A Computational Scheme for Evaluating the Stress Field of Thermally and Pressure Induced Unconventional Reservoir

The ﬂ uid ﬂ ow connecting the hydraulic fracture and associated unconventional gas or oil reservoir is of great importance to explore such unconventional resource. The deformation of unconventional reservoir caused by heat transport and pore pressure ﬂ uctuation may change the stress ﬁ eld of surrounding layer. In this paper, the stress distribution around a penny-shaped reservoir, whose shape is more versatile to cover a wide variety of special case, is investigated via the numerical equivalent inclusion method. Fluid production or hydraulic injection in a subsurface resource caused by the change of pore pressure and temperature within the reservoir may be simulated with the help of the Eshelby inclusion model. By employing the approach of classical eigenstrain, a computational scheme for solving the disturbance produced by the thermally and pressure induced unconventional reservoir is coded to study the e ﬀ ect of Biot coe ﬃ cient and some other important factors. Moreover, thermo-poro transformation strain and arbitrarily orientated reservoir existing within the surrounding layer are also considered.


Introduction
The International Energy Agency projections indicate that an increasing supply gap for liquid fuels will arise if discovery and investment in new fields are unable to meet the urgent requirements with the rapidly increasing consumption of resources. Although the unconventional sources will lower energy costs and allow the substitution of natural gas for coal for electrical power production, the lack of equipment and experience for exploring the complicated geological structure may inhibit or slow down development. The prospect of higher prices accompanying a supply crunch may lead to an acceleration of investment in the unconventional resources, such as shale gas and other frontier oils [1][2][3]. Hydraulic fracturing, whose process involves the high-pressure injection of fracking fluid into a wellbore to create cracks in the deep-rock formation, is a well stimulation technique commonly utilized in unconventional reservoirs.
Solving the thermally and pressure induced reservoir is of great importance to petroleum industry and the associ-ated applications [4][5][6]. The fluid flow connecting the fracture and surrounding rocks, heat transport, and temperature changes within the resource may significantly disturb the stress fields and deformation of the unconventional oil reservoirs including the shale and tight sand reservoir. The well-known Eshelby solution has been widely utilized to determine the effects of inhomogeneities which differ from the regional material [7,8]. Eshelby's technique can be also applied to deal with the elastic fields of stress, strain, displacement, deformation gradient around the intrusions, reservoirs, and underground structures due to the pore pressure or temperature changes [9,10]. Zhou et al. investigate the background, applications, and future research directions of inclusion problem, as well as the classical methods including the analytical and numerical approach to treat associated problems [11].
Extensively studies on the geostructures and applications of hydraulic fracture have been explored over the decades. An analytical solution considering the porosity and permeability of gas storage is presented and may be applied to study the corresponding properties of shale gas reservoirs with noncircular nanopores [12]. Cai et al. explore the mechanism and effect of shale bedding of supercritical carbon dioxide hydro-jet fracturing using the experimental method [13][14][15]. The problems of two interacting ellipsoidal inclusions and inhomogeneities with arbitrary orientation are solved to determine the elastic stress and strain components using Eshelby's technique, and such method can be extended to discuss the pressure induced fracture [16,17]. The failure mechanism, crack initiation, and propagation of brittle granite are evaluated on the basis of discrete element method [18], as well as the thermally induced microcracking [19]. The interpolations of the experimental data are proposed to obtain the characterization of fractured medium and corresponding crack porosity according to Eshelby-Cheng model [20]. Kanin et al. adopt the Carter's leak-off model to assess the effect of pore fluid leak-in, near-tip cavity, and the penny-shaped crack [21]. Mauludin et al. investigate the effect of mismatch fracture properties and simulate the circular capsule with different types surrounding by the mortar media utilizing the cohesive-zone model [22]. The hydrodynamic flow model is proposed to investigate the properties and distributions of heterogeneous fluid, the interaction between the liquid and wall, and liquid-liquid slip [23]. Eghbalian et al. use the approach of representative elementary volume to estimate the failure mechanism of brittle rocks containing the microcracks and nanopores [24]. Pang et al. study the mechanism of kaolinite swelling produced by CO2 and water with the help of molecular simulation method [25]. Cai et al. develop a numerical method to evaluate hydraulic fracture propagation with the help of the flow-coupled cohesive interface coupled method [26].
The problem of deformation and stress disturbance around a reservoir due to the heat transport and fluid migration may be solved on the basis of the Eshelby model, which involves the concept of eigenstrain to handle the inelastic issue [27][28][29]. In this study, the stress fields around a penny-shaped reservoir with thermo-porous eigenstrains in a full space are investigated via the numerical equivalent inclusion method. Section 2 introduces the thermal and porous eigenstrains and the classical equivalent inclusion method, as well as the corresponding Eshelby tensor. The computational scheme to treat the arbitrarily orientated reservoir surrounding by the diffident soil and rock materials is proposed in Section 2. Parametric studies and examples on the Biot coefficient, distribution of a reservoir, effect of linear thermal expansion, temperature changes, pore pressure fluctuation, and elastic stress disturbance are provided. The applications of proposed solution and computational scheme are given in Section 3. Conclusions are presented in Section 4.

