Bed form-induced hyporheic exchange and geochemical hotspots

This study follows the methodology of [7] to create synthetic 3D bed forms for use in coupled numerical surface water-groundwater flow and reactive transport models. Bed form morphology is generated using the geometric model from a widely used set of codes originally developed by Rubin [39] and later published by Rubin and Carter [40] as a set of MATLAB routines. More than eighty parameters are used in their model, which generates geomorphologically plausible, realistic, and visually satisfying bed form morphologies that represent a wide range of dynamic environmental conditions (see examples in Figure 1). The focus is on relating the depositional and geomorphological processes and resulting characteristics to different scales of bed forms. Specifically, the models synthesize the bed form geometry by merging two sets of surfaces generated using sine functions. Here, nine representative bed forms are used to represent typical types of sand dunes under unidirectional river flow, such as barchanoid or linguoid forms with dune crest lines that are sinusoidal and with pronounced lobes and saddles.

The mathematical formulation for reproducing the bed form geometry may be written as

where $m_{i}$ is defined as:

The bed form topographies (i.e., elevations) in the x-y plane are calculated as the maximum value of each set as:

where $i=1,2$ represent the first and second set of bed forms, respectively. The x and y are the Cartesian coordinates, and $L$ is the length of the repeating block. $A_{i}^{f}$ and $A_{i}^{s}$ are the planform amplitudes of the first and second sine curves of each set of bed forms, respectively. The $\varphi_{i}^{bf}$ is the bed form phase (in degrees), $\eta$ is the mean steepness (height-to-wavelength ratio), $\lambda_{L}$ and $\lambda_{T}$ are the longitudinal and transverse wavelengths of the bed form, respectively, and $\alpha_{i}$ and $\beta_{i}$ are the phases of the first and second planform sinuosity of each set of bed forms, respectively.

In this study, bed form shapes are specified by two phase shift parameters:

with $\beta_{2}=\beta_{1}$. Figure 1 and Table 1 present details of the nine different bed forms. Other parameters are held constant throughout the computations (Table 2). The overall bed form geometry changes based on variations in the phase shift parameters (Figure 1).For instance, a $\Delta\varphi_{1}$ value of $90^{\circ}$ represents linguoid bed forms while a $\Delta\varphi_{1}$ value of $270^{\circ}$ represents lunate (crescent-shaped) bed forms. There are also more subtle differences between bed forms that are better observed in 2D cross-sections. For example, while linguoid ripples produce an angle to the flow and have a random shape rather than a W-type shape, the stoss sides of lunate ripples are curved rather than having a lee slope. Photos of such bed forms exposed in aquatic environments are shown in Figure 2.

A modeling framework named SimSGI (simulation of surface water and groundwater interactions) is developed to more realistically simulate bed form-driven hyporheic exchange processes. SimSGI integrates two advanced models of flow and mass transport: a surface water solver (OpenFOAM) [43] and a subsurface simulator (PFLOTRAN) [44, 45, 46]. SimSGI also uses LaGriT [47], an advanced and open-source mesh generation tool, to create numerical meshes for the different bed form geometries discussed above. After the streambed topographies are geometrically simulated, LaGrit generates computational grids based on the resulting surfaces and the SimSGI framework provides interface scripts for using them in both the surface water (OpenFOAM) and groundwater (PFLOTRAN) domains (Figure 3). After solving the turbulent surface flow in OpenFOAM, the resulting pressure distribution over the bed surface is assigned as the top boundary condition for the groundwater domain, which is then used to solve the groundwater flow and reactive solute transport in PFLOTRAN.

Turbulent surface water flow generates non-uniform pressure distributions across the sediment-water interface in response to bed form shape, with more complex patterns associated with increasing $\Delta\varphi_{2}$. For example, as $\Delta\varphi_{2}$ increases for both linguoid ($\Delta\varphi_{1}$ = $90^{\circ}$) and lunate ($\Delta\varphi_{1}$ = $270^{\circ}$) bed forms, the flow regime shifts as adjacent bed forms become staggered, resulting in larger high-pressure zones on the lee faces of the bed forms. In-phase ($\Delta\varphi_{2}$ = $0^{\circ}$) and out-of-phase ($\Delta\varphi_{2}$ = $180^{\circ}$) bed forms display symmetrical pressure distributions with aligned low- and high-pressure zones, while the pressure distribution across intermediate bed forms ($\Delta\varphi_{2}$ = $90^{\circ}$) is more erratic. In contrast, as $\Delta\varphi_{1}$ increases the high and low pressure zones becomes more uniformly distributed across the bed form geometry, trending away from focused high pressure zones along bed form troughs that drive downwelling conditions. Spatial changes in the pressure distribution alter the pressure gradient, and thus control hyporheic exchange and the flux of solutes between the river and groundwater. Further, larger pressure zones ultimately affect hotspot formation by increasing both the volume of the reaction zone and the overall reaction rates (Table 5).

During the first few days of simulation solutes are only infiltrating until the amounts of infiltration and exfiltration are stabilized and the system reaches an equilibrium. As shown in Figures 6a and 6b, the net exchange flux of solute A is positive meaning that infiltration of solute A at the streambed is more than its exfiltration. It should be noted that characteristics of solute A and solute B are identical and thus they show the same behavior. Net flux of solute C through the streambed is negative, which indicates that solute C is mostly leaving the subsurface and going back to the surface water.

