# Enhancing Mixing During Groundwater Remediation via Engineered Injection-Extraction: The Issue of Connectivity

Oriol Bertran, D Fernàndez-garcia, Guillem Sole-mari, P Rodríguez-escales

► To cite this version:

Oriol Bertran, D Fernàndez-garcia, Guillem Sole-mari, P Rodríguez-escales. Enhancing Mixing During Groundwater Remediation via Engineered Injection-Extraction: The Issue of Connectivity. Water Resources Research, 2023, 59 (7), pp.e2023WR034934. 10.1029/2023wr034934. insu-04158302

HAL Id: insu-04158302

https://insu.hal.science/insu-04158302

Submitted on 11 Jul 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Enhancing Mixing During Groundwater Remediation via Engineered Injection-Extraction: The Issue of Connectivity

O. Bertran<fnref n="1" />,<fnref n="2" />, D. Fernàndez-Garcia<fnref n="1" />,<fnref n="2" />, G. Sole-Mari<fnref n="3" />, P. Rodríguez-Escales<fnref n="1" />,<fnref n="2" />

## Key Points:

* The spreading and dilution of a treatment solution is tested under rotating-dipole Engineered Injection-Extraction (EIE)

* EIE is similarly effective in both Multigaussian log-conductivity fields and well-connected non-Multigaussian fields

* EIE reduces the remediation outcome uncertainty associated to porous medium heterogeneity and connectivity

<footnote n="1">¹Department of Civil and Environmental Engineering (DECA), Universitat Politècnica de Catalunya (UPC), Jordi Girona 1-3, 08034 Barcelona, Spain</footnote>

<footnote n="2">²Associated Unit: Hydrogeology Group (UPC-CSIC)</footnote>

<footnote n="3">³Géosciences Rennes, Université de Rennes, France</footnote>
**Abstract**

In the context of in-situ groundwater remediation, mixing is vital for a successful out-come. A slow mixing rate between the contaminated groundwater and the injected treat-ment solution can severely weaken the effective degradation rate. Engineered Injection-Extraction (EIE) has been proposed as a means to accelerate dilution within the porous medium. However, existing studies on the subject have not considered the potential im-pact of connectivity and preferential flow-paths. Neglecting connectivity can lead to an overestimation of EIE's capabilities, since the fluid may in reality be carried mainly through a few high-permeability channels, thus hampering mixing and reaction. Due to the fact that channeling can be found in many actual sites, in this work we aim to evaluate EIE methods in both poorly-connected (represented as Multigaussian fields) and well-connected fields (represented as non-Multigaussians). The approach is to identify, for each given medium, a stirring protocol – defined by a specific combination of rotation angle and ro-tation rate – which maximizes mixing. To that end, metrics are proposed in order to (1) quantify both the mixing and the containment of the treatment solution within a given remediation volume, and (2) characterize the particle trajectories to explicitly evaluate if preferential paths are broken. The results obtained from these metrics are quite sim-ilar for both types of fields, proving that the enhancing of mixing by means of EIE is ef-fective regardless of the presence of preferential flow paths. This study demonstrates that EIE via rotating dipoles diminishes the remediation outcome uncertainty induced by medium heterogeneity.

# 1 Introduction

Most *in-situ* groundwater remediation schemes are based on the injection of a treat-ment solution into the subsurface in order to promote reactions that remove or trans-form the contaminants into harmless products. For instance, the success of both biore-mediation and *in-situ* chemical oxidation is based on the mixing between contaminated groundwater and nutrients or oxidant agents, which can promote the occurrence of tar-get reactions and facilitate pollution removal. These systems can be enhanced by favor-ing increased mixing, that is, by promoting mass transfer between the different solutions. The physical processes that boost mass transfer in porous media are mechanical disper-sion and molecular diffusion (Tartakovsky, 2010). The former spreads a solute through-out the porous medium, increasing contact surface and concentration gradients, and the latter homogenizes it leading the solute to occupy a larger volume in the porous medium. However, these processes are typically slow in natural conditions, and thus are often un-able to provide sufficient mixing – and hence degradation.

Mixing in porous media has received much attention in the last few decades, largely because chemical reactions in the subsurface are often limited by mixing (Sturman et al., 1995; Cirpka et al., 1999; Cirpka, 2002; de Simoni et al., 2005; Rolle et al., 2008). In the context of groundwater remediation, some authors have proposed Engineered Injection-Extraction (EIE) as a method to promote mixing and, consequently, reactions in porous media (e.g., Zhang et al., 2009; Lester et al., 2010; Trefry et al., 2012; Piscopo et al., 2013). EIE involves a sequence of injections and extractions through wells with time-dependent flow rates, resulting in transient velocity fields which stretch and fold the solute plume (Bagtzoglou & Oates, 2007). It has been demonstrated that EIE can generate chaotic advection in porous media (Metcalfe et al., 2008; Lester et al., 2009, 2010; Trefry et al., 2012). This phenomenon, first coined in such terms by Aref (1984), is typically charac-terized by particle trajectories which experience exponential stretching over time, result-ing in complex patterns and a virtually unpredictable evolution from the initial condi-tion (Lester et al., 2018; Turuban et al., 2019). EIE schemes increase the active spread-ing and the area occupied by the injected fluid in aquifers, favoring mixing and accel-erating the natural attenuation and degradation of contaminants. The benefits of this strategy has been proven for many configurations of extraction/injection wells through

Accepted Article
numerical simulations (e.g., Piscopo et al., 2013; Neupauer et al., 2014; Di Dato et al., 2018; Speetjens et al., 2021), laboratory experiments (e.g., Zhang et al., 2009; Sather et al., 2023) and field applications (Cho et al., 2019). Most recently, EIE has also been demonstrated to enhance NAPL recovery in multiphase flow systems (Wang et al., 2022).

Several schemes can be used to create time-dependent velocity fields via EIE. Bagtzoglou and Oates (2007) and Zhang et al. (2009), used an oscillatory pumping scheme consisting of three randomly located wells with some realistic constraints to avoid dewatering and preserve mass. Mays and Neupauer (2012), Piscopo et al. (2013) and Neupauer et al. (2014), employed a four-well configuration operating sequentially, alternating between injection and extraction and only using one well at a time. Using a purely mathematical approach, Lester et al. (2009) studied a temporally rotating 2D dipole scheme. Although recirculation of the treatment solution from the extraction to the injection well can exist in EIE systems, most of the previous literature has not addressed this issue. Recirculation can also enhance mixing and thus the effectiveness of the *in situ* remediation system (Cirpka & Kitanidis, 2001), since it extends the residence time of the treatment solution within the reactive zone (Luo et al., 2007). However, frequent recirculation could generate bioclogging near the injection wells during bioremediation schemes (McCarty et al., 1998; Cirpka et al., 1999; MacDonald et al., 1999a, 1999b; Gandhi et al., 2002). Yet, little is known about the role of recirculation in EIE applications.

