Gas Flow Tightly Coupled to Elastoplastic Geomechanics for Tight and Shale Gas Reservoirs: Material Failure and Enhanced Permeability

Unconventional natural gas such as tight and shale gas has become an increasingly important source of natural gas in the United States over the past decade and the estimates of the natural gas resource potential for shale gas range from 14.16×10¹² m³ to 28.3×10¹² m³ (500~1000 Tcf) (Arthur et al. 2008; Jenkins and Boyer 2008). Geological formations of tight and shale gas reservoirs exhibit extremely low permeability within which gas is trapped (Hill and Nelson 2000), and accordingly hydraulic fracturing is performed in order to enhance permeability and increase production rate (Dean and Schmidt 2008; Ji et al. 2009; Nassir et al. 2012). The success of gas production from the Barnett Shale is based on horizontal wells and hydraulic fracturing techniques and has led to gas production of other shale reservoirs such as the Marcellus Shale, one of the largest natural gas resources in the United States (Arthur et al. 2008; Cipolla et al. 2010).

Several studies have been made on production of gas from tight and shale gas reservoirs as well as hydraulic fracturing. Freeman et al. (2011) investigated non-Darcy flow without geomechanics, examining the effects of Knudsen diffusion on gas composition, because flow mechanism of shale gas is significantly different from that of conventional oil and gas reservoirs due to extremely low permeability and small pore throat. Vermylen and Zoback (2011) studied different scenarios for hydraulic fracturing along horizontal wells, and found significant differences in stimulation robustness for different fracturing procedures, for example, two alternatively fractured (zipperfrac) and simultaneously fractured (simulfrac) wells. Ji et al. (2009) developed a numerical model for hydraulic fracturing, accounting for coupled flow and geomechanics, and Nassir et al. (2012) incorporated shear failure to tensile fracturing. Dean and Schmidt (2008) employed the same fracturing algorithm, using different criteria of fracture propagation. Fisher and Warpinski (2011) analyzed the fracture growth induced by hydraulic fracturing with real data, and concluded that the fracture propagation was limited in the vertical direction, compared with the horizontal direction. Yet, coupled flow and geomechanics that can consider dynamic interrelations among pore volume, permeability, and material failure during production are little investigated. Geomechanics with material failure may change permeability and porosity followed by flow patterns significantly, and, in turn, redistribution of pore-pressure can affect geomechanics (e.g., Armero (1999) and Rutqvist and Stephansson (2003)), Thus, rigorous modeling for coupled flow and geomechanics is necessary for gas production in shale and tight gas reservoirs.

In this study, we perform numerical simulations of coupled flow and geomechanics in shale and tight gas reservoirs, focusing on geological failure and secondary fracturing of intact rock, which can increase permeability. To this end, we consider variations in not only permeability but also pore volume, which lead to dynamic permeability and porosity of the reservoirs, respectively. In particular, compaction in undrained condition such as low permeability induces the increase of pore-pressure, which can change effective stress fields significantly and may result in geological failure (Kim et al., 2012a).

Here, we employ the Mohr-Coulomb failure model for elastoplasticity, which is widely used to model failure in cohesive frictional materials. We also use dynamic permeability, which is a function of the failure status and geomechanical variables. Single-porosity and double-porosity models are both investigated. One may introduce double porosity and permeability in flow for better accurate modeling of flow because fractures occupy a much smaller fraction of bulk volume of a gridblock than rock matrix (Barenblatt et al.1960; Pruess and Narasimhan 1985; Bagheri and Settari 2008). On the other hand, it is required to use a slightly different geomechanics, because return mapping for elastoplasticity that models rock failure is performed at the gridblock level, not subelements. Note that the geomechanics still accounts for both geomechanical moduli of the fracture and intact rock. Then, in the case of the double porosity model in flow and geomechanics, we modify the constitutive relations proposed by Kim et al. (2012b), which generalized the constitutive model by Berryman (2002) to non-isothermal multiple porosity coupled flow and geomechanics systems.

