A Study of Nonlinear Elasticity Effects on Permeability of Stress Sensitive Shale Rocks Using an Improved Coupled Flow and Geomechanics Model : A Case Study of the Longmaxi Shale in China

Gas transport in shale gas reservoirs is largely affected by rock properties such as permeability. These properties are often sensitive to the in-situ stress state changes. Accurate modeling of shale gas transport in shale reservoir rocks considering the stress sensitive effects on rock petrophysical properties is important for successful shale gas extraction. Nonlinear elasticity in stress sensitive reservoir rocks depicts the nonlinear stress-strain relationship, yet it is not thoroughly studied in previous reservoir modeling works. In this study, an improved coupled flow and geomechanics model that considers nonlinear elasticity is proposed. The model is based on finite element methods, and the nonlinear elasticity in the model is validated with experimental data on shale samples selected from the Longmaxi Formation in Sichuan Basin China. Numerical results indicate that, in stress sensitive shale rocks, nonlinear elasticity affects shale permeability, shale porosity, and distributions of effective stress and pore pressure. Elastic modulus change is dependent on not only in-situ stress state but also stress history path. Without considering nonlinear elasticity, the modeling of shale rock permeability in Longmaxi Formation can overestimate permeability values by 1.6 to 53 times.


Introduction
Shale gas has become an important source of natural gas.Successful shale gas extraction requires accurate and detailed characterization of shale reservoir properties.Permeability, porosity, and elastic properties of shale reservoir rocks are critical parameters that affect the fluid flow and rock deformation in shale reservoirs.These parameters are stress sensitive, as their values can vary significantly with in-situ stress changes and consequently affect the fluid flow and geomechanics effects in shale [1][2][3].Therefore, it is critical to apply adequate techniques to the study of effects of stress sensitivity in shale rocks.Specifically, experimental analysis and mathematical models are often used to study stress sensitive effects.
Experimental studies on shale cores provide correlations between rock property changes and stresses changes.Based on experiments, many correlations between permeability, porosity, rock elastic properties, pore structures, anisotropy, stress, and pore pressure were established, and some of the correlations were used to validate theoretically proposed mathematical functions that depict property changes in stress sensitive shale rocks [1][2][3][4][5][6][7][8][9][10][11][12][13].These experimental results help to establish validated correlations describing shale rock stress sensitivity, and help to improve the accuracy and reliability of theoretically proposed analytical models.However, in these studies, the analytical solutions usually ignore the volume of the analyzed sample and treat the sample as a point.These solutions assume that properties are uniformly distributed in the analyzed volume, which limits the studies' capability of addressing spatial characteristics of shale rock stress sensitivity.Besides, exact analytical solutions for nonlinear behaviors in fluid flow and geomechanics problems are not available, while nonlinearity (e.g., nonlinear elasticity and nonlinear shale gas flow behaviors) is an important effect that cannot be neglected in flow and geomechanics problems, as proved in many studies [7,[14][15][16].
Many studies based on numerical models were also carried out to study the stress sensitivity in reservoir rocks.As a complement to analytical approaches, numerical models are known for their capability of handling complex temporal and spatial evolutions of reservoir characteristics, as numerical methods can efficiently address the heterogeneity in space and the nonlinear behaviors of flow and geomechanics problems [17,18].In the study of stress sensitivity related problems, the involved numerical simulators usually implement the coupling between the fluid flow problem and the geomechanics problem.Wu et al. incorporated an analytical correlation into their numerical simulator coupling flow and geomechanics problems to model the stress distribution in a field-scale problem considering stress sensitive reservoir rock properties [17].This lab-validated analytical correlation represents the stress sensitive permeability hysteresis curve depicting permeability changes due to stress changes.Cao et al. also used a coupled flow and geomechanics model to study the effects of stress sensitivity on shale permeability [18].The stress sensitive permeability was incorporated in their model using a lab-validated correlation.An et al. used a coupled flow and geomechanics simulator to test the effect of stress-dependent permeability on hydrocarbon production in organic-rich shale reservoirs [19].They indicated that the selection of the correlation for the stress-dependent permeability significantly affects the field production performance.Shovkun and Espinoza applied an empirical correlation between horizontal permeability, compressibility, and stress to their coupled flow-geomechanics simulator for shale reservoir, and pointed out that the degree of the permeability alteration due to its stress sensitivity varies with location in the reservoir: near-well locations have stronger stress changes due to production while far-field locations have smaller stress changes due to production [20].While these researches provide insights in numerical investigations of stress sensitive reservoir rocks, limitations of these previous numerical studies are observed.In An et al. [19] and Shovkun and Espinoza [20], the numerical methods for the geomechanics modeling only consider linear elasticity, implying that these numerical methods could not describe the nonlinear elastic behaviors as rocks are being compacted and consolidated due to in-situ stresses changes and hydrocarbon depletion in the reservoir.The correlations between rock deformation and stress sensitive rock properties were theoretically derived and they lack the validation with field/experimental results, indicating that their methods work for theoretical analysis but have limitations when it comes to realistic shale rocks [19,20].
Based on the review of previous works, it is meaningful to combine the numerical modeling with lab-validated analytical correlations so that the stress sensitive shale rock properties can be properly addressed.In addition, it is important to take nonlinear elasticity into account to improve the effectiveness of the modeling study.Thus, the mechanism of gas production in realistic shale gas reservoirs can be better studied.
In this work, an improved coupled flow and geomechanics model considering lab-validated nonlinear elasticity is proposed in a case study.The lab validation is based on laboratory experiments on shale samples taken from the realistic shale gas play in the Lower Silurian Longmaxi Formation in Southwest China.Using this model, the characteristics of spatial distribution of stress and stress sensitive rock properties are analyzed.Also, quantification of the effect of considering nonlinear elasticity on numerical results is presented.In the paper, first, the target shale gas reservoir in the Lower Silurian Longmaxi Formation and its core samples are described.The stress sensitive permeability hysteresis curves obtained from lab measurements are also presented.Second, the finite element based numerical model that couples fluid flow and geomechanics is introduced, and nonlinear elasticity is incorporated in the numerical model based on elastic modulus alteration due to nonlinear rock consolidation.The incorporation of nonlinear elasticity is validated with the experimental results.Finally, using the calibrated numerical model, sensitivity analyses for stress, permeability, pore pressure, and porosity are conducted to improve the understanding of stress sensitivity in the shale gas play of the Lower Silurian Longmaxi Formation.Results from numerical modeling only considering linear elasticity are compared with numerical results considering nonlinear elasticity so that the effect of nonlinear elasticity is quantified.

