New Pump-In Flowback Model Verification with In-Situ Strain Measurements and Numerical Simulation

<jats:p>This study presents an analytical model for estimating minimum horizontal stress in hydraulic fracturing stimulations. The conventional Diagnostic Fracture Injection Test (DFIT) is not practical in ultra-tight formations, leading to the need for pump-in/flowback tests. However, ambiguities in the results of these tests have limited their application. The proposed model is based on the linear diffusivity equation and material balance, which is analytically solved and verified using a commercially available numerical simulator. The model generates a linear graph in which the pressure drop and its derivative are plotted versus the developed solution time function. The closure pressure is determined when the slope of the derivative deviates from linearity. The model was applied to several cycles of field flowback tests and found to eliminate the ambiguity associated with identifying the fracture closure. Furthermore, the minimum In-situ stresses estimated using this approach are verified via downhole strain measurement and synthetic data from a fully 3D commercial fracturing simulator. The proposed technique outperformed other conventional methods in analyzing challenging injection/shut-in tests, showing improved results and reducing uncertainty in estimated fracture parameters. This model is expected to scale down the need for multiple field trials and provide a reliable estimation of minimum stress.</jats:p>

section_content: In-situ stress describes the local forces acting on lithologic layers in subsurface structures to be exploited. Stresses acting on a formation can be decomposed into two principal stresses, vertical and horizontal. Vertical stress is caused by the overburden weight of the upper geological layers, which usually makes it the largest stress acting on the rock in deep formations. On the other hand, horizontal compressive stresses result from the poroelastic deformation of the rock and external tectonic forces. The parameters that affect the magnitude of In-situ stresses include formation porosity, natural fractures, fluid pore pressure, overburden weight, and mechanical rock properties (such as Young's modulus and Poisson's ratio). The estimation of In-situ stresses is of paramount importance in many fields of subsurface engineering and earth geoscience problems, such as earthquake evaluation [1], enhanced geothermal reservoirs [2, 3], caprock effectiveness in CO 2 storage systems [4] [5] [6] [7] [8] [9] [10], drilling optimization and wellbore-stability [11]. In hydrocarbon reservoirs, the magnitude of In-situ stress is a critical parameter that plays a vital role in hydraulic fracturing treatment design [12] [13] [14] [15], fracture geometry design [16] [17] [18], and proppant placement design applications [19, 20].
section_name: Introduction

section_content: The extended leak-off test is a widely used technique to estimate In-situ stress essential for various engineering aspects of drilling operations. Generally, when a new section of the wellbore is drilled, the drilling fluid is pumped at a constant rate to pressurize the wellbore and fracture the formation. Pressure and rate data are observed during the test through surface or downhole gauges. Figure 1 illustrates a typical pressure response to injected fluids in XLOT. As fluid is pumped into a fixed volume (wellbore), the pressure increases proportionally to the pumping rate until reaching a noticeable deflection on the pressure curve at the leakoff-point (LOP). The inflection point is referred to as fracture initiation pressure (FIP). LOP is generally accepted as an approximation of fracture closure pressure (FCP, which represents minimum horizontal stress). However, it is not recommended in tests run with cased hole completions because of the near-wellbore tortuosity that dominates the pressure response [21, 22]. Alternatively, XLOT is recommended to estimate FCP because it provides far-field stress information as fracture growth increases with continued pumping beyond breakdown pressure (FBP). In permeable formations, the fracture may close as fluids leak-off during the shut-in period. However, in tight formations, the leak-off rate is meager, and a subsequent flowback period is desirable to estimate minimum horizontal stress within a reasonable time period [23].
section_name: A Brief Background of Flowback in Drilling Operations

