Flowback Test Analyses at the Utah Frontier Observatory for Research in Geothermal Energy (FORGE) Site

In 2017 and 2019, injection testing was carried out in three zones in a vertical well in granite at the Frontier Observatory for Research in Geothermal Energy site near Milford, Utah, USA. In several injection cycles, flowback was implemented rather than shut-in. The goal was to explore an alternative to prolonged shut-in periods for inferring closure stress, formation compressibility, and formation permeability (permeability thickness product). The flowback procedures involved a cyclic flowback/shut-in, while pressure decreased. The flowback data are presented, and analyses are shown. The inferred closure stress(es) from flowback analyses are lower than for equivalent injection cycles that were strictly shut-in. Relatively high formation compressibility obtained from the flowback analyses indicates an extensive, fractured system. This study also includes numerical simulation of the flowback events. The numerical model shows that the rebound pressure is not necessarily the lower bound of the minimum principal stress. The signature of stiffness change can be identified as the process when the depletion mainly transitions from hydraulic fracture to natural fractures from numerical analysis. Overall, flowback potentially has advantages over shut-in because of the reduced time to closure.


Introduction
Enhanced geothermal systems (EGS) offer the potential to bring low-cost geothermal energy to locations that lack natural permeability through hydraulic stimulation (Moore et al. 2019). The U.S. Department of Energy selected a location near Milford, Utah, as the site for the Frontier Observatory for Research in Geothermal Energy (FORGE). The FORGE program aims to develop the techniques required for creating, sustaining, and monitoring EGS reservoirs. In 2017 and 2019, multiple injections were carried out in Well 58-32 at the Utah FORGE site (see Xing et al. 2020b). Well 58-32 was drilled to 2294 m total depth with 46 m of openhole below the production casing shoe. Post-injection measurements were undertaken under shut-in conditions or while flowing back the well.
Flowback closure stress measurement has been used in the petroleum sector for in-situ stress inference. A historical overview for flowback measurements from the petroleum industry is provided in Xing et al. (2020c). Ideally, pumpin/flowback methods can reduce the operational time, and therefore serve as a substitute for unreasonably long shutin periods in Diagnostic Fracture Injection Testing (DFIT). In addition to inferring in-situ stresses, pump-in/flowback methods can also be used to assess formation permeability (permeability thickness product), formation compressibility, and reservoir pressure (Zanganeh et al. 2020). This study presented pump-in/flowback data for Well 58-32 at the FORGE site. Various analytical methods have been applied to infer the closure stress (often treated as minimum in-situ stress), permeability (permeability thickness product), and formation compressibility. Numerical simulation of flowback has also been carried out. A distinct element method simulator, 3DEC, is used to create the numerical model. This paper first describes the injection activities in Well 58-32. Then, the theory underpinning flowback methods is reviewed, especially the so-called stiffness method. After that, FORGE pump-in/flowback data are evaluated to infer the relevant in-situ properties. Finally, numerical flowback simulations are presented and the results are discussed.