For numerical simulation, we use a sequential implicit method, employing the fixed-stress split, which can provide unconditional stability and high accuracy (Settari and Mourits 1998; Kim et al. 2011). Specifically, flow is solved first, fixing the total stress fields, and then geomechanics is solved from the variables obtained at the previous flow step. We employ finite volume and finite element methods for flow and geomechanics in space discretization, respectively, and the backward Euler method in time discretization.

We perform numerical simulations for given original horizontal and vertical primary single fractures, under horizontal well production scenarios. During production, shear failure occurs not only near the tips of the primary fractures but also away from the fractures, generating secondary fractures, because production of gas causes redistribution of the effective stress fields, increasing the effective shear stress and resulting in plasticity. Secondary fractures increase permeability dramatically, which is a welcome development as it enhances gas production. The risk of potential detrimental environmental effects (e.g., through creation of a conduit to overlying aquifers) appears low because of the limited extent of the secondary fractures, but a definitive answer will require more focused study of this problem. This physics cannot be captured accurately by uncoupled or decoupled simulation, and thus tightly coupled flow and geomechanics are highly recommended for gas production in shale gas reservoirs.

We describe the governing equations for fluid and heat flow, as explained in Pruess et al. (1999). The governing equation for multiphase and multi-component flow comes from mass balance as

where the superscript k indicates the fluid component. d(·)/dt means the time derivative of a physical quantity (·) relative to the motion of the solid skeleton. m^k is mass of component k.f^k and q^k are its flux and source terms on the domain Ω with the boundary Γ, respectively, where n is the normal vector of the boundary.

The mass of the component k is written as

where the subscript J indicates fluid phases. φ is the true porosity, defined as the ratio of the pore volume to the bulk volume in the deformed configuration. S_J, ρ_J, and X_J^k are saturation and density of phase J, and the mass fraction of component k in phase J, respectively.

The mass flux term is obtained from

where w^k_J and f^k_J are the convective and diffusive mass flows of component k in phase J. For the liquid phase, J=L, w^k_J is given by Darcy's law as

where k is the absolute (intrinsic) permeability tensor. μ_J, k_rJ, p_J are the viscosity, relative permeability, and pressure of the fluid phase J, respectively. g is the gravity vector, and Grad is the gradient operator. For the gaseous phase, J = G,w^k_G is given by

where k_K is the Klinkenberg factor. The diffusive flow can be written as

where D^k_J and τ_G are the hydrodynamic dispersion tensor and gas tortuosity.

The governing equation for heat flow comes from heat balance, as

where the superscript H indicates the heat component. m^H, f^H, and q^H are heat, its flux, and source terms, respectively. m^H in the heat accumulation term is expressed as

where ρ_R, C_R, T, and T₀ are the density and heat capacity of the porous medium, temperature and reference temperature. The heat flux is written as

where K_H is the composite thermal conductivity tensor of the porous media. The specific internal energy, e_J, and enthalpy, h_J, of components k in phase J become, respectively,

The governing equations for geomechanics is based on the quasi-static assumption (Coussy 1995), written as

where Divis the divergence operator, σ is the total stress tensor, andp_b is the bulk density. The infinitesimal transformation is used to allow the strain tensor, ε, to be the symmetric gradient of the displacement vector, u,

We focus on single phase flow (i.e., gas flow) with elastoplastic geomechanics in this study.

We employ the constitutive relations for coupled multiphase non-isothermal flow and elastoplastic geomechanics, described in Coussy (1995). For elastoplastic mechanics and nonisothermal single phase flow, the constitutive relations are written as

where the subscript e indicates elasticity. σ', K_dr,C_e are the effective stress tensor, drained bulk modulus, and drained isothermal elasticity tensor at the level, respectively. 1 is the second order identity tensor. δζ =δm/ρ, where ρ is fluid density. α and M are the Biot coefficient and modulus in single phase fluid, written as