Experimental Procedure
A set of gas permeability measurement experiments with varying confining pressure are conducted to obtain correlations between shale permeability and stress.Core samples were taken from the Longmaxi Formation in Sichuan Basin in Southwest China as this formation exhibits a great potential of shale gas production [21], and the effect of stress sensitivity on shale properties is significant in this formation [11,12].The shale sample description is in Table 1.Apparent permeability was measured with nitrogen flooding as in Figure 1.Nitrogen was used as it has similar properties to shale gas.It is also safe to use nitrogen.The constant inlet pressure (p in ) is 3 MPa and the constant outlet pressure (p out ) is 2 MPa.The initial confining stress is 5 MPa, and it was increased to 22 MPa in the stress loading process.It was then decreased back to 5 MPa in the stress unloading process.In stress loading and unloading, permeability values were measured at certain discrete confining stress states.Once the permeability was measured, the confining stress was then adjusted to move to the next discrete stress state.Since the lab was based on gas measurement, no liquid saturation was involved.Thus, the fragile rock samples did not experience high pressure which would damage the samples.sensitive rock properties are analyzed.Also, quantification of the effect of considering nonlinear elasticity on numerical results is presented.In the paper, first, the target shale gas reservoir in the Lower Silurian Longmaxi Formation and its core samples are described.The stress sensitive permeability hysteresis curves obtained from lab measurements are also presented.Second, the finite element based numerical model that couples fluid flow and geomechanics is introduced, and nonlinear elasticity is incorporated in the numerical model based on elastic modulus alteration due to nonlinear rock consolidation.The incorporation of nonlinear elasticity is validated with the experimental results.Finally, using the calibrated numerical model, sensitivity analyses for stress, permeability, pore pressure, and porosity are conducted to improve the understanding of stress sensitivity in the shale gas play of the Lower Silurian Longmaxi Formation.Results from numerical modeling only considering linear elasticity are compared with numerical results considering nonlinear elasticity so that the effect of nonlinear elasticity is quantified.

Experimental Procedure
A set of gas permeability measurement experiments with varying confining pressure are conducted to obtain correlations between shale permeability and stress.Core samples were taken from the Longmaxi Formation in Sichuan Basin in Southwest China as this formation exhibits a great potential of shale gas production [21], and the effect of stress sensitivity on shale properties is significant in this formation [11,12].The shale sample description is in Table 1.Apparent permeability was measured with nitrogen flooding as in Figure 1.Nitrogen was used as it has similar properties to shale gas.It is also safe to use nitrogen.The constant inlet pressure ( ) is 3 MPa and the constant outlet pressure ( ) is 2 MPa.The initial confining stress is 5 MPa, and it was increased to 22 MPa in the stress loading process.It was then decreased back to 5 MPa in the stress unloading process.In stress loading and unloading, permeability values were measured at certain discrete confining stress states.Once the permeability was measured, the confining stress was then adjusted to move to the next discrete stress state.Since the lab was based on gas measurement, no liquid saturation was involved.Thus, the fragile rock samples did not experience high pressure which would damage the samples.