The Eshelby Model
2.1. Thermo-Poro Transformation Strain. The formulation of poroelasticity is a branch of the classical elastic theory and can be extended to study the thermal problem for a continuous porous material. Mura [30] extends the celebrated work by Eshelby [7,8] and uses the term "eigenstrain" to represent the inelastic strain, including the thermal expansion, pore pressure change, misfit, and dislocation. The distribution of uniform or nonuniform eigenstrains in a material will lead to the corresponding eigenstresses and also affect the mechanical properties of the surroundings. The well-known eigenstrain method can be widely utilized to handle the void, crack, and other material defects as strain released during geometry changes. Therefore, changes due to heat expansion or cold contraction in a porous medium may be computed with choosing the proper eigenstrain via the Eshelby model [31]. The eigenstrain caused by the pore pressure fluctuation and temperature change may be expressed as The terms μ and ν refer to the shear modulus and Poisson's ratio, Kronecker delta δ ij is 1 if the variables i, j are equal, and 0 otherwise. The Biot parameter α = 1 − C s /C, and the compressibility of the solid phase or called the unjacketed coefficient of compressibility is denoted as C s , which is usually considered identical to the compressibility of the grains [32]. The corresponding jacketed compressibility coefficient C means the drained compressibility of the porous medium. Furthermore, the constitutive equation of thermo-poro eigenstrain (1) also depends on the unstrained state of pore pressure P 0 and temperature T 0 and the strained state of pressure P and temperature T. The term λ represents the coefficient of linear thermal expansion, which characterizes the ability of an object to expand under the effect of temperature elevation. It is worth noting that if the porous reservoir is compliant, which is compared to the solid constituents, α = 1 − C s /C~1, whereas if the resource is stiff α = 1 − C s /C~0.

Formulation.
Eshelby's work [7,8] has shaped the fields of micromechanics and provided some outstanding and fundamental solutions to address the problems in solid with defects, especially for the fracture mechanics, plasticity and engineered composite materials, etc. It is noted that Eshelby's equivalent inclusion method in elasticity is one of the most famous and well-known micromechanical approach to analyze the deformation of arbitrarily shaped inhomogeneity, which is restrained by its surroundings with different elastic constants. The general ellipsoidally shaped inhomogeneity is versatile enough to cover a wide variety of special shape including the penny-shaped reservoir with thermal and porous transformation strain. In this case, an isotropic, linear, and full-space medium applied with remote external load σ 0 ij contains an ellipsoidal inhomogeneity, whose elastic constants C * ijkl differ from those of the rest of the surrounding matrix C ijkl . In this case, when letting three axes of the ellipsoid a 1 ≠ a 2 ≫ a 3 , the general shape of ellipsoid transforms into the penny-shaped reservoir, as shown in Figure 1.

