Biogeochemical and Hydrological Drivers of Heterogeneous Nutrient Exports From Subterranean Estuaries

Submarine groundwater discharge (SGD) to coastal zones contributes terrestrial freshwater and nutrients that may support harmful algal blooms (HABs). The magnitude of nutrient exports via SGD depends on volumes of fresh groundwater discharge, its chemical composition, and modifications by biogeochemical processing within subterranean estuaries. Thus, the ability to upscale SGD exports requires knowing the range of chemical composition of inland groundwater and how those compositions may be transformed as fresh and saltwater mix within subterranean estuaries. These processes may create heterogeneous magnitudes of solute exports, even at small spatial scales, and such heterogeneities have rarely been assessed for regional or global SGD nutrient export estimates. To evaluate heterogeneity in subterranean estuary processes and nutrient export, we collected seasonal pore water samples in 2015–2016 at three proximal (<20 km) subterranean estuary sites in Indian River Lagoon, FL. Sites have homogenous hydrogeological settings, but differ in land use and coastal features, and include a mangrove site, an urban site, and a site offshore of a natural wetland. All sites exhibit little variation through time in nutrient concentrations and modeled SGD rates. In contrast, each site exhibits significantly different nutrient concentrations of potential fresh groundwater sources, fresh groundwater discharge volumes, and nutrient transformations within subterranean estuaries. Groundwater specific discharge correlates with nutrient concentrations, suggesting that higher residence times in the subterranean estuary increase biogeochemical transformations that reduce anthropogenic nutrient loads but increase in situ nutrient sources derived from organic matter remineralization. The differences in transformations lead to SGD nutrient contributions that differ by orders of magnitude between sites and have N:P ratios that are greater than the Redfield ratio (15) for the mangrove (29) and urban sites (28), but less than the Redfield ratio for the wetland site (8). These results indicate that heterogeneity of both absolute and relative nutrient export via SGD complicates integration of nutrient fluxes across regional coastal zones and evaluations of its impacts to coastal ecosystems. A better understanding of the drivers of heterogeneity, including subterranean estuary processes, land use, coastal topography, and vegetation dynamics could improve assessments of regional nutrient loading and upscaling for estimates of global solute cycles.

Submarine groundwater discharge (SGD) to coastal zones contributes terrestrial freshwater and nutrients that may support harmful algal blooms (HABs). The magnitude of nutrient exports via SGD depends on volumes of fresh groundwater discharge, its chemical composition, and modifications by biogeochemical processing within subterranean estuaries. Thus, the ability to upscale SGD exports requires knowing the range of chemical composition of inland groundwater and how those compositions may be transformed as fresh and saltwater mix within subterranean estuaries. These processes may create heterogeneous magnitudes of solute exports, even at small spatial scales, and such heterogeneities have rarely been assessed for regional or global SGD nutrient export estimates. To evaluate heterogeneity in subterranean estuary processes and nutrient export, we collected seasonal pore water samples in 2015-2016 at three proximal (<20 km) subterranean estuary sites in Indian River Lagoon, FL. Sites have homogenous hydrogeological settings, but differ in land use and coastal features, and include a mangrove site, an urban site, and a site offshore of a natural wetland. All sites exhibit little variation through time in nutrient concentrations and modeled SGD rates. In contrast, each site exhibits significantly different nutrient concentrations of potential fresh groundwater sources, fresh groundwater discharge volumes, and nutrient transformations within subterranean estuaries. Groundwater specific discharge correlates with nutrient concentrations, suggesting that higher residence times in the subterranean estuary increase biogeochemical transformations that reduce anthropogenic nutrient loads but increase in situ nutrient sources derived from organic matter remineralization. The differences in transformations lead to SGD nutrient contributions that differ by orders of magnitude between sites and have N:P ratios that are greater than the Redfield ratio (15) for the mangrove (29) and urban sites (28), but less than the Redfield ratio for the wetland site (8). These results indicate that heterogeneity of both absolute and relative nutrient export via SGD complicates

