Do Simple Analytical Models Capture Complex Fractured Bedrock Hydraulics? Oscillatory Flow Tests Suggest Not

Fractured sedimentary bedrock aquifers represent complex flow systems that may contain fast, fracture‐dominated flow paths and slower, porous media‐dominated flow paths. Thus, characterizing the dynamics of flow and transport through these aquifers remains a fundamental hydrogeologic challenge. Recent studies have demonstrated the utility of a novel hydraulic testing approach, oscillatory flow testing, in field settings to characterize single bedrock fractures embedded in low‐porosity sedimentary bedrock. These studies employed an idealized analytical model assuming Darcian flow through a nondeforming, constant‐aperture, nonleaky fracture for data interpretation, and reported period‐dependent effective fracture flow parameters. Here, we present the application of oscillatory flow testing across a range of frequencies and inter‐well spacings on a fracture embedded in poorly cemented sedimentary bedrock with considerable primary porosity at the Field Site for Research in Fractured Sedimentary Rock. Consistent with previous studies, we show an apparent period‐dependence in returned flow parameters, with hydraulic diffusivity decreasing and storativity increasing with increasing oscillation period, when assuming an idealized fracture conceptual model. We present simple analyses that examine non‐Darcian flow and borehole storage effects as potential test design artifacts and a simple analytical model that examines fluid leakage to the surrounding host rock as a potential hydraulic mechanism that might contribute to the period‐dependent flow parameters. These analyses represent a range of conceptual assumptions about fracture behavior during hydraulic testing, none of which account for the measured responses during oscillatory flow testing, leading us to argue that other hydraulic processes (e.g., aperture heterogeneity and/or fracture hydromechanics) are necessary to accurately represent pressure propagation through fractured sedimentary bedrock.


Introduction
Owing to their areal coverage and significant thickness, sedimentary bedrock aquifers contain significant storage volumes that represent ideal targets for wastewater injection and carbon capture and storage (CCS).Though they represent a small volume of sedimentary bedrock aquifers, secondary porosity features (i.e., fractures) provide conduits that dominate flow and transport throughout these aquifer systems (Viswanathan et al. 2022), especially in deep sedimentary bedrock systems where subhorizontal bedding plane fractures may represent the primary conduit for pressure propagation and fluid flow.Seen through the lens of a carbon-constrained world, developing a better understanding of the flow processes occurring within these fractured sedimentary bedrock aquifers remains critical for reducing atmospheric greenhouse gas concentrations through alternative energy development and CCS (Viswanathan et al. 2022).
Due to the presence of multiscale heterogeneity, characterizing the physical properties that govern groundwater flow through fractured bedrock-specifically transmissivity and storativity-remains a fundamental hydrogeologic challenge and active area of research that has generated a significant body of literature describing multiple hydraulic characterization approaches (e.g., Butler 2005;Cardiff et al. 2012;Illman 2014).Pressure-based characterization approaches, where pressure responses are measured at discrete points in response to introduced hydraulic stresses, are by far the most common characterization method and provide the most direct local information on fracture flow properties.Due to minimal equipment requirements and ease of implementation, the single-well slug test and cross-well constant-rate pumping test-which tests fractured bedrock at distinctly different scales-are the most prominent pressure-based characterization approaches.Characterizing fractured bedrock with constant-rate pumping tests can involve the extraction of significant volumes of water from the subsurface, which risks dewatering wells due to low fracture storativity and has the potential to significantly alter the ambient stress and flow fields (Schweisinger et al. 2011;Cardiff et al. 2018).Furthermore, pressures measured under constant-rate pumping may incorporate sensor drift and/or hydrologic noise (e.g., recharge events, evapotranspiration, and barometric loading) into the recorded signal, complicating analysis.Single-well slug tests are smaller-scale and faster hydraulic tests that are less likely to significantly alter hydraulic stresses along the fracture and less likely to incorporate sensor drift and/or hydrologic noise into the recorded responses.However, slug tests are sensitive primarily to hydraulic properties immediately adjacent to the well, making them susceptible to near-wellbore skin effects that may occur as a result of drilling-induced formation damage or inadequate well completion activities (Butler Jr and Healey 1998;Butler 2005).
Oscillatory flow interference testing-where water is alternatively pumped from and injected into a borehole in a periodic manner-is an actively evolving pressurebased hydraulic characterization strategy designed to overcome the limitations with the more common hydraulic characterization approaches discussed above.Oscillatory flow interference testing provides many benefits-which are well documented in the literature-over other pressure-based characterization approaches (Fokker and Verga 2011;Cardiff et al. 2013;Bakhos et al. 2014;Guiltinan and Becker 2015;Rabinovich et al. 2015).For example, oscillatory flow tests are capable of characterizing fractured bedrock at multiple scales, with lowfrequency pumping signals testing fractures at scales similar to a constant-rate pumping test and high-frequency pumping signals providing near-field hydraulic information similar to the single-well slug test (Cardiff et al. 2013).Furthermore, oscillatory flow testing emphasizes the effects of fracture storativity at the scale of interest because recorded pressure responses remain in a state of constant transience (i.e., they do not achieve the steady-shape condition associated with constant-rate pumping tests) during hydraulic testing (Guiltinan and Becker 2015).This can eliminate the need for earlytime and late-time drawdown analysis, which can be impacted by near-borehole conditions and hydraulic conditions beyond the well field under constant-rate pumping conditions, respectively.Also, oscillatory flow tests can be designed to yield a net-zero drawdown resulting in minimal alterations to the ambient flow and stress fields, minimizing recovery time between individual tests and allowing for nearly instant follow-on testing.Finally, using standard signal processing techniques the introduced pressure signal is readily separated from instrument drift, instrument noise, and/or other nontesting-related hydrologic signals that commonly contaminate collected hydrologic data (Bakhos et al. 2014).
Previous studies highlight the utility of natural and introduced oscillatory pressure signals to characterize alluvial aquifers (Ferris 1952;Rabinovich et al. 2015;Sobolevskaia et al. 2021), unconsolidated sedimentary aquifers (Rasmussen et al. 2003), and petroleum reservoirs (Johnson et al. 1966;Kuo 1972;Hollaender et al. 2002).Here, we focus on studies that use oscillatory pressure signals to characterize effective flow properties in fractured sedimentary bedrock.These previous field experiments share two common characteristics: (1) they are conducted in well-cemented sedimentary bedrock with less than 5% primary porosity (Maineult et al. 2008;Guiltinan and Becker 2015); and (2) they utilized an idealized analytical model that conceptualizes the tested fracture as a nondeforming, parallel-plate fracture bounded by impermeable host rock, while neglecting testing non-idealities such as borehole storage effects.Using a simple analytical model that relates the amplitude ratio (i.e., attenuation) and phase delay between the observation and stimulation wells to fracture flow properties, oscillatory pressure signals have been used to estimate the effective hydraulic diffusivity of isolated bedrock fractures (Renner and Messar 2006;Guiltinan and Becker 2015;Sayler et al. 2018), infer inter-well connectivity (Becker and Guiltinan 2010;Guiltinan and Becker 2015), and investigate the flow geometry of a complex fracture network (Sayler et al. 2018).Assuming an idealized fracture during analysis, the above-mentioned studies all noted a decrease in hydraulic diffusivity estimates with increasing oscillation period.That is, the estimated flow parameters appear to be dependent on the period of the pumping signal.The observed decrease in diffusivity appears to be associated with increasing storativity at longer pumping periods as opposed to major changes in effective transmissivity (Renner and Messar 2006;Guiltinan and Becker 2015).
The observed period-dependence-when using these idealized analytical models to characterize fractured bedrock-suggests that unmodeled flow conditions occurring within the fracture may be altering observed pressure responses.For example, testing "non-idealities" such as non-Darcian flow and/or borehole storage effects and unmodeled hydraulic conditions such as fluid exchange (i.e., leakage) with the surrounding host rock may each impact measured oscillatory pressure signals at different oscillation periods (Figure 1).The impact of each of these non-idealities can be individually tested with modified analytical or simple numerical modeling approximations.Exploring other unmodeled hydraulic processes-such as the impacts of multiscale heterogeneous flow and storage occurring within the fracture, as suggested by Cardiff et al. (2013), or the fracture strain response to applied hydraulic stresses (i.e., fracture hydromechanical effects), as hypothesized by Guiltinan and Becker (2015)-requires more detailed numerical and/or stochastic modeling strategies, which are beyond the scope of this work.
To date, an accepted, and yet untested, hypothesis states that flow through the fracture plane is concentrated along a highly conductive backbone and at longer pumping periods water interacts with flow restricted areas within the fracture that do not contribute to flow, leading to an apparent increase in estimated storativity and period-dependent effective flow parameters (Renner and Messar 2006;Guiltinan and Becker 2015).In the most general sense, this hypothesis states that highpermeability channels are bounded by low-permeability regions of fluid storage that do not contribute to flow, which can be a result not only of heterogeneous storage within the fracture, but also fluid exchange between the fracture and host rock.This highlights the complexity of flow processes occurring within an individual fracture and motivates the importance of choosing the correct conceptual model when characterizing flow properties in fractured sedimentary bedrock.
In this work, we present the application and analysis of oscillatory flow testing to an isolated fracture embedded within poorly cemented sedimentary bedrock, which to our knowledge represents the first analysis of its kind.We describe oscillatory flow tests conducted across a range of oscillation periods, oscillation amplitudes, and inter-well spacings at a field site intentionally designed to investigate issues such as potential flow anisotropy and scale effects under oscillatory flow conditions.Utilizing this rich dataset, we investigated the presence of period-dependent effective fracture flow parameters when using an analytical model that assumes a nondeforming, parallel-plate fracture bounded by impermeable host rock (henceforth referred to as an ideal fracture), consistent with previous studies (Renner and Messar 2006;Guiltinan and Becker 2015;Sayler et al. 2018).Furthermore, using simple analytical modeling approaches we explore test design artifacts (such as non-Darcian flow and borehole storage effects) and unmodeled hydraulic processes (such as fracture leakage) as potential mechanisms that might contribute to the apparent period-dependence in estimated effective fracture flow parameters (Figure 1).Wellbore skin has been shown to not effect observation signals collected during oscillatory flow testing (Ahn and Horne 2011); therefore, we do not consider these wellbore skin effects as a potential mechanism that might contribute to the observed period-dependence.Finally, in addition to performing our analyses, this work provides collected data, data processing routines, and code used for analysis as Supporting Information so that other researchers may investigate this rich dataset and utilize our approaches for aquifer characterization efforts that employ oscillatory flow testing.