where K_s is the intrinsic solid grain bulk modulus, and c_fl is the fluid (gas) compressibility.

For plasticity in this study, we use the Drucker-Prager and Mohr-Coulomb models, which are widely used to model failure of cohesive frictional materials. Fig. 1 shows the comparison between the Drucker-Prager and Mohr-Coulomb models.

The Drucker-Prager model of this study is expresses as

where f and g are the yield and plastic potential functions, respectively. Ψ_f and Ψ_d are the friction and dilation angles, respectively. θ is the Lode angle, written as

I₁, J₂, and J₃ are invariants of the effective stresses as follows

where tr is the trace of a second order tensor.

The Mohr-Coulomb model is given as

where σ₁',σ₂',σ₃' are the maximum, intermediate, and minimum principal effective stresses. Note that the coefficients of Drucker-Prager model used in this study (i.e., Eqs. (18) and (19)) are adjusted to the Mohr-Coulomb model, and the two models have the same yield and potential functions. However, since the Drucker-Prager model is a I₁-J₂ plasticity model, the coefficients of the Drucker-Prager model in Eq. (17) are calculated explicitly while the coefficients associated with J₃ in the Mohr-Coulomb model are calculated implicitly. The different treatments of J₃may cause slight differences between the two plastic models. We will show validation of the numerical implementation and compare the results between the two models. When material failure occurs, micro-fractures are created. As a result, permeability can increase dramatically and discontinuously. We consider the change in permeability, using a permeability (or transmissibility) multiplier ω_ρ, as

