# Subsurface Mixing Dynamics Across the Salt-Freshwater Interface

K. de Vriendt, Tanguy Le Borgne, M. Pool, Marco Dentz

► To cite this version:

K. de Vriendt, Tanguy Le Borgne, M. Pool, Marco Dentz. Subsurface Mixing Dynamics Across the Salt-Freshwater Interface. Geophysical Research Letters, 2022, 49 (7), pp.e2021GL097094. 10.1029/2021gl097094. insu-03707425

HAL Id: insu-03707425

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

Submitted on 28 Jun 2022

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
# Geophysical Research Letters®

## RESEARCH LETTER

10.1029/2021GL097094

**Key Points:**

* The width of the saltwater-freshwater interface and mixing rate varies along the interface

* Interface width and mixing dynamics are related through flow deformation

* We derive an analytical model that predicts the mixing width and rate based on the interface profile and transverse dispersivity

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

M. Dentz,
marco.dentz@csic.es

**Citation:**

De Vriendt, K., Le Borgne, T., Pool, M., & Dentz, M. (2022). Subsurface mixing dynamics across the salt-freshwater interface. *Geophysical Research Letters*, *49*, e2021GL097094. https://doi.org/10.1029/2021GL097094

Received 17 NOV 2021

Accepted 14 MAR 2022

© 2022. The Authors.

This is an open access article under
the terms of the Creative Commons

Attribution-NonCommercial-NoDerivs

License, which permits use and
distribution in any medium, provided the
original work is properly cited, the use is
non-commercial and no modifications or
adaptations are made.

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

## Subsurface Mixing Dynamics Across the Salt-Freshwater Interface

K. De Vriendt¹ ID, T. Le Borgne² ID, M. Pool³ ID, and M. Dentz¹ ID

¹Institute of Environmental Assessment and Water Research, IDAEA-CSIC, Barcelona, Spain, ²Université de Rennes 1, CNRS, Géosciences Rennes UMR 6118, Rennes, France, ³AMPHOS 21 Consulting S. L., Barcelona, Spain

**Abstract** Mixing along the salt-freshwater interface is critical for geochemical reactions, transport, and transformation of nutrients and contaminants in coastal ecosystems. However, the mechanisms and controls of mixing are not well understood. We develop an analytical model, based on the coupling between flow deformation and dispersion, which predicts the mixing dynamics along the interface for steady-state flow in coastal aquifers. The analytical predictions are compared with the results of detailed numerical simulations, which show that nonuniform flow fields, inherent to seawater intrusion in coastal aquifer, result in a non-monotonic evolution of mixing width and mixing rates along the interface. The analytical model accurately captures these dynamics over a range of freshwater flow rates and dispersivities. It predicts the evolution of the mixing width and mixing rates along the interface, offering a new framework for understanding and modeling mixing and reaction processes in coastal aquifers.

**Plain Language Summary** Density differences between salt and freshwater lead to the formation of a convection cell in coastal aquifers, in which seawater intrudes inland along the aquifer bottom. Fresh and mixed waters flow upward along the salt-freshwater interface and are forced to accelerate before being discharged along the ocean seabed. The resulting nonuniform flow alters the concentration of the mixed waters along the interface, which in turn enhances mixing rates and creates local mixing hotspots. Our results show how nonuniform velocity fields result in enhanced local mixing dynamics and elucidate the mechanisms and controls of mixing processes along salt-freshwater interfaces in coastal aquifers.

### 1. Introduction

Coastal aquifers are some of the most vulnerable groundwater resources sustaining dense coastal populations globally (Ferguson & Gleeson, 2012). These subsurface environments are subject to significant anthropogenic pollutants that negatively impact ocean ecosystems (Kroeger & Charette, 2008; Moore, 2010; Slomp & Cappellen, 2004). Moreover, their inherently nonstationary flow dynamics on different temporal scales (tides, seasons, and glacial cycles) leads to a range of geochemical processes across coastal landscapes. A notable example is mixing-enhanced carbonate dissolution and karstification processes in coastal zones (Back et al., 1986). Over large time scales, seawater intrusion has acted as a primary mechanism to observable land features, such as the formation of “Flank Margin Caves” near the mixing discharge zone (Back et al., 1979; Mylroie & Carew, 1990) or cave and conduits formation in Bermudas (A. Palmer, 1992), Bahamas (R. Palmer & Williams, 1982), and Yucatán (Back et al., 1986). Freshwater discharge in coastal aquifers has also been associated with a variety of other biogeochemical reactions in beach environments. A well-known example is the enhanced iron oxide precipitation in Waquiot Bay (termed “iron curtain”) (Charette & Sholkovitz, 2002; Spiteri, Slomp, Charette, et al., 2008; Spiteri, Slomp, Tuncay, et al., 2008), which attenuates contaminants, such as phosphates and arsenic. Such reactions may hold a strong propensity in regulating the flux of terrestrial pollutants toward coastal marine ecosystems.