Overview of Injection Program at FORGE Site
During 2017 and 2019, injections were carried out in three zones in FORGE Well 58-32. Xing et al. (2020b) summarizes the reservoir characteristics and the injection operations . The three zones are described as follows: • Zone 1 is the barefoot (openhole) section of the hole, extending from the casing shoe at 2248 m measured depth (MD) to the plug back total depth (TD) at 2294 m MD. All depths are reported as 6.6 m above ground level (GL) to be consistent with the rotary kelly bushing (RKB) in the Sept 2017 injection program. For this zone, all gradient calculations were carried out at a depth of 2262 m true vertical depth (TVD) relative to RKB Sept 2017. This zone was stimulated in both 2017 and 2019. • Zone 2 was perforated over 3 m from 2123 to 2126 m MD. The perforating guns were loaded with 30-g charges at 6 shots per foot and 60° phasing. Gradients were calculated using a depth of 2122 m TVD relative to RKB Sept 2017. This zone was picked, because it contained an abundance of pre-existing fractures (determined from the Formation MicroScanner Image log run before casing in 2017) that were anticipated to be near-critically stressed and prone to shear, dilation, or re-opening. • Zone 3 was perforated over 3 m from 2001 to 2004 m MD. The perforating guns were fired at six shots per foot with 30-g charges and 60° phasing. Gradients were calculated at a depth of 2000 m TVD relative to RKB Sept 2017. This zone contained few fractures or they were expected to be not prone to shear and dilation. The conse-quent anticipation was that breakdown would be difficult. This proved to be true.
In each zone, up to nine injection cycles were carried out (seven cycles in Zone 3). These injection activities included pump-in/shut-in, pump-in/flowback, step rate tests, and combinations. These cycles were designed for injection at different rates from 0.001 to 0.039 m 3 /s and to test different injection protocols (microhydraulic fracturing, DFIT measurements, and step rate-step down testing). Among these injection activities, five cycles in Zone 1 and five cycles in Zone 2 incorporated flowback. Xing et al. (2020b) shows the entire chronology (pressure and pumping rate) in each zone. In this paper, the entire injection data of Zone 2 (including pressure and pumping rate in each cycle) are reproduced in Fig. 1, since all the pump-in/flowback case studies are in this zone. Inferred closure stress from pumpin/shut-in and step rate tests for Zone 1 and 2 are tabulated in Fig. 2. The closure stress interpretation in Fig. 2 does not include pump-in/flowback tests. Those interpretations from flowback methods are described in Sect. 4. Some of those gradients may not be representative of the minimum insitu stress. Possible factors affecting the gradients include near-wellbore features and/or losses, self-induced porothermo-elastic back stress in a fractured environment, and dilation/contraction of fractures not aligned with the major horizontal total principal stress (Xing et al. 2020b Xing et al. (2020b) 2017, excluding lower values measured in cycles with lower injection rates in 2017. • The inferred "apparent" stress gradients in Zone 2 are greater than those in Zone 1. • A higher pumping rate/volume gives higher closure stress. Plahn et al. (1995) summarized the use of flowback as a closure stress diagnostic. These authors provided excerpts from relevant publications-those are reproduced here, with attribution. They also stated: "The pump-in/flowback (PIFB) test is frequently used to estimate the closure stress magnitude. The test is attractive because bottomhole pressures during flowback develop a distinct and repeatable signature. This is in contrast to the pump-in/shut-in test where strong indications of fracture closure are rarely seen." Earlier, Nolte and Smith (1981), observed: "If the flow back rate is within the correct range, the resulting pressure decline will show a characteristic reversal in curvature (from positive to negative) at the closure pressure. The accelerated pressure decline at the curvature reversal is due to the flow restriction introduced when the fracture closes." Shlyapobersky et al. (1988) provided a different line of reasoning that is reminiscent of the compliance method in G-function analysis: "The distinct flowback pressure character is due to the increase of frictional pressure in the fracture and/or the decrease of fracture compliance during continuous fracture aperture reduction before the complete mechanical closure occurs. The mechanical fracture closure is the moment at which the fracture storage, V f ∕ p , equals 0. Therefore, this definition of closure suggests using the lower inflection point as an indication of mechanical closure [sic, the point at which wellbore pressure begins a more or less linear decline following the first inflection point]. At mechanical closure, the hydraulic fracture may still retain significant permeability because an incomplete hydraulic fracture closure caused by released formation particles or mismatched fracture faces. This hypothetical fracture behavior is supported by the fact that the slope of the linear pressure decline after fracture closure may be smaller than the slope estimated from the compressibility relation caused by enhanced flow from the fracture into the wellbore." As will be seen later, when actual data are considered, a shut-in following a flowback period leads to a rebound in the pressure (the examples shown later have multiple flowback-shut-in cycles). Nolte (1982) observed about the value of stabilized rebound pressures: "The rebound pressure is the near constant pressure which occurs (following a short period of increasing pressure) after shut-in of the flowback test. This pressure is an important confirmation for the lower bound of the closure pressure, being nearly equal to the closure pressure if the flowback is ended shortly after closure" (see also, Soliman and Daneshy 1991).

Overview of Flowback Methods
Other early references include Tan et al. (1988) and Hsiao and Tsay (1990). Like Shlyapobersky et al. (1988), Wang and Sharma (2020) and  considered the evolution of system stiffness during flowback: "The system stiffness is the response of the well pressure due to fluid content changes resulting from leak-off to the formation or flowback at the surface. It was shown that the pump-in flowback test gives a robust and attractive method for the estimation of the minimum in-situ stress. Also, it was shown that the flowback can be performed with a constant choke rather than a constant flow rate, which simplifies test procedures." Contemporary work, by Raaen and Brudy (2001), also suggested that flowback measurements provide an improved (and lower) measurement of in-situ stress relative to the Fig. 2 Compilation of reasonable stress gradients determined from pump-in/shut-in and step rate tests, not including pumpin/flowback tests. Reproduced from Xing et al. (2020b). Each cycle may have multiple interpretations by different methods labeled as a-d, respectively (e.g., G-function, pressure vs. square root of time, step rate test, and diagnostic plot) measurements of the shut-in type. A highly relevant paper, with excellent field observations and recommendations, is shown by Savitski and Dudley (2011). They used a fixed choke size rather than a constant flowback rate. In hindsight to the measurements described below, their recommendations of reduced inflow rate are significant. In the FORGE program, the smallest available orifice corresponded to a 1/64 inches choke. Even that may have been too aggressive, at least at early times.

Stiffness Analysis of Flowback
During the flowback process, the fracture volume is related to the flowback rate and leak-off into the formation where V f is the fracture volume, Q b is the flowback rate, and Q l is the leak-off rate. Presumably, during the flowback operations in this low permeability formation, the flowback rate Q b is much greater than the leak-off rate. Then, we can simplify Eq. (1) to be Hence, the volume change of the fracture is approximately equal to the flowback volume. The change of bottomhole pressure p w is directly related to the fracture volume change neglecting near-wellbore restrictions where S t is the system stiffness. Combining Eqs.
(2) and (3) gives where V b is the flowback volume. The system stiffness can be expressed as where C w is the wellbore storage coefficient, A f is the fracture area, and s f is the fracture stiffness (derivative of pressure with respect to the average fracture aperture according to, McClure et al. 2019). That is to say that the system stiffness is composed of two components: wellbore storage and fracture stiffness (Shlyapobersky et al. 1988;). The wellbore storage, C w , is assumed as constant. Then, the system stiffness increases as the fracture stiffness increases at fracture mechanical closure. After fracture closure, the fracture stiffness is so large that the system compliance is controlled by the wellbore storage. Therefore, the system stiffness after closure should be constant and close to the stiffness during the injection (Raaen and Brudy 2001;Savitski and Dudley 2011). (1) (3) dp w ∕dV f = S t , (4) dp w ∕dV b = S t ,

Flowback Data Evaluation at FORGE Site
Recognizing the insights of earlier researchers, it was decided to try flowing back-rather than shutting-inon some of the injection cycles that were pumped in 2019 at the FORGE site. Operationally, there was some trial and error, and consequently, the flowback data in all zones evaluated may not be suitable. The data and possible interpretation methods are presented to demonstrate the potential viability of this expedited measurement technique.
As with shut-in data at a minimum (as can be seen from the historical perspective of flowback measurements presented earlier), flowback data can be used to evaluate the closure pressure, formation permeability (permeability thickness product), and formation compressibility. Five cycles in Zone 1 and five cycles in Zone 2 included flowback. Not all of these data are interpretable for closure stress measurements, either because flowback was not started soon enough after shutdown or volumetric flowback rate procedures had not yet been adequately refined on the site for some of the early measurements. When the flowback is started too late after shutdown, the corresponding pressure could already be lower than the closure pressure, preventing inference of the closure stress. Xing et al. (2020a, c) show preliminary interpretations of flowback data at FORGE. Further detailed interpretations of closure stress, permeability (permeability thickness product), and compressibility are carried out here by analyzing three pump-in/flowback cases.

Description of Flowback Tests at FORGE Site
Three cycles with flowback operation for Zone 2 are investigated: Cycle 5, Cycle 7, and Cycle 9 (refer to Fig. 1). In all of the FORGE test cases, flowback was performed with a fixed choke size rather than a constant flow rate. However, the choke size was not kept fixed constant at all times through the test. Instead, the choke size started with 1/64″ and was progressively increased to 4/64″ in 1/64″ increment. The purpose of changing choke size was to study the performance at different choke size, not to maintain a constant flow rate. During this process, each choke size was maintained for some time. In hindsight, constant choke size would have been desirable. In some cases (e.g., Zone 2, Cycle 9), flowback/shut-in sub-cycles were implemented, where a short period of flowback was followed by a short period shut-in, and the process was repeated several times before a supplementary change in choke size. The flowback was recorded in Zone 2 with a turbine meter.
As was indicated, Zone 2 was perforated from 2123 to 2126 m MD. The perforating guns were loaded with 30-g charges at 6 shots per foot and 60° phasing. Stress gradients were calculated for a true vertical depth of 2122 m TVD RKB Sept 2017. For Cycle 5, the treatment entailed pumping Milford city water at 0.013 m 3 /s for about 5 min; 5.2 m 3 of fluid were pumped. After a 10-min shut-in, the well was flowed back through a 1/64″ choke. The choke was beaned up in 1/64-inch increments from 1/64″ to 4/64″ (Fig. 3). There were no shut-in period during the flowback. After 1 h, the flowback rate was too small to measure. A total of 2.8 m 3 were recovered .
Zone 2, Cycle 7 was a step rate/step down cycle before shut-in. 30.2 m 3 were pumped in Cycle 7. After shut-in for 19 min, flowback started through a 1/64″ choke. The choke was beaned up in 1/64″ increments from 1/64″ to 4/64″. After 16.7 m 3 of fluid were recovered, the flow was too small to measure. The pressure and rate data are shown in Fig. 4. Similar to Cycle 5, there was no shut-in period during the flowback.
Cycle 9 was the final injection cycle when treating Zone 2 in 2019. For this injection cycle, Milford city water was pumped at 0.039 m 3 /s for 10 min. In the oil and gas domain, rates of this magnitude may be used for diagnostic fracture injection testing. The well was then shut-in and the pressure dropped (see Fig. 5). After 28 min of shut-in, a controlled flowback program was initiated. At the beginning of the flowback, the choke size was kept at 1/64″, with cyclic flowback and shut-in, as can be seen in Fig. 5. After 16 min, the 1/64″ choke was changed to 2/64″, also with cyclic flowback and shut-in. Then, the choke size was increased to 4/64″ with 1/64″ increment. Finally, the pressure was bled to be zero. About 14.3 m 3 of fluid were recovered out of 29.9 m 3 pumped.

Closure Stress Evaluation Using the Method of Pressure Versus Returned Volume
Following Savitski and Dudley (2011), the closure pressure can be inferred from a plot of pressure versus returned volume, as shown in Fig. 6 for Cycle 9. The closure pressure corresponds to a deviation from linearity in the plot indicating the change of system stiffness. The surface pressure corresponding to apparent closure is 10.3 MPa (blue circle in Fig. 6). A hydrostatic pressure gradient of 9.8 MPa/km is assumed and the total hydrostatic pressure is calculated to be 20.8 MPa (TVD is 2122 m). Therefore, the closure pressure is 31.1 MPa and the closure stress gradient is 14.7 MPa/km. If the intersection point (red circle) of the two linear sections is selected, the surface pressure at closure is 11.0 MPa and the stress gradient is 14.9 MPa/km. As seen in Fig. 6, the system stiffness is high initially and then decreases until the end of the flowback test. However, this is puzzling, because the system stiffness should increase with the fracture closure (refer to Eq. 5). This particular trend of changing stiffness was also seen in Savitski and Dudley (2011), where they hypothesized that the lower stiffness is caused by choked flow when the fracture gets pinched near the wellbore. From Fig. 5, we notice that the rebound pressure did not reach maximum value, indicating that the pressure is not in quasi-equilibrium state and has a transient signal. The transient signal was mitigated through the flowback/shut-in sub-cycles conducted in Cycle 9. If flowback time was shorter and shut-in time was longer in each subcycle, it is anticipated that decreasing stiffness would not be evident.
Although the stiffness decreases with flowback in each sub-cycle with the 2/64″ choke, the initial stiffness increases with more sub-cycles, and it is almost unchanged in the last three sub-cycles (refer to Fig. 7a). There is a pronounced change of initial stiffness between sub-cycle 2 and sub-cycle 3 with the 2/64″ choke, which suggests the occurrence of fracture closure. Also, final sub-cycle's initial stiffness with the 2/64″ choke is 19.7 MPa/m 3 , which is very close to the stiffness during injection (20.4 MPa/m 3 in Fig. 7b). This similarity in stiffness indicates that the fracture was closed. Although the initial stiffness is affected by transient flow response, it still can be used as an indicator for fracture closure, because the near-wellbore restriction is coupled with fracture closure and the transient signal is not dominant. Therefore, the surface pressure at closure could be interpreted as 10.3 MPa (blue circle in Fig. 6) based on the initial stiffness analysis. The corresponding stress gradient is 14.7 MPa/km. Similar to Zone 2, Cycle 9, the flowback stiffness in Zone 2, Cycle 7 decreases gradually as the returned volume increases for the same choke size. This is shown in Fig. 8a. The initial flowback stiffness did not change significantly for the three choke sizes used and they are all smaller than Pressure versus returned volume for Zone 2, Cycle 9. The surface pressure at closure is around 10.3 MPa and the stress gradient is 14.7 MPa/km, presuming that the pressure (blue circle) corresponding to deviation from a linear trend is chosen. If the intersection point (red circle) of the two linear sections is selected, the surface pressure at closure is 11.0 MPa and the stress gradient is 14.9 MPa/km. Learnings include starting the flowback immediately following the shutdown and using shorter flowback/shut-in cycles. This operational precautions would ensure not missing early closure and having a more definitive plot of pressure versus returned volume Stiffness analysis for Zone 2, Cycle 5. a Pressure with returned volume during flowback: comparison of the apparent stiffness for each sub-cycle during flowback with a 2/64″ choke in Zone 2, Cycle 9. The initial stiffness of the first sub-cycle is 9.0 MPa∕m 3 . The initial stiffness for the last cycle is 19.7 MPa∕m 3 . b Pressure with injection volume during pumping: the stiffness during injection is 20.4 MPa∕m 3 the injection stiffness (24.3 MPa/m 3 in Fig. 8b). There is no indication of fracture closure from the stiffness analysis.
For Cycle 5, as shown in Fig. 9a, during the flowback, the stiffness is almost constant on a 1/64″ choke. While flowing back through a 2/64″ choke, the stiffness decreases as the flowback continued. There could be a fracture closure event when changing chokes from 1/64″ to 2/64″, because the initial stiffness of 2/64″ choke is much greater than for flow through 1/64″ choke. However, the initial stiffness (23.6 MPa/m 3 ) is still much smaller than the system stiffness during the injection (35.5 MPa/m 3 , see Fig. 9b), which makes the in-situ stress interpretation from stiffness analysis challenging in this case.
The stiffnesses of the three cycles during injection are 35.5 MPa/m 3 for Cycle 5, 24.3 MPa/m 3 for Cycle 7, and 20.4 MPa/m 3 for Cycle 9. Interestingly, earliest Cycle 5 has the largest value and latest Cycle 9 has the smallest. The system stiffness decreased as more injection cycles were conducted. The stiffness reduction could be due to permanent reservoir permeability enhanced resulted by the previous injections.

Closure Stress Evaluation from a Plot of Rate Normalized Pressure Versus Material Balance Time
There are different flow regimes during flowback. A log-log plot of rate normalized pressure (RNP) with respect to material balance time (MBT) can be used to identify the different flow regimes (Zanganeh et al. 2020).
In an RNP versus MBT log-log plot, Zanganeh et al. (2020) identified sequential regimes: flow regime 1 (FR1) corresponds to wellbore depletion with a unit slope, flow where Q(t) is the cumulative recovered volume at time t. Figure 10 shows a log-log plot of RNP versus MBT for the three investigated cycles in Zone 2. For Cycle 9, three flow regimes can be identified in Fig. 10a. At the start of the flowback, the unit slope indicates that there is still wellbore storage during this period. Then, it transitions into a 1/2 slope, which is the signature of intra-fracture transient linear flow. Finally, the slope returns to a unit value, which would correlate with a boundary-dominated flow for conventional reservoir analysis. Fracture closure is picked at the end of FR2 (1/2 slope) and the beginning of FR3 (unit slope). The corresponding closure pressure is 30.7 MPa, and the closure stress gradient is 14.5 MPa/km. Different flow regimes are also seen in an RNP versus MBT log-log plot for Zone 2, Cycle 7 (refer to Fig. 10b). The initial slope is unity in the diagram, indicating that there is still a wellbore storage effect at the beginning of the flowback. After this period, it transitions to a 1/2 slope, which indicates transient linear flow. Finally, the slope becomes unit again, which suggests a boundary-dominated flow. Similar to Cycle 9, fracture closure is picked at the end of FR2 (1/2 slope) and the beginning of FR3 (unit slope). The corresponding closure pressure is 30.2 MPa, and the closure stress gradient is 14.3 MPa/km.
The RNP versus MBT plot for Zone 2, Cycle 5 is shown in Fig. 10c. From this plot, only a late unit slope is observed, which is an indication of a boundary-dominated flow. Therefore, we are not able to identify the fracture closure for this cycle.

Fig. 10
Flow regime transitions during the flowback test for the three investigated cycles in Zone 2. a Cycle 9: the initial flow regime is FR1 (wellbore depletion with a unit slope), then transitions to FR2 (linear flow with a 1/2 slope), and finally changes to FR3 (boundarydominated flow with a unit slope). Fracture closure is picked at the end of FR2 (1/2 slope) and the beginning of FR3 (unit slope). The corresponding closure pressure is 30.7 MPa, and closure stress gradient is 14.5 MPa/km. b Cycle 7: fracture closure picking is similar to Cycle 9. The corresponding closure pressure is 30.2 MPa, and the closure stress gradient is 14.3 MPa/km. c Cycle 5: we only observe a late unit slope, which is an indication of boundary-dominated flow. Therefore, we are not able to identify the fracture closure for Cycle 5 ▸

Reservoir Compressibility Evaluation Through Flowing Material Balance Analysis
If a boundary-dominated flow is observed, the flowing material balance method can be used to estimate the "initial water-in-place in the fluid bank region" (Zanganeh et al. 2020). For the boundary-dominated flow, the flowing material balance equation can be written in the form of a straight line as (Clarkson and Williams-Kovacs 2019) where PNR denotes pressure normalized rate, m f is the slope of the flow material balance plot, b f is the intercept, and c t is the compressibility including both the water compressibility and the formation compressibility. Water compressibility at 200 °C and 20.7 MPa is 8.7×10 −10 1/Pa; rock compressibility is approximately 1.82×10 −11 1/Pa. The flowing material balance diagram is generated by plotting PNR versus . Then, the initial water-in-place of the fluid bank region, I WF , and productivity index, J f can be obtained as In this low permeability, low porosity, granitic reservoir, the calculated initial water-in-place should be comparable to the injected volume. The formation compressibility is not precisely known and requires an initial guess to calculate "initial" water in place. The two end members of the formation compressibility are 1.82×10 −11 1/Pa, which is the compressibility of granite, and 2.2×10 −8 1/Pa, which is a highend estimate of unpropped fracture compressibility obtained from core plug measurements performed on reservoir samples in the laboratory (Zhang et al. 2018). If the "predicted" initial water-in-place equals the injected volume, it gives a reasonable upper-limit for the formation compressibility. Figure 11a shows the flowing material balance diagram for sub-cycle 5 on the 2/64″ choke for Zone 2, Cycle 9. In this plot, formation compressibility c FR was taken as 1.5×10 −8 1/Pa, which is a high-side estimate of unpropped Fig. 11 Flowing material balance diagram for the three investigated cycles for Zone 2. For all the three cycles, formation compressibility c FR is assumed as 1.5×10 −8 1/Pa. a Cycle 9: sub-cycle 5 with a 2/64″ choke. The initial water-in-place is calculated as 27.5 m 3 , which is close to the injected volume of 29.9 m 3 . b Cycle 7: the section of flowback is taken from section on a 2/64″ choke. The calculated initial water-in-place is 21.6 m 3 , which is somewhat comparable to the injected volume of 30.2 m 3 . The discrepancy in these volumes suggests that the assessed formation compressibility varies from the assumption of Zone 2, Cycle 9. c Cycle 5: the processed data are taken from section with a 2/64″ choke. The calculated initial waterin-place is 8.9 m 3 , which is of the same order and comparable to the injected volume of 5.2 m 3 ▸ fracture compressibility. The initial water-in-place is calculated as 27.5 m 3 , which is very close to the injected volume of 29.9 m 3 . This flowing material balance diagram shows that the formation compressibility is at the high end. The implication is that there are abundant natural fractures that exist and have been "stimulated" in the reservoir. This is consistent with the Formation MicroScanner Image (FMI) log that suggests Zone 2 has abundance of natural fractures (Xing et al. 2020b).
The flowing material balance diagram for Zone 2, Cycle 7 is shown in Fig. 11b. In this diagram, the formation compressibility is also assumed to be 1.5×10 −8 1/Pa, the same for Zone 2, Cycle 9. The calculated "initial" water-in-place is 21.6 m 3 , which is comparable to the injected volume of 30.2 m 3 . Matching a calculated "initial" water-in-place with the injected volume, the formation compressibility is 1.0×10 −8 1/Pa, still at the high end. This elevated inferred formation compressibility for Cycle 7 also suggests the reservoir has abundant natural fractures.
The flowing material balance plot for Zone 2, Cycle 5 is shown in Fig. 11c. The formation compressibility is assumed to be 1.5×10 −8 1/Pa to be consistent with Zone 2, Cycle 9. The calculated initial water-in-place is 8.9 m 3 , which is in the same order and comparable to the injected volume of 5.2 m 3 . By matching the "initial" water-in-place with the injected volume, the formation compressibility should be 2.5×10 −8 1/Pa, at the high end.

Permeability Thickness Product Evaluation with the Method of Two-Rate Flow
The flowback data can also be used to calculate permeability thickness product using two-rate superposition concepts. Figure 12 shows a two-rate example taken from the flowback period for Zone 2, Cycle 9. A two-rate superposition for this period is realized by a plot of pressure p w versus log t+Δt � Δt � + q 2 q 1 logΔt � (see Fig. 13). Here, q 1 is the injection rate prior to rate change, q 2 is the rate after rate change, t is the time duration of q 1 , and Δt � is the time measured from the instant of the rate change. The permeability thickness product kh can be calculated as (Eq. 6.9 in Matthews and Russell 1967) where k is the permeability and h is the formation thickness. m is the slope, as shown in Fig. 13. The formation volume factor B is defined as the ratio of reservoir volume to surface tank volume. The compressibility of water at surface condition (65 • C and 1 atm) is 5.0×10 −10 1/Pa, while the compressibility at the reservoir condition (200 °C and 20.7 MPa) is 8.7×10 −10 1/Pa. Hence, the formation volume factor is calculated as 0.6. The viscosity is approximated as 1.5×10 −4 Pa ⋅ s at 200 °C and 20.7 MPa. This method offers potential and can presumably be refined by considering partial completion skin and fracture skin effects. logΔt � for the two flow rate tests in Zone 2, Cycle 9. The slope m is 3.6 MPa. Several representative data points from Fig. 12 are used to construct this plot. q 1 , the pressure prior to rate change, is 0.0031 m 3 /s, and q 2 is 0 m 3 /s The pressure build-up before fracture closure is not controlled by matrix permeability but by the near-wellbore connectivity, while it is controlled by matrix permeability after fracture closure, as shown in Eq. 11. Although the pressure build-up before fracture closure is not controlled by matrix permeability, the permeability calculated by Eq. 11 reflects the wellbore connectivity coupled with the fracture closure process. The permeability calculated by Eq. 11 can capture the change of the trend before and after fracture closure and hence identify the fracture closure.
Following the same method described above, we calculated each sub-cycle's permeability thickness product for the 2/64″ choke of Cycle 9. The results are shown in Fig. 14. The permeability thickness product decreases with the decreasing pressure. There is a bi-linear relation between permeability thickness product and pressure drop, which also indicates possible fracture closure between sub-cycle 2 and sub-cycle 3, similar to the stiffness analysis method.
Unlike the case of Cycle 9, there is no shut-in period during the flowback in Cycle 7. Alternatively, flowback period on a 1/64″ choke was divided into small sections (every 50 s). Thus, each small section has a different flow rate and the two-rate method is still applicable. The camdlculated permeability thickness product versus pressure drop is shown in Fig. 15. We can see that the permeability thickness product is around 32.9 md ⋅ m and there is no obvious trend change in this plot for Zone 2, Cycle 7, which could be due to starting flowback late.
Similar to Cycle 7, the two-rate method is also used to calculate the permeability thickness product by dividing the flowback period on the 1/64″ choke into small sections for Cycle 5. From the permeability thickness product versus pressure drop curve in Fig. 16, there is not a clear trend for Cycle 5. The permeability thickness product is around 9 md ⋅ m.

Permeability Evaluation with the Multiple-Rates Flow Method for Cycle 9
It is also possible to do a multiple cycle analysis to obtain the permeability thickness product using a crossplot of the RNP and the Odeh-Jones time function (Odeh and Jones 1965) Fig. 14 Calculated permeability thickness product versus pressure drop for each sub-cycle with a 2/64″ choke in Zone 2, Cycle 9. Each "dot" here represents one sub-cycle. The bi-linear behavior indicates the fracture closure between sub-cycle 2 and sub-cycle 3 where q k is the flowback rate for the kth step, and t k is the time of the kth step rate since the initiation of flowback. However, for Cycle 9, there were shut-in periods between each flowback rate, which makes both the RNP and the Odeh-Jones time infinite. Hence, a very small flowback rate is assumed during the shut-in period. Figure 17 demonstrates a multiple rate analysis for Zone 2, Cycle 9.
The slope of the multiple rate analysis is obtained as m x = 1236.5 MPa/(m 3 ∕s) from Fig. 18. The permeability thickness product can be calculated as The formation volume factor is also taken as 0.6 here. This calculated permeability thickness product value is smaller than that calculated using Matthew and Russell's two-rate method. This could be due to the mathematical difficulties of handling the shut-in period in the multiple rates method.

Summary of the Case Studies and Discussion
From the analyses of the three cases, the closure stress gradient is estimated at 14.3-14.9 MPa/km, the formation permeability thickness product at 5.5-36.5 md ⋅ m, and the formation compressibility is bracketed between 1.0×10 −8 and 2.5×10 −8 1/Pa for Zone 2. The high formation compressibility suggests that this zone has abundant compressible natural fractures, which is in consistent with the FMI log.
As discussed above, an estimate of the closure stress was obtained from two of the three cases (Zone 2, Cycles 9 and 7). For Zone 2, Cycle 9, closure stress is estimated using a stiffness analysis, a RNP versus MBT plot, and permeability change, while for Zone 2, Cycle 7, it was only possible to identify the fracture closure from an RNP versus MBT plot. For Cycle 9, all the interpretations from flowback data give similar closure stress, indicating that the dataset still gives us reasonable estimation of fracture closure due to the mitigated transient pressure signal.
It is appropriate to compare the closure stress obtained from flowback with the estimations from extended shut-in and step rate tests. Figure 19 shows the pressure-time data for Zone 2, Cycle 4, during 2019, one of the typical extended shut-in (DFIT) tests at FORGE; conventional closure stress gradient interpretation suggests a gradient of 18.1 MPa/km. The DFIT data (Cycle 4) in this zone do not suggest a severe choking during the fracture closing. The closure stress gradient from shut-in and step rate tests for Zone 2 ranges from 17.2-21.5 MPa/km (refer to Fig. 2), which is substantially higher than the inferred values for flowback tests, 14.3-14.9 MPa/km. This could suggest that 1) when analyzing flowback data (Fig. 6, for example), an artificial gradient is being picked due to the fact that the flowback started late, or 2) closure is masked by early near-wellbore closure, or 3) flowback offers a very useful method for closure stress interpretation in naturally fractured reservoirs where there is tortuous communication between the wellbore and a natural fracture system. The pump-in/flowback method minimizes the effects of natural fractures, because the flowback rate is much larger than after shut-in leak-off and has a much (13) kh = B 4 m x = 5.9 md ⋅ m. Fig. 17 Plot of a "multiple flow rate" test, taken from the 7590-8690 s during injection in Zone 2, Cycle 9, 2019. The first flowback rate q 1 is 0.0031 m 3 /s and the second flowback rate q 2 is 0.0 m 3 /s and the third flowback rate q 3 is 0.0028 m 3 /s Fig. 18 RNP versus Odeh-Jones time for the multiple rate tests in Zone 2, Cycle 9. The slope m x is used to infer the permeability thickness product in a conventional radial flow relationship stronger signature. In addition, the role of temperature is still unknown during the interpretation of pump-in/shut-in data or pump-in/flowback data, recognizing the 200 °C static formation temperature. As for the first scenario, it is possible that the flowback was not started soon enough in the studies presented. If that is the case, the closure stress picked from a pressure versus returned volume curve may not adequately represent the whole trend. This could result in an underestimation of the closure stress.

