Copyright © the author(s). This work is licensed under a Creative Commons Attribution 4.0 International License. DOI: 10.14800/IOGR.1185 Received February 23, 2021; revised March 15, 2021; accepted April 6, 2021. *Corresponding author: dongzhenzhen1120@hotmail.com 1 Optimization of Flowback Operations for Shale Gas Wells Cheng Dai, Research Institute of Petroleum Exploration And Development, SINOPEC Group, Beijing, China; Zhenzhen Dong*, Xi’an Shiyou University, Xi’an, China; Xiang Li, Ennosoft Ltd, Beijing, China; Weidong Tian, Xi’an Shiyou University, Xi’an, China Abstract Hydraulic fracturing has become one of the most important aspects of well completion because of the ever- increasing development of unconventional oil and gas reservoirs. Uncontrolled fluid flowback after the fracture treatment has the potential to impact well productivity and profitability, as aggressive flowback strategies can damage a well’s completion while conservative flowback procedures may hinder near-term economic performance. The goal of this study was to develop an integrated flowback model of fracturing fluid that can accelerate production while minimizing the risk of a decrease in the estimated ultimate recovery or a damage to the productivity of the well. First, by considering the complex situations occurring during fracturing fluid flowback, a fluid-solid coupling numerical model corresponding to fracture closure was established and solved with the finite-element method. A proppant transport model was established to investigate the time and displacement of proppant migration within the fracture. Fluid velocity in the fracture created as a result of flowback was considered, along with its effects on proppant movement and localized fracture closure. Finally, by integrating the two-phase tubing flow model, the proppant transport model, and the fluid-solid coupling model, an integrated reservoir-fracture-wellhead model was established, and the principles and methods for designing the flowback scheme were determined. The work considered actual data and information on flowback in shale gas wells from the available literature, as well as our own experiences. The results reveal how choke management, or flowback strategies, can impact potential damage mechanisms in the reservoir, completion, and wellbore. The proposed model can provide optimum flowback design and therefore lead to maximized near-wellbore fracture conductivity and maximum-attainable conductive length in communication with the wellbore. Introduction Multistage fracturing technology for horizontal wells has been applied to shale gas development for many years and has achieved great success within the history of the oil and gas industry. To create hydraulic fractures in horizontal wells, a high volume of fracturing fluid, potentially including slick water and various additives, is injected into the shale formation, along with proppant. The oil industry has developed a suite of tools or models to mitigate proppant flowback, including forced closure, resin-coated sands, deformable proppants, mechanical screens, fracture packs, resin injection, and surface modification agents (Shor and Sharma 2014). However, there are no standard industry practices or recommendations for well flowback to maximize well productivity and minimize proppant flowback. Current practices for flowback operations are, in general, based on rules-of-thumb and are embodied in mailto:dongzhenzhen1120@hotmail.com 2 the confidential flowback operational procedures of various well operators. Such rule-of-thumb flowback practices can cause extensive tensile rock failure, excessive proppant flowback, fines migration, and scale formation (Barree and Mukherjee1995). In tight sandstone gas reservoirs usually produced by single-wing hydraulically induced fractured wells, the forced closure flowback process is used for fracturing fluid flowback within 20 to 30 minutes after the treatment of the fracturing pump (Canon et al. 2003). While for shale gas wells, most multistage fracturing horizontal wells undergo more than 3 days of soaking, coiled tubing is used to drill and grind the bridge plug of the wellbore, and then multistage mixed flowback testing is carried out. Fracturing fluid loss, microfracture communication, and shale immersion in artificial fractures during soak affect reservoir properties near the wellbore and then affect gas well productivity. According to the fracturing technology, string characteristics and formation, the flowback process of shale gas wells usually goes through three stages: closed control, maximum production, and stable production. According to our experience, the whole flowback period lasts around 2 weeks, and flowback rate could be up to 30% to 40%. Thus, unlike tight gas wells, procedure and time for flowback operations are dictated by economic considerations and reservoir properties for shale gas wells. Sometimes, it is desirable to conduct flowback operations immediately after the fracturing treatment so that the well can benefit from nondissipated reservoir pressure. There are also reservoirs where wells show better production performance after “seasoning,” when fracturing fluid is allowed to dissipate in the formation for several weeks before initiating flowback procedures. In all cases after the flowback is initiated, it is desirable to flow the well back at the maximum technically and operationally allowable rate so that the well can be put into production quickly. At the same time, the flowback rate should not exceed certain limits defined by the formation and type of injected materials, as exceeding these limits may result in excessive flowback of propping material, formation destabilization, and, as a result, poorer well production performance. Thus, the optimization of fracturing fluid flowback strategy is mainly embodied in optimizing flowback timing and flowback rate. Early modeling work to understand proppant transport behavior was based on fluid flow velocity, proppant transport, and settling velocity predicted by broad correlations (Bratli and Risnes 1981). Later work (Robinson et al. 1988; Ely et al. 1990; Biezen et al. 1994) used probabilistic methods to calculate two-dimensional (2D) grid block failure of fracture proppant packs. A critical flow rate depended on proppant size, and closure stress was the key cause of catastrophic pack failure. Shortly after, a different approach based on Novotny’s work was introduced that described the movement of individual proppant grains. Rather than directly modeling the effects of fluid, these studies (Asgian et al. 1995; Gidley et al. 1995; Andrew and Kjorholt 1998; Parker et al. 1999; Crafton 2010; Wang et al. 2020) added point-drag forces to grains to calculate the viscous effect of fluid. Now, a fluid-solid coupling numerical model corresponding to fracture closure has been established and solved with the finite-element method (FEM). This paper presents this model, along with a coupled proppant transport model to investigate the stability of proppant packs, and makes recommendations on flowback procedures. Flowback Mechanism Model The flowback model is a system of mathematical equations that characterize the pressure distributions and fluid flow in hydraulically fractured reservoir rock and well piping over time as a function of wellhead pressure or choke size. The mathematical model of reservoir fluid flow showing fracturing fluid flowback is the same as the model showing fracturing fluid injection. The proppant particles suspended in the fracturing fluid enter the hydraulic fracture system during the course of fracturing fluid injection, but they cannot enter the matrix pore given the size scale (nanoscale matrix pore versus millimeter diameter proppant). During flowback, proppant particles in the hydraulic fracture system flow under natural 3 convection, as does the fluid. If reservoir pressure drops, the fluid leaked into the matrix during the fracturing fluid injection period may gradually release into the fracture during the flowback period. The flow of fluid through the fracture system is driven by viscous and gravity forces, simultaneously carrying suspended proppant particles. Proppant particle flow is propelled by both vertical forces, such as gravity and buoyancy, and horizontal forces, such as drag, inertia, and collision. Fluid flow in the matrix is propelled by viscous and gravity forces, and does not carry proppant particles. During the flowback of the fracture fluid, the fluid first flows to the fracture, flows next to the horizontal wellbore along the fracture, and then arrives at the wellhead along the wellbore (horizontally and vertically). Thus, to simulate the whole process of fracturing fluid flowback, one should simulate the flow in the well, as well as through the reservoir and/or along the hydraulic fracture. Figure 1 illustrates an embodiment of the flowback model including submodels. It incorporates a horizontal well model, reservoir model, and fracture flow model. 1) The horizontal well model is a two-phase pressure loss correlation that calculates the fracturing fluid and gas flow rates in the well over time as a function of the wellhead pressure. The horizontal well model receives input (bottomhole pressure and fluid velocity) from a proppant transport model and a fluid-solid coupling flow model. 2) The proppant transport model simulates the movement of particles along the fracture. 3) The fluid-solid coupling flow model incorporates the output of a fluid-fluid displacement model within the fractures, a model of the geomechanical behavior of the formation rock, and a reservoir model that models the inflow of fluid from the reservoir into the fractures.  The geomechanical model models the interaction among the stresses, pressures, and temperatures in the reservoir rock and the hydraulic fracture.  The fluid-fluid displacement model models the displacement of gas by hydraulic fracturing fluid in the reservoir rock and also the displacement of the hydraulic fracturing fluid by the resident reservoir fluids.  The reservoir model models the physical space of the reservoir by an array of discrete cells delineated by an irregular grid. Due to the high conductivity of hydraulic fractures, the flow of the fracturing fluid can be considered Darcy’s flow. The interaction of the fracturing fluid and the proppant flow is coupled by the equivalent viscosity. Figure 1—Schematic representation of the component parts of a flowback model. Mathematical Model and Discretization Assumptions. The following assumptions were made: (1) The flowback model comprises four interconnected domains, including wellhead, horizontal well, Proppant transport model Horizontal well model Fluid-solid coupling flow model Geomechanical model Fluid-fluid displacement model Reservoir model C o u p li n g 4 fractures, and matrix. (2) A set of fine grids with high permeability was used to characterize the primary fracture, which has a half-length (Lf), width (wf), and height (hf). (3) Proppant transport includes the drag force, inertia force, and collision force. (4) The proppant settlement and the change in proppant bed height occur due to hydraulic fracture closure. Horizontal Well Flow Model. Fluid flowing in a wellbore will experience pressure loss. Pressure loss can be categorized as hydrostatic, frictional, or kinetic. For wellbores, kinetic loss is often minimal and can be ignored. Many fluid correlations can be derived empirically that account for hydrostatic and frictional fluid losses in a wellbore under various flow conditions, including modified Beggs and Brill, Petalas and Aziz, Flanigan, and modified Flanigan. In this study, the modified Beggs and Brill method was used to calculate the pressure loss along the wellbore. Fluid-solid Coupling Model. Geomechanical Governing Equation. We treated the porous medium as the superimposition of two continua: skeleton and fluid. The physical model was based on the governing equations of quasistatic poroelasticity (Biot 1941, 1956; Rice and Cleary 1994). We assumed the porous medium to be of isotropy and of infinitesimal transformation and the fluid in the porous medium to be isothermal, single-phase, and compressible. The governing equations for fluid flow and mechanics were taken from mass and linear-momentum balance, respectively. Constitutive relationship for solid rock can be written as (Cook et al. 1974): (𝜎 − 𝜎0) + b(𝑝− 𝑝0)𝑰 = 𝑫휀̃,…………………………………………………………………………….….(1) where σ is the Cauchy total stress tensor; b is the Biot’s coefficient; p is the pore fluid pressure; I is thee unit matrix, I = [1 1 1 0 0 0]𝑇; D is the elastic tensor matrix; and 휀̃ is the linearized strain tensor under the assumption of infinitesimal transformation. The strain tensor, 휀̃ , is related to the displacement vector, �̃�, via the kinematic compatibility relations as follows: 휀̃ = 1 2 (∇�̃� + (∇�̃�)𝑇),……………………………………………………………………………………(2) The governing equation for mechanical deformation of the rock mass can be written as follows. ∇𝛔 + 𝜌𝑏𝒈 = 0,………………………………………………………………………………………….(3) where g is the gravity vector; 𝜌𝑏 = 𝜙𝜌𝑓 + (1 − 𝜙)𝜌𝑠 is the bulk density; 𝜌𝑓 is the fluid density; 𝜌𝑠 is the solid density; and 𝜙 is the true porosity. True porosity is defined as the ratio between current void volume and current bulk volume. Substituting Eqs. 1 and 2 in Eq. 3, the governing equation for rock can be expressed as: ∇ [ 1 2 𝐷(∇�̃� + (∇�̃�)𝑇) + 𝜎0 − b(𝑝 − 𝑝0)𝑰] + 𝜌𝑏𝒈 = 0…………………………….…………………...(4) Fluid Flow Governing Equation. Fluid flow in the fractured porous medium during the flowback period experiences three stages: flow in the porous matrix, flow in the fracture, and cross-flow between the matrix and fracture. The governing equations for the three stages of fluid flow are described in the following paragraphs. The fluid flow constitutive equation can be expressed as: 𝜕(𝜌𝑓𝑉𝑝) 𝑉𝑏𝜕𝑡 + ∇(𝜌𝑓𝒗𝑓) = 𝑞 𝑉𝑏 ,………………………………………………………………………………..(5) where 𝑉𝑝 is the pore volume; 𝑉𝑏is the bulk volume; 𝜌𝑓 is the fluid density; q is a sink/source term for 5 fluid; and 𝒗𝑓is the fluid velocity vector, which can be expressed by Darcy’s law: 𝒗𝑓 = − 𝐾 𝜇 (∇𝐩 − 𝜌𝑓𝒈),………………………………………………………………………………….(6) where 𝜇 is the fluid dynamic viscosity; K is the permeability vector; and p is the pressure. Based on poroelastic theory and Darcy’s law, the governing equation for single-phase fluid flow (Eq. 5) in porous media can be derived as (see derivative details in Attachment 1): 𝑏 𝜕𝜀𝑏 𝜕𝑡 + 1 𝑀 𝜕𝑝 𝜕𝑡 = 𝑞 𝜌𝑓𝑉𝑏 + ∇ [ 𝐾 𝜇 (∇𝐩 − 𝜌𝑓𝒈)],………………………………………………………………(7) where 휀𝑏is the bulk stress tensor; M is Biot’s modulus, 1 𝑀 = ϕ𝑐𝑓 + 𝑏−𝜙 𝐾𝑠 ; 𝑐𝑓 is the fluid compressibility; K is the skeleton modulus; and 𝑏 = 1 − 𝐾𝑑 𝐾𝑠 . Eqs. 4 and 7 constitute the governing equations for the rock and fluid, respectively. After discretizing the equations, the pressure and displacement fields of the model could be determined. Numerical Discretization. Time discretization was applied using a backward first-order and a fully implicit finite-difference scheme. We adopted FEM to discretize the objective domain, both for fluid flow and geomechanics (Nassir 2013). The benefit of adopting the same method is that FEM can easily handle general boundary conditions, complex geometry, and variable material properties. For most heterogeneous grid blocks, permeability is understood to be a 3×3 second-order tensor with stress-dependent elements. Because in FEM, the governing differential equation for flow is integrated over each discretized domain, permeability can be considered in its full tensor form. In FEM, we rewrote Eqs. 4 and 7 in an equivalent weak variation form. We partitioned the domain into nonoverlapping elements, Ω = 𝑈𝑗=1 𝑛𝑒𝑙𝑒𝑚Ω𝑗,…………………………………………………………………………………………..(8) where nelem is the number of elements. A linear space was defined as V={Functions v: v is a continuous function in the domain Ω and has a piecewise continuous and bounded first partial derivatives in Ω, and ν(Γ) = 0} The discrete approximation of the continuum problem of Eqs. 4 and 7 could then be given as: ∫ ∇𝑁𝑢: 𝜎𝑑Ω = ∫ 𝑁𝑢𝜌𝑏𝑔𝑑Ω + ∫ 𝑁𝑢𝑡̅𝑑Γ ΓΩΩ ,…………………………………………………………..(9) ∫ 𝑁𝑝 𝜕 𝜕𝑡 ( 𝑝 𝑀 + 𝑏휀𝑏)𝑑Ω + ∫ 𝑁𝑝∇𝒗𝑓𝑑Ω = ∫ 𝑁𝑝�̅�𝑑Γ. ΓΩΩ ..…………………………………………….(10) When the weight coefficients of solid and fluid were discretized by FEM, the variables of displacement field and pressure field were placed on the element vertex. The form of the coefficient matrix is shown in Figure 2. It can be concluded that the dimension of coefficient matrix for a model with N nodes is 4N×4N, with the coefficient matrix being a sparse matrix. To reduce the storage space, the coefficient matrix was stored in compressed sparse row (CSR) format. The calculation unit used in this study was hexahedron element, as shown in Figure 2. 6 Figure 2—Stiffness matrix form for FEM. Proppant Transport Model. Proppant Transport Grid. In a quasicontinuous medium, the flow passage of proppant does not really exist. To visualize the flow passage of proppant, a flow passage of proppant was made equivalent to the fracture element generated by the above algorithm, and the topological relationship between these flow passages was established, as shown in Figure 3. Figure 3—A schematic illustration of the topological structure of proppant. Proppant-Governing Equation. In the framework of a quasicontinuous medium, the proppant transport process satisfies the mass conservation equation (Boronin and Osiptsov 2014): 𝜕(𝐶𝜙𝐹) 𝜕𝑡 + ∇(𝐶𝜙𝐹�⃗� 𝑝) = 𝑞𝑖𝑛𝑗,……………..…….………………………………………………………(11) where C is the concentration of proppant in the fracturing fluid, cm3/cm3; ∅𝐹 is the fracture porosity (percentage of the fracture volume in the grid); �̅�𝑝 is the proppant migration velocity, cm/s; and 𝑞𝑖𝑛𝑗 is the proppant injection rate, cm3. Proppant migration velocity in a fracture is determined by fracturing fluid flowback velocity (Blyton et al. 2015): 7 �̅�𝑝 = 𝑘𝑟𝑒𝑡�̅�𝑓 ,...…………………………………………………………………………………………(12) where 𝑘𝑟𝑒𝑡 is the retardation factor derived from an experiment of fluid-particle flow considering the drag force, inertia force, and particle-to-particle and particle-to-wall collision forces. 𝑘𝑟𝑒𝑡 = 1 + ( 𝑑𝑝 𝑤𝑒 ) − 2.02 ( 𝑑𝑝 𝑤𝑒 ) 2 ,…………..…………………………………………………………..(13) where 𝑑𝑝 is the diameter of proppant, cm; and we is the effective fracture width, cm, which is related to the proppant diameter and the concentration. 1 𝑤𝑒 2 = 1.411 ( 1 𝑑𝑝 2 − 1 𝑤2 ) 𝐶0.8,………………….……………………………………………………......(14) where w is the fracture width, cm. There is an additional settling velocity of the proppant in the vertical direction: { 𝑢𝑝,𝑥 = 𝑘𝑟𝑒𝑡𝑢𝑓,𝑥 𝑢𝑝,𝑦 = 𝑘𝑟𝑒𝑡𝑢𝑓,𝑦 𝑢𝑝,𝑧 = 𝑘𝑟𝑒𝑡𝑢𝑓,𝑧 + 𝑢𝑠 ,.……….………………………………………………………………………(15) where 𝑢𝑠 is the proppant-settling velocity, cm/s, which is generally used as the corrected velocity based on Stokes’ law: 𝑢𝑠 = 𝑢𝑠𝑡𝑜𝑘𝑒𝑠𝑓(𝑅𝑒𝑝, 𝐶) = 𝑔𝐷𝑝 2(𝜌𝑝−𝜌𝑓) 18𝜂𝑓 𝑓(𝑅𝑒𝑝, 𝐶, 𝑤𝐹),……….……………………………………….(16) where Dp is the proppant diameter, cm; 𝜌𝑝 and 𝜌𝑓 are the particle and fluid densities, g/cm3, respectively; 𝜂𝑓is the hydrodynamic viscosity, mPa.s; and 𝑓(𝑅𝑒𝑝, 𝐶, 𝑤𝐹) is the correction factor, which is related to the particle Reynolds number, concentration, and fracture width. It is commonly believed that the Reynolds number, particle concentration, and fracture width are independent, so 𝑓(𝑅𝑒𝑝, 𝐶, 𝑤𝑓) can be expressed as: 𝑓(𝑅𝑒𝑝, 𝐶, 𝑤𝑓) = 𝑓1(𝑅𝑒𝑝)𝑓2(𝐶)𝑓3(𝑤𝑓),………………………………………………………………(17) 𝑓2(𝐶) = 𝑒−5.9𝑐,……………………………………………………………………………………….(18) 𝑓3(𝑤𝑓) = ∑ 𝑅𝑖 ( 𝐷𝑝 𝑤𝑓 ) 𝑖 𝑛 𝑖=0 ,……………………………………………………………………………..(19) where Ri is a calibrating model parameter. 𝑓1(𝑅𝑒𝑝) has no implicit expression. It is usually calculated iteratively based on the traction model and gravity balance as: 𝑑�⃗⃗� 𝑝 𝑑𝑡 = 𝐷𝑝(�⃗� 𝑓 − �⃗� 𝑝) + (1 + 𝜌𝑓 𝜌𝑝 )𝑔 = 0,………………………………………………………………(20) where: 𝐷𝑝 = 3𝜌𝑓 8𝜌𝑝 𝐶𝑑 |�⃗⃗� 𝑓−�⃗⃗� 𝑝| 𝑟 ,…………………………………………………………………………………..(21) 𝐶𝑑 = { 24 𝑅𝑒𝑝 (1 + 1 6 𝑅𝑒𝑝 2/3) 𝑅𝑒𝑝 < 103, 0.44 103 ≤ 𝑅𝑒𝑝 ≤ 3 × 105. ,……………………………….……………….(22) Thus, by combining Eq. 20 with Eqs. 18 and 19, the particle settlement velocity due to inertial effects can be calculated by Eq. 16. Numerical Discretization. We applied finite volume methodology to discretize proppant-governing 8 equation Eq. 11. Before presenting the discretization of Eq. 11, we present a discussion here. Note that in the quasicontinuous medium method, the real fracture structure is not analyzed. The existence of the fracture in the mesh and the direction of fracture propagation are determined by current stress conditions. Correspondingly, when calculating proppant migration, the macroscopic migration behavior should be considered with the grid as the basic control body, and the proppant flow between the grids depends on the fluid flow in the fracture. No matter how complex the fracture network is, the change in proppant concentration can be calculated by calculating the interposition equation of the channeling flow in the fracture between different grids. In the framework of a quasicontinuous medium, the total flow between grids can be obtained by Darcy's law, and the flow through the fracture can be separated by the average permeability law. If using the explicit central second-order scheme to discretize the flow term in the proppant-governing equation, the discrete form can be written as follows (Li et al. 2016): (𝐶𝜙𝐹)𝐶 𝑛+1−(𝐶𝜙𝐹)𝐶 𝑛 Δ𝑡 Δ𝑉 + ∑ [𝑤𝑐(𝐶𝜙𝐹)𝐶 𝑛 + 𝑤𝑛(𝐶𝜙𝐹)𝑁 𝑛 ](�⃗� 𝑝,𝐼 𝑛+1�⃗� Δ𝑆)𝑠 = 𝑞𝑖𝑛𝑗,𝐶Δ𝑉,………....……………..(23) where subscriptions C and N are the current grid and neighboring grid, respectively; �⃗� 𝑝,𝐼 𝑛+1 is the particle transport velocity at the interface in time step n+1, which is precalculated; and 𝑤𝑐 and 𝑤𝑛 are the weights of the first-order approximation at the interface between the current grid and the adjacent grid. The explicit discretization mentioned above is obviously unstable (referring to the explicit second-order central difference scheme). In order to ensure accuracy, some variables must be treated implicitly, and a set of linear equations must be constructed to solve them. The first-order upwind scheme is absolutely stable, but its accuracy is low. Combining the advantages and disadvantages of the two schemes, we implicitly dealt with the convective upwind scheme and modified the right hand of equation (RHE) by using the second-order margin. Thus, Eq. 23 can be discretized as: (𝐶𝜙𝐹)𝐶 𝑛+1−(𝐶𝜙𝐹)𝐶 𝑛 Δ𝑡 Δ𝑉 + ∑ [𝑤𝑐(𝐶𝜙𝐹)𝐶 𝑛+1 + 𝑤𝑛(𝐶𝜙𝐹)𝑁 𝑛+1](�⃗� 𝑝,𝐼 𝑛+1�⃗� Δ𝑆)𝑠 = ∑ [𝑤𝐷(𝐶𝜙𝐹)𝑢 𝑛 + 𝑤𝑛(𝐶𝜙𝐹)𝑁 𝑛 ](�⃗� 𝑝,𝐼 𝑛+1�⃗� Δ𝑆)𝑠 + 𝑞𝑖𝑛𝑗,𝐶Δ𝑉……………………(24) The relationship among the upwind grid, dead-wind grid, current grid, and adjacent grid should be determined according to the flow velocity. The above analysis is based on the suspending region. If the case of sand heap is considered, the first- order dead scheme must be adopted; otherwise, the result will be nonphysical. For this reason, when the concentration of the upwind grid in the direction of particle settlement approaches the limit accumulation concentration, the interface between the two grids must be corrected to adopt the first-order dead scheme. Critical Flow Rate of Proppant Flowback. The maximum stable pressure drop gradient for the proppant flowback within the fracture after closure can be obtained with the sum of the pressure drop gradients that the proppant adhesion force and fracture closure stress can resist. ( 𝑑𝑝 𝑑𝑥 ) 𝑚𝑎𝑥 = ( 𝑑𝑝 𝑑𝑥 ) 𝑠𝑡𝑎 +( 𝑑𝑝 𝑑𝑥 ) 𝐹𝑉 ,……………...…………………………………………………………..(25) where ( 𝑑𝑝 𝑑𝑥 ) 𝑠𝑡𝑎 is the pressure drop gradient that the fracture closure stress can resist (MPa/m); and ( 𝑑𝑝 𝑑𝑥 ) 𝐹𝑉 is the pressure drop gradient that the proppant adhesion force can resist (MPa/m), ( 𝑑𝑝 𝑑𝑥 ) 𝐹𝑉 = 3𝐹𝑛 106×4𝜋𝑅3.……………………………………………………………………………………(26) The effective cohesion force can be expressed as a function of cohesion, gravity, and buoyance. 9 𝐹n = 𝑎𝑐 ( 𝛾0 𝛾0∗ ) 2.5 𝜋 2 𝜌𝑓휀𝑑 5 3 √3 − (𝜌𝑠 − 𝜌𝑓) 1 6 𝜋𝑑3g,……………………………………….……………(27) where 𝜌𝑠 is the proppant density, kg/m3; 𝜌f is the fluid density, kg/m3; and 휀 is a constant, 1.75 cm3/s2. The pressure that the fracture closure stress can resist is a function of fracture width and fracture closure stress, ( 𝑑𝑝 𝑑𝑥 ) 𝑠𝑡𝑎 = 𝑊𝑇𝑒𝑥𝑝 [−0.5 ( 𝑙𝑛𝑝𝑛𝑒𝑡−𝑎′ 𝑆𝑇 ) 2 ],………………………………………………………………(28) 𝑊𝑇 = 1422.5exp(−1.0483𝑊𝑟),…………………………………………………………………….(29) 𝑊𝑟 = 𝑊𝑓 𝑑 ,……………………………………………………………………………………………..(30) where Wf is the fracture width, cm; d is the proppant diameter, cm; pnet is the fracture net closure pressure, Pa; 𝑎′ is a constant, 7.7172; and ST is a function of apparent strength of proppant. 𝑆𝑇 = 3 × 10−5S𝑀𝑎𝑥 + 0.22368,……………………………………………………………………..(31) where Smax is the apparent proppant strength. According to the stress and deformation of the element, the fracture width is calculated by the constitutive model of the fracture. In this study, we adopted Barton et al.’s (1985) constitutive model to calculate fracture width. Barton et al. (1985) proposed an empirical correlation between mechanical aperture, am, and hydraulic aperture, wf, on the basis of experimental data: 𝑤𝑓 = { 𝑎𝑚 2 𝐽𝑅𝐶−2.5 𝛿ℎ ≤ 0.75𝛿𝑝𝑒𝑎𝑘 √𝑎𝑚𝐽𝑅𝐶𝑚𝑜𝑏 𝛿ℎ > 𝛿𝑝𝑒𝑎𝑘 , ,……….……………………………………………….....(32) where JRC is the joint roughness coefficient; JRCmob is the mobilized value of JRC; 𝛿 is the shear displacement in the horizontal direction; h and v are horizontal and vertical, respectively; and 𝑤𝑓is the fracture width of the smooth fracture wall, which is the fracture width used in this study. Note that a linear interpolation determines the hydraulic aperture value when 0.75 < 𝛿ℎ/𝛿𝑝𝑒𝑎𝑘 < 1.0, and both am and wf are in the unit of mm. Given the fracture stress state, the mechanical aperture can be calculated as (Li et al. 2017): 𝑎𝑚 = 𝑎0 − ∆𝑎𝑛 + 𝛿𝑣 ……………………………………………………………………………….....(33) Fracture closure stress refers to the force on the fracture wall required to close the fracture after shut-in. Fracture closure stress is an important factor that affects fracture conductivity and can be calculated from the instantaneous shut-in pressure at the wellhead, the fluid column pressure, and the reservoir pressure. In addition, the minimal horizontal stress of the formation affects the fracture closure stress. On the basis of the study of triaxial stress and mechanical parameters of rock, an equation for closure stress was derived as 𝑝𝑛𝑒𝑡 = 𝛾 1−𝛾 S𝑣+𝑆ℎ𝑖+𝐴𝑝𝑒× 𝑝𝑖 2 1− 𝐴𝑝𝑒 2 ………………………...……………………………………………………...(34) The closure stresses in different regions can be obtained in field tests. If conditions permit, the best methods for obtaining closure stress are stepped injection, flowback, and equilibrium tests prior to hydraulic 10 fracturing. Otherwise, the closure stress can be estimated with Eq. 33. If the stability of the proppant filling layer is characterized with the critical flow rate, the Darcy equation can be introduced. 𝑁𝑅𝑒 < 10 𝑑𝑝 𝑑𝑥 = 𝜇 𝐾 𝑣,…………………………………………………………………….(35) 𝑁𝑅𝑒 > 10 𝑑𝑝 𝑑𝑥 = 𝜇 𝐾 𝑣 + 𝛽𝜌𝑣2,………..………………………………………………………(36) where μ is the fluid viscosity (Pa·s); v is the seepage velocity (m/s); ρ is the fluid density (kg/m3); p is the pressure (Pa); β is the seepage velocity (m–1); and K is the permeability (m2). Eqs. 35 and 36 were transformed to obtain the critical fluid velocity for zero proppant production after fracture closure. 𝑁𝑅𝑒 < 10 𝑣 = 𝐾 𝜇 ( 𝑑𝑝 𝑑𝑥 ) 𝑚𝑎𝑥 ,……………………………..…………………………………(37) 𝑁𝑅𝑒 > 10 𝑣 = − 𝜇 𝐾 +√( 𝜇 𝐾 ) 2 +4𝜌𝛽( 𝑑𝑝 𝑑𝑥 ) 𝑚𝑎𝑥 2𝜌𝛽 ……..…………………………………………………….(38) According to Eqs. 37 and 38, the fluid within the closed fracture exhibits mainly a gas-liquid flow for shale gas wells. Thus, calculating the critical fluid velocity for the zero proppant production after the fracture closure requires the maximum stable pressure drop gradient of the proppant-filling layer, mixed fluid viscosity, and mixed fluid density within the fracture. Full Model Coupling. By coupling the proppant transport model, fluid-solid coupling model, and gas- liquid two-phase tubing flow models, the optimal flowback velocity for the stable proppant-filling layer and minimal fluid loss could be calculated with the optimization algorithm; in addition, the choke size for the fluid flowback was adjusted by observing the wellhead pressure. Explicit iteration was used to couple the fracture model and the proppant migration model (Figure 4). Firstly, FEM was used to simulate the current time-step evolution process of fracture width, and then the proppant migration was simulated to obtain the concentration distribution. Then, the proppant distribution was used to update the viscosity of sand-carrying fluid and the fracture-equivalent conductivity, and it was substituted into the model of the average permeability tensor of a quasicontinuous medium to achieve coupling. 11 Figure 4—Schematic of the fully coupled model. Optimization of Fracturing Fluid Flowback in a Shale Gas Well As an example, the flowback scheme of Well A in the Fuling shale gas field (a horizontal multistage fractured gas well) was established with the reservoir and fracturing treatment parameters. Fracturing fluid flowback was simulated using the proposed integrated model for a hydraulically fractured horizontal well with a 1,000-m horizontal wellbore and 10 fracture stages located in the center of a shale reservoir. Each stage created transverse primary fractures along the horizontal wellbore with a fracture spacing of 25 m. All primary hydraulic fractures were assumed to be identical and to penetrate the whole reservoir. A specific fracturing fluid flowback scheme was planned for a horizontal shale gas well in Fuling. The reservoir and fracturing treatment parameters of the gas well are listed in Table 1. 12 Table 1—Reservoir and fracturing treatment parameters of a shale gas well (Well A) in Fuling. Parameter type Parameter Unit Value Reservoir parameter Effective reservoir thickness m 12 Poisson's ratio Dimensionless 0.23 Young's modulus MPa 28,000 Fracture closure pressure MPa 52 Reservoir pressure MPa 64.2 Shut-in pressure at wellhead MPa 69 Reservoir permeability 10–3μm2 0.003 Filter-cake permeability 10–3μm2 1.6 Relative permeability 10–3μm2 0.02 Reservoir porosity Dimensionless 0.09 Comprehensive compressibility MPa–1 0.001 Fracturing treatment parameter Proppant density kg/m3 2,650 Proppant particle size mm 0.8 Fracturing fluid density kg/m3 1020 Fracturing fluid viscosity mPa·s 2 Fracture half-length m 200 Fracture height m 40 Fracture dynamic maximum width mm 6 Adhesion force coefficient dyn/cm 2.6 Filtrate viscosity mPa·s 1.6 Borehole vertical height m 2300 Borehole radius m 0.05 Horizontal section length m 500 Stage number Dimensionless 5 Nozzle outlet pressure MPa 0.1 Local resistance coefficient Dimensionless 0.5 Treatment displacement m3/min 3 Total leak-off coefficient m/min0.5 0.00029 Absolute roughness of the borehole mm 2 Model Validation. To validate the proposed fracturing flowback model, we matched the measured casing pressure with the predicted pressure from the injection period to the flowback period, as shown in Figure 5. They were in good agreement. Consequently, the fracture width was predicted and is illustrated in Figure 6. The fracture width clearly declined as pressure dropped during the flowback period. 13 Figure 5—History match of bottomhole pressure for Well A during fracture treatment. Figure 6--Predicted fracture width for single fracture of Well A during fracture treatment. In the first couple of days, Well A flowback showed a relatively high flow rate and a choke size of 8 mm, which caused sand production (orange dots in Figure 7a and 7b). The proposed model was run to calculate critical flow rate and corresponding critical choke size, which is the maximum flow rate that will not cause proppant flowback (blue curves in Figure 7 and 7b, respectively). As actual flow rate is greater than critical flow rate and actual choke size is larger than the estimated choke size, proppant flowed back with fracturing fluid. After controlling flow rate by reducing choke size to be smaller than critical choke size, there was no sand production anymore. Thus, the proposed model was proved to be robust enough to optimize the flow rate and choke size. 0 20 40 60 80 0 50 100 150 200 P re ss u re , M P a t,min Field Data Simulated Result 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0 50 100 150 200 F ra ct u re w id th , m t, min 14 (a) (b) Figure 7—Flowback performance versus (a) flow rate and (b) choke size. Proppant Transport. In addition to flow pressure and flow rate, the proposed model can also display the progress of proppant migration from injection period to flowback period. Figure 8 shows proppant concentration profiles perpendicular to the fracture. The x-axis is the position along the fracture, which refers to the center of the fracture. The proppant concentration during the fracturing period ranged from 1.5% to 2.5%. From the concentration change curve in Figure 8, an area of high proppant concentration was found to form at the fracture tip due to the severe filtration. After stopping pumping (“End of fracturing” curve in Figure 8), the proppant concentration in the fracture front decreased gradually, but fracturing fluid in the fracture continued to advance because of high pressure in the fracture. At this time, the velocity of fluid in the fracture decreased and the distribution of fluid in the fracture gradually became uniform. However, it should be noted that the fluid continued to move forward. Figure 9 is the 2D distribution of proppant concentration in a single fracture in the course of fracturing to flowback. 0 50000 100000 150000 200000 250000 1 7 13 19 25 31 F lo w b ak r at e, m 3 /d Flowback time, days Critical flow rate Field flow rate Flow with sand 0 1 2 3 4 5 6 7 8 9 1 7 13 19 25 31 C h o k e si ze , m m Flowback time, days Critical choke size Field choke size 15 Figure 8—Proppant concentration profiles perpendicular to the fracture in the course of injection and flowback. (a) 48 min (b) 72 min (c) 113 min (d) End of fracturing (e) Start of flowback (f) Flowback for 20 min Figure 9—2D proppant distribution in the fracture in the course of injection and flowback. Fracturing Fluid Flow. The model was run to investigate the factors that affect fracture closure time. The pressure profiles of the wellhead during the fracturing fluid flowback period were established with the basic data in Table 1 by changing the matrix permeability (Figure 10) and the choke size (Figure 11). When the wellhead pressure was lower than the fracture closing pressure (53 MPa in this case), the fracture was considered closed. It was obvious that the fracture closure time gradually increased with decreasing matrix permeability and choke size. The lower the matrix permeability, the slower the fracture pressure drop and the longer the fracture closure time. And with decreased permeability, the decreasing trend of fracture 16 pressure was slower. In other words, for shale gas reservoirs, it takes longer for a fracture to close because of the extremely low permeability. But by increasing choke size in the wellhead, it is possible to obtain an optimal flow rate that can reduce fracture closing time and prevent proppant flowback. Figure 10—Wellhead pressure in the fracture with various matrix permeabilities. Figure 11—Wellhead pressure in the fracture with various choke sizes. Optimization of Flowback Scheme.On the basis of the above analysis, we applied the proposed model to optimize the flowback mode and choke size corresponding to measured wellhead pressure. Traditional Flowback Mode. The traditional flowback mode refers to the choke size (radius of 2 mm) not vary during fracturing fluid flowback. The wellhead pressure profile obtained by the proposed model is shown in Figure 12. The result shows that the fracture closed at 1,440 minutes and the pressure decreased by 32.4 MPa in 30 days. 20 30 40 50 60 70 80 0 10 20 30 40 50 W el lh ea d p re ss u re , M P a Time, Days k=0.001mD k=0.003mD k=0.005mD k=0.01mD Fracture closure pressure (52 MPa) 20 30 40 50 60 70 80 0 10 20 30 40 50 W el lh ea d p re ss u re , M P a Time, Days r=2mm r=3mm r=4mm r=5mm r=6mm Fracture closure pressure (52 MPa) 17 Figure 12—Wellhead pressure curve in the traditional flowback scheme. Natural Closure Mode. In natural closure mode, the well is shut in until the fracture is closed after fracturing, and then the flowback initiates. The predicted wellhead pressure is presented in Figure 13. The fracture closed at 1,584 minutes, and then choke size increased from 4 to 8 mm at 7.6 days. The pressure decreased by 34.5 MPa in 30 days. It was observed that when the fracture closed and the choke size increased, the pressure decreased suddenly. Figure 13—Wellhead pressure in natural fracture closure mode. Forced Closure Mode. In forced closure mode, flowback is performed immediately with a choke size of 2 mm after fracturing. The estimated wellhead pressure is presented in Figure 14. The fracture closed and the choke size increased to 4 mm at 1,044 minutes. After 7.1 days, the choke size increased to 8 mm. The pressure dropped by 35.8 MPa in 30 days. 0 10 20 30 40 50 60 70 80 0 5 10 15 20 25 30 35 W el lh ea d p re ss u re , M P a Time, Days Fracture closure pressure (52 MPa) 0 10 20 30 40 50 60 70 80 0 5 10 15 20 25 30 35 W el lh ea d p re ss u re , M P a Time, days No choke Choke size=4mm Choke size=8mm 18 Figure 14—Wellhead pressure curve in forced closure mode. The comparison of the above three schemes shows that in the forced closure mode of the fracture, the pressure decreased more rapidly, the fracture closed earlier, and the final stress equilibrium time was shorter, causing less damage to the reservoir. Optimization of Choke Size. According to the above analysis, the optimal choke size under forced closure mode without proppant production corresponding to wellhead pressure is presented in Figure 15. The fracture was closed when the wellhead pressure was 53 MPa. Before fracture closure, the critical flowback rate was 2.5 m3/h and the safe choke size was 2 mm; after fracture closure, the critical flowback rate was 4 to 14.6 m3/h, and the choke size without proppant production was 4.2 to 7.8 mm. One can refer to Figure 15 to adjust the choke size based on the measured wellhead pressure during the fracturing fluid flowback period. Figure 15—Maximum choke size and flow rate without proppant production for different wellhead pressures. Conclusions An integrated reservoir-fracture-wellbore model was successfully applied to evaluate the stability of the proppant pack and optimize fracturing fluid flowback procedures based on wellhead pressure in multistage fractured horizontal shale gas wells. The model accuracies were verified with an example, and the following conclusions were drawn. 0 10 20 30 40 50 60 70 80 0 5 10 15 20 25 30 35 W el lh ea d p re ss u re , M P a Time, Days Choke Size=2mm Choke Size=4mm Choke Size=8mm 0 4 8 12 16 20 0 2 4 6 8 10 20 30 40 50 60 70 80 F lo w r at e, m 3 /h C h o k e si ze , m m Wellhead pressure, MPa Choke Size Critical Flow Rate 19 1. Through simulations, the effects of the different parameters on wellhead pressure, fracture closure time, and final equilibrium time during fracturing fluid flowback were determined: larger choke size and higher reservoir permeability result in faster wellhead pressure drop, shorter fracture closure time, and shorter equilibrium time (Figures 10 and 11). 2. During fracturing fluid flowback, several parameters influence the critical flow rate and choke size required for proppant pack stability: as the fracture is closed, the critical flowback rate and choke size increase rapidly because of the stress acting on the proppant. Higher proppant density allows a higher critical flow rate and larger choke size before the fracture closure, but the choke size required for stable proppant pack after the fracture closure does not vary. The larger particle size of the proppant results in a higher critical flow rate and larger choke size for a safe blow-off before fracture closure. The opposite is the case after the fracture closure. The higher viscosity of the fracturing fluid corresponds to a lower critical flow rate and a smaller choke size for safe blow-off. 3. The proposed model for fracturing fluid flowback in a shale gas well considers the reservoir, fracture, and wellbore. Thus, the model provides the theoretical basis for the optimization of postfracture flowback time, fracture closure mode, and choke management. On the basis of the model, the flowback process of a horizontal shale gas well with multistage fracturing in Fuling was optimized. According to the results, the choke size should be increased from 2 to 4 and 8 mm in forced closure mode (Figure 15). Acknowledgements We would like to thank the sponsors of the project entitled “Study on Fracturing Fluid Flowback Model and Flowback Scheme for Shale Gas Well” (G5800-18-ZS-KFGY002) from State Key Laboratory of Shale Oil and Gas Enrichment Mechanisms and Effective Development, SINOPEC Group, for financial support. We also would like to thank Prof. Zhenzhen Dong from Xi’an Shiyou University and Dr. Xiang Li from Ennosoft Ltd. for technical support. Conflicts of Interest The author(s) declare that they have no conflicting interests. Reference Andrew, J.S. and Kjorholt, H. 1998. Rock Mechanical Principles Help to Predict Proppant Flowback from Hydraulic Fractures. Paper presented at SPE/ISRM Rock Mechanics in Petroleum Engineering, Trondheim, Norway, 8-10 July. SPE-47382-MS. Asgian, M.I., Cundall, P. A., and Brady, B.H.G. 1995. The Mechanical Stability of Propped Hydraulic Fractures: A Numerical Study. Journal of Petroleum Technology 4(3):75-78. SPE-28510-PA. Asadollahi, P. and Tonon, F. 2010. Constitutive Model for Rock Fractures: Revisiting Barton's Empirical Model. Engineering Geology 113(4):11-32. Barton, N., Bandis, S., and Bakhtar, K. 1985. Strength, Deformation and Conductivity Coupling of Rock Joints. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts 22(3):121-140. Barree, R.D. and Mukherjee, H. 1995. Engineering Criteria for Fracture Flowback Procedures. Paper presented at Low Permeability Reservoirs Symposium, Denver, Colorado. 19-22 March. SPE-29600-MS. Biezen, E.N., Bruining, J., and Molenarr, J. 1994. An Integrated Model for Underground Coal Gasification. Paper presented at the SPE Annular Conference and Exhibition. New Orleans, 25-28 September. SPE-28583-MS. Bratli, R.K. and Risnes, R. 1981. Stability and Failure of Sand Arches. Society of Petroleum Engineers Journal 21(2): 65-72. SPE-8427-PA. 20 Boronin, S.A. and Osiptsov, A.A. 2014. Effects of Particle Migration on Suspension Flow in a Hydraulic Fracture. Fluid Dynamics 49(2): 208-221. Biot, M.A. 1941. General Theory of Thress Dimensional Consolidation. J Appl. Phys. 8: 155-164. Biot, M. A. 1956. General Solution of the Equation of Elasticity and Consolidation for a Porous Material. J. Appl. Mech. 78:91-96. Blyton, C.A., Gala, D.P., and Sharma, M.M. 2015. A Comprehensive Study of Proppant Transport in a Hydraulic Fracture. Papre presented at the SPE Annual Technical Conference and Exhibition. Houston, Texas, USA, 28-30 September. SPE-174973-MS. Cook, R.D., Malkus, D.S., Plesha, M.E., et. al. 1974. Concepts and Applications of Finite Element Analysis (2nd ed.). New York: Wiley. Canon, J.M., Romero, D.J., Pham, T.T., et. al. 2003. Avoiding Proppant Flowback in Tight-Gas Completions with Improved Fracture Design. Paper presented at SPE Annual Technical Conference and Exhibition, Denver, Colorado, 5-8 October. SPE-84310-MS. Crafton, J.W. 2010. Flowback Performance in Intensely Naturally Fractured Shale Gas Reservoirs. Paper presented at SPE Unconventional Gas Conference, Pittsburgh, USA, 23-25 February. SPE-131785-MS. Ely, J.W., Arnold, W.T., and Holditch, S.A. 1990. New Techniques and Quality Control Find Success in Enhancing Productivity and Minimizing Proppant Flowback. Paper presented at SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, 23-26 September. SPE-20708-MS. Gidley, J.L., Penny, G.S., and McDaniel, R.R. 1995. Effect of Proppant Failure and Fines Migration on Conductivity of Propped Fractures. SPE Production & Facilities 10(1): 74-78. SPE-24008-PA. Parker, M., Weaver, J., and Van Batenburg, D. 1999. Understanding Proppant Flowback. Paper presented at SPE Annual Technical Conference and Exhibition, Houston, Texas, 3-6 October. SPE-56726-MS. Jaeger, J.C. and Cook, N. 1979. Fundamentals of Rock Mechanics (4th ed. Vol. volume 9). London: Chapman & Hall. Robinson, B.M., Holditch, S.A., and Whitehead, W.S. 1988. Minimizing Damage to a Propped Fracture by Controlled Flowback Procedures. Journal of Petroleum Technology 40(6): 753-759. Li, S., Li, X., and Zhang, D. 2016. A Fully Coupled Thermo-Hydro-Mechanical, Three-Dimensional Model for Hydraulic Stimulation Treatments. Journal of Natural Gas Science and Engineering 34(4): 64-84. Nassir, M. 2013. Geomechanical Coupled Modeling of Shear Fracturing in Non-Conventional Reservoirs. PhD dissertation. University of Calgary, Calgary, Canada. Shor, R. J. and Sharma, M.M. 2014. Reducing Proppant Flowback From Fracture: Factors Affecting the Maximum Flowback Rate. Paper presented at SPE Hydraulic Fracture Technology Conference, Woodlands, Texas, 4-6 February. SPE-168649-MS. Wang, F., Chen, Q., Liu, X., et al. 2020. Fracturing-Fluid Flowback Simulation with Consideration of Proppant Transport in Hydraulically Fractured Shale Wells. ACS Omega 5: 9491-9502. 21 Appendix A 𝑑휀𝑏 = 𝑑𝑉𝑏 𝑉𝑏 = 𝑐𝑏𝑐𝑑𝜎 + 𝑐𝑏𝑝𝑑𝑝,………………………………………………………………………..(A1) 𝑑휀𝑝 = 𝑑𝑉𝑝 𝑉𝑝 = 𝑐𝑝𝑐𝑑𝜎 + 𝑐𝑝𝑝𝑑𝑝,………………………………………………………………………..(A2) 𝑑𝜙 𝜙 = 𝑑( 𝑉𝑝 𝑉𝑏 ) 𝑉𝑝 𝑉𝑏 = 𝑑𝑉𝑝 𝑉𝑏 − 𝑉𝑝 𝑉𝑏 2𝑑𝑉𝑏 𝑉𝑝 𝑉𝑏 = 𝑑휀𝑏 − 𝑑휀𝑝,……………….……………………………………………..(A3) where 𝑏 = 1 − 𝐾𝑑 𝐾𝑠 ,…………………………………………………………………………………………..(A4) 𝑐𝑏𝑐 = 𝑏 𝐾𝑑 + 1 𝐾𝑠 = 1− 𝐾𝑑 𝐾𝑠 𝐾𝑑 + 1 𝐾𝑠 = 𝐾𝑠−𝐾𝑑 𝐾𝑠𝐾𝑑 + 1 𝐾𝑠 = 1 𝐾𝑑 ,………………………………..……………………..(A5) 𝑐𝑏𝑝 = 𝑏 𝐾𝑑 ,…..…………………………………………………………………………………………(A6) 𝑐𝑝𝑐 = 𝑏 𝜙𝐾𝑑 ,……………………………………………………………………………………………(A7) 𝑐𝑝𝑝 = 𝑏 𝜙𝐾𝑑 − 1 𝐾𝑠 ,………………………………………………………………………………………(A8) where 휀𝑏 is the bulk stress tensor; 휀𝑝 is the pore stress tensor; 𝜙 is true porosity; Ks and Kd are the skeleton modulus. Fluid compressibility is given as 𝑐𝑓 = 1 𝜌𝑓 𝑑𝜌𝑓 𝑑𝑝 ,……………………………………………………………………………………………(A9) 𝜕(𝜌𝑓𝑉𝑝) 𝑉𝑏𝜕𝑡 = 𝜌𝑓 𝜙 𝑉𝑝 𝜕𝑉𝑝 𝜕𝑡 + 𝜙 𝜕𝜌𝑓 𝜕𝑡 = 𝜌𝑓 𝜙 𝑉𝑝 𝜕𝑉𝑝 𝜕𝑡 + 𝜙𝜌𝑓 1 𝜌𝑓 𝜕𝜌𝑓 𝜕𝑝 𝜕𝑝 𝜕𝑡 .…….…...………..….…………………….(A10) Substituting Eqs. A2 and A9 into Eq. A10 𝜕(𝜌𝑓𝑉𝑝) 𝑉𝑏𝜕𝑡 = 𝜙𝜌𝑓 ( 𝜕𝜀𝑝 𝜕𝑡 + 𝑐𝑓 𝜕𝑝 𝜕𝑡 )…..………………………………..…………..……………………….(A11) From Eq. A2 𝜕𝜀𝑝 𝜕𝑡 + 𝑐𝑓 𝜕𝑝 𝜕𝑡 = 𝑐𝑝𝑐𝜕𝜎+𝑐𝑝𝑝𝜕𝑝 𝜕𝑡 + 𝑐𝑓 𝜕𝑝 𝜕𝑡 = 𝑐𝑝𝑐 𝜕𝜎 𝜕𝑡 + (𝑐𝑓 + 𝑐𝑝𝑝) 𝜕𝑝 𝜕𝑡 ……….....…………………………….(A12) Substituting Eqs. A11 and A12 into Eq. A10 𝜕(𝜌𝑓𝑉𝑝) 𝑉𝑏𝜕𝑡 = 𝜙𝜌𝑓 (𝑐𝑝𝑐 𝜕𝜎 𝜕𝑡 + (𝑐𝑓 + 𝑐𝑝𝑝) 𝜕𝑝 𝜕𝑡 ),………………………...…………………………………(A13) Substituting Eq. A13 into the fluid constitutive equation (Eq. 5) and dividing 𝜌𝑓 on both sides of Eq. 5, Eq. 5 can be simplified as 𝜙𝑐𝑝𝑐 𝜕𝜎 𝜕𝑡 + 𝜙(𝑐𝑓 + 𝑐𝑝𝑝) 𝜕𝑝 𝜕𝑡 + ∇(𝜌𝑓𝒗𝑓) 𝜌𝑓 = 𝑞 𝜌𝑓𝑉𝑏 .………………....……………………………………..(A14) From Eq. A1 𝜕𝜎 𝜕𝑡 = 𝜕𝜀𝑏−𝑐𝑏𝑝𝜕𝑝 𝑐𝑏𝑐𝜕𝑡 ……………………………………………….………………………………………(A15) Thus, 22 𝜙𝑐𝑝𝑐 𝜕𝜎 𝜕𝑡 = 𝑏 𝐾𝑑 𝜕𝜎 𝜕𝑡 = 𝑏 𝐾𝑑 𝜕𝜀𝑏− 𝑏 𝐾𝑑 𝜕𝑝 1 𝐾𝑑 𝜕𝑡 = 𝑏 𝜕𝜀𝑏 𝜕𝑡 − 𝑏2 𝐾𝑑 𝜕𝑝 𝜕𝑡 ……………......……………………………………(A16) From Eq. A8 𝜙(𝑐𝑓 + 𝑐𝑝𝑝) 𝜕𝑝 𝜕𝑡 = 𝜙 (𝑐𝑓 + 𝑏 𝜙𝐾𝑑 − 1 𝐾𝑠 ) 𝜕𝑝 𝜕𝑡 ,…..………………….…………………………………….(A17) ∇(𝜌𝑓𝒗𝑓) 𝜌𝑓 = ∇𝒗𝑓 + 𝒗𝑓 ∇(𝜌𝑓) 𝜌𝑓 = ∇𝒗𝑓 + 𝒗𝑓𝑐𝑓∇𝑝…………………………………………………………(A18) Substituting Eq. A16 through A18 into Eq. A14 𝑏 𝜕𝜀𝑏 𝜕𝑡 + [𝜙 (𝑐𝑓 + 𝑏 𝜙𝐾𝑑 − 1 𝐾𝑠 ) − 𝑏2 𝐾𝑑 ] 𝜕𝑝 𝜕𝑡 + ∇𝒗𝑓 + 𝒗𝑓𝑐𝑓∇𝑝 = 𝑞 𝜌𝑓𝑉𝑏 ………....………………………….(A19) In Eq. A19 𝜙 (𝑐𝑓 + 𝑏 𝜙𝐾𝑑 − 1 𝐾𝑠 ) − 𝑏2 𝐾𝑑 = 𝜙𝑐𝑓 + 𝑏(1−𝑏) 𝐾𝑑 − 𝜙 𝐾𝑠 = 𝜙𝑐𝑓 + 𝑏 𝐾𝑑 𝐾𝑠 𝐾𝑑 − 𝜙 𝐾𝑠 = 𝜙𝑐𝑓 + 𝑏−𝜙 𝐾𝑠 ...…..………………(A20) Let 1 𝑀 = 𝜙𝑐𝑓 + 𝑏−𝜙 𝐾𝑠 ….…………………………………………………………………………………...(A21) Eq. A19 can be finally simplified as 𝑏 𝜕𝜀𝑏 𝜕𝑡 + 1 𝑀 𝜕𝑝 𝜕𝑡 + ∇𝒗𝑓 + 𝒗𝑓𝑐𝑓∇𝑝 = 𝑞 𝜌𝑓𝑉𝑏 ………...………………………..……………………………(A22) Based on Darcy’s law 𝒗𝑓 = − 𝐾 𝜇 (∇𝐩 − 𝜌𝑓𝒈)………………………………………………………………………………(A23) Eq. A22 can be expressed as 𝑏 𝜕𝜀𝑏 𝜕𝑡 + 1 𝑀 𝜕𝑝 𝜕𝑡 = 𝑞 𝜌𝑓𝑉𝑏 + ∇ [ 𝐾 𝜇 (∇𝐩 − 𝜌𝑓𝒈)] + 𝐾 𝜇 (∇𝐩 − 𝜌𝑓𝒈)𝑐𝑓∇𝑝,...……………..………………….(A24) 𝐾 𝜇 (∇𝐩 − 𝜌𝑓𝒈)𝑐𝑓∇𝑝 is negligible as it is a second derivative. Thus, the governing equation of fluid flow can be expressed as 𝑏 𝜕𝜀𝑏 𝜕𝑡 + 1 𝑀 𝜕𝑝 𝜕𝑡 = 𝑞 𝜌𝑓𝑉𝑏 + ∇ [ 𝐾 𝜇 (∇𝐩 − 𝜌𝑓𝒈)]…………...……………………..………………………...(A25) Cheng Dai is a senior engineer in petroleum exploration and development research institute, SINOPEC. Dr. Dai holds a Bachelor's degree and a Ph.D. degree from Peking University. His research interest includes numerical simulation and shale gas development. Zhenzhen Dong is a professor in the Department of Petroleum Engineering at Xi’an Shiyou University. She worked as a reservoir engineer at Schlumberger from 2012 to 2016. Dong holds a bachelor’s degree in mathematics from Northeast Petroleum University, China; a master’s degree in petroleum engineering from the Research Institute of Petroleum Exploration and Development, China; and a Ph.D. degree in petroleum engineering from Texas A&M University. Xiang Li is the technical director of Beijing Energy Innovation Software Company. His work includes scientific computing, fracture modeling, and reservoir simulator development. Dr. Li holds a bachelor’s degree and a PhD degree from Peking Univerisy. 23 Weidong Tian is a master's candidate in the Department of Petroleum Engineering at Xi’an Shiyou University. He has focused his research in areas involving reservoir simulation, well testing, and production analysis. Tian holds a BS degree in Petroleum Engineering from Xi’an Shiyou University.