Subsurface Mixing Dynamics Across the Salt-Freshwater Interface

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.


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;, 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).
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.
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;. 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., 2014Pool et al., , 2015 as well as heterogeneity (Abarca Cameo, 2006;Kerrou & Renard, 2010; 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;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.

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 where q is the specific discharge, K is the hydraulic conductivity, h f the equivalent freshwater head, ρ the fluid density, ρ f the density of freshwater, and e z is the unit vector in y direction. Fluid mass conservation in the absence of sources and sinks implies ∇⋅ρq = 0. The fluid density is assumed to be linearly dependent on the salt mass fraction ω (mass of salt dissolved per unit mass of fluid) given by ρ = ρ f [1 + ϵ′c], where ϵ′ is the buoyancy factor given by ϵ′ = (ρ s − ρ f )/ρ f with ρ s the density of seawater and c is the normalized salt concentration defined as c = ω/ω s with ω s the salt mass fraction of seawater. The concentration c evolves according to the advection dispersion equation, which in the steady state reads as where D is the dispersion tensor (Bear, 1988), D m is the molecular diffusion coefficient, and ϕ 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, , which compares the advection and dispersion times, and , which compares the advection and diffusion times, and the gravity number, , which compares the viscous q f /K and buoyancy forces ϵ′ (see Supplementary Information) (see Abarca et al., 2007), where b defines the domain thickness, α t is the transverse dispersivity, q f is the specified fresh water flux, and ϕ is the porosity.

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 . 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 and 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 × 10 −2 m/d to 3 × 10 −2 m/d. Thus, the simulated scenarios are characterized by a of 500 and ranging between 17.3 and 7.2. Since we vary only the flow rate, the range of considered is equivalent to the one of . Therefore, in the following, the scenarios are characterized by their 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).

Mixing Measures
The variability of mixing along the SFI can be characterized by the local scalar dissipation rate, which is defined using 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, where ‖⋅‖ denotes the L 2 -norm. The salt concentration gradient at the SFI can be approximated by θ ≈ c s /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, Θ ζ , 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. 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.

Mixing Mechanisms and Mixing Model
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.

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 by changing the freshwater flow rate. The local mixing widths along the interface for different 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 (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, ∝ , while s m decreases as ∝ 1∕2 (see Supporting Information S1). Figure 2c shows the evolution of the concentration gradient θ along the interface for different . All θ collapse on a single curve when rescaled with their respective minima θ 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 θ/θ m can be well represented by the inverse interface width ( ∕ ) −1 . We observe the same behavior for the mixing rates (see Figure 2d), which are rescaled by their minima χ m . Their evolution is well represented by χ ≈ α t vθ 2 normalized by its minimum. This highlights the central role of the interface width on mixing along the interface.

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), where γ is the stretching (or compression) rate and D t /s 2 is the dispersive expansion rate with D t = D m + α t v 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γ where = 2 0 ∕ (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 γ, For times larger than t s , the mixing width given by Equation 7 is expected to converge to the Batchelor scale = √ ∕ . 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 = √ ∕ . 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 ≪ z m , expression (Equation 7) implies where we set the transverse dispersion coefficient D t ≈ α t v a . This explains the increase in the mixing width observed in Figures 2a and 2b. The dependence of s on α t is confirmed by additional numerical simulations for variable α t (see Supporting Information S1). For increasing distance along the interface, the acceleration and thus v and γ increase notably along the interface. Assuming that v and γ 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 where v(z) and γ(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 γ is given by the derivative of the flow velocity v(z) along the interface, γ(z) = dv(z)/dz. Thus, we obtain for the interface width in terms of v(z) the expression 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 γ. The velocity along the interface can be approximated by v(z) = q f b/ξ(z), where ξ(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 (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) 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 ξ(z) into (11), we obtain the compact expression The analytical solution explains the scaling behavior of s observed in Figure 2b. Note that the Glover solution predicts the toe length = ′ ∕2 . In fact, we can write Equation 12 as 7 of 10 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 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 ( ) = ∕ √ (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 . 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, γ 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 α 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.

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.