INTRODUCTION
Groundwater is an important component of coastal water budgets and contributes terrestrial solutes to coastal zones through submarine groundwater discharge (SGD). SGD includes both fresh groundwater (herein; fresh SGD), whose discharge is driven by inland groundwater hydraulic head, and recirculated seawater (herein: saline SGD), whose discharge is driven by processes including tidal pumping and bioirrigation (Martin et al., 2006;Taniguchi et al., 2019). Both fresh and saline SGD are important for coastal biogeochemical processes; however, only fresh SGD contributes "new" terrestrial solutes to coastal zones while saline SGD contributes recycled marine solutes to coastal zones. Terrestrial solutes transported by fresh SGD are increasingly recognized to affect coastal chemical budgets (Luijendijk et al., 2020) and SGD may be their principal or sole source where surface water runoff is limited (Burnett et al., 2006;Pain et al., 2020). Both terrestrial and recycled marine nutrients transported by SGD have been implicated as drivers of benthic primary productivity (Carruthers et al., 2005) and harmful algal blooms (Phlips et al., 2002;Hu et al., 2006).
Estimates of nutrient delivery by SGD are necessary to evaluate the potential impacts to coastal ecosystems. These estimates rely on accurate assessments of both water fluxes and solute compositions. Methods to estimate water fluxes range from water mass balance approaches to the measurement of conservative geochemical tracers (typically Ra and Rn isotopes; methods extensively reviewed in Taniguchi and others, 2019). While these methods provide a means to estimate total contributions of freshwater to the coast, they can only be used to estimate total nutrient inputs when coupled with an appropriate estimate of average discharging groundwater composition. Inland groundwater composition is frequently used as a freshwater end member to estimate nutrient fluxes via SGD (e.g., Paytan et al., 2006;Rodellas et al., 2015). However, the assumptions made in estimating groundwater composition are challenged by heterogeneity in inland groundwater composition in space and time, which can be related to hydrogeological factors (Zamrsky et al., 2020), vegetation type and coverage (McGowan and Martin, 2007;Humphries et al., 2011;Alaghmand et al., 2014), and anthropogenic nutrient sources including fertilizers and wastewater.
Once fresh groundwater flows offshore, it enters a biogeochemically active subsurface freshwater-saltwater mixing zone known as a subterranean estuary (Moore, 1999). The subterranean estuary is characterized by steep salinity and redox gradients that drive a variety of reactions and nutrient transformations. The types and magnitudes of reactions depend on factors such as inflowing groundwater composition (Slomp and Van Cappellen, 2004), rock:water ratio and hydrogeological characteristics (Pain et al., 2019b), and characteristics of sediments and organic matter available to drive reactions (Moore, 1999). The net impact of these reactions may be site-specific and result in net increases or decreases in solutes of interest (Kroeger et al., 2007;Kroeger and Charette, 2008;Santos et al., 2009;Knee and Paytan, 2011;Erler et al., 2014). Estimates for SGD nutrient delivery therefore must include both heterogeneity in the composition of groundwater sources as well as the complex and heterogeneous biogeochemical reactions that take place in subterranean estuaries prior to discharge to surface waters. These factors, and their interactions, are rarely considered in SGD nutrient loading estimates.
Here, we evaluate discharge dynamics and nutrient concentrations at three proximal, but distinct, subterranean estuary sites to evaluate the extent of heterogeneity in SGD nutrient loading to the Indian River Lagoon (Florida, United States). We couple concentration information with groundwater specific discharge, estimated through advectiondiffusion models, to evaluate the variability in point fluxes. Total nutrient loading from each site is estimated by integrating nutrient point fluxes with distance offshore. We additionally evaluate the relationships between groundwater specific discharge and nutrient concentrations in order to assess hydrological controls of the observed heterogeneity. We propose that nutrient transport and upscaling of nutrient fluxes can be improved based on concentration-discharge relationships, which are frequently used to evaluate hydrologic and biologic controls on solute fluxes in rivers (e.g., Moatar et al., 2013). These assessments will provide new information regarding the extent and drivers of heterogeneity in SGD nutrient fluxes, which may allow a better understanding of how fluxes may respond to changing hydrological conditions resulting from factors such as changing groundwater recharge patterns and sea level rise.

Study Site
The Indian River Lagoon (Florida, United States; Figure 1A) is a large (approximately 5,670 km 2 ) microtidal (<10 cm) lagoon (Figure 1). In the past several decades, the lagoon water quality and clarity has substantially deteriorated due to recurring harmful algal blooms (Phlips et al., 2004;Kamerosky et al., 2015) that are increasing in both magnitude and duration. Increasing external nutrient loads from lawns and agriculture draining to the lagoon, as well as groundwater contributions, have been implicated in driving blooms (Phlips et al., 2002).
The hydrogeologic setting of the lagoon is relatively homogenous. The surficial aquifer contributes most groundwater to the lagoon and is composed of Pliocene-aged deposits of quartz sand and coquina with some thin confining clay layers (Martin et al., 2004) with hydraulic conductivity of approximately 1 × 10 −3 cm/sec (Zimmermann et al., 1985). Freshwater contributions are estimated to range from 0.01 to 0.3 × 10 6 m 3 per kilometer of shoreline per year (m 3 km −1 yr −1 ; Martin et al., 2007), the magnitude of which is largely driven by interannual variation in groundwater recharge (Roy et al., 2013). Seasonal variation in lagoon water salinity and fresh groundwater head causes fluctuations in seepage face width (Roy et al., 2013), and storm-driven saltwater intrusion events can alter seepage face salinity for several months (Smith et al., 2008a). Saline SGD contributes approximately 43 × 10 6 m 3 km −1 yr −1 , and because of the small tidal range, fetch, and wind waves, saline SGD is largely driven by the large number of bioturbating organisms in the Indian River Lagoon (Martin et al., 2006).
We evaluate SGD dynamics and potential nutrient loading from three seepage faces offshore of regions of distinct inland land use: Eau Gallie North (EGN; Figure 1B) is offshore of a highly developed commercial area, Riverwalk Park (RWP; Figure 1C) is offshore of the last remaining natural wetland bordering Indian River Lagoon, and Banana River Lagoon (BRL; Figure 1D) is offshore of moderately developed residential properties. The shoreline of BRL is lined with mangroves while the other two locations have sandy beaches at the shoreline. Permanent multilevel wells (multisamplers; Martin et al., 2003) were installed perpendicular to the shoreline across each seepage face site. Multisamplers were constructed with multiple (4-8) well screenings at depths ranging from 7 cm to variable depths up to 2.5 m below the sediment-water interface (Figure 2A). Tubing led from the screened intervals to the surface and was sampled by pumping pore water using a peristaltic pump. Multisamplers were installed in 2004 at EGN and between May 2014 and September 2015 at RWP and BRL. Multisamplers are located 0, 10, 20, and 22.5 m offshore at EGN, 10, 20, and 35 m offshore at RWP, and 1, 11, 21, and 45 m offshore at BRL.
Freshwater end member selection is critical to evaluate compositional changes caused by reactions within the subterranean estuary and determine whether biogeochemical processing is a net source of sink of terrestrial nutrients within the subterranean estuary. We define three potential freshwater Regression lines include all data available except for the most distal piezometer at EGN (17.5 m offshore) in September 2015, as both the 15-and 17.5-m offshore piezometers had zero groundwater specific discharge. RWP-20 is not included in September 2015 discharge estimates because its inclusion forces the x-intercept to extend beyond RWP-35, where zero advective flow was calculated based on porewater chloride profiles. P-values are only indicated where the correlation between variables is significant. end members to evaluate the possible range of groundwater compositions flowing to each seepage face sites. One potential end member was sampled from inland wells near each seepage face (Figures 1B-D). This water may be modified by reactions as it flows along freshwater flow paths to the offshore seepage face and the flow paths may not connect wells with our sampled seepage face. The second and third potential end members are hereby collectively referred to as offshore end members: one offshore end member is the sample collected from the deepest port of the multisampler closest to the shoreline. These samples represent groundwater that is farthest from the lagoon water and is assumed to have the least interaction with overlying saltwater and lagoon sediments. The other offshore end member is the freshest groundwater sampled at each seepage face during each sampling time. The freshest sample is not necessarily located nearest the shoreline and thus may have interacted with lagoon sediments to a greater degree than the deep shoreline end member but represents the seepage face sample with the least contribution from saline pore waters. We differentiate these three end members to assess their impacts on nutrient dynamics in subterranean estuaries, as well on estimates of nutrient fluxes.

