Center for Sustainable Engineering of Geological and Infrastructure Materials
[0.1in] Department of Civil and Environmental Engineering
[0.1in] McCormick School of Engineering and Applied Science
[0.1in] Evanston, Illinois 60208, USA
Multiphysics Lattice Discrete Particle Modeling (MLDPM) for the Simulation of Shale Fracture Permeability
[0.5in] Weixin Li, Xinwei Zhou, J. William Carey, Luke P. Frash, Gianluca Cusatis
[0.75in] SEGIM INTERNAL REPORT No. 183/976M
[0.75in]
Submitted to Rock Mechanics and Rock Engineering March 2018
Abstract: A threedimensional Multiphysics Lattice Discrete Particle Model (MLDPM) 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 explicitimplicit integration scheme and a staggered coupling method that makes use of the dual lattice topology enabling the seamless twoway 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 postpeak 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.
Keyword: Fracture permeability; Triaxial direct shear; Hydromechanical coupling; Lattice Discrete Particle Model; Dual lattice
1 Introduction
Gas shale and shalelike rock play a crucial role in several engineering applications such as hydraulic fracturing and CO 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 damageinduced 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 directshear 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, 2016, 2017).
Integrating experimental studies on shale fracturepermeability 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.
Continuumbased approaches, in which the material domain of interest is treated as a homogeneous continuum body or a multiphase continuum mixture if the materials are fluidinfiltrated, 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 sizeeffect. 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, b, 2003a, 2003b; Smith 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 threedimensional 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 hydromechanical coupling effect can be simulated at different length scales. At the microscopic level, fluid flow can be solved by NavierStokes equations with certain assumptions and the hydromechanical coupling effect can be captured via forcebased interaction across the fluidsolid 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 tradeoff 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 (MLDPM) framework is formulated to study shale fracture permeability behavior by simulating triaxial directshear experiments.
2 Multiphysics Lattice Discrete Particle Model (MLDPM) for shale
MLDPM 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 mesoscale, which is the length scale of major material heterogeneities.
2.1 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 as one can see in Fig. 1a with reference to the Toarcian shale (Akono and Kabir, 2016). In the MLDPM formulation, the grainscale heterogeneity of shale is replicated by discretizing the material domain with a system of randomly distributed polyhedral cells where a polyhedral cell represents a bindercoated 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 tryandreject 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): where , , are particle size, minimum particle size, maximum particle size, respectively. In this work was used in all calculations. Particle diameters are calculated by sampling the CDF. New particles are generated until the total volume of the generated spherical particles, , exceeds , where is the total volume of the simulated particles, can be computed from the particle volume fraction , and V denotes the specimen volume. Note that , , and 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 (tetpoint), one on each edge of the tetrahedron (edgepoint), and one on a face of the tetrahedron (facepoint). 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 (tetpoints), which are also vertices of the corresponding polyhedral cells, as shown in Fig. 2a for the 2D representation and Fig. 2c for the 3D system. The fluid mass transport between these two tetrahedra is described by the FLE. Hence, the fluid flow across the entire computational domain can be represented by a network of FLE connecting all tetrahedra, as illustrated in Fig. 2b and Fig. 2d. In addition, a thin layer of elements orthogonal to the external surface of the domain is generated to facilitate application of the boundary conditions needed for solving flow problems.
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 dualgraph 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 for 2D, and , , for 3D. This makes possible the seamless twoway coupling of the mechanical and flow behaviors as will be presented later in this paper.
2.2 LDPM for shale mechanical behavior
2.2.1 Kinematics
In the LDPM formulation, adjacent cells interact through shared triangular facets. Rigidbody 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
(1) 
where 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; is a unit vector normal to each facet, and and are two mutually orthogonal unit vectors contained in a plane orthogonal to . 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 et al, 2017; Cusatis and Zhou, 2013; Rezakhani and Cusatis, 2016).
2.2.2 LDPM constitutive laws
The facet stress vector applied to the solid phase, , 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:
(2) 
where is the normal modulus, is the shear modulus, and shearnormal 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 . Following the work of Cusatis et al (2003a, b), one can define an effective strain and an effective stress as , , respectively, and compute the normal and shear stresses as follows:
(3) 
The effective stress is incrementally elastic () and must satisfy the inequality , in which
(4) 
where , , and . The function represents the strength limit for the effective stress, and is given by
(5) 
where is the ratio between the mesoscale shear strength (or cohesion) and the mesoscale tensile strength . After the maximum effective strain reaches its elastic limit , the stress boundary decays exponentially with a softening modulus, , which provides a smooth transition from softening behavior under pure tension (), , to perfectly plastic behavior under pure shear (), . In most cases, . For the correct energy dissipation during mesoscale damage localization to be preserved, the softening modulus in pure tension is expressed as , where is the characteristic length, is the mesoscale fracture 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
(6) 
The second source of nonlinearity is related to frictional behavior under compression/shear for . This can be simulated effectively through a non associative incremental plasticity formulation in which, during plastic flow, , the incremental shear stresses are computed as , and the normal stress is simply elastic . In the previous expressions, is the yielding function, and are the shear plastic strain increments, and is the plastic multiplier, is the plastic potential. The yielding function can be formulated according to a MohrCoulomb criterion as
(7) 
in which 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
(8) 
where is the Biot coefficient, and 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.