Results and Discussion
Correlations between gas permeability and confining stress for all five samples are shown in Figure 2a-e.Figure 2f conceptually shows the processes of stress loading and unloading.In all

Results and Discussion
Correlations between gas permeability and confining stress for all five samples are shown in Figure 2a-e.Figure 2f conceptually shows the processes of stress loading and unloading.In all Energies 2018, 11, 329 4 of 16 samples, permeability rapidly decreases as stress is being loaded, indicating permeability is damaged by increased stress.Then, permeability slowly recovers as the confining stress is unloaded.However, it is noted that the stress unloading process is not able to fully recover the damaged permeability to its original value.by increased stress.Then, permeability slowly recovers as the confining stress is unloaded.However, it is noted that the stress unloading process is not able to fully recover the damaged permeability to its original value.
This irreversible permeability damage is the hysteresis phenomenon for stress-sensitive permeability.A conceptual hysteresis model of stress-sensitive permeability in the literature also depicted similar phenomenon [17].This hysteresis phenomenon is attributed to the nonlinear consolidation in the stress loading/unloading cycling [14].Sone and Zoback [14] indicated that, as stress is changed, the stress-strain relationship in rock deformation is not always linear, and deformation exhibits nonlinear elastic behaviors [14].The stress sensitive permeability hysteresis correlations obtained in this lab represent the nonlinear elastic rock deformation.It is used to improve the coupled flow and geomechanics model in this study.

Improved Finite Element Numerical Model
An improved numerical model considering nonlinear elasticity is presented to address the stress sensitive effects in shale rocks.The improvement is achieved by adding nonlinear elasticity to a commonly used poroelastic model as in [22].The model is validated with the aforementioned experimental data and can effectively address stress sensitive reservoir characteristics specifically in the Lower Silurian Longmaxi Formation.

Numerical Solution to the Coupled Flow and Geomechanics Problem
Accurate modeling of multiphysics problems in shale reservoirs is critical for shale gas extraction.The physics in shale reservoirs can be effectively depicted by two problems: fluid flow (shale gas flow) in porous media and reservoir geomechanics.Extending from a set of diffusivity and displacement based equations [22], the basic fluid flow and geomechanics formulation in this study is shown in Equations ( 1)-(3).Equation ( 1) is based on the Darcy's equation as it is used to describe the fluid transport problem in shale reservoirs [23,24].Equation (2) represents the gas phase pressure This irreversible permeability damage is the hysteresis phenomenon for stress-sensitive permeability.A conceptual hysteresis model of stress-sensitive permeability in the literature also depicted similar phenomenon [17].This hysteresis phenomenon is attributed to the nonlinear consolidation in the stress loading/unloading cycling [14].Sone and Zoback [14] indicated that, as stress is changed, the stress-strain relationship in rock deformation is not always linear, and deformation exhibits nonlinear elastic behaviors [14].The stress sensitive permeability hysteresis correlations obtained in this lab represent the nonlinear elastic rock deformation.It is used to improve the coupled flow and geomechanics model in this study.

Improved Finite Element Numerical Model
An improved numerical model considering nonlinear elasticity is presented to address the stress sensitive effects in shale rocks.The improvement is achieved by adding nonlinear elasticity to a commonly used poroelastic model as in [22].The model is validated with the aforementioned experimental data and can effectively address stress sensitive reservoir characteristics specifically in the Lower Silurian Longmaxi Formation.

