{'fraction_non_alphanumeric': 0.029357269321035946, 'fraction_numerical': 0.01694322065302425, 'mean_word_length': 4.782049755981431, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multiphysics Lattice Discrete Particle Modeling (M-LDPM) for the Simulation of Shale Fracture Permeability

A three-dimensional Multiphysics Lattice Discrete Particle Model (M-LDPM) framework is formulated to investigate the fracture permeability behavior of shale. The framework features a dual lattice system mimicking the mesostructure of the material and simulates coupled mechanical and flow behavior. The mechanical lattice model simulates the granular internal structure of shale and describes heterogeneous deformation by means of discrete compatibility and equilibrium equations. The network of flow lattice elements constitutes a dual graph of the mechanical lattice system. A discrete formulation of mass balance for the flow elements is formulated to model fluid flow along cracks. The overall computational framework is implemented with a mixed explicit-implicit integration scheme and a staggered coupling method that makes use of the dual lattice topology enabling the seamless two-way coupling of the mechanical and flow behaviors. The proposed model is used for the computational analysis of shale fracture permeability behavior by simulating triaxial direct shear tests on Marcellus shale specimens under various confining pressures. The simulated mechanical response is calibrated against the experimental data, and the predicted permeability values are also compared with the experimental measurements. Furthermore, the paper presents the scaling analysis of both the mechanical response and permeability measurements based on simulations performed on geometrically similar specimens with increasing size. The simulated stress-strain curves show a significant size effect in the post-peak due to the presence of localized fractures. The scaling analysis of permeability measurements enables prediction of permeability for large specimens by extrapolating the numerical results of small ones.


Introduction
Gas shale and shale-like rock play a crucial role in several engineering applications such as hydraulic fracturing and CO 2 sequestration, requiring a basic understanding and characterization of their material properties including mechanical and hydraulic behavior. Particularly, it is essential to investigate the fracture-and damage-induced permeability evolution of the material for two opposite purposes: (1) utilizing the low permeability of shale in subsurface energy storage; (2) enhancing gas recovery from unconventional reservoirs. This involves quantifying local damage and fracturing of material and correlating permeability evolution with damage by taking into account the multiscale character of shale (Hyman et al, 2016;Ilgen et al, 2017;Li et al, 2016).
Previous experimental studies have shown that fracture permeability involves multiple physical processes, and is influenced by various factors such as fracture aperture and roughness (Witherspoon et al, 1980;Zimmerman and Bodvarsson, 1996;Yeo et al, 1998), chemical precipitation and dissolution (Detwiler, 2010;Noiriel et al, 2010;Elkhoury et al, 2013), and thermal effects (Rutqvist et al, 2008). In particular, open and connected cracks have a strong influence on permeability. A large amount of previous work on fracture permeability has been conducted on prefractured (split or sawn) specimens and has shown that rock permeability is strongly related to density, spacing, orientation, width, and length of fractures within specimens (Gutierrez et al, 2000;Davy et al, 2007;Jobmann et al, 2010;Cho et al, 2013).
However, researchers suggested that permeability studies of rocks fractured at in situ conditions through triaxial compression or direct shear tests can better represent fracture properties in the subsurface than prefractured or artificial specimens (Zhang et al, 2013;Nygård et al, 2006;Carey et al, 2015a;Frash et al, 2016;Su et al, 2017;Kluge et al, 2017). Indeed, some factors such as the influence of stress state during fracture, the transition from brittle to ductile failure, effective stress and hysteresis effects can be only investigated under in situ conditions. An experimental campaign on shale fracture permeability using triaxial direct-shear coreflood methods provides a unique opportunity to continually measure permeability during the entire loading phase under in situ conditions (Carey et al, 2015a,b;Frash et al, 2016Frash et al, , 2017. Integrating experimental studies on shale fracture-permeability behavior into a numerical model is meaningful from both theoretical and practical points of view, considering that characterization of fracture permeability is difficult and costly in the field and in the laboratory. Modeling of fracture permeability remains a challenging problem in the following three aspects: (1) description of cracking and fracturing events considering material heterogeneity; (2) proper selection of a fluid solver and an efficient coupling scheme to simulate coupled fluid flow in the damaged material; (3) capability of dealing with scaling effects.
Continuum-based approaches, in which the material domain of interest is treated as a homogeneous continuum body or a multiphase continuum mixture if the materials are fluid-infiltrated, are widely used, but are inherently incapable of dealing with material heterogeneity and they fail to a large extent to simulate appropriately damage localization and the resulting size-effect. In these approaches, the fracture surfaces are usually represented by various methods, including but not limited to, embedded discontinuities (Mohammadnejad and Khoei, 2013;Salimzadeh and Khalili, 2015;Jin and Zoback, 2017), phase field models (Mikelic et al, 2015;Miehe et al, 2015;Choo and Sun, 2018), or smeared crack bands (Chau et al, 2016;Li et al, 2017a).
On the contrary, approaches forfeiting the continuum assumption (also called discrete models), such as the discrete element method (Cundall and Strack, 1979;Potyondy and Cundall, 2004) and lattice/particle models (Cusatis et al, 2011a(Cusatis et al, ,b, 2003aSmith et al, 2014;Li et al, 2017e;Ashari et al, 2017;Kozicki and Donzé, 2008;Belheine et al, 2009;Bolander Jr and Saito, 1998;Nikolic et al, 2016), considers the material domain as a three-dimensional system of interacting cells or particles often connected through lattice elements. These approaches come with a natural representation of fractures and cracks by discontinuities between adjacent cells or particles and, if formulated at the appropriate length scale, are able to reproduce critical features of material heterogeneity.
Furthermore, the fluid transport and hydro-mechanical coupling effect can be simulated at different length scales. At the microscopic level, fluid flow can be solved by Navier-Stokes equations with certain assumptions and the hydro-mechanical coupling effect can be captured via force-based interaction across the fluid-solid interface (Catalano et al, 2014;Robinson et al, 2014). Although conceptually simple, the microscopic approach requires a knowledge of the pore network structure and places a high demand on computational resources. An alternative to the costly microscopic approach is a mesoscopic/macroscopic approach (Bolander and Berton, 2004;Grassl et al, 2015;Wang and Sun, 2016;Li et al, 2017c,b) which adopts diffusion equations to describe fluid flow in a porous medium and relies on phenomenological laws such as Darcy's law, Fick's law, and lubrication theory. By introducing diffusion coefficients, such as hydraulic conductivity, permeability, and diffusivity in the sense of homogenization over the multiphase mixture, this approach benefits from affordable computational cost. The choice of the length scale and the corresponding flow solver involves a trade-off between computational cost and physical realism.
Finally, discrete models can reproduce well size effect phenomena. It has been widely demonstrated that various localization events (tensile strain softening, shear band, compaction band, etc.) result in a strong size effect, that is the dependence of macroscopic material properties (such as, for example, strength and fracture toughness) on sample size (Bažant, 1984;Bažant and Planas, 1997;Li et al, 2017d;Jin et al, 2017b). However, not enough attention has been paid to the effect of sample size on rock fracture permeability. Indeed, this aspect is extremely important for models that are calibrated and validated on the basis of laboratory measurements and then used for field applications.
In this work, a three dimensional (3D) numerical model based on a Multiphysics Lattice Discrete Particle Model (M-LDPM) framework is formulated to study shale fracture permeability behavior by simulating triaxial direct-shear experiments.
2 Multiphysics Lattice Discrete Particle Model (M-LDPM) for shale M-LDPM is formulated in a discrete poromechanics setting by adopting two coupled dual lattices simulating mechanical and transport behaviors, respectively. The model adopts an "a priori" discretization of the internal structure of the material at the meso-scale, which is the length scale of major material heterogeneities.

Dual lattice topology
Typical rock materials feature an internal granular structures. For shale, grains consist of agglomerated clay particles and other minerals such quartz and feldspar, and in most cases are characterized by a maximum grain size in the order of 50 µm as one can see in Fig. 1a with reference to the Toarcian shale (Akono and Kabir, 2016). In the M-LDPM formulation, the grain-scale heterogeneity of shale is replicated by discretizing the material domain with a system of randomly distributed polyhedral cells where a polyhedral cell represents a binder-coated grain. Only the coarser grains are considered in order to reduce computational cost while preserving the main source of heterogeneity.
The granular lattice system is constructed through the following steps: (1) An artificial supporting system with spherical particles placed at the center of shale grains inside a certain volume of material is generated by a try-and-reject procedure devised to prevent mutual overlapping of particles and overlapping with the external boundary (see Fig. 1b for 2D representation and Fig. 1e for a 3D system). Details of the procedure can be found in Cusatis et al (2011b). The size distribution of the spherical supports follows a fractal Cumulative Distribution Function (CDF): are particle size, minimum particle size, maximum particle size, respectively. In this work q = 2.5 was used in all calculations. Particle diameters d i are calculated by sampling the CDF. New particles are generated until the total volume of the generated spherical particles,Ṽ a0 = i (πd 3 i /6), exceeds V a0 , where V a0 = v a0 V is the total volume of the simulated particles, v a0 = [1 − (d 0 /d a ) 3−q ]v a can be computed from the particle volume fraction v a , and V denotes the specimen volume. Note that d 0 , d a , and v a need to be calibrated against the grain size distribution of the simulated granular rock. The set of particles inside the interior volume is also augmented through surface nodes with zero radius and a given spacing to describe the external surface.
(2) A constrained Delaunay tetrahedralization connects the given particle centers and fills the space without gaps. For 2D illustration, the tetrahedralization corresponds to a triangular mesh, as shown in Fig. 1c, whereas in an actual 3D material domain, a mesh of tetrahedra is created, as shown in Fig. 1f. (3) A 3D domain tessellation, anchored to the Delaunay tetrahedralization but not coinciding with the classical Voronoi tessellation, subdivides the domain into a system of polyhedral cells (see Fig.   1d for 2D representation and Fig. 1g for the 3D system). A generated polyhedral cell, as illustrated in Fig. 1h, represents a rock grain. The surface of each polyhedral cell is subdivided into triangles, called hereinafter "facets" and describing potential crack locations (see Fig. 1e). Each facet has one vertex inside the tetrahedron (tet-point), one on each edge of the tetrahedron (edge-point), and one on a face of the tetrahedron (face-point). The tessellation procedure is described in detail in Cusatis et al (2011b) and it is formulated in a way to minimize the intersection between each facet and the grains. This is done so that, as often verified in practice, cracks are simulated to occur at the grain interface and in the embedding matrix as opposed to cutting through the grains.
While the LDPM cell system provides a geometrical characterization of the solid portion at the mesoscale and it serves as a structure for the description of the mechanical behavior, a dual lattice system allows for the coupled simulation of fluid flow at the same length scale. The dual lattice system is constructed as follows. Let's consider the 3D mesh of tetrahedra generated in step (2) as described above and its 2D representation as depicted in Fig. 2a. A Flow Lattice Element (FLE) connects two points inside two adjacent tetrahedra (tet-points), which are also vertices of the corresponding polyhedral cells, as shown in It is worth pointing out that the term "dual lattice" represents the topological relation between the LDPM cell system for the solid phase and the network of FLE for the fluid phase. This is similar to the dual-graph concept adopted by Grassl et al (2015); Grassl and Bolander (2016); Ulven and Sun (2017).
Under such an arrangement, the flow/transport elements are aligned in the general direction of potential crack paths, and thus the influence of local cracks on the fluid flow, and vice versa, the influence of the fluid pressure on the solid, can be captured naturally. With reference to Fig. 2a (for 2D) and 2c (for 3D), one can see that each FLE is associated to a certain characteristic number of facets (2 in 2D and 6 in 3D) belonging to the two LDPM cells associated to the mechanical lattice elements P 1 − P 2 for 2D, and P 1 − P 2 , P 1 − P 3 , P 2 − P 3 for 3D. This makes possible the seamless two-way coupling of the mechanical and flow behaviors as will be presented later in this paper.

Kinematics
In the LDPM formulation, adjacent cells interact through shared triangular facets. Rigid-body kinematics is adopted to describe the heterogeneous deformation of the cell system. Three strain measures, one normal component and two shear components, are defined at each facet (see Fig. 1h) as where u C is the displacement jump vector calculated by the displacements and rotations of the nodes adjacent to the selected facet; = tetrahedron edge associated with the facet; n is a unit vector normal to each facet, and m and l are two mutually orthogonal unit vectors contained in a plane orthogonal to n. It is worth pointing out that (1) to avoid unrealistic shear locking phenomena (Cusatis et al, 2011b), the local system of reference is defined on the projection of the LDPM facets in a plane orthogonal to the lattice connecting the two cells; and (2) the strain definitions in Eq. 1 correspond to the projection into the facet system of reference of the classical strain tensor Cusatis and Zhou, 2013; Rezakhani and Cusatis, 2016).

LDPM constitutive laws
The facet stress vector applied to the solid phase, t s = t N n + t M m + t L l, is calculated through appropriate constitutive laws, which describes various phenomena of the solid grain interaction (Fig. 1g).
In the elastic regime, the mechanical facet stress components are proportional to the corresponding strain components: where E N is the normal modulus, E M = E L = αE N is the shear modulus, and α = shear-normal coupling parameter.
Beyond the elastic limit, the vectorial constitutive laws are formulated to reproduce two distinct nonlinear phenomena.
The first source of nonlinearity is related to fracturing and cohesive behavior under tension and tension/shear that occur for N > 0. Following the work of Cusatis et al (2003a,b), one can define an effective strain and an effective stress as = compute the normal and shear stresses as follows: The effective stress t is incrementally elastic (ṫ = E N˙ ) and must satisfy the inequality 0 ≤ t ≤ σ bt ( , ω), in which The function σ 0 (ω) represents the strength limit for the effective stress, and is given by where r st = σ s /σ t is the ratio between the mesoscale shear strength (or cohesion) σ s and the mesoscale tensile strength σ t . After the maximum effective strain reaches its elastic limit 0 (ω) = σ 0 /E N , the stress boundary σ bt decays exponentially with a softening modulus, H 0 (ω) = H t (2ω/π) nt , which provides a smooth transition from softening behavior under pure tension (ω = π/2), H 0 = H t , to perfectly plastic behavior under pure shear (ω = 0), H 0 = 0. In most cases, n t = 0.2. For the correct energy dissipation during mesoscale damage localization to be preserved, the softening modulus in pure tension is expressed energy, and is the length of the tetrahedron edge (lattice element) associated with the facet.
Since relevant to the the effect of cracking on permeability discussed later, it is worth pointing out that the mesoscale crack opening components (normal and shear) can be calculated as The second source of nonlinearity is related to frictional behavior under compression/shear for N < 0.
This can be simulated effectively through a non associative incremental plasticity formulation in which, and the normal stress is simply elasticṫ N = E N˙ N . In the previous expressions, is the yielding function,˙ P M =λ∂ψ/∂t M and˙ P L =λ∂ψ/∂t L are the shear plastic strain increments, and λ is the plastic multiplier, ψ = ψ 0 [t 2 M + t 2 L ] 1/2 is the plastic potential. The yielding function can be formulated according to a Mohr-Coulomb criterion as in which µ 0 is the internal friction coefficient.
The facet stresses calculated through the constitutive laws described above represent the stresses carried by the solid phase. Equilibrium considerations at the facet level allow for the reasonable assumption of a parallel coupling between the stresses carried by the solid phase and those by the fluid phase. In this work, the effective stress concept is adopted, and the total stress vector on each facet can be computed as where b is the Biot coefficient, and t w = pn is the fluid pressure vector. The negative sign in Eq. 8 comes from the pressure sign convention, which is positive for the fluid and negative for the solid.
At this point a few comments are in order.
1. The constitutive equations described above depend on 6 parameters that must be identified by fitting macroscopic experimental data relevant to elastic behavior (for E N and α), tensile fracturing (for σ t and t ), unconfined compression and triaxial compression (for σ s and µ 0 ).
2. The resulting macroscopic mechanical response is statistically isotropic. Material anisotropy, typical of many shale formations, can be accounted for by postulating the dependence of the above parameters on the facet orientation (Jin et al, 2017a) and by explicitly simulating bedding planes. This was already done by Li et al (2017e) and it is not adopted here simply because the particular shale samples analyzed in this paper do not show appreciable anisotropy.
3. The constitutive equations do not include nonlinearity for high confining pressures under hydrostatic loading, pressure dependence on the elastic behavior, and creep behavior. These aspects, albeit important, are outside the scope of this study.

Equation 8
is formulated based Biot's theory of poroelasticity (Biot, 1941) and for intact shale one can assume b ≈ 0.6 (Suarez-Rivera and Fjaer, 2013). Beyond the elastic regime, the effective stress concept still holds, but the Biot coefficient needs to evolve with the level of damage (Chau et al, 2016) and must be equal to one (b = 1) at the location of a stress free crack (t s = 0). Such formulation, however, requires calibration with data not available in the current study. Hence the constant value b = 0.6 will be used in all simulations. It must be mentioned that this assumption does not reduce the validity of the results because in the analyzed experiments the fluid pressure is much smaller than the stresses on the solid phase (|t s | p), making the dependence of the response on b negligible.

Equilibrium
The equilibrium equations are obtained through the force and moment equilibrium of each cell subject to the force resultants obtained by multiplying the tractions in Eq. 8 times the facet areas for all facets belonging to the given cell. In absence of body forces, one can write where F I is the set of facets surrounding the node P I (located inside cell I) and obtained by collecting all the facets associated with each adjacent cell pairs (I, J); A f = projected facet area; c I = vector connecting the node P I to the facet centroid; M I u , M I θ = mass and moment of inertia of cell I; andÜ I , Θ I are acceleration and rotational acceleration, respectively, of cell I. In the current implementation, an explicit dynamic algorithm (based on a central difference scheme) is adopted to solve the equations above (Cusatis et al, 2011b). This offers the advantage of avoiding the convergence problems that implicit schemes often have in handling softening behaviors.

Discrete formulation of fluid flow
This section discusses the formulation of flow phenomena in a discrete setting consistent with the LDPM framework for shale introduced earlier. Attention is restricted to the flow of water under full saturation conditions, constant room temperature, and assuming that water behaves as a slightly compressible Newtonian fluid.
Let's consider two adjacent tetrahedra, each connecting four LDPM cells in the undeformed configuration: P 1 , P 2 , P 3 , P 4 and P 1 , P 2 , P 3 , P 5 in Fig. 3a. They have a common triangular face 123 across which a FLE connects the tet-points T 1 and T 2 located inside the two tetrahedra. The FLE is labeled E 12 and it is associated to two pyramidal volumes, V 1 and V 2 , identified by the points P 1 , P 2 , P 3 , T 1 and P 1 , P 2 , P 3 , T 2 .
The volumes can be computed as V 1 = A n L 1 /3, and V 2 = A n L 2 /3, where A n = Aa T e is the area calculated by projecting the original area of the face 123 , A, into a plane orthogonal to E 12 ; a is the unit vector orthogonal to the face 123 ; e is the unit vector oriented as E 12 ; L 1 , L 2 are the distance between point F and T 1 , T 2 , respectively; point F is the point where the triangular face 123 intersects E 12 . It is worth noting that L 1 + L 2 = L, where L is the length of E 12 . If one defines the fraction coefficients g 1 = L 1 /L and g 2 = L 2 /L (g 1 + g 2 = 1), one obtains represents the total volume of material associated with E 12 .

Water content and flux for uncracked shale
With reference to the control volume V i (i=1,2), the mass of water can be written as M i wu = m wi V i in which m wi is the water mass content, defined to be the water mass per unit reference volume (see, among others, Rice and Cleary (1976); Wang (2000)). The change in water mass content can be related to the increment of water content ζ i = (m wi − m w0 )/ρ w0 , where m w0 and ρ w0 are water mass content and density in the reference state, respectively. For slightly compressible fluids, the water density in the current state can be related to ρ w0 by defining the bulk modulus K w ; one can write where p i is the fluid pressure in V i , and p 0 is the initial/reference pressure.
According to the classic theory of poromechanics (Biot, 1941;Rice and Cleary, 1976;Detournay and Cheng, 1993), the increment of water content, ζ i , can be expressed as a linear combination of the volumetric strain, ε i , of the solid phase defined as the relative variation of the solid volume, and the water pressure, where M b denotes the Biot modulus (also defined as the reciprocal of the so-called storage coefficient). It is worth pointing out that b and M b may vary due to material hetereogeneity. The effect of this variation is insignificant in the context of this paper and will be neglected thereinafter. Also, the previous discussion on b and its evolution with damage also holds for M b . Again, in absence of relevant experimental data a constant value is actually used in this paper.
One can write the time variation of the water mass in the control volume aṡ The mass flux through the face 123 from V 1 into V 2 can be obtained by using Darcy's law, which can be written as where κ 0 and µ w denote the intrinsic permeability of the material and the water viscosity, respectively; p 1 , p 2 are the values of fluid pressure at points T 1 and T 2 , respectively; andρ w = g 1 ρ w1 + g 2 ρ w2 is an estimate of the average density in the volume V .
Although Darcy's law is classically adopted to calculate flow rate through a porous material, as pointed out by Falk et al (2015); Obliger et al (2018), it might not be completely appropriate to predict transport in shales due to the strong adsorption in kerogen and the breakdown of hydrodynamics at the nanoscale.
Analysis of non-darcy behavior of shale is beyond the scope of the current paper. Moreover, this aspect does not significantly affect the results since shale fracture-permeability is mainly governed by the flow along cracked surfaces, as will be discussed below.

Water content and flux for cracked shale
The influence of local cracking can be addressed by introducing the contribution of the water mass stored in the cracked domain and the water flux along the crack surfaces. Let's consider again E 12 and the six associated LDPM facets of areas A i f 1 , A i f 2 , and A i f 3 belonging to V i (i=1,2) as it is illustrated in Fig.   3b. The water mass stored in the cracks is M i wc = ρ wi V ci , where the cracked volume can be expressed as The time variation of the water mass in the cracks can be written aṡ The water mass flux, Q 12 wc , from V 1 into V 2 associated with the cracks, can be approximated by assuming a steady laminar flow between two crack surfaces with a cross section of length f j and width δ i N j with j = 1, ..., 3, where f j represents the intersection of the j th facet with the tetrahedron face 123 (see Fig. 3c and d). In this case, the solution of a two dimensional Poiseuille flow in a channel (Massey and Ward-Smith, 1998), known as Poiseuille's formula, can be adopted.
One can write where κ c = 1 12A and is derived assuming that the cracked permeabilities in V c1 and V c2 are coupled in series.

Mass balance equations
The total water mass and total water flux can be obtained by adding the contributions from the uncracked and cracked domains. In this case, the mass balance equations for volume V 1 and V 2 can be written as and bε 2 +ṗ where both equations are normalized with the reference water density ρ w0 , ξ =ρ w (κ 0 + κ c )/(ρ w0 µ w ), and the superimposed dot represents time derivative. Equations 16 and 17 can be rewritten in matrix form as where C and D are the capacity (or storage) matrix and the conductance matrix of the FLE E 12 , respectively; S is the source term; and p = [p 1 , p 2 ] T . One has It is worth observing that, to reduce the memory requirements of the calculations, it is possible, without significant difference in the results, to substitute c i and s i (i = 1, 2) with the volume averagesc = g 1 c 1 + g 2 c 2 ands = g 1 s 1 + g 2 s 2 . The global matrices of the overall governing equations of the flow problem are obtained by assembling the matrix contributions of all flow elements. It is worth mentioning that the contributions of the flow elements associated with the boundary layers are also included. Dirichlet type boundary conditions can be imposed on the external nodes of the boundary layers through the penalty method. The fluxes of the boundary layer elements can be also calculated.
Thanks to their perpendicularity to the external surfaces, the sum of the fluxes (Riemann sum) provides an estimate of the volumetric flow rate over the boundaries. In this work, the time integration of the flow problem was performed by means of the Crank-Nicolson method (Di Luzio and Cusatis, 2009; Jirásek, 2017). The staggered coupling scheme along with the mixed explicit-implicit integration scheme provides a simple yet efficient approach for two-way coupling simulations of mechanical and flow behaviors.
The proposed framework is implemented into the MARS software (Pelessone, 2006), which is a structural analysis computer code with an object-oriented architecture that makes the implementation of new computational technologies very effective.

Numerical simulations of triaxial direct shear tests
The M-LDPM framework is applied here to the numerical analysis of shale fracture permeability.

Simulation setup
The simulated test setup is the one of recent triaxial direct shear tests performed at Los Alamos National Laboratory (LANL) to study shale fracture permeability behaviors at in situ stress conditions and considering various effects such as confining pressure, in situ stress, time-dependent factors, and stress cycling (Carey et al, 2015a;Frash et al, 2016Frash et al, , 2017. The newly improved triaxial direct shear apparatus developed by Frash et al (2017) provides better stress control by limiting specimen rotation in direct shear loading, and thus improves the accuracy and sensitivity of permeability measurements.
A cylindrical specimen with length H, diameter D, and a 1:1 H to D ratio is placed between two  opposing semi-circular disks, as depicted in Fig. 4a. The disks are stiff enough so that can be simulated as rigid bodies with only vertical displacement allowed and are in direct contact with the upper-half and lower-half surface of the specimen. The load is applied through the upper disk, while the lower one is fixed. The specific loading configuration creates a plane of direct shear extending vertically through the specimen. In addition, a confining pressure is applied laterally on the external surface of the specimens and on the two upper and lower half surfaces not in contact with the loading platens. The resultant boundary stresses are qualitatively illustrated in Fig. 4b as per the description provided by Frash et al (2017). As one can see there is stress concentration, σ p , near the direct shear plane and part of the upper and lower surfaces are not loaded. This is due to the fact that the specimens ends tend to rotate under load and they detach partially from the loading platens.
In the simulations, this is obtained by using a unilateral, frictional penalty contact model applied between the rigid disks and the specimen. The penalty contact model prevents compenetration of the specimen and platens but allows specimen to platen detachment. Hence, it is able to reproduce the experimental conditions and to replicate the stress state illustrated in Fig. 4b. The penalty contact algorithm also includes simulation of Coulomb frictional effects. The static and kinetic coefficients of friction used in the simulations are 0.2 and 0.15, respectively.
In the experiments, water was pumped onto the entire end face of the specimen for permeability measurement. In the simulations, inlet and outlet pressures, p in and p out , respectively, were imposed on the nodal variables of the top and bottom boundary layers of the flow lattice model to create a differential pressure, ∆p = p in − p out , across the core, as shown in Fig. 4c. The corresponding upstream and downstream volumetric water flow rates were recorded continually during the entire test run. The apparent macroscopic permeability across the core area was calculated using Darcy's law applied to the entire sample. Due to dynamic processes occurring during the experiments (deformation, fines migration, etc.), permeability values changed slowly. The reported values represent pseudo-steady state where upstream and downstream flow rates were nearly equal (see Frash et al, 2016Frash et al, , 2017. Three different effective confining pressures, σ c = 4.4, 9.3, and 14.5 MPa, were considered in the experiments and simulations. Similar to the experiments, the effective confining pressure was calculated using the effective stress concept, i.e. σ c = σ c − 0.5b(p in + p out ), where σ c is the actual confining pressure applied on the specimen surfaces. A constant value of 0.6 was used for the Biot coefficient b in this work.
The simulations were performed following the experimental sequence provided by Frash et al (2017).
Firstly, the confining pressure was gradually increased to a target value while a minimal positive direct shear stress (up to 1.0 MPa) was maintained. Secondly, the axial displacement applied to the top semicircular disk was steadily increased while maintaining constant confining pressure until major fracture creation was identified. Finally, loading ceased when the axial stress reached a plateau (residual stress).

Mechanical behavior
The first step in the simulations is related to the LDPM parameter calibration.
The generation of the shale mesostructure depends on the minimum and maximum size of the support particles (see Sec. 2.1), which, without specific information on the shale under study, were estimated (Li et al, 2017e) as d 0 =15 µm and d a =35 µm. These values along with v a = 0.55 lead to a LDPM cell size distribution characterized by minimum, average, and maximum cell size equal to 15, 30, and 51 µm, respectively. In addition, the shale under study was reported to be statistically isotropic (Frash et al, 2017) since no anisotropy, pre-exiting flaws or bedding planes were observed in the tested specimens.
The experimental data reports an elastic modulus of 74 ± 6 GPa and a Poisson ratio of 0.30 ± 0.02.
Since the samples were saturated with water and the tests were completed in a relatively short time period, the reported elastic properties must be considered to be undrained values. However, since the LDPM parameters govern the mesoscale mechanical behaviors of the solid phase, they must be related to the macroscopic properties measured under drained condition. Given a reasonable assumption of Biot and Skempton coefficients equal to 0.6 for typical shale (Suarez-Rivera and Fjaer, 2013), the drained elastic properties can be obtained by using the classic theory of poroelasticity (Detournay and Cheng, 1993;Wang, 2000  Finally, the shear strength, σ s = 24 MPa, and the internal friction coefficient, µ 0 = 0.18, were calibrated by simulating the direct shear tests with various level of confinement. In order to save computational time, simulations were performed on reduced size samples. This was possible because under this loading condition the peak stress is basically size-independent. In order to demostrate this statement, specimens geometrically similar to the experimental ones but varting sizes were simulated: they had diamater equal to D = 0.5, 0.8, 1.2, and 1.6 mm. Numerical results of the stress-strain curves for the geometrically similar specimens with increasing size are shown in Fig. 5a, b and c for various confining pressures. In these plots, similarly to Frash et al (2017), the net direct-shear stress acting on the direct-shear plane is calculated as where F P is the force applied on the rigid disk, and F c is the resultant force of the confining pressure applied on half of the specimen cross section. Also, the nominal strain is calculated as γ = u/D, where u is the displacement of the upper loading platen.
For each confining pressure, the curves for different sizes coincide in the elastic and nonlinear regime before reaching the peak. The peak stress is only minimally influenced by the sample size. The post-peak regime, however, is markedly different and the post-peak stress vs. strain curve tends to be steeper as the specimen size increases. This behavior is typical of quasi-brittle materials and can be explained by the size effect induced by damage localization and fracture (Bažant, 1984;Bažant and Planas, 1997). Indeed, as τ DS approaches the peak, coalescence of microcracking occurs and finally leads to a localized fracture along the direct shear plane as one can see in Fig. 4d and 4e for a typical numerical response.
If one extrapolates the numerical scaling results, snap-back at peak load with a loss of stability even under displacement control is expected for sufficiently large specimens. This indeed was observed in the experiments where the samples featured a catastrophic failure with a sudden drop of the applied load after attaining the maximum value (Frash et al, 2017). In addition, one can see in Fig. 5a, b, and c that, similarly to typical experimental data, the specimens under higher effective confining pressure exhibit a more ductile response and higher residual stress after failure. and c M C = 26 ± 2 MPa reported in Frash et al (2017). The LDPM parameters were calibrated to slightly overestimate the experimental data in order to compensate for the size-effect discussed above.
The tested and simulated samples feature lateral expansion as a result of shear-dilation behavior.
This can be characterized by the so-called mechanical aperture calculated by the difference of the lateral displacements between two diametrically opposite points located at half height of the specimens. The mechanical aperture results are shown in Fig. 6. It is clear that the mechanical aperture depends significantly on the confining pressure but its evolution versus the loading platen displacement is only minimally affected by the specimen size. The latter observation is key for the size scaling of the fracture permeability that will be discussed in the next section.

Fracture permeability behavior
The numerical predictions of the permeability measurements for the simulated specimens during the triaxial direct shear tests are presented in this section.
The dynamic viscosity of water, µ w = 8.9 × 10 −4 Pa, the water bulk modulus, K w = 2.15 GPa, and the water reference density, ρ w0 = 1000 kg/m 3 , are values at room temperature and atmospheric pressure. In addition, the flow model depends on the intrinsic permeability and the Biot modulus. The value κ 0 = 10 −20 m 2 =10.13 nD was estimated from the typical permeability of intact shale samples measured in laboratory (Soeder, 1988), and M b = 61, 670 MPa was estimated from the experimentally reported elastic properties by using the classic theory of poroelasticity (Detournay and Cheng, 1993;Wang, 2000).
The typical upstream and downstream volumetric flow rate (discharge) data collected during one generic simulation run are plotted against the loading platen displacement in Fig. 7a, and are used to calculate the macroscopic apparent permeability through Darcy's law as follows where Q is the absolute value of the downstream flow rate (downstream = -upstream when attaining steady state). Figure 7b shows the typical permeability calculated using Eq. 22 (solid line) from the numerical The results in Fig. 7 are relevant to a sample with D = 1.2 mm.
In the elastic regime and before the peak the apparent permeability is very small (in the order of 10 nD). The mechanical aperture is negligible as well. As the direct shear stress reaches the peak, a rapid rise of the mechanical aperture is observed, followed by a substantial increase of permeability. This is strongly associated with the development of localized fractures along the shear plane, as illustrated in Fig. 4e.
The localized fractures offer wide open channels permitting fluid flow. Correspondingly, the main flow is concentrated within the region of fracture localization, as shown in Fig. 4f.
The simulated permeability as well as the mechanical aperture keeps increasing gradually in the postpeak softening regime. This behavior is, of course, different from what observed in the experiments due to size effect. As explained earlier, the experiments feature a sudden drop of the stress as opposed to a gradual post-peak. As a result, in the experiments, the permeability spiked sharply upward as soon as failure occurred. This difference must be taken into account when comparing the computed fracture permeability values with the experimental ones.
In the presence of a single fracture, the fracture permeability is usually considered to be proportional to the square of mechanical aperture according to Poiseuille's formula. As shown in Fig. 7b, the variations of the apparent permeability and the square of the measured mechanical aperture do not match. This indicates that global geometrical measurements, such as mechanical aperture, and fracture permeability are not closely correlated. In order to provide a sufficiently accurate prediction, a fine scale coupled model, such as M-LDPM presented in this work, is necessary to track the evolution of permeability as a consequence of interaction between local failure events and flow in three-dimension. Figure 8 presents the numerical results of the apparent permeability as a function of the loading platen displacement for the simulated specimens of various sizes and with different effective confining pressures.
If one compares the numerically calculated apparent permeability values at a certain level of displacement, it can be seen that for the specimens with the same size, the permeability values decrease with an increase of effective confining pressure σ c . The same trend was observed in the experiments and many other rock fracture permeability studies (Gutierrez et al, 2000;Davy et al, 2007;Jobmann et al, 2010;Zhang et al, 2013;Zhou et al, 2016). This is consistent with the previously discussed trend of decreased mechanical aperture for increased effective confining pressure.
Furthermore, one may note that the numerical results of the apparent permeability also depend on specimen size. A comparison of permeability for the specimens with different sizes yet under a same level of confinement shows that permeability tends to decrease with an increase of specimen size. The size dependence of the permeability measurements is caused by the fact that the damaged zone is localized and does not scale proportionally to the sample size.

Apparent permeability scaling analysis and predictions
Due to the size dependence of the apparent permeability, the comparison between numerical predictions and experimental results must be performed with reference to an appropriate scaling law.
A simple scaling law correlating the fracture permeability measurements obtained from different sized specimens can be formulated if one assumes that (1) the width of the permeable band in which fracture localizes, h, is independent of specimen size; (2) The average opening or aperture of the permeable band is also size independent. In a direct shear test, such a permeable fracture band is indeed a shear band, and the opening or aperture of the shear band results from the shear-dilation behavior. For granular materials such as soil, these two assumptions are widely accepted: the shear band width depends mainly on grain size, and the shear-dilation is related to material constitutive behavior.
With reference to the simulation performed in the current study, the first assumption is supported by comparing the normal cracking profiles of the geometrically similar specimens with increasing size at the same level of displacement, as shown in Fig. 10a, b, c and d. Note that the cracks are localized into a narrow band, i.e. shear band, with a statistically constant width for all investigated specimens. The second assumption is verified by the results on the mechanical aperture discussed in the previous section.
The initial macroscopic permeability can be calculated approximately as D in ≈ κ 0 ; D in deviates from κ 0 due to the intrinsic heterogeneity of LDPM and the dual FLE system. However, since D in is calculated on a certain volume, the effect of the heterogeneity is filtered out for a large enough volume. Moreover, it worth pointing out that since K b is related to the mechanical aperture, it can be assumed not to change with size. Hence, since also the width of the permeable shear band h is size independent, the scaling law in Eq. 23 predicts a linear increase of the of the permeability with the inverse of the specimen radius.
Let's consider now apparent permeability values sampled at two levels of displacement, u 1 = 0.016 mm and u 2 = 0.024 mm. In both cases, the corresponding stresses are in the post-peak regime after fracture localization occurred. The corresponding mechanical apertures are about 0.013 mm and 0.017 mm for σ c = 14.5 MPa, which are in the range of the experimentally measured mechanical aperture right after failure.
In order to take into account the effect of the randomness of the internal structure (captured by the proposed model), the apparent permeability values were calculated for three specimens with different internal structures randomly generated by the proposed model for each size and for each effective confining pressure. The results are plotted in Fig. 9a, and b for the two displacement levels, respectively, as function of the inverse of the specimen radius. A consistent trend of decreased permeability for increased specimen radius can be clearly identified.
By performing a linear regression analysis of the generated numerical results, the constant B in Eq.
23 can be identified.
Once B is identified, the scaling law in Eq. 23 enables an extrapolation from the small size numerical simulations to the prediction of larger size specimens. Hence, considering the actual specimen radius of 12.7 mm in the experiments, one can predict the corresponding permeability values. The prediction results are listed in Table. 1.  (Frash et al, 2017) Comparison of the predicted permeability values with the initial fracture permeability values measured in the experiments right after failure (see 4 th column of Table 1) shows that the predicted values are of the same order of magnitude as the experimental data and follow the correct trend for increasing effective confining pressures. In general, the predicted values tend to be greater than the experimental values. This is probably due to the fact the formulated model simply formulates the meso-scale permeability to depend on the mesoscale normal crack opening but neglects the effect of the shear crack sliding which is known to reduce permeability (Kluge et al, 2017). Also, effects such as creep, chemical precipitation, particle mobilization, mylonitization and grain crushing induced during shearing, also play a role in the fracture permeability experiments and they are not included in the proposed model. However, a more sophisticated formulation requires specifically designed experiments which are not available at the moment.

Conclusions
A Multiphysics-Lattice Discrete Particle Model (M-LDPM) framework was formulated and used for the computational analysis of fracture permeability behavior of shale. This was obtained by simulating triaxial direct shear tests on Marcellus shale specimens of different sizes and under various confining pressures.
Based on the obtained results, the following conclusions can be drawn.
1) The proposed model can simulate the mechanical response of shale in the triaxial direct shear test, and successfully captured the poroelastic properties, the pressure-dependent strength properties, and the failure pattern. Benefiting from the model's capability of addressing grain-scale heterogeneity, the localized fracture and the size effect of the stress-strain responses can be captured accurately.
2) The dual lattice topology enables the seamless coupling of the mechanical and flow behavior of the material at grain-scale. As a consequence, the variation of the apparent permeability induced by fracturing processes can be simulated by relating the hydraulic properties of flow elements with local cracking events.
The numerical results demonstrated that a consistent trend of reduced fracture permeability for increased confining pressure can be captured properly by the proposed model in agreement with experimental evidence.
3) The scaling analysis of the permeability measurements obtained from the simulations performed on geometrically similar specimens enables a prediction of permeability for actual sized specimens by extrapolating the results of reduced scale specimens, which are required to make the computational cost of the simulations manageable. The predicted permeability values are in the same order of magnitude as the initial fracture permeability data measured in the experiments.