section_content: Analogous to XLOT, Diagnostic Fracture Injection Test (DFIT) is also used to estimate minimum horizontal stress by injecting one to two barrels of low-viscosity fracturing fluid (without proppant) into the formation to create small fractures ranging from 10 to 20 feet deep into the formation. After the pumps shut down, the pressure is monitored and recorded for hours or days (depending on rock properties) using downhole or surface pressure gauges. One of the limitations of DFIT that limits its applicability in some cases is that the test run time may extend to several days to estimate the closure pressure. In addition, in ultra-tight formations, the leakoff is extremely low; hence, the fracture will require longer to close. Various analytical and numerical models have been developed to analyze fracture injection tests. One of the earliest attempts to estimate In-situ stress from hydraulic fracturing was found in [24, 25]. During that time, Instantaneous shut-in pressure (ISIP) was accepted as the upper approximation for In-situ stress if a small volume of fluid was injected to create a relatively short fracture. An alternative approach was deemed to calculate the closure pressure of larger fractures by plotting falloff pressure versus the square root of time. The most distinguished contributor to the falloff analysis, Nolte in 1979, introduced the material balance-based approach to obtain fracture closure, commonly known as G-function [26]. Two years later, Nolte and Smith 1981, highlighted the importance of the evolution of net pressure (the difference between the pressure in the fracture and closure pressure) in fracture propagation and proppant placement. Their contribution demonstrated the significance of accurately estimating closure pressure before the primary fracturing treatment. In the same work, Nolte and Smith used flowback analysis to determine fracture closure independent of the mentioned earlier decline analysis [27]. The test procedures employ the same tubing/annuals configuration and injection at a sufficient rate to initiate the fracture, then flowback at roughly one-quarter of the injection rate. The pressure plotted against flowback time on a Cartesian plot shows an increase in the decline rate followed by a characteristic reversal of curvature. The closure pressure by Nolte and Smith's approach identifies the pressure at point (A) in Figure 2, preceding the downward curvature. In 1988, Shlyapobersky suggested a repeated sequence of reopening, propagation, shut-in, and the flowback procedure is depicted in Figure 3. The test procedure is similar to the conventional mini-frac and differs only in the fluid volume/rate injected scale. The injection rate is sufficiently maintained during the first cycle to ensure formation breakdown. The test can then be followed by a period of shut-in or constant-rate flowback. The purpose of the flowback is to accelerate the fracture closure and make the dynamic fracture closure more pronounced. The distinct pressure shape during the flowback period is primarily affected by fracture storage/stiffness ( ∂V f / ∂p ) and frictional pressure (at lower flowback rates ∆P friction is minimal). Shlyapobersky in 1988 used the lower inflection point to indicate the mechanical fracture closure shown by point (B) in Figure 3. Soliman and Daneshy, in 1991, proposed a technique based on a compressibility equation coupled with a mass balance equation to calculate closure pressure during the flowback test. In addition, they presented a graphical representation of the system volume (fracture and wellbore volumes), allowing a better understanding of the closure mechanism with a simple yet robust physics-based approach. For instance, from Equation (1), a plot of total flowback volume against corresponding pressure difference would exhibit a straight line with a slope of (C f luid V) as illustrated in Figure 4 [30]. The early part indicates a system volume change, i. e., fracture starts to close, while the straight line represents a constant volume behavior, which can be accepted as the lower bound of closure pressure. As indicated in Figure 5, fluid inside the fracture will stabilize and move to dry areas at the fracture tip when the injection stops. Consequently, the fluid redistribution may provide excess energy that forces the fracture to propagate after shut-in. This combined effect will cause a rapid initial decline in pressure. Following the rapid initial decline, one would expect the decline rate to increase with flowback and fluid leakoff; however, the pressure drop is somehow compensated with fluid expansion and reduction in fracture volume (start of closure) corresponding to point (A). Beyond point (A), the pressure decline accelerates after the fracture tips start to close. Finally, after the fracture is completely closed at point (C), the decline rate drops proportionally to the flowback rate and the formation's ability to produce fluids. Plahn et al. in 1995 proposed another approach to estimate closure pressure based on a numerical model used to simulate the pump-in/flowback test [31]. The approach disagrees with the techniques mentioned above, as the closure pressure from their simulations lies somewhere between the upper and the lower bounds of the Soliman and Daneshy approach. Nevertheless, Plahn et al. concluded that using the intersection of double stiffness lines (fracture stiffness and wellbore stiffness) is the closest approximation of closure pressure obtained from the simulation model, as illustrated by point (C) in Figure 6. However, observing the straight fracture stiffness line may be ambiguous, and it varies from one to another depending on formation and fluid properties as well as flowback rate. Figure 7 shows the summary sketch of techniques mentioned above to calculate closure pressure from the injection flowback test. However, there is no consensus regarding which point represents the definite closure pressure on the stiffness plot. Furthermore, attaining straight-line stiffness and a characteristic curvature requires an appropriate flowback rate which may lead to multiple field trials before obtaining the desired curvature. In this work, we propose an analytical model that not only provides an accurate calculation of closure pressure but may also give an appropriate estimate for pore pressure which is a crucial parameter in hydraulic fracture design and evaluation. Our model is based on a well-testing fundamental equation to model fluid leakoff deployed with suitable boundary conditions to calculate closure pressure without requiring a stiffness straight-line plot.
section_name: A Brief Background of Flowback Tests in Fracture Diagnostics

section_content: In most fracture injection modeling approaches, it is assumed that the injection rate into one wing is constant up to the end of the injection time. After shutting down the pumps, the pressure at the wellbore declines as fluid leaks off from the fracture to the formation. Because the fracture closure process is controlled by leakoff, pressure falloff analysis has been a primary source of obtaining the fracture parameter and reservoir permeability based on the adopted leakoff model for the analysis method. Fundamentally, the models are different in their assumption of leakoff behavior. There are two leakoff approaches adopted by hydraulic fracturing scholars. The first approach was introduced dated back to [33]. This approach fundamentally considers the fluid loss to be a property of a fluid-rock interaction system in which fluid loss velocity u L is given by Carter's equation:

where C L is the leakoff coefficient, and t is the time elapsed since the start of the leakoff. Unfortunately, this approach neglects both the effect of fluid viscosity and solid mechanics, rendering it inadequate to describe the leakoff behavior, especially if Howard and Fast's experiment injection process is not strictly met. Equation (2) can be integrated over time, and it describes the flow of volume V L through surface area A L :

Carter introduced simple mass balance that considers that the fracture has a constant width both in space and in time given by the following equation:

Numerous efforts have been made to combine Carter's model, which assumes a constant fracture width to obtain fracture length as a function of time, and Perkins and Kern's model to estimate fracture width. Nordgren, in 1972, fortified attempts by adding leakoff and evolving fracture volume (resulting from width increase) to the fracture model. The continuity equation used by Perkins, Kern, and Nordgren was based on the mass balance Equation (5) [34] :

where q i (x, t) is the volume flow rate through a cross-section of fracture,q L (x, t) is the volume rate of fluid loss to formation per unit length of fracture, and A c is the crosssectional area of the fracture. Nolte in 1979 and 1986 used Equation (5) to develop the well-known G-function with considering the fluid loss term slightly different than in this paper, which is given by the Equation ( 6 ) [26, 35] :

where u L is given by Equation (1). 

The second approach considers that fluid loss is controlled by fluid flow in porous media. One advantage of this approach is that it depicts the leakoff behavior effectively since it is based on rigorous fluid flow equations. However, the solutions may become complicated for evolving fracture geometry as time progresses. Mayerhofer and Economides, in 1993, presented a model that overcomes the limitation of Carter's "bulk" leakoff model. In their technique, they describe the fracturing fluid flow into a formation with the presence of filter-cake at the fracture face. Although the fluid filtrate extends a few centimeters into the formation, the resulting pressure disturbance extends further into the reservoir. The fluids leakoff linearly and perpendicular to the fracture face governed by linear flow partial differential Equation [36]. The total pressure drop in their model is the summation of three pressure drop components: pressure drop across the fracture face dominated by filter-cake, pressure drop across the invaded zone, and pressure drop across the reservoir. The transient pressure drop is given by:

The advantage of this approach is that it adequately captures the leakoff phenomena; it can also be used to estimate fracture-face resistance and formation permeability. However, it is an iterative technique that requires an initial guess of permeability and the fracture area. Mayerhofer and Economides, in 1997, presented a modified version of the original Mayrhofer model [37]. Their solution suggests using a graphical technique that overcomes the original method's limitations and does not require an initial guess of permeability and fracture-face resistance. Valkó and Economides, in 1999, presented an update of the Mayerhofer model with a dimensionless pressure function, which enabled them to estimate the leakoff rate based on actual fracture length at the end of the time segment [16]. The authors show that the linear leakoff in their proposed model may differ significantly from the leakoff in high permeability formation [38].
section_name: Review of Fluid Leakoff Models

section_content: Our model fundamentally honors the Mayerhofer model using a linear flow Equation (Equation ( 8 )) as the governing equation for the fluid leakoff; however, it differs in how the problem is formulated and solved. The problem is solved continuously rather than discrete pressure drop elements, as depicted by Equation (7).
section_name: Model Formulation

section_content: The boundary condition was formulated based on the mass balance given by Equation ( 9 ). The mass balance implies that fluid is injected to create storage volume and then flowed back as this volume dwindles as fluids are produced and leaked off. 

The fluid loss formation leaks off linearly from the fracture, perpendicular to the fracture faces following Darcy's law as illustrated by Equation (10):

where λ = k/µ. The area of the fracture in this model comprises four fracture faces of bi-wing fracture [39] given by:

The inner boundary condition (Equation ( 12 )) describes the pressure gradient at the fracture walls at each time step. The width change in the model honors the pressure change inside the fracture with constant fracture compliance (Equation ( 13 )) [35]. As the solution is considered a pre-closure model, it is valid to assume constant fracture compliance, which may be calculated using Table 1 for different fracture geometries. Refer to Appendix A for the complete derivation of Equation (12). 

Table 1. Fracture compliance expressions for different fracture geometry models. 

Fracture Geometry PKN KGD Radial
section_name: Boundary Conditions

section_content: The model honors the underlying assumptions of the linear diffusivity equation. The model formulation assumes an injection of slightly compressible fluid at a constant rate into an infinite conductivity fracture (i. e., negligible pressure drop inside the fracture) whose area is already in place (i. e., fracture face area does not change with time). Therefore, the injection and flowback process only causes a change in fracture width. Thus, this assumption limits the model's applicability to the flowback period. The flowback rate is assumed to be constant and starts immediately after stopping the injection. As the flowback starts, the fracture begins to close along the hinge at the tip while the fracture area (Equation ( 11 )) remains constant. Since the fracture area is small and the fluid compressibility is small, their product can be neglected. Therefore, fluid storage/compressibility effects in the wellbore and fracture can be dropped, which is a reasonable assumption.
section_name: Model Assumptions

section_content: The linear diffusivity Equation (Equation ( 8 )) and inner boundary conditions (Equation ( 12 )) were solved analytically by using Laplace transform and given by Equation (14). 