Numerical Solution to the Coupled Flow and Geomechanics Problem
Accurate modeling of multiphysics problems in shale reservoirs is critical for shale gas extraction.The physics in shale reservoirs can be effectively depicted by two problems: fluid flow (shale gas flow) in porous media and reservoir geomechanics.Extending from a set of diffusivity and displacement based equations [22], the basic fluid flow and geomechanics formulation in this study is shown in Equations ( 1)-(3).Equation ( 1) is based on the Darcy's equation as it is used to describe the fluid transport problem in shale reservoirs [23,24].Equation (2) represents the gas phase pressure and the coupling between the gas flow problem and the geomechanics problem.Equation ( 3) is the geomechanics equilibrium describing the force balance in the domain of the analyzed volume and at the boundary of the domain.Nonlinear elasticity is described in Equation ( 4), where the Young's modulus E is a function of confining stress σ c .The detailed function is correlated in Equation ( 13) in Section 3.2 based on experimental results.
The full coupling between the flow problem and the geomechanics problem is fulfilled by incorporating the terms b−∅ K s and ρ g b ∂ε v ∂t in the pressure equation Equation ( 2): In Equation ( 1), v is the gas velocity.k is permeability.µ g is gas viscosity.∇p g is pressure gradient of the gas phase.In Equation ( 2), ρ g is gas density.K s is solid grain modulus.b is Biot's coefficient.
∅ is rock porosity.z g is real gas factor.M g is gas molar mass.R is gas constant.T is temperature.t is time.ε v is volumetric strain.In Equation ( 3), σ 0 is initial total stress.λ is first Lame's constant.µ is second Lame's constant.ε is strain tensor.tr(ε) is trace of strain tensor.p g0 is initial gas pressure.I is the second order identity tensor.Some effects in stress sensitive shale rocks are also considered in the numerical model.The Klinkenberg gas slippage effect is incorporated in Equation (5).Shale porosity and permeability are formulated as functions of rock deformation and pore pressure alterations in Equations ( 6) and ( 7): In Equation ( 5), k g is the lab-measured permeability or apparent permeability.b K is a constant coefficient of a certain rock type.p avg is the average pressure.In Equation ( 7), ∅ 0 and k 0 are reference porosity and permeability.γ is a correlation factor [25]. b K in this study is from published data on shale rocks in Longmaxi Formation [26].
In fact, for the given reservoir depth where rock samples were taken in this study, the Klinkenberg effect is not significant and it can be ignored.However, the consideration of this effect in the formulas improves the numerical simulator's modeling capacity for other numerical studies where this effect is significant.
Initial conditions and boundary conditions of the flow and geomechanics problems represent the confined core plug as in the lab setup in Figure 1.Boundaries are shown in Figure 3. Equations ( 8)-(12) show the initial and boundary conditions.n is the normal vector.u is the displacement tensor.t is the boundary traction of the geomechanics problem: Finite element methods are used for the numerical solution of this coupled problem based on the DEAL.II library developed by Bangerth et al. [27].Equation ( 1) is solved with the Raviart-Thomas method [28].Equation ( 2) is solved with the discontinuous Galerkin method [27].Equation ( 3) is solved with the continuous Galekin method [27].The numerical system is iteratively solved by the Newton-Raphson method at each time step, which addresses the nonlinearity of the flow and geomechanics problems.A direct solver UMFPACK with direct sparse LU factorization is applied in the matrix solution [29].
Energies 2018, 11, x 6 of 16 Finite element methods are used for the numerical solution of this coupled problem based on the DEAL.II library developed by Bangerth et al. [27].Equation ( 1) is solved with the Raviart-Thomas method [28].Equation ( 2) is solved with the discontinuous Galerkin method [27].Equation ( 3) is solved with the continuous Galekin method [27].The numerical system is iteratively solved by the Newton-Raphson method at each time step, which addresses the nonlinearity of the flow and geomechanics problems.A direct solver UMFPACK with direct sparse LU factorization is applied in the matrix solution [29].

Modeling of Nonlinear Elasticity
Instead of assuming that the rock has constant elastic modulus as in elasticity in the original model [22], the nonlinear elastic theory considers nonlinear stress-strain relationships [7,14,15].Alteration of elastic modulus during deformation was discussed as a method to describe the nonlinear relationship between stress and strain [30].The alteration of elastic modulus during stress changes can be used to describe the nonlinear stress-strain relationship in rock deformation.

Elastic Modulus Alteration
In this study, following the concept of changing elastic modulus, based on experimental results in Figure 2, Young's modulus (as an important form of elastic modulus) is treated as a nonlinear function of confining stress to incorporate nonlinear elasticity in the numerical model.For each of the five core samples, Young's modulus value calibration is conducted at each discrete confining stress value for both stress loading and unloading.The lab measured apparent permeability is used as the calibration criterion, and the Young's modulus value is being adjusted in the numerical simulation until the simulated average permeability is equal to the measured permeability.The measured permeability is corrected by Equation (5) for the gas slippage effect.Once these two permeability values match, the Young's modulus is recorded and associated with the corresponding confining stress and corresponding stress loading/unloading state.The calibration workflow is shown in Figure 4.The calibrated correlations between Young's modulus and confining stress are shown in Figure 5. Curve fittings with cubic equations are also provided.Equation (13) shows the cubic equations for Young's modulus as functions of stress for all cores.Table 2 shows the cubic coefficients for all five cores.In cubic equations, confining stress is in MPa and calibrated Young's modulus is in Pa.

