# Iterative Analytical Solutions for Nonlinear Two-Phase Flow with Gas Solubility in Shale Gas Reservoirs

Authors: Xiaoji Shang; J. G. Wang; Zhizhen Zhang
Journal: Geofluids
 (2019)

---

## Abstract

The governing equations of a two-phase flow have a strong nonlinear term due to the interactions between gas and water such as capillary pressure, water saturation, and gas solubility. This nonlinearity is usually ignored or approximated in order to obtain analytical solutions. The impact of such ignorance on the accuracy of solutions has not been clear so far. This study seeks analytical solutions without ignoring this nonlinear term. Firstly, a nonlinear mathematical model is developed for the two-phase flow of gas and water during shale gas production. This model also considers the effects of gas solubility in water. Then, iterative analytical solutions for pore pressures and production rates of gas and water are derived by the combination of travelling wave and variational iteration methods. Thirdly, the convergence and accuracy of the solutions are checked through history matching of two sets of gas production data: a China shale gas reservoir and a horizontal Barnett shale well. Finally, the effects of the nonlinear term, shale gas solubility, and entry capillary pressure on the shale gas production rate are investigated. It is found that these iterative analytical solutions can be convergent within 2-3 iterations. The solutions can well describe the production rates of both gas and water. The nonlinear term can significantly affect the forecast of shale gas production in both the short term and the long term. Entry capillary pressure and shale gas solubility in water can also affect shale gas production rates of shale gas and water. These analytical solutions can be used for the fast calculation of the production rates of both shale gas and water in the two-phase flow stage.


---

## Body

## 1. Introduction

Hydraulic fracturing is a key technology to the development of shale gas reservoirs with ultralow permeability [1]. After hydraulic fracturing, there are fracturing fluid, initial water, and gas in the gas reservoir [2, 3]. The gas flow in the shale gas reservoir usually experiences three stages as shown in Figure 1: the single-phase water flow stage, the unsaturated bubble flow stage, and the gas-water two-phase flow stage [4]. At the last two stages, both gas and water simultaneously flow to the wellbore. The two-phase flow in these stages has significant influence on the gas production life because of complex interactions between gas and water [5]. Therefore, transport behaviors and mechanisms should be carefully studied for gas production forecast.Figure 1
Process of shale gas production.Numerous numerical models have been proposed to describe the complex process of the two-phase flow in porous media [6–8]. For instance, Mahabadi et al. [6] combined the techniques of microfocus X-ray computed tomography (CT) and pore-network model simulation to identify key parameters and their values in characterizing immiscible gas-water two-phase flow during hydrate dissociation. For the safety assessment on the CO2 geological storage, Wang and Peng [7] proposed a fully coupled mathematical model to describe the combined effects of the two-phase flow, deformation, and carbon dioxide (CO2) sorption on caprock sealing efficiency. Also, Wang et al. [8] developed a two-phase numerical model with multiscale diffusion mechanisms for both the flowback and long-term production stages. All these numerical models can efficiently analyze the process of the two-phase flow. However, there exist two unsolved issues. One issue is to include the gas solubility in water. Some experiments observed that the organic matter in water, high temperature, and pressure can lead to the increase of the methane solubility in water [9–11]. How this gas solubility affects the gas production rate has not been studied. The other issue is to obtain analytical solutions of water and gas pressures if all nonlinear interactions in the two-phase flow are included. Analytical solutions can forecast the field production data in a more convenient and faster way.Considerable analytical solutions have been proposed for the single-phase flow in oil and gas reservoirs [12–15]. For example, Nobakht and Clarkson [14] developed an analytical method to analyze the linear flow when the constant gas production rate is considered. They got a linear relationship between the square-root-of-time and production rate. In practice, the flow pattern may not be simplified into a linear flow or single-phase flow. Further improvements were made by introducing some factors related to the two-phase flow into the governing equations of the single-phase flow [16, 17]. For example, Behmanesh et al. [17] introduced two-phase viscosity and two-phase compressibility into the single-phase flow. They established a semianalytical solution of this single-phase flow for rate-transient analysis of gas condensate reservoirs. However, they did not directly solve the governing equations of the two-phase flow with a nonlinear term.Several analytical solutions have been obtained for the two-phase flow after some simplifications. For example, Wang et al. [18] developed a semianalytical solution for the two-phase flow after the linearization of the mathematical model. This solution can quickly predict production performance. In order to estimate the residual saturation of the displaced fluid, Adibifard [19] decomposed the total pressure drop into Poiseuille and Young-Laplace terms; thus, an analytical solution was obtained for a two-phase flow in a capillary tube. Yang et al. [20] investigated the effects of multifracture networks on the gas production rate and developed a semianalytical solution for the two-phase flow during the flowback period with complex fracture networks. They found that increasing the fracture network complexity is favorable to gas production enhancement. For calculation convenience, these abovementioned analytical solutions of two-phase flow models are all obtained after omitting or linearizing the nonlinear term. The nonlinear term in governing equations comes from the interaction between gas and water such as capillary pressure, water saturation, and gas solution. Thus, ignoring the nonlinear term may induce some error in the forecast of field gas production.A partial differential equation for a physical process can be transformed into an ordinary differential equation through the travelling wave method, a variable substitution method. In this method, the differential equation is firstly transformed into the variable integral form. Then, a general solution is obtained from this integral form. Finally, definite solution conditions are applied to obtain the particular solution. The travelling wave method has been successfully applied to seek the analytical solutions of seepage problems [21]. It has also been applied to the coupling problem of temperature change with the gas-water two-phase flow in the thermal stimulation-enhanced coalbed methane recovery [22]. The analytical solutions were obtained for a fractal-hydrothermal model. Moreover, the variational iteration method is a useful tool in analytically solving a nonlinear flow problem. This method has the flexibility and ability in seeking analytical solutions of nonlinear equations through iteration. Generally, a nonlinear flow problem can be solved with a highly accurate solution or a series of exact solutions with fast convergence in the variational iteration method. He [23] successfully solved the Duffing equation and fractional partial differential equations of percolation through the variational iteration method. After converting differential equations into a dimensionless form, Kamran et al. [24] used the variational iteration method to solve the non-Newtonian third-grade fluid flow problem between two parallel plates. However, the above-solved problems are all one-phase fluid flow problems. No two-phase flow problem in porous media has been solved with this variational iteration method so far.This study develops iterative analytical solutions of a nonlinear two-phase flow after considering the influence of shale gas solubility in water. Firstly, a nonlinear mathematical model of pore pressures is developed based on the two-phase flow interactions during shale gas production. Secondly, this nonlinear mathematical model is transformed through the travelling wave method. The transformed nonlinear equations on water and gas pressures are analytically solved by the variational iteration method. The analytical solutions of water and gas pressures and the water and gas production rate with time are obtained. The production rate of gas or water or both is validated by the field data of both the flowback stage and long-term production with a China shale gas reservoir [20] and Barnett shale well 314 [25]. Finally, the effects of the nonlinear term in the two-phase flow, shale gas solubility in water, and entry capillary pressure on the shale gas production rate are investigated. Based on these studies, conclusions are drawn.