The constitutive equations described above depend on 6 parameters that must be identified by fitting macroscopic experimental data relevant to elastic behavior (for and ), tensile fracturing (for and ), unconfined compression and triaxial compression (for and ).

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.

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 (SuarezRivera and Fjær, 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 () at the location of a stress free crack (). Such formulation, however, requires calibration with data not available in the current study. Hence the constant value 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 (), making the dependence of the response on negligible.
2.2.3 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
(9) 
where is the set of facets surrounding the node (located inside cell ) and obtained by collecting all the facets associated with each adjacent cell pairs ; = projected facet area; = vector connecting the node to the facet centroid; , = mass and moment of inertia of cell ; and , are acceleration and rotational acceleration, respectively, of cell . 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.
2.3 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: , , , and , , , in Fig. 3a. They have a common triangular face across which a FLE connects the tetpoints and located inside the two tetrahedra. The FLE is labeled and it is associated to two pyramidal volumes, and , identified by the points , , , and , , , . The volumes can be computed as , and , where is the area calculated by projecting the original area of the face , , into a plane orthogonal to ; is the unit vector orthogonal to the face ; is the unit vector oriented as ; , are the distance between point and , , respectively; point is the point where the triangular face intersects . It is worth noting that , where is the length of . If one defines the fraction coefficients and (), one obtains , and , in which represents the total volume of material associated with .
2.3.1 Water content and flux for uncracked shale
With reference to the control volume (=1,2), the mass of water can be written as in which 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 , where and 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 by defining the bulk modulus ; one can write
(10) 
where is the fluid pressure in , and 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, , can be expressed as a linear combination of the volumetric strain, , of the solid phase defined as the relative variation of the solid volume, and the water pressure, , as , where denotes the Biot modulus (also defined as the reciprocal of the socalled storage coefficient). It is worth pointing out that and 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 and its evolution with damage also holds for . 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 as
(11) 
The mass flux through the face from into can be obtained by using Darcy’s law, which can be written as
(12) 
where and denote the intrinsic permeability of the material and the water viscosity, respectively; , are the values of fluid pressure at points and , respectively; and is an estimate of the average density in the volume .
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 nondarcy behavior of shale is beyond the scope of the current paper. Moreover, this aspect does not significantly affect the results since shale fracturepermeability is mainly governed by the flow along cracked surfaces, as will be discussed below.
2.3.2 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 and the six associated LDPM facets of areas , , and belonging to (=1,2) as it is illustrated in Fig. 3b. The water mass stored in the cracks is , where the cracked volume can be expressed as and are the normal crack openings.
The time variation of the water mass in the cracks can be written as
(13) 
The water mass flux, , from into associated with the cracks, can be approximated by assuming a steady laminar flow between two crack surfaces with a cross section of length and width with , where represents the intersection of the facet with the tetrahedron face (see Fig. 3c and d). In this case, the solution of a two dimensional Poiseuille flow in a channel (Massey and WardSmith, 1998), known as Poiseuille’s formula, can be adopted.
One can write
(14) 
where
(15) 
and (=1,2). Equation 15 is derived assuming that the cracked permeabilities in and are coupled in series.
2.3.3 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 and can be written as
(16) 
and
(17) 
where both equations are normalized with the reference water density , , and the superimposed dot represents time derivative. Equations 16 and 17 can be rewritten in matrix form as
(18) 
where and are the capacity (or storage) matrix and the conductance matrix of the FLE , respectively; is the source term; and . One has
(19) 
(20) 
(21) 
where and
.
It is worth observing that, to reduce the memory requirements of the calculations, it is possible, without significant difference in the results, to substitute and () with the volume averages and .
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 CrankNicolson method (Di Luzio and Cusatis, 2009; Bažant and Jirásek, 2017).
Finally, the numerical implementation is completed by the application of a staggered coupling scheme with the solver of the mechanical problem. At the selected simulation time step, the mechanical model makes use of the solution (fluid pressure) from the fluid flow model, whereas the flow problem properties are updated by using the information from the mechanical model (crack opening and volumetric strain). The staggered coupling scheme along with the mixed explicitimplicit integration scheme provides a simple yet efficient approach for twoway 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 objectoriented architecture that makes the implementation of new computational technologies very effective.
3 Numerical simulations of triaxial direct shear tests
The MLDPM framework is applied here to the numerical analysis of shale fracture permeability.
3.1 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, timedependent factors, and stress cycling (Carey et al, 2015a; Frash et al, 2016, 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 , diameter , and a 1:1 to ratio is placed between two opposing semicircular 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 upperhalf and lowerhalf 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, , 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, and , respectively, were imposed on the nodal variables of the top and bottom boundary layers of the flow lattice model to create a differential pressure, , 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 pseudosteady state where upstream and downstream flow rates were nearly equal (see Frash et al, 2016, 2017).
Three different effective confining pressures, 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. , where is the actual confining pressure applied on the specimen surfaces. A constant value of 0.6 was used for the Biot coefficient 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).
3.2 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 and . These values along with lead to a LDPM cell size distribution characterized by minimum, average, and maximum cell size equal to 15, 30, and , respectively.
In addition, the shale under study was reported to be statistically isotropic (Frash et al, 2017) since no anisotropy, preexiting flaws or bedding planes were observed in the tested specimens.
The experimental data reports an elastic modulus of GPa and a Poisson ratio of . 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 (SuarezRivera and Fjær, 2013), the drained elastic properties can be obtained by using the classic theory of poroelasticity (Detournay and Cheng, 1993; Wang, 2000). One has GPa and for the drained Young’s modulus and the drained Poisson’s ratio, respectively. With these values, the LDPM normal modulus, MPa, and shearnormal coupling parameter, , can be estimated with the microtomacro formulas reported by Cusatis et al (2011b). The calculation was verified by simulating an undrained uniaxial compression test: the obtained undrained macroscopic Young’s modulus and undrained macroscopic Poisson’s ratio are 71 GPa and 0.29, respectively, which are in a good agreement with the experimental data. The LDPM tensile strength, MPa, was estimated based on previous work by the authors (Li et al, 2017e). The characteristic length, mm, was estimated (Cusatis and Schauffert, 2009) according to measurements of the effective Fracture Process Zone (FPZ) length for Marcellus shale (Li et al, 2017d).
Finally, the shear strength, MPa, and the internal friction coefficient, , 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 sizeindependent. In order to demostrate this statement, specimens geometrically similar to the experimental ones but varting sizes were simulated: they had diamater equal to , and mm.
Numerical results of the stressstrain 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 directshear stress acting on the directshear plane is calculated as where is the force applied on the rigid disk, and is the resultant force of the confining pressure applied on half of the specimen cross section. Also, the nominal strain is calculated as , where 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 postpeak regime, however, is markedly different and the postpeak stress vs. strain curve tends to be steeper as the specimen size increases. This behavior is typical of quasibrittle 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 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, snapback 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.
Figure 5d shows the comparison between the experimentally recorded and numerically calculated peak stresses vs. effective confining pressure. The numerical values are relevant to the simulations of samples with mm.
The macroscopic internal friction angle, , and cohesion, , of the simulated material, were calculated from the maximum values according to the classic MohrCoulomb criterion: and MPa. These values are in a good agreement with the experimental values , and MPa reported in Frash et al (2017). The LDPM parameters were calibrated to slightly overestimate the experimental data in order to compensate for the sizeeffect discussed above.
The tested and simulated samples feature lateral expansion as a result of sheardilation behavior. This can be characterized by the socalled 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.
3.3 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, Pa, the water bulk modulus, GPa, and the water reference density, kg/m, are values at room temperature and atmospheric pressure. In addition, the flow model depends on the intrinsic permeability and the Biot modulus. The value m= nD was estimated from the typical permeability of intact shale samples measured in laboratory (Soeder, 1988), and 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. 9a, and are used to calculate the macroscopic apparent permeability through Darcy’s law as follows
(22) 
where is the absolute value of the downstream flow rate (downstream = upstream when attaining steady state). Figure 9b shows the typical permeability calculated using Eq. 22 (solid line) from the numerical results as a function of the loading platen displacement. For comparison, the same figure reports curves relevant to the square of the mechanical aperture (dashed line) and the direct shear stress (dotted line). The results in Fig. 9 are relevant to a sample with 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 postpeak. 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. 9b, 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 MLDPM presented in this work, is necessary to track the evolution of permeability as a consequence of interaction between local failure events and flow in threedimension.
Figure 10 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 . 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.
3.4 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, , 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 sheardilation behavior. For granular materials such as soil, these two assumptions are widely accepted: the shear band width depends mainly on grain size, and the sheardilation 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. 14a, 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.
Macroscopically, the flow across the specimen core can be decomposed into one component through the permeable shear band and the other one through the intact matrix. Correspondingly, the total flow rate can be approximated by the sum of the flow rate through the shear band, , and the one through the intact matrix, , i.e. , where , , , , and , , are the apparent macroscopic permeability, the macroscopic permeability of the intact material, the macroscopic permeability of the shear band, respectively. This results in the following scaling law
(23) 
The initial macroscopic permeability can be calculated approximately as ; deviates from due to the intrinsic heterogeneity of LDPM and the dual FLE system. However, since 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 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 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, 0.016 mm and = 0.024 mm. In both cases, the corresponding stresses are in the postpeak regime after fracture localization occurred. The corresponding mechanical apertures are about 0.013 mm and 0.017 mm for 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. 13a, 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 in Eq. 23 can be identified.
Once 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.
, [MPa]  Apparent permeability, , [mD]  

Simulation 1  Simulation 2  Experiments  
2.8  14.1  45.4  6.6 
9.3  7.5  28.1  3.4 
14.5  5.1  21.5  0.57 
Numerical results sampled at displacement level = 0.016 mm  
Numerical results sampled at displacement level = 0.024 mm  
Initial fracture permeability (Frash et al, 2017) 
Comparison of the predicted permeability values with the initial fracture permeability values measured in the experiments right after failure (see 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 mesoscale 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.
4 Conclusions
A MultiphysicsLattice Discrete Particle Model (MLDPM) 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 pressuredependent strength properties, and the failure pattern. Benefiting from the model’s capability of addressing grainscale heterogeneity, the localized fracture and the size effect of the stressstrain responses can be captured accurately.
2) The dual lattice topology enables the seamless coupling of the mechanical and flow behavior of the material at grainscale. 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.
Acknowledgements
The work of WL and GC was supported by the Center for Sustainable Engineering of Geological and Infrastructure Materials (SEGIM) and the Quest high performance computing facility at Northwestern University. The Work of XZ was supported with ES3 R&D resources. The work of WC and LF was supported by the Department of Energy’s Basic Energy Sciences program under grant DEAC5206NA25396/FWPLANL20171450.
References
 Akono and Kabir (2016) Akono AT, Kabir P (2016) Microscopic fracture characterization of gas shale via scratch testing. Mechanics Research Communications 78:86–92
 Ashari et al (2017) Ashari SE, Buscarnera G, Cusatis G (2017) A lattice discrete particle model for pressuredependent inelasticity in granular rocks. International Journal of Rock Mechanics and Mining Sciences 91:49–58
 Bažant (1984) Bažant ZP (1984) Size effect in blunt fracture: concrete, rock, metal. Journal of Engineering Mechanics 110(4):518–535
 Bažant and Jirásek (2017) Bažant ZP, Jirásek M (2017) Creep and Hygrothermal Effects in Concrete Structures. Springer
 Bažant and Planas (1997) Bažant ZP, Planas J (1997) Fracture and size effect in concrete and other quasibrittle materials. CRC press
 Belheine et al (2009) Belheine N, Plassiard JP, Donzé FV, Darve F, Seridi A (2009) Numerical simulation of drained triaxial test using 3d discrete element modeling. Computers and Geotechnics 36(12):320–331
 Biot (1941) Biot MA (1941) General theory of threedimensional consolidation. Journal of applied physics 12(2):155–164
 Bolander and Berton (2004) Bolander J, Berton S (2004) Simulation of shrinkage induced cracking in cement composite overlays. Cement and Concrete Composites 26(7):861–871
 Bolander Jr and Saito (1998) Bolander Jr J, Saito S (1998) Fracture analyses using spring networks with random geometry. Engineering Fracture Mechanics 61(56):569–591
 Carey et al (2015a) Carey JW, Lei Z, Rougier E, Mori H, Viswanathan H (2015a) Fracturepermeability behavior of shale. Journal of Unconventional Oil and Gas Resources 11:27–43
 Carey et al (2015b) Carey JW, Rougier E, Lei Z, Viswanathan H (2015b) Experimental investigation of fracturing of shale with water. In: 49th US Rock Mechanics/Geomechanics Symposium, American Rock Mechanics Association
 Catalano et al (2014) Catalano E, Chareyre B, Barthélemy E (2014) Porescale modeling of fluidparticles interaction and emerging poromechanical effects. International Journal for Numerical and Analytical Methods in Geomechanics 38(1):51–71
 Chau et al (2016) Chau VT, Bažant ZP, Su Y (2016) Growth model for large branched threedimensional hydraulic crack system in gas or oil shale. Phil Trans R Soc A 374(2078):20150,418
 Cho et al (2013) Cho Y, Ozkan E, Apaydin OG (2013) Pressuredependent naturalfracture permeability in shale and its effect on shalegas well production. SPE Reservoir Evaluation & Engineering 16(02):216–228
 Choo and Sun (2018) Choo J, Sun W (2018) Coupled phasefield and plasticity modeling of geological materials: From brittle fracture to ductile flow. Computer Methods in Applied Mechanics and Engineering 330:1–32
 Cundall and Strack (1979) Cundall PA, Strack OD (1979) A discrete numerical model for granular assemblies. geotechnique 29(1):47–65
 Cusatis and Schauffert (2009) Cusatis G, Schauffert EA (2009) Cohesive crack analysis of size effect. Engineering Fracture Mechanics 76(14):2163–2173
 Cusatis and Zhou (2013) Cusatis G, Zhou X (2013) Highorder microplane theory for quasibrittle materials with multiple characteristic lengths. Journal of Engineering Mechanics 140(7):04014,046
 Cusatis et al (2003a) Cusatis G, Bažant ZP, Cedolin L (2003a) Confinementshear lattice model for concrete damage in tension and compression: I. Theory. Journal of Engineering Mechanics 129(12):1439–1448
 Cusatis et al (2003b) Cusatis G, Bažant ZP, Cedolin L (2003b) Confinementshear lattice model for concrete damage in tension and compression: II. Computation and Validation. Journal of Engineering Mechanics 129(12):1449–1458
 Cusatis et al (2011a) Cusatis G, Mencarelli A, Pelessone D, Baylot J (2011a) Lattice discrete particle model (LDPM) for failure behavior of concrete. II: Calibration and validation. Cement and Concrete composites 33(9):891–905
 Cusatis et al (2011b) Cusatis G, Pelessone D, Mencarelli A (2011b) Lattice discrete particle model (LDPM) for failure behavior of concrete. I: Theory. Cement and Concrete Composites 33(9):881–890
 Cusatis et al (2017) Cusatis G, Rezakhani R, Schauffert EA (2017) Discontinuous cell method (dcm) for the simulation of cohesive fracture and fragmentation of continuous media. Engineering Fracture Mechanics 170:1–22
 Davy et al (2007) Davy CA, Skoczylas F, Barnichon JD, Lebon P (2007) Permeability of macrocracked argillite under confinement: gas and water testing. Physics and Chemistry of the Earth, Parts A/B/C 32(8):667–680
 Detournay and Cheng (1993) Detournay E, Cheng AHD (1993) Fundamentals of poroelasticity. Chapter 5 in Comprehensive Rock Engineering: Principles, Practice and Projects, II pp 113–171
 Detwiler (2010) Detwiler RL (2010) Permeability alteration due to mineral dissolution in partially saturated fractures. Journal of Geophysical Research: Solid Earth 115(B9)
 Di Luzio and Cusatis (2009) Di Luzio G, Cusatis G (2009) Hygrothermochemical modeling of highperformance concrete. II: Numerical implementation, calibration, and validation. Cement and Concrete composites 31(5):309–324
 Elkhoury et al (2013) Elkhoury JE, Ameli P, Detwiler RL (2013) Dissolution and deformation in fractured carbonates caused by flow of co 2rich brine under reservoir conditions. International Journal of Greenhouse Gas Control 16:S203–S215
 Falk et al (2015) Falk K, Coasne B, Pellenq R, Ulm FJ, Bocquet L (2015) Subcontinuum mass transport of condensed hydrocarbons in nanoporous media. Nature communications 6:6949
 Frash et al (2016) Frash LP, Carey JW, Lei Z, Rougier E, Ickes T, Viswanathan HS (2016) Highstress triaxial directshear fracturing of utica shale and in situ xray microtomography with permeability measurement. Journal of Geophysical Research: Solid Earth 121(7):5493–5508
 Frash et al (2017) Frash LP, Carey JW, Ickes T, Viswanathan HS (2017) Caprock integrity susceptibility to permeable fracture creation. International Journal of Greenhouse Gas Control 64:60–72
 Grassl and Bolander (2016) Grassl P, Bolander J (2016) Threedimensional network model for coupling of fracture and mass transport in quasibrittle geomaterials. Materials 9(9):782
 Grassl et al (2015) Grassl P, Fahy C, Gallipoli D, Wheeler SJ (2015) On a 2d hydromechanical lattice approach for modelling hydraulic fracture. Journal of the Mechanics and Physics of Solids 75:104–118
 Gutierrez et al (2000) Gutierrez M, Øino L, Nygård R (2000) Stressdependent permeability of a demineralised fracture in shale. Marine and Petroleum Geology 17(8):895–907
 Hyman et al (2016) Hyman J, JiménezMartínez J, Viswanathan H, Carey JW, Porter M, Rougier E, Karra S, Kang Q, Frash L, Chen L (2016) Understanding hydraulic fracturing: a multiscale problem. Phil Trans R Soc A 374(2078):20150,426
 Ilgen et al (2017) Ilgen AG, Heath JE, Akkutlu IY, Bryndzia LT, Cole DR, Kharaka YK, Kneafsey TJ, Milliken KL, PyrakNolte LJ, SuarezRivera R (2017) Shales at all scales: Exploring coupled processes in mudrocks. EarthScience Reviews
 Jin et al (2017a) Jin C, Salviato M, Li W, Cusatis G (2017a) Elastic microplane formulation for transversely isotropic materials. Journal of Applied Mechanics 84(1):011,001
 Jin and Zoback (2017) Jin L, Zoback M (2017) Fully coupled nonlinear fluid flow and poroelasticity in arbitrarily fractured porous media: A hybriddimensional computational model. Journal of Geophysical Research: Solid Earth 122(10):7626–7658
 Jin et al (2017b) Jin Z, Li W, Jin C, Hambleton J, Cusatis G (2017b) Elastic, strength, and fracture properties of marcellus shale. arXiv preprint arXiv:171202320
 Jobmann et al (2010) Jobmann M, Wilsnack T, Voigt HD (2010) Investigation of damageinduced permeability of opalinus clay. International Journal of Rock Mechanics and Mining Sciences 47(2):279–285
 Kluge et al (2017) Kluge C, Blöcher G, Milsch H, Hofmann H, Nicolas A, Li Z, Fortin J (2017) Sustainability of fractured rock permeability under varying pressure. In: Poromechanics VI, pp 1192–1199
 Kozicki and Donzé (2008) Kozicki J, Donzé F (2008) A new opensource software developed for numerical simulations using discrete modeling methods. Computer Methods in Applied Mechanics and Engineering 197(4950):4429–4443
 Li et al (2017a) Li C, Caner FC, Chau VT, Bažant ZP (2017a) Spherocylindrical microplane constitutive model for shale and other anisotropic rocks. Journal of the Mechanics and Physics of Solids 103:155–178
 Li et al (2016) Li W, Jin C, Cusatis G (2016) Integrated experimental and computational characterization of shale at multiple length scales. In: New Frontiers in Oil and Gas Exploration, Springer, pp 389–434
 Li et al (2017b) Li W, Bousikhane F, Carey JW, Cusatis G (2017b) Computational analysis of the fracturepermeability behavior of shale. In: Poromechanics VI, pp 1200–1207
 Li et al (2017c) Li W, Bousikhane F, Carey JW, Cusatis G (2017c) Discrete modeling of the fracturepermeability behavior of shale. In: 51st US Rock Mechanics/Geomechanics Symposium, American Rock Mechanics Association
 Li et al (2017d) Li W, Jin Z, Cusatis G (2017d) Characterization of marcellus shale fracture properties through size effect tests and computations. arXiv preprint arXiv:171007366
 Li et al (2017e) Li W, Rezakhani R, Jin C, Zhou X, Cusatis G (2017e) A multiscale framework for the simulation of the anisotropic mechanical behavior of shale. International Journal for Numerical and Analytical Methods in Geomechanics 41(14):1494–1522
 Massey and WardSmith (1998) Massey BS, WardSmith J (1998) Mechanics of fluids, vol 1. CRC Press
 Miehe et al (2015) Miehe C, Mauthe S, Teichtmeister S (2015) Minimization principles for the coupled problem of darcy–biottype fluid transport in porous media linked to phase field modeling of fracture. Journal of the Mechanics and Physics of Solids 82:186–217
 Mikelic et al (2015) Mikelic A, Wheeler MF, Wick T (2015) A phasefield method for propagating fluidfilled fractures coupled to a surrounding porous medium. Multiscale Modeling & Simulation 13(1):367–398
 Mohammadnejad and Khoei (2013) Mohammadnejad T, Khoei A (2013) An extended finite element method for hydraulic fracture propagation in deformable porous media with the cohesive crack model. Finite Elements in Analysis and Design 73:77–95
 Nikolic et al (2016) Nikolic M, Ibrahimbegovic A, Miscevic P (2016) Discrete element model for the analysis of fluidsaturated fractured poroplastic medium based on sharp crack representation with embedded strong discontinuities. Computer Methods in Applied Mechanics and Engineering 298:407–427
 Noiriel et al (2010) Noiriel C, Renard F, Doan ML, Gratier JP (2010) Intense fracturing and fracture sealing induced by mineral growth in porous rocks. Chemical Geology 269(3):197–209
 Nygård et al (2006) Nygård R, Gutierrez M, Bratli RK, Høeg K (2006) Brittle–ductile transition, shear failure and leakage in shales and mudrocks. Marine and Petroleum Geology 23(2):201–212
 Obliger et al (2018) Obliger A, Ulm FJ, Pellenq RJM (2018) Impact of nanoporosity on hydrocarbon transport in shalesâ organic matter. Nano letters
 Pelessone (2006) Pelessone D (2006) MARS, modeling and analysis of the response of structures – userââs manual. ES3, Inc, San Diego, CA, USA
 Potyondy and Cundall (2004) Potyondy D, Cundall P (2004) A bondedparticle model for rock. International journal of rock mechanics and mining sciences 41(8):1329–1364
 Rezakhani and Cusatis (2016) Rezakhani R, Cusatis G (2016) Asymptotic expansion homogenization of discrete finescale models with rotational degrees of freedom for the simulation of quasibrittle materials. Journal of the Mechanics and Physics of Solids 88:320–345
 Rice and Cleary (1976) Rice JR, Cleary MP (1976) Some basic stress diffusion solutions for fluidsaturated elastic porous media with compressible constituents. Reviews of Geophysics 14(2):227–241
 Robinson et al (2014) Robinson M, Ramaioli M, Luding S (2014) Fluid–particle flow simulations using twowaycoupled mesoscale sph–dem and validation. International journal of multiphase flow 59:121–134
 Rutqvist et al (2008) Rutqvist J, Freifeld B, Min KB, Elsworth D, Tsang Y (2008) Analysis of thermally induced changes in fractured rock permeability during 8 years of heating and cooling at the yucca mountain drift scale test. International Journal of Rock Mechanics and Mining Sciences 45(8):1373–1389
 Salimzadeh and Khalili (2015) Salimzadeh S, Khalili N (2015) A threephase xfem model for hydraulic fracturing with cohesive crack propagation. Computers and Geotechnics 69:82–92
 Smith et al (2014) Smith J, Cusatis G, Pelessone D, Landis E, O’Daniel J, Baylot J (2014) Discrete modeling of ultrahighperformance concrete with application to projectile penetration. International Journal of Impact Engineering 65:13–32
 Soeder (1988) Soeder DJ (1988) Porosity and permeability of eastern devonian gas shale. SPE Formation Evaluation 3(01):116–124
 Su et al (2017) Su K, Sanz Perl Y, Onaisi A, Pourpark H, VidalGilbert S (2017) Experimental study of hydromechanical behavior of fracture of vaca muerta gas shale. In: 51st US Rock Mechanics/Geomechanics Symposium, American Rock Mechanics Association
 SuarezRivera and Fjær (2013) SuarezRivera R, Fjær E (2013) Evaluating the poroelastic effect on anisotropic, organicrich, mudstone systems. Rock mechanics and rock engineering 46(3):569–580
 Ulven and Sun (2017) Ulven OI, Sun W (2017) Capturing the twoway hydromechanical coupling effect on fluiddriven fracture in a dualgraph lattice beam model. International Journal for Numerical and Analytical Methods in Geomechanics
 Wang (2000) Wang HF (2000) Theory of linear poroelasticity with applications to geomechanics and hydrogeology. Princeton University Press
 Wang and Sun (2016) Wang K, Sun W (2016) A semiimplicit discretecontinuum coupling method for porous media based on the effective stress principle at finite strain. Computer Methods in Applied Mechanics and Engineering 304:546–583
 Witherspoon et al (1980) Witherspoon PA, Wang JS, Iwai K, Gale JE (1980) Validity of cubic law for fluid flow in a deformable rock fracture. Water resources research 16(6):1016–1024
 Yeo et al (1998) Yeo Id, De Freitas M, Zimmerman R (1998) Effect of shear displacement on the aperture and permeability of a rock fracture. International Journal of Rock Mechanics and Mining Sciences 35(8):1051–1070
 Zhang et al (2013) Zhang J, Kamenov A, Zhu D, Hill A (2013) Laboratory measurement of hydraulic fracture conductivities in the barnett shale. In: IPTC 2013: International Petroleum Technology Conference
 Zhou et al (2016) Zhou T, Zhang S, Feng Y, Shuai Y, Zou Y, Li N (2016) Experimental study of permeability characteristics for the cemented natural fractures of the shale gas formation. Journal of Natural Gas Science and Engineering 29:345–354
 Zimmerman and Bodvarsson (1996) Zimmerman RW, Bodvarsson GS (1996) Hydraulic conductivity of rock fractures. Transport in porous media 23(1):1–30