Modeling of Nonlinear Elasticity
Instead of assuming that the rock has constant elastic modulus as in elasticity in the original model [22], the nonlinear elastic theory considers nonlinear stress-strain relationships [7,14,15].Alteration of elastic modulus during deformation was discussed as a method to describe the nonlinear relationship between stress and strain [30].The alteration of elastic modulus during stress changes can be used to describe the nonlinear stress-strain relationship in rock deformation.

Elastic Modulus Alteration
In this study, following the concept of changing elastic modulus, based on experimental results in Figure 2, Young's modulus (as an important form of elastic modulus) is treated as a nonlinear function of confining stress to incorporate nonlinear elasticity in the numerical model.For each of the five core samples, Young's modulus value calibration is conducted at each discrete confining stress value for both stress loading and unloading.The lab measured apparent permeability is used as the calibration criterion, and the Young's modulus value is being adjusted in the numerical simulation until the simulated average permeability is equal to the measured permeability.The measured permeability is corrected by Equation (5) for the gas slippage effect.Once these two permeability values match, the Young's modulus is recorded and associated with the corresponding confining stress and corresponding stress loading/unloading state.The calibration workflow is shown in Figure 4.The calibrated correlations between Young's modulus and confining stress are shown in Figure 5. Curve fittings with cubic equations are also provided.Equation (13) shows the cubic equations for Young's modulus as functions of stress for all cores.Table 2 shows the cubic coefficients for all five cores.In cubic equations, confining stress σ c is in MPa and calibrated Young's modulus E c is in Pa.

Discussion
Based on Figure 5, as confining stress increases, for the stress loading curve, Young's modulus increases since rock stiffness increases as the rock becomes more compacted.For the stress unloading curve, Young's modulus decreases as the confining stress decreases and the rock becomes less compacted.In Figure 5c,e, at the maximum confining stress of 22 MPa, the Young's modulus is not at its largest value.This is because, when the permeability was measured at the maximum confining stress in the lab, although the confining stress was at its greatest value, the stress unloading was initiated and the rock sample was no longer in a fully compressed state, resulting in decreased rock stiffness and decreased Young's modulus.Also, it is noted that for the same confining stress value, the elastic modulus associated with stress loading is greater than the elastic modulus associated with stress unloading.It indicates that the shale rock elastic modulus is not only affected by the confining stress value but also by the stress history path.This observation of elastic modulus's dependence on both stress history path and in-situ stress state is an extension of previous observations that permeability is dependent on both in-situ stress state and reservoir stress history path [31].Based on Equations ( 6) and (7), greater strain change leads to greater permeability decrease.For the same confining stress value, permeability associated with stress unloading is smaller than that with stress loading, which indicates greater strain change for stress unloading permeability.Thus, the numerical model needs to adopt a smaller elastic modulus value for stress unloading to allow for greater strain change under the same confining stress state.Note that these observation and explanation are only based on shale samples from Longmaxi Formation, and changing to other shale samples from new sampling locations can lead to different results.As indicated by Sone and Zoback, Young's modulus values in stress loading/unloading processes for shale samples can vary in a rather large range, and the Young's modulus measurements are affected by many factors such as measurement methods for elastic modulus, the strategy for shale core extraction, decompaction of shale samples, and in-situ over-pressurization state [14].This means that elastic modulus alteration in stress loading/unloading is a complicated problem and the patterns of elastic modulus alteration in this section cannot be generalized.
In Figure 5, as a reference, the constant Young's moduli in all samples assuming linear elasticity are also presented by green dashed flat lines.These moduli do not change with confining stress changes.The constant Young's modulus in each sample is taken from the initial point of the stress loading curve, since the rock is not affected by the stress path at the beginning and its elastic modulus can be used as the value not affected by any nonlinear stress-strain relationship during rock deformation.Equation ( 4) is now replaced with Equation (13) to enable the coupled flow and geomechanics simulator to model nonlinear elasticity:

Numerical Study
Using the improved numerical model considering nonlinear elasticity, a set of numerical results is obtained for detailed characterization of pore pressure, permeability, and porosity within the analyzed core volume to study the shale stress sensitivity.The numerical simulation models the stress and pressure evolution in the core plugs flooded as Figure 1.The 1D and 2D spatial distributions of the results are reported as in Figure 6.Numerical simulation is conducted for all five shale samples.

Numerical Study
Using the improved numerical model considering nonlinear elasticity, a set of numerical results is obtained for detailed characterization of pore pressure, permeability, and porosity within the analyzed core volume to study the shale stress sensitivity.The numerical simulation models the stress and pressure evolution in the core plugs flooded as Figure 1.The 1D and 2D spatial distributions of the results are reported as in Figure 6.Numerical simulation is conducted for all five shale samples.

