ABSTRACT
Fiber-optic sensing technologies allow petroleum engineering teams to detect hydraulic fracture interaction with boreholes during unconventional reservoir stimulation. In combination with high-repeatability seismic sources, the same distributed acoustic sensors (DASs) enable vertical seismic profiling (VSP) of the fracture evolution away from the boreholes. We discovered clear signatures of seismic scattering on activated fractures during nine days of continuous seismic monitoring of the fracturing stages at the Austin Chalk/Eagle Ford Field Laboratory. The present study applies a novel approach for quantitative analysis of the scattering events in terms of the evolution of the geometry and elastic stiffness of individual fractures. Our characterization strategy sequentially refines the fracture models: from a stack of 1D soft layers to 3D rectangular inclusions. First, we estimate the number of fracture locations and reflectivity using a modified sparse-spike deconvolution of the stacked VSP traces. The fracture set consists of five fractures spaced by 15–30 m with a reflectivity of approximately 1%. Then, we develop a scattering integral method to refine these estimates along with an inversion of the fracture top and bottom for each monitoring vintage. We find that, initially, some of the fractures are located above the monitoring fiber with the height of approximately 100 m. Then we integrate the seismic interpretation with the low-frequency DAS and pressure and microseismic monitoring to reconstruct the activation process of the fractures. Most likely, some of the natural fractures slowly grew downward to the monitoring fiber as a result of fluid injections in the stimulated well. This led to bright strain anomalies but did not trigger seismicity. The top of the fractures remained almost constant and were limited by a lithologic boundary/stress barrier. To our knowledge, this is the first time VSP data enabled tracking of the fracture evolution with such high spatial and temporal resolution, which was previously only available for crosswell surveys and at a much smaller scale.
INTRODUCTION
Hydraulic fracturing is the enabling technology for unconventional petroleum reservoir production (Zoback and Kohli, 2019b). Despite its ubiquitous application, control and optimization of hydraulic fracturing remain areas of active research. And imaging of the fracture growth and other hydromechanical reservoir processes remains one of the key aspects of this research, which was prioritized in the focused field experiments (Birkholzer et al., 2021; Ciezobka, 2021). Although microseismic monitoring is one of the key technologies to monitor hydraulic fracturing (Mawell, 2014; Eaton, 2018; Zoback and Kohli, 2019a), phenomenologically, it attributes the spatial extent of microseismic clusters and their density to a stimulated reservoir volume (SRV). The SRV consists of activated natural fractures and matrix rocks that are connected to newly created hydraulic fractures (Mayerhofer et al., 2010). This approach oversimplifies the sophisticated mechanics of earthquakes (Eisner et al., 2010). It is known to result in misleading interpretations when aseismic fault slip and stress transfer play an important role (Scholz, 2019).
Deformation sensors provide a more direct observation of hydraulic fracture growth. Originally, operators deployed tiltmeters in the available vertical boreholes around the fracturing stages to estimate fracture height (Fisher and Warpinski, 2012). Rapid advancement of fiber-optic sensing technology (Hartog, 2017) enabled in- and crosswell imaging of the fracture growth (Jin and Roy, 2017), which allowed petroleum engineering teams to detect fracture interaction with the boreholes. However, the strain measurements may not be interpreted in terms of the shape of the entire fracture (Sherman et al., 2019; Liu et al., 2021): the recorded anomalies are mainly sensitive to the fracture edges closest to the fiber, where the strain is concentrated, and thus these measurements are limited to characterization of only a small fraction of the fracture.
Borehole seismic may use the same optical fiber cables as distributed acoustic sensors (DASs) to map the fractures in the subsurface. Several vertical seismic profiling (VSP) surveys reported seismic anomalies attributed to hydraulic fracturing operations (Meek et al., 2019; Binder et al., 2020; Zhao et al., 2021). These studies reported a low-velocity zone along stimulated boreholes. Similar to the concept of SRV, the mapped velocity zones were continuous: a combined effect of pressurized fractures and microcracks in the reservoir rocks smeared by band-limited seismic signals. However, Titov et al. (2021, 2022) report discrete scattering events from several activated fractures and propose a method for estimation of the fracture height based on the scattering amplitude distribution along the fiber.
Our study proposes a new approach for advanced quantitative interpretation of the VSP monitoring of hydraulic fracturing using a comprehensive monitoring data set from the Austin Chalk/Eagle Ford Field Laboratory (ACEFFL) (Hill et al., 2020). The experiment featured fiber-optic static strain and temperature sensing as well as DAS in two horizontal boreholes (Saw et al., 2023) along with several in-well pressure gauges and microseismic monitoring using a surface geophone array. We use these data sets to support interpretation of nine days of continuous seismic monitoring using two permanent seismic sources (Correa et al., 2024). The hourly seismic vintages contain scattering events throughout the monitoring well. A peculiar feature of the data set is that seismic scattering and low-frequency DAS anomalies (Zhu et al., 2023) show that the fractures became activated for two to three stages adjacent to the stimulated interval (Figure 1).
This paper aims to explain the observed seismic signatures in terms of parameters of the activated natural fractures. First, we summarize the observations to develop a conceptual model of the stimulated reservoir: a set of discrete natural fractures. Then we prepare a rock-physics basis for the estimation of fracture locations and reflectivity using linear slip formalism (Schoenberg, 1980). To resolve individual fractures in the interfering fracture scattering events, we first apply a 1D modification of the sparse-spike deconvolution, which is valid for infinite fractures. To estimate the fractures’ vertical extent and position relative to the DAS cable, we incorporate the diffraction effects using the Born approximation. Finally, we integrate the seismic estimates with temporal and spatial distribution of strain anomalies, pressure measurements, and induced seismicity to reconstruct the hydromechanical processes underlying the observed geophysical response. To our knowledge, this is the first time VSP data enabled tracking of the fracture evolution with such high spatial and temporal resolution, which was previously only available for crosswell surveys and at a much smaller scale.
FIELD DATA SET
Outline of the ACEFFL experiment
ACEFFL is a research test-bed funded by the U.S. Department of Energy and SM Energy (Hill et al., 2020). Similar to prior pilot studies, Hydraulic Fracture Test Site 1 (Birkholzer et al., 2021), and Hydraulic Fracture Test Site 2 (Ciezobka, 2021), ACEFFL deployed a comprehensive subsurface monitoring system at a test site in Webb County, Texas, to image the process of unconventional reservoir stimulation and hydraulic fracture growth. The experiment involved stimulation of three boreholes, as shown in Figure 2a. Two of these boreholes, well 3 and well 5, were instrumented with fiber-optic cables that enabled strain and temperature monitoring (Saw et al., 2023) as well as detection of fracture activation by low-frequency DAS (Zhu et al., 2023) and acquisition of continuous seismograms (Figure 3).
The stimulation produced more than 2000 microearthquakes, which were detected and located by an array of geophones in shallow boreholes. The microseismic clusters are aligned orthogonal to the estimated minimum horizontal stress. In addition, one may clearly identify large clusters of strong microearthquakes that extend farther away from the stimulated boreholes, which correspond to natural faults and fractures. These produced numerous strong scattering events that were identified during wavefield separation for microseismic imaging by Ma et al. (2024). The fracture intersections also were clear in the borehole imaging logs in wells 5 and 7 (Figure 4a). The orientation of the conductive natural fractures tends to be close to orthogonal to the borehole trajectories, similar to the hydraulic fractures.
Continuous active seismic data were acquired throughout the first nine days of stimulation of well 5 using permanently deployed surface orbital vibrators (SOVs) and DASs in wells 3 and 5 (Correa et al., 2024). SOVs and DASs proved successful in many subsurface monitoring projects with remarkable sensitivity to small changes in the properties of target reservoirs (Kasahara et al., 2013; Pevzner et al., 2021, 2022; Tsuji et al., 2022). Most of the data processing was automated and done on site so that the data analysis focused on characterization of the time-lapse changes. During the stimulation, SOV-3 and SOV-5 (Figure 2a) generated a monitoring vintage every hour for each. The noise level in well 5 precluded quantitative analysis of the SOV data, so we focused on well 3.
The fluid injections in well 5 induced very few microseismic events in the vicinity of well 3 during nine days of seismic monitoring (Figure 2b and 2c). If new hydraulic fractures had grown from well 5 to well 3, we would have detected some seismicity produced by the strong overpressure and stress perturbation concentrated around the fracture edges. The subsequent stimulation of these same intervals in well 3 generated relatively strong seismicity.
Despite the lack of seismic events, Correa et al. (2024) find strong fracture scattering events in the processed VSP seismograms (Figure 3), which were clearly visible prior to the stimulation and then persisted in the data throughout the nine days of monitoring. Evolution of the amplitudes of the scattering events is somewhat consistent with the fracture interactions recorded by low-frequency DASs (see Figure 1). We discuss selected seismic signatures in the next subsection.
Seismic signatures of fractures
Figure 3b identifies all of the scattering modes in the processed seismogram for SOV-3 to well 3. “Reflected” PS waves, reflected PP waves, and “transmitted” PS waves are visible for multiple fracture sets at different measured depths (MDs): MD 4300–4400 m, MD 4400–4600 m, MD 4700–4800 m, and MD 4850–5000 m. Some dimensions of the fracture set (fracture height and spacing between individual fractures) are comparable to the wavelengths of the seismic signals, and thus the use of the terms “reflection” and “transmission” may be inaccurate. The strongest signals correspond to the reflected PS waves.
Strong P-to-S reflection and transmission coefficients at the surface of weak fractures are anticipated for the steep incidence angles that we have for SOV-3 and the toe region of well 3, approximately 60° (Fehler, 1982; Oelke et al., 2013). However, P-to-P reflection should have considerable strength as well. Part of the reason that PP waves are barely visible is the sensitivity of DAS. The sensitivity is linearly proportional to , projection of the wave vector on the fiber, and polarization of the body wave (Kuvshinov, 2015; Glubokovskikh et al., 2021). For the given incidence angles and seismic properties of the subsurface (Table 1), the receivers are almost three times more sensitive to PS waves than PP waves. For individual fractures, we may realistically have , so the reflected PP wave is six times weaker than the PS wave and may be well below the noise level.
Two other features of our data set are somewhat puzzling: negligible time delays of the direct P wave and the absence of a strong transmitted PS wave. The former one has been the main seismic characteristic of the SRV produced by hydraulic fracturing (Meek et al., 2019; Binder et al., 2020; Zhao et al., 2021). The P-wave slowdown is interpreted as the overpressure, which opens microcrack and fractures and expands pores in the formation (Binder et al., 2020). Binder et al. (2020) and Titov et al. (2021, 2022) report strong transmitted PS waves in their VSP data, which is caused by discrete low-velocity zones approximately 20 m thick. These zones, adjacent to major fractures, experienced stronger stiffness reduction than the rest of the SRV. Neither SOV-3 nor SOV-5 data contain strong . This leads us to a hypothesis that the fluid injections in well 5 pressurized natural fractures, which have very low communication with the surrounding rocks, causing extremely slow pressure relaxation.
To verify this hypothesis, we look at the evolution of apparent . Figure 4b shows the corridor stacks for each monitor vintage: traces computed by stacking the reflected PS waveforms along the traveltime curves picked for each event. We see that the reflected signals contain several wavelet phases, up to five for some events. This is likely caused by interference of signals scattered by multiple fractures. The amplitude trend for reflected wavelets remains stable for days, whereas the strain-rate anomalies and bottom-hole pressure transients occur in a matter of hours. In addition, we see that pressure gauges in a legacy vertical borehole between wells 3 and 5 record no pressure change precursors until certain abrupt events occur, and then the pressure remains constant (Figure 4c).
Therefore, our seismic analysis relies on the assumption of natural fractures isolated from the host rock. In the next sections, we detail how we quantify the parameters of the fracture set between MD 4400 and 4600 m, such as the number of fractures, their location, and their reflectivity (fracture compliance).
FRACTURE REFLECTIVITY AND LOCATIONS
Seismic scattering caused by fractures depends on spatial parameters, such as shape; orientation; position with respect to the monitoring fiber; and fracture “stiffness,” which is controlled by aperture, the rheology of the filling material, fluid pressure, and stress traction. A rigorous inversion of these parameters would require sophisticated modeling of wave propagation as well as a rock-physics model relating the fracture stiffness to coupled flow and mechanical processes. We lack prior information to constrain such models. Furthermore, DAS samples the wavefield along a single borehole and records only one component of the strain-rate tensor, which restricts the potential for a complete seismic inversion. Therefore, we implemented a more practical and robust approach, in which we sequentially refine the fracture models. This section treats fractures as infinite thin layers. This enables a relatively robust estimation of the fracture locations and reflectivity.
Rock-physics model of fracture reflectivity
For realistic cases, where the fracture has rough contacting surfaces and the stimulation fluid contains a high concentration of proppant grains, the compliance is difficult to estimate without a directly observed seismic response. However, for activated fractures, we may anticipate that an extremely large will have a much larger impact on the fracture reflectivity than a large but finite . A fracture filled by ideal fluid with bulk modulus has compliance parameters and . For a more realistic case of granular packing inside the fractures, we use the Hertz-Mindlin theory (Mavko et al., 2009), outlined in Appendix A. Such a model may capture the case of a proppant slurry inside hydraulic fractures or analogously a friable pressurized fault gauge material.
Reflection and transmission of a plane seismic wave illuminating a stack of parallel plane fractures may be computed using the matrix propagator method of Schoenberg and Protazio (2005). Figure 5a shows the angle dependence of and for three types of fracture filling material (see Table 1). We see that the reflection becomes stronger as the tangential component in the incident wave increases, due to a larger role of . Our field data correspond to an incidence angle of 60°. Among the three modeled scenarios, the field data resemble the case of a dense slurry: is of the order of 1%, which significantly exceeds and . In the long-wavelength limit, the scattering coefficients—, , and —are linearly proportional to a combination of angular frequency and tangential fracture compliance to . For realistic subsurface parameters, other terms in the power series are negligible. Figure 5b shows the frequency dependence of these parameters for a fracture filled by the dense slurry. All three of them change linearly with frequency as expected. This observation has two important implications for future analysis. First, the fracture (thin layer) acts as a differentiator filter. Second, the scattering strength scales linearly with fracture compliance . Thus, for activated fractures, we always observe a combined effect of an increased fracture aperture and reduced contact area or softer pore filling material.
Deconvolution of the fracture parameters
Now, we can apply the conceptual rock-physics model of the previous section to a quantitative analysis of the seismic data. Similar to the standard convolutional model of seismic trace (Mendel, 1990a), we can treat each trace in the stacked seismograms in Figure 3b as an interference of P-S wavelets reflected from the infinite parallel fractures. We detail the seismic deconvolution algorithm in Appendix B.
First, we need to estimate the wavelet. We stack direct P-wave arrivals along a lateral segment of the fiber between MD 4000 and 4900 m in the processed prestack seismic data for each monitoring vintage (Figure 3a). Note that in Figure 6a and 6b, the correlated waveform from SOV-3 contains a clear negative peak offset from the main one by 20 ms. The interference causes a clearly recognizable notch in the spectrum at a frequency of approximately 40 Hz (Figure 6c), similar to marine seismic data (Egorov et al., 2018). Most likely, the signal is a surface-related multiple produced by a strong near-surface reflection.
Second, we analyze the scattered wavelets. As shown in the previous subsection, a fracture acts as a differentiator filter on the scattered signals. Figure 6d and 6e compares the waveforms and amplitude spectra for the differentiated incident wavelet and stacked PS-wave reflection for the same seismic vintage (Figure 4b). The agreement is very good. We confirm again that the impact of the fracture pressurization on the matrix rock is minimal, and the reservoir changes are confined to the vicinity of the fractures, which means that the thickness of the perturbed zone is small compared to the fracture spacing. If not minimal, the propagating wave would experience significant scattering loss and attenuation within the fractured interval. Instead, the wavelet remains almost unchanged within the characterized interval, and the convolutional model should be valid.
Then we use the extracted wavelet along with a sparse-spikes deconvolution (Kormylo and Mendel, 1982) to determine the location of the fractures. This deconvolution refers to a group of inverse problem algorithms, which simultaneously estimate the location of each isolated reflector along with its reflectivity. Figure 7 summarizes the results of the deconvolution.
Figure 7 shows a diagram for the two snapshots of the fracture network. We determined that the scattered signal was largely produced by five fractures with spacing between 18 and 45 m. Note that the fracture locations are fixed to be the same among the vintages. For the dominant frequency of the PS wavelet (60 Hz), the wavelengths for P and S waves are 83 and 50 m, correspondingly. Thus, the reflections from fractures separated by approximately 20 m or less will effectively merge, whereas the more distant reflectors may be distinguishable. Figure 6d shows that the theoretical wavelet deviates from the observed wavelet at the very end of the waveform. The deconvolution algorithm picked up this feature by placing relatively strong reflectors closer to the toe. Appendix B provides more insights into the interference of the scattered events.
This section only considered the corridor stacks, meaning that the signals become stacked along the fiber length. Depending on the fracture height and lateral extent, the amplitude of PS reflections may change drastically along the fiber. Thus, the deconvolution results, including apparent reflectivity and fracture locations for abstract infinite fractures, are affected by the shape factor of the fractures. The next section examines this factor.
SEISMIC ESTIMATION OF FRACTURE EXTENT
Infinite 1D reflectors are just a first-order approximation of the seismic wave interaction with real natural fractures. When the fracture height is comparable to seismic wavelengths, approximately 100 m, diffraction effects play a major role (Born et al., 1999). Titov et al. (2021, 2022) use the amplitude distribution of the transmitted PS wave versus fiber length to estimate the height of the SRV. As mentioned previously, is weak in our data set. However, Figure 8 shows that the reflected PS waves feature a similar amplitude variation with distance. This section proposes a novel approach to efficiently model the seismic response from fractures. When combined with the results of the previous section, it provides robust estimates of the fracture location, their shapes, and reflectivity.
Scattering integral modeling
Numerical modeling of seismic wave interaction with fractures is known to be a formidable task for grid-based methods. To represent the fracture geometry and heterogeneous high-contrast properties, one would have to refine the grid size and decrease the time steps to an intractable limit even for modern supercomputers (sim, 2013). In addition, enforcing desired boundary conditions requires special treatment that requires larger models and increases the computational costs (Higdon et al., 2013). As a result, the fractures are often approximated by 2D flat rectangles of exaggerated thickness of approximately 1 m compared with a realistic thickness of approximately 1 mm (Wu et al., 2005; Titov et al., 2022). This approximation ignores potentially crucial effects associated with 3D wave propagation and fracture configuration, such as the proximity to lateral edges, irregular top and bottom edges, orientation relative to the fiber-optic cable, and variable fracture compliance.
We implement a modeling approach that uses the Born approximation (Snieder, 2002), and it is described in detail in Appendix C. This approximation is often called a small-contrast approximation, which may sound incorrect for fractures. However, the essential requirement relates to the magnitude of phase shift accumulated inside each inhomogeneity compared with the baseline wavefield (Hudson and Heritage, 1981). For thin low-reflectivity fractures (Figure 5), this requirement is fulfilled.
The Born approximation considers each element of fractures as a source of scattered wavefield excited by an incident wave that would exist in a baseline (without fracture) medium. Thus, we may add the contributions of individual elements to estimate the scattering from the entire fracture. Furthermore, we may switch on and off different types of scattering, so we focus on P-to-S conversion and target certain frequency bands by choosing an appropriate source function for the scattering elements.
For an incident plane wave, the scattering strength of each element is directly proportional to contrasts of density and stiffness inside the fracture. Thus, for a set of flat parallel fractures, we need to compute the scattered wavefield only once for an elementary fracture patch (Figure 9a). The total wavefield is then a linear combination of such patches with appropriate complex amplitudes: the modulus depends on the contrast between the fracture and background, and the phase depends on the phase shift of the incident plane wave at the location of the patch (Figure 9b).
Figure 10 shows the application of Born modeling for a P wave at 50 Hz and an incident angle of 60° on a 200 m high fracture and filled with a dense slurry (Table 1). The scattering is computed as interference of 10 fractures of 20 m height and 400 m width and 2 cm thickness, orthogonal to the monitoring fiber and the plane of incidence for the seismic signal. Figure 10b shows the strain amplitude distribution along the three vertical locations of the DAS receivers, as shown in Figure 10a. For the same incidence signal and fracture shape, the distance to the top of the fracture changes the observed trends drastically. Note how a bell-shaped curve for the DASs located at the top of the fracture transitions into a two-peak curve for the DASs at the bottom of the fracture.
Best-fit fracture shapes
Now, we can analyze prestack seismic data quantitatively based on the Born modeling approach. Ideally, we would iterate through the locations, orientations, and reflection strength of individual fracture patches until we match the scattered waveforms at different offsets from the fractures. Then we need to restrict the inversion to a simpler estimation problem. First, we sampled the axial strain only along a single borehole. Second, the data are noisy with relatively extensive intervals of signal-to-noise ratio that precludes interpretation, for instance, MD 4600–4700 m. We may only rely on a narrow frequency band of approximately 60 Hz, where the stacked wavelet (Figure 6e) has the strongest signal levels.
Thus, each monitoring vintage provides one curve of PS-wave amplitude along the fiber to estimate several parameters of multiple fractures. To constrain the search, we make the following five basic assumptions:
The fracture set consists of five fractures (the same locations in each vintage).
All fractures are rectangular.
The fracture width is sufficient to neglect the diffraction from the lateral edges (>70 m from fiber).
The orientation is orthogonal to the borehole.
Fracture compliance is constant within each fracture but may differ between fractures.
These assumptions are necessary if we would like to use the output of 1D deconvolution as an initial model for the search. Essentially, we refine the 1D model and augment it by locations of the top and bottom of each fracture. Effectively, the inversion is equivalent to a 2.5D model (infinite fracture length). However, the proposed framework is flexible and can incorporate irregular fractures and nonuniform fracture stiffness.
Figure 11 shows the results for three monitoring snapshots. For each of them, we performed a grid search over the fracture location along the fiber, vertical distance to the top and bottom, and fracture reflectivity in a frequency range between 58 and 65 Hz. Interestingly, the best-fit fracture shapes also predict weak transmitted PS wave, a puzzling field observation, as we mentioned in the preamble to this section. Figure 11d and 11f shows the predicted and observed amplitude curves at 25 Hz as a qualitative confirmation of the estimated fracture parameters. However, the noise level is too high to make any quantitative estimates based on these curves.
Fracture locations and reflectivity remained similar for the finite fractures compared with the infinite fractures in Table 2. The vertical extent of the fractures changed during the nine monitoring days. Initially, a majority of the fractures were located above the monitoring fiber, but they extended downward and intersected the borehole later (see diagrams in Figure 11d and 11f). It is not possible to estimate the depth of the fracture bottom below the fiber because recorded PS-wave amplitudes are almost insensitive to the signals produced by the parts of the fractures located under the fiber. Figures 9a and 10a show that the back-scattered wavefield is very weak when the fiber is located above the fracture.
DISCUSSION
Using the observed scattering events and other observations, we are able to explore the processes associated with fracture activation. First, we believe that the scattering events in the interval MD 4400–4600 m were largely caused by activation/relaxation of the same five natural fractures. Figure 9b shows the initial appearance of the fracture set. We think that no new fractures propagated to well 3 during the nine days of monitoring; otherwise, we would see clear fracture hits and induced seismicity associated with this propagation. Instead, these natural fractures produced scattering events prior to the stimulation of the corresponding stages in well 5; they were then activated and remained visible for five days after the conclusion of fluid injections (Figure 4b). The changes in reflectivity were caused by fracture dilation and subsequent closure along with some moderate growth. Next, we analyze some peculiar features in the evolution of fracture reflectivity.
Strain anomalies first extended to the interval of interest on days 3.5 and 3.9 (Figure 4b). The anomalies do not appear as fracture hits. Nevertheless, the fracture reflectivity abruptly started increasing (Figure 7). Most likely, the natural fracture set became hydraulically connected to a new cluster of hydraulic fractures created by stage 5. When the new fractures close at the end of injection, the natural fractures take a long time to release pressure. This is why the reflectivity trends are so smooth.
Then stages 6 and 7 injected much more fluid under higher overpressure, which boosted the reflectivity increase. At the same time, we detected a few microearthquakes northeast of well 3, where we expect to have fracture 2, which became the brightest fracture in the set up to this point. In addition, stage 7 produced high overpressure at a depth of 1970 m, far exceeding the fracture pressure. We believe that fracture 2 was strongly activated and grew laterally. However, we did not detect any fracture hits in the target interval, suggesting that the vertical extent did not change. Qualitatively, this conclusion agrees with the fact that the pressure gauge at 2010 m remained at ambient conditions. In addition, no microseismic events were detected close to well 5.
During the long waiting period between stage 7 and stage 8, all of the fracture reached their maximum activation state and started closing down. Then stage 9 led to an extremely extensive strain anomaly at 1000 m, which resulted in several fracture hits. Our Born-based interpretation suggests that fracture 1 intersected the monitoring fiber. The reflectivity of this fracture also changed the trend and stayed relatively high for the next three days. Stage 9 produced a strong strain anomaly, which included two interpreted fracture hits in the target interval. The Born modeling indicates that fracture 4 intersects the borehole and its reflectivity trend changed at the same time. In agreement with the seismic interpretation, the deeper pressure gauge, located at 2010 m, recorded noticeable overpressure.
We do not have reliable estimates of the extent of the natural fractures, except for the microseismic catalog. Microearthquakes can be triggered outside of the activated fractures as well as some parts of the faults and fractures may slip aseismically (Eaton, 2018). However, for the ACEFFL project, the vertical extent of the microseismic clouds (Figure 2c), approximately 100 m, matches the results of the analysis of seismic scattering events. In addition, the orientation of the fracture relative to the boreholes agrees with the borehole image logs (Figure 4) and the microseismic imaging by Ma et al. (2024).
Although the ACEFFL data set seems very comprehensive, we may not confidently estimate the parameters of the fractures and their permeability. Our joint interpretation involved a few assumptions that we must make to reconstruct the process of activation. Yet we may not determine with confidence the microscopic processes that reduced the stiffness of the fracture during activation: dilation, contact area, and removal of cohesive material. We are limited by the fact that fracture compliance and depends on stiffness of the fracture fill and its aperture. Thus, we may not distinguish various processes that may change these values: fracture opening, changes in the contacts between the fracture surfaces, vertical and lateral growth of the existing fractures, and formation of new fractures in-between the existing ones. Seismic resolution does not allow us to go beyond the assumption of constant and .
Potentially, an in-depth analysis of the strain anomalies and pressure measurements could provide additional insights. We are still unclear about the origin of the strain anomaly observed throughout the entire borehole on day 1. Also unclear is the fracture interaction with the monitoring fiber that manifested itself in extended strain anomalies concurrent with the changes in reflectivity trends. These anomalies creep slowly away from the stimulated interval.
Ideally, such analysis of the strain and pressure response would be combined with history-matching hydromechanical simulations. This may explain the characteristic temporal and spatial scales of the formation stimulation and relaxation. For now, we assume that the matrix permeability is extremely low, and connectivity to the hydraulic fractures disappears when the injection pressure drops to low values that explain the observed seismic features: negligible traveltime delays in the stimulated interval, the absence of the transmitted PS wave, and large relaxation times for the fracture reflectivity. Modeling of such aspects is beyond the scope of this study, which attempted to summarize the geophysical evidence supporting or opposing various subsurface scenarios.
CONCLUSION
This paper analyzed nine days of high-repeatability continuous borehole seismic data to reconstruct the process of reactivation of preexisting fractures in an unconventional reservoir, the ACEFFL unconventional reservoir. Our analysis relied on seismic scattering on fracture planes, more specifically, incident P-wave signals from permanent surface seismic sources to S waves recorded by horizontal DAS receivers in a monitoring well.
A seismic corridor stack for the P-S scattering events reveals how the fluid injections pressurized different natural fractures and how slowly the pressure dissipates between the injections. After the injections start, reflectivity increases for the fractures located two to three stages behind and three to four stages ahead of the hydraulic fracture corridor expected for the current stimulation stage. Typically, the intervals of increased reflectivity coincides with strain-rate anomalies recorded by low-frequency DASs. Such behavior strongly suggests that the scattering corresponds to reactivated natural fractures interacting with the monitoring fiber rather than new hydraulic fractures.
The scattered waveforms represent an interference of scattered signals for multiple fractures with spacing either smaller or comparable with the seismic resolution of approximately 25 m. To separate the contributions of individual fractures, we applied a deconvolution with fracture sparsity constraints. The deconvolution resulted in a set of five fractures spread over 120 m with slowly varying reflectivity for each individual fracture. One of the fractures exhibited significantly greater activation than the other fractures and dominated the seismic response over the nine day period. The pressure and microseismic data corroborate that this fracture grew laterally over that period of time.
Finally, we estimate the fracture height using amplitude variation of the reflected PS-waves along the fiber. We iterate through various configurations of the fracture network to find a fit to the observed amplitudes at two dominant frequencies of 60 and 25 Hz. The modeling algorithm implements a simple and fast integral formulation of single scattering on the fracture planes, the Born approximation. Although the 1D deconvolution assumes infinite plane fractures, it provides a robust initial model for the analysis of the prestack seismic waveforms; locations and reflectivity estimates changed only slightly.
We find that the fracture planes extended to 100–125 m above the monitoring well. Such a height agrees with the vertical extent of the microseismic cloud produced by the fracturing of the monitoring well later on. In addition, we find that the lateral edges of the fractures have negligible effect on the scattered events, which again supports the hypothesis of reactivation of extensive natural fractures. We do not find evidence of significant growth of the fractures, apart from some downward advancement of the fracture bottom. The fracture growth was probably driven by a gradual increase in the aperture of the preexisting fracture planes. Such a slow process may explain the atypical strain anomalies coinciding with the changes in apparent seismic reflectivity over time. Only rarely do these fracture interaction appear as fracture hits in the low-frequency DAS data. Previously, Tan et al. (2022) suggested that a bedding-plane slip may explain the complex spatial and temporal relationships between the low-frequency anomalies and seismicity. This might be a potential explanation for our observations that requires future investigations.
Overall, our study clearly shows the value that continuous active seismic monitoring using DASs can bring to understand the stimulation of unconventional reservoirs. The dimensions of the fracture set, height, and spacing are right at the edge of seismic resolution. The estimated fracture reflectivity, approximately 1%, is below the noise level. However, the high repeatability and high frequency of the seismic snapshots along with a meticulous analysis of the seismic scattered amplitudes enable high-precision tracking of the fracture reactivation. Without these subtle observations, it would not be possible to decipher the atypical strain anomalies that extend beyond the stimulated intervals as well as the patchy distribution of induced seismicity.
ACKNOWLEDGMENTS
This project was funded by the U.S. Department of Energy, assistant secretary for Fossil Energy and Carbon Management, Office of Fossil Energy and Carbon Management, Resource Sustainability Program through the National Energy Technology Laboratory under contract no. DE-AC02-05CH11231. We are grateful to SM Energy for providing the distributed acoustic sensor data and permission to publish this work. We thank D. Hill and D. Zhu for their guidance in this project. We are very grateful to R. Bohn for his insights and guidance regarding subsurface operations and engineering data. We thank Silixa for assistance in data acquisition.
DATA AND MATERIALS AVAILABILITY
Data associated with this research are confidential and cannot be released.
APPENDIX A HERTZ-MINDLIN THEORY
Figure A-1 shows how the fracture stiffness changes with more proppant getting into the fracture as well as relaxation of the fluid pressure inside the fracture. We see that when the normal stress at the grain contacts disappears, the shear modulus becomes zero, which corresponds to the Wood’s model, equation A-1. In contrast, when the fluid pressure drops below the fracturing pressure, the compacted proppant pack becomes stiff enough to effectively seal thin fractures, so seismic scattering is very weak.
APPENDIX B SPARSE-SPIKE DECONVOLUTION
The objective function in equation B-3 leads to a strongly nonlinear optimization to determine the locations and a relatively stable subsequent determination of (Velis, 2008). Fortunately, the dimensionality of our problem is low: and . Thus, the optimization does not require tedious development of computationally efficient methods. We used an optimization module of Wolfram Mathematica (Wolfram Research, 2008b). As part of the optimization process, the module selects an optimal strategy for the problem, which resulted in a simulated annealing for the location of the fracture and subsequent gradient-based search for the optimal , similar to the methodology used by Velis (2008), although the objective functions were slightly different.
The deconvolution determined fracture reflectivity for each seismic vintage independently, but must remain the same. This is because we assume that the fluid injections in well 5 did not propagate new hydraulic fractures that intersect well 3. Figure B-1 shows that the seismic likelihood and spacing constraint become flat when the number of fracture is equal to five. In addition, seismic likelihood suggests that we need to have at least three fractures to achieve a reasonable fit between the inverted and observed seismic reflectivity. Figure B-2 shows the agreement between the observed and estimated reflectivities. To facilitate the comparison, we added a model curve that accounts for the limited bandwidth of the seismic data, SOVs operate in the frequency range of 8–90 Hz.
APPENDIX C COMPUTATION OF THE BORN INTEGRAL
We used a fracture patch of 25 m × 25 m × 0.004 m. We computed a 3D distribution of the scattered PS wavefield only at the origin of the coordinate system. According to equation C-3, to model a patch at a different location , we only need to multiply the computed field by the account for the phase difference of the incident P wave between the origin and specified location.
For numerical integration of equation C-1, we used Wolfram Mathematica (Wolfram Research, 2008a). This module implements adaptive numerical integration strategies that account for the rate of variation of the integrand and symmetries of the integration domain.
Biographies and photographs of authors are not available.