## 2. Governing Equations for Shale Gas and the Water Two-Phase Flow

The pores in a fractured shale gas reservoir are assumed to be fully filled with gas or water or both. The two-phase flow in the fractured shale gas reservoir is described by following governing equations and the constitutive model.


### 2.1. Capillary Pressure Model

The capillary pressure,pc, is defined as the pressure difference between the gas pressure and the water pressure [7, 26]:
(1)pc=pg−pw.The pores are fully saturated by gas and water:(2)sw+sg=1.A normalized saturation of the water phasesw∗ is defined as
(3)sw=sw∗1−srg−srw+srw,where pg is the pore pressure of gas and pw is the pore pressure of the water phase. sw and sg are the saturations of water and gas, respectively. srw denotes the irreducible saturation of water and srg denotes the irreducible saturation of gas.The normalized water-phase saturation,sw∗, is a function of capillary pressure as [7, 27]
(4)sw∗=pepcλ,where pe is the entry capillary pressure and λ is the pore size distribution index.Integrating equation (4) into equation (3) yields
(5)sw=1−srg−srwpepcλ+srw.Therefore, the derivative of water saturation with respect to time is(6)∂sw∂t=1−srg−srwpeλ⋅∂pc−λ∂t.Similarly, for the gas phase,(7)sg=1−srw−1−srg−srwpepcλ,(8)∂sg∂t=−1−srg−srwpeλ∂pc−λ∂t.


### 2.2. Governing Equations of the Two-Phase Flow

In the two-phase flow, the continuity equation is expressed for each phase. The continuity equation for gas is(9)∂∂tϕsgρg+∇⋅ρgVg+RswρgVw=ρgqg.The continuity equation for water is(10)∂∂tϕρwsw+∇⋅ρwVw=ρwqw,where ϕ is the porosity and Rsw is the gas solubility in water. ρg is the density of gas and ρw is the density of water under formation conditions. Vg represents the velocity of gas and Vw represents the velocity of water. qg and qw are the gas and water sources, respectively. ∇ is the Hamiltonian operator.Water is slightly compressible but its densityρw can still be regarded as a constant. The gas density ρg follows the equation of the state:
(11)ρg=MZRTpg,where M is the molecular weight of gas, Z is the gas compressibility factor, R is the universal gas constant, and T is the gas temperature.Integrating equation (11) into equation (9) yields
(12)∂∂tϕsgpg+∇⋅pgVg+RswpgVw=pgqg.Hence, equation (10) is simplified into
(13)∂∂tϕsw+∇⋅Vw=qw.BothVg and Vw in porous media are described by the Darcy law as [7]
(14)Vg=−kkrgμg∇pg+ρgg∇H,(15)Vw=−kkrwμw∇pw+ρwg∇H,where k is the absolute permeability of porous media, krg denotes the relative permeability of gas, and krw denotes the relative permeability of water. μg and μw are the viscosities of gas and water, respectively. g is the gravity acceleration and H refers to the vertical elevation.Ifqw and qg are all zeros and the second-order phase is ignored, the governing equations are obtained from equations (12)-(15).For the gas phase,(16)ϕ∂sgpg∂t−∇kkrgpgμg∇pg+Rswkkrwpgμw∇pw=0.For the water phase,(17)ϕ∂sw∂t−kμw∇krw∇pw=0.The following constitutive laws are used for relative permeabilitieskrw and krg [28, 29]:
(18)krw=krwmaxsw∗Nw,(19)krg=krgmaxsg∗Ng,where sw∗+sg∗=1 and sg∗=1−pe/pcλ.Integrating equation (4) into equation (18) yields
(20)krw=krwmaxpeλNw⋅pc−λNw,(21)krg=krgmax1−pepcλNg.Integrating equations (1), (6), and (20) into equation (17) gets the governing equation for the water flow.
(22)ϕ1−srg−srwpeλ∂pg−pw−λ∂t−kkrwmaxpeλNwμw∇pg−pw−λNw∇pw=0.Similarly, integrating equations (1), (8), (21), and (16) obtains the governing equation for the gas flow as
(23)ϕ1−srw−1−srg−srwpeλ⋅pg−pw−λ∂pg∂t−ϕ1−srg−srwpeλ⋅pg∂pg−pw−λ∂t−kkrgmaxμg∇1−pepg−pwλNg⋅pg⋅∇pg−RswkkrwmaxpeλNwμw∇pg−pw−λNwpg∇pw=0.


## 2.1. Capillary Pressure Model

The capillary pressure,pc, is defined as the pressure difference between the gas pressure and the water pressure [7, 26]:
(1)pc=pg−pw.The pores are fully saturated by gas and water:(2)sw+sg=1.A normalized saturation of the water phasesw∗ is defined as
(3)sw=sw∗1−srg−srw+srw,where pg is the pore pressure of gas and pw is the pore pressure of the water phase. sw and sg are the saturations of water and gas, respectively. srw denotes the irreducible saturation of water and srg denotes the irreducible saturation of gas.The normalized water-phase saturation,sw∗, is a function of capillary pressure as [7, 27]
(4)sw∗=pepcλ,where pe is the entry capillary pressure and λ is the pore size distribution index.Integrating equation (4) into equation (3) yields
(5)sw=1−srg−srwpepcλ+srw.Therefore, the derivative of water saturation with respect to time is(6)∂sw∂t=1−srg−srwpeλ⋅∂pc−λ∂t.Similarly, for the gas phase,(7)sg=1−srw−1−srg−srwpepcλ,(8)∂sg∂t=−1−srg−srwpeλ∂pc−λ∂t.


## 2.2. Governing Equations of the Two-Phase Flow