Quantification of Effects of Nonlinear Elasticity
Numerical results obtained by the improved model are compared with numerical results simulated by the basic model [22] which only considers linear elasticity so that effects of nonlinear elasticity can be quantified.
In hydrocarbon production, effective stress generally increases as the hydrocarbon in the pores are depleted, and initial reservoir condition before hydrocarbon depletion has the lowest effective stress.In this study, as shown by flat dashed lines in Figure 5, the confining stress of 5 MPa at the beginning of the stress loading curve is used to represent the intact/initial shale reservoir condition.In the reference case, the Young's modulus at this initial state is used in the modeling that only considers linear elasticity, and this value is not changing with confining stress changes.
In contrast, the improved modeling that considers nonlinear elasticity has a changing Young's modulus as confining stress is altered.The increase of confining stress is used to represent the increase of effective stress caused by shale gas production in the reservoir.Several confining stress values in the stress loading process (13, 16 and 19 MPa) are used for the numerical modeling that considers nonlinear elasticity.Young's modulus values used in this quantification study are in Table 3, corresponding to the calibrated correlations in Equation (13).Numerical simulations of both the reference linear elasticity scenarios and the nonlinear elasticity scenarios are carried out.The spatial distribution results averaged over the 1D axis as in Figure 6 are shown in Tables 4-6 with values averaged over the 1D profiles.Additionally, following the 2D profiling scheme, Figures 7 and 8 are plotted.The reported results are average permeability, average porosity, average minimum principal stress, and average pressure.The numerical results for stress presented in this section are the effective stress excluding pore pressure.
The difference between linear elasticity simulation results and nonlinear elasticity simulation results (Equation ( 14)) becomes larger as confining stress increases.This trend is observed for all results of , ∅ , , and ̅ , which indicates that, as effective stress caused by hydrocarbon depletion becomes larger during gas extraction, the effects of nonlinear elasticity on permeability of rock parameters and stress/pressure spatial evolution become more significant.This is because in the nonlinear elasticity model, the increment of Young's modulus is caused by stress increase, and the increased Young's modulus leads to greater stress magnitude in the analyzed volume.The stress then significantly affects values of the stress sensitive rock properties.However, in the linear elasticity model, the initial Young's modulus is kept as constant and does not increase as stress increases, which leads to underestimated in-situ stress magnitudes.

Quantification of Effects of Nonlinear Elasticity
Numerical results obtained by the improved model are compared with numerical results simulated by the basic model [22] which only considers linear elasticity so that effects of nonlinear elasticity can be quantified.
In hydrocarbon production, effective stress generally increases as the hydrocarbon in the pores are depleted, and initial reservoir condition before hydrocarbon depletion has the lowest effective stress.In this study, as shown by flat dashed lines in Figure 5, the confining stress of 5 MPa at the beginning of the stress loading curve is used to represent the intact/initial shale reservoir condition.In the reference case, the Young's modulus at this initial state is used in the modeling that only considers linear elasticity, and this value is not changing with confining stress changes.
In contrast, the improved modeling that considers nonlinear elasticity has a changing Young's modulus as confining stress is altered.The increase of confining stress is used to represent the increase of effective stress caused by shale gas production in the reservoir.Several confining stress values in the stress loading process (13, 16 and 19 MPa) are used for the numerical modeling that considers nonlinear elasticity.Young's modulus values used in this quantification study are in Table 3, corresponding to the calibrated correlations in Equation (13).Numerical simulations of both the reference linear elasticity scenarios and the nonlinear elasticity scenarios are carried out.The spatial distribution results averaged over the 1D axis as in Figure 6 are shown in Tables 4-6 with values averaged over the 1D profiles.Additionally, following the 2D profiling scheme, Figures 7 and 8 are plotted.The reported results are average permeability, average porosity, average minimum principal stress, and average pressure.The numerical results for stress presented in this section are the effective stress excluding pore pressure.
The difference between linear elasticity simulation results and nonlinear elasticity simulation results (Equation ( 14)) becomes larger as confining stress increases.This trend is observed for all results of k, ∅, σ, and p, which indicates that, as effective stress caused by hydrocarbon depletion becomes larger during gas extraction, the effects of nonlinear elasticity on permeability of rock parameters and stress/pressure spatial evolution become more significant.This is because in the nonlinear elasticity model, the increment of Young's modulus is caused by stress increase, and the increased Young's modulus leads to greater stress magnitude in the analyzed volume.The stress then significantly affects values of the stress sensitive rock properties.However, in the linear elasticity model, the initial Young's modulus is kept as constant and does not increase as stress increases, which leads to underestimated in-situ stress magnitudes.In the results above, it is noted that the in-situ stress and pore pressure are altered as confining stress changes.In the stress loading process, as the confining stress increases, the stress increases and the pore pressure decreases.The decreasing pore pressure is due to the decreased porosity as shown in the numerical results: the increased confining stress leads to decreased pore volume which facilitates pore fluid depletion.The increased principal stress is because of the increased compression caused by the increased confining stress.