Field Methods
We collected water samples during four sampling campaigns in September 2015, and May, August, and December 2016. Samples were collected by pumping water to the surface through 0.5 cm diameter flexible poly(vinyl chloride) tubing. Sample tubing connected the multisampler ports to an in-line overflow cup containing a YSI Pro-Plus sensor, calibrated daily, that measured salinity, temperature, pH, dissolved oxygen concentration, and oxidation-reduction potential (ORP) of the pumped water. Once these parameters were stable, water samples were filtered through 0.45-µm trace-metal grade Geotech medium capacity disposable canister filters and collected and preserved in the field using methods appropriate for each solute. Anion samples for chloride analysis were collected in 20-mL HDPE vials. No preservative was added and samples were kept chilled (4 • C) in the dark before analysis. Nutrient samples were collected in 15-mL polypropylene vials and frozen until analysis.

Laboratory Methods
Prior to analyses for nutrient concentrations, frozen samples were thawed at room temperature for 24 h and chilled samples were warmed to room temperature. Concentrations of NO 2 + NO 3 , NH 4 , and PO 4 were measured using a Seal AA3 HR Autoanalyzer. PO 4 concentrations were measured using the ammonium molybdate spectrophotometric method. NO 2 + NO 3 concentrations were measured using a cadmium coil reduction of NO 3 to NO 2 followed by sulfanilamide reaction in the presence of N-(1-naphthyl)ethylenediamine dihydrochloride. NH 4 concentrations were measured using the Berthelot reaction with alkaline phenate and hypochlorite and sodium nitroprusside. Detection limits were 0.12, 0.06, and 0.009 µM for NO 2 + NO 3 , NH 4 , and PO 4 , respectively, and values below the detection limit are reported as zero values. The sum of the NO 2 + NO 3 and NH 4 concentrations is reported as dissolved inorganic nitrogen (DIN), and the ratio of DIN to PO 4 is taken as the molar N:P ratio. Anion and cation concentrations were measured on an automated Dionex ICS-2100 and ICS-1600 Ion Chromatograph, respectively. Error on replicates was less than 5%.

Advection-Diffusion Modeling for Groundwater Specific Discharge
For conservative components such as chloride (Cl), vertical concentration gradients can be modeled as a balance between downward diffusion of the solute into the freshwater seeping to the lagoon and upward advection of fresh water into the lagoon . Where fresh SGD occurs, depth profiles of Cl concentrations asymptotically approach Cl concentrations of the fresh groundwater end member with depth in the sediment. The specific discharge decreases offshore resulting in no decrease in Cl concentrations on the lagoon side of the seepage face. The upper portion of profiles within the seepage face exhibit Cl concentrations that are altered by changes in lagoon water Cl concentrations caused by evaporation and precipitation as bioirrigation circulates lagoon water into shallow sediment (Martin et al., 2006). We used advection-diffusion modeling following methods outlined in Martin et al. (2007) to calculate the advection rate based on curvature of the porewater Cl profile and multiplied by average sediment porosity, assumed to be 0.45, to calculate specific discharge (q; cm/d). We then plotted values of q versus distance offshore and integrated under a linear regression line between distance offshore and q, which linearly decreases offshore (Glover, 1959). The regression provides an estimate for discharge from the seepage face for an infinitely thin transect perpendicular to the shoreline where Q is the total discharge (m 2 /day). These values are converted to discharge (L 3 /t) per meter of shoreline by assuming that discharge does not vary over a 1-m wide section of the shoreline centered on the infinitely thin transect. With this assumption, multiplying Q by 1 m provides an estimate of discharge from the transect in units of volume/time per m of shoreline.