In the two-phase flow, the continuity equation is expressed for each phase. The continuity equation for gas is(9)∂∂tϕsgρg+∇⋅ρgVg+RswρgVw=ρgqg.The continuity equation for water is(10)∂∂tϕρwsw+∇⋅ρwVw=ρwqw,where ϕ is the porosity and Rsw is the gas solubility in water. ρg is the density of gas and ρw is the density of water under formation conditions. Vg represents the velocity of gas and Vw represents the velocity of water. qg and qw are the gas and water sources, respectively. ∇ is the Hamiltonian operator.Water is slightly compressible but its densityρw can still be regarded as a constant. The gas density ρg follows the equation of the state:
(11)ρg=MZRTpg,where M is the molecular weight of gas, Z is the gas compressibility factor, R is the universal gas constant, and T is the gas temperature.Integrating equation (11) into equation (9) yields
(12)∂∂tϕsgpg+∇⋅pgVg+RswpgVw=pgqg.Hence, equation (10) is simplified into
(13)∂∂tϕsw+∇⋅Vw=qw.BothVg and Vw in porous media are described by the Darcy law as [7]
(14)Vg=−kkrgμg∇pg+ρgg∇H,(15)Vw=−kkrwμw∇pw+ρwg∇H,where k is the absolute permeability of porous media, krg denotes the relative permeability of gas, and krw denotes the relative permeability of water. μg and μw are the viscosities of gas and water, respectively. g is the gravity acceleration and H refers to the vertical elevation.Ifqw and qg are all zeros and the second-order phase is ignored, the governing equations are obtained from equations (12)-(15).For the gas phase,(16)ϕ∂sgpg∂t−∇kkrgpgμg∇pg+Rswkkrwpgμw∇pw=0.For the water phase,(17)ϕ∂sw∂t−kμw∇krw∇pw=0.The following constitutive laws are used for relative permeabilitieskrw and krg [28, 29]:
(18)krw=krwmaxsw∗Nw,(19)krg=krgmaxsg∗Ng,where sw∗+sg∗=1 and sg∗=1−pe/pcλ.Integrating equation (4) into equation (18) yields
(20)krw=krwmaxpeλNw⋅pc−λNw,(21)krg=krgmax1−pepcλNg.Integrating equations (1), (6), and (20) into equation (17) gets the governing equation for the water flow.
(22)ϕ1−srg−srwpeλ∂pg−pw−λ∂t−kkrwmaxpeλNwμw∇pg−pw−λNw∇pw=0.Similarly, integrating equations (1), (8), (21), and (16) obtains the governing equation for the gas flow as
(23)ϕ1−srw−1−srg−srwpeλ⋅pg−pw−λ∂pg∂t−ϕ1−srg−srwpeλ⋅pg∂pg−pw−λ∂t−kkrgmaxμg∇1−pepg−pwλNg⋅pg⋅∇pg−RswkkrwmaxpeλNwμw∇pg−pw−λNwpg∇pw=0.


## 3. Analytical Solutions

### 3.1. Travelling Wave Transform

In order to obtain the analytical solution for the gas flow more easily, the travelling wave method is applied to get the ordinary differential equation. The travelling wave method works with the following procedure. When a given time and particular location are concerned, where the impact of the border is too late to reach, the wave always moves forward and forms a travelling wave [30–32]. A general form of nonlinear equations is
(24)Hx,t,u,ux,ut,uxx,uxt,utt,⋯=0,where x is a variable to represent the space coordinate and t represents the time coordinate. u is the function of x and t, while H is a suitable function of u and its derivative.If the solutionϕξ of equation (24) just depends on x and t in the form ξ=x+ct, among which c is a constant to represent the wave velocity, ϕξ is the travelling wave solution.Let(25)ξ=x+y+z+ct,then
(26)pgξ=pgx+y+z+ct=pgx,y,z,t,pwξ=pwx+y+z+ct=pwx,y,z,t,∂pg∂t=c∂pg∂ξ,∂pw∂t=c∂pw∂ξ,∇=∂∂xi→+∂∂yj→+∂∂zk→=∂∂ξn→,∇2=∂2∂x2+∂2∂y2+∂2∂z2=∂2∂ξ2.The above procedure is applied to the partial differential equations of equations (22) and (23). At this time, equation (22) is transformed into
(27)cϕ1−srg−srwpeλ∂pg−pw−λ∂ξ−kkrwmaxpeλNwμw∂∂ξpg−pw−λNw⋅∂pw∂ξ=0.Equation (27) can be further simplified into
(28)pg=pw+Apw′1/λNw−1,where
(29)A=kkrwmaxpeλNw−1cμwϕ1−srg−srw,pw′=∂pw∂ξ.Equation (23) is transformed into
(30)ϕ1−srw−1−srg−srwpeλ⋅pg−pw−λ∂pg∂ξ−ϕ1−srg−srwpeλ⋅pg∂pg−pw−λ∂ξ−kkrgmaxμg⋅∂∂ξ1−pepg−pwλNg⋅pg⋅∂pg∂ξ−RswkkrwmaxpeλNwμw⋅∂∂ξpg−pw−λNw⋅pg⋅∂pw∂ξ=0.For the convenience of calculation,uξ is used as a substitute of pwξ as
(31)pwξ=uξ.Integrating equations (28) and (31) into equation (30) yields
(32)ϕc1−srwu′ξ−ϕc1−srg−srwpeλA1/1−Nwu′ξ2−Nw/1−Nw+fuξ,u′ξ,u″ξ=0,where fuξ,u′ξ,u″ξ is given by equation (55) in Appendix.This is our nonlinear ordinary differential equation for the water flow and is analytically solved by the variational iteration method.


### 3.2. Variational Iteration Method

The variational iteration method [23, 24] is applied to the nonlinear problem of equation (32) to get an iterative analytical solution. The basic procedure for the variational iteration method is briefly described below. A general partial differential equation is
(33)Lξu+Nu=gξ,where Lξ is a linear operator of ξ and N is a nonlinear operator of u.A correction functional is constructed as(34)un+1ξ=unξ+∫0ξγLξun+Nun−gds,where γ is a Lagrange multiplier to be identified by the variational theory. un is a restricted variational [33] under δun=0.Therefore, the correction functional for equation (32) is
(35)un+1ξ=unξ+∫0ξγLξunω+Nunωdω.The Corey parameter for water and the pore size distribution index can be taken from publications [7, 8]. In this study, Nw=2 and λ=2; thus, the linear operator in equation (32) is
(36)Lξunω=−ϕc1−srg−srwpe2A−1−kkrgmaxμg+Rswkkrwmaxμwpe4A−2+2kkrgmaxpe2A−1μgunω+ϕc1−srwun′ω.The nonlinear operatorNunω is given by equation (56) in Appendix.Making the above correction functional of equation (35) stationary gets
(37)δun+1ξ=δunξ+δ∫0ξγLξunω+Nunωdω=0.The nonlinear term in equation (37) is considered as a restricted variational which is
(38)δNunω=0.Equation (37) can be translated into
(39)δun+1ξ=δunξ+δ∫0ξγ2kkrgmaxpe2A−1μgunω+ϕc1−srwun′ω−ϕc1−srg−srwpe2A−1−kkrgmaxμg+Rswkkrwmaxμwpe4A−2+Nunωdω=0,where
(40)δ−ϕc1−srg−srwpe2A−1−kkrgmaxμg+Rswkkrwmaxμwpe4A−2=0.Substituting equations (38) and (40) into equation (39) yields
(41)δun+1ξ=δunξ+δ∫0ξγ2kkrgmaxpe2A−1μgunω+ϕc1−srwun′ωdω=0.Integrating equation (29) into equation (41) yields
(42)δun+1ξ=δunξ+δ∫0ξγ2ϕcμwkrgmax1−srg−srwμgkrwmaxunω+ϕc1−srwun′ωdω=δunξ+ϕc1−srwγωδunωω=ξ+∫0ξ−ϕc1−srwγ′ω+2ϕcμwkrgmax1−srg−srwμgkrwmaxγωδundω=0.The following stationary conditions are then obtained:(43)1+ϕC1−srwγωω=ξ=0,−ϕC1−srwγ′ω+2ϕCμwkrgmax1−srg−srwμgkrwmaxγω=0.The general solution of the Lagrange multiplierγω is
(44)γω=Ce2μwkrgmax1−srg−srw/μgkrwmax1−srwω,where C is a constant to be determined by the first formula of equation (43) as
(45)C=−1ϕc1−srwe−2μwkrgmax1−srg−srw/μgkrwmax1−srwξ.Finally, the multiplier is expressed as(46)γ=−1ϕc1−srwe2μwkrgmax1−srg−srw/μgkrwmax1−srwω−ξ.The iteration solution in theξ direction has the following form:
(47)un+1ξ=unξ−∫0ξe2μwkrgmax1−srg−srw/μgkrwmax1−srwω−ξgnωdω,where gnω is given by equation (57) in Appendix.