Mixing aside, solute transport in porous media is highly dependent on the spatial variability of hydraulic conductivity (heterogeneity) and also the presence of preferential flow paths. The observations of anomalous transport and scale effects can often be associated with the omission of connectivity features (Boggs & Adams, 1992; Schulze-Makuch & Cherkauer, 1998). Connectivity in hydrogeology refers to the prevalence of higher-conductivity (K) paths that can carry and accelerate the transport of solutes from one region to another, largely bypassing the lower-conductivity areas (Knudby & Carrera, 2005; Renard & Allard, 2013). High connectivity is thus a true hindrance in the application of a remediation scheme since it can dominate the transport of the injected solution (McGregor & Benevenuto, 2021). Yet, specific information about medium connectivity is usually omitted due to: (1) the lack of exhaustive geological/hydrogeological data and (2) the selection of inadequate models to characterize the heterogeneity of the porous media.

Regarding the latter point, the most common approach for characterizing the spatial variability of the log-conductivity fields is the use of two-point – also called variogram-based – geostatistics (Hashemi et al., 2014) and Multigaussian random fields, which are characterized by their mean, variance and covariance function. These multigaussian models are mathematically simple and easy to treat analytically, yet they cannot reproduce high-connectivity patterns. Some non-Multigaussian models, on the other hand, are capable of emulating high connectivity by arranging the high log-K values as channels (Gmez-Hernndez & Wen, 1998). Although several studies conclude that EIE can enhance mixing in Multigaussian heterogeneous porous media (Rodrguez-Escales et al., 2017), the applicability of EIE in real settings featuring strong preferential flow paths is still not clear.

Remediation efficiency depends on the degree of mixing of the treatment solution within the target treatment zone, typically defined based on the areal extent of contamination. In this context, transport of this solution beyond the treatment zone represents a waste of resources which should be taken into consideration. Surprisingly, dilution metrics for permeable control volumes have not been proposed in the literature of solute transport in porous media, which has essentially focused on infinite domains (Kitanidis, 1994).

This work is aimed at studying how EIE systems based on rotating dipoles (Metcalfe et al., 2008) enhance dilution in a remediation scenario and help break the preferential flow paths that often exist in highly heterogeneous aquifers. For that, we evaluate how

Accepted Article
different EIE configurations improve (1) the dilution of an injected solution in an active remediation volume, and (2) the ability of treatment solution particles to fully explore the medium. Through a set of Monte Carlo simulations we explore different EIE configurations within both Multigaussian and well-connected non-Multigaussian log-conductivity fields.

## 2 The Engineered Injection-Extraction Method

We study an Engineered Injection-Extraction (EIE) setup based on a rotating dipole, emulating the idealized Rotated Potential Mixing (RPM) system (e.g. Metcalfe et al., 2008; Lester et al., 2009; Trefry et al., 2012). The procedure's goal is to enhance the mixing, within a circular treatment zone of radius $R_r$, of a treatment solution initially delivered at its center $(x_0, y_0)$. The treatment zone encompasses the contaminated region of the aquifer, and the determination of its size in relation to the positioning of wells is decided case-by-case attending the properties of the media and the distribution of pollutants. As a general rule, its radius ($R_r$) is larger than the radius of the dipoles ($R$) itself to prevent mass for immediately leaving the treatment area upon reinjection.

At $t=0$, the treatment solution is assumed to be uniformly distributed within a concentric circular region with radius $r_i < R_r$. We wish that the injected treatment solution: (1) effectively mixes with the resident groundwater; and (2) stays within the treatment zone during remediation. That is to say that we wish to transform the treatment zone into a well-mixed reactor. The treatment solution is stirred by a 2D dipole system that consists of two active wells (one injector and one extractor) separated by a distance $2R$ that operate at the same time with a constant flow rate $Q$. The injection and extraction wells are symmetrically placed around the center of the domain and periodically re-oriented (emulating a larger circular array where only two wells are active at a given time). The pumped water from the extraction well is recirculated instantaneously into the injection well. Figure 1(a) shows a sketch diagram of the setup. The parameters adopted are summarized in Table 1.

During EIE operations, the dipole wells are periodically activated for a rotation period $\tau$, after which the active dipole position rotates counter-clockwise around the origin by an angle $\Theta$ and operates again for a time $\tau$, and so on. We define the characteristic dipole-zone emptying time $t_c$ as the time needed to renew the entire volume of water occupied within the wells zone at a constant flow rate, that is, $t_c = \pi R^2 b \phi / Q$, where $\phi$ is the porosity and $b$ is the aquifer thickness. The dimensionless time is then defined in this work as $t^* = t/t_c$. A range of values for the rotation angle $\Theta$ (11 values from $\pi/5$ to $\pi$) and the rotation period $\tau$ (25 values from $\tau=0.05t_c$ to $\tau=10t_c$) are tested, resulting in a total of 275 combinations (stirring protocols). In the following, the rotation period and the angle of rotation will be presented in normalized form, defined as $\tau^* = \tau/t_c$ and $\Theta^* = \Theta/\pi$.

## 3 Flow and transport simulations

We solve flow and transport in a two-dimensional confined square aquifer of size $L$. The control volume and the treatment solution are located at the center of the domain, $x_0 = L/2$ and $y_0 = L/2$. In all outer boundaries a no-flow condition is imposed. Regional flow is neglected for simplicity. The hydraulic conductivity is assumed locally isotropic but spatially heterogeneous. To simplify the system, we consider that dipole flows are time-dependent but piecewise steady. Consequently, groundwater flow is described by the following equation,

$$ \nabla \cdot [\mathbf{K}(\mathbf{x})b\nabla h(\mathbf{x}, t)] + Q\delta(\mathbf{x} - \mathbf{x}^+(t)) - Q\delta(\mathbf{x} - \mathbf{x}^-(t)) = 0, \quad (1) $$

-4-

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

Accepted Article
<figure><img src="image_3.png" /><figcaption>Figure 1: (a) Conceptualization of the model setup. The left figure shows 6 wells forming 3 dipoles separated by a radius R. The circular treatment area is determined by a radius R<sub>r</sub>, and r<sub>i</sub> delineates the area of injection. Periodically, every time lapse τ, the active EIE dipole rotates by an angle Θ; (b) Three representative scenarios are shown, each one being the result of a different combination of Θ and τ. Two of them, the scenarios 1 and 3, display undesirable outcomes of insufficient dilution and ineffective containment, respectively. The scenario 2, in green, represents the desired outcome, featuring: (i) high dilution, and (ii) effective containment of the solute within the reactor.</figcaption></figure>