Methodology Field Site
The Field Site for Research in Fractured Sedimentary Rock ([FSR] 2 )-an easily accessible research site located approximately 24 km northeast of Madison, Wisconsin, USA (Figure 2A)-was constructed in 2017.The (FSR) 2 lies within the Elmer and Edna Culver Springs Conservancy, which is centered around a system of natural springs fed by groundwater discharge through streambed sediments and the surface expression of bedrock fractures, forming the headwaters of Token Creek (Bahr and Parent 2001).The (FSR) 2 contains four 0.15 m diameter vertical boreholes installed to a depth of 122 m below land surface (bls), and one 0.15 m diameter vertical borehole installed to a depth of 62 m bls.All boreholes were steel cased to a depth of 10 m bls and left open to the formation below that depth.We arranged the boreholes in a variable-distance 5-spot configuration, with inter-well spacings from 4 to 16 m (Figure 2A), allowing us to investigate spatial heterogeneity across these scales.Additionally, we oriented two borehole pairs (A1 to B2 and B1-B4)-both with an inter-well spacing of 10 m-in an orthogonal manner to allow the investigation of potential fracture flow anisotropy (Figure 2A).
The bedrock geology at the (FSR) 2 consists of Cambrian-aged sandstone, which underlies a layer of unconsolidated glacial till approximately 4 m in thickness (Figure 2B).The sandstone bedrock contains primarily horizontal bedding with a minor regional dip of approximately 2 to 3 m/km to the East (Cline 1965;Bahr and Parent 2001;Gellasch et al. 2013).Locally absent at our field site, the Eau Claire shale forms a regional aquitard that separates the Wonewoc and Mt.Simon Formations into shallow and deep aquifer systems throughout south-central Wisconsin (Bahr and Parent 2001;Bradbury et al. 2013;Parsen et al. 2016).The absence of this confining unit at the (FSR) 2 makes the Wonewoc to Mt. Simon transition indistinguishable, leading us to consider the sandstone bedrock a single hydrostratigraphic unit with similar hydraulic properties.Core analysis of samples throughout northern Illinois and south-central Wisconsin in the geologic formations that comprise the sandstone bedrock at the (FSR) 2 indicate effective primary porosities ranging from 11% to 20% (Wisconsin Geological and Natural History Survey [WGNHS] 2015), which is approximately one order of magnitude larger than fractured sedimentary bedrock units previously characterized using oscillatory flow testing (Renner and Messar 2006;Guiltinan and Becker 2015).The Mount Simon Formation, at the base of our well installations, represents the same geologic unit being considered for large-scale CCS operations in Illinois, where more saline groundwater is found (Frailey et al. 2011).Locally, the Mount Simon formation serves as the main source of drinking water for Madison owing to its high transmissivity and high water quality (Gellasch et al. 2013).
Geophysical logging at the (FSR) 2 showed multiple laterally extensive, flat-lying bedrock fractures that appear to follow bedding planes, which is consistent with regional geologic descriptions (e.g., Parsen et al. 2016).The shallowest bedrock fracture occurs at 13 m depth, which we believe to be the surface exposed fracture contributing groundwater discharge to the nearby Token Creek springs.Identified fractures extend to a depth of approximately 45 m bls, with larger fracture spacing at increasing depths (Figure 2B).Geophysical logging showed no evidence of high-angle or vertical fractures throughout the boreholes.In this work, we focus our discussion on a fracture at 32 m depth, shown as the dashed red line in Figure 2B.