While reaction kinetics and redox conditions are strong precursors to these reactive hotspots, their interplay with the nonuniform velocity field and mixing dynamics in coastal aquifers remains poorly understood. Sanford and Konikow (1989) and Rezaei et al. (2005) demonstrated numerically that the mixing of salt and freshwater in coastal aquifers induces local dissolution hotspots at both the discharge zone as well as at the toe of the saltwater wedge. Studies have since also highlighted the importance of heterogeneity across the salt-freshwater interface (SFI) in generating local reaction hotspots (De Vriendt et al., 2020).
A key challenge for capturing mixing and reaction hotspots is to quantify the size of the mixing zone between freshwater and saltwater, which sets concentration gradients and thus mixing rates across the interface. Under steady-state and homogeneous conditions, mixing across the SFI is dominantly controlled by density effects and transverse dispersion (Abarca et al., 2007; Paster & Dagan, 2007). Laboratory-scale experiments (e.g., Abarca et al., 2007; Goswami & Clement, 2007; Robinson et al., 2015; Yoshihiro et al., 2010) and some field observations (Paster, 2010) have shown relatively sharp mixing zones with small widths compared to the aquifer scale. On the other hand, large-scale field studies have observed mixing zones ranging from tens to hundreds of meters (Barlow, 2003; Kim et al., 2007; Kroeger & Charette, 2008; Langevin, 2003; Price et al., 2003; Spiteri, Slomp, Charette, et al., 2008; Spiteri, Slomp, Tuncay, et al., 2008). Widening of the mixing zones in real-world coastal aquifers has mainly been attributed to transient effects, such as tides (e.g., Ataie-Ashtiani et al., 1999; Pool et al., 2014, 2015) as well as heterogeneity (Abarca Cameo, 2006; Kerrou & Renard, 2010; Lu et al., 2013) or kinetic mass transfer (Lu et al., 2009). However, while all these investigations provide valuable insight into water resources management and general mixing dynamics, in these studies, the width of the mixing zone has been addressed mainly through averaging across and along the saltwater-freshwater interface (e.g., Abarca et al., 2007; Kerrou & Renard, 2010; Lu et al., 2013; Pool et al., 2014). Therefore, how the mixing widths vary along the interface and what are the mechanisms driving the formation of mixing and reaction hotspots are outstanding questions. Recent theoretical developments have demonstrated that fluid stretching in nonuniform flow fields can lead to increased local mixing and reactions (e.g., Bandopadhyay et al., 2018; Le Borgne et al., 2014). Here, we apply these concepts to investigate the impact of flow deformation, driven by velocity gradients inherent to salt-freshwater interfaces on mixing dynamics across the SFI. We quantify the evolution of the mixing width along the SFI for a range of freshwater flow rates and dispersivities and relate these dynamics to the stretching rate driven by nonhomogeneous flow along the interface. We derive an analytical solution, which provides accurate predictions of the mixing dynamics along the SFI and allows understanding and modeling the development of mixing hotspots. We discuss the implications of our findings regarding their impact on mixing and reaction rates in coastal aquifers.

## 2. Methods

### 2.1. Flow and Transport

We study mixing under steady variable density flow in a two-dimensional cross section of a coastal aquifer. Density-dependent flow is described by the Darcy equation

$$ \mathbf{q} = -K \left( \nabla h_f + \frac{\rho - \rho_f}{\rho_f} \mathbf{e}_z \right), \qquad (1) $$