Effect of Nonlinear Elasticity on Permeability
The permeability presented in these numerical results are the corrected liquid permeability instead of the apparent permeability.Thus the value is generally smaller than the lab measured permeability data.Comparing with results from the linear elasticity model, consideration of nonlinear elasticity decreases the shale permeability value .This is because the Young's modulus in the nonlinear elasticity model is increased due to confining stress increase, which leads to greater effective stress magnitudes.Greater in-situ stress indicates larger damage to permeability.In linear elasticity modeling, the overestimated permeability will lead to inaccurate modeling of shale gas production and largely decrease modeling effectiveness.

Effect of Nonlinear Elasticity on Effective Stress
Compared with linear elasticity model, consideration of nonlinear elasticity increases the effective stress within the analyzed volume.This is because Young's modulus values are larger in the nonlinear elasticity model, resulting in greater stress simulation results.
This observation indicates that, as linear elasticity modeling does not consider the Young's modulus increase due to stress increase, the assumption of linear elasticity leads to underestimated effective stress during the shale gas depletion process.

Ranking of Effects
Based on the differences in percentages, the effects of nonlinear elasticity on numerical results are ranked from the greatest to the smallest as permeability > stress > porosity > pore pressure.Regarding numerical results of all five samples with all three confining stresses (13 MPa, 16 MPa, and 19 MPa), the consideration of nonlinear elasticity can make the numerical result for permeability differ by one order of magnitude, with the difference ranging from 1.6 times (160.19%difference in sample LS1-2-1-X for 13 MPa ) to 53 times (5343.93%difference in sample LS2-4-2 for 19 MPa ).The consideration of nonlinear elasticity makes the numerical result for effective stress differ by 76% to 85%.The consideration of nonlinear elasticity makes the numerical result for porosity differ by 4% to 39%.The consideration of nonlinear elasticity makes the numerical result for pore pressure differ by 0.03% to 1.73%, which is very small: In the results above, it is noted that the in-situ stress and pore pressure are altered as confining stress changes.In the stress loading process, as the confining stress increases, the stress increases and the pore pressure decreases.The decreasing pore pressure is due to the decreased porosity as shown in the numerical results: the increased confining stress leads to decreased pore volume which facilitates pore fluid depletion.The increased principal stress is because of the increased compression caused by the increased confining stress.

Effect of Nonlinear Elasticity on Permeability
The permeability presented in these numerical results are the corrected liquid permeability instead of the apparent permeability.Thus the value is generally smaller than the lab measured permeability data.Comparing with results from the linear elasticity model, consideration of nonlinear elasticity decreases the shale permeability value k.This is because the Young's modulus in the nonlinear elasticity model is increased due to confining stress increase, which leads to greater effective stress magnitudes.Greater in-situ stress indicates larger damage to permeability.In linear elasticity modeling, the overestimated permeability will lead to inaccurate modeling of shale gas production and largely decrease modeling effectiveness.

Effect of Nonlinear Elasticity on Effective Stress
Compared with linear elasticity model, consideration of nonlinear elasticity increases the effective stress σ within the analyzed volume.This is because Young's modulus values are larger in the nonlinear elasticity model, resulting in greater stress simulation results.
This observation indicates that, as linear elasticity modeling does not consider the Young's modulus increase due to stress increase, the assumption of linear elasticity leads to underestimated effective stress during the shale gas depletion process.

Ranking of Effects
Based on the differences in percentages, the effects of nonlinear elasticity on numerical results are ranked from the greatest to the smallest as permeability > stress > porosity > pore pressure.Regarding numerical results of all five samples with all three confining stresses (13 MPa, 16 MPa, and 19 MPa), the consideration of nonlinear elasticity can make the numerical result for permeability differ by one order of magnitude, with the difference ranging from 1.6 times (160.19%difference in sample LS1-2-1-X for 13 MPa σ c ) to 53 times (5343.93%difference in sample LS2-4-2 for 19 MPa σ c ).The consideration of nonlinear elasticity makes the numerical result for effective stress differ by 76% to 85%.The consideration of nonlinear elasticity makes the numerical result for porosity differ by 4% to 39%.The consideration of nonlinear elasticity makes the numerical result for pore pressure differ by 0.03% to 1.73%, which is very small: Di f f erence = |Property f rom nonlinear elasticity − Property f rom linear elasticity| Property f rom nonlinear elasticity (14)

