Complex Fracture Closure Pressure Analysis During Shut-in: A Numerical Study

Conventional fracture closure models typically assume homogeneous formation. This assumption renders complex fracture closure pressure analysis during shut-in insufficient at present. In this paper, we use a cohesive zone method (CZM)-based finite element model to obtain the fracture closure pressure and minimum horizontal principal stress of the major fracture and the branch fracture based on pressure fall-off analysis. Key factors including leak-off rate, injection time, and stress anisotropy are discussed in detail. A quantitative relationship between fracture closure pressure and these factors is plotted to reduce or eliminate the errors caused by conventional fracture injection diagnostic models. The results show that leak-off rate, injection time, and stress anisotropy have a significant effect on the fracture closure process. Fracture closure in naturally fractured formations is a slow process, and the early closure pressure represents the closure of the branch fracture, which is much higher than the minimum horizontal stress. This investigation provides new insight into how to estimate the in-situ stress in naturally fractured reservoirs.


Introduction
Hydraulic fracturing is a well stimulation technology that is effective in increasing hydrocarbon production in unconventional reservoirs such as shale gas, coal-bed methane, and tight sandstone (Montgomery and Smith, 2010;Dahi Taleghani et al., 2020;Wang et al., 2020b). A large volume of viscous fluid is continuously injected into such a formation to promote hydraulic fracture propagation, which occurs at a depth of several thousand metres and accordingly cannot be directly observed. The hydraulic fracturing process is very complex (Nolte, 1986;van Dam et al., 2000;Smith and Montgomery, 2015). However, the pressure data recorded at the time of an injection via a diagnostic fracture injection test (DFIT) can provide a wealth of information on fracture geometry, fluid leak-off characteristics, and fracture closure time (Nolte, 1986(Nolte, , 1991Mayerhofer et al., 1995;Plahn et al., 1997). Among DFITs, fracture pressure decline analysis is a common method of determining the fracture closure pressure and closure time due to the elimination of wellbore friction after shut-in (de Pater et al., 1996;Kholy et al., 2019). Thus, fracture closure analysis has been one of the trending topics in research on hydraulic fracturing over the past decades.
A laboratory experiment was designed to investigate hydraulic fracture closure behaviour in different materials, including cement, plaster, soft plaster, and diatomite (van Dam et al., 2000). The results showed that the radius changes of penny-type fractures after shut-in significantly influence fluid leak-off characteristics, mechanical behaviour, and fracture closure pressure. The closure time increases as the fractured radius decreases in size, while the calculated leak-off coefficient is overestimated when the fractured radius increases in size.
Over the last few decades, different fracture closure models have been developed to determine key parameters such as closure pressure, formation permeability, leak-off coefficient, and fluid efficiency. These models include the square root plot, log-log plot, G-function plot, and empirical equation (Nolte, 1986(Nolte, , 1991Barree et al., 2007;Belyadi et al., 2019;Kholy et al., 2019). The square root plot suggests that closure pressure occurs at the inflection point on the second derivative of pressure (Belyadi et al., 2019). The straight line between the origin and the inflection point deviates from the second derivative curve. The log-log plot is employed to determine fracture closure and various flow regimes such as the linear flow regime, bilinear flow regime, and pseudoradial flow regime (Mayerhofer et al., 1995;Barree et al., 2007). Closure occurs at the point where the slope of the second derivative curve changes from positive to negative. G-function derivative analysis requires a plot that includes the bottom-hole pressure, the pressure derivative (dP/dG), and the superimposed pressure derivative (GdP/dG). The characteristic curved shape of the pressure derivative and the superposition pressure derivative is used to identify the leak-off behaviour, of which there are four common types, i.e, normal leak-off behaviour, pressure-dependent leak-off from opened natural fractures, fracture-height recession during shut-in, and fracture-tip extension (McLellan and Janz, 1991;Barree and Mukherjee, 1996;Craig et al., 2000;Chipperfield, 2006;Cai and Dahi Taleghani, 2021). Plahn et al. (1997) proposed a numerical model of the fracture pump-in and flow-back test, which gives a better estimate of fracture closure pressure than the aforementioned conventional methods. Based on the signal processing approach, Eltaleb et al. (2021) proposed a novel methodology for DFIT analysis in geothermal reservoirs, which is independent of rock mechnical properties and fracture geometry. Wang et al. (2021b) anaylzed the characteristic of pressure decline in fractured horizontal wells, and build up a new model which can quantitatively estimate closure stress and average pore pressure. This model can also infer the uniformity of proppant distribution in the hydraulic fractures. Carpenter (2021) pointed out that a slow increase in the dP/dG plot indicates that long-time fracture closure is caused by early near-wellbore fracture closure. However, it is not practical to use the aforementioned fall-off analysis methods to determine fracture closure pressure in tight formations because the closure time is very long (Barree et al., 2014). Kholy et al. (2019) developed an empirical equation to predict hydro-fracture closure pressure based on the instantaneous shut-in pressure (ISIP) using subsurface solids-injection data. The newly proposed method can be used to obtain the fracture closure pressure with a relative error of less than 3% compared to the closure pressure measured in the oilfield. Noting that a hydraulic fracture retains a limited width after asperities come into contact, McClure et al. (2016) proposed a fracture compliance method that picks up the fracture closure pressure from DFITs. Afterwards, the proppant transport with gravitational settling and fracture closure are integrated into a threedimensional (3D) hydraulic fracture simulator (Shiozawa and McClure, 2016). A new constitutive equation is introduced into the 3D simulator and can seamlessly describe the fracture closure process during shut-in. This model includes the nonlinear joint closure law denoting fracture compliance and fracture surface roughness, as well as the effect of proppant accumulation in the hydrofractures. He et al. (2017) considered the different gas-rate production values of fractures and developed a semi-analytical method to diagnose the locations of underperforming hydro-fractures through pressure-transient analysis in tight gas reservoirs. This method was derived from the Green function and was verified by comparison using the commercial software KAPPA Ecrin. Using post-closure pulse testing, Upchurch (2003) developed a method to determine the fracture closure pressure in soft formations. This method was likewise verified with several field examples. No reduction in the accuracy of the leak-off coefficient when compared to the recorded pressure data was found. Cai and Dahi Taleghani (2021) numerically investigated the inititation and closure behaviour of an axial fracture during DFIT, and the results show that the closure time of the axial fracture from fracture compliance based DFIT analysis is earlier than that of the transverse fracture. They recommend that the combined method of the G × dP/dG plot and the dP/dG plot is helpful for DFIT interpretation.
Taken together, the current fracture closure models do not consider the long closure process of complex fracture networks including branch fractures. Furthermore, the effect of geomechanical behaviour on fracture closure is not fully integrated into the calculations of the fluid flow and leak-off in the current fracture closure models. Which fracture closes first after shut-in: the main hydraulic fracture or the branching hydraulic fracture? Errors may occur when field technicians employ a conventional closure model to identify the closure pressure or minimum horizontal principal stress (σ h ) in tight formations. According to Barree et al. (2014), such poor or nonexistent results stem from either problematic data acquisition or the incorrect analysis of the acquired data.
The goal of this paper is to numerically investigate complex fracture closure pressure analysis during shut-in. We employ a cohesive zone method (CZM)-based finite element model to identify the fracture closure pressure and minimum horizontal principal stress of major and branch fractures from pressure fall-off analysis. Key factors including leak-off rate, injection time, and stress anisotropy are discussed in detail. A quantitative relationship between fracture closure pressure and these factors is plotted to reduce or eliminate the errors caused by conventional fracture injection diagnostic models. This investigation provides new insight into how to estimate the in-situ stress in naturally fractured reservoirs.
The paper is organized as follows: The "Method and Theory" section presents the governing equations of hydraulic fracturing problems, including the stress equilibrium equation for rock deformation and the lubrication equation for fluid flow in hydraulic fractures. The "Numerical Simulation" section discusses the effects of key factors including leak-off coefficient, injection time, and stress anisotropy on the fracture closure pressure and in-situ stress estimation. The quantitative relationship between fracture closure pressure and these factors is also presented. The main conclusions of this study are summarized in the last section.

Method and theory
We consider the hydraulic fracture closure process to be a two-dimensional (2D) coupled stressseepage field problem in porous media. This problem consists of a rock matrix that includes predefined zero-thickness cohesive elements representing the hydraulic fracture propagation paths. CZM is utilized to simulate the damage initiation and evolution of hydraulic fracturing problems (Barenblatt, 1962;Wang et al., 2018b;Wang et al., 2020a).
The coupled problem consists of rock deformation and fluid flow in the rock matrix and fracture networks. Fluid exchange exists between the porous rock matrix and fractures. The pore space of the rock matrix is considered fully saturated. We assume that only a single-phase flow exists in the rock matrix, as this type of flow usually dominates the early flow-back process. We assume the fluid flow in the rock matrix obeys Darcy's flow. According to the theory of linear poroelasticity (Biot, 1941), the governing equation for the rock matrix is expressed as where σ ij are the geo-stress components; σ 0 ij are the initial geo-stress components; E is the elastic modulus of rock; ν is Poisson's ratio; α is the Boit constant; ϵ ij are the strain components; p is the pore pressure; and p 0 is the initial pore pressure. It is worth noting that we use initial stress to obtain a more efficient convergence of our numerical model.
According to the principle of mass balance, the fluid flow in fractures is expressed as where w f is fracture width; t is injection time; x is fracture location; v n is the leak-off velocity normal to fracture surface (Carter, 1957); and v t is the fluid velocity along the fracture. The nominal strains are expressed as where δ n and δ s are the normal and shear separation, respectively; ϵ n and ϵ s are the nominal strain components in the normal and shear directions, respectively; and T 0 is the initial thickness of the cohesive elements. According to traction-separation law (TSL), the constitutive equation of the cohesive elements is expressed as where t n and t s are the two components of nominal stress traction vector t in the normal and shear directions, respectively; and E nn , E ns, and E ss are the components of elastic constitutive matrix E.
The maximum nominal stress criterion is utilized to represent damage initiation. The criterion is expressed as where t n 0 and t s 0 are the peak values of the nominal stress traction vector in the normal and shear directions, respectively. The bilinear TSL, which is suitable for sedimentary rocks in practical applications (Dahi Taleghani et al., 2013), is used to model the fracture zone process in this study. The traction linearly increases with the increase in separation before damage initiation. The cohesive elements are damaged once the traction reaches its peak value. After that, the traction linearly decreases with the increase in separation. When the traction reaches zero, complete damage occurs in the cohesive elements.
As shown in Figure 1, a 2D plane strain model is used to represent the coupled stress-seepage field problem in ABAQUS. The size of the square domain is 200 m × 200 m. The symmetric finite element model is used to reduce the computation cost. The injection point is located at the centre of the left edge. Two groups of cohesive elements (i.e. a hydraulic fracture and a natural fracture, denoted by the red lines in Figure 1) are inserted into the unstructured meshes to represent the possible fracture propagation and closure paths during flow-back. The approach angle between the hydraulic fracture and the natural fracture is denoted as β. Although we focus on flow-back in this study, we include the injection stage in our numerical simulations to provide a realistic fracture pattern for an initial condition of the flow-back model.
We conduct a single planar fracture propagation as the initial simulation to provide an accurate initial condition during flow-back for verification and validation ends. The input parameters for the flow-back simulation model are listed in Table 1, which are from the MaHu tight oil reservoir in Xinjiang Autonomous Region, China. In this planar fracture model, the minimum horizontal stress σ h , the maximum horizontal stress σ H, and the overburden stress σ v are equal to 77 MPa, 82 MPa, and 124 MPa, respectively, and pore pressure is equal to 54 MPa. The rock matrix is considered to be a linear elastic medium in which E = 25 GPa and ν = 0.23. The CZM parameters of E nn and K n are set to 25 GPa and 87.5 GPa, respectively. A linear elastic TSL is used to represent the damage evolution during hydraulic fracturing. Complete cohesive element failure occurs when the fracture energy satisfies the Benzeggagh-Kenane fracture criterion (Benzeggagh and Kenane, 1996). The calculated initial fracture width and the net pressure within the hydro-fracture before the simulation of flow-back presents a non-uniform distribution. Therefore, the analytical step of modelling fracture propagation in ABAQUS plays a very important role in simulating fracture closure during flow-back.

Numerical simulation
In this section, we use a CZM-based finite element model to simulate the major and branch fracture closure processes. Key factors including approach angle, injection time, leak-off rate, and horizontal stress difference are analysed in detail. As shown in Figure 1, the major fracture is a part of the hydraulic fracture in the horizontal direction, and the branch fracture is a part of the natural fracture in the inclined direction.

Model verification
To verify the CZM-based finite element model, the numerical solutions are compared to the analytical solutions. The analytical solutions are derived from the famous KGD hydraulic fracturing model (Zheltov, 1955;Geertsma and De Klerk, 1969). The fracture length L frac (t) and fracture width at the injection point w 0 h (t) are respectively expressed as in Eqs. (7) and (8).
w h 0 (t) = 2.36 where t is the injection rate; q 0 is the injection rate; μ f is the fluid viscosity; E' = 2G/(1-ν) is A rectangular 100 m × 120 m domain is built up in the finite element model, which is divided into 21,280 quadrilateral elements with pore pressure. The injection point is located on the midpoint of the left edge. The injection rate is equal to 2 L/s per metre thickness of the reservoir. The shear modulus of rocks is 8.3 GPa, and the Poisson's ratio is 0.2. The fluid viscosity is 0.001 Pa·s, and the bulk modulus of fluid is 2.2 GPa. A total of 133 cohesive elements with pore pressure is inserted along the central axis of the domain in the horizontal direction. The initial fracture opening is equal to 0.002 mm as an initiation path during the numerical simulation. The result comparison between numerical solutions and analytical solutions is shown in Figure 2. We observe that the numerical solutions match well with the analytical solutions. It shows that the hydraulic fracturing finite element model is comparatively accurate.
To evaluate the model performance, the relative error curves of fracture length and width to their analytical solutions are plotted in Figure 3. We observe that the value of the relative error is less than 0.5%, which indicates that our finite element model is reliable.

Effect of injection time on fracture closure
The contour plots of the maximum principal stress at different injection times are shown in Figure 4. Figure 4(a) corresponds to the time of instantaneous shut-in (t ISIP ), and Figure 4(b) corresponds to the time of branch fracture closure (t branch ). Thus, Figure 4(b) provides the initial conditions for flow-back analysis. Next, we use the instantaneous shut-in time t ISIP as the starting time during the flow-back process. The time of major fracture closure is denoted as t major . Figure 5 presents an elapsed time comparison of the major fracture closure and the branch fracture closure at different approach angles, including 30 degrees, 60 degrees, and 75 degrees. The figure shows that the longer the injection time is, the higher the shut-in pressure is. This is because more fluids are filtrated to the rock matrix normal to the fracture surface, which in turn increases the local pore pressure in the vicinity of the hydraulic fracture. The variation in Figure 5(b) is different from that in Figure 5(b) and 5(c) because the hydraulic fracture patterns are different. Only the upper part of the natural fracture in Figure 5(b) is opened during the interaction between hydraulic fractures and natural fracture. We also observe that there is an inflection point on the curve. The inflection point reflects the fact that the branch fracture is closed during the early stage and the major fracture is closed during the later stage. Thus, fracturing fluid flows from the branch fracture to the major fracture during the flow-back process. Figure 6 presents an elapsed time comparison of the major fracture closure and the branch fracture closure at different approach angles. We observe that the longer the injection time is, the higher the ISIP and the closure pressure of the branch fracture are. We also observe that the closure times of the major fracture and the branch fracture, i.e., t branch -t ISIP and t major -t branch , increase as the injection time increases. Thus, the closure effect of the hydraulic fracture is a slow, gradual process. The closure pressure in the early stage represents the closure of the branch fracture, while the closure  pressure in the late stage represents the closure of the major fracture. In addition, the value of t major--t branch (red curve) in Figure 4(c) is higher than that of t branch -t ISIP , which is different from Figure 4(a) and 4(b). When the intersection angle between hydraulic fracture and natural fracture is high, the opened natural fracture is easy to close during shut-in process under the action of maximum horizontal principal stress.

Effect of leak-off rate on fracture closure
Leak-off rate is an important factor during hydraulic fracturing, which impacts the fluid pressure distribution in hydraulic fractures. Figure 7 presents an elapsed time comparison of the major fracture closure and the branch fracture closure at different approach angles, including 30 degrees, 60 degrees, and 75 degrees. The figure shows that the smaller the leak-off rate is, the higher the shut-in pressure is. This is because less fluid is filtrated to the rock matrix normal to the fracture surface, which in turn increases the fluid pressure within the hydraulic fracture. When the value of the leak-off rate is small, the closure time of the branch fracture is obviously delayed and the closure of the major fracture occurs even later.
As shown in Figure 8, an elapsed time comparison of the major fracture closure and the branch fracture closure at different approach angles is plotted. We observe that the leak-off rate has a noticeable effect on the closure time of the branch fracture or the ISIP. The smaller the leak-off rate is, the longer the closure times of the major fracture and the branch fracture are. The smaller the leak-off rate is, the larger the difference between the closure time of the branch fracture and the major fracture is. When the leak-off rate is high, there is little difference between the closure time of the major fracture and that of the branch fracture.
Given the slow leak-off rate and long fracture closure time in tight formations, the fracture closure pressure obtained on site mainly reflects the closure of the branch fracture. Accordingly, the fracture closure pressure cannot be used to determine the minimum horizontal stress. Moreover, the closure pressure of the branch fracture reflects the complexity of the fracture. The larger the difference between the closure pressure of the branch fracture and the minimum stress, the more complex the fracture is considered to be.

Effect of horizontal stress difference on fracture closure
Horizontal stress difference impacts the interaction between hydraulic fracture and natural fractures, and thus effect of horizontal stress difference on fracture closure is investigated in this section. Figure 9 presents an elapsed time comparison of the major fracture closure and the branch fracture closure at different approach angles, including 30 degrees, 60 degrees, and 75 degrees. The figure shows that the larger the horizontal stress difference is, the lower the shut-in pressure is. We also observe that the larger the horizontal in-situ stress difference is, the longer the closure time of the branch fracture is. Figure 10 presents an elapsed time comparison of the major fracture closure and the branch fracture closure at different approach angles. We observe that as the horizontal principal stress difference increases, the ISIP decreases and the branch fracture closure pressure increases. We also     observe that the larger the horizontal stress difference is, the shorter the closure time of the branch fracture and the longer the closure time of the major fracture are. In the case of a 75-degree approach angle (Figure 8(c)), there is a certain critical stress difference (3.5 MPa). To the left of the critical point, the closure time of the branch fracture is longer than the closure time of the major fracture; to the right of the critical point, the closure time of the major fracture is longer than the closure time of the branch fracture.

Discussion
Natural fractures widely exist in unconventional reservoirs such as shale gas, shale oil and tight sandstone reservoirs. Natural fractures mainly include cemented and partially uncemented natural fractures. Our investigation aims at the formation with fully cemented natural fractures. Our previous investigation show that the reservoir with an uncemented or partially cemented natural fracture have stronger formatiom capcity of complex fracture than a cemented natural fracture (Wang et al., 2018a;Wang et al., 2021a). Thus, the pressure decline response during hydraulic fracturing of cemented natural fracture is different from that of uncemented or partially cemented natural fractures. Futher investigation of the effect of uncemented or partially cemented natural fractures on fracture closure can be done in the next step.

Conclusion
This paper numerically investigated complex fracture closure pressure during shut-in. The CZM-based finite element model was used to determine the fracture closure pressure and minimum horizontal principal stress of the major hydro-fracture and the branch fracture based on pressure fall-off analysis. The effects of key factors including leak-off rate, injection time, and stress anisotropy on fracture closure process were discussed in detail. The main conclusions are as follows: The closure times of the major fracture and the branch fracture increase as the injection time increases. Fracture closure is a slow process, and the early closure pressure represents the closure of the branch fracture, which is much higher than the minimum horizontal stress.
The smaller the leak-off rate is, the longer the closure times of the major fracture and the branch fracture are. When the leak-off rate is high, there is little difference between the closure time of the major fracture and that of the branch fracture.
Due to slow reservoir filtration and long fracture closure times, the fracture closure pressure obtained on site mainly reflects the closure of the branch fracture. Accordingly, it cannot be used to determine the minimum horizontal stress. Moreover, the closure pressure of the branch fracture reflects the complexity of the fracture.
The larger the horizontal stress difference is, the shorter the closure time of the branch fracture and the longer the closure time of the major fracture are.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.