Oscillatory Flow Interference Testing
In this section, we describe a 2019 field campaign where we conducted 149 oscillatory flow interference tests using oscillation periods across 1.5 orders of magnitude-from 10 to 480 s-and stimulation amplitudes across more than three orders of magnitude-from 0.005 to 11 m-to characterize the fractured sedimentary bedrock at the (FSR) 2 .Throughout the field campaign we tested a single hydraulically isolated fracture, and we also explored our ability to propagate pressure signals from a hydraulically isolated fracture to the surrounding host rock.For the purposes of this analysis, we focus our discussion on a subset of 85 oscillatory flow tests conducted in the identified fracture at 32 m bls.
We installed inflatable packer strings into each borehole, where a packer string consisted of two individual packers connected by a perforated stainless-steel pipe 0.9 m in length, to hydraulically isolate the tested fracture.The packer strings were lowered into the borehole using 0.05 m diameter steel pipe so that the perforated interval was centered at a depth of 32 m bls.Each individual packer was then inflated to 700 kPa above the ambient hydrostatic pressure to ensure an adequate seal along the borehole wall and limit fluid leaks around the packers.With the packers fully inflated, the static water level within the riser pipe was approximately 8 m bls in all wells, providing a water column 24 m in length that could be forced into/out of the fracture using the pneumatic testing approach described below.
To conduct a single oscillatory flow interference test, we used an in-house developed, software-controlled, pneumatic well-head apparatus (the "Osc-a-Lot") to generate flow through the tested fracture.We connected Osc-a-Lot to the riser pipe at the stimulation well, creating an air-tight seal during testing.We passed two fiber optic transducers in air-tight cable glands through Osca-Lot, into the stimulation well riser pipe.One transducer extended to the center of the perforated interval at 32 m depth, which recorded pressure changes due to both the height of the water column and air pressure NGWA.org above the transducer.The second transducer was placed above the water column to record air pressures within the riser pipe.This transducer configuration allowed us to precisely determine water level changes within the riser pipe during testing by removing the effects of air pressure applied to the water column from the recorded pressure signals.To monitor observation head changes, we placed fiber optic transducers (Roctest FOP-MicroPZ Piezometer) approximately two meters below the water column in the riser pipe-which is hydraulically connected to the perforated interval-at each observation well (Figure 1).We selected transducers with a calibration range of 350 kPa and a 125 Hz native sampling rate, which provides sub-cm measurement resolution while allowing us to investigate a broad range of stimulation amplitudes.
The Osc-a-Lot apparatus consists of an inlet valve connected to a gas compressor, and an outlet valve that vents to the atmosphere; it is otherwise airtight.During testing, pressure command signals from a computer were sent to a pressure control valve that opened the proper valves to match the requested pressure-i.e., the inlet valve was opened if pressures were too low, or the outlet valve was opened if pressures were too high.The computer specified temporal variations in pressure such that a sinusoidal pressure signal was produced at the stimulation well with prescribed amplitude and oscillation period.Meanwhile, the pressure change is recorded at the remaining four observation wells.
Prior to testing, Osc-a-Lot applied a prescribed air pressure into the stimulation well to depress and "charge" the water column, which facilitates water level increases during the "pumping" portion of the oscillation cycle.We started each oscillatory flow test by initializing the data acquisition system and utilizing the live visualization capabilities to monitor pressure stabilization in each well.Our analysis approach requires we only know pressure changes from an established baseline; therefore, once we observed pressure signal stability in all wells, we reset all pressure offsets to an initial head change value of zero.
Following a brief period of quiescence of approximately 5 s, Osc-a-Lot began cycling fluid through the tested fracture by adding and removing air from the riser pipe to generate the prescribed sinusoidal pressure signal.Individual test files for each well, with the associated transducer scaling data, are automatically saved when signal recording begins.Using the water level changes recorded in the stimulation well we calculated the pumping rate as: where r pipe (L) is the radius of the riser pipe and dh/dt (L/T) is the rate of change of the water level in the riser pipe.
Figure 3 shows the recorded head signals and the calculated flow rate using Equation 1 for a single oscillatory flow test with an oscillation period of 20 s, and the pressure stimulation at well A1.The head change amplitude at the stimulation well is 0.33 m and the peak flow rate is 2.1 × 10 −4 m 3 /s (Figure 3).We highlight that flow rates of 0 m 3 /s correspond with the positive and negative peak head amplitudes in the stimulation well (Figure 3).We also note that the amplitude of the observation signals decrease, and timing of the amplitude arrival (i.e., phase delay) increases with increasing distance from the pumping well (Figure 3).
The range of oscillation periods available for testing was constrained on the low end by the software controlling Osc-A-Lot and on the high end by the length of the water column in the riser pipe-i.e., the volume of water available for cycling (Sayler et al. 2018).We discovered that Osc-A-Lot could not produce the desired input signal for oscillation periods shorter than 10 s or longer than 480 s during our field experiments.With oscillation periods less than 10 s, the pressure control valve on Osc-a-Lot did not release air fast enough to match the prescribed signal.With oscillation periods longer than 480 s the low volumetric flow rates resulted in pressure equilibration along the fracture and little to no detectable signals at observation wells.The constraint on longer oscillation periods could be overcome by having a larger volume of water available for cycling that would allow larger amplitude pressure stimulations, either by testing deeper fractures to increase the length of the water column or using a dual-pump system similar to Rasmussen et al. (2003), as opposed to the pneumatic approach presented in this work.We assessed the robustness and repeatability of oscillatory flow testing as a fracture characterization technique by rotating Osc-A-Lot across all five wells, generating multiple reciprocal tests (Figure 2A).We implemented test design strategies that seek to minimize uncertainty in estimated flow parameters.First, we recorded pressure signals at a high sampling frequency of 125 Hz, which has a significant effect on reducing parameter uncertainty by allowing us to more accurately identify amplitude arrival and phase delay at the observation wells (Patterson and Cardiff 2022).Furthermore, we collected a time series that is at least 10 full cycles in length for each oscillatory flow test, which generated at least five full cycles available for analysis and avoided any transience associated with testing "ramp-up" or "ramp-down."This approach balanced the trade-off in minimizing uncertainty in the estimated parameters and duration of the field experiments (Bakhos et al. 2014;Patterson and Cardiff 2022).