Lithosphere
The nonlinear inhomogeneity problem may be converted into linear inclusion problem via the equivalent inclusion method. The constitutive equation of equivalent inclusion method is wherein the term ε 0 ij denotes the remote strain field, which can be calculated in absence of the inhomogeneity according to Hook's law, σ 0 ij = C ijkl ε 0 kl . In addition, ε P ij is the preliminary eigenstrain distributed inside the reservoir and is given in Equation (1). The equivalent eigenstrain ε * * ij may be expressed as The fictitious eigenstrain ε * ij can be solved according to the constitutive equation in Equation (2). Moreover, Equation (2) can be rewritten in terms of the classical Eshelby strain tensors S ijkl . The detailed formulation of Eshelby tensors is listed in [30].
Solving above equation can obtain the fictitious eigenstrain ε * ij and equivalent eigenstrain ε * * ij . Equation (4) is enforced at a set of particularly chosen points, resulting in a system of algebraic equations, which can be numerically solved. Thus, the stress field of thermally and pressure induced unconventional reservoir may be determined by adopting Eshelby's equivalent inclusion method, as shown in Figure 2.
When the penny-shaped reservoir and surrounding rock matrix are both isotropic materials, Equation (4) may be expressed as [30] 2μ wherein the parameters λ * and μ * refer to the Lame constants of the penny-shaped reservoir and λ and μ denote the corresponding Lame constants for the rock matrix. In order to obtain the shear components of ε * * ij , the deviatoric strains are introduced.
By adopting the above expressions, Equation (5) arrives at Arbitrary orientation After the algebraic operation, the shear parts of ε * * ij may be calculated as In this special case, the reservoir may be assumed be a penny-shaped inhomogeneity, whose geometric shape is a 1 = a 2 and a 3 = 0. Then, Equation (5) becomes and the equivalent eigenstrain can be solved as To make the case a bit more general, the following examples and discussions of the geometry for the resource inhomogeneity are set to be a 1 ≠ a 2 ≫ a 3 . Due to the arbitrarily orientated reservoir existing within the surrounding rock after hydraulic injection, it is necessary to apply the Euler angles to evaluate this problem. According to Euler's rotation theorem, any rotation can be measured by utilizing three angles. Thus, in conjunction with the equivalent inclusion method, the stress distributions of thermally and pressure induced reservoir with arbitrary orientation can be numerically solved.

