Analysis of flow characteristics of the dual media shale gas reservoir with the elastic outer boundary

In order to more accurately describe the seepage characteristics of shale gas reservoirs, in this paper, an elastic outer boundary condition is introduced, and a new dual-media shale gas seepage model is established to describe the seepage characteristics of shale gas reservoirs more accurately, while considering the adsorption and desorption process. Combining Laplace transformation and Similar Structure Method, the solution of the percolation model is obtained in Laplace space. Furthermore, the solutions are inverted into real space with the Stehfest numerical inversion method. Type curves of dimensionless pressure and dimensionless pressure derivative can be used to evaluate the reservoir characteristics. The results show that the conventional three kinds of outer boundary conditions (infinite, constant pressure and closed) are actually three special cases of elastic outer boundary. The introduction of elastic outer boundary conditions not only expands the scope of data interpretation, but also closer to the actual situation of the outer boundary of the reservoir. The theory of similar structure greatly simplifies the complex process of solving the model and will have a far-reaching impact on the development of well test analysis software in the future. Elastic outer boundary conditions are introduced to the seepage model of shale gas reservoir for the first time. Elastic outer boundary conditions include three special outer-boundary conditions (closed, constant pressure, infinite), and it is more general. The introduction of elastic outer-boundary conditions provides a new perspective and method to make a more complete theoretical chart. Elastic outer boundary conditions are introduced to the seepage model of shale gas reservoir for the first time. Elastic outer boundary conditions include three special outer-boundary conditions (closed, constant pressure, infinite), and it is more general. The introduction of elastic outer-boundary conditions provides a new perspective and method to make a more complete theoretical chart.