Field Data Analysis
In this section we describe our largely automated data processing and analysis workflow used to analyze the head change signals collected during an individual oscillatory flow test to generate estimated fracture flow parameters with uncertainty.To maintain clarity throughout the manuscript we adopt the language presented by Patterson and Cardiff (2022), where "observation signal" refers to the timeseries of recorded head changes with noise at an observation well, "data" refers to the extracted Fourier coefficients estimated through linear signal processing techniques, "error" refers to the uncertainty in the estimated Fourier coefficients as a result of noise in the observational signal, and "uncertainty" refers to the 95% confidence interval associated with fracture flow parameter estimates determined through linearized error propagation.We applied the gradient-based inversion approach described below to optimize parameters associated with all conceptual models described later: the ideal fracture analysis, the leaky fracture analysis, and borehole storage-impacted analysis.
Initial processing of recorded head change observation signals proceeded in the following manner.The original text files are loaded and converted into MAT-LAB binary files that we use to manually visualize each pressure signal and eliminate any observation signals that show unexpected jumps in head change data.While manually visualizing the individual signals, we identify the steady-periodic (i.e., constant amplitude and phase) portion of the signal to be trimmed for further analysis.We automatically excluded the first two cycles and the last two cycles of the observation signal to reduce the potential impacts of transient effects from test initiation or termination on the recorded signal.
Following initial data processing and trimming, we resample the observation signal to a constant time step of 0.008 s to remove any inconsistencies related to the sampling rate of the data acquisition system.Under steady-periodic conditions, each recorded observation signal can be represented as a linear combination of sinusoids with their associated noise (Figure 4): where P (T) is the stimulation period, r and i (L) are the real and imaginary Fourier coefficients, respectively, and ε (L) is the measured noise.Note that these Fourier coefficients are directly related to amplitude (A) and angular phase (ϕ) of measured signals as: where we define artan() as the "four-quadrant" arctangent with two inputs, as implemented in several computational languages.
Taking advantage of the linearity of Equation 2, we can write this as a matrix system of equations (Equation 5), with the Fourier coefficients (Φ) representing the only unknown.We take a least-squares approach to relate the noisy observation signal to the estimated Fourier coefficients (Equation 6), which serves as the data used during inversion (Bakhos et al. 2014;Patterson and Cardiff 2022).Figure 4 highlights that this workflow can reliably extract the noise-free signal from observation signals with amplitudes as small as 5 mm with 0.2 mm amplitude noise (i.e., low signal-to-noise ratio).
We used the mean squared difference between the observation signal and the reconstructed noise-free observation signal (Figure 4) as a measure of the observation signal measurement noise σ 2 .Using the estimated signal measurement noise, we calculated the data error covariance matrix (R) given by Equation 7, which quantifies the error associated with our estimated Fourier coefficients (i.e., data error).For our field experiments, R is an n × n-where n is the number of Fourier coefficients-diagonal matrix with zeros for the off-diagonal components, indicating independent data errors.In the case of multi-frequency inversion, R is a block-diagonal matrix, where the main diagonal is the data error covariance for each individual frequency component.
After extracting the Fourier coefficients and calculating the data error covariance matrix, we implemented a numerical inversion approach to estimate the fracture flow properties and their associated uncertainty.Specifically, we applied a Levenberg-Marquardt inversion algorithm under a Bayesian framework-as described by Patterson and Cardiff (2022)-to minimize the objective function: where Φ is the n × 1 vector of extracted Fourier coefficients and h(s) represents the chosen forward model that takes fracture flow parameters as inputs and returns Fourier coefficients as outputs.Minimizing Equation 8is equivalent to maximizing the likelihood of fracture flow parameters given the extracted Fourier coefficients for the observation signal (Aster et al. 2018).For our analysis, we selected an uninformative prior-that is all parameter sets are equally likely-and assuming the estimated parameters follow a log-normal distribution, the mean value of the posterior distribution represents the maximum likelihood of the posterior distribution.We imposed a nonnegativity constraint during inversion by perturbing log-transformed parameters.
Our inversion algorithm iteratively moves toward the optimal parameters by numerically approximating the Jacobian, calculating the objective function gradient at each iteration, and taking steps toward the minimum objective function value.We declare convergence at the optimal values when relative change in parameters and relative change in objective function value is less than or equal to 1e−6 between consecutive iterations.This numerical inversion approach differs from the analytical inversion algorithm developed by Rasmussen et al. (2003); however, our methodology can help overcome challenges with parameter nonuniqueness that has been shown to occur in the presence of phasewrapping (Cardiff and Sayler 2016), that is, observation signal phase delays greater than 2π .
Last, assuming local linearity at the optimal parameters, we used linearized error propagation to estimate the parameter uncertainty.We calculated the parameter covariance matrix (C), given by Equation 9: where J is the m × n numerically approximated Jacobian (i.e., sensitivity) matrix at the optimal flow parameters (s * ) and m is the length of the parameter vector.We extracted the diagonal elements of the parameter covariance matrix, which represents the variance for each parameter and calculated the 95% confidence interval for the optimal fracture flow parameters as s * ± 1.96 √ tr[C(s * )] following Aster et al. (2018).