Nutrient Point Fluxes and Total Nutrient Loads
We estimate nutrient point fluxes at each multisampler location and for each sampling period by multiplying the groundwater specific discharge calculated with advection-diffusion models by the average concentration of nutrients across the sediment porewater profiles for each sampling period. We estimate uncertainty by multiplying the calculated specific discharge plus or minus one standard deviation of the concentration of nutrients in the pore water profile (Supplementary Figures 1-3).
To convert individual nutrient point discharges to total seepage face nutrient loads, we integrate between multisamplers, assuming a constant rate of change between each pair of multisamplers (Supplementary Figure 4). Unlike groundwater specific discharge, offshore changes in nutrient point fluxes are not a priori expected to follow a linear or exponential decay because nutrient concentrations may vary with local sediment composition, biogeochemical reactions mechanisms and rates, and groundwater residence time. Stepwise integration between each multisampler to estimate total nutrient load, rather than with a regression as for total groundwater discharge estimates, better accounts for heterogeneous nutrient point fluxes with distance offshore. Because little seasonal variation is observed compared to variation between sites, we report average yearly fluxes for each site. Average yearly fluxes are calculated by integrating point nutrient fluxes with distance offshore for each sampling and averaging over the four sampling campaigns.

Discharge Modeling Results
Chloride modeling yields values of q that generally are highest at the shoreline multisamplers and systematically decrease to no modeled flow with distance offshore. At EGN, the maximum value of q occurs at the shoreline for three sample times and overall q is about an order of magnitude greater than at the other two sample sites (Figure 2A). No information at the shoreline was available for December 2016 because the multisampler was destroyed. However, a linear extrapolation between the two multisamplers located at 10 and 20 m offshore yielded a maximum value of q at the shoreline that was slightly lower (0.08 cm/d) than other time points. Because the true q value may have been higher than estimated considering higher values at other sample times, we consider our December 2016 discharge assessments for EGN to be conservative. At RWP, values of q generally decreased with distance offshore, but exhibited maxima in September, May, and August 20 m offshore, potentially reflecting differences in the sediment hydraulic conductivity ( Figure 2B). Average values of q for BRL were intermediate between RWP and EGN and generally decreased offshore with maxima typically observed at 11 m offshore ( Figure 2C).
The maximum values of q differ significantly between sites, with the highest time-averaged maximum q observed at EGN (0.11 cm/d) compared to BRL (0.03 cm/d) and RWP (0.01 cm/d; Figure 3A). Seepage face widths vary significantly between sites, with lowest time-averaged value at EGN (17.6 m), followed by RWP (35.7 m) and BRL (67.8 m; Figure 3B). Seepage face widths at EGN and RWP are relatively invariable through time compared to BRL, which exhibits a higher degree of temporal variability (from approximately 50 to 110 m). The time-averaged amount of water discharged over each sampling campaign also differs significantly between sites. EGN and BRL FIGURE 3 | Box plots of the distribution of seepage face hydrologic characteristics over the four sampling campaigns including (A) distribution of maximum groundwater specific discharge (q; cm/day) measured at each sampling time, (B) distribution of seepage face widths measured at each sampling time, and (C) distribution of total seepage face discharge measured at each sampling time. The average value is indicated with a number above each boxplot (n = 4 in all cases). P-values indicate significant differences among the three sites resulting from one-way ANOVA tests. discharge similar volumes of water, at 8.0 and 9.5 L m −2 d −1 , respectively, while RWP discharges significantly less (2.7 L m −2 d −1 ; Figure 3C). One way ANOVA comparison between sites of time-averaged maximum q, seepage face width, and total discharge were significantly different between the three sites at p < 0.0001 (Figure 3).

Potential Freshwater End Members and Subterranean Estuary Nutrient Content
With a few exceptions, potential freshwater ORP values and end member nutrient concentrations range widely at individual sites as well as between sites (Figure 4). In general, freshwater and saltwater end members have ORP values elevated above the reducing values (as low as to −350 mV) in the subterranean estuary ( Figure 4A) with the exception of BRL, for which the deep shoreline groundwater value is similarly reducing. Saltwater end member ORP values are consistently higher than subterranean estuary values and freshwater end member values except at EGN, where freshwater end member values range from 0 to 150 mV and are higher than the slightly reducing saltwater end member (average value of −35 ± 90 mV). EGN freshwater end members are all low for NH 4 concentrations while NO 3 and PO 4 are low at the inland well but elevated for the two offshore end members (Figures 4A-C). All freshwater end members from EGN have molar N:P ratios that are elevated above the Redfield ratio (N:P of 15:1) with a range from 19 to 38. RWP end members are less variable than EGN, though span approximately 50 µM for NH 4 (Figure 4A), compared to NO 3 and PO 4 , which are within 2-3 µM in concentration (Figures 4B,C). All freshwater end member N:P ratios from RWP are lower than the Redfield ratio (Figure 4D), with a range from 4 to 11. BRL freshwater end members span a wide range for NH 4 and PO 4 , with lowest concentrations observed in the inland well, followed by the freshest offshore groundwater sample and deep shoreline groundwater. All freshwater end member N:P ratios from BRL are higher than the Redfield ratio ranging from 40 to 160, with the highest ratio observed at the inland well.
Compared to potential freshwater and saltwater end members, nutrient concentrations in the subterranean estuary are highly variable (see Supplementary Figures 1-3 for pore water depth profiles). NH 4 concentrations at EGN and RWP are generally higher than the fresh and saltwater end members, while NH 4 concentrations at BRL span the range encompassed by the potential freshwater end members. NO 3 concentrations at EGN are all near zero, except for freshest samples in the shoreline multisampler ( Figure 4B). NO 3 concentrations are relatively low and invariable with salinity for RWP and BRL. Like NO 3 , PO 4 concentrations at EGN are lower than offshore freshwater end members except for the freshest samples, though some elevated PO 4 concentrations are observed in samples with a salinity of approximately 20. PO 4 concentrations are uniformly higher than freshwater end members at RWP and BRL ( Figure 4C). N:P ratios of some subterranean estuary samples are elevated above freshwater end members for EGN and RWP, while ratios fall within the range spanned by freshwater end members at BRL (Figure 4D).