<table><caption>Table 1: Parameters that define the conceptual EIE model</caption><thead><tr><th>Parameters</th><th></th><th>Values</th></tr></thead><tbody><tr><td>Treatment solution radius</td><td>r<sub>i</sub></td><td>10 m</td></tr><tr><td>Dipole radius</td><td>R</td><td>25 m</td></tr><tr><td>Reactor radius</td><td>R<sub>r</sub></td><td>37.5 m</td></tr><tr><td>Aquifer thickness</td><td>b</td><td>10 m</td></tr><tr><td>Aquifer length</td><td>L</td><td>300.25 m</td></tr><tr><td>Characteristic time</td><td>t<sub>c</sub></td><td>10 d</td></tr><tr><td>Porosity</td><td>φ</td><td>0.25</td></tr><tr><td>Injection-Extraction rate</td><td>Q</td><td>491 m<sup>3</sup>/d</td></tr><tr><td>Amount of treatment solution particles</td><td></td><td>30000</td></tr><tr><td>Normalized rotation angle</td><td>Θ*</td><td>{<sup>1</sup>&frasl;<sub>5</sub> ; <sup>1</sup>&frasl;<sub>4</sub> ; <sup>1</sup>&frasl;<sub>3</sub> ; <sup>2</sup>&frasl;<sub>5</sub> ; <sup>1</sup>&frasl;<sub>2</sub> ; <sup>3</sup>&frasl;<sub>5</sub> ; <sup>2</sup>&frasl;<sub>3</sub> ; <sup>3</sup>&frasl;<sub>4</sub> ; <sup>4</sup>&frasl;<sub>5</sub> ; <sup>9</sup>&frasl;<sub>10</sub> ; 1}</td></tr><tr><td>Normalized rotation period</td><td>τ*</td><td>{0.05 ; 0.1 ; 0.15 ; 0.2 ; 0.3 ; 0.4 ; 0.6 ; 0.8 ; 1 ;<br>1.2 ; 1.4 ; 1.6 ; 1.8 ; 2 ; 2.4 ; 2.8 ; 3.2 ; 3.6 ;<br>4 ; 4.8 ; 5.6 ; 6.4 ; 7.2 ; 8 ; 10}</td></tr></tbody></table>
where $K(\mathbf{x})$ is the hydraulic conductivity at the $\mathbf{x}$ location, $b$ is the aquifer thickness (assumed constant), $h(\mathbf{x},t)$ is the hydraulic head, $\delta(\mathbf{x})$ is the Dirac delta function, and $\mathbf{x}^+$ and $\mathbf{x}^-$ are the injection and extraction positions of the rotating dipole,

$$x^{\pm}(t) = x_0 \pm R \cos \left( f \left( \frac{t}{\tau} \right) \Theta \right), \quad y^{\pm}(t) = y_0 \pm R \sin \left( f \left( \frac{t}{\tau} \right) \Theta \right), \quad (2)$$

where $\Theta$ is the rotation angle, $\tau$ is the rotation period and $f$ is the floor function defined as $f(x) = \max\{n \in \mathbb{Z}|n \le x\}$. The transport of the treatment solution is described by the Advection-Dispersion equation (ADE):

$$\phi \frac{\partial c}{\partial t}(\mathbf{x}, t) = -\nabla \cdot (\mathbf{q}(\mathbf{x}, t)c(\mathbf{x}, t)) + \nabla \cdot (\phi \mathbf{D}(\mathbf{x}, t)\nabla c(\mathbf{x}, t)), \quad (3)$$

where $c(\mathbf{x}, t)$ is the solute concentration, $\mathbf{q}(\mathbf{x}, t)$ is the Darcy velocity, $\phi$ is the porosity (assumed constant), and $\mathbf{D}(\mathbf{x}, t)$ is the local dispersion tensor. For simplicity and since mixing is typically controlled by local transverse (rather than longitudinal) dispersion, we assume that $\mathbf{D}(\mathbf{x}, t)$ is locally isotropic,

$$\mathbf{D}(\mathbf{x}, t) = \alpha \|\mathbf{v}(\mathbf{x}, t)\| \mathbf{I}_d, \quad (4)$$

where $\alpha$ is the constant dispersivity, $\|\mathbf{v}(\mathbf{x}, t)\|$ is the norm of the fluid velocity vector at any $\mathbf{x}$ location and time $t$ (with $\mathbf{v}(\mathbf{x}, t) = \mathbf{q}(\mathbf{x}, t)/\phi$), and $\mathbf{I}_d$ is the identity matrix. A value of $\alpha = 0.25$ m is set. This yields a Péclet number (defined here as the ratio between dipole radius and dispersivity) of 1000.

To solve Eq. (1) we use the finite-difference code MODFLOW 2005 (Harbaugh et al., 2017), with the Preconditioned Conjugate-Gradient method (PCG). The model domain is discretized into 1201 × 1201 square cells of size $\Delta = 0.25$ m. The resulting cell-interface fluxes are then used to solve the advection-dispersion equation (ADE) by means of the Random-Walk Particle Tracking method (RWPT). The RWPT simulations are performed by initially distributing 30000 particles uniformly within the circular treatment solution. The results have been verified to have converged for this number of particles. The RWPT method used here is based on the approach of LaBolle et al. (2000), which is suitable for simulating transport through porous media in systems with abrupt changes in dispersion. The formal derivation and details of this RWPT method are not presented here; the interested reader may refer to the source itself or other works such as Delay et al. (2005), Salamon et al. (2006) or Sole-Mari et al. (2021).

The scattered data points of the particle positions are then reconstructed into concentrations using an open-source MATLAB code developed by Sole-Mari et al. (2019). With this approach, concentrations are efficiently reconstructed by combining histogram methods with locally adaptive Kernel Density Estimation methods. This makes the estimation of concentration efficient, robust and less restrictive on the number of particles.

## 4 Stochastic approach

Monte Carlo simulations are performed to represent multiple equiprobable realities of the natural logarithm of the hydraulic conductivity field $Y(\mathbf{x}) = \ln K(\mathbf{x})$, which is described as a random function of space. We generate 50 Multigaussian and 50 non-Multigaussian random fields (Fig. 2). The difference between Multigaussian (MG) and the particular type of non-Multigaussian (nMG) fields is that the former have the extreme values of $Y(\mathbf{x})$ isolated from the rest and the latter have the higher values connected. Both MG and nMG fields share similar two-point statistical properties: (1) identical log-normal hydraulic conductivity distribution and (2) near-identical isotropic spatial covariance function.

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

Accepted Article
<figure><img src="image_4.png" /><figcaption>Figure 2: Central region of two selected Multigaussian and connected non-Multigaussian fields used in the simulations. the wells zone is represented with a solid line, and the treatment zone with a dashed line. Below, the corresponding histograms of $Y = \ln K$.</figcaption></figure>

MG fields are generated using the SGeMS software (Remy, 2005) via the sequential Gaussian simulation subroutine SGSIM with zero mean, variance $\sigma_Y^2 = 2$, and an isotropic Gaussian covariance function model with a range of $a=5$ meters, which is 10% of the distance between injection and extraction wells in a dipole. This relatively low range results in a somewhat poor hydraulic connection between the injection and extraction wells. nMG fields, on the other hand, are specifically generated to display high connectivity of high values that favor the formation of preferential flow paths. We use the method developed by Zinn and Harvey (2003), which transforms any MG field into a nMG field by (1) normal scoring the MG field; (2) taking the absolute value; and (3) transforming the distribution of absolute values into a log-normal distribution with zero mean and $\sigma_Y^2$ variance:

$$Y' = -\sigma_Y \sqrt{2 \operatorname{erf}^{-1} \left(2 \operatorname{erf} \left(\frac{Y}{\sqrt{2}}\right) - 1\right)}, \quad (5)$$