Exploring Fracture Flow Conceptual Models
In this section, we explore a range of conceptual models used to model flow through fractured bedrock and their potential impacts on recorded observation signals and the resulting fracture flow parameter estimates under oscillatory flow conditions.First, we replicate the analysis approach seen in previous studies assuming an ideal fracture to assess for the presence of period-dependent fracture flow parameters at the (FSR) 2 .Then, we employ a range of alternative conceptual models to analyze and understand the impacts of these conceptual models on observation signals collected during our field experiments.

Ideal Fracture Analysis
First, we estimate the effective flow properties for the tested fracture assuming it behaves like an ideal fracture, which is consistent with previous approaches to estimate effective fracture flow properties under oscillatory flow conditions (Renner and Messar 2006;Guiltinan and Becker 2015;Sayler et al. 2018).Under this ideal fracture conceptualization, we assumed the tested fracture can be modeled as a fully confined aquifer and invoked assumptions identical to those of the Theis solution, with the exception of an oscillatory flow rate as opposed to a constant pumping rate.More specifically, we assume a rigid, parallel-plate fracture of infinite areal extent bounded by impermeable fracture walls.Furthermore, we assumed steady-state hydraulic conditions prior to testing with zero head change at infinite distance.Under these assumptions, we selected the analytical model described in Rasmussen et al. (2003)-originally developed by Black and Kipp (1981)-as our forward model which takes fracture flow parameter inputs and outputs the real and imaginary Fourier coefficients: where Q 0 (L 3 /T) is the peak pumping rate, T (L 2 /T) is transmissivity, D (L 2 /T) is hydraulic diffusivity, r (L) is the radial distance between the pumping and observation well, K 0 is the zero-order modified Bessel function of the second kind, i is the imaginary number (i.e., √ −1), and ω (T −1 ) is the angular frequency (2π/P ), where P (T) is the stimulation period.This expression differs from Black and Kipp (1981) in that it applied a flux boundary condition at the pumping well, allowing the separation of hydraulic diffusivity into its component parts, transmissivity and storativity (Rasmussen et al. 2003).We employed the gradient-based inversion approach described above in conjunction with the described forward model to invert for transmissivity (T ) and storativity (S) of the tested fracture.
Using this approach, we generated a set of effective flow parameters for each observation signal, yielding four sets of optimal flow parameters for a single oscillatory flow test, one for each observation well (Figure 2A).As a quality control measure, we first evaluated the parameter estimates for reciprocal and repeated flow tests, and we observe that parameter estimates for a specified stimulation period and radial distance cluster around a consistent value (Figure 5).This result highlights the repeatability and provides confidence in the returned fracture flow parameter estimates.
As seen in previous fracture characterization efforts, our results show a period-dependence in the fracture flow parameter estimates that follow a logarithmic trend.Specifically, we observe that diffusivity estimates decrease and that storativity estimates increase with increasing oscillation period (Figure 5), which is consistent with previous studies (Renner and Messar 2006;Guiltinan and Becker 2015;Sayler et al. 2018).There is no obvious period-dependent trend in transmissivity estimates; however, we note a large variability in transmissivity at the shortest oscillation periods that converges to a consistent value at the longest oscillation periods (Figure 5B).Previously unreported in the literature, our analysis also showed correlations between estimated fracture flow parameters and inter-well spacing.Specifically, for a given oscillation period we note higher diffusivity and lower storativity values at the longest inter-well spacings (Figure 5).Additionally, our analysis indicates a lack of flow anisotropy when comparing the parameter estimates of the well pairs with an orthogonal orientation at equal distance (Figure S1).
Figure 6 shows the recorded observation signals with the best fit observation signal for four representative oscillatory flow tests at a single observation well.The best fit signal was reconstructed using the returned optimal parameters at multiple oscillation periods and a constant inter-well spacing of 8 m.We highlight the presence of a secondary frequency component in the recorded signal at 120 and 300 s oscillation periods that leads to reduced model fit to the observation signal at longer oscillation periods (Figure 6).These secondary frequency components result from the Osc-a-Lot system adjusting pressures within the riser pipe mid-cycle to match the sinusoid prescribed in the controlling software.We could improve the fit to these observation signals by considering the secondary frequency components when extracting the signal Fourier coefficients and during inversion; however, we opted to focus our characterization efforts on the single frequency analysis using the prescribed oscillation period for each test.Patterson and Cardiff (2022) demonstrated that if the selected forward model used for inversion captures the processes present in the tested system during oscillatory flow testing, then a multi-frequency inversion approach can return a single set of flow parameters (i.e., consistent T and S) that fits the observed data within the expected data measurement error.Here, we test the hypothesis that we cannot adequately fit observed data using a homogeneous analytical model (Equation 10).We used a multi-frequency approach-described in Patterson and Cardiff (2022)-to simultaneously invert multiple oscillatory flow tests for a single T and S that provide the best fit.This multi-frequency analysis examines our assumptions that the tested fracture can be represented as two smooth walls with a constant separation distance bounded by impermeable host rock.To reject our hypothesis that the observed data cannot be adequately fit with a homogeneous model, we would expect the root mean square error (RMSE) between the modeled and observed data to be the same order of magnitude as the estimated data error given in Table 1, which would indicate that any modeled data misfit due to hydraulic processes cannot be separated from misfit due to data measurement error.
First, we simultaneously inverted all oscillatory flow tests at the largest inter-well spacing of 16 m and noted significant misfit (>1 cm) for many of the data points   (Figure 7A).Furthermore, our results show that the modeled data misfit-given by the RMSE-is three orders of magnitude larger than the misfit expected due to data measurement error (Table 1).Next, we simultaneously inverted all oscillatory flow tests with a stimulation period of 300 s, and we note that the modeled data misfit provided by the returned flow parameters are generally overestimated (Figure 7B).Consistent with the previous results, this analysis shows that the RMSE is three orders of magnitude larger than the misfit expected due to data measurement error (Table 1).

