Continuous Dissolved Gas Tracing of Fracture‐Matrix Exchanges

Transport in fractured media plays an important role in a range of processes, from rock weathering and microbial processes to contaminant transport, and energy extraction and storage. Diffusive transfer between the fracture fluid and the rock matrix is often a key element in these applications. But the multiscale heterogeneity of fractures renders the field assessment of these processes extremely challenging. This study explores the use of dissolved gases as tracers of fracture‐matrix interactions, which can be measured continuously and highly accurately using mobile mass spectrometers. Since their diffusion coefficients vary significantly, multiple gases are used to probe different scales of fracture‐matrix exchanges. Tracer tests with helium, xenon, and argon were performed in a fractured chalk aquifer, and resulting tracer breakthrough curves are modeled. Results show that continuous dissolved gas tracing with multiple tracers provides key constrains on fracture‐matrix interactions and reveal unexpected scale effects in fracture‐matrix exchange rates.

The signature of fracture-matrix exchanges on tracer breakthrough curves is a strong late-time tailing, characterized by a power law decay of late time concentrations (Haggerty et al., 2000;Haggerty & Gorelick, 1995;Hyman et al., 2019;Kirchner et al., 2000;Maloszewski & Zuber, 1993;Małoszewski & Zuber, 1985;Tsang, 1995;Weede & Hötzl, 2005). For idealized parallel plate fractures and a homogeneous matrix, the expected late time decay is (t)~t −b , with b ¼ 1.5 (Carrera et al., 1998;Haggerty et al., 2000). However, a broad range of power law exponents b have been observed, depending on the context, and significance in terms of flow heterogeneity within preferential flowpaths and multiple rates of mass transfer in the matrix is still debated and difficult to assess (Becker & Shapiro, 2003;Carrera et al., 1998;Dentz & Berkowitz, 2003;Haggerty et al., 2001;Hyman et al., 2019;Kang et al., 2015;. The use of tracers with multiple diffusion coefficients (Becker & Shapiro, 2000) and combining different forced gradient configurations such as convergent and push-pull (Becker & Shapiro, 2003;Kang et al., 2015) have been shown to considerably reduce the uncertainty for the inference of relevant transport physics from tracer tests. But the range of diffusion coefficients of fluorescent dyes that can be continuously monitored in situ is limited, and other solute tracers require heavy sampling and analysis. To overcome these challenges, this study investigates the potential of in situ continuous dissolved gas tracing to quantify fracture-matrix interactions.
The potential use of dissolved gases as tracers in groundwater studies was established around 20 years ago (Sanford et al., 1996;Solomon et al., 1998). Due to their higher molecular diffusion coefficients, interactions of dissolved gases (e.g., D ¼ 58.8 × 10 −10 m 2 s −1 for helium 12°C; Jähne et al., 1987) with the rock matrix are stronger compared to classical tracers (e.g., D ¼ 4.5 × 10 −10 m 2 s −1 for uranine; Skagius & Neretnieks, 1986). Sanford et al. (2002) combined bromide (D ¼ 25 × 10 −10 m 2 s −1 ) and helium in a tracer test performed in fractured rock to calibrate a fracture-matrix model. Although uncertainty was large due to low mass recovery, the larger diffusion coefficient of helium was shown to delay and dilute the breakthrough curve peak compared to bromide, highlighting the effect of matrix diffusion. While these early dissolved gas studies relied on field sampling and laboratory measurements, the recent development of portable mass spectrometer now allows on-site continuous measurements of multiple gas tracers with high accuracy (Brennwald et al., 2016;Chatton et al., 2017). Innovative techniques based on membrane inlet mass spectrometry (MIMS) enable the in situ measurement of reactive and noble dissolved gases in oceanic and continental waters (Bell et al., 2007;Brennwald et al., 2016;Camilli & Hemond, 2004;Cassar et al., 2009;Chatton et al., 2017;Mächler et al., 2012;Roques et al., 2020;Tortell, 2005).
Here we explore the potential of this technology in active mode. Forced gradient tracer test experiments were performed in a dual porosity porous/fractured chalk aquifer using helium, argon, and xenon as dissolved gases and uranine as fluorescent dye, characterized by diffusion coefficients of in total one order of magnitude difference. After transport in the fracture network, tracer recovery was measured continuously on-site using a portable mass spectrometer for dissolved gases and a fluorimeter for uranine. This technique first allows the possible measurement of multiple tracers with different diffusion coefficients probing different time scales of fracture-matrix interactions. Second, the high accuracy of measurements at low concentration allows a fine analysis of the late time breakthrough curve tailing, which provides critical information on the physics of fracture-matrix exchanges. Using analytical solutions of tracer transport in fractures and diffusion into the matrix, the study shows that continuous high accuracy measurement of tracers of different diffusion coefficients allows inferring unexpected transport dynamics in fracturematrix systems.

Test Site and Tracer Properties
Injection of a fluorescent dye tracer (uranine) jointly with dissolved gases was performed in a fractured chalk aquifer (supporting information Text S1 and Figure S1). This aquifer is considered as a typical dual porosity, dual permeability reservoir with a total porosity of the porous matrix estimated to about 40%, including about 1% of fractures (Batlle Aguilar et al., 2007;Brouyère, 2001;Dassargues, 2018;Goderniaux et al., 2018;Rorive & Goderniaux, 2014). While fractures allow fast preferential flows and explain most of the aquifer transmissivity, the matrix porosity is characterized by a lower hydraulic conductivity but enables the storage of large quantities of groundwater (Dassargues, 2018;Dassargues & Monjoie, 1993;Downing et al., 1993).

CF-MIMS-Dissolved Gas Measurement
A membrane inlet mass spectrometry (MIMS) technique is used for monitoring dissolved gas concentration on-site. For this purpose, an EQuilibrium (EQ-) Membrane system (Brennwald et al., 2016) or a Continuous Flow (CF-) Membrane Inlet Mass Spectrometer (MIMS) is suitable (Chatton, 2017;Chatton et al., 2017). The system uses an inlet system developed to enhance its sensitivity and response time to "light" gases (<50 amu especially <20 amu), and previous tests showed that this selectivity has a price for "heavy" gases (>100 amu) which sensitivity decrease and response times increase with the gas mass. Response times are typically below 30 s for light gases and reach up to a few minutes for detectable heavy gases. In this study, the heavy gas xenon is therefore affected by a delay in measurements changing the amplitude for small observation times, while the light gases argon and helium are not affected. With more time, these effects are less significant, so that the tailing of xenon can be still used.

Experimental Setup
The two wells Pz1 and Pz2 (distant 7.55 m) are connected by open horizontal fractures, and two experiments were performed with joint injections of dissolved gases and uranine in a single fracture located at 34.8 m deep in Pz2, using an inflatable double packer system: (1) a convergent forced gradient experiment with recovery in Pz1 and (2) a push-pull experiment in Pz2 ( Figure S1). In the convergent experiment, 89.46 L of water, containing 2.0 g of uranine, were partially saturated with helium, xenon, and argon corresponding to concentrations of 2.0 × 10 −5 mol L −1 (0.007 g), 1.7 × 10 −5 mol L −1 (0.20 g), and 1.8 × 10 −3 mol L −1 (6.37 g), respectively. All tracers were injected jointly in Pz2 at a flow rate of 0.57 m 3 hr −1 and extracted continuously in Pz1, at a constant flow rate of 7.2 m 3 hr −1 . Uranine and dissolved gases were continuously measured using a field fluorimeter and the CF-MIMS, respectively. In the push-pull experiment, 92.63 L of water, containing 0.2 g of uranine, were partially saturated with helium, xenon, and argon corresponding to concentrations of 2.7 × 10 −5 mol L −1 (0.01 g), 9.0 × 10 −6 mol L −1 (0.11 g), and 6.4 × 10 −4 mol L −1 (2.2 g). All tracers were injected jointly in Pz2 at a flow rate of 0.57 m 3 hr −1 , followed by a "resting time" period of 0.33 hr, and a pumping period of 19.70 hr in the isolated fracture at the same flow rate than injection. Pumped water is continuously conveyed to the CF-MIMS and a field fluorimeter for continuous analysis. During this experiment, no pumping is operated in Pz1.

Breakthrough Curves Normalization
Measured breakthrough curves are normalized by the integrated concentration, which provides a residence time distribution (RTD) p(t) (s −1 ), as described by Equation 1

Geophysical Research Letters
where C(t) is the concentration measured at time t (s) in the recovery well.

Analysis of Peak Time and Concentration
In the convergent test, the peak times increase, and the normalized peak values decrease with increasing diffusion coefficient (Figures 1a, 3a, 3c, and 3e and Table 1). For example, compared to uranine, the helium breakthrough is 35% more diluted and has a 15% larger peak time (Figure 1a). The peak time of xenon, a heavy gas, is however delayed, compared to the other light gases (Figure 1a) as discussed in section 2.2.
For the push-pull test, the tracer peak times decrease with increasing molecular diffusion coefficients (Figures 1b, 3b, 3d, and 3f and Table 1). For example, the peak time value of uranine (0.76 hr) is delayed of 5% compared to helium (0.72 hr). As in the convergent test, xenon is delayed (11% and 5% compared to helium and uranine, respectively) ( Figure 1b). The trend of the normalized peak values is less clear, and helium and xenon values are unexpectedly high, compared with uranine and argon (Figures 1a, 1b, and 3f and Table 1). Note however that some uncertainty affecting these values may arise from background concentration estimation issues related to remnant tracer from previous experiments. Measured xenon concentration values are also characterized by some visible oscillations, due to strong night and day temperature variations affecting the accuracy of measurements and the assessment of the background for the push-pull experiment.
The different trends observed for peak time and amplitude in convergent and push-pull tests may be understood as follows. For the convergent test, an increasing diffusion coefficient leads to enhanced exchange with the matrix leading to both delayed peak time and smaller peak amplitude (Sanford et al., 2002). For the push-pull test, increasing the diffusion coefficient enhances solute diffusion into the rock matrix and limits transport of the tracer from the injection point to the fracture network. In the pull phase, a larger mass of tracer is close to the injection point, and therefore the peak time decreases with the diffusion coefficient, as shown by Klepikova et al. (2016) for heat diffusion. The peak amplitude decays with diffusion coefficient like in the convergent test, due to enhanced diffusion. In the following, these trends are confirmed from analytical solutions.

Analysis of Late-Time Tailing
Late time tailing could be analyzed with high precision for helium and uranine, while xenon and argon were affected by the resolution limit ( Figure 1c). In the case of xenon, the background measurements are fluctuating around 4 × 10 −10 mol L −1 , which is a too low-resolution limit for accurate late time tailing observations (Text S1 and Figures S2 to S4). As argon is one of the major dissolved gas in waters and the tracer injection amount is limited, it remains difficult to induce a contrast which is large enough ( Figure S2).
Helium and uranine display power law behavior both for the convergent and push-pull tests. The exponent value is calculated considering concentration values where the slope is constant (from 2.5 to 36 hr for the convergent test and from 6 to 16 hr for the push-pull tests). For the convergent test, the tailing of helium has a significantly lower exponent (b ¼ 1.10) compared to uranine (b ¼ 1.53) (Figure 1c). In the push-pull test, the late time tailing of helium and uranine is similar (Figure 1d) and characterized by an exponent approximately equal to b ¼ 1.5 (Table 1).
The exponent b ¼ 1.5 is typical of diffusion perpendicular to a planar fracture into the adjacent rock matrix (Haggerty et al., 2000). So all measured late time tailings except for helium in the convergent test suggest the relevance of a simple model of fracture-matrix exchange. However, the convergent test tailing of helium, the tracer with the largest diffusion coefficient, is significantly different (b ¼ 1.10). Klepikova et al. (2016) showed that for highly channelized flow, diffusion from fractures to matrix becomes radial instead of effectively one-dimensional for the parallel plate model. This leads to a late time tailing as c(t)~t −1 /log(t 2 ), which is close to c(t)~t −1 . Hence, the late time tailing of helium suggests that tracers that penetrate more deeply into the matrix are influenced by a radial diffusion behavior, possibly indicating strong channeling of flow in fractures. Over the experimental time (39 hr), the characteristic scale of penetration of uranine into the matrix is 2 cm, while that of helium is 7 cm. In the investigated aquifer, fractures are regionally observed in outcropping chalks, including numerous bedding joints, each of them being characterized by a highly variable aperture, inducing channelization of flow. Based on that, it is hypothesized that for small matrix penetration depths, solutes only probe the fringes bordering the channels. At this scale, the effect of larger-scale channels may be negligible and mass transfer may be effectively perpendicular to the local  (Table 1).

Geophysical Research Letters
fracture-rock interface (Figure 2). At larger penetration depths, it is likely that the effect of radial diffusion from channels starts to influence fracture-matrix mass transfer ( Figure 2).
The change of tailing exponent with tracer diffusion is not observed for the push-pull tests. This may be due to the smaller involved rock volumes, time intervals, and distances in the push-pull experiment. The characteristic scale of the convergent test may be estimated as the distance between the boreholes d ¼ 7.55 m. It is significantly smaller for the push-pull test where it is estimated as the average radius of injected fluid around the injection point, r ¼ (3V/4φπ) 1/3 , with V the total injected volume. This radius remains below 1.4 m if considering a minimum porosity φ of 1% (Table 1). Note however that the uncertainty related to the tailing exponent calculation is higher for this push-pull experiment than for the convergent test. This is mainly due to the combined challenge related to the duration of the monitoring after the peak time and the correct evaluation of the background concentration, which was here influenced by previous experiments and not strictly constant during the push-pull monitoring, affecting the normalization and slope calculation procedures.

Breakthrough Curve Modeling
The analysis of peak time and amplitude, as well as late time tailing provides some qualitative insights into the transport dynamics and fracture-matrix interactions. To extract quantitative information on these dynamics, we interpreted the breakthrough curves with analytical solutions of advection and dispersion in fractures and diffusion into the matrix. Following the observations of breakthrough curve tailing, two different conceptual models for the fracture geometry are tested: the parallel plate model and the channel model ( Figure S5). The analytical solutions described by  and  based on Laplace transforms (Text S2) are used. These Laplace transform solutions are inverted using the algorithms from Hoog et al. (1982) and Hollenbeck (1998). Manual parameter adjustments are simultaneously operated on the experimental uranine and helium breakthrough curves in RTD using Equation 1 and considering equal weight on the peak value, peak time, and slope in RTD and minimizing the difference between observations and simulations according to the diffusion coefficient ( Figure 3).
Optimized parameters include the number of fractures or channels, the fracture aperture or channel radius, and the dispersivity in fractures/channels.
As expected from the different power law tailing exponents of the uranine and helium breakthrough curves in the convergent tracer test (Figure 1c), the uranine breakthrough curve is well captured by the parallel plate model, while the helium breakthrough curve is well captured by the channel model. The uranine breakthrough curve is simulated using n f ¼ 40 fractures with an aperture a ¼ 1 mm and a dispersivity of α ¼ 2.4 m. The helium recovery is simulated using n c ¼ 700 channels characterized by a radius of a ¼ 1 mm and a dispersivity of α ¼ 1.5 m. This suggests that, in the area sampled by the tracer, fractures contained on average 17 channels per fracture. Using the parallel plate model for the uranine breakthrough curve, fails to predict the helium breakthrough curve and vice versa (Text S4 and Figure S10), which confirms the transition in the fracture-matrix interaction mode with different tracer penetration depths into the matrix. This transition tends to be confirmed by the modeling of peak

Geophysical Research Letters
time and amplitude, which shows a transition from the parallel plate trend to the channel trend (Figures 3c and 3e). As discussed in Text S3 and in Figures S6 to S9, the dependency of the peak time and amplitude on the diffusion coefficient shows a strong sensitivity to model parameters. Hence, the analysis of the dependency of peak time and amplitude on the diffusion coefficient offers an alternative way to characterize the dynamics of fracture-matrix interactions.
Push-pull experiments are generally consistent with the parallel plate fracture model, which successfully captures the uranine and helium breakthrough curves with a single set of parameters (Figure 3b). Note that the slight mismatch for the tailing of uranine may be due to more uncertain background concentrations  (Table 1).

10.1029/2020GL088944
Geophysical Research Letters related to previous convergent experiment. The best fit aperture is a ¼ 1 mm, as for the convergent test. The number of fractures for the best fit model was 4, that is, much less than for the convergent test, which reflects the much smaller rock volume investigated by the push-pull test. The best fit dispersivity was α ¼ 0.1 m, which is much lower than for the convergent test. This is consistent with the fact that a large part of dispersion, the reversible dispersion, is canceled in push-pull tests (Kang et al., 2015). The model also captures the dependency of the peak time and peak value to the diffusion coefficient, although the trend is less clear for the peak (Figures 3d and 3f).

Conclusions
In this study, we have explored the use of continuously monitored dissolved gases as tracers of fracture-matrix interactions using recently developed field mass spectrometers. Highly accurate measurements at low concentrations allowed for a refined analysis of late time breakthrough curve tailing, a key for capturing the physics of fracture-matrix interactions and bringing new insights about the debate on power law decay. The experimental results suggested that the model derived from a conventional fluorescent dye, with relatively low matrix penetration depth, could not be used to predict the transport behavior of helium, characterized by a deeper penetration depth into the matrix. The study showed that the effects of channeled fracture flow started to be dominant at larger scales, highlighting the challenge to describe the scale dependency of fracture-matrix interactions and the interest of using tracers of different diffusion coefficients for conceptual model definition.
Tracer tests results bring complimentary data supporting the reconstruction of the porous/fractured medium. Interpretation with analytical solutions showed that these innovative data are sensitive to hydraulic and geometrical parameters, which are also better constrained in comparison with classical experiments. Results thus support characterization and further modeling of dual media toward more robust simulations and predictions.