Identification of dominant gas transport frequencies during barometric pumping of fractured rock

We demonstrate that although barometric pressures are complicated signals comprised of numerous frequencies, it is a subset of these frequencies that drive the overwhelming majority of gas transport in fractured rock. Using an inverse numerical analysis, we demonstrate that a single barometric component with seasonally modulated amplitude approximates gas transport due to a measured barometric signal. If past barometric tendencies are expected to continue at a location, the identification of this frequency can facilitate accurate long term predictions of barometrically induced gas transport negating the need to consider stochastic realizations of future barometric variations. Additionally, we perform an analytical analysis that indicates that there is a set of barometric frequencies, consistent with the inverse numerical analysis, with high production efficiency. Based on the corroborating inverse numerical and analytical analyses, we conclude that there is a set of dominant gas transport frequencies in barometric records.


Identification of dominant gas transport frequencies during barometric pumping of fractured rock Dylan R. Harp , John P. ortiz & Philip H. Stauffer
We demonstrate that although barometric pressures are complicated signals comprised of numerous frequencies, it is a subset of these frequencies that drive the overwhelming majority of gas transport in fractured rock. Using an inverse numerical analysis, we demonstrate that a single barometric component with seasonally modulated amplitude approximates gas transport due to a measured barometric signal. If past barometric tendencies are expected to continue at a location, the identification of this frequency can facilitate accurate long term predictions of barometrically induced gas transport negating the need to consider stochastic realizations of future barometric variations. Additionally, we perform an analytical analysis that indicates that there is a set of barometric frequencies, consistent with the inverse numerical analysis, with high production efficiency. Based on the corroborating inverse numerical and analytical analyses, we conclude that there is a set of dominant gas transport frequencies in barometric records.
Predictions of gas transport can be improved for many applications by better understanding the effect of barometric variations on gas transport in fractured rock. Barometric variations push gases deeper into the fractured rock during barometric highs and pull gases upward during barometric lows [1][2][3][4][5] . Figure 1 contains a schematic conceptualization of this process with details that will be referred to throughout this paper. Fractures provide fast pathways for gas transport while the rock matrix provides relatively immobile gas storage in between barometric cycles, allowing for a ratcheting mechanism that greatly enhances gas transport. Barometric variations can drive leakage from CO 2 sequestration sites 6-10 , leakage of methane from hydraulic fracturing operations 11 , radon gas entry into buildings 12 , and radionuclide gas seepage from underground nuclear explosions and waste storage [13][14][15][16][17][18][19] .
One of the primary uncertainties in making predictions of barometrically-induced gas transport is that the future barometric variations are highly uncertain and cannot be accurately forecasted past a few weeks. Barometric pressure fluctuations are complex signals composed of many frequencies driven by multiple processes including atmospheric tides, weather patterns, and seasonal and annual cycles and are highly dependent on latitude and elevation. Many previous researchers have assumed that barometric variations due to roughly weekly weather patterns are important frequencies to consider for gas transport in fractured rock. Nilson et al. 1 identified that the barometric period with the highest amplitude at the National Nuclear Security Site (NNSS) was approximately 7.2 days based on historical data. Mourzenko et al. 20 use a synthetic sinusoidal barometric signal with period of approximately 7.3 days to represent weather patterns at Roselend Natural Laboratory in France for their numerical investigations. Neeper 3 also assumed that weekly weather patterns are important for gas transport in fractured rock. However, in our investigations, a spectral analysis of a barometric record results in many frequencies with similar amplitude making it difficult to identify a single dominant barometric frequency. We also demonstrate that the barometric sinusoidal component with the largest amplitude will not necessarily be associated with the dominant gas transport frequency.
This paper presents the first detailed analysis of the identification of the dominant gas transport frequency through corroborating inverse numerical and analytical analyses. Using barometric pressures measured in Anchorage, AK from 2014 through 2017, we first present the decomposition of a barometric pressure record into sinusoidal components (amplitude/frequency pairs). Anchorage was chosen since it is representative of a location (2019) 9:9537 | https://doi.org/10.1038/s41598-019-46023-z www.nature.com/scientificreports www.nature.com/scientificreports/ with large barometric variations. Then, we present an inverse analysis using a numerical model of gas transport within a fractured domain demonstrating that a dominant transport frequency can be identified that reproduces gas transport simulated using the measured barometric pressures. Next, we perform an analytical analysis of gas production efficiency using the decomposed barometric pressure record. From this we identify a set of highly efficient frequencies for producing gas that are clustered around the frequency identified in the inverse numerical analysis. Finally, we provide conclusions based on the corroborating evidence from the inverse numerical and analytical analyses.