Non-Darcian Flow Analysis
Non-Darcian flows are a common concern when conducting hydraulic tests in fractured bedrock because of the increased fluid velocities occurring within fractures compared to porous media.We applied the smallest possible stimulation amplitudes needed to propagate measurable signals to the surrounding observation wells to minimize the potential of non-Darcian flow regimes during our field experiments.In this section, we test our assumption that Darcy's Law is valid during our field experiments by determining the Reynolds numbers for each oscillatory flow test and by assessing for the expected positive linear correlation between observation head amplitude and peak volumetric flow rate at a constant oscillation period.
First, we calculated the Reynolds number (Equation 11) at the pumping well-where we expect fluid velocities to be the greatest-for each oscillatory flow test to determine if the flow regime within the fracture is dominated by viscous or inertial forces during testing: where v (L/T) is flow velocity through the fracture, which we calculated as the peak flow rate divided by the cross-sectional area of the riser pipe v = Q max A pipe , and d (L) is the geometry specific characteristic length, which we assumed to be the effective hydraulic aperture of the fracture (d = 2b), consistent with previous studies (Nicholl et al. 1999;Zimmerman et al. 2004;Ranjith and Darlington 2007).We used the transmissivity estimates returned through our ideal fracture analysis to estimate the effective hydraulic aperture (b) of the tested fracture (Equation 12), assuming validity of the local cubic law (Zimmerman and Bodvarsson 1996): where T (L/T 2 ) is fracture transmissivity, μ (M/[LT]) is the dynamic viscosity, ρ (M/L 3 ) is the fluid density and g (L/T 2 ) is acceleration due to gravity.Our analysis returned effective hydraulic aperture estimates from 0.1 to 0.4 mm (Figure 8), which is consistent with previously reported aperture distributions (Hakami and Larsson 1996).Using the effective aperture estimates, we calculated a maximum Reynolds number of 0.7 (Figure 8), which is more than one order of magnitude lower than the typical threshold of 10, where turbulent flow effects are commonly observed (e.g., Ranjith and Darlington 2007).
As another approach to investigate the presence of non-Darcian flow effects during our hydraulic testing, we conducted a series of oscillatory flow tests at a constant oscillation period with increasing peak flow rates and monitored the change in observation signal amplitudes.If our assumption holds that Darcy's Law is valid during our field experiments, we expect to see a positive linear correlation between the applied peak flow rate and measured observation signal amplitude.As an illustrative example, we present analysis from oscillatory flow tests conducted at oscillation periods of 10, 30, and 60 s.For each oscillation period, we plotted the peak flow rate (Q max ) versus observation head amplitude (| h|) at the observation wells and used linear regression to determine the line of best fit that goes through the origin for each individual well (Figure 9).This analysis shows head amplitudes increasing linearly with the peak flow rate at all wells across the range of explored oscillation periods (Figure 9).These results support our Reynolds number analysis above and taken together support our assumption that Darcy's law is valid during our field experiments.These results suggest that the presence of non-Darcian flow is not the mechanism causing the period-dependent fracture flow parameters returned during our ideal fracture analysis.

Leaky Fracture Analysis
To investigate the presence of fluid exchange between the fracture and surrounding host rock, we conducted oscillatory flow tests with a pressure stimulation in the tested fracture, while monitoring for a pressure response in the surrounding bedrock.To accomplish this, we located four packer strings at the level of the tested fracture, and we placed one packer string in competent bedrock (well A1), approximately 1.5 m below the tested fracture (Figure 10).We introduced the pressure stimulation in the isolated fracture (well B2) while monitoring for the presence of a measurable pressure response in the surrounding host rock.Analyzing the recorded pressure signals at all wells, we note a largely attenuated observation signal in the host rock-relative to observation signals measured in the isolated fracture-with an amplitude of approximately 5 mm that is readily identifiable above the recorded instrument noise (Figure 10).These results illustrate that we can propagate a pressure signal outside of the fracture and provides the motivation to employ a leaky fracture analytical model for initial exploration of fluid exchange as a potential mechanism leading to the observed period dependence.
Under a leaky fracture conceptual model, we invoked the same assumptions as the ideal fracture analysis described above, with the exception that we now allow for vertical flow into and out of the fracture along pressure gradients.We used the leaky aquifer analytical model (Equation 13) developed by Rasmussen et al. (2003) as our forward model, which outputs Fourier coefficients for a given set of flow parameter inputs: where B 2 (L 2 ) is T /L and L (T −1 ) is the fracture leakance, which represents the vertical flow of fluid into and/or out of the fracture.This model is identical to the Hantush and Jacob (1955) leaky confined aquifer model, with the exception of the oscillatory flow rate.The additional assumptions in the leaky fracture model are that flow is vertical into and out of the fracture due to a large hydraulic conductivity contrast between the fracture and host rock, and the surrounding host rock is assumed to be incompressible (i.e., no fluid storage in the host rock).
Using this leaky fracture model as our forward model during inversion, we minimized the objective function (Equation 8) to determine the effective transmissivity (T ), storativity (S), and leakance (L) of the isolated fracture that best fit the observed data.It is important to note that this forward model contains three unknown parameters, leading to an underdetermined inverse problem when using a single-frequency analysis.We overcame this nonuniqueness by using the multi-frequency inversion approach described in the ideal fracture analysis above.
We used all oscillatory flow tests with a pumping period of 300 s to determine a single set of effective flow parameters for the leaky fracture.This approach allowed for a direct comparison of effective flow properties and modeled data fit with the ideal fracture analysis (Figure 7B) to determine if the addition of a leakance term improved modeled data misfit relative to our ideal fracture analysis above.