Discussion
In order to calculate the stress distributions of thermally and pressure induced reservoir surrounding by the unconventional oil reservoir, a computational code is programmed Figure 2: The equivalent inclusion method to treat the thermally and pressure induced reservoir. 4 Lithosphere to investigate some fundamental parameters affecting the elastic stress fields. In this section, dimensionless factors are employed to study the penny-shaped reservoir caused by thermo-porous eigenstrains. The practical parameters including the Young's modulus, Poisson's ratio, remote applied loads, coefficient of linear thermal expansion, Biot's coefficient of the unconventional reservoir, and rock matrix are considered. Due to the randomness of reservoir, the Euler angles to represent its orientation are also involved to determine the Von Mises Stress of reservoir caused by the heat elevation or cold contraction and pressure raising or dropping before and after the hydraulic injection. To make the case a bit more general, the values of Biot's coefficient are on the basis of Bakken shale samples from Williston Basin; the text measurements of samples are detailed in [33]. Thus, the practical factors of Young's modulus, Poisson's ratio, pore pressure change, temperature fluctuation, coefficients of linear thermal expansion, remote applied loads, and rotation angles of reservoir can be input into the codes to evaluate this problem. For the purpose of investigating each factor affecting the final results, the general parameters are set as follows. Reservoir is embedded in an infinitely extended shale matrix with remote loads applied σ 0 11 = σ 0 22 = σ 0 33 = 1 at infinity, the size of reservoir is a 1 = 1, a 2 = 0:75, a 3 = 0:01, Young's modulus and Poisson's ratio for the reservoir matrix and surrounding rock E = 30000Mpa and ν = 0:3, the corresponding coefficient of linear thermal expansion λ = 3 × 10 -5 , Biot's coefficient for the unconventional reservoir α = 1 − C s /C = 0:7, pore pressure change ΔP = 15, and temperature fluctuation ΔT = −50. The target line is from -2 to 2, and the Von Mises stress distribution and the each stress components along this line is obtained.
The effect of parameter on the Young's modulus for the surrounding rock matrix with respect to those of reservoir E reservoir = 30000 Mpa may be computed according to the formulation detailed in Section 2. The calculated Von Mises stress for a penny-shaped reservoir starting from x 1 /a 1 = −2 to 2 is shown in Figure 3(a). Results are displayed as the ratio of E matrix /E reservoir ; the Von Mises stress reaches the maximum points at x 1 /a 1 = −1, 1 when E matrix = 5/3E reservoir and minimum points at x 1 /a 1 ≤ −1, or ≥ 1. That means the jump between the matrix and reservoir has its maximum value when E matrix is larger than E reservoir . The effect of Poisson's ratio on the Von Mises stress is plotted in Figure 3 Results are given for two parameters based on the equation of eigenstrain caused by the pore pressure fluctuation and temperature change, which is related to the coefficients of Biot and linear thermal expansion. Figures 4(a) and 4(b) show the normal stress components for Biot's coefficient α = 0:7 and the coefficient of linear thermal expansion λ = 1 × 10 -5 , 3 × 10 -5 , respectively. The difference between 4(a) and 4(b) is that the stresses σ 11 , σ 22 when λ = 1 × 10 -5 show smaller than for the λ = 3 × 10 -5 . It is noted that the stress σ 33 exhibits almost the same tendency inside and outside the reservoir due to the target route lying on the x-axis. In reality, the stress distribution inside the inhomogeneity would like to be disturbed by surrounding geostructures such as fault and cracked layer. In this case, the reservoir is assumed to be perfectly bounded by the rock matrix, without being disturbed by any other geostructures or neighbor resources. The total tendency of Figures 4(a) and 4(b) is seen to be almost equal; only the maximum and minimum stress values differ.
As expected, before or after the hydraulic fracturing, there is a temperature difference between the injected fluid and surrounding matrix underground. Figure 5 displays the predicted normal stress distributions of the penny-shaped reservoir due to temperature changes  5 Lithosphere produced by hydraulic injection. Figure 5(a) presents the normal stresses after injection leading to the temperature decreasing of reservoir matrix caused by the injected cold water. In like manner, Figure 5(b) describes the normal stress distributions when the temperature elevates. It can be seen that stress σ 33 keeps a constant as the temperature heating or decreasing as opposite to the stresses σ 11 and σ 22 . The results tell that the stress disturbance may be influenced by the injected cold fluid, and the associated stress field of reservoir would like to be strongly disturbed leading to the stress concentration. The stresses inside the reservoir show a big difference among each stress component if the shape of resource is irregular.
As we know, the pressure of reservoir formation would like to increase with the injection of fracturing fluid until it peaks at the breakdown pressure. In order to compare the pressure dropping and increasing, the last example describes the pore pressure fluctuation affecting the total stress field after or before the fracturing work. Figures 6(a) and 6(b), respectively, exhibit the predicted normal stress parts of the reservoir produced by the change of porous eigenstrains, which means the pressure dropping ΔP = P − P 0 = −30 and pressure elevating ΔP = P − P 0 = 30. In Figure 6(a), the stresses σ 11 and σ 22 reach their maximum value at the interface between the resource and formation when ΔP = −30, while the stresses σ 11 and σ 22 when ΔP = 30 in Figure 6(b) are lower than those of Figure 6(a). It is concluded that the pore pressure of reservoir matrix dropping or elevating may both influence the stress components σ 11 and σ 22 at this special case, but stress σ 33 in Figures 6(a) and 6(b) shows no change due to any pore pressure fluctuation along the target line from x 1 /a 1 = −2 to 2 because of the special shape of reservoir.

Conclusions
Fluid production or hydraulic injection in a subsurface resource would like to change the pore pressure and temperature within the reservoir, which generally induces the alteration of stress distribution within and around the reservoir. The general ellipsoid is more versatile to cover a wide variety of special case, including the penny-shaped inclusion. Using the classical Eshelby inclusion method to solve the induced stress changes within an ellipsoidal resource in a full space or half space is an effective way to deal with the geostructure problem. In this paper, a computational scheme for treating the stress distributions a penny-shaped reservoir with thermo-porous eigenstrains in an infinitely extended formation is explored with the help of numerical equivalent inclusion method. Parametric studies on the Biot coefficient, effect of linear thermal expansion, temperature changes, pore pressure fluctuation, Young's modulus, and Poisson's ratio of the unconventional reservoir are performed. Results are displayed as the ratio of E matrix /E reservoir , the Von Mises stress reaches the maximum points at x 1 /a 1 = −1, 1 when E matrix = 5/3E reservoir , and minimum points at x 1 /a 1 ≤ −1, or ≥ 1. That means the jump between the matrix and reservoir has its maximum value when E matrix is larger than E reservoir . Some other results and discussions are detailed in Section 3.

Data Availability
The processed data required to reproduce these findings cannot be shared at this time as the data also forms part of an ongoing study.

Conflicts of Interest
The author(s) declare(s) that there is no conflict of interest regarding the publication of this paper.