by which the intrinsic permeability, k₀, is expressed as k₀ = (O_pk_m, where k_m is the intact rock permeability. Thus, the permeability is a step function for a failure status. In this study, we simply use ω_ρ = 10^-2 when the intact rock faces failure. Other complicated models for ω_ρ may be used for high accuracy of certain specific reservoirs.

We use the finite volume and finite element methods for flow and geomechanics in space discretization, respectively, which are widely used in reservoir and geotechnical engineerings, respectively (Aziz and Settari 1979; Lewis and Schrefler 1998; Hughes 2000). This mixed space discretization can provide advantages such as local mass conservation and better numerical stability in space, compared with the finite element method for both flow and geomechanics (Jha and Juanes 2007; Kim et al. 2011).

For time discretization, we use the backward Euler method, which is typically used in reservoir simulation. We use TOUGH+ family codes as a fluid and heat flow simulator (Pruess et al., 1999) and ROCMECH for a geomechanics simulator, recently developed in the Lawrence Berkeley National Laboratory (Kim et al., 2012a, b).

When we solve coupled problems numerically, there are two approaches; fully coupled (monolithic) and sequential methods. The fully coupled methods solve the coupled problems simultaneously, using the Newton-Raphson method (Lewis and Sukirman, 1993a,b; Sukirman and Lewis, 1993; Pao et al., 2001; Gutierrez and Lewis, 2002; Pao and Lewis, 2002; Lewis et al., 2003; White and Borja, 2008). This approach typically provides unconditional stability, but requires a unified flow and geomechanics, which results in the enormous software development and huge computational cost. On the other hand, the sequential methods can offer the use of existing robust flow and geomechanics simulators only by constructing an interface between them (Armero and Simo, 1992, 1993; Settari and Mourits, 1998; Armero 1999). This provides wide flexibility and high efficiency for code management, but may be limited by numerical stability and accuracy, which are the main issues in using the sequential methods. According to Kim et al. (2011b), among several sequential methods, the fixed stress sequential method can provide unconditional stability and high accuracy, competitive with the fully coupled methods. The fixed stress method solves the flow problem, fixing the total stress field, where the strain and displacement fields can be changed. Then, it solves geomechanics, based on the solutions obtained at the previous flow problem. This sequential method can easily be implemented by the Lagrange porosity function Φ and its correction ΔΦ (Settari and Mourits, 1998; Kim et al., 2012a), as follows.

where c_p is the pore compressibility in conventional reservoir simulation (Aziz and Settari, 1979), and σ_v is the total volumetric mean stress. Φ is defined as the ratio of the pore volume in the deformed configuration to the bulk volume in the reference (initial) configuration. The porosity correction term, ΔΦ, is calculated from geomechanics, which corrects the porosity estimated from the pore compressibility.

When implementing the Drucker-Prager and Mohr-Coulomb models, we follow the algorithms proposed by Wang et al. (2004). For the Mohr-Coulomb model, we may face numerical instability at the discontinuous corners, in which case we use the Drucker-Prager model at the corners (Borja et al., 2003; Wang et al., 2004).

We perform two test cases for 2D gas flow simulation coupled to the plane strain geomechanics, which are aimed at gas production from tight gas or shale gas reservoirs with several uniform long horizontal wells. This production scenario can reduce a 3D system to 2D because of symmetry and the plane strain condition, as shown in Fig. 2. Cases 1 and 2 contain horizontal and vertical fractures, respectively.

For Case 1, the 2D reservoir domain is divided into 52×96 gridblocks in horizontal and vertical (x, z) directions. Horizontally, we have a uniform grid size, Δx = 1.0m. On the other hand, in the vertical direction, we have a uniform grid size that has Az = 1.0m from the first to the 40th horizontal grid layers, non uniform sizes that are 0.64m, 0.32m, 0.16m, 0.08m, 0.04m, 0.02m from the 41th to 46th horizontal grid layers, a horizontal fracture whose aperture is 0.01m at the 47th horizontal grid layer, non uniform sizes that are 0.02m, 0.04m, 0.08m, 0.16m, 0.32m, 0.64m from 48th and 53th horizontal grid layers, and a uniform grid size that has Az = 1.0m from the 54th to the last horizontal grid layers. The length of the horizontal fracture, where the production well is located at (x=0.5m, z=41.265m), is 10m, and the aperture is 0.01m.

For Case 2, the reservoir domain is divided into 56×80 gridblocks in horizontal and vertical (x, z) directions. Vertically, we have a uniform grid size, Δz = 1.0m. On the other hand, in the horizontal direction, we have a 0.01m grid size of first vertical grid layer that contains a vertical fracture whose aperture is 0.01m, non uniform sizes that are 0.02m, 0.04m, 0.08m, 0.16m, 0.32m, 0.64m from second and 7th vertical grid layer, and a uniform grid size that has Ax = 1.0m from the 8th to the last vertical grid layers. The length of the vertical fracture, where the production well is located at (x=0.005m, z=40.5m) is 20m.

For both cases we have the same flow and geomechanics conditions. For flow, we have no flow at the boundaries. The initial pressure is 68.95MPa with no gravity. The initial temperature is 176.67 oC. The heat capacities of the media, and the wet and dry thermal conductivities for all layers are 1000 Jkg^-1° C^-1, 3.1Wm^-1° C^-1, and 0.5Wm^-1° C^-1, respectively. We apply a constant pressure at the production well, 30MPa. The initial permeabilities for the fracture and rock matrix are 2.556 × 10¹³m² and 8.645 × 10¹⁸m², respectively, where 1 Darcy is 9.87 × 10¹³m². The initial porosities for them are 1.0 and 0.3, respectively. For geomechanics, we have no horizontal displacement at both sides, based on symmetry shown in Fig. 2, and no displacement at the bottom, the overburden at the top, σ = 68.95MPa. The initial vertical and horizontal principal stresses are -68.95MPa. Young's moduli for the fracture and rock matrix are 185MPa and 1.282GPa, respectively, and Poisson's ratios for them are 0.0 and 0.22, respectively. For both materials, the Biot coefficients, a, are 1.0, the bulk densities are 2600 kg·m^-3, and the thermal dilation coefficients are 1.0×10^–5°C^–1m. The cohesion, fraction and dilation angles for plasticity are 4.0MPa, Ψ_f = 28.60and Ψ_d = 28.60 for the both materials. We use a generalized reservoir model, rather than selecting geomechanical properties for a certain specific reservoir. The geomechanical properties used in this study are within a range of the properties of shale gas reservoirs (Eseme et al. 2007). We use low Young's modulus for the primary fracture because the fracture is typically much deformable than intact rock (e.g., Nguyen and Abousleiman (2010)).

We first perform numerical simulations for Cases 1 and 2 on the basis of elasticity. The simulations can indicate the potential areas of failure during gas production. Let us introduce λ in order to indicate the distance between the effective stresses and yield function, for an example of the Mohr Coulomb model,

High values indicate the potential areas of shear failure, when using the Mohr Coulomb model. In addition, failure can occur when λ > 0. From numerical simulation, Fig. 3 shows distributions of λ for Cases 1 and 2 when using elasticity. λ is calculated from the effective stresses obtained from elasticity. The shear stresses are concentrated at the right boundaries because of the horizontally constrained boundary condition, especially around the top and center of the right side boundary, resulting in high values of λ. This implies that gas production can cause failure away from the production well. This failure may be beneficial to increase permeability in tight gas or shale gas reservoirs, but it could be problematic if surface facilities exist, or weak geological sealing occurs, in the vicinity of the failure area. In this study, we focus on productivity of the gas reservoirs, investigating changes in permeability, failure, and flow and geomechanical variables.

We perform failure analyses for Cases 1 and 2, using the Drucker-Prager model adjusted for the Mohr-Coulomb model. For Case 1, failure occurs at the right boundary and it moves into the domain, shown in Fig. 4. The areas that have non-zero in the Fig. experienced plasticity during simulation. Specifically, the value indicates the number of Gauss points within a gridblock which faced plasticity. Most of the non-zero areas returned to elasticity, after facing plasticity, but they yielded the enhanced permeability because of failure. In this simulation, the areas with the enhanced permeability are not yet connected to the original horizontal fracture. However, considering the failure propagation, they may possibly be connected when further simulation proceeds.

We introduce three monitoring points for Case 1; P1, P2, and P3. P1, P2, and P3 are the points at (x=4.5m, z=-41.265m), (x=39.5m, z=-41.265m), and (x=49.5m, z=-41.265m), respectively. P1 is located within the fracture, while P2 and P3 are far from the fracture. In Fig. 5, P1, P2, and P3 enter plasticity at early time, then they re-enter elasticity. The effective stresses at P1 reach plasticity almost immediately, and propagate along the yield function. Then, they return to the elastic region. The effective stresses at P2 and P3 propagate vertically, having high values of shear effective stresses. Once reaching plasticity, they return to elasticity after short time. This behavior mainly results from enhancement of permeability, followed by redistribution in fluid pressure and effective stress. The Lode angle at P1 varies during simulation, which leads to softening at late times. We also observe that the pressures at P2 and P3 become significantly closer because of the enhanced permeability. Fig. 6 shows the effective stresses on the σ_m'-τ_m' plot in order to test the algorithm of the Drucker-Prager model adjusted for the Mohr-Coulomb model in ROCMECH (i.e., explicit treatment of J₃). In the Fig., the adjusted Drucker-Prager model properly simulates the failure of the Mohr-Coulomb model. From the results of Figs. 4 and 5, we find significant nonlinearity and complexity when using stress dependent permeability in conjunction with plasticity, changing permeability fields dynamically. This behavior cannot properly be simulated by the conventional flow-only or decoupled simulation.

For Case 2, we obtain similar results to Case 1. Fig. 7 shows the propagation of the region where permeability increases due to plasticity. The region propagates from the middle of the right boundary to the center of the domain, similar to Case 1. Additionally, failure begins to occur at the right top corner, and there are a few gridblocks near the original fracture tips, although they are invisible in Fig. 7 because of the tiny scale.

The three monitoring points for Case 2 are introduced; P4, P5, P6. P4, P5, and P6 are the morning points at (x=0.005m, z=-29.5m), (x=23.77m, z=-40.5m), and (x=43.77m, z=-40.5m), respectively. P4 is located right near the fracture tip, while P5 and P6 are far from the fracture. The effective stresses at P4, P5 and P6 propagate toward the yield function, entering plasticity and then returning to elasticity, similar to those for Case 1, as shown in Fig. 8. We observe non-smooth evolution in effective stress and pressure at P5. When failure occurs at P5, the permeability at P5 increases, connected to the enhanced permeability region. This can redistribute the pressure fields over the domain, followed by the effective stresses. Compared with Fig. 5, the evolution of effective stress at P5, shown in Fig. 8, is highly discontinuous. Discontinuity is also observed in evolution of pressure. Fluid flow for the enhanced permeability region, such as P6, can support pressure at P5 significantly from influx from P6 to P5, changing the effective stresses considerably. We plot σ_m' and τ_m' in Fig. 9 in order to test validity for the algorithm of the Drucker-Prager model for the Mohr-Coulomb model, again, and the Fig. still shows good validity for modeling the failure of the Mohr-Coulomb model.

We modify the friction and dilation angles, ψ_f= 10.00 and ψ_d =10.00, and apply a constant production rate of Q_p = 4 × 10^-3kg/s in order to account for various plastic cases. We use the different algorithm for the Mohr-Coulomb model with the implicit treatment of J₃ in ROCMECH. The initial mass of gas in place is M _ini = 3.105 x10⁵kg. The left of Fig. 10 shows the area of the enhanced permeability due to plasticity after simulation, at t_d = 0.04, where the dimensionless time for a constant production rate is defined as t_d = Q_pt/M_ini. In contrast with the previous results, failure occurs near the original horizontal fracture. Since the Mohr-Coulomb model is not significantly different from the Drucker-Prager model used in this study, the differences stem from the different friction and dilation angles.

The right of Fig. 10 shows evolution of effective stresses at P7 and P8, which are located at (x=10.5m, z=-41.265m) and (x=19.5m, z=-41.265m). P7 is near the original horizontal fracture tip, while P8 is away from the fracture tip. The region near the fracture faces failure during production because the friction angle is lower than the previous case. P8 has not yet faced plasticity, but it will possibly reach the yield function in the light of the trend of its propagation.

For Case 2, we use the friction and dilation angles, ψ =10.00 and ψ_d =10.00, the Mohr-Coulomb model, and a constant production rate of Q =4×10³kg/s, same as the previous case, while the total initial mass of gas in place is M_ini = 2.81 × 10⁵kg. The region of the enhanced permeability occurs near the original vertical fracture, which is similar to Case 1. From the left of Fig. 11, failure propagates into the reservoir slowly. We also find failure to propagate diagonally from the fracture tips. We introduce another monitoring point, P9, located horizontally next to center of the fracture (x=0.02m, z=-40.5m). The effective stresses at P4 and P9 propagate to the yield function, and enter the plastic regime. After plasticity, both follow the yield surface, and return to elasticity at late time. When the Lode angle approaches 0° or 60°, The Mohr Coulomb plasticity is switched into the Drucker-Prager model to avoid numerical instability at the edges of the Mohr Coulomb model. From the change in the plastic model, the effective stresses at P4 can propagate slightly over the Mohr Coulomb failure line, as shown in Figs. 11 and 12.

As failure proceeds during simulation, the permeability field changes discontinuously because permeability is a discontinuous function of the failure status. From this discontinuity, the responses in pressure are also discontinuous (Fig. 12). When the permeability of the regions away from P4 or P9 becomes enhanced by plastic failure, the flow of gas into P4 or P9 increases significantly, supporting pressures at the points suddenly. Thus, the pressures at P4 and P9 increase and decrease repeatedly.

Failure implies creation of small fractures, and the theory of single porosity has significant limitation not only in coupled flow and geomechanics also in flow-only simulations because the single porosity cannot properly represent two distinct materials such as the fracture and intact rock (Bai 1999; Barenblatt et. al. 1960; Pruess and Narasimhan 1985). To solve this issue, the double porosity (, or more generally multiple porosity model) has been introduced for more realistic simulation, representing local heterogeneity (Barenblatt et al.1960; Pruess and Narasimhan 1985; Berryman 2002; Bagheri and Settari 2008). The fracture-rock matrix systems based on the double porosity model are representative of any composite systems that consist of a high permeable material transporting fluid over the domain and the other materials storing fluid and conveying it to the high permeable material. Then, the double porosity model can be well suited for flow for fractured reservoirs.

We apply the constitutive relations for the multiple porosity model in thermoporomechanics proposed by Kim et al. (2012b), extended from the relations by Berryman (2002). We do not use the Einstein notation in this section for clarity. Then, the constitutive relations for the double porosity model such as a fracture-rock matrix system are written as

where the subscript l indicates a material (sub-element) within a gridblock. K_dr and C_up are the upscaled elastoplastic drained bulk and tangent moduli at the level of the gridbolck, respectively. $bl*$ and $b˜l*$ are the coupling coefficients. b_l, $b˜l$⁠, K_dr, $Ll,m−1$ are written as

where a_l, α_{T, l}, η_l, and K_l are the Biot coefficient, thermal dilation coefficient, volume fraction, drained bulk modulus for material l, respectively. L_{l, m}(≡L) represents the Biot modulus matrix of the double porosity model (e.g., fracture-rock matrix system), where N_f and N_M are the inverse of the Biot moduli, M_f and M_M, for the fracture and rock matrix media, respectively, (i.e., N_f =1/M_f and N_M =1/M_M). The subscripts f and M indicate the fracture and rock matrix.

For naturally fractured reservoirs, the double porosity model is used initially, while, in this study, we change the single porosity model into the double porosity during simulation when a material faces plasticity. Thus, for the naturally fractured reservoirs, C_up and K_dr at a gridblock are obtained from an upscaling from given properties of fracture and rock matrix materials. Accordingly, the return mapping for elastoplasticity is performed at all the subelements (Kim et al., 2012b). On the other hand, in this study, C_up and K_dr are directly obtained from the elastoplastic tangent moduli at a gridblock (global) level, not the subelements, while we need to determine the drained bulk moduli of fracture and rock matrix materials for the double porosity model, followed by the coupling coefficients. To this end, we assume that the rock matrix has the same drained bulk modulus as that of the single porosity material before plasticity (i.e., elasticity), because the rock matrix is undamaged. Then, from Eq. 29₃, the drained bulk modulus of the fracture can be determined as

Considering K_dr and Ä^to be positive for wellposedness, the volume fraction of the fracture, η_f, has the constraint as

Then, according to Kim et al. (2012b), Lagrange's porosity for material l can be written as

from which, the modified fixed stress split yields

where σ_v is the total volumetric mean stress at the gridblock. The permeability (or transmissibility) multiplier for the fracture $ω˜p$ in the double porosity system can be given as $ω˜p=ωp/ηf$⁠, based on conservation of fluid fluxes between the single and double porosity models.

We perform numerical simulations for Cases 1 and 2 with the Drucker-Prager model, using the aforementioned modified double porosity model when the material faces failure. In the case of the double porosity model, the volume fractions of the secondary fracture from failure and rock matrix are 0.1 and 0.9, respectively. The volume fraction of the fracture satisfies the constraint in Eq. (31) for the test cases.

Fig. 13 shows evolution of the enhanced permeability area at different times when using the double porosity model. The propagation of the enhanced permeability area is faster than that of the single porosity model, because the permeability at the secondary fractures from failure is higher than that of the single porosity model. Due to the fast propagation of failure, the enhanced permeability area is connected to the original horizontal fracture for the given simulation time. After they are connected, noticeable further failure does not occur.

Fig. 14 shows evolutions of effective stress and pressure at P1, P2, and P3. The original horizontal fracture where P1 is located faces failure almost at initial time, and then it returns to elasticity. At P2 and P3, away from the horizontal fracture, the effective stresses propagate vertically at early times until they reach the yield function. The vertical propagation indicates the significant increase of shear stress. After plasticity, the effective stresses at P2 and P3 re-enter elasticity.

At early times, pressure at P2 decreases lower than that at P3 because P2 is closer to the production well than P3. When P2 faces failure and has the secondary fracture with high permeability, the pressures at P2 and P3 become almost identical because of the enhanced permeabilities at P2 and P3, causing significantly low pressure gradient between the two points. The behavior of the effective stresses and pressure is the same as that of the single porosity model, even though the time scale and enhanced permeability area are different.

For Case 2 with the double porosity model, we find from Fig. 15 faster evolution of the enhanced permeability zone than the case with the single porosity model, which is the same observation as that of Case 1. But, unlike Case 1 with the double porosity model, in Fig. 15, the enhanced permeability is not connected to the original vertical fracture.

From Fig.16, the effective stresses at P4, P5, and P6 enter plasticity, and then they re-enter elasticity, which is similar to the behavior of Case 1 with the double porosity model. The pressures at P5 and P6 are identical because the two points face failure at the same time, yielding the high permeability area between them.

We investigated coupled flow and geomechanics in gas production from the extremely low permeability reservoirs such as tight and shale gas reservoirs. We accounted for dynamic changes in pore volume and permeability in the context of coupled flow and geomechanics. The fixed-stress sequential implicit method was adopted for simulation of coupled flow and geomechanics. During simulation, permeability is updated every time step, which is based on the status of material failure.

We performed numerical tests for two cases: the horizontal and vertical fractures, considering the single and double porosity models. The Drucker-Prager and Mohr-Coulomb models were employed and tested for cohesive frictional materials. For the single porosity model, we found that pressure drop at the well caused failure not only near the fracture but also away from the well, increasing shear stresses significantly. Even though the pressure support from reservoir compaction is small because of low coupling strength in pore volume from compressible fluid (i.e., gas), changes in permeability from material failure are significant, altering flow regimes dynamically. The changes in flow in turn affected geomechanics, causing further failure, followed by the propagation of the enhanced permeability area. We also observed that the enhanced permeability area highly depends on plastic models and their parameters (e.g., cohesion, friction and dilation angles).

We introduced the double porosity model for secondary fracturing with slight modification from the multiple porosity constitutive relations. From numerical simulation, we obtained different results from the single porosity model because the permeability of the secondary fracture is higher than the permeability of the single continuum. The different areas of the enhanced permeability were found, and the faster propagation of the areas was also observed.

In conclusion, from the numerical tests for tight and shale gas reservoirs, the physics in highly stress-sensitive reservoirs may not be captured accurately by flow only simulation or flow simulation decoupled from geomechanics, and thus tightly coupled flow and geomechanics in porosity and permeability are strongly recommended for gas production in stress-sensitive tight and shale gas reservoirs.

This study was supported by the US Environmental Protection Agency, Office of Water, under an Interagency Agreement with the U.S. Department of Energy at the Lawrence Berkeley National Laboratory through Contract No. DE-AC02-05CH11231, and by RPSEA (Contract No. 08122-45) through the Ultra-Deepwater and Unconventional Natural Gas and Other Petroleum Resources Research and Development Program as authorized by the US Energy Policy Act (EPAct) of 2005. The research described in this article has been funded wholly (or in part) by the U.S. Environmental Protection Agency through Interagency Agreement (DW-89-92235901-C) to the Lawrence Berkeley National Laboratory. The views expressed in this article are those of the author(s) and do not necessarily reflect the views or policies of the EPA.