Concentration-Discharge Relationships and Nutrient Point Fluxes
Significant relationships occur between the average of all nutrient concentrations in pore water from each multisampler (Supplementary Figures 1-3) and values of q at EGN and RWP, but no relationships are significant at BRL (Figure 5). At EGN, NH 4 is significantly negatively correlated with q, while NO 3 and PO 4 are significantly positively correlated with q. At RWP, NH 4 , and PO 4 are significantly negatively correlated with q, while NO 3 is significantly positively correlated with q.
Nutrient fluxes vary between each seepage face. Seasonal variability between nutrient point fluxes is low, with most variability occurring as a function of distance from the shoreline. In general, nutrient fluxes are highest at the shoreline multisampler and decrease offshore, though the decrease does not follow a strong linear trend. Specifically, the NH 4 fluxes are similar between EGN and RWP, ranging from 0 to 15 µmol m −2 d −1 and less than BRL with a flux of up to 90 µmol m −2 d −1 (Figure 6A). NO 3 fluxes are higher from EGN (up to 400 µmol m −2 d −1 ) compared to RWP and BRL (up to 0.5 µmol m −2 d −1 , respectively; Figure 6B). PO 4 fluxes are also higher for EGN (up to 10 µmol m −2 d −1 ) compared to RWP and BRL (up to 1.5 and 3 µmol m −2 d −1 , respectively; Figure 6C).

Total Nutrient Loads
By integrating beneath the linear extrapolations between each nutrient point flux (Supplementary Figure 4), we calculate a total seepage face load for each site. The time-averaged total nutrient loads, both within and between sites, exhibit significant variability (Figure 7). BRL contributes significantly (p < 0.001) more NH 4 (321 mmol m −1 yr −1 ) than EGN and RWP (14 and 31 mmol m −1 yr −1 , respectively; Figure 7A). EGN contributes significantly (p < 0.001) more NO 3 (540 mmol m −1 yr −1 ) than RWP and BRL (1 and 4 mmol m −1 yr −1 , respectively; Figure 7B). EGN also contributes significantly more PO 4 (18 mmol m −1 yr −1 ) compared to RWP and BRL (4 and 12 mmol m −1 yr −1 , respectively; Figure 7C). The N:P ratios of nutrient loads are similar between EGN (28) and BRL (29), while RWP is significantly lower (8). The N:P ratio of average surface water collected from all three sites at all four sampling periods is similar to the N:P ratio of RWP nutrient fluxes (Figure 7D).

DISCUSSION
The three seepage face sites sampled here display significant heterogeneity in discharge, nutrient concentrations, and nutrient loads despite the similar geological and hydrological setting between sites. We first discuss likely causes of the difference in flow between sites, followed by a discussion of differences in nutrient compositions of potential freshwater sources and biogeochemical transformations in subterranean estuaries. We then describe possible causes of the relationships between average multisampler nutrient concentrations and q and differences in the estimate nutrient point fluxes at each site, followed by a discussion of likely drivers of heterogeneity and the potential impact of heterogeneity on lagoon nutrient budgets and primary productivity. Our results reflect local heterogeneity of biogeochemical processes and nutrient cycling at proximal seepage sties, which implies these heterogeneities may be common in other subterranean estuaries. We suggest that evaluations of fundamental biogeochemical and hydrologic controls of the heterogeneity will improve our ability to make accurate estimates of regional SGD loads as well as possible   Figures 1-3).

Groundwater Specific Discharge and Total Seepage Face Discharge Estimates
We observe similar groundwater specific discharge over time at the three sampled seepage faces (Figure 2). However, significant differences occur between sites in the maximum groundwater specific discharge, seepage face width, and the total volume of water discharged over the integrated seepage face (Figure 3). Variations in these three parameters likely result from differences in the nearshore topography and land use that may impact the local water balance and groundwater discharge dynamics. For example, elevated EGN specific discharge may result from a nearshore elevation that is approximately 9 m higher than the lagoon level unlike the other sites that have gradual increase in elevation with distance inland. This topographic relief would increase the hydraulic gradient between the inland aquifer and lagoon surface water level, and contribute to the 4-10 times greater maximum values for q at EGN compared with other sites (e.g., Mulligan and Charette, 2006).
Despite high specific discharge at the shoreline, the average width of the seepage face at EGN is narrower (17.6 m) than RWP (35.7 m) and BRL (67.8 m; Figure 3B). The width of the seepage face is regulated by the specific discharge of groundwater, the total amount of water being discharged and the hydrogeological characteristics such as sediment porosity and hydraulic conductivity (Glover, 1959). Considering that hydrogeologic variables are relatively constant between sites, the narrow seepage face may result from the rapid discharge at EGN. However, there is not a perfect relationship between values of q and the seepage face width among the three sites, which suggests controls on discharge other than just rates of flow. Specifically, RWP has the slowest maximum q values and lowest total discharge, Q (2.7 L/m 2 /d) compared to EGN (8 L/m 2 /d) and BRL (9.5 L/m 2 /d; Figure 3C). But the seepage face width is intermediate between EGN and BRL (Figures 3A,B). While the cause of lower flow and total water discharge at RWP is uncertain, it could result from the natural wetland inland of the RWP seepage face, which would lower the hydraulic gradient and lead to more evapotranspiration than in the urbanized BRL or EGN seepage faces (e.g., Zheng et al., 2020). Some fresh SGD may also be captured by the wetland prior to flowing to the seepage face thereby diminishing the discharge offshore.