Introduction
Elasticity has been applied in all fields, but it has different meanings in different fields. In physics, elasticity shows how an object moves or deforms under the action of external forces; In economics, Alfred Marshall [1] firstly introduced the concept of elasticity to express the degree of the consumers and producers' react to changes in prices. The concept of elasticity can be applied to all variables with causal relationship. Furthermore, Woods [2] gave an expression of the elastic relationship between independent variable x and dependent variable y if y = f (x) is a differentiable function and f (x) ≠ 0: where is the elasticity of the function y = f (x) with respect to x , and y x 0 is the elasticity of the function y = f (x) with respect to x at point x 0 .
According to Eq. (1), in oil and gas reservoir engineering, the concept of the elastic outer boundary is proposed, and the elastic coefficient of pressure difference P = p 0 − p(r, t) with respect to radial radius r is defined as follows [3] P r reflects the influence of the change amplitude Δr r of radial radius r on the change amplitude ΔP P of pressure difference P, that is, it reflects the intensity of ln P relative to the change of ln r . Specifically, P r (r 1 , t 1 ) indicates that at r = r 1 , t = t 1 , when the time t = t 1 remains constant and ln r changes by 1%, ln P changes approximately by P r (r 1 , t 1 )% . And the same is true for P t . According to the definition of P r and P t , it is not difficult to find that: (1) The elastic coefficient (function) is a dimensionless, which is independent of the system of units used for pressure, radial radius and time. (2) The essence of the double logarithmic chart analysis commonly used in oil and gas reservoir engineering is the elastic analysis of pressure with respect to time ( log P̃log t).
According to Eq. (2), the elastic outer boundary condition of the pressure difference P with the outer boundary Γ ∶ r = R is defined as follows where, P Γ is the elasticity coefficient with the outer boundary Γ ∶ r = R and a dimensionless quantity. p 0 is the original reservoir pressure, and p is the reservoir pressure.
According to Eq. (3), it is easy to obtain that those three kinds of outer boundary conditions are special cases of elastic outer boundary conditions: (1) When the outer boundary is closed, P r = P r (r, t) = ln P ln r = r P ⋅ P r , P t = P t (r, t) = ln P ln t = t P ⋅ P t . (3) When the outer boundary is constant pressure, P D | |r D =R D = 0 corresponds to P D Γ → +∞. (3) When the outer boundary is closed, P D | |r D =+∞ = 0 corresponds to R D → +∞.
Therefore, the elastic outer boundary condition is the extension of closed, constant pressure and infinite outer boundary, which can describe the actual state of the reservoir more truly, and has more general meanings and generalities.
In recent three years, the elastic outer boundary conditions have been applied to the reservoir seepage model [3,4], replacing the three special outer boundary conditions commonly used before. The results show that the elastic outer boundary conditions can better reflect the real situation of the reservoir than the special outer boundary conditions. At the same time, the similar construction method used to solve the reservoir seepage model has also been greatly developed [3][4][5][6][7], and the simple algebraic method used to solve the seepage model has provided great convenience for scholars.
Nowadays, shale gas is a kind of biogenic or thermogenic unconventional natural gas that is preserved in organic-rich shale in absorption or free gas status. Shale gas resources are very abundant all over the world and have a greater potential of development than conventional gas. As a new filed of global oil and gas exploration and development, shale gas is attracting more and more attention.In 2012, Azom P. N. et al. [8] expressed the matrix -fracture interflow coefficient as a function of pseudo pressure, and established a dual media seepage model for naturally fractured shale gas reservoirs considering slippage effect and Knudsen diffusion. In 2012, Guo et al. [9] studied double-permeable seepage in horizontal well production of fractured-vuggy carbonate reservoirs for the first time, and the flow curves are analyzed deeply, but their analysis was just based on the three ideal outer boundary conditions. In 2014, Yu et al. [10] established and solved a shale gas reservoir mathematical model considering the desorption and adsorption process, they developed a single well productivity formula and bottom hole pressure formula for shale gas reservoirs. In 2015, Guo et al. [11] based on the theory of trilinear seepage in horizontal wells fractured in shale gas reservoirs to build a mathematical model considering the slip seepage, and they obtained the productivity equation of fractured horizontal well that can be used for predicting on-site production. In 2015, Ma et al. [12] established a trilinear seepage model for the finite conductivity of vertical fractures in shale gas | https://doi.org/10.1007/s42452-021-04787-y Research Article reservoirs by considering the dispersion and desorption of shale gas in the matrix, they obtained the model's solution by using the Laplace transform and Stehfest numerical inversion, and they also studied the influences of sensitive parameters on the dimensionless wellbore pressure to draw corresponding conclusions. In 2018, Zhao et al. [13] proposed a comprehensive mathematical model of the percolation mechanism of multi-stage fractured horizontal wells and studied the influences of sensitive parameters such as fracture number and formation permeability on the pressure profile. In 2019, Huang et al. [14] established a seepage model of horizontal well in double-medium shale gas reservoir with the infinite outer boundary, then they drew the diagrams of the seepage pressure curves of horizontal well in shale gas reservoir, they also analyzed the influences of sensitive parameters(such as reservoir and seepage) on the variation characteristics of horizontal well pressure. In 2019, Liu et al. [15] established a semianalytical model for finite conductivity multi-wing fractured wells in a dual-permeability reservoir at constant rates in the Laplace space, then they verified the solution of their model under the condition of dual-hole and dualpermeability hydraulic fractures, and they also analyzed the influences of important parameters on dimensionless pressure and its derivative in detail. It can be seen that although many previous scholars have established different shale gas reservoir seepage models basing on different factors from different perspectives and the scholars have solved the models with different methods. But all these models have one thing in common and that is the models established are all based on three special outer boundary conditions. As a result, the models can not be used to describe the real state of shale gas reservoir accurately.
Based on the above reasons, the seepage model of shale gas reservoir with dual media is established based on elastic outer boundary condition from three different perspectives( 1. matrix pore slippage effect is only considered; 2. considering matrix micropore slippage effect and gas diffusion factor; 3. considering the matrix stress sensitive, slippage effect and micropore gas diffusion) in the second part, and the model is solved by using the similar construction method. In the third part, we draw curves which can show the flow characteristic of bottom-hole pressure and analyze the influence of related parameters on the curves of flow characteristic. In the fourth part, we discuss the results from two perspectives of calculation results and graphics drawing, and illustrate the importance of elastic outer boundary conditions. Finally, the corresponding conclusions are drawn.