### 3.3. Final Solution of Water-Phase Pressure

The final analytical solution is obtained as follows.For the water-phase pressure,(48)pwξ=u1ξ.For the gas-phase pressure,(49)pg=u1ξ+Au1′ξ1/λNw−1.


### 3.4. Analytical Solution of the Gas Production Rate

The shale gas production rate is defined as [34]
(50)dGptdt=−∫ϕpadpgdtdv,where pa=101.325kPa.Integrating equation (49) into equation (50) yields the final form of the gas production rate as
(51)dGpgtdt=−∫Ωϕpadu1t+Au1′t1/λNw−1dtdΩ.


### 3.5. Analytical Solution of the Water Production Rate

Based on Darcy’s law, the relationship between shale gas production and water production is obtained by [35]
(52)GpgGpw=krgkrwΔpgΔpwμwμg.Integrating equation (51) into equation (52) yields the final form of the water production rate as
(53)dGpwtdt=−krwkrg⋅μgμw⋅ΔpwΔpg⋅∫Ωϕpadu1t+Au1′t1/λNw−1dtdΩ.


## 3.1. Travelling Wave Transform

In order to obtain the analytical solution for the gas flow more easily, the travelling wave method is applied to get the ordinary differential equation. The travelling wave method works with the following procedure. When a given time and particular location are concerned, where the impact of the border is too late to reach, the wave always moves forward and forms a travelling wave [30–32]. A general form of nonlinear equations is
(24)Hx,t,u,ux,ut,uxx,uxt,utt,⋯=0,where x is a variable to represent the space coordinate and t represents the time coordinate. u is the function of x and t, while H is a suitable function of u and its derivative.If the solutionϕξ of equation (24) just depends on x and t in the form ξ=x+ct, among which c is a constant to represent the wave velocity, ϕξ is the travelling wave solution.Let(25)ξ=x+y+z+ct,then
(26)pgξ=pgx+y+z+ct=pgx,y,z,t,pwξ=pwx+y+z+ct=pwx,y,z,t,∂pg∂t=c∂pg∂ξ,∂pw∂t=c∂pw∂ξ,∇=∂∂xi→+∂∂yj→+∂∂zk→=∂∂ξn→,∇2=∂2∂x2+∂2∂y2+∂2∂z2=∂2∂ξ2.The above procedure is applied to the partial differential equations of equations (22) and (23). At this time, equation (22) is transformed into
(27)cϕ1−srg−srwpeλ∂pg−pw−λ∂ξ−kkrwmaxpeλNwμw∂∂ξpg−pw−λNw⋅∂pw∂ξ=0.Equation (27) can be further simplified into
(28)pg=pw+Apw′1/λNw−1,where
(29)A=kkrwmaxpeλNw−1cμwϕ1−srg−srw,pw′=∂pw∂ξ.Equation (23) is transformed into
(30)ϕ1−srw−1−srg−srwpeλ⋅pg−pw−λ∂pg∂ξ−ϕ1−srg−srwpeλ⋅pg∂pg−pw−λ∂ξ−kkrgmaxμg⋅∂∂ξ1−pepg−pwλNg⋅pg⋅∂pg∂ξ−RswkkrwmaxpeλNwμw⋅∂∂ξpg−pw−λNw⋅pg⋅∂pw∂ξ=0.For the convenience of calculation,uξ is used as a substitute of pwξ as
(31)pwξ=uξ.Integrating equations (28) and (31) into equation (30) yields
(32)ϕc1−srwu′ξ−ϕc1−srg−srwpeλA1/1−Nwu′ξ2−Nw/1−Nw+fuξ,u′ξ,u″ξ=0,where fuξ,u′ξ,u″ξ is given by equation (55) in Appendix.This is our nonlinear ordinary differential equation for the water flow and is analytically solved by the variational iteration method.


## 3.2. Variational Iteration Method