Barometric Pressure Decomposition
We obtained hourly barometric pressure data from Anchorage, AK from 2014 through 2017 from Weather Underground (www.wunderground.com). In order to obtain a uniform one hour spacing, we linearly interpolated a few missing measurements and removed a few extra intra-hour measurements in the record (top plot in Fig. 2). We decomposed the data into the frequency domain using a Fast Fourier Transform (FFT) algorithm 21 in the bottom plot of Fig. 2, where the barometric period is T = 2π/ω, where ω is the barometric frequency. From this plot, it is apparent that the barometric signal is a complicated combination of many components (frequency/ amplitude pairs). While a period of around 24 days has the largest amplitude, our analysis below will demonstrate that this is not the dominant gas transport frequency. Therefore, unlike the case in Nilson et al. 1 , where the dominant gas transport frequency is identified as the frequency with the largest amplitude, it is not possible to identify the dominant gas transport frequency simply based a spectral analysis alone.

Numerical Modeling Analysis
We use a 2D numerical model to simulate air flow and gas transport and immobile pore-water storage of dissolved gas within a partially-saturated, fractured rock with 1 mm vertical fractures separated by 10 m. The  III I  I II I I I I  I II I I I  II I I I II I I I  II I I I I I I I  I I I I I I I I I I I I I I I I I I  www.nature.com/scientificreports www.nature.com/scientificreports/ is implemented in the PFLOTRAN simulator 22 with modifications to allow for kinetic dissolution and exsolution of gas between air and immobile pore-water. The model is described in detail and verified against a suite of analytical solutions in Harp et al. 5 .
Identifying the dominant gas transport frequency. We calibrate the synthetic pressures using a Levenberg-Marquardt optimization approach 23 implemented in the PEST software package 24 . An initial calibration using a single frequency synthetic barometric signal with constant amplitude resulted in a seasonal mismatch in concentrations (for details, refer to Section S1 in Supporting Material). By inspecting the top plot in Fig. 2, it is apparent that the barometric pressures during the summer months have lower amplitudes than during winter months. We therefore defined a seasonally modulated synthetic barometric pressure signal as where d  is the mean amplitude of the dominant gas transport frequency, s  is the amplitude of the seasonal modulation, ω d is the dominant gas transport frequency, ω s is the seasonal modulation frequency (T s = 1 year, where ω s = 2π/T s ), and γ d and γ s are the phase shift of the dominant gas transport frequency and seasonal modulation, respectively, and is a vector containing the calibration parameters. The s ) term captures the seasonal modulation about the mean of the dominant gas transport frequency. The objective function F minimized in the calibration is where C i m and C i s are the ith tracer concentrations at the ground surface driven by the measured and synthetic barometric signals, respectively. We plot the calibrated seasonally modulated barometric signal in the top plot of Fig. 3 as the red line along with the measured barometric record as the black line. In the bottom plot of Fig. 3, concentrations simulated using the synthetic barometric signal (red line) and measured barometric record (black line) are plotted, demonstrating that a seasonally modulated single frequency barometric signal is able to closely approximate concentrations driven by a measured barometric record. This calibration reduced the standard error of the concentration residuals by approximately half compared to the calibration with constant amplitude barometric frequency (Section S1 in Supporting Material).
The calibrated parameter values are  = 1310 d Pa, T d = 7.29 days,  = 422 s Pa, and γ s = 1.85 radians. We calibrated γ d , but the calibration was insensitive to its value because the gas transport in our example is a longer term process resulting from a succession of barometric cycles, and therefore identifying a phase shift is not critical to match the overall concentration trend.