The establishment of dual-media seepage model for shale gas reservoirs under the elastic outer boundary condition and its solution
Assumptions of the model: (1) The reservoir is of homogeneous, horizontal, isopachous and anisotropic; (2) The gas seepage is single-phase, gravity and capillary force are also ignored; (3) The temperature at each point in the shale gas reservoir remains unchanged, that is the seepage process is of isothermal; (4) Stress sensitivity of fractures are not considered; (5) The skin factor of matrix and fracture are the same; (6) Shale gas formation is slightly compressible, and its micro-compressibility coefficient is a constant; (7) Fractures are hydraulic fractures, and production is determined by the shale gas wells; (8) The matrix cannot supply shale gas directly to the wellbore, and the matrix can produce quasi-steady state interflow to the fracture ( Fig. 1). (9) The Langmuir isothermal adsorption equation is obeyed.
In view of the parse complex problems (such as adsorption, diffusion, slip and seepage etc.), we consider the following three angles respectively: (1) matrix pore slippage effect is only considered; (2) considering matrix micropore slippage effect and gas diffusion factor; (3) considering the matrix stress sensitive, slippage effect and micropore gas diffusion, the shale gas seepage Eqs. (54) and (55) are established on the basis of considering desorption and adsorption. At the same time, combining the initial condition (56) with the inner boundary conditions (57) and the elastic outer boundary condition (58), the shale gas seepage model is formed as follows (modeling process details are included in Appendix A): The established mathematical model of seepage can be turned into the mathematical model of dimensionless form as follows: Taking the Laplace transform of the above dimensionless mathematical model dimensionless time t D : The following can be obtained: By using the similarity construction method,the structure of the solution is as follows (derivation details are included in Appendix B): Substituting p D1 or p D1 into Eq. (9), let r D = 1 , the dimensionless bottom hole pressure in the Laplace space can be obtained: where m,n r D , R D , i , (m, n = 0, 1) is a linearly independent combination of the first and second order v-variant

Characteristic curves and sensitivity analysis
The model's solution in the real space can be obtained by the Stehfest numerical inversion [16], and all log-log type curves are drawn with dimensional variables in the real domain. At the same time, to understand the effect of the parameters on well performance and obtain more accurate results during well testing, the effect factors of characteristic curves are analyzed. In this section, five kinds of relevant parameters (including elasticity coefficient, skin factor, wellbore storage coefficient, outer boundary radius, permeability ratio) are analyzed. In addition, in the following figure, the following figure the solid line is the bottom hole pressure curve, and the dotted line is the derivative of the bottom hole pressure.
In the first stage, the shale gas from the fracture system and the pore medium flows into the wellbore simultaneously. At this time, the bottom hole pressure reflects the characteristics of homogeneous gas reservoirs, so those four bottom hole pressures and its derivative's curves coincide with each other respectively. The curves of the first stage indicates that the change of the elastic coefficient P Γ at this stage will not affect the characteristic curves of the bottom hole pressure and its derivative. The second stage is the transition stage. There is an obvious difference between the average pressure of pore medium and fracture. Under the action of this difference, the phenomenon of interporosity flow or overflow occurs between the two medias. Figure 1 shows that when the elastic coefficient is p Γ = inf , the bottom hole pressure curve is almost not influenced by channeling. However, and when P Γ satisfies inf → 0 , the bottom hole pressure is affected by interporosity flow for a longer time. In addition, The downward trend of the bottom hole pressure increases with the decrease of the elastic coefficient during the process of interporosity flow, and then slowly decreases when reaching a critical point of P Γ , and finally rises to the bottomhole pressure curve corresponding to P Γ = 0. In the third stage, the bottomhole pressure curve starts to go up normally and steadily increases after the interporosity flow process. As can be seen from in Fig. 1, When the elastic coefficient decreases from inf → 0 , the bottomhole pressure increases gradually, and in this process, the smaller the elastic coefficient, the longer the process takes.
In the fourth stage, the bottom hole pressure curve changes from steep to flat, then to horizontal, and the greater the elasticity coefficient, the earlier the process begins. In the first stage, the characteristic curves of the dimensionless bottom hole pressure and its derivative overlap respectively, and the dimensionless bottom hole pressure shows a steady increase, indicating that the change of skin factor will not affect the dimensionless bottom hole pressure in this stage, but the completion time of this stage is relatively short.

Effect of S on the characteristic curves
In the second stage, Fig. 3 shows that the dimensionless bottom hole pressure slowly rises in a short period of time, and then quickly go up to a normal level, and the smaller skin factor is, the greater the impact on this stage. During this stage, the dimensionless bottom hole pressure in the first half time roses more slowly, and rises more quickly in the second half time. However, no matter what the skin factor value is, the beginning and end times of this stage will not occur significantly. That is to say, the duration of this stage does not change with the changes of skin factor.
In the third stage, the characteristic curves of the dimensionless bottom hole pressure and its derivative overlap respectively, and the dimensionless bottom hole pressure returns to the normal rising trajectory and shows a steady growth, which indicates that the change of skin factor will not affect the dimensionless bottom hole pressure in this stage. But unlike the first stage, it takes longer to complete this stage.
In the fourth stage, the dimensionless bottom hole pressure begins to transition from a uniform growth rate to a straight horizontal line and finally to a horizontal line. If the skin factor S is larger, the transition time is much earlier, and the pressure value corresponding to the final dimensionless bottom hole pressure is smaller.

