Denitrification hotspots in intertidal mixing zones linked to geologic heterogeneity

The mixing between fresh and saline groundwater in beach aquifers promotes biogeochemical transformations that affect nutrient fluxes to the coastal ocean. We performed variable-density groundwater flow and reactive transport simulations with geostatistical representations of sedimentary structure to understand the influence of heterogeneity on groundwater dynamics and denitrification in intertidal mixing zones. Ensemble-averaged simulation results show that heterogeneity can enhance mixing between fresh and saline groundwater and increase residence time, resulting in up to 80% higher nitrate removal relative to equivalent effective homogeneous aquifer sediment. Denitrification hotspots form in high permeability structures where DOC and nitrate are readily supplied by convergent flow. The results provide a physical explanation for the formation of denitrification hotspots observed in beach aquifers and illustrate for the first time the influence of sediment heterogeneity on rates and spatial patterns of biogeochemical processes in intertidal aquifers that are critical mediators of land-sea solute fluxes along world coastlines.


Introduction
The intertidal zone is an energetic environment where interactions between physical, chemical, and biological processes lead to the exchange of fluid, chemicals, heat, and organic matter across the aquifer-ocean interface. Fluxes of these solutes in submarine groundwater discharge (SGD) can degrade water quality and constitute a large portion of global nutrient loads to the ocean (Burnett et al 2003). Tide and wave action across the beach are responsible for contributing large quantities of seawater and various constituents to the beach aquifer (Li et al 1999, Robinson et al 2007b, Santos et al 2009, Mcallister et al 2015. These constituents mix with solutes in terrestrially-derived fresh groundwater in the aquifer to drive redox and other reactions that alter groundwater chemistry, including concentrations of nitrate, phosphate, sulfate, iron, manganese, and uranium Sholkovitz 2002, Santoro 2010; see review in Robinson et al 2018).
In order to quantify and manage fluxes across the land-sea interface, it is critical to understand the fate of nutrients along groundwater discharge pathways leading to the coastal ocean.
The physical hydrological controls on beach groundwater flow processes have been well studied. Field and numerical modeling studies of beach aquifers in the past decade have demonstrated the prevalence of a region of elevated salinity within the beach aquifer that overlies a zone of discharging fresh groundwater (figure 1). This saltwater circulation cell is characterized by its distinct flow pattern whereby seawater infiltrates across the upper beachface and circulates downward until discharging near the low tide mark. The effects of hydrologic controls, such as tides (Robinson et al 2007b, Boufadel et al 2011, waves (Robinson et al 2006, Sous et al 2016, Geng et al 2017, and terrestrial fresh groundwater discharge (Robinson et al 2007c, Kuan et al 2019, on flow and transport indicates that flowpaths and salinity distribution in the circulation cell depends on tidal stage, beach slope, freshwater forcing, tidal amplitude, hydraulic conductivity (K), and the terrestrial freshwater hydraulic gradient. Flow and solute transport in beach aquifers, and resulting SGD patterns, can be highly dynamic in response to transient hydrologic conditions (Robinson et al 2007a, Greskowiak 2014. These physical hydrological and morphological factors can affect nitrogen transformations along coastal groundwater flowpaths (Anwar et al, 2014. Field observations have identified hotspots of nitrate processing within the intertidal zone (Schutte et al 2015, Kim et al 2017. While some of these hotspots correspond well with expected location based on salinity and solute distributions, others do not. While these hotspots could be explained by geochemical heterogeneity, we explore other physical factors that could also contribute to unanticipated distributions of beach aquifer reactivity. Coastal groundwater flowpaths and fluid and chemical fluxes to surface water are strongly affected by geologic heterogeneity (Mulligan et al 2007, Bokuniewicz et al 2008, Russoniello et al 2013. Preferential flow though high-K conduits at the shelfscale creates complex salinity distributions and spatially variable SGD patterns (Michael et al 2016). At the shoreline, high-K sediments act as preferred pathways of fresh groundwater and nutrient transport (Mulligan et al 2007, Russoniello et al 2013, and heterogeneity of shallow bay sediments can control nitrate flux distributions and redox chemistry of SGD (Sawyer et al 2014). Because geologic heterogeneity has a strong control on flow (e.g. Li et al Li, et al., 2016, Trglavcnik et al 2018 and has been shown to greatly enhance mixing (Pool et al 2015, Michael et al 2016, and initiate denitrification (Sawyer 2015), it has the potential to substantially modify nutrient cycling and other biogeochemical processes by modifying intertidal flow patterns (figure 1), with important implications for chemical ocean budgets.
In this study we investigate the influence of geologic heterogeneity on flow patterns, mixing, and denitrification in sandy beach aquifers under tidal influence. Variable-density groundwater flow and reactive transport models are used with geostatistical realizations of heterogeneous sediments to explore linkages between preferential flow paths resulting from geologic heterogeneity, and the formation and distribution of denitrification hotspots in beach aquifers. The models assume that denitrification is the primary mechanism of N removal and does not incorporate all N transformation pathways and speciation. We compare attenuation of terrestrially-derived nitrate in ensemble-averaged heterogeneous models to effective equivalent homogeneous models.

Geostatistical and numerical modeling
Geologic heterogeneity exists in beaches due to sediment source-switching and deposition, barrier island migration into marsh sediment, and re-working of intertidal sediments (e.g. Rine and Ginsburg 1985, Woodroffe 2002, Holland and Elmore 2008. Ten geologic realizations of sand and silt were generated geostatistically using the sequential Gaussian simulation algorithm in SGeMS (Remy et al 2009). The horizontal and vertical variogram model ranges were 10 m and 1.0 m, respectively, to capture the spatial correlation of sand and silt facies (figure S1 and additional details in Text S1(https://stacks.iop.org/ERL/15/084015/mmedia)). For each realization, the silt fraction was varied from 0.1 to 0.9 in increments of 0.1 (figure 2) by truncation of the Gaussian random fields using threshold values corresponding to proportions of sand and silt (Beucher and Renard 2016), totaling 90 K fields. Representative values of K were assigned to sand (1.0 × 10 −4 m s −1 ) and silt (5.0 × 10 −6 m s −1 ), representing a twentyfold difference in K (following Sawyer 2015).
Equivalent effective hydraulic conductivity values were calculated for the 9 silt fractions for each realization by simulating groundwater flow vertically and horizontally through the model domain. Darcy's law was then used with the simulated fluxes to calculate effective K (table S1), and the average of the effective horizontal and vertical K values for each silt fraction were used in the homogeneous models. The geostatistical and numerical models represent a 2-D shoreperpendicular cross-section of a beach aquifer with a shoreface slope of 0.1 (figure 1).
Variable-density groundwater flow and salt transport in the heterogeneous and equivalent homogeneous permeability fields were simulated with SEA-WAT v4.0 (Langevin et al 2008). A uniform sinusoidal time-varying hydraulic head was implemented at the shoreface boundary in all models using the Periodic Boundary Condition package to simulate tidal forcing (Post 2011) where h (m) is the tidal elevation at time t, A (1.0 m) is the tidal amplitude, and ω (12.567 rad d −1 ; 2π T ) is the tidal angular frequency. Along the shoreface a constant seawater concentration (35 ppt) was assigned for inflow and a zero concentration gradient was assigned for outflow. A constant head of 1.0 m was assigned to the landward vertical boundary for the heterogeneous models, and the simulated freshwater fluxes were assigned to that boundary for each corresponding equivalent homogeneous model. Reactive transport of labile dissolved organic carbon (DOC), dissolved oxygen (O 2 ), and nitrate (NO 3 − ) was simulated using the phase-resolved flow fields with PHT3D v2.13 (Prommer and Post 2002). The chemical reaction network for aerobic respiration and denitrification was adopted from Bardini et al (2012) and are based on the redox kinetics of Hunter et al (1998). DOC oxidation is assumed a linear first-order reaction with respect to the DOC concentration, with the reduction of O 2 and NO 3 − proceeding from the most (O 2 ) to least (NO 3 − ) thermodynamically favorable electron acceptor: where β i is the stoichiometric ratio between moles of transferred electrons per mole of oxidized DOC and the moles of electrons per mole of electron accepter. f i describes a modified linear Monod behavior where the reduction rate is inhibited by the availability of the electron acceptor: with f 0 = 0 and or where α i is a dimensionless parameter that considers the inhabitation of the reduction rate, r red, i due to the availability of the i-th electron acceptor, C i and C i, lim are respectively the molar concentration and the limiting molar concentration of oxygen and nitrate. If the oxygen or nitrate concentration is greater than their limiting concentration, then the reduction rate is independent of their respective concentration, while the reduction rate is linearly proportional to the concentration if the concentration falls below the limiting concentration. We chose to not include nitrification because mineralization of marine-derived DOC and subsequent nitrification represents a recycled marine source of nitrogen, while nitrate in the fresh groundwater endmember represents a new nitrogen source to coastal surface water ecosystems. The models are therefore most representative of systems where groundwater ammonium concentrations are low (e.g. Mcallister et al 2015, Couturier et al 2017. The source concentrations of DOC, O 2 , and NO 3 − at the shoreface boundary were specified according to values previously measured coastal sediments (Berner 1981, Van Cappellen andWang 1996). Seawater inflowing across the shoreface boundary contained DOC (2 × 10 -3 M) and O 2 (3.13 × 10 -4 M), and inflowing fresh groundwater along the left vertical boundary represented a terrestrial source of NO 3 − (1.29 × 10 -4 M) and was anoxic (0 M O 2 ). Model parameters are within the range of published values (table S2).
Two sets of flow and reactive transport simulations were performed. The first set incorporated the 90 heterogeneous permeability fields and the second set included the 9 corresponding equivalent effective homogeneous sediments. Simulations were run to 2000 d to achieve dynamic steady-state with respect to hydraulic heads and concentrations. Tidally-averaged results from the last tidal cycle on day 2000 are reported. Nitrate removal efficiency was computed as the difference between influx across the left boundary and outflux across the seabed divided by the influx.

Heterogeneity effects on intertidal salinity and nitrate discharge zone
Simulated intertidal salinity distributions in heterogeneous sediments for all silt fractions are more complex than in the equivalent homogeneous models; examples are shown in figure 2. The salinity distributions in the homogeneous cases are consistent with previous modeling studies of tidal effects on beach aquifers, showing the formation an intertidal saltwater circulation cell and fresh discharge zone landward of the traditional saltwater interface (Vandenbohede and Lebbe 2005, Abarca et al 2013, Greskowiak 2014. Salinity in the equivalent effective homogeneous fields exhibits smooth transitions between the fresh and saline groundwater with concentration gradients that shift gradually along the perimeter of the circulation cell and along the lower saltwater interface. However, salinity distributions for the heterogeneous model cases differ from those in equivalent effective homogeneous fields (figure 2). The salinity contours along the intertidal mixing zone and lower interface are more irregular and concentration gradients are shaper in some regions and more diffuse in others because high groundwater flow rates within sandy conduits transport groundwater of differing salinity into adjacent and less mobile water masses.
The size of the intertidal mixing zone in the heterogeneous models increases with increasing silt fraction and is consistently larger than in the equivalent effective homogeneous models (figures 2 and S2-S3). In contrast, the size of the mixing zone for the homogeneous models decreases from 200 m 2 to 150 m 2 for silt fractions between 0.1 and 0.5 and then increases to 423 m 2 between silt fractions 0.5 and 0.9. The decrease in size is the result of an increase in anisotropy (equivalent horizontal to vertical hydraulic conductivity; figure S4) from low to intermediate silt fractions; equivalent vertical hydraulic conductivity decreases as more horizontal silt structures are added to the system, which results in lower vertical effective hydraulic conductivities in the equivalent homogeneous models and reduces seawater infiltration. Beyond a silt fraction of 0.5, the lithology is reversed and silt comprises a majority of the sediment and hence any additional silt results in less anisotropy, allowing for greater infiltration relative to the opposing terrestrial freshwater input. The salinity distributions for the heterogeneous fields are largely unaffected by anisotropy. The ensemble-averaged size of the mixing zone in the heterogeneous models increases with silt fraction from 228 m 2 to 564 m 2 due to greater mixing (figures 2 and S2). Thus, heterogeneity is a more important control than anisotropy on mixing zone size and salinity patterns in beach aquifers.
Geologic structure strongly influences the spatial variability of nitrate fluxes across the seabed ( figure  S5). Nitrate fluxes for all heterogeneous and homogeneous models are focused in the fresh discharge zone, however nitrate fluxes in the heterogeneous models exhibit greater spatial variability, with fluxes differing up to an order of magnitude over a distance of 30 cm across the sediment-water interface.

Nitrate removal efficiency
Ensemble-averaged model results show that removal efficiency can be 80% higher in heterogeneous compared to homogeneous sediments due to enhanced mixing between nitrate-enriched fresh groundwater and DOC-enriched seawater (figure 3). The largest difference occurs at 0.5 silt fraction, where the size of the mixing zone in the heterogeneous models relative to the equivalent homogeneous model is largest, as discussed above.
Removal efficiency consistently increases from 0.0-1.0 silt fraction for both the heterogeneous and equivalent homogeneous models, however the rate of increase is greater for silt fractions above 0.5 for both sets of models. The trend in removal efficiency is due to the relative balance between saltwater-freshwater mixing intensity and residence time in the circulation cell. In heterogeneous sediment, the size of the mixing zone increases on average at a rate of 44 m 2 per 0.1 increase in silt fraction (figure S2) and the residence time increases rapidly on average from 62-561 d between 0.5 and 1.0 silt fraction (figure S6). The combination of a steadily growing mixing zone with increasing silt fraction and substantially longer residence times from 0.5-1.0 silt fraction result in more complete nitrate removal at higher silt fractions. In the equivalent homogeneous sediment, the size of the mixing zone decreases from 0.0-0.5 silt fraction, but residence times increase at approximately the same rate as the heterogeneous models. The effect of longer residence times in maintaining contact between DOC and nitrate is more important than the intensity of mixing within this lower range of silt fractions, resulting in a gradual increase in removal efficiency. Thus, removal efficiency increases more rapidly in heterogeneous sediment at lower silt fractions due to the combination of an increase in mixing zone size and longer residence times. These findings are consistent across a range of beach slopes (figure S7). Total denitrification decreases with increasing silt content in heterogeneous and homogeneous sediment because solute supply is limited in low-permeability Figure 2. Example heterogeneous fields and simulated salinity, nitrate, and dissolved oxygen concentrations, and equivalent homogeneous models. White contours are 10% and 90% seawater salinity. sediments ( figure 3(b)). The ensemble-averaged results indicate that the greatest benefit of heterogeneity occurs at moderate proportions of sand and silt.

Reactivity distribution
Heterogeneous sediment in the intertidal mixing zone promotes the formation of denitrification hotspots ( figure 4(a)). The hotspots are maintained by high fluxes of nitrate and DOC in sand, while low permeability silt facies restrict solute transport and are thus largely devoid of denitrification with the exception of small reactive areas in silt at sand-silt interfaces. Close to the sand surface and in the interior of the circulation cell where salinity is >90% seawater, aerobic respiration consumes dissolved oxygen along circulating flow paths. The landward fringe of the mixing zone hosts a denitrification hotspot that is confined vertically by the thickness of the sand facies in which it is located ( figure 4(a)). Denitrification within this first zone reduces nitrate concentrations as circulating flow transports pore water from the hotspot into the underlying silt facies where low nitrate concentrations inhibit denitrification. Farther along circulating flow paths beyond the silt facies, a horizontal conduit of high-K sediment transports nitrate-rich fresh groundwater to the base of the mixing zone. Here, contact between DOC in saline pore water and nitrate in fresh pore water form a second denitrification hotspot. Nitrate is fully removed along the circulating flow path until upward vertical flow of freshwater along the lower saltwater-freshwater interface provides a new source of nitrate, forming a third hotspot. Unlike the equivalent effective homogeneous models in which denitrification rates transition gradually along the mixing zone beneath the high tide line ( figure 4(b)), denitrification can The removal difference is the difference between total nitrate removed in the homogeneous models and the average of the total removed in the heterogeneous models.
occur as hotspots at the confluence of fresh and saline groundwater in high-K sand facies along the margin of the circulation cell. The results demonstrate that geologic heterogeneity in beach aquifers is a major control on the formation and distribution of denitrification hotspots by regulating the supply of reactive solutes and their contact time.

Discussion
While there are a growing number of combined field and modeling studies of groundwater dynamics in sandy beaches, few have considered the role of geologic structure on flow and transport (e.g. Trglavcnik et al 2018). Consequently, homogeneity is often assumed in groundwater models despite field observations suggesting strong hydrogeological influence on salinity distributions. Here, we show that geologic heterogeneity can lead to the complex salinity and flow patterns widely observed in beach aquifers (Robinson et al 2006, Abarca et al 2013. Similarly, SGD and associated nutrient fluxes to receiving coastal ecosystems are spatially heterogeneous (Kroeger et al 2007, Sawyer et al 2014. Previous field measurements along the interfluve of a paleochannel have shown that nitrate fluxes can vary by an order of magnitude over tens of meters (Sawyer et al 2014). The model results in this study demonstrate that ubiquitous heterogeneity on a smaller scale is also a strong factor affecting the spatial distribution of nitrate fluxes across the seabed, varying by an order of magnitude over a distance of 30 cm. In models with heterogeneity, nitrate removal increased despite a decrease in sand connectivity at higher silt fractions. While lower connectivity of sand reduces DOC and nitrate contact, higher silt fractions increase the area of the mixing zone and promote denitrification. The effect of mixing zone area on nitrate removal is thus more important than connectivity of high-K sediments in heterogeneous beach aquifers.
We demonstrate that geologic heterogeneity in beach aquifers can form biogeochemical hotspots within high-K zones that hydrologically connect terrestrially-derived fresh and marine-derived saline groundwater, providing a physical explanation for observed beach aquifer denitrification hotspots (e.g. Schutte et al 2015, Kim et al 2017. In other environments, such as the hyporheic zone, heterogeneity in pore throat size can create anoxic microzones owing to variability in residence time (Briggs et al 2015). At the centimeter scale, denitrification in streambed sediments is most efficient in low permeability zones where residence time and contact time between reactive solutes is longest (Flewelling et al 2012). However, flow processes in beach aquifers differ considerably from those in streambed sediments. For many coastlines, beach aquifers represent the terminal end of groundwater flow paths leading to the ocean, where interaction between waves, tides, and terrestrial fresh groundwater promote reactivity that attenuate nutrients, arsenic, and other contaminants that would otherwise discharge to coastal surface water ecosystems (Anwar et al 2014. The results of this study provide new insights into the interaction between geology and hydrologic processes within these important land-sea transition zones, and the effects of those interactions on spatial patterns of reactivity. Geochemical controls may also influence hotspot formation, such as DOC filtration dynamics (e.g. Kim et al 2019) and sediment geochemical heterogeneity (e.g. Sawyer 2015). While we did not consider organic matter filtration in this study, physical heterogeneity may also play a role in this process, as organic particles may become differentially entrapped in different sediment. Sawyer (2015) demonstrated that low permeability carbon aggregates (i.e. chemical heterogeneity) in the hyporheic zone can enhance local reactivity by serving as an in-situ organic carbon source. In our study, denitrification was greater in high-permeability sediment where nitrate and DOC were supplied by the flow system, since the reaction was dependent on the mixing between solutes, rather than contact between solutes and sediments. The reactive hotspots shown in figure 4 would likely be larger if geochemical heterogeneity were considered because hotspots formed by physical heterogeneities would merge with reaction zones supported by chemical heterogeneity. This in turn would enhance total removal, suggesting that our estimates are conservative. The localized zones of reactivity demonstrated in our study even in the absence of geochemical heterogeneity suggest that both physical and geochemical conditions are primary controls on biogeochemical reactions within beach aquifers. The results demonstrate that mixing-dependent reactions can be substantially affected by relatively mild (twentyfold) heterogeneity in hydraulic conductivity that is likely typical of beach aquifers, and are therefore applicable to a range of sandy coastlines worldwide.
Beach nourishment and scraping are likely to affect nitrate attenuation in the intertidal circulation cell. There is often a poor match between sediment size of the pre-existing beach and of the sand used for nourishment (Dean and Campbell 1999). In beach scraping, the coarsest sediment at the low tide line is transported to the backshore where sediments are finer (Stauble 2005). Thus, beach nourishment and scraping can form a geologically heterogeneous intertidal subsurface and likely will affect attenuation of nitrate in through-flowing fresh groundwater. The new understanding of the role of geologic structure on nitrate fluxes to coastal surface water ecosystems should be considered in beach management decisions in order to maximize the benefits or mitigate adverse effects of beach nourishment, scraping, and other artificial interventions that reshape the beach.
In addition to enhanced nitrate attenuation, beach heterogeneity affects the spatial distribution of nutrient fluxes from terrestrial aquifers to coastal surface water ecosystems. Nutrient loading via SGD is a key control on benthic habitat (Knee and Paytan 2011), thus the spatial variability of groundwater nitrate fluxes demonstrated in this study is likely to impact local benthic and seabed microbial communities. In addition, distributed pockets of fresh and saline SGD zones may influence other seabed organisms, as benthic invertebrates are sensitive to nearbed salinity (Metaxas and Young 1998). The effects of seabed geology on the spatial distribution of SGD and concentrations of associated constituents should therefore be considered in benthic habitat mapping in the shallow subtidal zone.
The model assumed saturated flow, but if unsaturated flow was considered the longer residence times in the intertidal zone (e.g. Geng et al 2014) could increase oxygen consumption of percolating seawater, promoting anoxia and denitrification at depth. However, oxygenation of shallow groundwater can occur due to air entrapment with a fluctuating water table (Williams and Oostrom 2000), inhibiting denitrification. We did not include DOC mineralization and nitrification because the objective of this study was to understand the effects of geologic heterogeneity on denitrification of land-derive nitrate, which is a major control of productivity. Nitrification would promote denitrification by consuming oxygen and producing nitrate, hence our estimates may underestimate total denitrification considering marine-and terrestrially-derived nitrate sources. The models were also based on 2-D heterogeneous permeability fields to maintain practical computational times and did not consider 3-D preferential transport paths in the along-shore direction (e.g. Siena and Riva 2018), which could enhance mixing and reactivity.

Conclusion
Simulations of groundwater flow and reactive transport in beach aquifers show that geologic heterogeneity can enhance mixing between fresh and saline groundwater, increase residence times, and create preferential transport pathways. These effects on physical flow can result in up to 80% higher nitrate removal relative to equivalent effective homogeneous aquifers for the model setup, reaction network, and geologic realizations considered. Reactive hotspots form within high-permeability sediments in the intertidal mixing zone where nitrate and DOC are supplied to maintain high denitrification rates. Geologic heterogeneity can result in order-of-magnitude spatial variability in nitrate fluxes at the submeterscale. The impact of heterogeneity on mixing and nitrate attenuation will likely vary with tidal forcing, pore water geochemistry, the scale of heterogeneity, and the extent of geologic continuity. The results of this study are thus an approximation of the role of geologic heterogeneity on affecting nutrient cycling in these systems. The results suggest that geologic heterogeneity may be responsible for the denitrification hotspots and complex saltwater-freshwater mixing zones observed in beach aquifers along a variety of coastlines. The findings have important implications for quantifying nutrient loads to coastal surface water ecosystems and indicate that geologic heterogeneity is an important consideration when interpreting geochemical signatures and biogeochemical reactivity in coastal aquifers.

Data Availability
The data that support the findings of this study are available from the corresponding author upon request.