where $\mathbf{q}$ is the specific discharge, $K$ is the hydraulic conductivity, $h_f$ the equivalent freshwater head, $\rho$ the fluid density, $\rho_f$ the density of freshwater, and $\mathbf{e}_z$ is the unit vector in y direction. Fluid mass conservation in the absence of sources and sinks implies $\nabla \cdot \rho \mathbf{q} = 0$. The fluid density is assumed to be linearly dependent on the salt mass fraction $\omega$ (mass of salt dissolved per unit mass of fluid) given by $\rho = \rho_f[1 + \epsilon'c]$, where $\epsilon'$ is the buoyancy factor given by $\epsilon' = (\rho_s - \rho_f)/\rho_f$ with $\rho_s$ the density of seawater and $c$ is the normalized salt concentration defined as $c = \omega/\omega_s$ with $\omega_s$ the salt mass fraction of seawater. The concentration $c$ evolves according to the advection dispersion equation, which in the steady state reads as

$$ \mathbf{q} \cdot \nabla c - \nabla \cdot (\mathbf{D} + \phi D_m \mathbf{1}) \nabla c = 0, \qquad (2) $$

where $\mathbf{D}$ is the dispersion tensor (Bear, 1988), $D_m$ is the molecular diffusion coefficient, and $\phi$ is the porosity. We consider here a uniform hydraulic conductivity and assume that the impact of subscale heterogeneity is captured by the dispersivity. For this particular problem, the key dimensionless numbers that emerge are two Péclet numbers, $Pe_I$, which compares the advection and dispersion times, and $Pe_{II}$, which compares the advection and diffusion times, and the gravity number, $Ng$, which compares the viscous $q/K$ and buoyancy forces $\epsilon'$ (see Supplementary Information) (see Abarca et al., 2007),

$$ Pe_I = \frac{b}{\alpha_t}, \quad Pe_{II} = \frac{q_f b}{\phi D_m}. \quad Ng = \frac{K \epsilon'}{q_f}, \qquad (3) $$
<figure><img src="image_3.png" /><figcaption>**Figure 1.** (a) Steady-state concentration map for $\text{Ng} = 17.3$. The figure illustrates the prescribed freshwater flux boundary on the left and hydrostatic head conditions on the right boundary. The inset image depicts a map of $c(1 - c)$ along with a local profile of $c(1 - c)$ perpendicular to the interface along the *n*-coordinate. (b) Map of the concentration gradient and (c) the scalar dissipation rate.</figcaption></figure>

where *b* defines the domain thickness, $\alpha_t$ is the transverse dispersivity, $q_f$ is the specified fresh water flux, and $\phi$ is the porosity.

## 2.2. Numerical Model

We consider a shallow coastal aquifer of constant thickness *b* and length *L* extended offshore with a specific freshwater discharge from inland $q_f$ (see Figure 1a). The connection with the sea is represented as a prescribed head along the offshore model top and the offshore vertical boundaries. Different values for the fresh water flux and for the longitudinal and transverse dispersivities have been considered to evaluate their impact on mixing along the interface. The base-case scenario used in this study is largely inspired from the study of Spiteri, Slomp, Charette, et al. (2008) and Spiteri, Slomp, Tuncay, et al. (2008). However, the general relationship between fluid stretching and mixing dynamics derived from this numerical example is expected to apply more generally over a large range of coastal aquifer systems.

The values used for longitudinal and transverse dispersivities are based on typical literature values where numerical simulations were calibrated to field measurements (see Table S2 in Supporting Information S1). The values chosen for $Pe_I$ and $Pe_{II}$ are consistently larger than unity as typically found in field studies and laboratory experiments (see Table S2 in Supporting Information S1). A summary of the parameters used in the numerical simulations is provided in Table S1 in Supporting Information S1. The freshwater flux ranges from $q_f = 1.25 \times 10^{-2}$ m/d to $3 \times 10^{-2}$ m/d. Thus, the simulated scenarios are characterized by a $Pe_I$ of 500 and $\text{Ng}$ ranging between 17.3 and 7.2. Since we vary only the flow rate, the range of $\text{Ng}$ considered is equivalent to the one of $Pe_{II}$. Therefore, in the following, the scenarios are characterized by their $\text{Ng}$ values. It should be noted that the gravity number in general plays a fundamental role in the movement of the wedge and has also been shown to play an important role in the mixing of stable stratification problems (Dell'Oca et al., 2018).
## 2.3. Mixing Measures

The variability of mixing along the SFI can be characterized by the local scalar dissipation rate, which is defined using

$$ \chi = \nabla c \cdot (\mathbf{D}\nabla c). \qquad (4) $$

For reversible mixing-limited reactions, this measure is directly proportional to the reaction rate (De Simoni, 2005). In order to separate the impact of (velocity-dependent) dispersion and concentration gradient in the scalar dissipation rate, we also consider the concentration gradient,

$$ \theta = \|\nabla c\|, \qquad (5) $$

where $\|\cdot\|$ denotes the L²-norm. The salt concentration gradient at the SFI can be approximated by $\theta \approx c/s$, where $c_s$ is the concentration of salt in the seawater and $s$ is the interface width. Accordingly, the evolution of the concentration gradient and thus the mixing rate is determined by the interface width. The interface width is therefore a crucial element toward understanding the mixing dynamics (Paster & Dagan, 2007; Abarca et al., 2007). The width of the mixing zone normal to the principal direction of flow is determined from the width of the auxiliary function $c(1 - c)$ as detailed in Section 1.2 of the supporting information. All quantities are evaluated along the curvilinear length of the interface, where the toe is located at $z = 0$. We compare the scalar dissipation rate and the gradient of concentration by evaluating their local maximum values at a given depth along the length of the interface. Finally, we evaluate the rate of strain to highlight zones of enhanced fluid strain, $\Theta_c$, across the interface, where flow deformation may compress the mixing zone and thus enhance concentration gradients (De Barros et al., 2012).

These concepts are illustrated in Figure 1, which shows the general mixing and flow features for a saltwater wedge at a steady state. The model setup and a visual description of the mixing width are also provided in Figure 1a.

# 3. Mixing Mechanisms and Mixing Model

Figure 1b shows the evolution of the concentration gradient, which is maximum at the toe and head. This evolution is also reflected in the mixing rate (see Figure 1c). This behavior indicates that the width, which is inversely proportional to the concentration gradient, is small at the toe and head and evolves non-monotonically in between. To illustrate the relation of the flow deformation, we also include a map of the rate of strain (De Barros et al., 2012; Okubo, 1970; Weiss, 1991) (Figure 1d). These dynamics are quantified in the following by deriving an analytical model for the evolution of the mixing width in response to dispersion and flow deformation.

## 3.1. Mixing Along the Interface

To investigate the impact of flow deformation on the interface width, concentration gradient, and mixing rate, we vary the gravity number $Ng$ by changing the freshwater flow rate. The local mixing widths along the interface for different $Ng$ are shown in Figures 2a and 2b. The SFI is initially narrowest at the toe where the two fluids initially mix. From here $s$ broadens to a maximum value, $s_m$ before narrowing again toward the discharge zone. While it has been speculated that under velocity-dependent dispersion, the mixing width should increase with increasing freshwater flux (Werner et al., 2012), we find that the overall interface width increases for decreasing freshwater flow, that is, increasing $Ng$ (see Figure 2a). We also show that all curves can be collapsed by scaling $s$ by $s_m$ and $z$ by the toe length, $L_t$ (Figure 2b). We find that $L_t$ grows proportional to the freshwater flux, $L_t \propto Ng$, while $s_m$ decreases as $s_m \propto Ng^{1/2}$ (see Supporting Information S1). Figure 2c shows the evolution of the concentration gradient $\theta$ along the interface for different $Ng$. All $\theta$ collapse on a single curve when rescaled with their respective minima $\theta_m$ and plotted against $z/L_t$. This behavior mirrors the evolution of the mixing width as it decays from the toe toward a minimum and again increases toward the discharge. In fact, the evolution of the concentration gradient $\theta/\theta_m$ can be well represented by the inverse interface width $(s/s_m)^{-1}$. We observe the same behavior for the mixing rates (see Figure 2d), which are rescaled by their minima $\chi_m$. Their evolution is well represented by $\chi \approx \alpha_v\theta^2$ normalized by its minimum. This highlights the central role of the interface width on mixing along the interface.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** (a) Mixing width along the interface for Ng = 17.3 (purple triangles), 14.4 (pink circles), 10.8 (beige squares), 8.6 (light red triangles), and 7.2 (red circles). (b) Mixing widths scaled by the respective maximum interface widths $s_m$ versus distance along the interface scaled by the toe length $L_t$. (c) Concentration gradients scaled by the respective minimum gradients $\chi_m$. The blue-dotted line denotes the inverse mixing width $\theta/\theta_m \approx (s/s_m)^{-1}$. (d) Scalar dissipation rates scaled by their respective minima $\chi_m$. The dashed blue line denotes $\alpha_v\theta^2$ normalized by its minimum.</figcaption></figure>

## 3.2. Interface Mixing Model

The evolution of the interface width can be understood from the interplay between transverse dispersion and flow deformation. Initially, near the toe, we observe enhanced mixing reflected by high concentration gradients and mixing rates. They are attributed to a local stagnation point resulting from opposing flow, which leads to enhanced interface compression. Moving away from the toe, flow velocities accelerate, which imply stretching along the interface and at the same time, interface compression is perpendicular to the stretching direction. Near the toe, the compression rates are so low that transverse dispersion dominates over compression, and the interface width grows diffusively with distance as $z^{1/2}$ (Figures 2a and 2b). Further up the interface, freshwater velocities increase faster due to a decrease in area between the confining unit and the interface. Eventually, at a characteristic depth $z_c$, the acceleration along the interface and the concurrent compression are large enough to overcome transverse dispersion. Thus, a maximum interface width is reached, followed by a succession of compression events of increasing rates that lead to a decrease in the mixing width. A similar behavior was observed by Eeman et al. (2011) when investigating upwelling of saline water across a freshwater lens into a ditch. The authors found that despite increasing velocities toward the outlet, the mixing width continued to narrow due to converging streamlines.

The competition between hydrodynamic compression and dispersive expansion can be understood more quantitatively by the following evolution equation for the mixing width $s$ (Villermaux, 2012),

$$ \frac{1}{s} \frac{ds}{dt} = -\gamma + \frac{D_t}{s^2}, \qquad (6) $$

where $\gamma$ is the stretching (or compression) rate and $D_t/s^2$ is the dispersive expansion rate with $D_t = D_m + \alpha_v\gamma$ the transverse dispersion coefficient. The mixing time $t_s$, that is the time at which dispersion and compression equilibrate, is defined by $t_s = \ln(1 + Pe_s)/2\gamma$ where $Pe_s = s_0^2\gamma/D_t$ (Villermaux, 2019). Although in our system the compression rate varies along the interface, it is useful to consider the solution to Equation 6 for a constant $\gamma$,

$$ s = \sqrt{\frac{D_t}{\gamma} [1 - \exp(-2\gamma t)] + s_0^2 \exp(-2\gamma t)}. \qquad (7) $$
For times larger than $t_s$, the mixing width given by Equation 7 is expected to converge to the Batchelor scale $s_B = \sqrt{D_t/\gamma}$. We define the mixing distance $z_m = v_a t_s$ with $v_a$ the average velocity along the interface as the distance over which the mixing width converges to the local Batchelor scale $s_B = \sqrt{D_t/\gamma}$. The evolution of the interface width along the interface as a function of $z$ is obtained from $s(t)$ by setting $t = z/v_a$ such that $s'(z) = s(z/v_a)$. For simplicity of notation, we drop the prime in the following. Close to the toe, $z < z_c$, the compression rate is small, which implies a large mixing distance $z_m$. For $z \ll z_m$, expression (Equation 7) implies

$$s(z) = \sqrt{\alpha_t z}, \quad \text{for } z < z_c, \qquad (8)$$

where we set the transverse dispersion coefficient $D_t \approx \alpha_t v_a$. This explains the increase in the mixing width observed in Figures 2a and 2b. The dependence of $s$ on $\alpha_t$ is confirmed by additional numerical simulations for variable $\alpha_t$ (see Supporting Information S1). For increasing distance along the interface, the acceleration and thus $v$ and $\gamma$ increase notably along the interface. Assuming that $v$ and $\gamma$ change on length scales larger than the corresponding mixing distance $z_m$, then $s$ evolves in a quasi-steady manner as a succession of Batchelor scales such that

$$s(z) = \sqrt{\frac{\alpha_t v(z)}{\gamma(z)}}, \quad \text{for } z > z_c, \qquad (9)$$

where $v(z)$ and $\gamma(z)$ are the local velocity and compression rates along the interface, respectively. This second, quasi-steady regime describes the recompression of the interface after it has reached its maximum width $s_m$. We notice that $\gamma$ is given by the derivative of the flow velocity $v(z)$ along the interface, $\gamma(z) = dv(z)/dz$. Thus, we obtain for the interface width in terms of $v(z)$ the expression

$$s(z) = \sqrt{\alpha_t \left[ \frac{d \ln v(z)}{dz} \right]^{-1}}. \qquad (10)$$

This means that the interface width can be estimated from the velocity profile. In summary, the transition between dispersive growth and compression corresponds to the crossover between two competing mechanisms. Dispersive growth is overcome by accelerating flow toward the discharge zone, which stretches the interface. This leads to a compression of the mixing width in a quasi-steady manner as expressed by Equation 9.

To derive an approximate analytical solution for the mixing width during recompression toward the discharge zone, we must find an expression for $\gamma$. The velocity along the interface can be approximated by $v(z) = q/b/\xi(z)$, where $\xi(z)$ is the interface height. Inserting these approximations in Equation 9, we obtain for the evolution of the interface width in the compression regime the expression

$$s(z) = \sqrt{-\alpha_t \left[ \frac{d \ln \xi(z)}{dz} \right]^{-1}}, \qquad (11)$$

(see Supplementary Information). This means that the interface width can be estimated directly from the interface profile. In order to test this expression, we approximate the interface height by the solution of Glover (1959) as $\xi(z) = \sqrt{b^2 - 2bz/\text{Ng}'}$ (see Supplementary Information). Note that $\text{Ng}' = \text{Ng}/[1 - (\alpha_t/b)^{1/4}]$ is a modified gravity number to correct for the impact of dispersion in the interface position in the Glover solution (Lu & Werner, 2013; Pool & Carrera, 2011). Inserting the expression for $\xi(z)$ into (11), we obtain the compact expression

$$s(z) = \sqrt{\alpha_t \text{Ng}' b \left( 1 - \frac{2z}{\text{Ng}' b} \right)}. \qquad (12)$$

The analytical solution explains the scaling behavior of $s$ observed in Figure 2b. Note that the Glover solution predicts the toe length $L_t = \text{Ng}'b/2$. In fact, we can write Equation 12 as

$$s(z) = s_m \sqrt{3 \left( 1 - \frac{z}{L_t} \right)}. \qquad (13)$$
<figure><img src="image_5.png" /><figcaption>**Figure 3.** (a) Saltwater interface defined by the 50% concentration isoline. Symbols denote numerical simulation results and solid lines denote the Glover solutions. (b) Numerically determined stretching rate and stretching rates determined by Glover solution (solid lines). (c) Numerical mixing width compared against the numerically derived Batchelor scale (solid lines). The solid black line denotes dispersive growth $s \sim \sqrt{z}$ prior to interface recompression. The asterisks denote the predicted crossover width and position.</figcaption></figure>

The crossover position $z_c$ between the expansion and compression regimes is obtained by matching the solution Equation 8 for the expansion regime and Equation 12 for compression. Thus, we obtain for crossover position $z_c$ and the maximum interface with $s_m = s(z_c)$, the explicit expressions

$$z_c = Ng'b/3, \quad s_m = \sqrt{\alpha_t z_c}. \qquad (14)$$

This means that the maximum interface width and its position can be estimated from the modified gravity number and the aquifer thickness. Note that inserting $z_c$ in the Glover solution for the interface height leads to $\xi(z_c) = b/\sqrt{3}$, which gives the depth above which mixing is most active due to recompression along the interface. It is interesting to note that this depth is simply a fraction of the aquifer thickness and is independent on other system properties.

Figure 3a confirms the match of the Glover solution with the interface height determined from the direct numerical simulations for different $Pe_{II}$. We also show the predicted stretching rate along $z$ together with the data from the direct numerical simulation (Figure 3b). Note that no fitting parameter is used. Discrepancies at the toe can be attributed to local deceleration due to the stagnation zone. In addition, since the Glover solution assumes flow is forced through an infinitely small outlet rather than a gap as in the numerical simulations, $\gamma$ is overestimated as it asymptotes near the outlet. Figure 3c shows the match between the analytical expressions for the Batchelor scale and numerically derived mixing widths. Note that we multiply $\alpha_t$ by a factor of 3/4 to match the evolution of the data at a short distance from the toe. This can be traced back to the fact that the concentration profile across the interface is not Gaussian (see inset in Figure 1a). We find that the transition between dispersive growth and recompression of the interface is slightly overestimated for interfaces with small freshwater fluxes. However, in general, there is good agreement between the numerical and analytical solutions. It should be emphasized that the Glover solution used in this study is a means to approximate the position and velocity along the interface for this given problem. Naturally, for problems with different boundary conditions, the interface position and velocity field may deviate from the idealized scenario studied here and therefore require further evaluation.

## 4. Conclusion

Our study has examined mixing dynamics for seawater intrusion under steady-state conditions. Evaluation of the mixing width along the salt-freshwater interface has highlighted several mixing processes that are influenced by nonuniform flow from the mixing of saline and freshwater bodies. We find that the mixing width initially grows due to transverse dispersion up to a characteristic location where it then recompresses due to accelerating flow toward the discharge zone. Interface compression near the outlet is accompanied by enhanced concentration gradients and mixing rates. We attribute stronger mixing rates near the interface toe to enhanced local compression resulting from opposing flow, which results in a stagnation point. The expansion and recompression of the interface can be understood in terms of the flow deformation along the interface and are quantified by a mixing model that accounts for the competition of dispersive expansion and hydrodynamic compression of the interface. We show that the mixing width can be estimated from the interface profile and transverse dispersivity. Using the Glover solution for a sharp interface, we propose an analytical model that describes the initial growth of the interface width near the toe and its subsequent recompression near the outlet.

The quantification of the mechanism of mixing across the SFI resulting from variable density-induced nonuniform flow may shed light on mixing-limited reactions in coastal landscapes. This is particularly relevant when evaluating the chemical composition of submarine groundwater discharge (SGD), which is often altered
by biogeochemical reactions resulting from the mixing of salt and freshwater (Moore, 1999). Given that high concentrations of nutrients in coastal groundwater have been associated with eutrophication and the onset of algal blooms (LaRoche et al., 1997; Valiela et al., 1990), understanding mixing dynamics that lead to the transformation of chemicals along the interface warrants careful consideration. Furthermore, the understanding of the mixing dynamics at the discharge zone is of particular interest as it has been linked to an array of geochemical activity (e.g., Charette & Sholkovitz, 2002; De Vriendt et al., 2020; Kroeger & Charette, 2008; Mylroie & Carew, 1990; Rezaei et al., 2005), such as the precipitation of iron oxide and the dissolution of calcite.

Our study links the mixing evolution along the interface and the resulting mixing patterns to the flow patterns via the mechanisms of interface compression and local-scale dispersion. Due to their fundamental nature, we expect these mechanisms to govern mixing also under more realistic aquifer conditions characterized by variable thickness and lithology. The mixing patterns are expected to change according to the heterogeneity and geometry-induced flow patterns. The analytical approach provides a basis for the estimation of mixing along the interface also under more complex aquifer conditions and a compact expression for the mixing width if the interface profile is known.

## Data Availability Statement

Open Research Data can be accessed at the open repository https://digital.csic.es under the permanent identifier http://hdl.handle.net/10261/264456.

## Acknowledgments

This project has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie Grant Agreement No. 722028 (ENIGMA ITN). Maria Pool acknowledges the support of MCIN/AEI/10.13039/501100011033 through the Torres-Quevedo program (PTQ2018-010081). MD acknowledges the support of MCIN/AEI/10.13039/501100011033 through the project HydroPore (PID2019-106887GB-C31) and the Grant CEX2018-000794-S.

## References

* Abarca, E., Carrera, J., Sánchez-Vila, X., & Dentz, M. (2007). Anisotropic dispersive Henry problem. *Advances in Water Resources*, 30(4), 913-926. https://doi.org/10.1016/j.advwatres.2006.08.005

* Abarca Cameo, E. (2006). Seawater intrusion in complex geological environments (Doctoral dissertation). Technical University of Catalonia, UPC. Retrieved from http://hdl.handle.net/10261/33570

* Ataie-Ashtiani, B., Volker, R. E., & Lockington, D. A. (1999). Tidal effects on sea water intrusion in unconfined aquifers. *Journal of Hydrology*, 216(1-2), 17-31. https://doi.org/10.1016/S0022-1694(98)00275-3

* Back, W., Hanshaw, B. B., Pyle, T. E., Plummer, L. N., & Weidie, A. E. (1979). Geochemical significance of groundwater discharge and carbonate solution to the formation of Caleta Xel Ha, Quintana Roo, Mexico. *Water Resources Research*, 15(6), 1521-1535. https://doi.org/10.1029/WR015i006p01521

* Back, W., Hanshaw, B. B., Herman, J. S., & Van Driel, J. N. (1986). Differential dissolution of a Pleistocene reef in the ground-water mixing zone of coastal Yucatan, Mexico. *Geology*, 14(February), 137-140.

* Bandopadhyay, A., Davy, P., & Borgne, T. L. (2018). Shear flows accelerate mixing dynamics in Hyporheic zones and Hillslopes. *Geophysical Research Letters*, 659-668. https://doi.org/10.1029/2018GL079914

* Barlow, P. M. (2003). Ground water in freshwater-saltwater environments of the Atlantic Coast (Tech. Rep.). Geological Survey (USGS).

* Bear, J. (1988). *Dynamics of fluids in porous media*. Dover publications, inc.

* Charette, M. A., & Sholkovitz, E. R. (2002). Oxidative precipitation of groundwater-derived ferrous iron in the subterranean estuary of a coastal bay. *Geophysical Research Letters*, 29(10), 85-1-85-4. https://doi.org/10.1029/2001gl014512

* De Barros, F. P. J., Dentz, M., Koch, J., & Nowak, W. (2012). Flow topology and scalar mixing in spatially heterogeneous flow fields. *Geophysical Research Letters*, 39(8), 1-5. https://doi.org/10.1029/2012GL051302

* Dell'Oca, A., Riva, M., Carrera, J., & Guadagnini, A. (2018). Solute dispersion for stable density-driven flow in randomly heterogeneous porous media. *Advances in Water Resources*, 111(April 2017), 329-345. https://doi.org/10.1016/j.advwatres.2017.10.040

* De Simoni, M. (2005). A procedure for the solution of multicomponent reactive transport problems. *Water Resource Research*, 41, 1-16. https://doi.org/10.1029/2005WR004056

* De Vriendt, K., Pool, M., & Dentz, M. (2020). Heterogeneity-induced mixing and reaction Hot Spots facilitate Karst propagation in coastal aquifers. *Geophysical Research Letters*, 47(10), 1-12. https://doi.org/10.1029/2020GL087529

* Eeman, S., Leijnse, A., Raats, P. A., & van der Zee, S. E. (2011). Analysis of the thickness of a fresh water lens and of the transition zone between this lens and upwelling saline water. *Advances in Water Resources*, 34(2), 291-302. https://doi.org/10.1016/j.advwatres.2010.12.001

* Ferguson, G., & Gleeson, T. (2012). Vulnerability of coastal aquifers to groundwater use and climate change. *Nature Climate Change*, 2(5), 342-345. https://doi.org/10.1038/nclimate1413

* Glover, R. E. (1959). The pattern of fresh-water flow in a coastal aquifer. *Journal of Geophysical Research*, 64(4), 457-459. https://doi.org/10.1029/jz064i004p00457

* Goswami, R. R., & Clement, T. P. (2007). Laboratory-scale investigation of saltwater intrusion dynamics. *Water Resources Research*, 43(4), 1-11. https://doi.org/10.1029/2006WR005151

* Kerrou, J., & Renard, P. (2010). A numerical analysis of dimensionality and heterogeneity effects on advective dispersive seawater intrusion processes. *Hydrogeology Journal*, 18(1), 55-72. https://doi.org/10.1007/s10040-009-0533-0

* Kim, K.-Y., Chon, C.-M., & Park, K.-H. (2007). A Simple method for locating the fresh water-salt water interface using pressure data. *Ground Water*, 45(6), 723-728. https://doi.org/10.1111/j.1745-6584.2007.00349.x

* Kroeger, K. D., & Charette, M. A. (2008). Nitrogen biogeochemistry of submarine groundwater discharge. *Limnology and Oceanography*, 53(3), 1025-1039. https://doi.org/10.4319/lo.2008.53.3.1025

* Langevin, C. D. (2003). Simulation of submarine ground water discharge to a marine estuary: Biscayne bay, Florida. *Ground Water*, 41(6), 758-771.
LaRoche, J., Nuzzi, R., Waters, R., Wyman, K., Falkowski, P., & Wallace, D. (1997). Brown Tide blooms in Long Island's coastal waters linked to interannual variability in groundwater flow. *Global Change Biology*, *3*(5), 397-410.

Le Borgne, T., Ginn, T. R., & Dentz, M. (2014). Impact of fluid deformation on mixing-induced chemical reactions in heterogeneous flows. *Geophysical Research Letters*, *41*(22), 7898-7906. https://doi.org/10.1002/2014GL062038

Lu, C., Chem, Y., Zhang, C., & Luo, J. (2013). Steady-state freshwater-seawater mixing zone in stratified coastal aquifers. *Journal of Hydrology*, *9*, https://doi.org/10.1016/j.jhydrol.2013.09.017

Lu, C., Kitanidis, P. K., & Luo, J. (2009). Effects of kinetic mass transfer and transient flow conditions on widening mixing zones in coastal aquifers. *Water Resources Research*, *45*(December 2008), 1-17. https://doi.org/10.1029/2008WR007643

Lu, C., & Werner, A. D. (2013). Advances in water resources timescales of seawater intrusion and retreat. *Advances in Water Resources*, *59*, 39-51. https://doi.org/10.1016/j.advwatres.2013.05.005

Moore, W. S. (1999). The subterranean estuary: A reaction zone of ground water and sea water. *Marine Chemistry*, *65*(1-2), 111-125. https://doi.org/10.1016/S0304-4203(99)00014-6

Moore, W. S. (2010). A reevaluation of submarine groundwater discharge along the southeastern coast of North America. *Global Biological Cycle*, *24*(December 2009), 1-9. https://doi.org/10.1029/2009GB003747

Mylroie, J. E., & Carew, J. L. (1990). The flank margin model for dissolution cave development in carbonate platforms. *Earth Surface Processes and Landforms*, *15*(5), 413-424. https://doi.org/10.1002/esp.3290150505

Okubo, A. (1970). Horizontal dispersion of floatable particles in the vicinity of velocity singularities such as convergences. *Deep-Sea Research and Oceanographic Abstracts*, *17*(3), 445-454. https://doi.org/10.1016/0011-7471(70)90059-8

Palmer, A. (1992). Origin and morphology of limestone caves. *Geological Society of America Bulletin*, *103*(January 1991), 1-21.

Palmer, R., & Williams, D. (1982). Cave development under Andros island, Bahamas. *Cave Science*, *11*(1), 50-52.

Paster, A. (2010). Mixing between fresh and salt waters at aquifer regional scale and identification of transverse dispersivity. *Journal of Hydrology*, *380*(1-2), 36-44. https://doi.org/10.1016/j.jhydrol.2009.10.018

Paster, A., & Dagan, G. (2007). Mixing at the interface between two fluids in porous media: A boundary-layer solution. *Journal of Fluid Mechanics*, *584*, 455-472. https://doi.org/10.1017/S0022112007006532

Pool, M., & Carrera, J. (2011). A correction factor to account for mixing in Ghyben- Herzberg and critical pumping rate approximations of seawater intrusion in coastal aquifers. *Water Resources Research*, *47*, 1-9. https://doi.org/10.1029/2010WR010256

Pool, M., Post, V. E. A., & Simmons, C. T. (2014). Effects of tidal fluctuations on mixing and spreading in coastal aquifers: Homogeneous case. *Water Resources Research*, *50*(8), 6910-6926. https://doi.org/10.1002/2014wr015534

Pool, M., Post, V. E. A., & Simmons, C. T. (2015). Effects of tidal fluctuations and spatial heterogeneity on mixing and spreading in spatially heterogeneous coastal aquifers. *Water Resources Research*, *51*(3), 1570-1585. https://doi.org/10.1002/2014WR016068

Price, R. M., Top, Z., Happell, J. D., & Swart, P. K. (2003). Use of tritium and helium to define groundwater flow conditions in Everglades National Park. *Water Resources Research*, *39*(9), 1-12. https://doi.org/10.1029/2002WR001929

Rezaei, M., Sanz, E., Raeisi, E., & Ayora, C. (2005). Reactive transport modeling of calcite dissolution in the fresh-salt water mixing zone. *Journal of Hydrology*, *311*, 282-298. https://doi.org/10.1016/j.jhydrol.2004.12.017

Robinson, G., Hamill, G. A., & Ahmed, A. A. (2015). Automated image analysis for experimental investigations of salt water intrusion in coastal aquifers. *Journal of Hydrology*, *530*, 350-360. https://doi.org/10.1016/j.jhydrol.2015.09.046

Sanford, W. E., & Konikow, L. F. (1989). Simulation of calcite dissolution and porosity changes in saltwater mixing zones in coastal aquifers. *Water Resources Research*, *25*(4), 655-667.

Slomp, C. P., & Cappellen, P. V. (2004). Nutrient inputs to the coastal ocean through submarine groundwater discharge: Controls and potential impact. *Journal of Hydrology*, *295*, 64-86. https://doi.org/10.1016/j.jhydrol.2004.02.018

Spiteri, C., Slomp, C. P., Charette, M. A., Tuncay, K., & Meile, C. (2008). Flow and nutrient dynamics in a subterranean estuary (Waquoit Bay, MA, USA): Field data and reactive transport modeling. *Geochimica et Cosmochimica Acta*, *72*(3), 3398-3412. https://doi.org/10.1016/j.gca.2008.04.027

Spiteri, C., Slomp, C. P., Tuncay, K., & Meile, C. (2008). Modeling biogeochemical processes in subterranean estuaries: Effect of flow dynamics and redox conditions on submarine groundwater discharge of nutrients. *Water Resources Research*, *44*, 1-18. https://doi.org/10.1029/2007WR006071

Valiela, I., Costa, J., Foreman, K., Teal, J. M., Howes, B., & Aubrey, D. (1990). Transport of groundwater-borne nutrients from watersheds and their effects on coastal waters. *Biogeochemistry*, *10*(7418), 177-197.

Villermaux, E. (2012). Mixing by porous media. *Comptes Rendus - Mecanique*, *340*(11-12), 933-943. https://doi.org/10.1016/j.crme.2012.10.042

Villermaux, E. (2019). Mixing versus stirring. *Annual Reviews*, *51*(August 2018), 245-273.

Weiss, J. (1991). The dynamics of entrophy transfer in two-dimensional hydrodynamics. *Physica D: Nonlinear Phenomena*, *48*(2-3), 273-294. https://doi.org/10.1016/0167-2789(91)90088-Q

Werner, A. D., Bakker, M., Post, V. E. A., Vandenbohede, A., Lu, C., Ataie-Ashtiani, B., et al. (2012). Seawater intrusion processes, investigation and management: Recent advances and future challenges. *Advances in water resources*. https://doi.org/10.1016/j.advwatres.2012.03.004

Yoshihiro, O., Tamio, T., Hisashi, H., Atsushi, S., & Takanori, W. (2010). Optical measurement of the salinity distribution by saltwater intrusion experiment. *Japan Society of Civil Engineers*, *67*(2), 186-197.

## References From the Supporting Information

Abarca, E., & Clement, T. P. (2009). A novel approach for characterizing the mixing zone of a saltwater wedge. *Geophysical Research Letters*, *36*(January), 1-5. https://doi.org/10.1029/2008GL036995

Abarca, E., Karam, H., Hemond, H. F., & Harvey, C. F. (2013). Transient groundwater dynamics in a coastal aquifer: The effects of tides, the lunar cycle, and the beach profile. *Water Resources Research*, *49*(May), 2473-2488. https://doi.org/10.1002/wrcr.20075

Heiss, J. W., & Michael, H. A. (2014). Saltwater-freshwater mixing dynamics in a sandy beach aquifer over tidal, spring-neap, and seasonal cycles. *Water Resources Research*, *50*(8), 6747-6766. https://doi.org/10.1002/2014wr015574

Masahiro, T., Momii, K., & Luyun Jr, R. (2018). Laboratory scale investigation of dispersion effects on saltwater movement due to curoff wall installation. In *25th salt water intrusion meeting* (pp. 323-328).

Neuman, P. (1990). Universal scaling of hydraulic conductivities and dispersivities in geologic media. *Water Resources Research*, *26*(8), 1749-1758.
Nick, H. M., Raoof, A., Centler, F., Thullner, M., & Regnier, P. (2013). Reactive dispersive contaminant transport in coastal aquifers: Numerical simulation of a reactive Henry problem. *Journal of Contaminant Hydrology*, *145*, 90–104. https://doi.org/10.1016/j.jconhyd.2012.12.005

Paster, A., Dagan, G., & Guttman, J. (2006). The salt-water body in the Northern part of Yarkon-Tanimim aquifer: Field data analysis, conceptual model and prediction. *Journal of Hydrology*, *323*, 154–167. https://doi.org/10.1016/j.jhydrol.2005.08.018