Effect of C D on the characteristic curves
In Fig. 4, the values of wellbore storage coefficient are C D = 10000, 20000, 40000, 80000 , respectively, while other  In the first stage, the initial position of the bottom hole pressure curve decreases with the increase of C D , the upward trend of the wellbore pressure curve is the same, and it can be seen from Fig. 4 that C D has the greatest influence on this stage.
In the second stage, affected by channeling, the rate of increase of the dimensionless bottom hole pressure slowed down in a short period of time, and then reached a stable rising state. Moreover, in this stage, the smaller the wellbore storage coefficient C D , the more obvious the channeling process, and the greater the variation degree of dimensionless bottom hole pressure will be.
In the third stage, those four dimensionless characteristic curves show the same trend from top to bottom with the increasing of C D .
In the fourth stage, the characteristic curves coincide respectively, and the dimensionless bottom hole pressure begins to transition from a steady increase to a horizontal straight line, and tends to a horizontal line finally, which shows that the change of C D will not affect this stage.

Effect of R D on the characteristic curves
In Fig. 5, the values of the outer boundary radius are R D = 500, 1000, 2000, 4000 , respectively, while other parameters are p Γ = 0 , r D = 1,S = 2 , k = 0.7 , w = 10 −3 , = 10 −4 , C D = 10000. In the first stage, the characteristic curves of the dimensionless bottom hole pressure and its derivative overlap respectively, and the dimensionless bottom hole pressure shows a slow growth trend, indicating that the change of R D will not affect the dimensionless bottom hole pressure in this stage, but the completion time of this stage is relatively short.
In the second stage, Fig. 5 shows that when R D is smaller, the time of interporosity flow process occurs earlier and shorter, and the corresponding dimensionless bottom hole pressure value after the completion of interfacial flow process is smaller.
In the third stage, the dimensionless bottom hole pressure returns to the normal rising trajectory and shows a steady growth after passing through the interporosity flow process. However, according to Fig. 5, the smaller R D is, the longer the stable growth will be.
In the fourth stage, the dimensionless bottom hole pressure gradually transitions from a steady increase to a straight horizontal line and to a horizontal line finally. If R D is smaller, the transition time is later, and the pressure value corresponding to the final dimensionless bottom hole pressure is larger.
In the first stage, the dimensionless bottom hole pressure shows a slow growth trend, and the greater k is, the later this stage starts and ends.
In the second stage, Fig. 6 shows that the smaller k is, the earlier and more obvious the interporosity flow process occurs, and the smaller the dimensionless bottom hole pressure after the completion of the interporosity flow process.  In the third stage, the dimensionless bottom hole pressure returns to the normal rising trajectory and shows a steady growth after passing through the interporosity flow process. However, according to Fig. 6, the greater k is, the longer the stable growth will be.
In the fourth phase, the dimensionless bottomhole pressure begins to grow steadily, gradually, to a horizontal line, and finally to a horizontal line. When k is smaller, the transition time is later, and the final pressure value corresponding to dimensionless bottomhole pressure is larger.

Discussion
According to the calculation results in the second part and the detailed description of elastic outer boundary conditions in the first part, we can get the following three special cases: (1) For the closed outer boundary condition, only let p Γ = 0 , and then the similar kernel function of the solution can be obtained: (2) For the constant pressure outer boundary condition, only let p Γ → ∞ , and then the similar kernel function of the solution can be obtained: (3) For the infinite outer boundary condition, only let R D → ∞ , and then the similar kernel function of the solution can be obtained: It can be seen from the above three conclusions that, in fact, the traditional outer boundary conditions can be regarded as three special cases of the elastic outer boundary conditions, because when the elastic coefficient changes, the structure of the solution does not change, only the similarity kernel function changes, that is, the elastic coefficient is only related to the similarity kernel function, and has nothing to do with the structure of the solution.
Meanwhile, it can be seen from Fig. 2 in the third part that the elastic outer boundary characteristic curves lies between the constant pressure outer boundary and the closed outer boundary characteristic curve.It is shown that the introduction of the elastic outer boundary can not only describe the real state of the reservoir boundary, but also expand the reasonable interpretation range of the measured data, and greatly improve the situation where the reservoir parameters cannot be explained well due to the data deviation.This is also the spirit of the introduction of elastic outer boundaries in this paper.