where u is the Laplace parameter, and t f b is the flowback time. The inversion of the Laplace space of Equation ( 14 ) was carried out analytically using the convolution theorem [40]. Equations ( 15 ) and ( 16 ) show the final solution form of the solution. Appendices A-C give detailed steps of the derivation of boundary conditions, Laplace transforms, and numerical inversion verification using Stehfest's algorithm [41], respectively. 

(p

The mathematical examination of the solution during injection suggests increased pressure during pumping. However, since the boundary condition does not include model fracture propagation, it is not useful to use the model during the injection period, as the change in pressure in the model yields only a change in fracture width according to Equation (13). 

Figure 8 shows the solution plotted in dimensionless variables for different values of α from Equation (15). Careful examination of Equation (15) indicates that the solution can be applied to injection/shut-in data when the flowback is set to zero. The resulting equation was reduced to the solution by Siddiqui et al. [39], which corresponds to curve (A) in Figure 8. As the value of α increases from zero to 1, the solution suggested an increase in decline rate of pressure as indicated by curves (B) to an extreme case of (F) where the flowback rate is equal to the injection rate. Further inspection of the physics robustness of the model is to be acquired by eliminating the fracture volume from the inner boundary condition (Equation ( 12 )) by setting fracture compliance to zero. The resultant equation describes a conventional square root of time solution of injection followed by production period, as illustrated by Equation (18).
section_name: Analytical Solution and Mathematical Validation

section_content: The analysis methodology of the proposed approach requires prior knowledge of mobility ratio (reservoir permeability to fluid viscosity), total compressibility, porosity, and fracture compliance which can be calculated based on the fracture geometry model from Table 1. The following steps summarize the analysis workflow:

• Convert time to F-function using Equation ( 16 ).
section_name: Analysis Methodology

section_content: Calculate the derivative of pressure with respect to F-function of Equation ( 15 ). • Plot FdP/dF and (pp r ) versus F on linear graph paper on the same axis.
section_name: •

section_content: Find the linear portion of FdP/dF, which superimposes on a straight line from the origin.
section_name: •

section_content: The pore pressure is calculated when (pp r ) overlays FdP/dF and exhibits a straightline from the origin.
section_name: •

section_content: The closure pressure is calculated when FdP/dF and (pp r ) deviate from linearity. • Equation (15) suggests that the slope of the straight line leads to calculating the fracture area using A = q i c f 2mλϕc t.
section_name: •

section_content: The proposed model is validated using a numerical simulator and against accurate gauge measurement deployed in experimental pump-in/flowback at Stanford underground research facility.
section_name: Model Validation and Field Applications

section_content: This section verifies the proposed approach against the ResFrac commercial simulator. ResFrac is a fully integrated 3D hydraulic fracturing, reservoir, and wellbore simulator; therefore, it is better suited for our purposes, as coupling the fracturing model with a reservoir model is crucial to our model verification. The verification process workflow includes creating a simulation model, mimicking a pressure decline response in the pump-in/flowback test, and then analyzing the pressure decline using the proposed approach. Finally, our analytical model's calculated closure and pore pressure values are compared to the numerically simulated results. Table 2 provides the input parameters for the base case simulation of a fracture injection test in a horizontal well. The F-function was calculated based on Equation ( 16 ) for a time period starting from flowback time. The derivative FdP/dFF is then computed and plotted against F-function as illustrated in Figure 10. The F-function is plotted in a reverse direction so that time increases to the right, and the origin's point lies on the right-hand side. A straight line is easily identified as the derivative line passes through the origin. Equation ( 15 ) suggests that plotting pressure difference against F-function on the same axis would create a straight line that superimposes the derivative line that extends to the origin. The reservoir pressure "pore pressure" can be a matching parameter as implied by Equation ( 15 ). The closure pressure is identified as 50. 6 MPa with an error of less than 1% from the simulator input for minimum horizontal stress (Sh min ), as shown in Table 2. The pore pressure was also calculated as a matching parameter with merely a 2% error from the model input. Table 3 summarizes the analysis results. In this section, we use experimental data to verify the validity of our approach in order to give an accurate estimation of closure pressure. Several injection tests are conducted in this series of experiments by injecting the fluid at sufficient pressure to beak the formation, followed by either a long shut-in period or immediate flowback. The fracture initiation and closure were instrumented with the Step-Rate Injection Method for Fracture In-situ Properties tool (SIMFIP). The tool setup consists of a double-packer hydrofracturing probe with a 6-component displacement Fiber Brag sensor. The high-resolution sensor is placed within the target interval sealed between the two packers and set in the center of the 2. 41 m long interval. The sensor is a 0. 24 m long and 0. 1 m diameter pre-calibrated aluminum cage connected to two 0. 58 m long elements that allow clamping of both ends of the cell on the borehole wall (clamps indicated by red arrows in Figure 11 ). Guglielmi et al., in 2022, presented a detailed description of the procedure adopted in the experiments [42]. During the injection, the tool captures the relative borehole displacement between the upper and lower clamps without influence from the packer system. The tool provides highresolution measurements in the micrometer range, sufficient to capture a few millimeters of borehole displacement. The pressure flow rates are continuously recorded during the test with a half-second resolution [42]. The tests were conducted in 2018 and 2019 during hydrofracturing experiments in boreholes drilled into the highly fractured metamorphic rock of Sanford underground research facility in South Dakota [43, 44]. The boreholes drilled from drifts at the mine were 1490 m (4850 ft) and 1250 m (4100 ft) below ground surface (BGS), respectively. The borehole measuring 1490 m was around 69 m deep, slightly inclined from the horizontal, while the 1250 m borehole was approximately 50 m long. In this study, we present an investigation of the 4100 ft borehole since it includes flowback cycles with corresponding measures of the SIMFIP tool. In total, eight injection tests were conducted to investigate the fracture displacement response at different depths to construct a stress profile at each depth along a 50 m deep borehole. All eight tests across 4100 intervals contained several natural fractures before the experiment, even though it is highly expected that those fractures will be reactivated under pressure during pumping.
section_name: Simulation and Analysis

section_content: In a standard SIMFIP measurement procedure, a series of injection cycles with a small volume of fluid induce fractures in the rock. Typically, the pump-in cycles are followed by several flowback cycles and conclude with a long shut-in to characterize the hydromechanical transient response of the rock. The SIMFIP clamps reflect the rock deformation response during injection and flowback/shut-in periods. The pressure recorded during the experiments usually correlates with the measurements of the SIMFIP tool. However, the correlation is sometimes more complex due to the fractured rock heterogeneity. Ideally, the borehole displacement measured by the tool during injection is reversible while flowback and shut-in. Figure 12 shows a sketch of the displacement measured by the SIMFIP tool versus the recorded pressure during experiments. Similar to the graphical approach to obtain (Sh min ) from normalized tiltmeter response by Economides and Nolte, as in Figure 13, the (Sh min ) is approximated as the linear portion of displacement during shut-in or reopening extends to zero displacements [45]. For fractures perpendicular to (Sh min ), the relationship between the deformation and the change in pressure is linear, as given by Equation ( 19 ) [45, 46]. In this study, we are only considering TV4100 borehole tests 4, 7, and 8. Test 7 is an excellent candidate to validate our proposed model, comprising three injection cycles, including an extended shut-in injection test followed by a flowback test. Figure 14 and Table 4 quantitatively describe the workflow of the test. Fracture Pressure of 3500 psi is reached, followed by a pressure drop to 3244 psi.
section_name: SIMFIP Tool Measurements Analysis

section_content: Increase to constant flowrate of 1. 54 L/min for 30 s.
section_name: 2.

section_content: Increase to constant flowrate of 2. 09 L/min for 120 s. 4. 

Increase to constant flowrate of 2. 79 L/min for 360 s.
section_name: 3.

section_content: Overnight shut-in (14 h 24 min 46 s). 6. 

Bleed off.
section_name: 5.

section_content: Fast pressure increase at a quasi-constant injection flow rate of 2. 78 L/min for 300 s Cycle 3

High step-pressure test Test 7 Cycle 2. 2

In cycle 2. 2, the injection pressure reached a fracturing pressure of 22. 7 MPa, followed by a pressure drop of 22. 2 MPa. Then, bleed off occurred at a target rate of 0. 5 L/min, around (18%) of the injection rate. However, challenges were encountered to maintain a constant flowback rate, as indicated in Figure 15. Therefore, it is anticipated that the fluctuation in pressure response will affect the analysis once the pressure derivative is calculated. Figure 17 shows the analysis of cycle 2. 2 using the proposed approach. The F-function is plotted in the reverse direction (F increases to the left). Since the pressure decay is highly disturbed when the derivative is plotted, (pp r ) is used to identify the slope of the model. The pore pressure for the 4100 borehole is 4. 23 MPa [48]. Table 5 exposes all the input parameters used in the analysis [47] [48] [49]. As illustrated in Figure 17, plotting the pressure difference (pp r ) and Fdp/dF versus F should reveal the straight-line region which extends to the origin. Despite the noisy derivative highlighted by three spikes in Fdp/dF due to fluctuating flowback rate, the straight-line region can be accurately identified once the derivative deviates from linearity. The closure pressure is estimated to be 18. 97 MPa (2752 psi) with less than 0. 1% error from the SIMFIP tool. 

Table 6 shows a summary of the analysis results of cycle 2. 2. To confirm the consistency and the validity of the proposed approach when applied to conventional DFITs, we analyzed injection cycle 2. 1 from test 7, which was concluded with a long overnight shut-in period, as illustrated in Figure 14. In this cycle, 2. 8 L/min of fluid was injected at a constant flow rate for around 12 min, followed by 14 h long shut-in. It was reported that several natural fractures had been reactivated between the SIMFIP tool clamps. During the falloff period, the SIMFIP tool shows a smooth variation of fracture opening and shearing [42]. In Figure 18, the x-intercept of the straight line during closure period indicates a closure pressure of 18. 7 MPa (2712 psi). 

The application of the model on injection/shut-in tests (i. e., DFIT, Minifrac) is obtained by setting α in Equation ( 16 ) (ratio of the flowback rate to the injection rate) to zero. The resulting equation agrees with the model offered by Siddiqui et al. in 2016 [39]. After computing F, Equation ( 17 ) is used to generate Fdp/dF versus F, as in Figure 19. In the analysis of this test, the slope manifests itself when it extends to the origin (the x-axis increases to the left). Since the model straight-line is independently evident, pore pressure can be used as a matching parameter by plotting on the same axis. The resulting pore pressure is in excellent agreement with the reported pore pressure by the EGS collab team [48, 49]. The derivative and pressure difference deviation from linearity indicates a closure pressure of 18. 9 MPa, with less than 2% error from the SIMFIP tool. Table A1 (Appendix D) highlights the analysis results from cycle 2. 1. 

The analysis of the injection cycle 2. 1 was compared to other techniques as further verification of the F-function model. Nolte introduced a dimensionless pressure function known as G-Function that describes the leakoff from pressure decline during injection tests [26, 35]. Later, Barree and Mukherjee developed a new diagnostic plot of GdP/dG versus G known as the tangent method [50, 51]. One of the most recent interpretation techniques of G-function is the compliance method which uses dP/dG behavior to estimate the "contact pressure" and closure pressure [52]. Ideally, the tangent method uses the deviation GdP/dG from linearity as indication closure pressure, while compliance uses dP/dG minimum as the contact pressure, then closure pressure is 75 psi lower than contact pressure as suggested by numerical simulations [53]. Figure 20 demonstrates the application of both the tangent method and the compliance method. The tangent method estimated the closure to be around 18. 8 MPa (2725 psi), while closure from the compliance method does not show a definitive closure because dP/dG does not reach the minimum. The upward deflection occurs at 19. 7 MPa (2895 psi); therefore, the closure would be around 19. 4 MPa (2820 psi), 100 psi higher than the closure by the tangent method. As seen from the previous analysis, the compliance method did not show the typical signature where dP/dG decreases to the minimum and then deflects upward; thereby, the closure pressure cannot precisely be determined. On the other hand, closure pressure from F-function and G-function is within excellent agreement with SIMFIP tool measurement.
section_name: Cycle 2.2

section_content: In this section (test 4), we show an injection test in which G-function fails to deliver a closure signature. In test 4, nine injection cycles were conducted, as shown in Figure 21. The SIMFIP tool anchors are set across (28. 6 m-31. 01 m) displacement intervals. This test is characterized by the reactivation of cemented natural fractures between the tool anchors that would reflect on fracture propagation and the pressure decline analysis. In cycle 4, a reasonably constant flow rate of 1. 2 L/min 10. 5 min was applied followed by 7 h shut-in. Figure 22 illustrates the tool measurement against zone pressure. At the start of the injection, a complex fracture initiation behavior was observed. The falloff period showed a smooth closure pressure at around 18. 94 MPa (2750 psi). In the G-function analysis indicated by Figure 23, the compliance method followed the typical behavior in which dP/dG decreases and reaches a minimum, yielding a contact pressure of 17. 58 MPa (2550 psi) and closure pressure of approximately around 17 MPa (2475 psi). Conversely, GdP/dG monotonically increases, implying that the fracture does not close during the transient period while the fracture is closed based on actual measurements from the SIMFIP tool. Figure 24 shows an analysis of the F-function of test 4 cycle 4. First, the F time was computed, similar to test 7 cycle 2. 1, using Equation ( 16 ) after letting the parameter α be zero. Next, Equations ( 15 ) and ( 17 ) are used to plot (pp r ) and FdP/dF, respectively. 

The pressure derivative FdP/dF decreases to a minimum value and then levels up to follow the model straight-line that extends to the origin and then deflects upward, leaving the straight line and indicating the fracture closure. The deviation from linearity occurs at a pressure of 18. 78 MPa (2725 psi) with less than 25 psi lower than the actual measurements from the SIMFIP tool. The following section will show our approach's application in analyzing multiple flowback cycles. Although no SIMFIP measurements were published for test 8 data from EGS Collab, showing the consistent results of the technique when applied to the same formation is of paramount importance to this paper.
section_name: Test 4 DFIT

section_content: In test 8, three injection cycles were conducted, starting with step pressure tests in the first and last cycles, as seen in Figure 25. Cycle 2 is divided into multiple injection/flowback cycles to evaluate the In-situ principal stress. In cycle 2. 1, the injection was followed by around 25 min of shut-in, then bleeding off at a flow rate of 0. 37 L/min. Cycles 2. 2, 2. 3, and 2. 4 nearly have the same injection volume followed by immediate flowback, making them good candidates for our analysis. In this part of the analysis, we will compare our presented model to mentioned techniques in Section 1. 2 and summarized in Figure 7. Figure 26 shows the analysis of cycles 2. 2, 2. 3, and 2. 4 using a plot of pressure versus time on a linear plot. Nolte and Smith describe the pressure decline during flowback by a characteristic reversal of curvature, and they proposed a point of closure at point (A) [27]. Shlyapobrsky suggested a lower point of closure at the start of the wellbore stiffness straight line (B) [29]. Soliman and Daneshy's model indicates a progressive fracture closure with point (B) as complete fracture closure [32]. The closure pressure, according to the most recent approach, is at point C, where the fracture stiffness straight-line intersects with wellbore stiffness [31]. The same injection cycles were analyzed using F-function, as shown in Figure 27. The data in Table 5 are used to compute F according to Equation (16). FdP/dF and are plotted according to Equations ( 15 ) and ( 16 ). The derivative behavior is similar to the simulation case in Figure 10, in which FdP/dF decreases and then follows the model straight-line that extends to the origin. Due to data granularity, the derivative oscillates and then deviates upwards. Since the flowback rates are relatively constant, the slope was distinct without plotting the pressure difference; pore pressure can be estimated as a matching parameter.
section_name: Test 8

section_content: The main advantage of the flowback test over conventional injection/falloff is the test run time, particularly in tight formations with low leakoff. Various authors have presented tools and techniques to analyze the pressure data from pump-in/flowback tests [27, 29, 31, 32]. However, all of these models adopt a plot of pressure versus time/cumulative flowback rates to identify the characteristics signature of fracture closure. Flowback rates must be appropriately designed and maintained for these models to work. Therefore, the flowback rate plays a significant role in a meaningful flowback test. Nolte and Smith [26, 27, 35] suggested a constant flowback rate no larger than 15% of the injection rate. On the other hand, Clarkson and Williams-Kovacs presented a flowback conceptual model with a medium rate of (1-10%) of the injection rates [54]. Zanganeh et al. recommended a flowback rate of (0. 1%) of the injection rate, followed by an extended shut-in period [55]. 

The model presented in this paper offers a novel diagnostic plot that limits the need for conventional pressure versus time plots. Therefore, it can be used in wide ranges of flowback rates. For instance, methods by Nolte and Plahn rely on the presence of fracture stiffness straight-line to identify the fracture closure. 

In Figure 9, the early sharp decline following the shut-in is due to the dissipation of frictional losses instantaneously because of the injection process. The first couple of minutes after shut-in, a pressure gradient inside the fracture (higher pressure at the fracture mouth) will force the fluid to redistribute, providing energy to propagate the fracture further. The excess energy is then dissipated after a lower pressure area is created at the wellbore because of the flowback. Consequently, the pressure gradient inside the fracture is reversed. This combined effect will limit the possibility of observing the stiffness's straight lines, thus limiting the applicability of Plahn's and Nolte's approaches. 

The analysis of the F-function of the simulation case, Figure 10, exhibits a clear closure signature as the FdP/dF decreases to reach linearity, which extends to the origin. The model assumes a reduction in width under constant area, leading to a masked fracture closure due to fracture growth after shut-in. This caused the straight-line region to last only briefly as the pressure's excess energy contributed to the extension of fractures. FdP/dF and (pp r ) merge to a straight line as fracture growth stops and the fracture closes according to model formulation and assumptions. 

The behavior of FdP/dF and (pp r ) following deviation from linearity is related to after-closure pressure change inside the fracture. Typically, in injection shut-in tests (DFITs), the derivative (FdP/dF) after closure may deflect upward or downward. However, in flowback tests, the derivative deviates upward, and the pressure difference shifts downward; as the substantial pressure drop near the wellbore masks any pressure gradient behavior far from the wellbore. 

Ideally, designing the flowback test requires the flowback rate to be constant. A constant flowback rate is essential to generate a smooth pressure falloff response. If the flowback rate is allowed to vary, the pressure response will be erratic, and the fracture closure signature may be impacted. 

In test 7 cycle 2. 2 (Section of Test 7 Cycle 2. 2), the flowback rate is approximately (15%) of the injection rate. The injection rate was fairly constant at 2. 78 L/min; however, the flowback rate is somewhat variable, as illustrated in Figure 15. Pressure variability between minutes 10 and 15 is reflected in the analysis once the pressure derivative is calculated, as shown in Figure 17. Despite the noisy spikes in the derivative, the closure signature was readily identified when the derivative and pressure difference deviated from linearity. 

In the same test, cycle 2. 1 (Section of Test 7 Cycle 2. 1), the F-function approach was applied to an injection test with a long shut-in period as shown in Figure 19. In this case, the near-wellbore pressure drop effect did not mask the after-closure derivative behavior; therefore, FdP/dF and (pp r ) developed a different after-closure signature. The pressure derivative declines downward with a changing slope, and the pressure difference declines gentler than the model linearity period. The model straight-line suggests that (pp r ) decreases as the fracture loses width under constant area assumption. Thus, the afterclosure behavior implies a remaining system conductivity contributing to the pressure drop after closure. Furthermore, the late (pp r ) slope would generate a larger fracture surface area. 

The tangent method in Figure 20 shows a clear signature of natural fracture opening and about the same closure as the F-function. However, the G-function fails to estimate closure in test 4 as GdP/dG monotonically increases. The compliance method, in contrast, displayed an ideal upward concave signature, but with a 300 psi lower closure pressure as compared to the closure pressure obtained from the SIMFIP tool, as depicted in the Table 7. The EGS collab team reported that the injection period of test 4 did not propagate a new fracture; instead, it reactivated the natural fracture at the wellbore. As a result, tangent and compliance performed poorly as G-function considers fracture propagation with its underlying assumptions. Unlike tangent and compliance methods, F-function showed a clear closure signature as the linearity region manifests itself when it extends to the origin, yielding nearly the same closure as the SIMFIP tool. 

In the 8 flowback cycles of the test, the analysis of the F-function produced consistent results and closure pressures that were closely aligned with those from other techniques. In contrast, the results from the different methods were not consistent. For example, even though the injection scheme for the three analyzed cycles was nearly identical, it was expected that the closure pressure would be similar, or potentially slightly higher in later cycles due to increasing pore pressure from multi-cycle injection. This increase acts as back pressure and helps to support the minimum horizontal stress [56]. As illustrated in Table 7, our new model estimates the closure difference between cycles 2. 2 and 2. 3 to be around 15 psi and between cycles 2. 3 and 2. 4 to be 30 psi. The last cycle closure pressure difference is double that of the earlier pressure because the flowback rate is (15%) larger than the first two cycles. Therefore, the fracture might pinch at the wellbore, leading to a larger closure pressure [2, 4].
section_name: Discussion

section_content: This paper provides a comprehensive review of the available tools and models for analyzing pump-in/flowback test data. Additionally, a new analytical model is introduced that utilizes a linear flow equation and a specially designed boundary condition. This allows the fracture width to decrease as fluids are withdrawn from the system, accelerating the fracture closure process. The model was then verified and validated using experimental and simulated data. Based on the findings in this study, the following conclusions were drawn:

1. 

The fluid leakoff after flowback problem was formulated and solved thoroughly based on a linear flow differential equation instead of relying on stiffness straight-line methods.
section_name: Summary and Conclusions

section_content: The model was verified against a fully 3D fracture model and strain gauge measurements, and exhibits a clear signature with an accurate and consistent assessment of closure pressure in both DFIT and flowback tests. Flowback rate, m 3 /s q L Leakoff rate, m 3 /s q L,j

Leakoff rate at time step j, m 3 /s q L,j-1

Leakoff rate at time step j-1, m At fracture face, the velocity across porous medium is given by:

where (λ = k/m) for a bi-wing fracture, the fracture face area is given by Siddiqui et al. in 2016:

Substituting in the mass balance equation:

Considering the fluid compressibility and chain rule: The mass balance equation then yields:

Thus, arranging the previous equation and omitting the compressibility term because the fracture volume is small relative to the total injected volume. 

The fracture width changes according to the change of net pressure as given by Nolte in 1979 [26] :

The average width is calculated based on Table 1 for different fracture geometries. To appropriately solve the equation continuously, the Heaviside step function is used as in the following figure:

Initial condition: p(x, t = 0) = p r

Removing the Heaviside step function and carrying out the integral:
section_name: 2.

section_name: Acknowledgments:

section_content: The authors recognize HESS Corporation, and Halliburton Energy Service, members of the hydraulic fracturing consortium at the University of Houston, for sponsoring this work.
section_name: Acknowledgments:

section_content: Funding: This research received funding from HESS Corporation, Shell PLC, and Halliburton Energy Service, members of the hydraulic fracturing consortium at the University of Houston. 

Institutional Review Board Statement: Not applicable. 

Informed Consent Statement: Not applicable.

section_content: Data Availability Statement: Data are publicly available on Geothermal Data Repository (https: //gdr. openei. org/submissions/1270).

section_name: Conflicts of Interest:

section_content: The authors declare no conflict of interest.
section_name: Conflicts of Interest:

section_content: The following abbreviations are used in this manuscript:

Let a = ϕc t /λ and b = 2λ √ a ⇒ 4λϕc t Applying Laplace transform to PDE and finding the exponential solution:

After substituting from outer boundary condition:

Using the convolution theorem to invert the above expression
section_name: Abbreviations

section_content: To confirm the inversion of Equation (A13) from Laplace space to real space, Stehfest's algorithm was used to find the numerical inversion that matched the analytical inversion as shown in Figure A2.
section_name: Appendix C. Verification of the Analytical INVERSION

section_name: Appendix D. Summary of Analysis Results of Cycle 2.1