Moreover, the comparison of exchange fluxes for all bed forms in Figures 6a and 6b shows that as $\Delta\varphi_{2}$ grows, the net flux of solute A transferring from surface into subsurface, as well as the net flux of solute C leaving the sediment increase. In other words, bed forms with larger $\Delta\varphi_{2}$ transfer higher amounts of solute. In particular, for lunate bed forms (c, f, i) the absolute value of solute transferred (both solute A and solute C) is growing as $\Delta\varphi_{2}$ is increasing from bed form (c) to bed form (i). Linguoid bed forms (a, d, g) shows similar trend. This is also shown in Table 5 that higher $\Delta\varphi_{2}$ results in higher average velocity, larger reaction rate and higher amounts of solute exchange. Bed form (i) has a phase shift ($\Delta\varphi_{2}$) of 180^∘ compared to bed form (a) which leads to 21% change in average velocity and average reaction rate. The solute exchange flux is also increased by 21% for bed form (i).

Furthermore, both streambed morphology and hydraulic conductivity ($K$) influence the solute exchange amount between surface water and groundwater. Figure 7 displays the equilibrated mass flux at time = 20 days for linguiod type bed forms and with different $K$. Other bed form types such as lunate with the same $\Delta\varphi_{1}$ shows similar trends. As shown, increasing $K$ will enhance the amount of entering mass of solute A from surface water to the hyporheic zone. However, the amount of solute C moving out of the hyporheic zone and into the stream decreases. It is also clear that if sediment is more conductive the differences between linguoid bed forms are amplified in terms of solute mass flux.

The residence time of solutes infiltrating from river water to riverbed sediments can provide critical information about biogeochemical reaction times [55, 56]. In fact, it is known that residence time distributions can determine the type and rate of many processes occurring in surface, near surface, and deep environments from solely physical to complex biogeochemical processes[57]. Therefore, there have been several recent studies directed at unraveling mechanisms controlling solute residence times for hyporheic exchange processes at multiple scales [58, 59, 34, 60]. PFLOTRAN does not explicitly simulate residence time distributions using particle tracking, but it does use the concept of groundwater age [61]. Here, this capability is used to quantify residence times for infiltrated (i.e., A and B) and produced (i.e., C) solutes to further study the impacts of bed form topographies on hyporheic processes.

Bed form geometry controlled the extent of solute transformations by regulating residence times within the streambed. At steady-state (t=20 days), residence times are lowest in the shallow streambed, associated with higher groundwater velocities and shorter hyporheic flowpaths than deeper regions (Figure 8). Longer residence times (up to 14 days) in the deeper streambed allowed for more complete solute transformation, however, resulting in lower concentrations of A and B and the highest concentrations of C. As a governing factor in the dynamics of groundwater flow patterns, the phase shift in the second planform sinuosity ($\Delta\varphi_{2}$) strongly controlled the distribution of residence times based on bed form geometry. Average residence times increased from $2.08\pm 0.01$ days to $2.52\pm 0.03$ days as phase shift increased from $\Delta\varphi_{2}$=0^∘ to $\Delta\varphi_{2}$=180^∘. Though solute concentration profiles followed a similar trend for all bed form geometries, with solutes A and B decreasing to roughly 0.1 M and solute C reaching up to 0.9 M as cell residence time increased, concentrations are much more variable throughout the domain when bed forms are out-of-phase ($\Delta\varphi_{2}$=90^∘). For example, the average standard deviation of the concentration of solute C is 568% and 351% larger in bed forms where $\Delta\varphi_{2}$=90^∘ than those where $\Delta\varphi_{2}$=0^∘ and $\Delta\varphi_{2}$=180^∘, respectively.

Temporal variation of concentration can also be described as time series at a particular depth. As shown in Figure 9, at 0.2 m deep in the center of the subsurface domain, solute concentrations for out-of-phase bed forms ($\Delta\varphi_{2}\neq 0$) experience sharp increases while they gradually increase for in-phase bed forms ($\Delta\varphi_{2}=0$). This outcome agrees well with the exchange flux results in section (3.2), which demonstrated that the out-of-phase bed forms allow for a greater amount of solute exchange.

The volume of the reaction zone, $V_{rz}$, is defined where reaction rates are greater than 1% of the maximum rate (k_max). For each bed form type (e.g., linguiod versus lunate), $V_{rz}$ is most sensitive to $\Delta\varphi_{2}$, and increases for out-of-phase bed forms ($\Delta\varphi_{2}\neq 0$) (Figures 10). Table 5 also shows that with a 180^∘ phase shift, $V_{rz}$ grows by 20%.

Hydraulic properties of sediments also influence solute transport and reaction processes. Hyporheic exchange volumes are lower in less conductive sediments, and the flux of solutes into the streambed is thus decreased. As a result, bed form morphology is negligible to $V_{rz}$ because the magnitude of hyporheic exchange is below the threshold of influence. Conversely, bed form shape is a primary control on $V_{rz}$ in highly conductive sediments (K $<$ $10^{-6}$ m/s) because flow and solute exchange are not hindered by sediment permeability (Figure 11b). However, given the different flow regimes imposed by varying bed form geometries, this relationship may not be consistent in all environments. This inverse relationship also represents another layer of complexity that should be evaluated in future studies: under otherwise analogous environmental conditions (i.e., flow depth, surface water velocity), streambed sediments of varying permeability will form bed forms of different size and shape. As a result, the relative influence of bed form geometry and sediment permeability on hyporheic processes will shift across time and space as sediments are transported and the complex streambed topography evolves in response to episodic flow events.

Depending on the bed form geometry, the maximum reaction rate for identical biogeochemical processes may change in space. Figure 12a depicts the plan view of locations with maximum reaction rates in different bed forms. Figures 12b, 12c, and 12d display the locations of maximum reaction rate from a cross-sectional view. Here, bed forms (a), (e), and (h) are presented to show distinct locations of maximum reaction rate spots.