Numerical Simulation
We have analyzed the field flowback data of FORGE site to infer the closure stress, permeability, and compressibility of the reservoir. It is also desirable to investigate fracture closure process and pressure response of the reservoir under flowback operation through numerical modeling. The numerical model considers the interaction between hydraulic fracture and discrete fracture network (DFN).

Theory and Background of the Numerical Method
The numerical simulation of flowback was carried out using Itasca Consulting Group's code 3DEC. 3DEC is a threedimensional numerical program based on the distinct element method for discontinuum modeling (Cundall 1971). The discontinuous medium is represented as an assemblage of discrete blocks. The discontinuities are treated as boundary conditions between blocks. Individual blocks behave as either rigid or deformable material. Deformable blocks are subdivided into a mesh of finite difference elements, and each element responds according to a prescribed linear or nonlinear stress-strain law. The relative motion of the discontinuities is also governed by linear or nonlinear force-displacement relations for movement in both the normal and shear directions. In 3DEC, the contacts between blocks are treated as joints. The basic joint constitutive model incorporated in 3DEC is the generalization of the Coulomb slip law. The Coulomb slip model provides a linear representation of joint stiffness (both normal and shear) and a yield limit, and is based on elastic stiffness, frictional, cohesive, tensile strength, and dilation. The model simulates displacement weakening of the joint by loss of cohesive and tensile strength at the onset of shear or tensile failure. Shear and normal stresses on a joint develop elastically until the stress reaches the peak strength. The peak shear stress is given by where c is the cohesion, ′ n is the effective normal stress, and is the friction angle. When the peak shear strength is exceeded, the shear strength drops instantaneously to the residual shear strength. The shear stress can then not exceed the residual strength. The residual strength is given by where c res is the residual cohesion, and res is the residual friction angle. By default in 3DEC, c res = 0 and res = .
A fluid flow model in the joints is an essential part of these simulation. The relation between fluid flow rate per unit width and the pressure in the joint (fracture) can be expressed as, using a conventional cubic relationship where u h0 is the joint aperture at initial state (in-situ stress state); and Δu n is the joint normal displacement (positive denoting opening). Joint aperture is also bound by a minimum value, u res . The mechanical closure does not affect the contact permeability. The continuity equation for a slightly compressible fluid in a deformable rock fracture has the form where K w is the bulk modulus of the fluid. Combining Eqs. (16) and (18), we obtain the expression which is the diffusion equation for the pressure inside the fracture. In the general case, the fracture aperture u h is not only a function of the local fracture stiffness and fluid pressure, but also the mechanical response of larger volume of the rock mass. Therefore, the term (K w ∕u)( u∕ t) is external, or a "source", which is determined by coupling between fluid flow and mechanical deformation of the solid model.
The fluid pressure in the fracture affects the fracture normal deformation. Normal deformation of the fracture is a function of effective stress, which is a linear combination of the confining stress (or total normal stress) in the rock, n , and the fluid pressure, p.
Pressures are calculated and stored in a data structure called a flow-knot. The flow-knot pressures are updated, taking into account the mass balance and possible changes in flow-knot volume due to the incremental motion of the surrounding blocks. The new flow-knot pressure becomes where p 0 is the flow-knot pressure in the previous timestep; Q s is the sum of flow rates into the flow-knot from all surrounding contacts; V and V 0 are the new and old flow-knot volumes, respectively.