where $Y'$ are the transformed values and $Y$ are the absolute values obtained in the previous step. This rearrangement reduces the integral scale of the random field by a factor of 1.86 (yellow line in Fig. 3). This is addressed by resizing the field by a factor of 1.86 and then cropping it also by a factor of 1.86 in order to keep the same domain size for both random field types.

203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218

Accepted Article
<figure><img src="image_5.png" /><figcaption>Figure 3: The different spatial semivariograms obtained during the transformation procedure from one selected MG to the nMG field</figcaption></figure>

# 5 Performance metrics

## 5.1 Volume-control reactor efficiency

We aim to quantify the degree of mixing within the treatment zone, a fixed control volume V defined previously based on remediation targets and the areal extent of contamination (Rr in Figure 1(a)). In order to reduce remediation costs and risks, the treatment solution should only spread and mix within the contaminated area. In this context, it is convenient to have a metric that quantifies the degree of mixing within the remediation target volume. Kitanidis (1994) defines the reactor ratio M as the ratio of the actual dilution index E of a solute plume in a volume V to its maximum theoretical value Emax,

$$M = \frac{E}{E_{max}}. \qquad (6)$$

The dilution index is defined as

$$E(t) = \exp \left[ - \int_V p(\mathbf{x}, t) \ln(p(\mathbf{x}, t)) dV \right], \qquad (7)$$

where p(x, t) is the probability density for a mass particle to be located on x at time t,

$$p(\mathbf{x}, t) = \frac{\phi c(\mathbf{x}, t)}{\int_V \phi c(\mathbf{x}, t) d\mathbf{x}}. \qquad (8)$$

The dilution index E is a well-established metric of dilution that quantifies the volume occupied by the solute plume in V (Kitanidis, 1994). Note for instance that when the solute is uniformly distributed in the aquifer over a volume W (W < V), assuming constant porosity, p approaches 1/W and E tends to W, which is the volume occupied by the solute. The reactor ratio M can thus be interpreted as the fraction of the volume V occupied by the solute plume, i.e., in this example we have that M = W/V. Thus, M is a dimensionless number that ranges between 0 and 1. The larger the value of M, the closer the system is to well-mixed conditions (complete mixing). Also note that, in our case, maximum dilution in V is achieved when

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239
$$E = E_{max} = V = \pi R_{r}^{2} b \phi. \qquad (9)$$

Kitanidis (1994) defines the reactor ratio M for two specific cases only: a bounded and an unbounded system. Neither of these cases accurately describes our system. The former would imply the assumption of a closed system where the injected solution cannot leave the control volume V. The latter describes the state of dilution of the plume as it is transported within an infinite aquifer. In our case the target volume is not isolated with impermeable boundaries nor is it infinite, but is a finite permeable volume defined by remediation targets. Depending on the pumping and injection sequence used, it is possible for a portion of the injected treatment solution to escape from the control volume V, and that portion should not count towards the well-mixed volume metric, since our goal is to maximize mixing within V. Hence we define the volume-control reactor efficiency as

$$\varepsilon_V = \omega M. \qquad (10)$$

Here, $\omega$ is the mass fraction of the injected treatment solution that remains in V at all times, i.e., $\omega = m_{in}/m_{tot}$, where $m_{in}$ is the mass of injected solution remaining inside of V, and $m_{tot}$ is the total injected mass. $\varepsilon_V$ can also be seen as a measure of the probability that the injected treatment solution is completely diluted in the control volume V. It is therefore a dimensionless number that ranges between 0 and 1. A value close to 0 indicates that mixing is not complete inside of V and/or is taking place outside, and a value close to 1 indicates that the total injected solution has perfectly mixed within the control volume V (see Scenario 2 in Fig. 1(b)).

## 5.2 Lagrangian semivariogram

An important objective of this work is to determine whether EIE can partially overcome the detrimental effect of preferential channels on solute mixing. In this context, the dilution index might not be sufficient in order to completely describe the interaction between solute particles and the porous medium, since it does not provide any information about individual particle trajectories but rather a snapshot of their global positions. Hence we also examine the Lagrangian semivariogram of $Y(\mathbf{X}_p(t))$ values visited by particles along their paths,

$$\gamma(t') = \frac{1}{2} E \left[ \left( Y(\mathbf{X}_p(t+t')) - Y(\mathbf{X}_p(t)) \right)^2 \right], \qquad (11)$$

where E is the expected value operator, and $t'$ is the time lag between any two particle positions. This Lagrangian semivariogram measures the average degree of dissimilarity (in terms of log-conductivity) between the zones visited by the particles.

Note, however, that during EIE the recirculation of particles from extraction to injection well can produce artificial jumps that do not reflect the temporal persistence of particles moving through similar conductivity zones. To avoid this, the Lagrangian semivariogram is computed by segmenting the particle trajectories into different paths. A particle path is defined as the concatenation of positions visited by a particle without being recirculated. Figure 4 illustrates this with an example of the Y values sampled by one particle and its corresponding displacement. For each realization and stirring protocol, the sample semivariogram is calculated as