The variational iteration method [23, 24] is applied to the nonlinear problem of equation (32) to get an iterative analytical solution. The basic procedure for the variational iteration method is briefly described below. A general partial differential equation is
(33)Lξu+Nu=gξ,where Lξ is a linear operator of ξ and N is a nonlinear operator of u.A correction functional is constructed as(34)un+1ξ=unξ+∫0ξγLξun+Nun−gds,where γ is a Lagrange multiplier to be identified by the variational theory. un is a restricted variational [33] under δun=0.Therefore, the correction functional for equation (32) is
(35)un+1ξ=unξ+∫0ξγLξunω+Nunωdω.The Corey parameter for water and the pore size distribution index can be taken from publications [7, 8]. In this study, Nw=2 and λ=2; thus, the linear operator in equation (32) is
(36)Lξunω=−ϕc1−srg−srwpe2A−1−kkrgmaxμg+Rswkkrwmaxμwpe4A−2+2kkrgmaxpe2A−1μgunω+ϕc1−srwun′ω.The nonlinear operatorNunω is given by equation (56) in Appendix.Making the above correction functional of equation (35) stationary gets
(37)δun+1ξ=δunξ+δ∫0ξγLξunω+Nunωdω=0.The nonlinear term in equation (37) is considered as a restricted variational which is
(38)δNunω=0.Equation (37) can be translated into
(39)δun+1ξ=δunξ+δ∫0ξγ2kkrgmaxpe2A−1μgunω+ϕc1−srwun′ω−ϕc1−srg−srwpe2A−1−kkrgmaxμg+Rswkkrwmaxμwpe4A−2+Nunωdω=0,where
(40)δ−ϕc1−srg−srwpe2A−1−kkrgmaxμg+Rswkkrwmaxμwpe4A−2=0.Substituting equations (38) and (40) into equation (39) yields
(41)δun+1ξ=δunξ+δ∫0ξγ2kkrgmaxpe2A−1μgunω+ϕc1−srwun′ωdω=0.Integrating equation (29) into equation (41) yields
(42)δun+1ξ=δunξ+δ∫0ξγ2ϕcμwkrgmax1−srg−srwμgkrwmaxunω+ϕc1−srwun′ωdω=δunξ+ϕc1−srwγωδunωω=ξ+∫0ξ−ϕc1−srwγ′ω+2ϕcμwkrgmax1−srg−srwμgkrwmaxγωδundω=0.The following stationary conditions are then obtained:(43)1+ϕC1−srwγωω=ξ=0,−ϕC1−srwγ′ω+2ϕCμwkrgmax1−srg−srwμgkrwmaxγω=0.The general solution of the Lagrange multiplierγω is
(44)γω=Ce2μwkrgmax1−srg−srw/μgkrwmax1−srwω,where C is a constant to be determined by the first formula of equation (43) as
(45)C=−1ϕc1−srwe−2μwkrgmax1−srg−srw/μgkrwmax1−srwξ.Finally, the multiplier is expressed as(46)γ=−1ϕc1−srwe2μwkrgmax1−srg−srw/μgkrwmax1−srwω−ξ.The iteration solution in theξ direction has the following form:
(47)un+1ξ=unξ−∫0ξe2μwkrgmax1−srg−srw/μgkrwmax1−srwω−ξgnωdω,where gnω is given by equation (57) in Appendix.


## 3.3. Final Solution of Water-Phase Pressure

The final analytical solution is obtained as follows.For the water-phase pressure,(48)pwξ=u1ξ.For the gas-phase pressure,(49)pg=u1ξ+Au1′ξ1/λNw−1.


## 3.4. Analytical Solution of the Gas Production Rate

The shale gas production rate is defined as [34]
(50)dGptdt=−∫ϕpadpgdtdv,where pa=101.325kPa.Integrating equation (49) into equation (50) yields the final form of the gas production rate as
(51)dGpgtdt=−∫Ωϕpadu1t+Au1′t1/λNw−1dtdΩ.


## 3.5. Analytical Solution of the Water Production Rate

Based on Darcy’s law, the relationship between shale gas production and water production is obtained by [35]
(52)GpgGpw=krgkrwΔpgΔpwμwμg.Integrating equation (51) into equation (52) yields the final form of the water production rate as
(53)dGpwtdt=−krwkrg⋅μgμw⋅ΔpwΔpg⋅∫Ωϕpadu1t+Au1′t1/λNw−1dtdΩ.


## 4. Model Verifications

The above two-phase flow model describes the interaction of water and gas in the shale gas production process, including the two-phase flow, capillary pressure, relative permeability, and solubility of gas. The convergence of the iterative analytical solutions is studied here. These iterative analytical solutions are also validated by following two sets of field data. One set is for the shale gas production at the flowback stage and the other set is for the long-term shale gas production.


### 4.1. Convergence of Analytical Solution

The convergence of iterative analytical solutions is checked here. The following mean relative tolerance is defined for this convergence study:(54)δ=dGpt/dtn+1−dGpt/dtndGpt/dtn+1×100%.The mean relative tolerance with the number of iterations is presented in Figure2. The mean relative tolerance is 0.241% for the first iteration, 0.139% for the second iteration, and 0.087% for the 3rd iteration. This indicates that the analytical solutions of the shale gas production rate are convergent with the number of iterations.Figure 2
Convergence check for the iterative analytical solution.


### 4.2. Model Verifications of Water and Gas Production Rates at the Flowback Stage

The early time-pressure-rate data immediately after hydraulic fracturing were collected as flowback data [36]. In the flowback stage, a set of field data of water and gas production rates was obtained from a China shale gas reservoir [20]. Table 1 lists the computational parameters used in analytical solutions. These parameters are taken from Wang and Peng [7] and Yang et al. [20]. Figure 3(a) compares the shale gas production rates among the field production data, the analytical solution of our model, and the semianalytical solution of Yang et al.’s model. Figure 3(b) presents the normalized mean square errors (NMSE) between our analytical solution of shale gas production rates and the field data. According to the fitting results, the NMSE is between 0.1% and 1%. The analytical solution of our model is in good agreement with the field data. Further, Figure 4(a) compares the water production rates of the field production data (black circles), the analytical solution of our model (black line), and the semianalytical solution of Yang et al.’s model (red double dash). Figure 4(b) presents the NMSE of the water production rate between our analytical solution and the field data. According to the fitting results, the NMSE is between 0.1% and 2%. The analytical solution of our model can well predict the field data. These figures indicate the following understandings.Table 1
Model parameters in the computation for the China shale gas reservoir.
Parameter
Unit
Value
Physical meanings
Srg
0.15
Residual saturation of gas
Srw
0.2
Residual saturation of water
μw
Pa∗s
3.6×10−4
Water viscosity
μg
Pa∗s
2.0×10−5
Shale gas viscosity
c
m/s
5.3×10−6
Travel wave viscosity
λ
2
Pore size distribution index
p0
MPa
27.4
Initial reservoir pressure
p1
MPa
19.67
Well pressure
k0
mD
0.1
Initial permeability in the fractured zone
ϕ0
0.18
Initial porosity
pe
MPa
2
Entry capillary pressure
krwmax
0.004
End-point relative permeability for water
krnwmax
1
End-point relative permeability for shale gas
Nw
2
Reference parameter for water
Ng
2
Reference parameter for shale gas
Rsw
5.42
Shale gas solubilityFigure 3
Model verifications of gas production in the flowback period of a China shale reservoir.
(a)
Gas production data
(b)
Normalized mean square error between analytical solution and field dataFigure 4
Model verifications of water production in the flowback period of a China shale reservoir.
(a)
Water production data
(b)
Normalized mean square error between the analytical solution and field dataFirstly, the production rate of shale gas increases to the maximum and then decreases. The gas production rate of our model goes up fast till the peak rate of1.69×105 m3/d at the 24th day and then goes down gently. On the other hand, the water production rate of our model goes down fast from 22 m3/d to 2.8 m3/d during the first 50 days and then keeps flat from 50 to 200 days. These results show that the analytical solutions are in good agreement with field production rates of both shale gas and water.Secondly, for the shale gas production rate, the peak rate of our model is1.69×105 m3/d at the 24th day, while the semianalytical solution of Yang et al.’s model is 1.91×105 m3/d at the 55th day as shown in Figure 3(a). The shale gas production rate of our model goes up faster than that of Yang et al.’s model in the rising period. At the 100th day, the shale gas production rates are 1.86×105 m3/d and 1.65×105 m3/d for the two cases, respectively. On the contrary, the shale gas production rate of our model goes down slower than that of Yang et al.’s model in the decline period. For the water production rate, the production rate of our model decreases from 21.86 m3/d to 7.35 m3/d but the production rate by Yang et al.’s model decreases from 29.14 m3/d to 6.14 m3/d after the first 24 days (see Figure 4(a)). The decrease percentages are 66.4% and 78.9%, respectively. The water production rate predicted by our model goes down slower than that by Yang et al.’s model in the whole 200 days. These differences between our model and Yang et al.’s model may be induced by the nonlinear term of the two-phase flow and the gas solubility in water. These two factors may significantly impact the production rates of shale gas and water. The nonlinear term stands for the complex interaction between shale gas and water. Meanwhile, the dissolved gas may change into free gas during depressurization, leading to the going up of the shale gas production rate by our model faster than that of Yang et al.’s model in the rising period.