Borehole Storage Impacts
Finally, we explore the impacts of borehole storage at the observation well as a potential mechanism leading to the observed period-dependence.We considered borehole storage at the observation well to occur within the riser pipe used to install packer strings and through which water is cycled during oscillatory flow testing, which is consistent with the borehole instrumentation used during our field experiments.We are unaware of any analytical models that account for borehole storage under oscillatory flow conditions; therefore, we developed a simple two-dimensional phase-domain numerical model to generate synthetic data (i.e., Fourier-coefficients) and then used the generated data to estimate the flow parameters following the approach described in our ideal fracture analysis.
We generated a square modeling domain 3000 m × 3000 m in size to simulate fluid flow through a semi-infinite, nondeforming, constant-aperture fracture bounded by an impermeable host rock.We designated the lateral boundaries as constant-head boundaries (h = 0 m), and we prescribed the model top and bottom boundaries as no-flow boundaries.We discretized the domain into 2.25 million grid cells with a constant discretization ( x, y) of 2 m in the x -and y-directions.We arranged wells throughout the modeling domain consistent with the (FSR) 2 well geometry (Figure 2A).For these modeling simulations, well A1 serves as the pumping well and wells B1, B2, B3, and B4 serve as the observation wells.
The model parameter inputs are described in Table 3. First, we selected a representative hydraulic fracture aperture of 0.3 mm, based on our field analysis at the longest oscillation periods (Figure 8A), and we then assigned a homogeneous transmissivity (Table 3) to the modeling domain through using the local cubic law (Equation 12).We also assigned a homogeneous storativity throughout the modeling domain, assuming that fluid storage along the fracture is a linear elastic response due to fluid compressibility alone (Equation 14) where ρ (M/L 3 ) is the fluid density, g (LT −2 ) is gravitational acceleration, and β (M −1 LT 2 ) is fluid compressibility.
To include the effects of borehole storage at the observation wells, we calculated an effective storativity that includes the effects of fluid storage within the riser pipe as well as the fracture.The change in water volume along the riser pipe in response to a unit head change is given by: where r pipe (L) represents the radius of the riser pipe through which water is being cycled.The change in water volume along the fracture due to a unit head change is given by: We averaged the volume of water released from the fracture and riser pipe to determine the effective storage coefficient (S eff ) and applied that to grid cells within the model domain that contain an observation well.
We simulated the steady-periodic-constant amplitude and phase-phase-domain response at the observation wells using OHT3DINV (Cardiff 2016), which is a MATLAB-based numerical forward model.OHT3DINV takes a finite volume approach to simulate the Fourier coefficients throughout the modeling domain by solving the steady-periodic formulation of the groundwater flow equation (Equation 18)-derived from the confined groundwater flow equation-subject to the boundary conditions described above.
where Q ω (T −1 ) represents the Fourier coefficients for the pumping source.As a first step, we ensured the accuracy of the developed homogeneous model by comparing the simulated Fourier coefficients at all observation wells to the Fourier coefficients output by our ideal fracture model (Equation 10) using oscillation periods from 10 to 500 s.We found the misfit between the numerical and analytical approaches to be consistent with the estimated data measurement error during our field experiments (Table 1), providing confidence in the accuracy of our developed numerical model.Now considering the impact of borehole storage effects, we simulated oscillatory flow tests with oscillation periods from 10 to 500 s using the described numerical model with input parameters given in Table 3 to generate synthetic data (i.e., Fourier coefficients) that we then inverted for effective fracture flow parameters.For this analysis, we considered only the well pair (A1 to B3) with the shortest inter-well spacing of 4 m (Figure 2), where we expect observation amplitudes to be least attenuated and thus borehole storage effects to be the greatest.Because we are most interested in the impact of borehole storage, we assumed we can perfectly measure a noisefree signal.Using the generated data, we individually estimated the effective fracture flow properties for each simulated oscillatory flow test using the ideal fracture analysis approach described above.
When considering the effects of borehole storage, we found period-dependent hydraulic diffusivity estimates that decrease with increasing oscillation period (Figure 12), consistent with our field analysis (Figure 5).However, our results show decreases in both estimated transmissivity and estimated storativity with increased pumping period (Figure 12), in contrast to the trends noted in our field analysis.These results indicate the decreasing diffusivity is driven by the large decrease in estimated transmissivity (Figure 12), as opposed to significant increases in storativity estimates seen in our field analysis (Figure 5).