Nonlinear Elasticity Modeling Results and Discussion
Numerical results based on the improved model considering nonlinear elasticity are discussed in this section.Core sample LS1-16-1 is specifically analyzed for profiles of numerical result distribution.The 1D and 2D characterizations indicate that the numerical model in this study can describe the heterogeneous rock characteristics in shale core samples, while the heterogeneity is usually neglected in analytical models.

Effective Stress
In Figure 7 effective stress distribution exhibits non-uniformity along both the axial and transverse directions.The transverse non-uniformity is the least significant in the middle of the studied area at around 2.2 cm along the axial direction.The effective stress monotonically increases along the axial direction from the inlet on the left to the outlet on the right.This is because, before core flooding, the entire volume has the initial pressure equal to inlet pressure.In core flooding, the pressure depletion near the outlet results in large rock deformation and large stress.Besides, the overall effective stress increases with confining stress increase.

Permeability
In Figure 8, permeability distribution is also non-uniform in the analyzed volume.The permeability is liquid permeability excluding the gas slippage effect and is smaller than lab measured apparent permeability.The non-uniform permeability distribution is caused by the non-uniform stress distribution, as permeability is significantly affected by stress and pore pressure.1D results show that larger confining stress leads to greater damage to permeability.Also, the permeability decreases from the inlet to the outlet, as the effective stress increases from the inlet to the outlet.The first derivative of permeability represents how fast the permeability decreases along the axial direction.A larger absolute value of permeability indicates a faster permeability decrease along the axis.As confining stress increases, the permeability decrease along the axis becomes slower.This indicates that, as confining stress increases, although permeability damage is increased, the effect of axial stress distribution on axial permeability distribution is decreased.

Conclusions
In this work, based on the calibration with lab data, an improved finite element model considering nonlinear elasticity is proposed in a case study of the Longmaxi shale in China.The quantification of effects of nonlinear elasticity on stress sensitive shale rock characteristics is also provided.The model presented in this work specifically focuses on the stress sensitive shale gas reservoir with high hydrocarbon production potential in the Lower Silurian Longmaxi Formation.Some conclusions of this study are drawn.
(1) Nonlinear elasticity should be considered in modeling of stress sensitive shale reservoir rocks as it largely affects the results.The concept of changing elastic modulus can be used to model nonlinear elasticity, and its usage is validated with the experimental study in this work.(2) Based on experimental validation of nonlinear elastic modulus change in the numerical model, in Longmaxi Formation, the shale elastic modulus is dependent on both in-situ stress state and the stress history path.
(3) The consideration of nonlinear elasticity affects the stress distribution in the analyzed volume, and consequently affects stress-related characteristics of permeability, effective stress, porosity, and pore pressure.Effects of nonlinear elasticity on the numerical results are ranked as permeability > stress > porosity > pore pressure.(4) Comparing with linear elasticity, the numerical result for permeability with the consideration of nonlinear elasticity can differ by one order of magnitude in the Longmaxi Formation shale samples.In the modeling of Longmaxi shale, without considering the nonlinear elasticity effect, the permeability can be overestimated by 1.6 to 53 times in the analyzed samples.(5) Comparing with nonlinear elasticity, the assumption of linear elasticity underestimates the magnitude of the effective stress caused by hydrocarbon depletion.This is because linear elasticity does not consider the increase of elastic modulus during the process of stress increase.(6) As the confining stress increases, although the damage to permeability increases, the effect of stress on the spatial distribution of permeability is decreased.(7) Shale rock properties of porosity and permeability are not uniformly distributed in the analyzed shale core volume that is flooded by gas in the experimental study.The heterogeneous characteristics was modeled by the improved numerical model, while heterogeneity is usually ignored in analytical models.

Figure 3 .
Figure 3. Boundaries in the numerical model.

Figure 3 .
Figure 3. Boundaries in the numerical model.

Figure 4 .Figure 5 .
Figure 4. Nonlinear elasticity calibration workflow based on adjustment of elastic modulus values.

Figure 4 .
Figure 4. Nonlinear elasticity calibration workflow based on adjustment of elastic modulus values.

Figure 4 .Figure 5 .
Figure 4. Nonlinear elasticity calibration workflow based on adjustment of elastic modulus values.

Table 2 .
Cubic equation correlations for numerical model calibration.

Table 3 .
Calibrated Young's modulus values used in the quantification numerical study.

Table 4 .
Average property values for confining stress of 13 MPa and comparison with the linear elasticity case.

Table 6 .
Average property values for confining stress of 19 MPa and comparison with the linear elasticity case.