### 4.3. Model Verifications for Long-Term Production

The analytical solution of the shale gas production rate is applied to the field production data forecast for the horizontal Barnett shale well 314 [37]. As the field water production rate is not available in the literature, thus, the model verification of the water production rate is not conducted here. Table 2 lists the model parameters for the Barnett shale. The reservoir parameters used in this analysis are determined after referencing the publications by Wang and Peng [7] and Yu and Sepehrnoori [25]. The travelling wave velocity is also decided by the reservoir length and the production time from theirs. The actual field data of shale gas production are presented as the black circle in Figure 5(a). The black line is obtained by the analytical solution of the shale gas production rate. The line is in reasonable agreement with the field production data for the whole production time. The gas production rate declines about 70% in the first 200 days. At the 200th day, the gas production rate reaches 6.3×104 m3/d and then slowly descends during the later period of production. As the bottom hole pressure is set as 3.45 MPa, the gas production rate also keeps about 1.6×104 m3/d in the production tail. Figure 5(b) shows the normalized mean square errors between our analytical solution of shale production rates and the field data. This NMSE is between 0.6% and 1%. Therefore, our analytical solution is feasible to predict the long-term gas production rate.Table 2
Model parameters in computation for the Barnett shale gas well.
Parameter
Unit
Value
Physical meanings
Srg
0.15
Residual saturation of gas
Srw
0.6
Residual saturation of water
μw
Pa∗s
3.6×10−4
Water viscosity
μg
Pa∗s
2.0×10−5
Shale gas viscosity
c
m/s
1.2×10−6
Travel wave viscosity
λ
2
Pore size distribution index
p0
MPa
20.34
Initial reservoir pressure
p1
MPa
3.45
Well pressure
k0
mD
5.0
Initial permeability in the fractured zone
ϕ0
0.2
Initial porosity
pe
MPa
2
Entry capillary pressure
krwmax
0.00035
End-point relative permeability for water
krnwmax
1
End-point relative permeability for shale gas
Nw
2
Reference parameter for water
Ng
2
Reference parameter for shale gas
Rsw
1.16
Shale gas solubilityFigure 5
Model verifications of long-term production of a horizontal Barnett shale reservoir.
(a)
Gas production data
(b)
Normalized mean square error between the analytical solution and field data


## 4.1. Convergence of Analytical Solution

The convergence of iterative analytical solutions is checked here. The following mean relative tolerance is defined for this convergence study:(54)δ=dGpt/dtn+1−dGpt/dtndGpt/dtn+1×100%.The mean relative tolerance with the number of iterations is presented in Figure2. The mean relative tolerance is 0.241% for the first iteration, 0.139% for the second iteration, and 0.087% for the 3rd iteration. This indicates that the analytical solutions of the shale gas production rate are convergent with the number of iterations.Figure 2
Convergence check for the iterative analytical solution.


## 4.2. Model Verifications of Water and Gas Production Rates at the Flowback Stage

The early time-pressure-rate data immediately after hydraulic fracturing were collected as flowback data [36]. In the flowback stage, a set of field data of water and gas production rates was obtained from a China shale gas reservoir [20]. Table 1 lists the computational parameters used in analytical solutions. These parameters are taken from Wang and Peng [7] and Yang et al. [20]. Figure 3(a) compares the shale gas production rates among the field production data, the analytical solution of our model, and the semianalytical solution of Yang et al.’s model. Figure 3(b) presents the normalized mean square errors (NMSE) between our analytical solution of shale gas production rates and the field data. According to the fitting results, the NMSE is between 0.1% and 1%. The analytical solution of our model is in good agreement with the field data. Further, Figure 4(a) compares the water production rates of the field production data (black circles), the analytical solution of our model (black line), and the semianalytical solution of Yang et al.’s model (red double dash). Figure 4(b) presents the NMSE of the water production rate between our analytical solution and the field data. According to the fitting results, the NMSE is between 0.1% and 2%. The analytical solution of our model can well predict the field data. These figures indicate the following understandings.Table 1
Model parameters in the computation for the China shale gas reservoir.
Parameter
Unit
Value
Physical meanings
Srg
0.15
Residual saturation of gas
Srw
0.2
Residual saturation of water
μw
Pa∗s
3.6×10−4
Water viscosity
μg
Pa∗s
2.0×10−5
Shale gas viscosity
c
m/s
5.3×10−6
Travel wave viscosity
λ
2
Pore size distribution index
p0
MPa
27.4
Initial reservoir pressure
p1
MPa
19.67
Well pressure
k0
mD
0.1
Initial permeability in the fractured zone
ϕ0
0.18
Initial porosity
pe
MPa
2
Entry capillary pressure
krwmax
0.004
End-point relative permeability for water
krnwmax
1
End-point relative permeability for shale gas
Nw
2
Reference parameter for water
Ng
2
Reference parameter for shale gas
Rsw
5.42
Shale gas solubilityFigure 3
Model verifications of gas production in the flowback period of a China shale reservoir.
(a)
Gas production data
(b)
Normalized mean square error between analytical solution and field dataFigure 4
Model verifications of water production in the flowback period of a China shale reservoir.
(a)
Water production data
(b)
Normalized mean square error between the analytical solution and field dataFirstly, the production rate of shale gas increases to the maximum and then decreases. The gas production rate of our model goes up fast till the peak rate of1.69×105 m3/d at the 24th day and then goes down gently. On the other hand, the water production rate of our model goes down fast from 22 m3/d to 2.8 m3/d during the first 50 days and then keeps flat from 50 to 200 days. These results show that the analytical solutions are in good agreement with field production rates of both shale gas and water.Secondly, for the shale gas production rate, the peak rate of our model is1.69×105 m3/d at the 24th day, while the semianalytical solution of Yang et al.’s model is 1.91×105 m3/d at the 55th day as shown in Figure 3(a). The shale gas production rate of our model goes up faster than that of Yang et al.’s model in the rising period. At the 100th day, the shale gas production rates are 1.86×105 m3/d and 1.65×105 m3/d for the two cases, respectively. On the contrary, the shale gas production rate of our model goes down slower than that of Yang et al.’s model in the decline period. For the water production rate, the production rate of our model decreases from 21.86 m3/d to 7.35 m3/d but the production rate by Yang et al.’s model decreases from 29.14 m3/d to 6.14 m3/d after the first 24 days (see Figure 4(a)). The decrease percentages are 66.4% and 78.9%, respectively. The water production rate predicted by our model goes down slower than that by Yang et al.’s model in the whole 200 days. These differences between our model and Yang et al.’s model may be induced by the nonlinear term of the two-phase flow and the gas solubility in water. These two factors may significantly impact the production rates of shale gas and water. The nonlinear term stands for the complex interaction between shale gas and water. Meanwhile, the dissolved gas may change into free gas during depressurization, leading to the going up of the shale gas production rate by our model faster than that of Yang et al.’s model in the rising period.


## 4.3. Model Verifications for Long-Term Production

The analytical solution of the shale gas production rate is applied to the field production data forecast for the horizontal Barnett shale well 314 [37]. As the field water production rate is not available in the literature, thus, the model verification of the water production rate is not conducted here. Table 2 lists the model parameters for the Barnett shale. The reservoir parameters used in this analysis are determined after referencing the publications by Wang and Peng [7] and Yu and Sepehrnoori [25]. The travelling wave velocity is also decided by the reservoir length and the production time from theirs. The actual field data of shale gas production are presented as the black circle in Figure 5(a). The black line is obtained by the analytical solution of the shale gas production rate. The line is in reasonable agreement with the field production data for the whole production time. The gas production rate declines about 70% in the first 200 days. At the 200th day, the gas production rate reaches 6.3×104 m3/d and then slowly descends during the later period of production. As the bottom hole pressure is set as 3.45 MPa, the gas production rate also keeps about 1.6×104 m3/d in the production tail. Figure 5(b) shows the normalized mean square errors between our analytical solution of shale production rates and the field data. This NMSE is between 0.6% and 1%. Therefore, our analytical solution is feasible to predict the long-term gas production rate.Table 2
Model parameters in computation for the Barnett shale gas well.
Parameter
Unit
Value
Physical meanings
Srg
0.15
Residual saturation of gas
Srw
0.6
Residual saturation of water
μw
Pa∗s
3.6×10−4
Water viscosity
μg
Pa∗s
2.0×10−5
Shale gas viscosity
c
m/s
1.2×10−6
Travel wave viscosity
λ
2
Pore size distribution index
p0
MPa
20.34
Initial reservoir pressure
p1
MPa
3.45
Well pressure
k0
mD
5.0
Initial permeability in the fractured zone
ϕ0
0.2
Initial porosity
pe
MPa
2
Entry capillary pressure
krwmax
0.00035
End-point relative permeability for water
krnwmax
1
End-point relative permeability for shale gas
Nw
2
Reference parameter for water
Ng
2
Reference parameter for shale gas
Rsw
1.16
Shale gas solubilityFigure 5
Model verifications of long-term production of a horizontal Barnett shale reservoir.
(a)
Gas production data
(b)
Normalized mean square error between the analytical solution and field data


## 5. Discussions

In this study, both shale gas solubility and capillary pressure are introduced into the two-phase flow model. The nonlinear term of governing equations is included for the iterative analytical solutions. How these two factors impact the gas production rate of shale gas is analyzed here. In these analyses, the field production data from the Barnett shale well [37] are used.


### 5.1. Impact of the Nonlinear Term on the Shale Gas Production Rate

The gas production rates with and without the nonlinear term are compared in Figure6. At the 200th day, the shale gas production rates are 5.7×104 m3/d with the nonlinear term and 4.4×104 m3/d without the nonlinear term. The gas production rate decreases from its initial production rate by 71.8% and 78.3% with and without the effects of the nonlinear term. Further, the shale gas production rates at the 1200th day are 1.6×104 m3/d and 1.1×103 m3/d, respectively. Without the impact of the nonlinear term, the gas production rate goes down faster and reaches the lower gas production rate. Combined with Figure 3, the gas production rate with the nonlinear term is in better agreement with field production data. Thus, the nonlinear term should not be ignored in the solving process of governing equations for the two-phase flow.Figure 6
Impact of the nonlinear phase of the governing equations for the two-phase flow on the gas production rate.


### 5.2. Impact of Shale Gas Solubility in Water on the Shale Gas Production Rate

The solubility of methane in water is very low under standard temperature and pressure. It may be big at the evolution stage of geological history due to in situ geological condition [9]. This section will study the impact of shale gas solubility on the shale gas production rate.The shale gas solubilities are taken as 1.16, 5.42, and 105. The values of 1.16 and 5.42 for gas solubility are taken from the experimental data obtained by Fu et al. [9]. We also take an extreme value of 105. Figure 7 presents the declines of the gas production rate with time. They all decline rapidly in the first 200 days. Little difference is observed for the history curves when the shale gas solubilities are 1.16 and 5.42. Difference is observed for the shale gas solubility value of 105. The shale gas production rates at the 200th day are 5.6×104 m3/d, 5.6×104 m3/d, and 5.2×104 m3/d, respectively. Moreover, the shale gas production rates at the 1200th day are 1.6×104 m3/d, 1.6×104 m3/d, and 1.1×104 m3/d, respectively. For the extreme case, the gas production rate goes down faster in the gas production stage and reaches the lower gas production rate in the production tail. This implies that the gas production rate is affected by the shale gas solubility when shale gas dissolves heavily in water. In that case, the large amount of dissolved gas will be released into free gas for production.Figure 7
Impact of gas solubility on the gas production rate.


### 5.3. Impact of Entry Capillary Pressure on the Shale Gas Production Rate