Numerical Modeling of Flowback at FORGE Reservoir
The domain of the model is a 300×300×300 m cube. A FORGE-derived deterministic DFN is incorporated in the model (see Fig. 20). The DFN contains three sets of natural fractures, about 2000 natural fractures in total (Finnila et al. 2019). In 3DEC, fractures can only be generated in the preexisting joints. Therefore, the hydraulic fracture path has In the numerical model, wellbore storage is not considered. Therefore, the stiffness/compliance transition from combination of fracture and wellbore storage to wellbore storage only is not captured. However, the process of fracture closure and influence of DFN is studied and presented. Also, the model was calibrated by the injection data. The modeling of flowback without wellbore storage will shed light on relation between fracture closure process and reservoir pressure response. The numerical model aims to investigate the reservoir pressure change under flowback operation considering the influence of DFN. It helps us to understand the pressure change to the fractures closure only and the process of fractures closure under flowback.
The initial aperture of the DFN and minimum aperture are both set as 1 ×10 −4 m. The injection point is at the center of this model. The fluid injection schedule is the same with Zone 2, Cycle 9: it was injected at 1.3×10 −2 m 3 ∕s for 210 sec, and then injected at 3.9×10 −2 m 3 ∕s for 670 sec. The simulated injection pressure match well with the field data of Zone 2, Cycle 9 (refer to Fig. 21). After injection, flowback immediately initiated at 5.5×10 −3 m 3 ∕s for 390 s and pressure dropped during this period. After that, the well was shut in for 270 s and the pressure rebounded. The entire pressure history of the numerical model is shown in Fig. 22.
We observed that the rebound pressure at the end of shutin is lower than h min . This supports the argument in Plahn et al. (1995) that the rebound pressure is not necessarily the lower bound of closure pressure. The pressure distribution and hydraulic aperture of the three stages are shown in Figs. 23 and 24. From these figures, we can see that, at the end of injection, the area around the injection point (the center of the model) has the maximum fluid pressure and hydraulic aperture. At the end of flowback, the pressure around the injection point is reduced to lower than the elements far away; the hydraulic aperture is also reduced to the same level as the other surrounding elements. At the end of shut-in, the pressure around the injection point recovered. We also realized that although the pressure at some region is reduced to below the h min at the end of flowback, the fractures are not fully closed there. That could be due to the tortuous connection between the hydraulic fracture and the natural fractures.
The pressure versus returned volume during flowback of the simulation is shown in Fig. 25. The inferred closure pressure is at the intersection of the two linear section, about 34.8 MPa, which is lower than h min , 37 MPa. The actual h min is between the two linear sections in plot of pressure with returned volume (between point A and C in Fig. 25). The stiffness decreases faster at the beginning, and then keeps almost constant when pressure drops below 34.3 MPa. As shown in the figure, the fractures (both the hydraulic fracture and DFN fractures) are only partially closed during the flowback, which could induce pressure transient signal. However, pressure is not evenly distributed between the main hydraulic fracture and the stimulated DFN fractures. The stiffness change could still reflect the difference between the main hydraulic fracture and the DFN fractures. From the simulation results, we can see that, initially, the main hydraulic fracture has a much higher pressure but a relatively small volume. Therefore, the signal of larger stiffness from the main hydraulic fracture dominates initially. Shortly afterward, the pressure of the main hydraulic fracture is reduced to the same level as the surrounding failed natural fractures, and at the same time, the main hydraulic fracture is not fully closed. From then on, the signal from DFN fractures dominates. Since the DFN fractures have smaller pressure and large volume, the stiffness of later stage of flowback is smaller than the initial stage of flowback dominated by the main hydraulic fracture. Therefore, the total stiffness decreases when the depletion mainly transitions from hydraulic fracture to natural fracture network. The signature of stiffness change could be different for the field data due to the wellbore storage effect.
From the numerical results, we can see the complexity of fracture closure process during flowback. There is even further complexity because of changes in temperature. The near-wellbore restrictions are not addressed in this model. Fig. 23 Pressure distribution of the numerical simulation at different stages. Top is in the plane perpendicular to H max ; bottom is in the plane perpendicular to h min Fig. 24 Hydraulic aperture of the numerical simulation at different stages. Top is in the plane perpendicular to H max ; bottom is in the plane perpendicular to h min It can be done through a fine grid and pressure-dependent fracture conductivity in the future work.