Discussion and Conclusion
This work represents the first published field application of oscillatory flow interference testing to characterize a single fracture embedded in poorly cemented sedimentary bedrock with considerable primary porosity, and the largest known dataset of its type.This work is also the first attempt to analyze a robust experimental dataset with multiple conceptual models to investigate multiple plausible hydrologic processes that might contribute to the apparent period-dependence of estimated hydraulic properties.The results presented in this work eliminated multiple candidate mechanisms that have been hypothesized to contribute to the period-dependent flow parameters and provide an important lens into unexplored hydraulic processes occurring within fractures that might be causing observed pressure responses and the resulting perioddependent effective fracture flow parameters.The results of this analysis highlight the necessity of choosing the appropriate conceptual model when characterizing fractured sedimentary bedrock, challenging the way typical idealized conceptual models are applied in many fracture characterization efforts.Increased understanding of the complex hydraulic processes occurring in fractured sedimentary bedrock can lead to improved characterization of reservoir features that are likely to be of outsized importance in deep aquifer applications such as CCS and fluid waste storage.
Hydraulic parameters in fractured bedrockspecifically hydraulic conductivity or permeabilityexhibit a well-described scale effect when determined in an effective, or homogeneous, sense (e.g., Neuman 1994;Rovey and Cherkauer 1995;Sánchez-Vila et al. 1996;Vesselinov et al. 2001;Illman 2006;Illman and Tartakovsky 2006).This scale effect has been primarily observed in complex fracture networks, and has been attributed to the increasing connectivity of high transmissivity fractures, which control flow through the fracture network, with increasing scale (Vesselinov et al. 2001;Martinez-Landa and Carrera 2005;Illman 2006;Illman and Tartakovsky 2006).While our storativity estimates show a prominent scale effect (Figure 5C), we do not observe a significant increase in our transmissivity estimates (Figure 5B).These results are consistent with results reported by Guiltinan and Becker (2015), which showed a factor of 2 decrease in estimated transmissivity with oscillation periods from 60 to 240 s.Neuman (1994) suggested that a lack of permeability scaling could result under purely radial flow conditions.However, our results showed the presence of fluid exchange between the fracture and host rock (Figure 10), suggesting that fluid flow is not likely to be purely radial during our field experiments.This leads us to conclude that the lack of transmissivity scaling in our analysis arises from the fact that we tested a single fracture and therefore we sampled aperture variability that likely varies over much smaller scales than network connectivity.
Though our analytical modeling approach does not explicitly allow the investigation of heterogeneity, we note multiple pieces of evidence pointing to the impact of heterogeneity on our recorded observation signals.First, using a multi-frequency inversion approach with our ideal fracture model we are not able to find a set of homogeneous or "effective" flow parameters that adequately explains multiple oscillatory flow tests within expected data measurement error across the range of tested oscillation periods and radial distances (Figure 7).Additionally, using the same multi-frequency approach with the addition of a leakance parameter (i.e., our leaky fracture model) does not significantly improve modeled data misfit (Figure 11).These results suggest that the assumption of aperture homogeneity in the chosen forward models is not valid for our analysis.The observed variability in transmissivity at a given pumping period (Figure 5B) further suggests the potential impact of heterogeneity on our recorded observation signals.If our homogeneity assumption was valid, we would expect consistent transmissivity estimates for multiple tests conducted at a consistent oscillation period and inter-well spacing.Taken together, these results indicate that homogeneous analytical models do not capture the necessary complexities occurring within an individual fracture under oscillatory flow conditions and motivate the exploration of more complex numerical modeling approaches using stochastically generated aperture realizations to understand the impacts of heterogeneity when using oscillatory flow tests to characterize bedrock fractures.A simple example of such modeling analysis creating the apparent period-dependent trends we observe in our ideal fracture analysis is included in Supporting Information (Figure S2); however, a more thorough exploration is necessary to fully understand the impact of heterogeneity on recorded oscillatory flow observation signals and the resulting parameter estimates.
Fluid exchange between the fracture and surrounding host rock remains an underexplored potential mechanism for the observed period dependence in effective fracture flow parameter estimates.Despite the volume-and thus oscillation period-limitations of our pneumatic testing approach, we successfully propagated a measurable pressure signal into the host rock using a pressure stimulation in the fracture during our field experiments, indicating that fluid exchange between the fracture and host rock does occur during oscillatory flow tests at the longest oscillation periods.We see the impacts of this fluid exchange in our returned storativity estimates (Figure 5C), which are two to three orders of magnitude larger than what would be expected due to fracture storativity alone (Rutqvist et al. 1998).Despite these observations, we found that the addition of a leakance parameter did not significantly improve the modeled data misfit (Figure 11) relative to the ideal fracture model (Figure 7).We interpret this to be an impact of the relative model insensitivity to the leakance parameter relative to T and S (Table 2), which is consistent with previous modeling studies (Patterson and Cardiff 2022).These results indicate that the leaky fracture analytical model does not adequately explain the period-dependent fracture flow parameters, and a more complex, three-dimensional numerical modeling approach is required to fully explore the impacts of fluid exchange on effective parameter estimates under oscillatory flow conditions.
The impact of non-ideal borehole conditions, such as borehole storage can represent a common source of parameter error during pressure-based hydraulic characterization efforts.Our modeling results indicated that the presence of borehole storage effects at the observation well can generate period-dependent fracture flow parameter estimates (Figure 12).Specifically, we demonstrated that borehole storage produces transmissivity estimates that decrease and storativity estimates that remain approximately constant with increasing oscillation period.These period-dependent trends are in contrast with the perioddependent trends returned in our ideal fracture analysis (Figure 5), leading us to eliminate borehole storage as a mechanism contributing to the trends in our field observations.This conclusion is consistent with and supported by previous field studies that demonstrated borehole storage effects are not a primary mechanism leading to the observed period-dependence (Renner and Messar 2006).
The results presented in this work demonstrate that oscillatory flow testing helps to illuminate important processes-not represented with idealized analytical models-that can influence pressure propagation (and thus fluid flows) through rock fractures.Proposed explanations for anomalous fracture behavior, including non-Darcian flow, fluid exchange with the surrounding host rock, or borehole storage effects, are deemed insufficient to explain the observed behavior under oscillatory flow conditions.Said another way, none of the common models that conceptualize fractures as infinite, uniform-aperture features appear sufficient to explain the observations from oscillatory flow tests conducted at multiple oscillation frequencies.
In addition to the potential role of fracture heterogeneity, fractured bedrock characterization efforts to date-including this one-do not explore the role of fracture hydromechanical behavior on the observed period-dependence.Our analysis shows observation signal amplitude increasing linearly with increasing peak pumping rate, indicating a lack of excess pressure losses due to hydromechanical behavior (Quinn et al. 2011); however, nanometer scale fracture displacement has been measured under oscillatory flow conditions with observation head amplitudes smaller than those used during our field experiments (Becker et al. 2017).Exploring the nonlinear relationship between head change and storage (Rutqvist and Stephansson 2003;Cappa et al. 2008) in deformable fractures as a mechanism leading to the observed perioddependence requires more complex numerical modeling approaches that are beyond the scope of this work, but represents an area of future research currently being explored.

Figure 1 .
Figure 1.Conceptual schematic of borehole instrumentation for oscillatory flow interference testing.The gray rectangle represents the tested fracture.Inset shows unmodeled flow conditions explored using simple analytical modeling approaches.
Figure 2. (A) Aerial view of the (FSR) 2 with insets showing the field site location within Wisconsin and well geometries with inter-well distances indicated.(B) Hydrostratigraphy at the (FSR) 2 .The blue dashed lines indicate laterally extensive, horizontal, bedding plane fractures (not all-inclusive), and the dashed red line represents the tested fracture at 32 m depth.Vertical black lines represent boreholes used for hydraulic testing.

Figure 3
Figure 3. (A) Recorded stimulation head signal at well A1 (blue) and volumetric flow rate for a representative oscillatory flow test with 20 s oscillation period.(B) Recording observation signals showing attenuated and delayed versions of the stimulation signal.

Figure 4 .Figure 5 .
Figure 4. Representative oscillatory flow test observation signal with 180 s oscillation period.The gray line represents the measured signal with associated noise, and the red line represents the extracted noise-free signal, highlighting the ability to extract noise-free pressure signals from noisy signals with low signal-to-noise ratios.

Figure 6 .
Figure 6.Measured (red) and modeled (black) observation head signals at 8 m inter-well spacing for oscillatory flow tests at 20, 60, 120, and 300 s oscillation periods highlighting the quality of data fit with fracture flow parameter estimates.

Figure 7 .
Figure 7. Modeled versus observed Fourier coefficients using a multi-frequency inversion approach for all tests conducted at (A) a distance of 16 m and (B) a stimulation period of 300 s.
Figure 8. (A) Effective hydraulic aperture estimates and (B) Reynolds numbers for each oscillatory flow test suggesting Darcian flow regimes during oscillatory flow testing.
Figure 11.Cross-plot of measured and modeled Fourier coefficients using a multi-frequency approach with a leaky conceptual model and all oscillatory flow tests conducted at an oscillation period of 300 s.

Figure 12 .
Figure 12.Estimated fracture flow parameters for the borehole storage model showing decreasing diffusivity (A), decreasing transmissivity (B), and decreasing storativity (C) with increasing stimulation period.