Conclusions
Although the matrix permeability has different meanings from three different angles, that is, combining those two factors(i.e. only consider the matrix pore slippage effect, the matrix micropore slippage effect and gas diffusion effect) with those three factors(i.e. the matrix micropore stress sensitivity, slippage effect and gas diffusion effect), the expression of the percolation model combined with different matrix permeability is similar, and the model's solution can be written in a unified expression. The similarity construction method is a simple algebraic method to solve complex seepage models. The dual-media seepage model of shale gas reservoir under the elastic outer boundary condition is solved by the similarity construction method, which greatly simplifies the solution process and obtains a unified expression form. The influences of different parameters on dimensionless bottom hole pressure and its derivative can be clearly analyzed, and the optimization idea is provided for the research and development of well test analysis software.
Elasticity is unique to anything, and the outer boundary of oil and gas reservoir is of also elastic. The elastic outer boundary condition is reasonably introduced in the dualmedia seepage model of shale gas reservoir. Based on the original three ideal outer boundary conditions (i.e. closed, constant pressure, infinite), a reasonable interpretation range of the measured data is expanded, and the actual situation of shale gas reservoir can be described better. long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.

Appendix A. (Modeling process details)
According to the definition of desorption and adsorption [17], when a well is opened, the free gas in the fracture first flows to the wellbore and is recovered. After that, as the pressure and concentration of the gas reservoir decrease, the adsorbed gas on the surface of the matrix undergoes desorption and adsorption and the free gas. Into free gas, which conforms to the following Langmuir isotherm adsorption equation.
Where, V is the adsorption capacity, V L and p L are Langmuir volume and pressure, and the p m is the matrix pressure.
Put the desorption adsorption term into the comprehensive compression coefficient C tm , which is defined as follows [18,19]: Where, the expression of desorption adsorption compression coefficient C d is Where, p sc , Z sc and T sc are standard pressure, gas compression factor under standard condition, and standard temperature.
The following first discusses the expression of matrix permeability and the expression of gas flow velocity in nano-scale pore of the matrix when considering several characteristics of shale gas reservoir 1) Only considering slippage effect in in the nano-scale pore matrix: The formula for matrix permeability affected by slippage effect is as follows [20] Where, k m0 is the initial permeability of the matrix, p m is the average pressure of the matrix, b is the slip factor, and (26) C tm = C m + C g + C d . (28) k app1 is the permeability in the nano-pore of the matrix only considering the slip effect. The gas seepage velocity in the nano-scale pore matrix is calculated by Where, g is gas viscosity, v m1 is the gas seepage velocity in the nano-scale pore matrix.
2) In the above case, gas diffusion [21] is also considered: Gas diffusion velocity is calculated by Where v sm is gas diffusion velocity, M g is the molecular weight of the gas, D m is the mass diffusion coefficient, and C m is the standard gas concentration. The seepage velocity in nano-scale pores of matrix can be composed via Darcy velocity v pm and gas diffusion flow velocity v sm : The equations of state in shale matrix can be expressed as follows The gas concentration equation is given by (33) g = p m M g ZR g T .
(34) C m = m M g = p m RTZ . (35) (36) v sm = −D m C gm ⋅ p m r . 3) Based on the above discussion, stress sensitivity of nano-pore [22], slippage effect and gas diffusion are considered comprehensively. The relationship between fracture and rock permeability and pressure is as follows Where, p 0 is the initial pressure, k i0 is the initial permeability and i is the permeability modulus.
Therefore, the combined Eq.(e37) and Eq.(e39), considering stress sensitivity, slippage effect and gas diffusion are used to express the gas seepage velocity in the modified nanoscale pores of the matrix as follows: In shale reservoirs, matrix permeability, slippage factor is defined as follows Where k app3 is Matrix permeability, b bm is slippage factor. Substituting Eqs. (41) and (42) into Eq. (40), the expression of gas seepage velocity in the nano porosity of the matrix can be obtained: Next, the following model is established.
Continuous equation of fracture system: