The flowback and production analysis in sub-saturated fractured shale reservoirs

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

Contents lists available at ScienceDirect

Journal of Petroleum Science and Engineering journal homepage: http://www.elsevier.com/locate/petrol

The flowback and production analysis in sub-saturated fractured shale reservoirs Shiming Wei, Yan Jin *, Yang Xia, Botao Lin State Key Laboratory of Petroleum Resource and Prospecting, China University of Petroleum, Beijing, China

A R T I C L E I N F O

A B S T R A C T

Keywords: Sub-saturated shale reservoirs Flowback Production Gas-water flow Discrete fractures

The initial water saturation of shale reservoirs is in general lower than the irreducible water saturation. However, the hydraulic fracture network region is highly saturated by fracturing fluid, which results in a two-phase flow inside fractured region. A long-time two-phase flow period is always observed in the fields even after a high flowback rate, which cannot be explained by existing models. In this paper, we develop a 2D gas-water flow model coupling multi-scale discrete fractures and continuum framework to describe the flowback process of a multi-stage fractured horizontal well in shale gas reservoirs. The model considers the non-linear flow mechanism of free gas, desorption of adsorbed gas and the co-existence of single gas flow and gas-water flow. The capillary pressure model and relative permeability model are both employed with different parameter values for matrix and fractures, respectively. The permeability of naturally fractured shale is described using a fractal model. To overcome non-convergence problem of finite element method for solving two-phase flow in cross fractures, the model is solved by finite element method and finite volume method. The model is verified through a simple twophase flow simulation and a single gas phase simulation with some simplifications. Then the model is used for history matching as well as production prediction of a well in Changning block, Sichuan, China. The simulation results show that the high water saturation induced “skin factor” in the vicinity of hydraulic fractures signifi­ cantly delays the gas production for a long time because some of the fracturing fluids become irreducible water and lower the gas phase permeability. Due to the great influence of fracture interference on water flow, the fracturing fluids concentrate on regions near hydraulic fractures and their junctions. The long-term gas pro­ duction of the complex fracture network is slightly higher than that of the simple fracture network, which is different with the simulation results of single gas flow model.

1. Introduction Shale gas reservoirs are characterized by ultra-low permeability, ultra-low porosity and well-developed natural fractures, thus the effi­ cient exploitation of shale reservoirs needs multi-stage hydraulic frac­ turing technology to extract the industrial gas flow (Carpenter. 2014; Warpinski et al., 2012; Denney. 2009). Due to the tight relation between the geological features of hydraulic fracture networks and shale gas production, a lot of works has been done to estimate the gas production of shale reservoirs with the consideration of construction and optimi­ zation of discrete fracture networks (Wei et al., 2019a,b; Wang et al., 2017; Rammay and Awotunde. 2016; Wang and Shahvali. 2016; Wilson and Durlofsky. 2013; Basquet et al., 2004). However, the early pro­ duction stage, i.e. flowback stage, exhibits a gas-water flow period. Even if the fracture description is precise, it is impossible to accurately

simulate the shale gas flow without considering the flow of fracturing fluid. On the other hand, flowback of fracturing fluids after the stimu­ lation of a well is a necessary step to condition the well for long-term production performance (Crafton and Gunderson. 2007). Robinson et al. (1988), Ely et al. (1990) and Barree and Mukherjee (1995) have proposed three different flowback modes, which aimed at keeping the flow conductivity of hydraulic fractures as high as possible (Barree and Mukherjee. 1995; Ely et al., 1990; Robinson et al., 1988). Flowback data after hydraulic fracturing gives the earliest opportu­ nity to characterize the properties of hydraulic fracture networks (Clarkson and Williams-Kovas, 2013a; 2013b). And researchers have realized the potential of the flowback data to extract key parameters such as the volume of hydraulic fracture networks, fracture length and fracture permeability, etc. (Cugnart et al., 2017; Williams-Kovacs and Clarkson, 2015; Ezulike et al., 2014; IIk et al., 2010), which are essential

* Corresponding author. E-mail address: [email protected] (Y. Jin). https://doi.org/10.1016/j.petrol.2019.106694 Received 20 May 2019; Received in revised form 11 November 2019; Accepted 13 November 2019 Available online 16 November 2019 0920-4105/© 2019 Elsevier B.V. All rights reserved.

Please cite this article as: Shiming Wei, Journal of Petroleum Science and Engineering, https://doi.org/10.1016/j.petrol.2019.106694

S. Wei et al.

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

to estimate the long-term gas production. There is no apparent flowback stage during shale gas exploitation in Sichuan, China, i.e. the early production stage is the flowback stage. Therefore, we need a method to simulate the flowback process and the production process continuously. Field data indicates that gas-water flow presents a long time in the production process, because a great amount of the fracturing fluid flows into the shale matrix through the hydraulic fracture networks that is controlled by various mechanisms (Liu et al., 2018), and most of the fracturing fluid becomes the irre­ ducible water. Then the shale matrix in the stimulated area, i.e. hy­ draulic fracture network region, is highly saturated with fracturing fluids. However, the area away from the hydraulic fractures is not saturated by the fracturing fluid, and its water saturation is still the initial water saturation, which is lower than the irreducible water saturation (Birdsell et al., 2015). The initial water saturation is lower than the irreducible water saturation, and it is neglected by previous literatures. Taken the commonly used exponential model of relative permeability as example (Yang et al., 2017; Brooks and Corey, 1964), when the initial water saturation is lower than its irreducible water saturation, the effective water saturation is negative, which results in the negative relative permeability for both the water phase and the gas phase. That is the reason why previous researchers always study the flowback process and the production process separately, and the initial water saturation is always set higher than its irreducible water satura­ tion in the flowback simulation (Yang et al., 2017). In order to simulate the flowback process and the production process continuously, we have to set a single gas flow away from hydraulic fractures and a gas-water flow near hydraulic fractures. In this paper, the relative permeability function is fitted with experimental data, using the water saturation as the variable, not the effective water saturation. The prerequisite for an accurate simulation on flowback is to reasonably consider some of the fracturing fluid as irreducible water after hydraulic fracturing treatment. The flow direction of water must be selected according to both the water pressure and water saturation, i.e. the transmissibility of water is decided by both the water pressure and water saturation. The finite volume method (FVM), employing the “upwind” scheme, has been widely used to solve the multi-phase flow in discrete fracture networks and continuous matrix pores. Reichenberger et al. (2006) presented a vertex-centered finite volume method for the gas-water flow in fractured porous media, which modeled fractures as lower dimensional elements (Reichenberger et al., 2006). Using the finite volume method, a three-dimensional model was presented by Alkhlaifat and Arastoopour, (2003). Li and Lee (2008) proposed a hybrid finite volume method to solve the gas-oil-water flow in discrete fracture networks and homoge­ nized media (Li and Lee. 2008). However, these works cannot be used to simulate gas-water flow in shale gas reservoirs due to the complex flow mechanisms of shale gas. Besides, the relative permeability curves and capillary curves in discrete fractures should be different with that in continuum media. The flow mechanisms of shale gas are complicated, including viscous flow, slip flow, Knudsen diffusion, surface diffusion, etc. (Wu et al., 2016; Sun et al., 2015; Sheng et al., 2015). And the shale gas can exist in two forms: free gas and adsorbed gas. The above phenomena enhance the nonlinearity of the gas flow equation. With the evolution of nu­ merical methods and computational speed of computers, different hi­ erarchical models coupling discrete fractures and multi-continuum (i.e. multi-porosity) were proposed, and the nonlinear flow mechanisms of shale gas are considered (Xia et al., 2019; Wei et al., 2019a,b; Xia et al., 2017; Akkutlu et al., 2017; Zhang et al., 2015). As the most important flow channels in shale reservoirs, hydraulic fractures (macro-fractures) are built as discrete fractures or embedded discrete fractures (Ding et al., 2017; Wang et al., 2017; Jiang and Younis, 2016; Fan et al., 2015; Moinfar et al., 2013). In order to simplify the calculation, the concept of SRV is introduced, and the effect of secondary fractures (intersecting with main hydraulic fractures) is replaced by an estimation of the

effective permeability in different flow regions (Wu et al., 2014; Fan et al., 2015). The naturally fractured shale matrix consists of matrix pores and natural fractures (micro-fractures). Due to the random dis­ tribution and large number of natural fractures, it is unwise to describe natural fractures with discrete fracture model. Warren and Root firstly introduced the dual porosity model to describe the oil and gas transport in the naturally fractured rocks (Warren and Root. 1963). Each contin­ uum has their porosity and permeability, and the flow in each contin­ uum is connected through a mass source term. The naturally fractured shale matrix can be subdivided into organic matters, organic pores, inorganic pores and natural fractures, the multi-porosity model can be deduced by referring to the dual porosity model (Wei et al., 2019a,b; Fung and Du. 2016; Zhang et al., 2015; Jiang and Younis, 2015). However, shape factors are needed to calculate the mass transfer among different continua for the multi-porosity model, and the shape factors are related to the matrix block size (Warren and Root. 1963), which are empirical parameters. For the dual/multi-porosity model, the shale matrix permeability is so small, thus it is hard for the mass transfer between matrix pores and natural fractures to achieve the pseudo-steady-state. Therefore, the multiple interacting continua (MINC) model is proposed to overcome this deficiency of dual/multi-porosity model, and the hybrid model coupling discrete fracture with MINC model in shale gas simulation gets better results than that with dual/multi-porosity model. The self-similar fractal structures of natural fracture networks were extensively studied (Sahimi. 2012; Bonnet et al., 2001). Numbers of researchers reported that the re­ lationships between the length and the number of natural fractures exhibit the power-law, exponential and log-normal types (Kolyukhin and Torabi. 2012; Andrade et al., 2009; Watanabe and Takahashi. 1995; Chelidze and Guguen. 1990; Chang and Yortsos. 1990). Also, the rela­ tionship between aperture and length of natural fractures can be described with an exponential law (Schultz et al., 2008; Hatton et al., 1994). The fractal nature of natural fractures allows the derivation of the permeability of naturally fractured rock, therefore, the overall perme­ ability can be expressed as the function of several statistical fracture parameters, like fractal dimension, fracture density, fracture azimuth and facture dip angle. (Miao et al., 2015). Proposed hybrid models are focused on the division of the pore structure, which are single gas flow model. Even though the hybrid models captured the multi-continuum nature of the shale and the interaction between shale matrix and hy­ draulic fractures, they failed to describe the gas and water production behaviors, especially at the early production stage. In our experiments, testing parameters (like porosity, relative permeability and capillary pressure curve) have already included the influence of natural fractures and matrix pores, hence the naturally fractured shale matrix is described using the fractal model (Miao et al., 2015). Summarizing previous literatures, a gas-water flow model consid­ ering complex gas transport mechanisms is needed to simulate the flowback process and long production process continuously. In this paper, the naturally fractured shale matrix is taken as a single contin­ uum, and its permeability is described by a fractal model. Both the viscous flow and Knudsen diffusion is considered in the naturally frac­ tured shale matrix. The hydraulic fractures are described with the discrete fracture model. The capillary curves and relative permeability curves are different in the naturally fractured shale matrix and hydraulic fractures. We combine the finite element method (FEM) and the finite volume method (FVM) to write the discrete equations of the mathe­ matical model for the purpose of overcoming the selection of water flow direction. The model and the program are verified through two special cases by simplifying our model. Then our model is applied to an actual well in Changning block, Sichuan, China. Finally, the influences of water saturation near hydraulic fractures, parameters of natural fractures, gas desorption and different hydraulic fracture networks are analyzed.

2

S. Wei et al.

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

2. Physical equations 2.1. Permeability of naturally fractured shale Shale is a kind of naturally fractured porous rock, which consists of multi-scale matrix pores and natural fractures (Shahid et al., 2015; Clarkson et al., 2013; Cho et al., 2013). It is neither efficient nor precise to model all the natural fractures as discrete fractures, especially when they co-exist with large scale hydraulic fractures (Wei et al., 2019a,b). However, the length and the number of natural fractures exhibit certain statistical laws, such as power-law, exponential and log-normal types (Miao et al., 2015). The relationship between the aperture of fracture and the fracture length is formulated as (Hatton et al., 1994; Schultz et al., 2008): Fig. 1. Water saturation in the shale reservoir.

(1)

a ¼ βln

where a and l are respectively the aperture and length of natural frac­ tures, m. β is the proportionality coefficient. n is a constant. Based on Eq. (1), with n ¼ 1, Miao et al. (2015) has proposed a fractal model of permeability for fractured rocks (Miao et al., 2015): � β2 D 1 Df l3max 1 cos2 α sin2 θ k¼ (2) 1 Df 128 4 Df 2 D 1 φf f Df ¼ d E þ

ln φ lnðlmax =lmin Þ

(3)

0 2 D¼

1

� B Df φf B @1 � Df βlmax 1

The relative permeability curve and capillary curve of hydraulic frac­ tures are both described as linear function of saturation. Fig. 2 and Fig. 3 are respectively the relative permeability curves and the capillary pressure curves of shale matrix and hydraulic fractures. When fitting the relative permeability curve with the core test of shale matrix, the gas saturation (i.e. 1 sw ) is chose as the horizontal axis in the relative permeability curve. It is can be seen from Fig. 2 that the irreducible water saturation of shale matrix is 0.6452, but the initial water saturation in shale matrix is only 0.4. Eq. (5) is the most used Brooks-Corey relation between water saturation and relative permeability (Brooks and Corey, 1964). 0 1

1 1 Df 2 Df

φf

φf

C C A �

2þ3λ

krw ¼ se λ ;

krg ¼ ð1

2B se Þ @1

2þλ C se λ A

(5)

(4) pc ¼ pd se 1=λ

where k is the permeability of fractured porous rocks, m2. Df is the fractal dimension of fracture length. D is the fracture density, 1/m. D is the liner density, which is the ratio of total fracture length to the area of reservoir. α and θ are respectively the mean facture azimuth and facture dip angle. ∅f is the porosity of natural fractures. dE is the dimension of simulated area. ∅ is the total porosity of fractured porous medium. lmax and lmin are respectively the maximum and minimum length of natural fractures, m. In this paper, Eq. (2) is used to give the overall permeability of shale matrix containing natural fractures.

(6)

where se is the effective water saturation, and se ¼ ðsw swr Þ=ð1 swr sor Þ. λ is related to the pore size distribution: materials with smaller variations in pore size have a larger λ value. Usually λ is in the range 0:2 � λ � 3 (Reichenberger et al., 2006). krw and krg are respectively the relative permeabilities of water and gas. pc is the capillary pressure, MPa. pd is the minimum capillary pressure, MPa. Even though the parameters in Eqs. (5) and (6) have very clear physical meanings, the effective water saturation is negative when the initial water saturation is lower than the irreducible water saturation. The relative permeability and capillary pressure will both be negative, which is against physical laws. Therefore, we directly use the water saturation rather than the effective water saturation to fit the experi­ mental data. The relative permeability function and capillary pressure function of hydraulic fractures are described by linear functions of water saturation (Yang et al., 2017). As we known that the hydraulic fracture aperture is much larger than the matrix pore diameter. Therefore, the capillary pressure in hydraulic fractures is much smaller than that in the shale matrix. Fig. 4 shows the capillary pressure curves of shale matrix and hydraulic fractures. The phase pressure in hydraulic fractures is continuous, thus the capillary pressure is continuous. However, two different capillary pressure func­ tions are selected to describe the relationship between capillary pressure and water saturation, which results in discontinuous water saturation in hydraulic fractures. In order to make sure the water saturation contin­ uous in hydraulic fractures, the water saturation in hydraulic fractures is transformed into the water saturation in shale matrix according to the continuous capillary pressure, i.e. Eq. (7). There is a critical water saturation s* that once the water saturation in shale matrix is lower than s*, the capillary pressure in shale matrix cannot be matched in hydraulic fractures. Therefore, the water saturation in hydraulic fractures is zero once the water saturation in shale matrix is lower than s*.

2.2. Two-phase flow area The initial water saturation of shale matrix is lower than conven­ tional reservoirs. Field observations in Sichuan province of China indi­ cate that the initial water saturation of shale is smaller than the irreducible water saturation. During the hydraulic fracturing process, water enters the shale matrix through hydraulic fractures. Therefore, the water saturation near the hydraulic fracture networks is higher than the irreducible water saturation. In the flowback process, the two-phase flow is only confined in the fracture network region. Fig. 1 shows the water distribution in the shale reservoir. The green area is the highly saturated area, which happens two-phase flow. Outside of the green area is single phase flow area. 2.3. Relative permeability and capillary pressure The relative permeability curves and capillary curves are different in shale matrix and hydraulic fractures. A linear function was selected to describe the relative permeability curve of hydraulic fractures (Yang et al., 2017). In this paper, the relative permeability curve and capillary pressure curve of shale matrix are both obtained from the experimental data. And shale cores used in experiments are from Changning oil field. 3

S. Wei et al.

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

Fig. 2. The relative permeability of the shale matrix (left) and hydraulic fractures (right).

Fig. 3. The capillary pressure of shale matrix (left) and Hydraulic fractures (right).

fractured shale matrix are respectively:

∂ðφρw sw Þ ∂t ∂ φρg sg ∂t

�

r⋅ðKρw λrw rpw Þ ¼ mw

(8)

� r⋅ Kρg λrg rpg ¼ mg

(9) (10)

sw þ sg ¼ 1

where φ is the total porosity of the naturally fractured shale matrix. λrw is the relative mobility of water, and λrw ¼ krw =μw . K is the permeability matrix of the naturally fractured shale matrix. pw and sw are respectively the pressure and saturation of water phase. pg and sg are respectively the pressure and saturation of gas phase. λrg is the relative mobility of gas, ~rg is the apparent permeability of gas phase when we and λrg ¼ ~ krg =μ . k g

consider different flow mechanisms. ρw and ρg are respectively the water density and gas density, kg/m3. mw and mg are respectively mass sources of water and gas, kg/(m3⋅s). The experimental results of the gas transport speed in a 10 nm capillary is comparable with that in a 100 nm capillary, which is a solid proof that there must be other flow mechanisms enhanced gas transport, and it cannot be defined by the conventional convection model (Roy et al., 2003). The gas transport enhancement is named as Knudsen diffusion, which plays an important role in nanopores (Wang and Li, 2003). Previous literatures indicated that Knudsen diffusion can contribute to the accumulated production up to 20% (Darabi et al., 2012). In this paper, we consider methane as an ideal gas, thus we have the following relations among gas concentration, gas density and gas pressure.

Fig. 4. Capillary pressure curves in shale matrix and a fracture.

� sw

F

¼

0; ðpc F Þ 1 pc m ðsc

where sw

F

m Þ;

sw sw

m m

� s* > s*

(7)

is the water saturation in hydraulic fractures.

2.4. Non-linear flow mechanism of shale gas

pg ¼ Cg RT;

The mass conservation equation of gas and water in the naturally 4

ρg ¼ MCg

(11)

S. Wei et al.

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

where pg , Cg and ρg are respectively the gas pressure, gas concentration and gas density. M is the molecular mass of methane, kg/mol. R is the gas constant, J/(mol⋅K). T is the formation temperature, K, which is Kelvin temperature and is a constant in this paper. Considering the viscous flow and Knudsen diffusion mechanisms, the gas flux through unit cross-sectional area of the shale matrix can be expressed as (Wei et al., 2019a,b): Jm ¼

~krg

μg

Fig. 5. Two nodes with different pressure and water saturation.

(12)

rpg

~krg ¼ krg þ φDk μ pg K

near the wellbore. So hydraulic fractures are built as discrete fractures in this paper. The aperture of hydraulic fractures is usually in the scale of mm ~ cm, however, the length of hydraulic fractures is over 100 m. Therefore, the fluid velocity profile across the hydraulic fracture section is seen as a straight line, and it is perpendicular to the hydraulic fracture. Hydraulic fractures are propped with proppants to maintain their con­ ductivity in flowback and production process. Therefore, hydraulic fractures can be seen as the porous medium. And Darcy’s law is employed to describe the gas and water flow in hydraulic fractures. In hydraulic fractures, the contribution of adsorbed gas and Knudsen diffusion to the free gas flow is neglected, thus the control equations in hydraulic fractures is: 8 > < � c s þ c s � 0 �� p_ � tg g tw w w wF � φF > s_w F ctw sw 1 :

(13)

~rg is the apparent permeability of shale gas considering viscous where k flow and Knudsen diffusion. Dk is the Knudsen diffusion coefficient, m2/ s. 2.5. Desorption/adsorption phenomenon As a kind of source rock, shale contains both free gas and adsorbed gas, which makes the gas desorption/adsorption process happen during production. Langmuir (1918) first proposed the isothermal Langmuir equation to describe the gas adsorption, and his remarkable research is republished in 2015 (Langmuir, 1918). The isothermal Langmuir equation is widely used to describe the state of adsorbed gas (Brohi et al., 2011): Cμ ¼

Ke Cμs C Ke Cμs pg ¼ Cμs þ Ke C Cμs RT þ Ke pg

2

(14)

Considering the desorption of adsorbed gas and Knudsen diffusion of the free gas, the control equations in the naturally fractured shale matrix is given below, and the detailed derivation of the mathematic model is shown in Appendix A.

φctw sw

Ke C2μs RT Cμs RT þ Ke pg

�2

2 3 0 7� p_ � 6 Kλt 5 w þ4 s_w Kλrw φ

Kλrg

0

9 3 �> � � 2 qg F þ qw r pw = 7 F5 ¼ qw F r2 sw F > ;

�

F

(17)

F

qw

Γ

¼

2K

qg

Γ

¼

2K

krw ðsw

mÞ

μ krg ðsw

μ

rpw

mÞ

rðpw

(18)

m ⋅n

m

þ pc

m Þ⋅n

(19)

where qw Γ and qg Γ are respectively the water and gas transfer between shale matrix and hydraulic fractures, m/s, and they exist only in the surfaces of hydraulic fractures. Therefore, Eq. (18) and Eq. (19) don’t appear in Eq. (16). n is the normal vector of hydraulic fractures. The specific mathematic model is shown in Appendix A.

2.6. Control equations in the matrix and hydraulic fractures

φÞM

dpc dsw

where φF is the porosity of the proppant pack inside of hydraulic frac­ tures. wF is the fracture aperture, m. qw F and qg F are respectively the water and gas source terms in hydraulic fractures, m/s. The mass transfer between hydraulic fractures and the shale matrix can be described using the normal pressure gradient of hydraulic frac­ tures, i.e. (Kazemi et al., 1969):

where mg is the change of adsorbed gas when gas pressure changes, kg/ (m3⋅s). And Eq. (15) is the mass source term in Eq. (9).

� 6 φ ctg sg þ ctw sw þ ð1 4

Kλrg

Kλrw

where Cμ and C are respectively the gas concentration of adsorbed gas and free gas in shale matrix, mol/m3. Cμs is the maximum adsorption capacity of shale, and Cμs ¼ pL Ke =R=T (Wei et al., 2019a,b). Ke is the equilibrium constant of adsorption/desorption. The change of adsorbed gas is the supplement of free gas, which can be expressed with the time derivative of gas pressure. # " Ke C2μs RT ∂Cμ ∂pg mg ¼ ð1 φÞM ¼ ð1 φÞM (15) �2 ∂t ∂t Cμs RT þ Ke pg

2

Kλt

6 þ4

3. Numerical treatment

3 dpc m � � � qg m þ qw r2 pw dsw 7 ¼ 5 2 qw m r sw 0

where ctg and ctw are respectively the total compressibility of gas and water phase, MPa 1, and their specific expression can be found in Ap­ pendix A. qg m and qw m are respectively the source terms of gas phase and water phase in the naturally fractured shale matrix, s 1. Compared with natural fractures, hydraulic fractures are larger scale fractures, but smaller in the quantity, and they are distributed in the area

� m

(16)

3.1. The treatment of the control equations with FEM and FVM method Fig. 5 shows two nodes with different pressure and water saturation. As we can see from Fig. 5 that the water is not supposed to flow from node 1 to node 2, i.e. the flow coefficient of water from node 1 to node 2 is zero, because the water saturation in node 1 is lower than the 5

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

S. Wei et al.

irreducible water saturation. If finite element method (FEM) were applied to process the water flow, the flow coefficient from node 1 to node 2 is the interpolation of node 1 and node 2, therefore, the flow coefficient is non-zero. Because the water pressure in node 1 is higher than that in node 2, there will be non-zero water flux from node 1 to node 2. Then the water saturation in node 1 will be lower than the irreducible water saturation, which is unreasonable. Therefore, the flow direction of water must be selected according to both the water pressure and water saturation, i.e. the transmissibility of water is decided by both

2 6 φ ctw sw þ ð1 4

sw Þctg þ ð1

Ke C2μs RT

φÞM

Cμs RT þ Ke pg

ctw sw M

! �2 M

node i to node j and node i to node m. Eq. (22) is the stiffness matrix we get using FVM, and the difference between it with FEM is that the relative mobility of water is perfectly compatible with “upwind” scheme. Thus, in our program, the relative mobility of water in the stiffness matrix is calculated with nodal water saturation that selected by nodal pressure. Other terms of Eq. (20) is discretized using FEM method, thus the discretized form of flow equation in shale matrix is given below, which neglects the effect of porosity change on water saturation.

2 3 � ~krg 0 7� p_ � 6 Tðλw Þ þ T λrg k w þ K4 rg 5 s_w Tðλ Þ rw M

the water pressure and water saturation. There are two reasons why we use the hybrid numerical method coupling the finite element method (FEM) and the finite volume method (FVM). Firstly, in order to better solve the flow direction of water, the “upwind” scheme is employed through the use of FVM, and the detailed derivation is given in Appendix B. For another reason is that it facilitates us to use the FEM gridding method and the already existed FEM program, if we only process the flow item with FVM and the others with FEM. Using the finite volume method (FVM), the mass conservation equation of water in matrix space can be written as: Z Z Z ∂ ðρw φsw ÞdΩ þ rðρw vw ÞdΩ ¼ ρw qw m dΩ (20) ∂ t Ωm Ωm Ωm

T λrg

Tii ¼

Tij þ Tim (

j→i λi→j rw ¼ λrw ¼

∂Ni ∂Nm ∂Ni ∂Nm þ ∂x ∂x ∂y ∂y

�

m

¼ Ωm

Z qg

m¼ Ωm

(27)

∂Ni ∂Nj ∂Ni ∂Nj þ ∂x ∂x ∂y ∂y

� (29) � (30)

ðmw = ρw ÞdΩ

(31)

� � mg ρg dΩ

(32)

The whole space of the hydraulically fractured shale consists of naturally fractured shale matrix and discrete hydraulic fractures, i.e.: Xfracnum Ω ¼ Ωm þ wF;i � ΩF;i (33) i¼1 where fracnum represents the number of hydraulic fractures. Eq. (33) shows the superposition law of shale matrix space and hydraulic fracture space, the treatment of mass conservation equation in discrete fractures is the same of that in matrix. In order to employ the superposition law, the variables of the dis­ cretized control equations must be consistent. Therefore, the discretized control equation of hydraulic fractures is represented with the water phase pressure and water saturation in shale matrix. 8 > �� � < � c s þ ð1 s Þc �M p_w 0 tw w w tg F wF � φF > s_w c s M a M tw w F f F : (34) 9 2 � dpc f 3� �> " # = Tðλ Þ T λ t rg qw f þ qg f pw 6 dsw 7 þ Kf 4 ¼ 5 qw f sw > ; 0 Tðλrw Þ

(23) � (24) (25)

where af ¼ dsw f =dsw m . The relation of water saturation in shale matrix and discrete fractures is derived from the capillary pressure functions, which are fitting functions of core testing data. Z MF ¼ NT NdΩ (35)

�

i λi→j rw sw ; ​ ​ pi > pj � i→j j λrw sw ; ​ ​ pi < pj

∂Ni ∂Nj ∂Ni ∂Nj þ ∂x ∂x ∂y ∂y

Z qw

�

� Tim ¼ Tmi ¼ λi→m rw

� ¼ λi→j rg ij �

where

∂Ni ∂Nj ∂Ni ∂Nj þ ∂x ∂x ∂y ∂y

�

Tðλrw Þij ¼ λi→j rw

where λrw is the relative mobility of water, and λrw ¼ krw = μw . Eq. (21) represents the water mass flux from control volume i to control volume j and m. After FVM discretization of Eq. (21), which is descripted detailedly in Appendix B, the mass flux vector of a triangular grid shown in Fig. A1 is: 0 1 10 1 0 Tii Tij Tim f pi ! @ iA 1 f ¼ fj ¼ ρw K⋅@ Tji Tjj Tjm A@ pj A (22) 2 Tmi Tmj Tmm fm pm

�

� m

Ωm

Γ

Tij ¼ Tji ¼ λi→j rw

3 � ~krg dpc m � � � q þ qg p krg dsw 7 5 w ¼ w m sw qw m 0

where M and T are respectively the mass matrix and stiffness matrix. pw and sw represent the water phase pressure and water saturation in shale matrix. The superscript dot of pw and sw means their time derivative. Z M¼ NT NdΩ (28)

where ρw and sw are respectively the water density and water saturation. vw is the flow velocity of water. Ωm represents the matrix space. The second item in the left side of Eq. (20) is the water flux. In the subdivided grid, the flow term experienced a dimension reduction with Gauss’s law of integration (Spiegel et al., 2009). Here we do not consider the anisotropy of matrix permeability. Z Z fi ¼ r⋅ðρw Kλrw rpw ÞdΩ ¼ ðρw Kλrw rpw Þ⋅ndΓ (21) Ωm

T λrg

(26)

i→m where λi→j rw and λrw are respectively the relative mobility of water from

ΩF

6

S. Wei et al.

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

Fig. 6. Assembly of Eq. (29) and Eq. (36). From Zhang et al., 2013.

Fig. 7. Treatment of discrete fractures. From Wei et al., 2019a,b.

Superpose Eq. (27) and Eq. (34) as Eq. (33) and Fig. 6 shows, we can obtain the final weak form of the model coupling shale matrix and discrete fracture. In this paper, LU decomposition is used to precondition the matrix equation (Bunch and Hopcroft, 1974), and Newton-Simpson iteration method (Kelley. 2003) is used to solve it.

Table 1 Simulation parameters.

3.2. The treatment of the mass transfer between matrix and hydraulic fractures

Parameters

Values

Explanations

km

10 D

Matrix permeability

wF

1 mm

Fracture aperture

kF

w2F =12

cw

Fracture permeability 10

Pa

1

Water compressibility coefficient

μw

1 mPa s

μo

5 mPa s

Oil viscosity

0.4

Matrix porosity

1

Fracture porosity

co

The numerical measurement of coupling the matrix and hydraulic fractures is carried out through the processing of Eq. (18) and Eq. (19). If no water or gas injection from the hydraulic fracture, we have qw Γ ¼ qw F and qg Γ ¼ qg F . When meshing, the discrete fracture is treated as line elements, without breadth, which are edges of area elements, seen Fig. 7b. As can be seen from Fig. 7a that the normal vector of the hy­ draulic fracture surfaces nΓ is the reverse of the normal vector of shale matrix inner surfaces nm . Therefore, Eq. (18) and Eq. (19) is eliminated when we superpose weak equations of shale matrix and hydraulic fractures, i.e. Eq. (27) and Eq. (34).

4 � 10

∅m ∅F

1.714 � 10

9

Pa

1

Water viscosity Oil compressibility coefficient

correctness of our model. Firstly, it can be seen from Eq. (16) and Eq. (17) that if the gas compressibility factor cg is replaced by oil compressibility factor co , and the Knudsen diffusion coefficient is neglected, the gas-water flow model is simplified to the oil-water flow model. Therefore, the experimental results of water flooding profile in previous literatures can be used to validate our model. Second, if the initial water saturation in the whole simulated area is set lower than the irreducible water saturation, the gas pressure solved by our gas-water flow model should be the same with the single gas flow model. Thus, we can fully verify our model from both the water saturation and pressure.

4. Model verification Unlike the gas-water flow in a pipe, such as the riser (Xu et al., 2019), it is hard to detect the gas-water velocity profile and water saturation in the nano to micro meter pores of shale. Therefore, two special cases, obtained by simplifying our model, are performed to verify the 7

S. Wei et al.

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

Fig. 8. The physical model and the comparation of meshing between MsFEM and FEM-FVM (our model).

Fig. 10. The physical model and meshing used in production simulation.

Fig. 9 shows the water saturation profiles after injecting different amount of water. It can be seen from Fig. 9 that our simulation results show a good agreement with the experimental results and the MsFEM results. In the experimental results, it is noticeable that the outer boundary of the water distribution is half circle (the dash line) from the tip of the fracture after injecting 0.5 pore volume (PV) water, which is consistent with our simulation results. However, the outer boundary of the water saturation distribution in MsFEM results is triangle (the dash line). Therefore, we know that the hybrid numerical method coupling FVM and FEM can solve the two-phase flow better than the simply modified FEM method.

Fig. 9. The comparation of water saturation profiles among experimental re­ sults (a), MsFEM results (b) and our FEM-FVM results (c). Both the experimental results and the MsFEM results are from Zhang et al. (2017).

4.1. Verification through oil-water flow simulation The water saturation solved by our model is compared with the Multi-scale Finite Element Method (MsFEM) method proposed by Zhang et al. (2013) and the experimental results (Zhang et al. 2013, 2017). The water flooding experiment in oil reservoirs is much easier to be carried out than that in gas reservoir. Besides, there are lots of numerical method in solving the oil-water flow in fractured porous media (Fumagalli and Scotti, 2013; Hauge and Aarnes, 2009; Hoteit and Fir­ oozabadi, 2008). In this section, the gas compressibility factor cg is replaced by oil compressibility factor co , and the Knudsen diffusion co­ efficient is zero. The simulation parameters are given in Table 1. In order to accelerate the injection process, Zhang’s experiment chooses a highly permeable media. The physical model and meshing results are shown in Fig. 8. And the simulation parameters are from Zhang et al. (2017), which is shown in Table 1. The relative permeability in matrix and fracture are respectively krw ¼ sw and kro ¼ 1 sw . The capillary is neglected both in the matrix and the fracture. The initial oil pressure is set to 1 MPa, and the initial oil saturation is 1. A constant injection rate 0.01 PV/s is applied to the water injection well, and the oil well is produced at a constant pressure of 0.1 MPa.

4.2. Verification through gas pressure As we have mentioned that once the initial water saturation in our gas-water flow model is set lower than the irreducible water saturation, the gas-water flow model should get the same result with single gas flow model, i.e. the single gas flow is the special case of gas-water flow. Our model is solved with our own program developed in MATLAB. And we just set the initial water saturation lower than the irreducible water saturation in our program. The mesh shown in Figs. 8c and 10 is generated using COMSOL. And the purpose is to avoid the effect of grids on the results of COMSOL and our program. The single gas flow model is given in Eq. (36), and it is solved by COMSOL (COMSOL, 2016). And Eq. (37) is employed to calculate the gas production in our simulation process.

8

S. Wei et al.

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

"

Table 2 Simulation parameters. Parameters km

kF

LF

μw

φctg 1 swr

Values 2

∅F w2f =12

Hydraulic fracture permeability

63.75 m

Half-length of hydraulic fracture

0.23 mPa s

Water viscosity

0.02 mPa s

Gas viscosity

0.03

∅f

0.004

Total porosity of the naturally fractured shale matrix Porosity of natural fractures

lmax

10 m

Maximum length of natural fractures

lmin

0.1 m

Minimum length of natural fractures

α

30�

Natural fracture azimuth

θ

30�

β

1 � 10

∅F

0.25

cw

4 � 10 10 Pa 1 1 � 10 11 Pa 1 5 � 10 11 Pa 1 8.314 J/mol/ K 333 K

cp

m

cp

F

R T Dk pL

swr

m

7

Kλrg swr r2 pg ¼qg

φð1 Ω

� sw Þ pini g

pg

�. . R T þ ð1

φÞ Cini μ

�i Cμ Vst dΩ

Formation temperature Knudsen diffusion coefficient Langmuir pressure

0.1

Equilibrium constant of adsorption/desorption

0.645

Irreducible water saturation in shale matrix Residual gas saturation in shale matrix

swr

F

0

Irreducible water saturation in hydraulic fractures

0

Residual gas saturation in hydraulic fractures

sini w

0.4

Initial water saturation

0.9

Water saturation in the highly saturated area

pini w

30.12 MPa

Initial water phase pressure

(37)

ini where Cini μ is the initial concentration of adsorbed gas. Cμ and Cμ are calculated using Eq. (14). Vst is the molecular volume of methane in the standard condition, and it is 0.0224 m3/mol in this paper. Wg is the accumulated gas production. The physical model used in this section is given in Fig. 10. And the constant pressure boundary is applied to the wellbore. The parameter values are from an actual well in Sichuan and are given in Table 2. The relative permeability and capillary pressure curves used in this section are shown in Figs. 3 and 4. Fig. 11 is the production comparation between our gas-water model solved by FEM-FVM method and the single gas flow model solved by COMSOL. Fig. 12 is the pressure comparation in the monitoring point. From Figs. 11 and 12, we can see that the gas production and gas

Pore compressibility coefficient of proppant pack inside of hydraulic fractures Gas constant

0

shigh

!

Pore compressibility coefficient of matrix

m

F

Wg ¼

The proportionality coefficient of natural fracture aperture Porosity of the proppant pack inside of hydraulic fractures Compressibility coefficient of water

sgr sgr

Z h

Natural fracture dip angle

13

9.906 � 10 m2/s 5 MPa

Ke

1 φ

# MKe C2μs R2 T 2 ∂pg �2 Cμs RT þ Ke pg pg ∂t

(36)

Overall permeability of the naturally fractured shale matrix Hydraulic fracture aperture

∅m

μg

m

!

Explanations

2.26 � 10 mD 1 mm

wF

!

Fig. 12. The comparation of gas pressure in the monitoring point.

Fig. 11. The production comparation between our model and single gas flow model. 9

S. Wei et al.

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

Fig. 13. The geometry of a quarter of the actual well.

pressure simulated using our gas-water flow model shows a good agreement with the single gas flow model solved by COMSOL. And the maximum relative error between our model and single gas flow model is 2.67%, which indicates that our model can be applied to simulate the distinctive phenomenon that there are single-gas phase flow and gaswater flow regions co-exist in the shale gas reservoirs. 5. Case study The amount of injected water per meter formation depth can be calculated by the following equation in this paper. Z � Inj ¼ φ shigh sini (38) w dΩ Ω

where Inj is the amount of injected water. shigh is the water saturation in the highly saturated region. Changning oil field, locating in the southern part of the Sichuan basin, is the demonstration area of shale gas exploitation in China (Chen et al., 2015). Changning oil field firstly achieves the commercial development of shale gas in China (Wang et al., 2016). A multi-stage fractured horizontal well from Changning oilfield is selected for the case study. The actual well, which is fractured 22 stages, have be injected into fracturing fluid for 48,356 m3, and the thickness of the reservoir is 52 m. The hydraulic fracture data is interpreted from micro-seismic data. The input data for the model simulation is given in Table 2. The relative permeability and capillary pressure curves are shown in Figs. 2 and 3, and the specific values are given in Table 2. Assuming that the fracturing fluid is equally spaced around these frac­ tures, we build the following physical model of a quarter of the actual well, seen Fig. 13. Assuming that shigh is a constant in the highly saturated region, which is set to 0.9 in this paper. Therefore, Eq. (38) can be written as:

Fig. 15. The history match and forecast of gas production.

Fig. 16. The history match and forecast of water production.

Z Inj ¼

φ shigh Ω

� sini w dΩ ¼ φ shigh

� sini w Ah

(39)

where Ah is the area of the highly saturated region. Thus, we calculated the area of the highly saturated region around each fracture is 700 m2. In this paper, the highly saturated area is assumed as an ellipse. Thus, the major axis and short axis are respectively 70 m and 10 m. The bottom hole pressure (BHP) of the actual well is shown in Fig. 14, and it is well fitted by a natural logarithmic function of pro­ duction time. The fitting function of BHP is the boundary condition in

Fig. 14. The bottom hole pressure of the actual well. 10

S. Wei et al.

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

Fig. 17. The distribution of gas phase pressure after 40 days production.

Fig. 18. The distribution of water saturation in matrix pores after 40 days production.

Fig. 19. The distribution of gas phase pressure after 340 days production.

Fig. 20. The distribution of water saturation in matrix pores after 340 days production.

the wellbore, and the stable BHP is set as 14 MPa. Therefore, the boundary condition in the wellbore is given in Eq. (40). And the other outer boundaries are set as nonflow boundaries. � ð 3:582 ln t þ 27:733ÞMPa; ​ ​ BHP > 14MPa BHP ¼ (40) 14MPa; ​ ​ ​ ​ ​ ​ BHP < 14MPa

correctness of our model and the program. In Fig. 16, because the fracturing fluid flow in the hydraulic fracture is accelerated by the closure of the hydraulic fracture, the simulated water production in the early flowback stage is smaller than the actual data. As the production time goes on, the effect of the fracture closure on the production weakens, and the difference between simulation results and the field data decreases. After 200 days production, the accumulated water pro­ duction is 5273.988 m3, which is 10.91% of the amount of the injected fracturing fluid. And the accumulated water production is almost stable. The water production rate is 1.66 m3/d after 200 days production and 0.0013 m3/d after 300 days production. However, the gas production still remains a high rate. The fracturing fluid is concentrated near

In Changning oilfield, China, there is no separate flowback process after the shale gas well is completed. We can obtain the gas and water production from the production record log. The comparison between the simulation results and the field data are shown in Fig. 15 and Fig. 16. It can be seen from Fig. 15 that the simulated gas production shows a good agreement with the actual gas production, which also verified the 11

S. Wei et al.

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

hydraulic fractures, and the flowback of the fracturing water enhances gas phase permeability near hydraulic fractures. Figs. 17–20 are distributions of gas phase pressure and water satu­ ration during the production. It can be seen from Fig. 17 that the gas phase pressure experiences a large pressure drop in the vicinity of hy­ draulic fractures. And Fig. 18 shows that the fracturing fluid is concentrated near hydraulic fractures, which indicates that the gas relative permeability in the vicinity of hydraulic fractures is much smaller than the areas away from hydraulic fractures. The highly satu­ rated area gets an extra “skin factor” for gas flow due to the high satu­ ration of the fracturing fluid, which increases the gas pressure gradient near the fracture. It can be seen from Fig. 19 that the gas phase pressure experiences a large pressure drop near the wellbore, especially in the vicinity of hy­ draulic fractures. Gas flow also occurs in the area 150 m away from the highly saturated region. The pressure drop in the highly saturated region reaches 16 MPa, which is obviously higher than the unsaturated region. In Fig. 20, the water saturation in the highly saturated area is almost stable. The water saturation in the shale matrix near the fracture is as high as 0.85, and the water saturation distribution decreases rapidly to irreducible water saturation from fractures to the outside. The reason is that most of the injected fracturing fluid flows from the shale matrix to the wellbore through hydraulic fractures. Before the water saturation exceeds the irreducible water saturation, the fracturing fluids injected dooms to be trapped in the shale matrix, which occupies most of the large pore channels in shale matrix. There­ fore, the gas relative permeability is only 0.113 under the irreducible water saturation. For the reason mentioned above, it is of great mea­ surement for improving the gas production to decrease the irreducible water saturation of shale matrix, such as changing the rock surface properties through injecting surface-active agent.

Fig. 22. The influence of highly saturation region on gas production.

6. Flowback and production sensitive studies In this section, a stage of the hydraulic fractures is analyzed, which consists of two clusters, i.e. two main fractures. The permeability and porosity of the shale matrix are respectively 1 � 10 2 mD and 0.03 in this section. Values of other parameters are given in Table 2. The relative permeability curve and capillary pressure curve are from Figs. 3 and 4. The fracture length is 80 m, and the fracture spacing between two fractures is 20 m. The initial water phase pressure is 30 MPa, and the borehole pressure of water phase is fixed into 20 MPa. The physical model and boundary conditions are shown in Fig. 21. For the sake of simplifying the grids and improving the convergence of our program, the highly saturated area is simplified from oval to rectangle.

Fig. 23. The influence of highly saturation region on water production.

Fig. 24. The influence of water injection volume on gas production.

Fig. 21. The geometry of this model used in production sensitive study. 12

S. Wei et al.

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

Fig. 27. The influence of matrix permeability on gas production. Fig. 25. The influence of water injection volume on water production.

production. If we only simulate the gas-water flow region or simulate the flowback process and production process separately, neither the flow­ back of fracturing fluid nor the gas production can be accurately predicted. For the scenario shigh ¼ 0:7, most of the injected water retains in the matrix pores as the supplement of irreducible water since the irreducible water saturation (0.645) is little smaller than 0.7, and the flowback water is only 0.051 percent of the injected water after 500 days pro­ duction. For other two scenarios shigh ¼ 0:8 and shigh ¼ 0:9, ratios of flowback water to injected water are respectively 7.1% and 21.76%. Thus, we known that the larger the amount of injected water, the larger the ratio of flowback water to injected water. In Fig. 24, the larger the amount of injected water, the lower the gas production. Comparing the accumulated gas production at 500 days, the accumulated gas produced increases by 6.22% and 14.62% for shigh from 0.9 (base case) to 0.8 and 0.7, respectively. Figs. 24 and 25 indicate the relation that the larger the ratio of flowback water to injected water, the lower the gas production. However, the retained water in the shale matrix for shigh ¼ 0:7, shigh ¼ 0:8 and shigh ¼ 0:9 are respectively 1853.28 m3, 2336.51 m3 and 2509.1 m3, which means the higher the shigh , the more the water retained in the shale matrix, and the lower the relative permeability of gas. That is the reason why the accumulated gas production decreases with the increase of shigh . Since the irreducible water saturation is higher than the initial water saturation, the retained water volume in the shale matrix is determined by the water injection volume and the irreducible water saturation of shale matrix. The superior limit of flowback ratio is ðshigh swm Þ=ðshigh sini w Þ, where swm is the irreducible water saturation. Due to

6.1. The influence of the highly saturated region Firstly, we analyzed the influence of highly saturation region. In this paper, we consider two scenarios of water saturation distribution with the same water injection volume. 1) Homogeneous water saturation: the whole simulated area is given a homogeneous initial water saturation 0.6911. 2) Inhomogeneous water saturation: the initial water saturation of the highly saturation region in Fig. 21 is set to 0.9, the water satu­ ration in other regions is set to 0.6452, which is equal to irreducible water saturation. Then, we analyzed three scenarios of the water satu­ ration in the highly saturated region namely 0.9 (base case), 0.8, 0.7 (little higher than the irreducible water saturation), and the initial water saturation is 0.4. In this section, the area of the simulated region is a constant. Therefore, the three scenarios of the highly saturated regions represent three scenarios of water injection volume, i.e. 2808 m3, 2246.4 m3 and 1684.8 m3. Fig. 22 and Fig. 23 compare the difference between whether considering highly saturated region (i.e. inhomogeneous water satura­ tion) or not (i.e. homogeneous water saturation). The amount of injected water is the same in the two scenarios, i.e. 1486 m3. In Fig. 22, without considering highly saturated region, the gas production increase greatly than that of considering highly saturated region. In Fig. 23, accumulated water productions of two scenarios are respectively 565.17 m3 and 0.96 m3. We can clearly see that it is unreasonable to give a homogeneous initial water saturation in the whole simulated area, of which the accumulated water production is nearly 0, which is conflict with the field observation, i.e. the long-time two-phase flow period during the

Fig. 26. Influences of fracture dip angle, fracture azimuth and fracture porosity on shale permeability. 13

S. Wei et al.

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

Fig. 28. The influence of matrix permeability on water production.

Fig. 30. The influence of Langmuir pressure on water production.

the heterogeneity of shale reservoir, the irreducible water saturation varies a lot among different wells, which is part of the reason why the flowback ratio varies a lot. In Fig. 23, for the scenario of inhomogeneous water saturation, the initial water saturation is 0.6452 (the irreducible water saturation), shigh ¼ 0:9, and the accumulated gas production is 1.84 � 106 m3. The scenario shigh ¼ 0:9 in Fig. 25 have the initial water saturation 0.4, and the accumulated gas production is 2.51 � 106 m3. Comparing the accumulative water production of these two scenarios, we know that the higher the initial water saturation, the lower the gas production.

several parameters. When we know the influence of shale permeability on gas and water production, influences of fractal parameters on gas and water production are also clear. We analyzed three scenarios of the matrix permeability, i.e. the absolute permeability of the naturally fractured shale, namely 1 � 10 3 mD, 1 � 10 2 mD (base case) and 5 � 10 2 mD. From Fig. 27, it is clearly seen that the matrix permeability has a great influence on the gas production. And the larger the matrix permeability, the larger the gas production. The influence of matrix permeability on water production is shown in Fig. 28. For the scenario km ¼ 5 � 10 2 mD, the water production rate (the slope of accumulated water production curve) declines to 0 m3/d after 100 days production. However, for the scenario km ¼ 1 � 10 3 mD, the water production rate is almost the same as the early production stage even after 500 days production. Form Fig. 28, we clearly see that there is an upper limit of flowback ratio. The max amount of flowback water when km ¼ 5 � 10 2 mD is 616 m3, and the flowback ratio is 21.89%. From Figs. 26–28, we can draw the conclusion: the larger the natural fracture length (both the minimum length and the maximum length), the larger the natural fracture porosity, the larger the gas production, and the longer the water production time.

6.2. The influence of natural fracture parameters The influence of natural fractures on shale permeability can be described using Eq. (2) – (4). In the fractal model, shale permeability is determined by the length and porosity of natural fractures, the mean facture azimuth and facture dip angle. The larger the maximum length and the minimum length of natural fractures, the larger the natural fracture permeability. Fig. 26 shows influences of fracture dip angle, fracture azimuth and fracture porosity on shale permeability. In Fig. 26, fractal parameters of the base case are: lmax ¼ 10m, lmin ¼ 0:1m, α ¼ � � 30 , θ ¼ 30 , ∅ ¼ 0:03, ∅f ¼ 0:004, dE ¼ 2. As we can see from Fig. 26, shale permeability is determined by

6.3. The influence of gas desorption In this paper, the adsorbed gas is considered using isothermal Langmuir equation, i.e. Eq. (14). In the base case, the Langmuir pressure is 5 MPa. We analyze three scenarios of Langmuir pressure namely 4 MPa, 5 MPa (base case) and 6 MPa. Fig. 29 shows the simulation results of accumulated gas production. Comparing the accumulated gas production at 500 days, the cumulative gas production increases by 7.57% and 7.75% for the Langmuir pressure from 4 MPa to 5 MPa (base case) and 6 MPa, respectively. The reason is that a larger Langmuir pressure indicates a larger maximum adsorption capacity of the shale matrix (see Eq. (11)), and the gas adsorption effect is weaker (Cao et al., 2016). Therefore, the free gas in matrix pores will get more supplement with a higher Langmuir pressure. It is shown in Fig. 30 that the higher the Langmuir pressure, the larger amount of the flowback water. The reason is that the larger the Langmuir pressure, the greater the mass supplement of the desorption of the adsorbed gas, thus the higher the gas and water pressure, and the larger the water pressure gradient at the same production time. There­ fore, the area of immobile water between the two fractures is smaller with larger Langmuir pressure.

Fig. 29. The influence of Langmuir pressure on gas production. 14

S. Wei et al.

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

Fig. 31. Geometries of fracture networks (left), and the coordinates of fracture ends (right); The size of the simulated area is the same as Fig. 21 shows.

Fig. 32. The influence of fracture network on gas production.

Fig. 33. The influence of fracture network on water production.

6.4. The influence of fracture networks In this section, six scenarios of the fracture networks are compared, and the physical models can be seen in Fig. 31. The fracture networks 1–3 show different cluster spacing of a stage of hydraulic fractures, which are respectively 20 m, 10 m and 5 m. The fracture networks 4–6 show three different complex fracture networks. In Fig. 32, the fracture network 1 has the smallest accumulated gas production, which is obviously smaller than other fracture networks. The fracture network 6 has the highest accumulated gas production. And the fracture network 6 has the cluster spacing 10 m and better con­ nectivity than fracture network 2. In our previous study of single gas flow (Wei et al., 2019a,b), the fracture density, i.e. the ratio of total fracture length to the simulated area, has a great influence on gas pro­ duction. However, when we consider the flowback of fracturing fluids, the gas production is slightly improved by the improvement of fracture density. The main reason is that lots of fracturing fluids retains in the vicinity of hydraulic fractures (see Fig. 32), especially near fracture junctions, which greatly reduces the gas phase permeability in these places. It can be seen from Fig. 33 that the larger the cluster spacing, the

Fig. 34. The influence of fracture network on water production.

15

S. Wei et al.

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

Fig. 35. The distribution of the water saturation after 500 days production in the 6 fracture networks.

larger the accumulated water production, i.e. the larger the flowback ratio. The accumulated water production of the three complex fracture networks are close. Fig. 34 shows the water production rate of 6 fracture networks, of which the time axis is logarithmic. It can be seen from Fig. 34 that complex fracture networks have the obviously advantage in the early time flowback rate. However, due to the fracture interference, the flowback rate decreases rapidly and the overall flowback ratio is smaller than the simple fracture network. In the actual hydraulic fracturing process, the hydraulic fracture networks are absolutely different in different wells, even with exactly the same fracturing measurements. Therefore, the fracture interference effect can contribute to the variation of the flowback ratio among different wells. In Fig. 35, water saturation in the corner of crossed fractures is up to 0.9, which is higher than other places near hydraulic fractures. The reason is that these corners are subject to greater fracture interference, which again emphasizes the importance of reducing the irreducible water saturation near hydraulic fractures. Combining results in Figs. 33 and 35, we know that the fracture interference of fracture network 1 is smaller than other fracture networks. The smaller the cluster spacing, the larger the fracture interference, and the more the fracturing fluids retained in the formation.

(2)

(3)

(4) (5)

7. Conclusion

formation into highly saturated region and initial water saturated region satisfies well with the field observation. The flowback ratio is limited by the irreducible water saturation and fracture interference. Due to the space variation of irreduc­ ible water saturation (higher than the initial water saturation) and the difference of hydraulic fracture networks among different wells, the flowback ratio varies a lot among different wells and fails to be linked with the gas production. The amount of retained water in the shale matrix is negatively related to the gas production. The retained water forms the pos­ itive “skin factor” for the gas flow in the vicinity of hydraulic fractures, which calls measurements to reduce the irreducible water saturation near hydraulic fractures. Due to the fracture interference, the fracturing fluid will be concentrated in the hydraulic fracture junctions, which will seriously reduce the gas production. The larger the natural fracture length, the larger the natural fracture porosity, the larger the gas production, and the longer the water production time. In the range of 0–90� , the larger the fracture dip angle and the smaller the fracture azimuth, the larger the gas production and the longer the water production time. The larger the Langmuir pressure, the higher the gas production, and the slightly bigger the flowback ratio.

8. Future works

In this paper, a 2D gas-water flow model coupling shale matrix and discrete fractures is developed, which is solved by the hybrid numerical method (finite element and finite volume method). The proposed model solves two problems which are not answered in previous studies. Firstly, the shale reservoir is constructed by a highly saturation region and an initial saturation region. The initial water saturation is lower than the irreducible water saturation. We find that the field observation of long water production time can be well explained by dividing the formation into highly saturated region (the water saturation is higher than the irreducible water saturation) and initial water saturated region (the water saturation is lower than the irreducible water saturation). Sec­ ondly, it is not the flowback ratio but the amount of retained fracturing fluids should be used to reflect the influence of flowback measures on gas production. The less the fracturing fluids retained in the formation, the larger the gas production. The irreducible water saturation and fracture interference are two main reasons influence the flowback ratio, which are the reasons why the flowback ratio varies a lot among different wells. Based on simulation results, the following conclusions can be drawn:

Due to the reason that the initial water saturation is lower than the irreducible water saturation, once the fracturing fluid is injected into the shale matrix, a great amount of the fracturing fluid is retained in the shale matrix. And the large matrix pores, near hydraulic fractures, are occupied by the fracturing fluid, which forms the positive “skin factor” in the vicinity of hydraulic fracture, and reduces the gas production. Hence the future works we mainly focus on are the following two as­ pects: 1) the gas and water production behaviors in the complex fracture networks; 2) the measurements to reduce the irreducible water satura­ tion near hydraulic fractures. The first work can be done by optimizing the efficiency of our program. The other forces us to add the consider­ ation of chemical reaction in our model. Acknowledgments The authors are grateful to the supports provided by the National Natural Science Foundation of China (Grand No. 51490650), and the Ministry of Science and Technology of the People’s Republic of China (2017ZX05037-004).

(1) The water saturation in the highly saturated region has a strong influence on the flowback ratio. The consideration of dividing the

16

S. Wei et al.

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

Nomenclature The area of the highly saturated region, m2 The gas concentration of adsorbed gas and free gas, mol/m3 The maximum adsorption capacity of shale, mol/m3

Ah Cμ ; C C μs

The initial concentration of adsorbed gas, mol/m3 Gas concentration, mol/m3 Coefficient of water compressibility, 1/Pa Coefficient of gas compressibility, 1/Pa The fractal dimension for fracture length, dimensionless Natural fracture density, m/m2 The dimension of simulated area, dimensionless The coefficient of Knudsen diffusion, m2/s The amount of injected water, m3 The water phase permeability, 1 The gas phase permeability, 1

Cini μ Cg cw cg Df D dE Dk Inj krw krg ~rg k

lmax ; lmin M mw ; mg pc pc m pc F pd pg pw R sw sw m sw F se shigh T

μw μg

Vst Wg

α θ ∅ ∅F

αij

Ke

ρw

λi→j rw

The apparent permeability of gas phase in shale matrix, 1 The maximum and minimum length of natural fractures, m The molecular mass of methane, kg/mol Mass sources of water and gas, kg/m3/s Capillary pressure, MPa Capillary pressure in shale matrix, MPa Capillary pressure in hydraulic fractures, MPa The minimum capillary pressure, MPa Gas pressure, MPa Water pressure, MPa Gas constant Water saturation, 1 Water saturation in matrix pores, 1 Water saturation in hydraulic fractures, 1 The effective water saturation, 1 The water saturation in the highly saturated region, 1 Formation temperature, K Viscosity of fracturing fluids, mPa⋅s Gas viscosity, mPa⋅s The molecular volume of methane in the standard condition, m3/mol The accumulated gas production, m3 Natural facture azimuth, Natural facture dip angle, � Total porosity of the naturally fractured shale matrix Porosity of the proppant pack inside of hydraulic fractures Shape factor The equilibrium constant of adsorption/desorption, 1 Density of fracturing fluids, kg/m3 The relative mobility of water from node i to node j, 1/(mPa⋅s)

Appendix A. Mathematic model of two-phase flow The mass conservation equations of water and gas in porous media are:

∂ðφρw sw Þ ∂t ∂ φ ρg s g ∂t

�

r⋅ðKρw λrw rpw Þ ¼ mw

(A-1)

� r⋅ Kρg λrg rpg ¼ mg

(A-2) (A-3)

sw þ sg ¼ 1

where φ is the total porosity of the naturally fractured shale matrix. λrw is the relative mobility of water, and λrw ¼ krw =μw . K is the permeability matrix. pw and Sw are respectively the pressure and saturation of water phase. pg and sg are respectively the pressure and saturation of gas phase. λrg is the ~rg =μ . k ~rg is the apparent permeability of gas phase when we consider different flow mechanisms. ρ and ρ are relative mobility of gas, and λrg ¼ k g

respectively the water density and gas density. mw and mg are respectively mass sources of water and gas. Considering the water and rock as slightly compressible material, thus we have: � �� cw ðpw pini Þ ¼ ρini ρw ¼ ρini pini w e w 1 þ cw pw 17

w

g

(A-4)

S. Wei et al.

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

φ ¼ φini ecp ðpw

pini Þ

� ¼ φini 1 þ cp pw

pini

��

(A-5)

where cp is the pore compressibility coefficient the superscript -ini means the initial state. With Eq. (A-4) and Eq. (A-5), the time derivative and convection terms of Eq. (A-1) can be expressed as:

∂ðφρw Sw Þ ∂pw ∂Sw ini ¼ ρini þ ρw φ w φ ctw Sw ∂t ∂t ∂t � cw ðpw r⋅ðKρw λrw rpw Þ ¼ Kr⋅ λrw ρini w e

(A-6) pini Þ

ini � ecw ðpw p Þ r rpw ¼ Kr⋅ λrw ρini w cw

!

�� � � � � 1 þ cw pw pini ini 2 ¼ Kr⋅ λrw ρini ¼ Kr⋅ λrw ρini w r w rpw ¼ Kρw rλrw ⋅rpw þ λrw r pw cw � � dλrw 2 ¼ Kρini rsw ⋅rpw þ λrw r2 pw ¼ Kρini w w λrw r pw dss

(A-7)

Where ctw ¼ cw þ cp . In Eq. (A-6), we neglect the quadratic infinitesimals containing cp cw . In Eq. (A-7), rSw ⋅rpw and a rpw ⋅rpw re also neglected. Substitute Eq. (A-6) and Eq. (A-7) into Eq. (A-1), we get the flow equation of the slightly compressible water in the slightly compressible porous media. ini ρini w φ ctw Sw

∂pw ∂Sw þ ρw φ ∂t ∂t

(A-8)

2 Kρini w λrw r pw ¼ mw

Divide ρini w on both sides of the equation at the same time. Compared with the time derivative of the pressure, the time derivative of the water saturation is an infinitesimal. Thus, we get: φini ctw sw

∂pw ∂sw þφ ∂t ∂t

(A-9)

Kλrw r2 pw ¼ qw

where qw ¼ mw =ρini w . The state equation of ideal gas is:

ρg ¼

pg M RT

(A-10)

Thus, we have the following relation:

ρg ¼ ρini g

p pini

(A-11)

With Eq. (A-10) and Eq. (A-11), the time derivative term of Eq. (A-2) can be expressed as: ! � ρini ∂ φ ρg s g ∂pw ∂sg g ini ¼ φ ini þ φ cp pg sg þ φρg ∂t pg ∂t ∂t In Eq. (A-12), the time derivative of capillary pressure is neglected. The convection term of Eq. (A-2) can be expressed as: � � � � r ⋅ Kρg λrg rpg ¼ K r ρg λrg ⋅ rpg þ ρg λrg r2 pg

(A-12)

(A-13)

It is known that ρg λrg is the function of gas pressure and gas saturation. rsg ⋅rpg and rpg ⋅rpg are also neglected. Therefore, Eq. (A-13) can be expressed as: � � � dpc 2 r ⋅ Kρg λrg rpg ¼ Kρg λrg r2 pg ¼ Kρg λrg r2 pw þ r sw (A-14) dsw Substitute Eq. (A-12) and Eq. (A-14) into Eq. (A-2), we get: ! � � ρini ∂pw ∂sg dpc 2 g φ ini þ φini cp ρg sg r sw ¼ mg þ φρg Kρg λrg r2 pw þ pg ∂t ∂t dsw

(A-15)

Divide the gas density on both sides of Eq. (A-15) at the same time. ! � � ρini ∂pw ∂sg dpc 2 g ini 2 φ s þ φ c p þ r s þ φ Kλ r ¼ qg p g rg w w ρg pini ∂t ∂t dsw g

(A-16)

where qg ¼ mg =ρg . From Eq. (A-10) and Eq. (A-11), we know that:

ρini d ρg 1 g ¼ ¼ ¼ cg pg ρg pg ρg pini g

(A-17)

18

S. Wei et al.

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

where cg is the coefficient of gas compressibility. Substitute Eq. (A-17) into Eq. (A-16), we obtain the gas flow equation in the same form as Eq. (A-9). � � ∂pw ∂sg dpc 2 φini ctg sg r sw ¼ qg þφ Kλrg r2 pw þ ∂t ∂t dsw

(A-18)

where ctg ¼ cg þ cp . Here, we make an approximation, i.e. φcg þ φini cp ¼ φini ðcg þ cp Þ. By adding Eq. (A-9) and Eq. (A-18), we can get the following equation: φini ctg sg þ ctw sw

� ∂pw ∂t

Kλt r2 pw

Kλrg

dpc 2 r sw ¼ qg þ qw dsw

(A-19)

Eq. (A-9) and Eq. (A-19) constitute the mathematical model of gas-water two-phase flow in porous media. And Eq. (A9) and Eq. (A17) can be assembled into the matrix form: 3 2 dpc m � � ini � � � �� � � p_w qg m þ qw m 0 φ ctg sg þ ctw sw 6 Kλt Kλrg ds 7 r2 pw þ4 ¼ (A-20) w 5 ini 2 s _ qw m φ ctw sw φ r sw w Kλrw 0 Because the width of hydraulic fractures is much smaller than the length of the crack and the proppant is filled in the crack, the Darcy flow equation is used in the hydraulic fractures. And the velocity is seen as a constant along the cross section of hydraulic fractures. Therefore, the flow equation in the hydraulic fractures is: 8 9 2 3 #� � " dpc f � > �> � = � < φini c s þ c s � 2 Kλ Kλ 0 t rg qg f þ qw f p _ tg g tw w p r 6 7 F w w dsw 5 wF � þ4 ¼ (A-21) 2 ini qw f s_w > r sw > φF ctw sw af φF : ; Kλrw 0 where φF is the porosity of the proppant pack inside of hydraulic fractures. wF is the fracture aperture. Appendix B. FVM discretization of the flow part of mass conservation equation Consider the mass conservation in the subdivided grid, Eq. (14) can be expressed as: Z fi ¼ ðρw Kλrw rpw Þ⋅ndΓ dcþcb

(A-22)

Fig. A-1 shows that one of the FEM grids can be subdivided into three control volume. The point c is the barycenter of the given triangular grid. And the coordinates of each node are given in the cartesian coordinate system. Thus, the unit normal vector of curve d→c and c→b is: nd→c ¼

ðyc

yd ; xd jdcj

xc Þ ðyb ; ​ nc→b ¼

yc ; xc jcbj

xb Þ

(A-23)

In order to couple FVM with FEM, the same shape function used in both FVM and FEM, i.e.: �� �T � pw ¼ Ni ; Nj ; Nm Pi ; Pj ; Pm Substitute Eq. (A-23) and Eq. (A-24) into Eq. (A-22), we get the discrete form of Eq. (A-22): � �T � �T �� �� ∂pw ∂pw 1 ∂pw ∂pw fi ¼ ρw Kλrw ¼ ρw Kλrw ; ⋅ yb yd ; xd xb ; ⋅ ym yj ; xj xm 2 ∂x ∂y ∂x ∂y � �1 0 0 1 ∂Ni ∂Ni ∂Nj ∂Nm ⋅ B ∂x ∂x ∂x ∂x C pi B CB C 1 CB C ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ¼ ρw Kλrw B B � � C@ pj A 2 @ ∂Ni ∂Ni ∂Nj ∂Nm A pm ⋅ ∂y ∂y ∂y ∂y As we know that Ni þ Nj þ Nm ¼ 1, we get another relation among the three shape functions: � rNi ¼ rNj þ rNk Substitute Eq. (A-26) into Eq. (A-25), we get � �� � � � 1 ∂Ni ∂Nj ∂Ni ∂Nj ∂Ni ∂Nm ∂Ni ∂Nm fi ¼ ρw Kλrw þ pj pi þ þ ðpm ∂x ∂x ∂y ∂y ∂x ∂x ∂y ∂y 2

� pi Þ

(A-24)

(A-25)

(A-26)

(A-27)

As we can see from Eq. (A-27) that the water mass flux among control volumes is decided by nodal pressure. And the water mass flux, form control volume i to control volume j and m, is split into two parts, which respectively expressed as the function of pressure difference between node i and node j, node i and node m. As we have mentioned before, the flow direction of water must be selected according to both the water pressure and water saturation. Therefore, the relative mobility in Eq. (A-27) is merged with nodal pressure difference.

19

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

S. Wei et al.

2 6 1 6 f i ¼ ρw K 6 4 2

�

�

� 3 pi þ 7 7 7 � � 5 ∂Ni ∂Nm i→m ∂Ni ∂Nm λrw þ ðpm pi Þ ∂x ∂x ∂y ∂y

λi→j rw

∂Ni ∂Nj ∂Ni ∂Nj þ ∂x ∂x ∂y ∂y

pj

(A-28)

In Eq. (A-28), the relative mobility of water is chosen using “upwind” scheme (Patankar, 1980). ( � i λi→j rw sw ; ​ ​ pi > pj j→i � ¼ λ ¼ λi→j rw rw j λi→j rw sw ; ​ ​ pi < pj

(A-29)

Similarly, we get the water mass flux from control volume j to control volume i and m, and the water mass flux from control volume m to control volume i and j. � � 2 � 3 ∂Nj ∂Ni ∂Nj ∂Ni λj→i þ pi pj þ rw 7 6 ∂x ∂x ∂y ∂y 1 7 6 f j ¼ ρw K 6 (A-30) 7 � � 4 5 2 � ∂ N ∂ N ∂ N ∂ N j m j m j→m λrw þ pm pj ∂x ∂x ∂y ∂y �

2 6 1 6 fm ¼ ρw K 6 4 2

λm→i rw

� λm→j rw

�

∂Nm ∂Ni ∂Nm ∂Ni ðpi þ ∂x ∂x ∂y ∂y ∂Nm ∂Nj ∂Nm ∂Nj þ ∂x ∂x ∂y ∂y

3 pm Þþ

� pj

pm

7 7 7 � 5

(A-31)

Combine Eq. (A-28 - A-31) into a vector, i.e. the mass flux vector of the given grid in Figure A-1: 0 1 10 1 0 Tii Tij Tim pi f ! @ iA 1 f ¼ fj ¼ ρw K⋅@ Tji Tjj Tjm A@ pj A 2 Tmi Tmj Tmm fm pm

(A-32)

Fig. A-1. The subdivided grid of FEM.

Appendix D. Supplementary data Supplementary data to this article can be found online at https://doi.org/10.1016/j.petrol.2019.106694.

References

In: SPE Western North American Regional Meeting. https://doi.org/10.2118/ 144057-MS. Brooks, R., Corey, T., 1964. In: Hydraulic Properties of Porous Media. Hydrology Papers, vol. 24. Colorado State University, pp. 1–27 (3). Bunch, J.R., Hopcroft, J.E., 1974. Triangular factorization and inversion by fast matrix multiplication. Math. Comput. 28 (125), 231–236. https://doi.org/10.2307/ 2005828. Cao, P., Liu, J., Leong, Y.K., 2016. A fully coupled multiscale shale deformation-gas transport model for the evaluation of shale gas extraction. Fuel 178, 103–117. https://doi.org/10.1016/j.fuel.2016.03.055. Carpenter, C., 2014. Characterization of hydraulic-fracture geometry in shale-gas reservoirs. JPT (J. Pharm. Technol.) 66 (11), 125–128. https://doi.org/10.2118/ 1114-0125-JPT. Chang, J., Yortsos, Y., 1990. Pressure transient analysis of fractal reservoirs. SPE Form. Eval. 5 (1), 31–38. https://doi.org/10.2118/18170-PA. Chelidze, T., Guguen, Y., 1990. Evidence of fractal fracture. Int. J. Rock Mech. Min. 27 (3), 223–225. https://doi.org/10.1016/0148-9062(90)94332-N. Chen, Z.Y., Moh, T.C., Phan, V.C., et al., 2015. Drilling evolution in Changning shale gas block of Central China. In: SPE/IATMI Asia Pacific Oil & Gas Conference and Exhibition, 20-22 October, Nusa Dua, Bali, Indonesia. https://doi.org/10.2118/ 176141-MS.

Akkutlu, I.Y., Efendiev, Y., Vasilyeva, M., et al., 2017. Multiscale model reduction for shale gas transport in a coupled discrete fracture and dual-continuum porous media. J. Nat. Gas Sci. Eng. 48, 65–76. https://doi.org/10.1016/j.jngse.2017.02.040. Alkhlaifat, A.L., Arastoopour, H., 2003. Simulation of water and gas flow in fractured porous media. Biol. Reprod. 209 (1), 201–212. https://doi.org/10.1144/GSL. SP.2003.209.01.17. Andrade, J.J., Oliveira, E., Moreira, A., et al., 2009. Fracturing the optimal paths. Phys. Rev. Lett. 103 (22), 225503. https://doi.org/10.1103/physrevlett.103.225503. Barree, R.D., Mukherjee, H., 1995. Engineering criteria for fracture flowback procedures. In: Low Permeability Reservoirs Symposium. https://doi.org/10.2118/29600-MS. Basquet, R., Jeannin, L., Lange, A., et al., 2004. Gas-flow simulation in discrete fracturenetwork models. SPE J. 7 (5), 378–384. https://doi.org/10.2118/88985-PA. Birdsell, D.T., Rajaram, H., Lackey, G., 2015. Imbibition of hydraulic fracturing fluids into partially saturated shale. Water Resour. Res. 51 (8), 6787–6796. https://doi. org/10.1103/physrevlett.103.225503. Bonnet, E., Bour, O., Odling, N.E., et al., 2001. Scaling of fracture systems in geological media. Rev. Geophys. 39 (3), 347–383. https://doi.org/10.1029/1999rg000074. Brohi, I., Darvish, M.P., Aguilera, R., 2011. Modeling fractured horizontal wells as dual porosity composite reservoirs-application to tight gas, shale gas and tight oil cases.

20

S. Wei et al.

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

Cho, Y., Ozkan, E., Apaydin, O.G., 2013. Pressure-dependent natural-fracture permeability in shale and its effect on shale-gas well production. SPE J. 16 (2), 216–228. https://doi.org/10.2118/159801-PA. Clarkson, C.R., Solano, N., Bustin, R.M., et al., 2013. Pore structure characterization of North American shale gas reservoirs; using USANS/SANS, gas adsorption, and mercury intrusion. Fuel 103 (1), 606–616. https://doi.org/10.1016/j. fuel.2012.06.119. Clarkson, C.R., Williams-Kovas, J., 2013. Modeling two-phase flowback of multifractured horizontal wells completed in shale. SPE J. 18 (4), 795–812. https://doi. org/10.2118/162593-PA. Clarkson, C.R., Williams-Kovas, J., 2013. A new method for modeling multiphase flowback of multi-fractures horizontal tight oil wells to determine hydraulic fracture properties. In: SPE Annual Technical Conference and Exhibition. https://doi.org/ 10.2118/166214-MS. COMSOL, 2016. Multi-physics Modeling Software. COMSOL, Burlington, MA. www. comsol.com. Crafton, J.W., Gunderson, D.W., 2007. Stimulation flowback management: keeping a good completion good. In: SPE Annual Technical Conference and Exhibition. Anaheim. https://doi.org/10.2118/110851-MS. Cugnart, R., Rasoanaivo, S., Williams-Kovacs, J.D., et al., 2017. Early time SRV characterization through flowback analysis: application of Clarkson/WilliamsKovacs technique to Vaca Muerta. In: SPE/AAPG/SEG Unconventional Resources Technology Conference. https://doi.org/10.15530/URTEC-2017-2689844. Darabi, H., Ettehad, A., Javadpour, F., et al., 2012. Gas flow in ultra-tight shale strata. J. Fluid Mech. 710, 641–658. https://doi.org/10.1017/jfm.2012.424. Denney, D., 2009. Evaluating implications of hydraulic fracturing in shale-gas reservoirs. JPT (J. Pharm. Technol.) 61 (8), 53–54. https://doi.org/10.2118/0809-0053-JPT. Ding, D.Y., Farah, N., Bourbiaux, B., et al., 2017. Simulation of matrix-fracture interaction in low-permeability fractured unconventional reservoirs. In: SPE Reservoir Simulation Conference. https://doi.org/10.2118/182608-MS. Ely, J.W., Arnold, W.T., Holditch, S.A.A., 1990. New techniques and quality control find success in enhancing productivity and minimizing proppant flowback. In: SPE Annual Technical Conference and Exhibition. https://doi.org/10.2118/20708-MS. Ezulike, D.O., Dehghanpour, H., Virues, C.J.J., Hawkes, R.V., 2014. A flowback-guided approach for production data analysis in tight reservoirs. In: SPE/CSUR Unconventional Resources Conference. https://doi.org/10.2118/171636-MS. Fan, D., Yao, J., Hai, S., et al., 2015. A composite model of hydraulic fractured horizontal well with stimulated reservoir volume in tight oil & gas reservoir. J. Nat. Gas Sci. Eng. 24, 115–123. https://doi.org/10.1016/j.jngse.2015.03.002. Fumagalli, A., Scotti, A.A., 2013. Numerical method for two-phase flow in fractured porous media with non-matching grids. Adv. Water Resour. 62 (12), 454–464. https://doi.org/10.1016/j.advwatres.2013.04.001. Fung, L.S., Du, S., 2016. Parallel-simulator framework for multipermeability modeling with discrete fractures for unconventional and tight gas reservoirs. SPE J. 21 (04), 1370–1385. https://doi.org/10.2118/179728-PA. Hauge, V.L., Aarnes, J.E., 2009. Modeling of two-phase flow in fractured porous media on unstructured non-uniformly coarsened grids. Transp. Porous Media 77 (3), 373–398. https://doi.org/10.1007/s11242-008-9284-y. Hoteit, H., Firoozabadi, A., 2008. An efficient numerical model for incompressible twophase flow in fractured media. Adv. Water Resour. 31 (6), 891–905. https://doi.org/ 10.1016/j.advwatres.2008.02.004. Hatton, C.G., Main, I.G., Meredith, P.G., 1994. Non-universal scaling of fracture length and opening displacement. Nature 367 (6459), 160–162. https://doi.org/10.1038/ 367160a0. IIk, D., Currie, S.M., Symmons, D., et al., 2010. A comprehensive workflow for early analysis and interpretation of flowback data from wells in tight gas/shale reservoir systems. In: SPE Annual Technical Conference and Exhibition. https://doi.org/ 10.2118/135607-MS. Jiang, J., Younis, R.M., 2015. A multimechanistic multicontinuum model for simulating shale gas reservoir with complex fractured system. Fuel 161, 333–344. https://doi. org/10.1016/j.fuel.2015.08.069. Jiang, J., Younis, R.M., 2016. Hybrid coupled discrete-fracture/matrix and multicontinuum models for unconventional-reservoir simulation. SPE J. 21 (03), 1009–1027. https://doi.org/10.2118/178430-PA. Kazemi, H., Seth, M.S., Thomas, G.W., 1969. The interpretation of interference tests in naturally fractured reservoirs with uniform fracture distribution. SPE J. 4 (9), 463–472. https://doi.org/10.2118/2156-B. Kelley, C.T., 2003. Solving Nonlinear Equations with Newton’s Method. Siam. Kolyukhin, D., Torabi, A., 2012. Statistical analysis of the relationships between faults attributes. J. Geophys. Res.: Solid Earth. 117 (B5), 15–20. https://doi.org/10.1029/ 2011JB008880. Langmuir, I., 1918. The adsorption of gases on plane surfaces of glass, mica and platinum. J. Am. Chem. Soc. 40 (9), 1361–1403. https://doi.org/10.1021/ ja02242a004. Li, L., Lee, S.H., 2008. Efficient field-scale simulation of black oil in a naturally fractured reservoir through discrete fracture networks and homogenized media. SPE Reserv. Eval. Eng. 11 (04), 750–758. https://doi.org/10.2118/103901-PA. Liu, S., Wang, J., He, H., Wang, H., 2018. Mechanism on imbibition of fracturing fluid in nanopore. Nanosci. Nanotechnol. Lett. 10 (1), 87–93. https://doi.org/10.1166/ nnl.2018.2594. Miao, T., Yu, B., Duan, Y., Fang, Q., 2015. A fractal analysis of permeability for fractured rocks. Int. J. Heat Mass Transf. 81 (81), 75–80. https://doi.org/10.1016/j. ijheatmasstransfer.2014.10.010. Moinfar, A., Varavei, A., Sepehrnoori, K., Johns, R.T., 2013. Development of a coupled dual continuum and discrete fracture model for the simulation of unconventional

reservoirs. In: SPE Reservoir Simulation Symposium, 18-20 February, the Woodlands, Texas, USA. https://doi.org/10.2118/163647-MS. Patankar, S.V., 1980. Numerical Heat Transfer and Fluid Flow. Taylor & Francis. Rammay, M.H., Awotunde, A.A., 2016. Stochastic optimization of hydraulic fracture and horizontal well parameters in shale gas reservoirs. J. Nat. Gas Sci. Eng. 36, 71–78. https://doi.org/10.1016/j.jngse.2016.10.002. Reichenberger, V., Jakobs, H., Bastian, P., et al., 2006. A mixed-dimensional finite volume method for two-phase flow in fractured porous media. Adv. Water Resour. 29 (7), 1020–1036. https://doi.org/10.1016/j.advwatres.2005.09.001. Robinson, B.M., Holditch, S.A., Whitehead, W.S., 1988. Minimizing damage to a controlled flowback procedure. SPE J. 40 (6), 753–760. https://doi.org/10.2118/ 15250-PA. Roy, S., Raju, R., Chuang, H.F., et al., 2003. Modelling gas flow through microchannels and nanopores. J. Appl. Phys. 93, 4870–4879. https://doi.org/10.1063/1.1559936. Sahimi, M., 2012. Flow and Transport in Porous Media and Fractured Rock. Wiley.Com, p. 167. Schultz, R.A., Soliva, R., Fossen, H., et al., 2008. Dependence of displacement–length scaling relations for fractures and deformation bands on the volumetric changes across them. J. Struct. Geol. 30 (11), 1405–1411. https://doi.org/10.1016/j. jsg.2008.08.001. Shahid, A.S.A., Wassing, B.B.T., Fokker, P.A., et al., 2015. Natural-fracture reactivation in shale gas reservoir and resulting microseismicity. SPE J. 54 (6), 450–459. https:// doi.org/10.2118/178437-PA. Sheng, M., Li, G., Huang, Z., et al., 2015. Pore-scale modeling and analysis of surface diffusion effects on shale-gas flow in Kerogen pores. J. Nat. Gas Sci. Eng. 27 (4), 979–985. https://doi.org/10.1016/j.jngse.2015.09.033. Spiegel, M.R., Lipschutz, S., Spellman, D., 2009. Vector Analysis. Schaum’s Outlines, second ed. McGraw Hill, USA. Sun, H., Chawathe, A., Hoteit, H., et al., 2015. Understanding shale gas flow behavior using numerical simulation. SPE J. 20 (1), 142–154. https://doi.org/10.2118/ 167753-PA. Wang, K., Liu, H., Luo, J., et al., 2017. A comprehensive model coupling embedded discrete fractures, multiple interacting continua and geomechanics in shale gas reservoirs with multi-scale fractures. Energy Fuels 31 (8), 7758–7776. https://doi. org/10.1021/acs.energyfuels.7b00394. Wang, M.R., Li, Z.X., 2003. Nonideal gas flow and heat transfer in micro- and nanochannels using the direct simulation Monte Carlo method. Phys. Rev. E. 68 (4), 046704 https://doi.org/10.1103/PhysRevE.68.046704. Wang, Y.H., Shahvali, M., 2016. Discrete fracture modeling using centroidal voronoi grid for simulation of shale gas plays with coupled nonlinear physics. Fuel 163, 65–73. https://doi.org/10.1016/j.fuel.2015.09.038. Wang, Y., Zheng, G., Wood, K., et al., 2016. Microseismic mapping improves understanding of a complex reservoir: a case study in a southern Sichuan shale gas field. In: SPE Asia Pacific Hydraulic Fracturing Conference, 24-26 August, Beijing, China. https://doi.org/10.2118/181872-MS. Warpinski, N.R., Du, J., Zimmer, U., 2012. Measurements of hydraulic-fracture-induced seismicity in gas shales. SPE J. 27 (3), 240–252. https://doi.org/10.2118/151597PA. Warren, J.E., Root, P.J., 1963. The behavior of naturally fractured reservoirs. SPE J. 3 (3), 245–255. https://doi.org/10.2118/426-PA. Watanabe, K., Takahashi, H., 1995. Parametric study of the energy extraction from hot dry rock based on fractal fracture network model. Geothermics 24 (2), 223–236. https://doi.org/10.1016/0375-6505(94)00049-I. Wei, S.M., Xia, Y., Jin, Y., et al., 2019. Study on the 3D fluid-solid coupling model with multi-pressure system of shale (in Chinese). Sci Sin Phys. Mech. Astron. 49, 1–13. Wei, S.M., Xia, Y., Jin, Y., et al., 2019. Quantitative study in shale gas behaviors using a coupled triple-continuum and discrete fracture model. J. Pet. Sci. Eng. 174, 49–69. https://doi.org/10.1016/j.petrol.2018.10.084. Wilson, K.C., Durlofsky, L.J., 2013. Optimization of shale gas field development using direct search techniques and reduced-physics models. J. Pet. Sci. Eng. 108, 304–315. https://doi.org/10.1016/j.petrol.2013.04.019. Williams-Kovacs, J.D., Clarkson, C.R., 2015. A modified approach for modelling 2-phase flowback from multi-fractured horizontal shale gas wells. In: Unconventional Resources Technology Conference. https://doi.org/10.15530/URTEC-20152149183. Wu, K., Li, X., Guo, C., et al., 2016. A unified model for gas transfer in nanopores of shalegas reservoirs: coupling pore diffusion and surface diffusion. SPE J. 21 (5), 1583–1611. https://doi.org/10.2118/2014-1921039-PA. Wu, Y.S., Li, J., Ding, D., et al., 2014. A generalized framework model for the simulation of gas production in unconventional gas reservoirs. SPE J. 19 (05), 845–857. https:// doi.org/10.2118/163609-PA. Xia, Y., Jin, Y., Chen, K.P., et al., 2017. Simulation on gas transport in shale: the coupling of free and adsorbed gas. J. Nat. Gas Sci. Eng. 41, 112–124. https://doi.org/ 10.1016/j.jngse.2017.02.024. Xia, Y., Jin, Y., Chen, M., Chen, K.P., 2019. An enriched approach for modeling multiscale discrete-fracture/matrix interaction for unconventional-reservoir simulations. SPE J. 24, 1349–1374. https://doi.org/10.2118/194012-PA. Xu, Y., Guan, Z., Xu, C., et al., 2019. Numerical method and analysis of ultrasonic detection of gas kick in deepwater risers during Offshore drilling. Int. J. Heat Mass Transf. 136, 1311–1326. https://doi.org/10.1016/j.ijheatmasstransfer.2019.03.102. Yang, R., Huang, Z., Li, G., et al., 2017. A semianalytical approach to model two-phase flowback of shale-gas wells with complex-fracture-network geometries. SPE J. 22 (6), 1808–1833. https://doi.org/10.2118/181766-PA. Zhang, N., Yao, J., Huang, Z., Wang, Y., 2013. Accurate multiscale finite element method for numerical simulation of two-phase flow in fractured media using discrete-

21

S. Wei et al.

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

fracture model. J. Comput. Phys. 242 (6), 420–438. https://doi.org/10.1016/j. jcp.2012.12.006. Zhang, Q.F., Huang, Z.Q., Yao, J., et al., 2017. Two-phase numerical simulation of discrete fracture model based on multiscale mixed finite element method (in Chinese). Chin. Sci. Bull. 62, 1392–1401.

Zhang, M., Yao, J., Sun, H., et al., 2015. Triple-continuum modeling of shale gas reservoirs considering the effect of kerogen. J. Nat. Gas Sci. Eng. 24, 252–263. https://doi.org/10.1016/j.jngse.2015.03.032.

22