Three entry capillary pressures of 1.5 MPa, 2.0 MPa, and 2.5 MPa are assumed to investigate the impact of entry capillary pressure on the gas production rate. The initial reservoir pressure is 20.34 MPa. Figure8 presents the gas production rate with time. All three cases have the same initial shale gas production rate of 2.02×105 m3/d. At the 200th day, the shale gas production rates become 4.4×104 m3/d, 5.6×104 m3/d, and 9.5×104 m3/d, respectively. It decreases from its initial production rates by 77.9%, 71.8%, and 53.1% for entry capillary pressures of 1.5 MPa, 2.0 MPa, and 2.5 MPa, respectively. Further, at the 1200th day, the shale gas production rates are 3.7×103 m3/d, 1.8×104 m3/d, and 6.7×104 m3/d, respectively. These results indicate that the entry capillary pressure has a significant impact on the production of shale gas. In the same production period, the shale gas production rate with higher entry capillary pressure declines faster and reaches the lower production rate. This is because higher entry capillary pressure causes more gas to flow to the wellbore faster. In fact, entry capillary pressure affects the nonlinear term. The capillary entry pressure in relation to the applied pressure difference determines whether the gas will move or get trapped in water [38]. When the interface pressure achieves the entry capillary pressure, the gas begins a relative motion with water, thus starting the interaction. In addition, the water saturation can be expressed by the relationship between entry capillary pressure and capillary pressure through a normalized saturation of the water phase. The nonlinear term is produced after the coupling of water or gas saturation and pressure of water or gas in the governing equations.Figure 8
Impact of entry capillary pressure on the gas production rate.


## 5.1. Impact of the Nonlinear Term on the Shale Gas Production Rate

The gas production rates with and without the nonlinear term are compared in Figure6. At the 200th day, the shale gas production rates are 5.7×104 m3/d with the nonlinear term and 4.4×104 m3/d without the nonlinear term. The gas production rate decreases from its initial production rate by 71.8% and 78.3% with and without the effects of the nonlinear term. Further, the shale gas production rates at the 1200th day are 1.6×104 m3/d and 1.1×103 m3/d, respectively. Without the impact of the nonlinear term, the gas production rate goes down faster and reaches the lower gas production rate. Combined with Figure 3, the gas production rate with the nonlinear term is in better agreement with field production data. Thus, the nonlinear term should not be ignored in the solving process of governing equations for the two-phase flow.Figure 6
Impact of the nonlinear phase of the governing equations for the two-phase flow on the gas production rate.


## 5.2. Impact of Shale Gas Solubility in Water on the Shale Gas Production Rate

The solubility of methane in water is very low under standard temperature and pressure. It may be big at the evolution stage of geological history due to in situ geological condition [9]. This section will study the impact of shale gas solubility on the shale gas production rate.The shale gas solubilities are taken as 1.16, 5.42, and 105. The values of 1.16 and 5.42 for gas solubility are taken from the experimental data obtained by Fu et al. [9]. We also take an extreme value of 105. Figure 7 presents the declines of the gas production rate with time. They all decline rapidly in the first 200 days. Little difference is observed for the history curves when the shale gas solubilities are 1.16 and 5.42. Difference is observed for the shale gas solubility value of 105. The shale gas production rates at the 200th day are 5.6×104 m3/d, 5.6×104 m3/d, and 5.2×104 m3/d, respectively. Moreover, the shale gas production rates at the 1200th day are 1.6×104 m3/d, 1.6×104 m3/d, and 1.1×104 m3/d, respectively. For the extreme case, the gas production rate goes down faster in the gas production stage and reaches the lower gas production rate in the production tail. This implies that the gas production rate is affected by the shale gas solubility when shale gas dissolves heavily in water. In that case, the large amount of dissolved gas will be released into free gas for production.Figure 7
Impact of gas solubility on the gas production rate.


## 5.3. Impact of Entry Capillary Pressure on the Shale Gas Production Rate

Three entry capillary pressures of 1.5 MPa, 2.0 MPa, and 2.5 MPa are assumed to investigate the impact of entry capillary pressure on the gas production rate. The initial reservoir pressure is 20.34 MPa. Figure8 presents the gas production rate with time. All three cases have the same initial shale gas production rate of 2.02×105 m3/d. At the 200th day, the shale gas production rates become 4.4×104 m3/d, 5.6×104 m3/d, and 9.5×104 m3/d, respectively. It decreases from its initial production rates by 77.9%, 71.8%, and 53.1% for entry capillary pressures of 1.5 MPa, 2.0 MPa, and 2.5 MPa, respectively. Further, at the 1200th day, the shale gas production rates are 3.7×103 m3/d, 1.8×104 m3/d, and 6.7×104 m3/d, respectively. These results indicate that the entry capillary pressure has a significant impact on the production of shale gas. In the same production period, the shale gas production rate with higher entry capillary pressure declines faster and reaches the lower production rate. This is because higher entry capillary pressure causes more gas to flow to the wellbore faster. In fact, entry capillary pressure affects the nonlinear term. The capillary entry pressure in relation to the applied pressure difference determines whether the gas will move or get trapped in water [38]. When the interface pressure achieves the entry capillary pressure, the gas begins a relative motion with water, thus starting the interaction. In addition, the water saturation can be expressed by the relationship between entry capillary pressure and capillary pressure through a normalized saturation of the water phase. The nonlinear term is produced after the coupling of water or gas saturation and pressure of water or gas in the governing equations.Figure 8
Impact of entry capillary pressure on the gas production rate.


## 6. Conclusions

In this study, a mathematical model of the two-phase flow in shale gas production was developed with extra consideration of gas solubility. This mathematical model was then analytically solved with travelling wave and variational iteration methods. Being different with previous analytical solutions, the nonlinear term in the governing equations of the two-phase flow was not ignored in the solutions of our model. The analytical solution of the shale gas production rate was compared with the field production data at both the flowback stage and long-term production from a horizontal Barnett shale gas reservoir and a China shale gas reservoir, respectively. The analytical solution of the water production rate was also compared with the field production data from the China shale gas reservoir. The effects of shale gas solubility and entry capillary pressure on the shale gas production rate were finally analyzed. Based on these studies, the following conclusions can be made.Firstly, the nonlinear term of governing equations is an important factor in the analytical solutions for the two-phase flow. With the impact of the nonlinear term, the predicted shale gas production rate falls slower and the total predicted shale gas production is higher at the same period. The analytical solutions with the nonlinear term can better predict the production rates of gas and water in the flowback stage.Second, shale gas solubility has little influence on the Barnett shale gas production rate when its value is not so large. When the shale gas solubility value is high enough, such as a gas solubility value of105, its impact on the shale gas production rate becomes bigger in the later stage of shale gas production. The shale gas production rate would decline faster and reach the lower gas production rate in the production tail.Finally, entry capillary pressure has a significant impact on the production rate of shale gas. In our computation cases, higher entry capillary pressure leads to a faster decline of the shale gas production rate and a lower production rate at the same time period.



---

*Source: 4943582-2019-05-08.xml*