Discussion and Conclusions
Several field cases with flowback data were analyzed. These data were from injections in Zone 2 of Well 58-32 at the FORGE site. Through various methods, the reservoir in-situ information inferred were closure stress gradient 14.3-14.9 MPa/km, permeability thickness product of 5.5-36.5 md ⋅ m, and formation compressibility of 1.0×10 −8 -2.5×10 −8 1/Pa. An evaluation of the closure stress was possible for two of the case studies Zone 2, Cycle 9 and Cycle 7. For Zone 2, Cycle 5, there was no significant indication of fracture closure using either the stiffness analysis or an RNP versus MBT plot, because the flowback was started too late in this cycle. The horizontal minimum stress gradient determined by flowback analysis for Zone 2 (14.3-14.9 MPa/km) is notably smaller than the values from the extended shutin and step rate analyses which ranged from 17.2-21.5 MPa/km. That could be attributed to either starting flowback late or the impact of natural fracture networks on the extended shut-in interpretations, or near-wellbore closure masking pressure equilibrium in the bulk of the fracture system.
The permeability thickness product obtained from the flowback data is on the order of 10 md ⋅ m, which is consistent with the value inferred using after closure analysis following conventional DFIT shut-in practices. The high formation compressibility inferred from flowing material balance method indicates that this zone has abundant compressible natural fractures, which is consistent with the FMI log. These evaluated values of permeability thickness product and formation compressibility are still under evaluation.
There may be alternative interpretations-especially for the closure stress-if the flowback had been started earlier.
Regardless, flowback seems to be a promising methodology with significant operational advantages in terms of time. The measurements are slightly more complicated than simple shut-ins because some form of continuous recording of flowback rate is necessary. Lessons learned were: (1) shorter duration flowback/shut-in cycles with a smaller choke size would be desirable; (2) it may be prudent to start flowback as soon as feasible after shutdown; (3) ensure the rebound pressure reach the peak value. In addition, the smallest practical choke size should be used to quickly achieve a quasi-equilibrium state and the choke size should be kept unchanged during the flowback.
Finally, numerical simulation of the flowback in Zone 2 was carried out using a distinct element method simulator. In the numerical model, the rebound pressure is lower than the h min , which demonstrates that the rebound pressure should be used very cautiously as a closure stress indicator. During the flowback, the system stiffness decreases when depletion transitions from the hydraulic fracture to the natural fracture network.