$$\gamma(t') \approx \frac{1}{2N(t')} \sum_{\substack{(i,j)|t'=|t_j-t_i|, \\ (i,j) \in path}} (Y_{t_i} - Y_{t_j})^2, \qquad (12)$$

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

Accepted Article
<figure><img src="image_6.png" /><figcaption>Figure 4: Values of lnK sampled by a selected particle, and its discrete displacements during the simulation, versus cumulative dimensionless time $t^* = t/t_c$. Whenever a jump from the extraction to the injection well occurs, a new path is generated.</figcaption></figure>

where $Y_{t_i}$ is the Y value at the $\mathbf{X}_p(t_i)$ position, and $N(t')$ is the number of pairs $(Y_{t_i}, Y_{t_j})$ such that $t' = |t_j - t_i|$, and $\mathbf{X}_p(t_i)$ and $\mathbf{X}_p(t_j)$ belong to the same path and particle trajectory.

In contrast to the dilution index, which is based on a snapshot at a given time, the Lagrangian semivariogram is an integrated metric which encapsulates the movement of particles during the entire simulation, giving valuable information about the hydrological response to connectivity.

The sample Lagrangian semivariogram of Y approaches a stable plateau at increasing time lags (shown later on in section 6), suggesting that, over the temporal and spatial scales of this problem, we deal with a stationary regionalized variable. Based on this, we characterize the Lagrangian semivariogram of Y by its sill ($\sigma^2$) and integral scale ($I_t$). The sill is estimated by the asymptote of the semivariogram and represents the range of Y values sampled by the particles. When the most frequented locations within paths have similar Y values, we have that $\sigma^2 \ll \sigma_Y^2$. This can occur for instance when particles concentrate their movement through preferential channels. When, on the contrary, the random field is fully sampled (complete mixing), we have that $\sigma^2 \approx \sigma_Y^2$. The integral time scale $I_t$ quantifies the temporal persistence of Y values along particle paths. The greater its value, the longer the particles tend to remain in the same conductivity zone during the simulation. The integral time-scale is defined as

$$I_t = \frac{1}{\sigma^2} \int_0^\infty C(t') dt', \qquad (13)$$

where $C(t')$ is the Lagrangian covariance function of $Y(\mathbf{X}_p(t))$. Since, for a stationary random process, $C(t') = \sigma^2 - \gamma(t')$, the integral time scale can be calculated from the sample Lagrangian semivariogram by numerically solving the following limit,

$$I_t \approx \lim_{t_f \to \infty} \left[ t_f - \frac{1}{\sigma^2} \int_0^{t_f} \gamma(t') dt' \right]. \qquad (14)$$

In practice, this limit is attained for $t_f$ larger than the time lag at which the correlation vanishes.

278

279

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

Accepted Article
During EIE, we wish to maximize the sill $\sigma^2$ (particles' proneness to move through different zones) and reduce the integral time scale $I_t$ (typical residence time within zones of similar conductivity). To quantify this in a synthesized manner, we define the following dimensionless metric of domain sampling effectiveness,

$$ \Omega = \left( \frac{\sigma^2}{\sigma_Y^2} \right) \left( \frac{t_c I_Y}{R I_t} \right), \qquad (15) $$

where $I_Y$ is the integral scale of the porous medium, $I_Y = a\sqrt{\pi}/6$, and $a$ is the spatial range of $Y(\mathbf{x})$.

In order to characterize treatment solution recirculation during EIE, We also measure the frequency of recirculation, defined as the average number of instances per unit of time that particles recirculate through wells during operation.

# 6 Results and discussion

## 6.1 Reactor efficiency

We start the analysis of the results by describing the temporal evolution of the volume-control reactor efficiency $\varepsilon_V$ during EIE. Figure 5 shows the ensemble average of $\varepsilon_V$ over all realizations of (a) MG and (b) nMG random fields as a function of time for all stirring protocols (combinations of $\Theta^*$ and $\tau^*$). Results show the existence of a set of optimal EIE pumping sequences which achieve the highest $\varepsilon_V$ values with a nearly stable plateau after a few characteristic times. These optimal protocols are highlighted with dotted lines. On the other hand, $\varepsilon_V$ exhibits a weak or slow response, sometimes followed by a late-time decline when a sub-optimal protocol is used.

<figure><img src="image_7.png" /><figcaption>Figure 5: Temporal evolution of the volume-control reactor efficiency $\varepsilon_V$ for every combination of $\Theta^*$ and $\tau^*$, within (a) MG fields, and (b) nMG fields. Optimal stirring protocols are highlighted with dotted lines.</figcaption></figure>

A similar temporal evolution of $\varepsilon_V$ is observed in both MG and nMG random fields, which are found to share the same optimal sequences of $\Theta^*$ and $\tau^*$. However, nMG fields exhibit relatively smaller optimal reactor efficiencies ($\varepsilon_V$ decreases by about 10%). This

299

300

301

302

303

304

305

306

307

308

309

310

311

312

313

314

315

316

317

318

319

320

Accepted Article
<figure><img src="image_8.png" /><figcaption>Figure 6: Ensemble average of the temporal mean over $t^* \in [0, 20]$ of (a) the Reactor Ratio of Kitanidis $M$, (b) the mass fraction inside the treatment zone $\omega$, and (c) the volume-control reactor efficiency $\varepsilon_V = \omega M$ as a function of $\Theta^*$ and $\tau^*$, for MG and nMG fields. The dashed line in the Reactor Efficiency figures follows the maximum values of dilution found for each $\Theta^*$.</figcaption></figure>

difference in efficiency between MG and nMG fields indicates that even though EIE can significantly enhance mixing, preferential flow-paths do exert some negative impact on its effectiveness in nMG fields. Figure 6, displays maps of the ensemble average of (a) $M$, (b) $\omega$ and (c) $\varepsilon_V = \omega M$, averaged over $t^* \in [0, 20]$, as a function of $\Theta^*$ and $\tau^*$, in both MG and nMG fields. These results indicate that the connectivity of nMG fields cause reduced dilution ($M$), as well as a slightly lower ability to contain the treatment solution ($\omega$) within the designated area. For a given value of $\Theta^*$, an exceedingly slow rotation period, $\tau^* \ll 1$, will not stir the injected solution properly and, as a consequence, much of the solution will stay near the injection location and will not mix substantially. On the other hand, if the rotation period is too high, $\tau^* \gg 1$, a large amount of the solution will leave the treatment zone. This is illustrated in Figure 7, where the particle positions at time $t^* = 8$ are shown for $\tau^* \in [0.1, 1, 10]$, and two values of $\Theta^*$. We finally note from Figure 6 that, generally, the higher the rotation angle, the higher the corresponding optimal rotation period. Actually, for $\Theta^* < 3/5$ the optimal rotation period can be estimated as $\tau^* \approx 2\Theta^*$. A singular behavior occurs at $\Theta^* = 2/3$, where the volume-control reactor efficiency is rather low and quite insensitive to the rotation period. The explanation for this singular behavior is not trivial, but it is related to the Hamiltonian of the flow field for $\tau \to 0$ and the associated mode locking, as shown by Lester et al. (2010).

Figure 8 displays the coefficient of variation (CV) of $\varepsilon_V$, as a function of $\Theta^*$ and $\tau^*$, for MG and nMG fields. Results show that the region of maximum reactor efficiency is also that of minimum uncertainty (CV). Hence, the EIE system proposed here to enhance mixing-driven reactions in predefined treatment areas not only enhances mixing, but is also capable of reducing its uncertainty, which makes the system more reliable and less dependent on heterogeneity. From a practical perspective, this means that EIE can enhance remediation efficiency and reduce the risk of not reaching target remediation goals at the same time.

Accepted Article
<figure><img src="image_9.png" /><figcaption>Figure 7: Particle positions at time $t^* = 8$ with $\tau^* \in [0.1, 1, 10]$, for 2 configurations of $\Theta^*$ ($\frac{3}{5}$ and $\frac{2}{3}$), in a representative Multigaussian field (upper two) and non-Multigaussian field (lower two). The dotted line encircles the wells area and the continuous line delimits the treatment zone.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 8: Coefficient of variation (CV) of ε<sub>V</sub> as a function of Θ* and τ* for MG and nMG fields.</figcaption></figure>

Similar results were obtained by Wang et al. (2022) for multiphase flow systems, where the authors demonstrated that EIE can not only enhance the removal of non-aqueous phase liquids (NAPLs) but also reduce its uncertainty. Yet, important differences exist between NAPL removal and contaminant transport. NAPL removal efficiency should first undergo an early stage with detrimental effects in order to enhance removal in the long-term. Here, we found that the optimal stirring protocols are capable of enhancing the volume-control reactor efficiency at all times.

## 6.2 Breaking the preferential flow paths

Preferential flow paths could have a detrimental effect on in-situ treatment tech-
nologies during remediation, because the injected treatment solution might concentrate
in high-permeability regions, leaving a large portion of the system without treatment.
We wish to understand whether the enhancement of mixing produced by the EIE sys-
tem deactivates this channeling effect. For this, we examine the Lagrangian semivari-
ogram of the *Y* values sampled along particle paths. While the sill (σ²) of the semivar-
iogram quantifies the variability of *Y* along particle paths, the integral time-scale (*I*_t)
indicates how quickly these *Y* values are visited. We assume that if particles tend to visit
a large variety of *Y* values in a short time, they can easily escape from preferential flow
paths.

Figure 9 shows the ensemble average of the Lagrangian semivariogram over all re-
alizations of MG (red lines) and nMG fields (blue lines), for all rotation periods and 4
representative rotation angles Θ*. In equivalent scenarios, MG fields always exhibit higher
semivariogram values, indicating more effective spatial sampling than in nMG fields.

Accepted Article
<figure><img src="image_11.png" /><figcaption>Figure 9: Ensemble average of Lagrangian semivariograms for all rotation periods $\tau^*$ and four chosen rotation angles $\Theta^*: \Theta^* = \frac{2}{3}; \Theta^* = \frac{3}{5}; \Theta^* = \frac{1}{2}$; and $\Theta^* = 1$</figcaption></figure>

The sill ($\sigma^2$) and the integral time-scale ($I_t$) of the Lagrangian semivariogram of $Y$, together with the sampling effectiveness ($\Omega \propto \sigma^2 I_t^{-1}$), are depicted in Fig. 10 as a function of $\Theta^*$ and $\tau^*$. For comparison purposes, this Figure also overlaps the contour line of $\varepsilon_V = 0.55$. This region of optimal dilution displays considerable overlap with that of reduced integral time-scale ($I_t$), both for MG and nMG fields. Note that a low value of $I_t$ suggests high propensity of zone change by particles. Hence, this result underpins the notion, established in section 6.1, that the optimal protocols are capable of partially overcoming the detrimental effect of preferential flow paths. However, nMG fields do exhibit a noticeably lower sill ($\sigma^2$) within this optimal dilution region, which is a persistent consequence of their higher connectivity, and which explains the slightly lower values of $\varepsilon_V$ displayed by nMG fields with respect to their MG counterparts (see section 6.1 above).

The stirring protocols leading to maximum volume-control reactor efficiency and sampling effectiveness are somewhat different. Sampling effectiveness requires slightly higher rotation periods. Nevertheless, there is a region of overlap between high reactor dilution and high Lagrangian sampling effectiveness. Note that the latter does not measure whether the particles stay within the remediation zone.

Accepted Article
<figure><img src="image_12.png" /><figcaption>Figure 10: Zonation maps of the parameters that summarize the information given by the Lagrangian semivariograms: (a) the sill $\sigma^2$, (b) the integral time scale $I_t$ and (c) the sampling effectiveness $\Omega \propto \sigma^2 I_t^{-1}$, as a function of $\Theta^*$ and $\tau^*$, overlapped by the contour of Fig. 6(c), corresponding to $\varepsilon_V \ge 0.55$, considered as the optimal dilution zone.</figcaption></figure>

## 6.3 On recirculation

In this section we examine the role of treatment solution recirculation between ex-traction and injection wells during the application of EIE. Figure 11 shows the ensemble average of the frequency of recirculation over all realizations of MG and nMG fields as a function of $\Theta^*$ and $\tau^*$ after all simulation time $t^* = 20$.

Results show that, while recirculation does clearly contribute to mixing, the overlap with the high-dilution region $\varepsilon_V \ge 0.55$ (delimited by dotted lines) is only partial. The higher $\tau^*$ values within this optimized region exhibit only moderate frequency of recirculation. These are the same optimal stirring protocols that display high Lagrangian sampling efficiency (which is not affected by recirculation events, see sections 5.2 and 6.2). This suggests that each of these two mechanisms – (1) heterogeneity sampling by moving through the medium, and (2) recirculation-induced mixing – contributes, separately, to the overall picture of dilution. In cases where a high recirculation rate should preferably be avoided (such as when dealing with clogging risks), this upper region of high $\tau^*$ might be the adequate choice since it can provide high mixing efficiency with relatively low recirculation frequency. On the other hand, the singular case of $\Theta^* = 2/3$, which has the potentially advantageous feature of low sensitivity of $\varepsilon_V$ to $\tau^*$, also displays a low frequency of recirculation. Hence, this configuration could be beneficial for the aforementioned type of field applications where clogging is a potential concern, although dilution efficiency will not be as high.

Accepted Article
<figure><img src="image_13.png" /><figcaption>Figure 11: Ensemble average of the frequency of recirculation over $t^* \in [0, 20]$ over all realizations of MG and nMG fields as a function of $\Theta^*$ and $\tau^*$. For comparison purposes, the contour of Fig. 6(c) corresponding to the values of $\varepsilon_V \ge 0.55$ is overlapped.</figcaption></figure>

# 7 Summary and conclusions

This work assesses the use of Engineered Injection-Extraction (EIE) for enhancing the mixing between an injected treatment solution and the contaminants emplaced in a specific remediation treatment zone. Chaotic advection is generated through periodically activating well dipoles separated by a rotation angle $\Theta$ for a period $\tau$. To evaluate the method, flow and transport simulations are performed for multiple realizations of randomly heterogeneous log-conductivity fields specifically designed to compare the performance of the method in MultiGaussian (MG) and high-connectivity non-MultiGaussian (nMG) fields. These nMG fields display well-connected high permeability structures that concentrate the flow in preferential paths and complicate remediation efforts.

To evaluate the performance of EIE, we propose the use of a new volume-control reactor efficiency metric $\varepsilon_V$ which aims at maximizing both the mixing of the injected treatment solution and its containment within a designated remediation zone. In addition, Lagrangian semivariograms of $Y$ along particle trajectories are used to explicitly evaluate whether preferential flow paths are broken during EIE. The main findings are listed below:

1. Optimal EIE stirring protocols ($\Theta$, $\tau$) exist which maximize the volume-control reactor efficiency within the treatment zone. These optimal stirring protocols are independent of the type of random field and its connectivity structure. For small rotation angles, $\Theta < 3\pi/5$, the optimal rotation period can be approximately estimated as $\tau = 2t_c\Theta/\pi$. From a practical point of view, there are no clear benefits in using a large number of dipoles, since the gain in volume-control reactor efficiency is not relevant compared to the increment in the number of wells and complexity of the EIE system; 2 or 3 dipoles with $\Theta = \pi/2$ and $\Theta = \pi/3$ seem sufficient to reach a high reactor efficiency in practical applications. An EIE system with $\Theta = 2\pi/3$ achieves a slightly lower reactor efficiency than these optimal stirring protocols but it is less dependent on the rotation period. This can be an adequate choice in cases where the rotation period is not meant to be optimized. This singular case ($\Theta = 2\pi/3$) can also be advantageous when dealing with clogging risks, since it features a low recirculation rate of the treatment solution.

Accepted Article
2. The optimization of the stirring protocol is important for designing an adequate EIE operation. With optimal stirring protocols, the volume-control reactor efficiency increases monotonically to a maximum value after approximately 8 characteristic times, $t = 8Q/\pi R^2 b\phi$, where Q is the pumping rate, R is the distance between injection and extraction wells, b is the aquifer thickness, and $\phi$ the porosity. For $\Theta < 3\pi/5$, that is to say that EIE approaches maximum strength after 2 rotations. Suboptimal stirring protocols exhibit a weak response and/or a late-time decline of the reactor efficiency. When $\tau$ is too small ($\tau \ll t_c$), the injected treatment solution is trapped near the injection location, and when $\tau$ is too high ($\tau \gg t_c$), a large amount of the injected solution is lost outside the treatment zone.

3. Heterogeneity typically complicates remediation efforts and makes field applications highly uncertain. Optimized EIE not only enhances the mixing between the injected treatment solution and the contaminants, but also reduces its uncertainty, making the remediation outcome more reliable and less dependent on heterogeneity. This means that EIE can improve remediation effectiveness and at the same time reduce the risk of not reaching remediation targets.

4. The presence of preferential channels in nMG fields can reduce the optimal reactor efficiency of the EIE system. In our case, when $\sigma^2=2$, we obtain a factor of 1.14 compared to MG fields. Optimal EIE protocols can generate reduced temporal correlations of visited aquifer zones along particle paths, both for MG and nMG fields. The results lead us to conclude that EIE helps to reduce the impact of preferential flow channels on remediation but cannot completely eliminate the tendency of contaminants to enter into preferential channels.

5. Maximum (Eulerian) volume-control reactor efficiency and maximum (Lagrangian) sampling effectiveness involve different stirring protocols. Sampling effectiveness requires slightly larger rotation periods to force the particles to escape from preferential channels. This effect is minor and, when wisely chosen, optimal stirring protocols (i.e., with high volume-control reactor efficiency) also exhibit large sampling effectiveness.

## Acknowledgments

This work was financially supported by GRADIENT (PID2021-127911350-OB-100), the Catalan Research Project RESTORA (ACA210/18/00040), and the European research project MixUP (HORIZON-MSCA-2021-PF-01-101068306). Additional funding was obtained from the Generalitat de Catalunya (2017 SGR1485). The authors thankfully acknowledge the computer resources at TITANI and the technical support provided by the CaminsTECH. The output data used to elaborate the figures in section 6 is available on Zenodo (https://doi.org/10.5281/zenodo.7753809).

## References

Aref, H. (1984). Stirring by chaotic advection. *Journal of fluid mechanics*, **143**, 1-21.

Bagtzoglou, A. C., & Oates, P. M. (2007). Chaotic advection andenhanced groundwater remediation. *J. Mat. Civil Eng.*, **19**, 75-83.

Boggs, J. M., & Adams, E. E. (1992). Field study of dispersion in a heterogeneous aquifer. 4. investigation of adsorption and sampling bias. *Water Resources Research*, **28** (12), 3325-3336.

Cho, M. S., Solano, F., Thomson, N. R., Trefry, M. G., Lester, D. R., & Metcalfe, G. (2019). Field trials of chaotic advection to enhance reagent delivery. *Ground-

437

438

439

440

441

442

443

444

445

446

447

448

449

450

451

452

453

454

455

456

457

458

459

460

461

462

463

464

465

466

467

468

469

470

471

472

473

474

475

476

477

478

479

480

481

482

483

484

Accepted Article
water Monitoring & Remediation, 39(3), 23-39.

Cirpka, O. A. (2002). Choice of dispersion coefficients in reactive transport calculations on smoothed fields. *Journal of Contaminant Hydrology*, 58, 261282. doi: https://doi.org/10.1016/S0169-7722(02)00039-6

Cirpka, O. A., Frind, E. O., & Helmig, R. (1999). Numerical simulation of biodegradation controlled by transverse mixing. *Journal of Contaminant Hydrology*, 40, 159-182. doi: https://doi.org/10.1016/S0169-7722(99)00044-3

Cirpka, O. A., & Kitanidis, P. K. (2001). Travel-time based model of bioremediation using circulation wells. *Groundwater*, 39(3), 422-432.

Delay, F., Ackerer, P., & Danquigny, C. (2005). Simulating solute transport in porous or fractured formations using random walk particle tracking: a review. *Vadose Zone Journal*, 4(2), 360-379.

de Simoni, M., Carrera, J., Snchez-Vila, X., & Guandagnini, A. (2005). A procedure for the solution of multicomponent reactive transport problems. *Water Resour. Res.*, 41, W11410. doi: doi:10.1029/2005WR004056

Di Dato, M., de Barros, F. P., Fiori, A., & Bellin, A. (2018). Improving the efficiency of 3-d hydrogeological mixers: Dilution enhancement via coupled engineering-induced transient flows and spatial heterogeneity. *Water Resources Research*, 54(3), 2095-2111.

Gandhi, R. K., Hopkins, G. D., Goltz, M. N., Gorelick, S. M., & McCarty, P. L. (2002). Full-scale demonstration of in situ cometabolic biodegradation of trichloroethylene in groundwater 2. comprehensive analysis of field data using reactive transport modeling. *Water Resources Research*, 38(4), 11-1.

Gmez-Hernndez, J. J., & Wen, X.-H. (1998). To be or not to be multi-gaussian? a reflection on stochastic hydrogeology. *Advances in Water Resources*, 21, 47-61.

Harbaugh, A., Langevin, C., Hughes, J., Niswonger, R., & Konikow, L. (2017). Modflow-2005 version 1.12.00, the us geological survey modular groundwater model: *Us geological survey software release, 03 february 2017*. F7RF5S7G.

Hashemi, S., Javaherian, A., Ataee-pour, M., & Khoshdel, H. (2014). Two-point versus multiple-point geostatistics: the ability of geostatistical methods to capture complex geobodies and their facies associationsan application to a channelized carbonate reservoir, southwest iran. *Journal of Geophysics and Engineering*, 11. doi: https://doi.org/10.1088/1742-2132/11/6/065002

Kitanidis, P. (1994). The concept of the dilution index. *Water Resour. Res.*, 30 (7), 20112026.

Knudby, C., & Carrera, J. (2005). On the relationship between indicators of geo- statistical, flow and transport connectivity. *Advances in Water Resources*, 28, 405-421. doi: 10.1016/j.advwatres.2004.09.001

LaBolle, E. M., Quastel, J., Fogg, G. E., & Gravner, J. (2000). Diffusion processes in composite porous media and their numerical integration by random walks: Generalized stochastic differential equations with discontinuous coefficients. *Water Resources Research*, 36(3), 651-662.

Lester, D., Dentz, M., Le Borgne, T., & de Barros, F. P. J. (2018). Fluid deformation in random steady three-dimensional flow. *Journal of Fluid Mechanics*, 855, 770803. doi: 10.1017/jfm.2018.654

Lester, D., Metcalfe, G., Trefry, M., Ord, A., Hobbs, B., & Rudman, M. (2009). Lagrangian topology of a periodically reoriented potential flow: Symmetry, optimization, and mixing. *Physical Review E*, 80(3), 036208.

Lester, D., Rudman, M., Metcalfe, G., Trefry, M., Ord, A., & Hobbs, B. (2010). Scalar dispersion in a periodically reoriented potential flow: Acceleration via lagrangian chaos. *Physical Review E*, 81(4), 046319.

Luo, J., Wu, W.-M., Carley, J., Ruan, C., Gu, B., Jardine, P. M., ... Kitanidis, P. K. (2007). Hydraulic performance analysis of a multiple injection-extraction well system. *Journal of Hydrology*, 336(3-4), 294-302.

<figure><img src="image_14.png" /></figure>
540

541

542

543

544

545

546

547

548

549

550

551

552

553

554

555

556

557

558

559

560

561

562

563

564

565

566

567

568

569

570

571

572

573

574

575

576

577

578

579

580

581

582

583

584

585

586

587

588

589

590

591

592

593

594

MacDonald, T. R., Kitanidis, P. K., McCarty, P. L., & Roberts, P. V. (1999a). Effects of shear detachment on biomass growth and in situ bioremediation. *Groundwater, 37*(4), 555-563.

MacDonald, T. R., Kitanidis, P. K., McCarty, P. L., & Roberts, P. V. (1999b). Mass-transfer limitations for macroscale bioremediation modeling and implications on aquifer clogging. *Groundwater, 37*(4), 523-531.

Mays, D. C., & Neupauer, R. M. (2012). Plume spreading in groundwater by stretching and folding. *Water Resources Research, 48*(7).

McCarty, P. L., Goltz, M. N., Hopkins, G. D., Dolan, M. E., Allan, J. P., Kawakami, B. T., & Carothers, T. (1998). Full-scale evaluation of in situ cometabolic degradation of trichloroethylene in groundwater through toluene injection. *Environmental Science & Technology, 32*(1), 88-100.

McGregor, R., & Benevenuto, L. (2021). The effect of heterogeneity on the distribution and treatment of pfas in a complex geologic environment. *Frontiers in Environmental Chemistry, 2*, 729779.

Metcalfe, G., Lester, D., Trefry, M., & Ord, A. (2008). Transport in a partially open porous media flow. In *Complex systems ii* (Vol. 6802, p. 680201).

Neupauer, R. M., Meiss, J. D., & Mays, D. C. (2014). Chaotic advection and reaction during engineered injection and extraction in heterogeneous porous media. *Water Resour. Res.*, 50, 14331447. doi: doi:10.1002/2013WR014057

Piscopo, A. N., Neupauer, R. M., & Mays, D. C. (2013). Engineered injection and extraction to enhance reaction for improved in situ remediation. *Water Resour. Res.*, 49, 3618-3625. doi: doi:10.1002/wrcr.20209

Remy, N. (2005). S-gems: the stanford geostatistical modeling software: a tool for new algorithms development. In *Geostatistics banff 2004* (pp. 865-871). Springer.

Renard, P., & Allard, D. (2013). Connectivity metrics for subsurface flow and transport. *Advances in Water Resources, 51*, 168-196. doi: doi:10.1016/j.advwatres.2011.12.001

Rodrguez-Escales, P., Fernndez-Garcia, D., Drechsel, J., Folch, A., & Sanchez-Vila, X. (2017). Improving degradation of emerging organic compounds by applying chaotic advection in managed aquifer recharge in randomly heterogeneous porous media. *Water Resour. Res.*, 49, 43764392. doi: doi:10.1002/2016WR020333

Rolle, M., Clement, T. P., Sethi, R., & Molfetta, A. D. (2008). A kinetic approach for simulating redox-controlled fringe and core biodegradation processes in groundwater: model development and application to a landfill site in piedmont, italy. *Hydrol. Process.*, 22, 49054921. doi: 10.1002/hyp.7113

Salamon, P., Fernández-Garcia, D., & Gómez-Hernández, J. J. (2006). A review and numerical assessment of the random walk particle tracking method. *Journal of contaminant hydrology, 87*(3-4), 277-305.

Sather, L. J., Roth, E. J., Neupauer, R. M., Crimaldi, J. P., & Mays, D. C. (2023). Experiments and simulations on plume spreading by engineered injection and extraction in refractive index matched porous media. *Water Resources Research, 59*(2), e2022WR032943.

Schulze-Makuch, D., & Cherkauer, D. (1998). Variations in hydraulic conductivity with scale of measurement during aquifer tests in heterogeneous, porous carbonate rocks. *Hydrogeology Journal, 6*, 204-215. doi: https://doi.org/10.1007/s100400050145

Sole-Mari, G., Fernndez-Garcia, D., Sanchez-Vila, X., & Bolster, D. (2019). Particle density estimation with grid-projected and boundary-corrected adaptive kernels. *Advances in Water Resources, 131*. doi: https://doi.org/10.1016/j.advwatres.2019.103382

Sole-Mari, G., Riva, M., Fernndez-Garcia, D., Sanchez-Vila, X., & Guadagnini, A. (2021). Solute transport in bounded porous media characterized by generalized

Accepted Article
sub-gaussian log-conductivity distributions. *Advances in Water Resources*, 147, 103812.

Speetjens, M., Metcalfe, G., & Rudman, M. (2021). Lagrangian transport and chaotic advection in three-dimensional laminar flows. *Applied Mechanics Reviews*, 73(3).

Sturman, P., Stewart, P., Cunningham, A., Bouwer, E., & Wolfram, J. (1995). Engineering scale-up of in situ bioremediation processes: a review. *Journal of Contaminant Hydrology*, 19, 171-203. doi: https://doi.org/10.1016/0169-7722(95)00017-P

Tartakovsky, A. M. (2010). Langevin model for reactive transport in porous media. *Physical Review E*, 82(2), 026302.

Trefry, M. G., Lester, D. R., Metcalfe, G., Ord, A., & Regenauer-Lieb, K. (2012). Toward enhanced subsurface intervention methods using chaotic advection. *Journal of Contaminant Hydrology journal*, 127, 15-29. doi: 10.1016/j.jconhyd.2011.04.006

Turuban, R., Lester, D. R., Heyman, J., Borgne, T. L., & Mheust, Y. (2019). Chaotic mixing in crystalline granular media. *Journal of Fluid Mechanics*, 871, 562594. doi: 10.1017/jfm.2019.245

Wang, Y., Fernández-Garcia, D., Sole-Mari, G., & Rodríguez-Escales, P. (2022). Enhanced napl removal and mixing with engineered injection and extraction. *Water Resources Research*, e2021WR031114.

Zhang, P., Devries, S. L., Dathe, A., & Bagtzoglou, A. C. (2009). Enhanced mixing and plume containment in porous media under time-dependent oscillatory flow. *Environ. Sci. Technol.*, 43, 62836288. doi: 10.1021/es900854r

Zinn, B., & Harvey, C. F. (2003). When good statistical models of aquifer heterogeneity go bad: A comparison of flow, dispersion, and mass transfer in connected and multivariate gaussian hydraulic conductivity fields. *Water Resources Research*, 39(3).

Accepted Article

595

596

597

598

599

600

601

602

603

604

605

606

607

608

609

610

611

612

613

614

615

616

617

618

619

620

621

622