Sensitivity of calibrated barometric parameters to subsurface domain scenarios. Based on an
analysis of the sensitivity of tracer concentrations at the top of the fracture to subsurface domain parameters, we identified that the concentrations are most sensitive to the depth to impermeable layer L and matrix permeability k m (refer to Section S2 in Supporting Material for details). Therefore, since L and k m have the largest potential to alter the calibration of the barometric parameters ( T , , ,

Analytical Analysis
The mass discharge efficiency of a barometric frequency is dependent on its ability to push and retract atmospheric air into and out of fractured rock and allow sufficient time for the diffusive exchange of gas into this air while at the depth of the tracer, and retain the gas during the return trip to the atmosphere (i.e., not lose too much gas due to diffusion into matrix without tracer). The mass discharge efficiency (η M ) can be defined as the mass of tracer removed from the subsurface during one barometric cycle (ΔM c ) relative to the original mass of tracer in the subsurface (M 0 ) as where, assuming a linear concentration gradient with depth, M 0 can be approximated as where C B is the concentration of the tracer gas and V 0 is the volume of air with tracer 1 . The dominant gas transport frequency will not only depend on η M but will also depend on the frequency of occurrence of the barometric component. For example, a barometric component with lower η M may remove more gas from the subsurface over time than a less frequent component with higher η M simply because it occurs more often. We define the production efficiency (η P ) as a metric for identifying dominant gas transport frequencies by multiplying η M by the frequency as www.nature.com/scientificreports www.nature.com/scientificreports/ In order to approximate ΔM c analytically, the exchange of fresh air with air containing tracer can be conceptualized as occurring within a packet of air with volume ΔV that travels down the fracture and a portion of the matrix wall that has been invaded by incoming fresh air 1 . Figure 1 illustrates this conceptualization and many of the quantities used below. The thickness of this packet of air in half-space considering both the half fracture thickness and the invaded region of the associated matrix half block (d) is = +  Fig. 1) as where W a is the Womersley number modified to account for the porosities of the given problem defined as , where φ f is one in our case, and = − i 1 . If the subsurface were in perfect equilibrium with the surface, i.e., the pneumatic diffusivity were extremely large so that pressure variations at the ground surface are immediately transferred throughout the subsurface, the ratio of ΔV to the total volume of subsurface air with tracer V 0 would be proportional to the ratio of the change in pressure (amplitude) Δp to the mean static pressure p 0 as Based on this relationship, the conceptual maximum volume of air that could be extracted from the subsurface due to a single cycle of a barometric component is www.nature.com/scientificreports www.nature.com/scientificreports/ where the variables have been defined earlier. However, due to finite pneumatic diffusivity in the fracture and matrix, perfect pressure equilibrium is never achieved throughout the subsurface. Using a double porosity (fracture/matrix) analytical solution, Nilson et al. 1 derive an analytical solution for ΔV as where the λ's are dimensionless Fourier numbers that define the ratio of the diffusive transport rate to the storage rate defined as Inserting Eq. 10 into Eq. 7, and then Eqs 5 and 7 into Eq. 6, η P can be expressed as This equation can also be derived from the breathing efficiency (η B ), which quantifies the volume of air that a barometric cycle is able to extract (ΔV) relative to the maximum volume that could be removed if the subsurface were in perfect equilibrium with the atmosphere (ΔV max ), defined as Refer to Fig. 1 for depictions of η B and η D . Nilson et al. 1 combine η B and η D into the overall transport efficiency η as (17) B D η η η = to describe the efficiency of a single cycle of a barometric component to extract gas from the subsurface (i.e., considering a single cycle in isolation, not the ability of the barometric component to extract gas from the subsurface per unit of time). As such, production efficiency can be equivalently expressed as (18) P B D η η η ω ηω = = .
In Fig. 5, we plot the breathing efficiency (η B ), diffusive exchange efficiency (η D ), overall transport efficiency (η), and production efficiency (η P ) calculated for the components of the measured barometric record as a function of period. In the top plot, it is apparent that η B increases monotonically with barometric period (decreases with frequency). It is not dependent on amplitude, which cancels out during its derivation. Therefore, η B effectively quantifies the ability of lower frequency components to more effectively penetrate into fractured rock irrespective of their amplitude. In general, η D decreases with increasing period, however, the combination of period and amplitude leads to more nuanced (non-monotonic) behavior. The combination of these conflicting efficiencies in η leads to a maximum efficiency at around 24.3 days (refer to the vertical dashed line in Fig. 5). Referring to the bottom plot of Fig. 2, this is the period with the largest amplitude. As discussed previously, η quantifies the gas transport efficiency of a single cycle of a barometric component and does not account for the frequency of occurrence of the barometric component. Therefore, it is not appropriate for identifying the dominant gas transport barometric component.
www.nature.com/scientificreports www.nature.com/scientificreports/ In the bottom plot of Fig. 5, we plot the production efficiency (η P ; Eq. 14) as a function of barometric period. In this case, when the frequency of occurrence of the barometric component is taken into account, it is apparent that there is a cluster of high efficiency amplitude/frequency pairs around a period of 7.3 days (refer to the inset in Fig. 5), consistent with the inverse numerical analysis above. However, the results are complicated by the highest production efficiency occurring at 0.5 days. To investigate this, in Fig. 6, we plot the cumulative mean period of the amplitude/frequency pairs sorted in order of decreasing production efficiency. This demonstrates that although the highest production efficiency occurs at 0.5 days, the average period of the highest efficiency amplitude/frequency pairs hovers around 7.3 days. Therefore, it is likely that for the Anchorage, AK data that we analyzed, 7.3 days is the average period of a set of high production efficiency periods that lead to the vast majority of the gas transport.  www.nature.com/scientificreports www.nature.com/scientificreports/

Conclusions
The combined inverse numerical and analytical analyses presented in this manuscript support the following conclusions: 1. There is a set of barometric frequencies responsible for the vast majority of gas transport in fractured rock. 2. A single barometric frequency with seasonally modulated amplitude can be used to accurately predict the barometrically-induced gas transport from a measured barometric signal. 3. The dominant gas transport frequency is the average of the high production efficiency amplitude/frequency barometric sinusoidal components.
For practical applications, these conclusions indicate that as long as future barometric pressures at a location have similar characteristics to the past, it is possible to predict future gas transport using a single seasonally modulated barometric frequency. This eliminates the need to consider an ensemble of possible future barometric variations, significantly simplifying predictions of gas transport in fracture rock and subsequent breakthrough times.

Data Availability
The data used in this research is freely available from Weather Underground (www.wunderground.com).