Potential Freshwater End Members
Identifying the composition of the fresh water source to subterranean estuaries is critical to evaluate the magnitude of nutrient delivery to the subterranean estuary, how subterranean estuary reactions may alter nutrient concentrations, and the absolute flux of nutrients from subterranean to lagoon water. Potential freshwater end members defined here vary in composition resulting from differences in inland land use and shoreline characteristics. Both EGN and BRL are located offshore of developed areas (commercial and residential, respectively) and their potential freshwater end member compositions vary widely, unlike RWP, located near an undisturbed wetland, where the three potential end member compositions are similar. We discuss below possible controls of the differences in compositions and evaluate those effects relative to processing within the subterranean estuary and magnitude of nutrient fluxes to the lagoon.
Inland wells at EGN and BRL have lower concentrations of the dominant inorganic N species (NH 4 for BRL and NO 3 for EGN) and PO 4 than the two offshore potential freshwater end members (Figure 4). These differences in concentrations indicate water sampled at the inland wells do not represent freshwater compositions entering the seepage face as the terrestrial groundwater composition is altered along the flow path from the wells to the coastline. Likely causes for compositional alteration are from contamination and/or biogeochemical reactions. However, at both EGN and BRL, nutrient concentrations also vary between the two potential freshwater end members collected at the freshwater side of the seepage face. Elevated NO 3 and NH 4 concentrations in the nearshore freshwater end members at EGN and BRL, respectively, may reflect fertilizers applied in these developed settings. Alternatively, the differences in nitrogen speciation may reflect the presence of mangroves along the shoreline at BRL, which have large impacts on biogeochemical reactions of pore water at seepage faces, including organic carbon remineralization (McGowan and Martin, 2007). Remineralization reactions consume terminal electron acceptors, including NO 3 , and available NO 3 may be reduce to NH 4 through denitrification. Higher nearshore remineralization at BRL compared to EGN is supported by the lower ORP values of BRL nearshore groundwater, including the deep shoreline groundwater freshwater end member ( Figure 4A). Regardless of the exact pathway, anthropogenic contamination and biogeochemical reactions associated with mangroves appear to elevate DIN concentrations at both EGN and BRL.
In contrast with the large variability observed for EGN and BRL, the NO 3 and PO 4 concentrations are similar and low for all three potential freshwater end members at RWP but the nearshore end members have elevated NH 4 concentrations compared with the inland well. This similarity may reflect a lack of point nutrient sources resulting from land use at RWP compared to EGN and BRL. Alternatively, the elevated NH 4 concentration may reflect reducing conditions in the nearshore wetland, similar to elevated NH 4 concentrations caused by biogeochemical reactions in the shoreline mangroves at BRL.
The variations in absolute nutrient concentrations between sites are important to nutrient delivery from the seepage face, but assimilatory demand in the subterranean estuary and lagoon will depend on the relative delivery of nutrients. Although freshwater end members at all three sites have similar PO 4 concentrations, the low N:P ratios of the freshwater end member at RWP compared with EGN and BRL and Redfield ratio reflects low DIN concentrations. These low DIN concentrations may result from limited anthropogenic contamination in the natural wetland and/or assimilation in the wetland. These results suggest that coastal wetlands are important to the regulation of both absolute and relative concentrations of inflowing fresh groundwater and reflect the importance of nearshore land use to the nutrient input to subterranean estuaries, processing within the subterranean estuaries, and ultimately the delivery of nutrients to coastal waters.
The variability among potential freshwater end members, as demonstrated by our sites (Figure 4), complicates evaluation of changes in nutrient compositions within the subterranean estuary through conservative mixing models, which are commonly used in surface water estuaries (Sanders et al., 2012;Yang et al., 2015;Pain et al., 2019aPain et al., , 2020. Although lagoon water composition is clearly defined as the saline end member of mixing models, our results reflect difficulties in defining the freshwater end member composition for mixing models because freshwater composition depends on coastal land use and modification of water composition from flow inland to the subterranean estuary. Though our inland wells are within a few hundred meters of the coast, the large deviation in their compositions from the offshore freshwater end members indicate such samples are not accurate representations of the compositions of fresh groundwater flowing to the subterranean estuary. In some cases, variability among potential freshwater end member samples precludes development of conservative mixing models. For example, the range of NH 4 concentration of potential freshwater end members at BRL is greater than the NH 4 concentrations of the pore water compositions within the subterranean estuary ( Figure 4A). The choice of an end member thus results in either net production or net consumption of NH 4 . In many cases, however, similarity of freshwater end member compositions allows construction of conservative mixing models to assess how nutrient concentrations are modified by biogeochemical reactions in subterranean estuaries and ultimately the delivery via SGD of nutrients to surface water.

Variations With Salinity
Some nutrients deviate systematically from conservative mixing regardless of complications that stem from differing compositions of the potential freshwater end members. For instance, NH 4 concentrations of EGN and RWP pore waters are higher than any potential freshwater or saltwater end member ( Figure 4A). Elevated NH 4 concentrations coincide with elevated PO 4 concentrations at a salinity of ∼20 at EGN and a salinity of ∼5 to 10 at RWP (Figures 4A,C). The correspondence of elevated N and P concentrations suggests a common source for both nutrients, possibly from organic matter remineralization. Considering differences of salinity of the elevated NH 4 and PO 4 concentrations, remineralization at EGN occurs on the saline side of the subterranean estuary (salinity ∼20) but remineralization at RWP occurs near equal mixtures of lagoon and fresh water (salinity ∼5 to 10). These inferences are supported by ORP values, which reach an average minimum value of −348 mV at EGN at a salinity of 19, while low ORP values are observed at RWP between average salinities of 5 and 20 and reach a minimum of −190 to −200 mV ( Figure 4A).
Some nutrients in the nearshore freshwater end members appear to be lost from pore waters through processing in the subterranean estuary. Elevated NO 3 and PO 4 concentrations occur for offshore potential freshwater end members and low salinity water at EGN. As salinity increases above 10, the NO 3 concentrations become negligible, suggesting denitrification reactions occur as water flows into reducing sediments of the subterranean estuary (e.g., Roy et al., 2011). Although NO 3 loss occurs through multiple microbial pathways, including biological assimilation and denitrification, denitrification is cited as the predominant consumptive pathway in most subterranean systems (Canfield et al., 2010). Dissimilatory reduction of NO 3 , which would depend on activities of the sedimentary microbial communities can produce NH 4 and could partially contribute to the elevated NH 4 concentrations observed in higher salinity EGN pore waters ( Figure 4A).
In contrast with microbial redox controls of nitrogen concentrations, controls of PO 4 concentrations more likely depend on interactions with sedimentary phases. Sediments at EGN are dominated by quartz sand grains and where salinity is low have extensive iron oxide coating . Iron oxide coatings decrease with increasing salinity as sulfate in the saline water is reduced to sulfide and sequesters iron as iron sulfide minerals. Unlike iron sulfide minerals, iron oxide phases can sorb and retain dissolved PO 4 . The switch from predominantly iron oxide to iron sulfide phases with salinity may result in the elevated PO 4 concentrations observed with elevated salinity.
Changes in nutrient concentrations within the subterranean estuary that are greater or less than conservative mixing indicate that nutrient loads to the lagoon cannot be determined based on the composition of any potential freshwater end member. For example, production of NH 4 and PO 4 in the subterranean estuary will generally result in an underestimation of nutrient delivery based on potential freshwater end member concentrations. Similarly, loss of nutrients through biogeochemical reactions and sequestration with authigenic sedimentary phases such as iron oxides suggest that estimates of SGD nutrient fluxes could be overestimated by not accounting for sinks within the subterranean estuary. In addition, changes to the absolute concentrations of individual nutrients by distinct processes (e.g., microbial processing or inorganic reactions with sedimentary phases) may alter pore water N:P ratios. For instance, EGN samples exhibit clear elevations in pore water N:P ratios in samples with salinity approaching 20 (Figure 4D), where both NH 4 and PO 4 are elevated, likely indicating decoupling of the processes controlling these two nutrient concentrations.

Concentration-Discharge Relationships
The extent to which subterranean estuary pore water may be modified by reactions depends on its residence time in the reactive freshwater-saltwater mixing zone. For example, residence time was found to regulate the proportion of CH 4 that is oxidized prior to discharge, thus impacting coastal greenhouse gas fluxes from groundwater (Schutte et al., 2016). Residence time also exerts an important control on nitrogen transformations because increased residence time results in more reducing environments (Addy et al., 2005;Kroeger and Charette, 2008;Santos et al., 2008;Smith et al., 2008b;Gonneea and Charette, 2014). Residence time in the subterranean estuary depends on the flow rates through the sediments, which is controlled by hydraulic conductivity and gradients. Flow rates also regulate the rate of delivery of terminal electron acceptors and organic matter, the primary electron donor, to the reaction zone. Therefore, flow rates should impact the extent and type of biogeochemical reactions in subterranean estuaries and thus SGD nutrient loads from any individual site (e.g., Pain et al., 2019b).
Concentration-discharge relationships between sites illustrate the effects of flow on nutrient concentrations and may indicate processes controlling compositions (Figure 5). The inverse relationships of NH 4 concentrations with q at both EGN and RWP suggest that NH 4 concentrations increase with higher residence times within the seepage face. Both sites also exhibit lower NH 4 concentrations in potential freshwater end members than an extrapolation of concentration as q approaches zero ( Figure 5A). These relationships suggest in situ sources, such as remineralization of organic matter, more strongly modify the composition of water with a higher residence time in the subterranean estuary to allow increased accumulation of remineralization products including NH 4 . In contrast, both EGN and RWP exhibit positive relationships between NO 3 concentrations and discharge, suggesting that slow flow decreases NO 3 concentrations ( Figure 5B) as would be expected if the NO 3 source is derived from fresh groundwater and is consumed through dissimilatory reduction of nitrate to ammonia (DNRA) or denitrification to N 2 as water flows through the subterranean estuary. These reactions suggest that delivery of NO 3 from fresh groundwater to lagoon water would require elevated flow rates and likely would be confined to nearshore settings rather than extending across the width of the seepage face. The elevated PO 4 concentrations at high nearshore q values at EGN (Figure 5C), similar to NO 3 , suggests the two nutrients may have a similar source, possibly from anthropogenic contamination. Similarly, their low concentrations at low q values indicate both are lost in the subterranean estuaries, although by distinct mechanisms described above. In contrast, elevated PO 4 and NH 4 concentrations at RWP at low values of q suggest that both have sources from reactions that require long residence times such as organic matter remineralization (Pain et al., 2019b).
Unlike EGN and RWP, no significant relationships exist between nutrient concentrations and discharge at BRL. The lack of significant relationships suggests that flow rates may be less important than biogeochemical reactions for nutrient concentrations at this site where mangroves may increase the extent of evapotranspiration and change redox chemistry (e.g., McGowan and Martin, 2007). The effect of mangroves is show by the NH 4 concentrations in the shoreline multisamplers, which are greater by several hundred micromoles over the trend line established by non-shoreline multisampler samples ( Figure 5A).
In contrast, the NO 3 and PO 4 concentrations versus discharge (q) are highly variable at BRL with no clear difference between shoreline and offshore multisamplers. The cause of limited relationships between nutrient concentrations and q at BRL is uncertain but could result from a number of factors including sediment heterogeneity and preferential flow paths induced by bioturbation that complicate the observed nutrient-discharge relationship. BRL also has the widest seepage face of any of the locations and thus may have more distinct reaction zones as water flows through the subterranean estuary.
Significant relationships between discharge and nutrient concentrations at EGN and RWP may facilitate the estimates of nutrient fluxes from these sites, and areas analogous to these sites, if discharge can be measured and concentration-discharge relationships established. However, complex biogeochemical and hydrological drivers of variability at the BRL site in particular highlight the challenges in upscaling nutrient fluxes from such systems with inland nutrient point sources, variations in inland land use activities, intense shoreline modification due to vegetation, and offshore sediment chemical or physical heterogeneity.

Integrated Nutrient Loads and Implications for Surface Water Processes
The heterogeneity in groundwater discharge volumes, nutrient concentrations, and concentration-discharge relationships leads to a wide range in estimates of individual nutrient loads from each of the three seepage faces examined here, as well as the ratio of these nutrients (Figure 7). The estimated nutrient loads also exhibit a range for each of our four sampling campaigns, but in most cases the inter-site variations are larger than the variations through time at any one site. BRL contributes approximately 6-10 times more NH 4 per year than EGN or RWP likely due to rapid biogeochemical cycling by mangroves. EGN delivers orders of magnitude more NO 3 and approximately 2-10 times more PO 4 than RWP or BRL possibly from anthropogenic nutrient contamination at this heavily urbanized site. Both the anthropogenic (EGN) nutrients and mangrove (BRL) nutrient sources result in N:P ratios of nutrient loads (28:1 and 29:1, respectively) that are nearly twice the Redfield ratio (15:1) and these ratios exhibit a high degree of interannual variability compared to RWP. RWP contributes low quantities of nutrients that exhibit relatively low interannual variability both in total nutrient load as well as N:P ratio (8:1), which is consistent through time (Figure 7D). This low N:P ratio appears to result from low DIN concentrations, which could indicate preferential uptake and/or loss of N as water flows through the wetlands.
The impact of SGD nutrient loads on surface water processes depends on both the lagoon water nutrient concentrations and ratios as well as biological nutrient demands. Surface water N:P ratios throughout the four sampling campaigns are relatively invariable and close to the ratios delivered by nutrient loads from RWP ( Figure 7D). Therefore, nutrient loads from RWP, and likely other natural wetlands, would not serve to greatly alter Indian River Lagoon surface water nutrient ratios. EGN and BRL, on the other hand, deliver much more N relative to P than is present in surface waters (Figure 7D), and nutrient loads from these sites would therefore serve to increase surface water N:P ratios. If the lagoon ecosystems are N-limited, delivery of excess N to P could result in enhanced primary productivity due to groundwater discharge from anthropogenically impacted sites, but not from naturally vegetated sites such as RWP. A further impact may occur due to the dominant N species delivered from sites: while they have similar N:P ratios, EGN exports most N as NO 3 while BRL exports most N as NH 4 . As there is some evidence that NH 4 is more readily assimilated by phytoplankton than NO 3 (Domingues et al., 2011), the similar N:P ratios yet different N species further complicate assessments of impacts of SGD nutrient delivery on surface water biogeochemical processes.

CONCLUSION
This study reveals significant heterogeneity in nutrient loading to the Indian River Lagoon from hydrogeologically similar subterranean estuaries. Heterogeneities occur both in the composition of fresh groundwater as well as the impact of subterranean estuary processes on nutrient concentrations and fluxes to surface waters. Some controls on the heterogeneity can be constrained, for example through concentration-discharge relationships at EGN and RWP. However, the heterogeneity at BRL likely results from multiple factors, including intense biogeochemical processing of pore waters by mangroves that line the shoreline. Inland land use appears to exert an important control on the quantity, composition, and variability of SGD and subterranean estuary processes. Natural wetlands, such as inland of RWP, attenuate groundwater fluxes and lead to lower concentrations of nutrients in potential freshwater end members and subterranean estuary samples through time. Anthropogenically impacted sites, such as EGN and BRL, are more variable both in total discharge volumes as well as the composition of potential freshwater end members. Land use differences, potential freshwater end member composition, and disparate subterranean estuary nutrient processes lead to orders of magnitude differences in NH 4 and NO 3 fluxes, and four-fold differences in PO 4 fluxes from the individual sites. This variability among geographically proximal and hydrogeologically homogeneous subterranean estuaries illustrates the complications in upscaling to lagoon-wide fluxes based on observations of discrete sections of the along shore seepage face. Such upscaling may be enhanced through surveys of land use and shoreline features linked to observations of seepage face groundwater composition. Despite these challenges, SGD is an important source of nutrients in many locations and its quantification essential to fully depict coastal nutrient budgets. Improved estimates of SGD nutrient loads and drivers, including the extent and drivers of heterogeneity, will therefore be a critical step to the inclusion of SGD in Indian River Lagoon biogeochemical budgets as well as potential extrapolation of SGD to global coastal nutrient budgets.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
JM participated in conceptualization, data collection, data interpretation, reviewing and editing the manuscript, and acquired the funding for this project. AP and CY participated in conceptualization, data collection, analysis, and interpretation. AP took the lead on writing the manuscript with contributions from JM and CY. All authors contributed to the article and approved the submitted version.