Global subterranean estuaries modify groundwater nutrient loading to the ocean

The main nutrient sources to the ocean include atmospheric deposition, rivers, and groundwater. Of these sources, groundwater-borne nutrients transported to the ocean via submarine groundwater discharge have remained the most uncertain at the global scale. We quanti ﬁ ed global nutrient loading via groundwater by compiling the largest meta-dataset of coastal groundwater nutrient concentrations available. Dissolved organic nitrogen was identi ﬁ ed as a key component of the ground-water nutrient pool and salinity and land cover were important drivers of nutrient concentrations. We provide evidence that nutrients behave non-conservatively in subterranean estuaries resulting in increases in groundwater inorganic nitrogen and phosphorus but decreases in organic nitrogen. Lastly, estimates of groundwater nutrient loading suggest submarine groundwater discharge deliver a similar amount of nutrients to the global ocean as rivers and nitrogen ﬁ xation. Our ﬁ ndings indicate that submarine groundwater discharge is an important source of nitrogen and phosphorus to the ocean that should be accounted for in nutrient budgets. Abstract Terrestrial groundwater travels through subterranean estuaries before reaching the sea. Groundwater-derived nutrients drive coastal water quality, primary production, and eutrophication. We determined how dissolved inorganic nitrogen (DIN), dissolved inorganic phosphorus (DIP), and dissolved organic nitrogen (DON) are transformed within subterranean estuaries and estimated submarine groundwater discharge (SGD) nutrient loads compiling > 10,000 groundwater samples from 216 sites worldwide. Nutrients exhibited complex, non-conservative behavior in subterranean estuaries. Fresh groundwater DIN and DIP are usually produced, and DON is consumed during transport. Median total SGD (saline and fresh) ﬂ uxes globally were 5.4, 2.6, and 0.18 Tmol yr (cid:1) 1 for DIN, DON, and DIP, respectively. Despite large natural variability, total SGD ﬂ uxes likely exceed global riverine


Scientific Significance Statement
The main nutrient sources to the ocean include atmospheric deposition, rivers, and groundwater.Of these sources, groundwater-borne nutrients transported to the ocean via submarine groundwater discharge have remained the most uncertain at the global scale.We quantified global nutrient loading via groundwater by compiling the largest meta-dataset of coastal groundwater nutrient concentrations available.Dissolved organic nitrogen was identified as a key component of the groundwater nutrient pool and salinity and land cover were important drivers of nutrient concentrations.We provide evidence that nutrients behave non-conservatively in subterranean estuaries resulting in increases in groundwater inorganic nitrogen and phosphorus but decreases in organic nitrogen.Lastly, estimates of groundwater nutrient loading suggest submarine groundwater discharge deliver a similar amount of nutrients to the global ocean as rivers and nitrogen fixation.Our findings indicate that submarine groundwater discharge is an important source of nitrogen and phosphorus to the ocean that should be accounted for in nutrient budgets.

Abstract
Terrestrial groundwater travels through subterranean estuaries before reaching the sea.Groundwater-derived nutrients drive coastal water quality, primary production, and eutrophication.We determined how dissolved inorganic nitrogen (DIN), dissolved inorganic phosphorus (DIP), and dissolved organic nitrogen (DON) are transformed within subterranean estuaries and estimated submarine groundwater discharge (SGD) nutrient loads compiling > 10,000 groundwater samples from 216 sites worldwide.Nutrients exhibited complex, nonconservative behavior in subterranean estuaries.Fresh groundwater DIN and DIP are usually produced, and DON is consumed during transport.Median total SGD (saline and fresh) fluxes globally were 5.4, 2.6, and 0.18 Tmol yr À1 for DIN, DON, and DIP, respectively.Despite large natural variability, total SGD fluxes likely exceed global riverine nutrient export.Fresh SGD is a small source of new nutrients, but saline SGD is an important source of mostly recycled nutrients.Nutrients exported via SGD via subterranean estuaries are critical to coastal biogeochemistry and a significant nutrient source to the oceans.
Globally, coastal waters receive large anthropogenic inputs of nitrogen and phosphorus, resulting in widespread water quality issues.Coastal eutrophication modifies biological communities, creates hypoxic, anoxic, or acidic conditions, and harms marine life (Basu et al. 2022).Nutrient inputs to the coastal ocean originate from rivers, atmospheric deposition, and

Wilson et al.
Subterranean estuaries modify groundwater nutrients 2 submarine groundwater discharge (SGD).The broad definition of SGD includes any flow of groundwater to the ocean, including both fresh terrestrial groundwater and seawater circulating through coastal aquifers on spatial scales greater than meters (Burnett and Dulaiova 2003;Moore 2010;Santos et al. 2021).
Seawater circulation delivers organic matter, oxygen, and other electron acceptors to sediments, accelerating biogeochemical reactions and the release of nutrients from aquifers to coastal waters.Seawater circulation can also be referred to as saline SGD or advective porewater exchange and is often traced by radon and radium isotopes (Taniguchi et al. 2019).
Resolving global nutrient budgets is critical for understanding marine biogeochemistry, productivity, and predicting future conditions.SGD is often overlooked in global ocean nutrient budgets (Gruber and Galloway 2008) and models (Shan et al. 2023).However, SGD-derived nutrient loads may exceed those from rivers and atmospheric deposition in some areas (Cho et al. 2018).Most existing estimates of SGD nutrient loads to the coastal ocean have large, unquantified uncertainties, often focus on local-scale observations, and frequently overlook transformations in subterranean estuaries.
Before discharge to the coastal ocean, SGD flows through the subterranean estuary, the subsurface transition zone between land and ocean (Moore 1999).Microbial processes either remove (e.g., denitrification), transform (e.g., nitrification), or produce (e.g., remineralization) inorganic nitrogen within subterranean estuaries (Ruiz-Gonz alez et al. 2021).Phosphorus is attenuated by sorption or released by desorption from particles like metaloxides (Charette and Sholkovitz 2002).Subterranean estuaries thus determine the speciation and concentration of SGDderived nutrients transported to the ocean.Quantifying SGD loads requires estimating flows and nutrient concentrations in fresh and saline SGD.
Here, we compiled a global meta-dataset with > 10,000 samples from 216 subterranean estuaries (Fig. 1) to characterize coastal groundwater nutrients, resolve drivers of nutrient distributions, quantify biogeochemical processing in subterranean estuaries, and finally, estimate SGD-derived nutrient loads to the ocean.We hypothesize that SGD is a major source of nitrogen and phosphorus to the global ocean and that biogeochemical transformations within subterranean estuaries modify SGD loads.Our global compilation builds on recent reviews focusing on tracer approaches to quantify SGD (Garcia-Orellana et al. 2021), SGD driving forces (e.g., Robinson et al. 2018;Taniguchi et al. 2019), and flux estimates from study cases (Santos et al. 2021).

Data compilation
Data were compiled from Web of Science searches and contacting authors, resulting in > 10,000 samples from 216 subterranean estuaries within 1-km of the coastline of 6 continents and 42 countries (Supporting Information Figs.S1, S2).Season and time of sampling, representing more than two decades of effort, was different for individual sites.While repeat surveys are not available for most sites, this represents the most comprehensive coastal groundwater dataset to date.The specific location, time of sampling, and all available data for each sample are reported in the Pangaea open data repository (https://doi.pangaea.de/10.1594/PANGAEA.9 55032).We defined samples with a salinity < 10 as low salinity and those with salinity > 10 as saline groundwater (Cho et al. 2018) because nutrient concentrations, speciation, and drivers are expected to depend on salinity.We also grouped samples into aquifer lithology types, continents, climate zones, and land use (see Supporting Information).To characterize the groundwater nutrient pool and resolve potential drivers of nutrient flux and speciation, a Random Forest analysis (Fig. S3) and linear regressions were used (Figs.S4-S10) as explained in Supporting Information.

Subterranean estuary transformations
To assess nutrient transformations, we applied the standard estuarine mixing model initially developed for surface estuaries (Boyle et al. 1974;Officer and Lynch 1981) and later applied to subterranean estuaries (Ullman et al. 2003;Santos et al. 2009a).This analysis assumes steady-state, homogeneous mixing of fresh groundwater and seawater, and the presence of only two water sources with constant endmembers.Each observed nutrient concentration was compared to theoretical conservative mixing to determine whether nutrients were produced or consumed along the salinity gradient.The result is an estimate of net transformations that provides insight to biogeochemical cycling within the subterranean estuary.The data (concentration vs. salinity) used in the mixing models developed for each subterranean estuary can be visualized in our online open-access tool (https://marineresearch.shinyapps.io/Gobal_STE_Nutrients/),including examples of sites where the model cannot be applied due to invalid assumptions.Further information regarding endmember determination and site selection can be found in the Supporting Information.Here, we focus on large scale patterns of nutrient behavior in subterranean estuaries.

Calculating SGD endmembers and nutrient fluxes to the global ocean
With our focus on subterranean estuary concentrations and transformations, we relied on previously estimated SGD water fluxes to estimate nutrient loads.The best available total (Cho and Kim 2016;Cho et al. 2018) and fresh (Kwon et al. 2014;Zhou et al. 2019;Luijendijk et al. 2020) SGD water fluxes were used to represent the potential range of nutrient loads.Total SGD fluxes were derived from global 228 Ra budgets (Cho et al. 2018) that can neither resolve the contribution of fresh and saline SGD nor the relative contribution of different aquifer types and driving forces (Taniguchi et al. 2019).These 228 Ra-derived SGD rates exclude diffusive benthic fluxes.These estimated water fluxes are thus naturally variable and have high uncertainties.Monte Carlo simulations were used to define endmember nutrient concentrations; probability density functions of various subsets of the data, including the total dataset and subsets grouped by continents, oceans, aquifer lithologies, salinities, and latitude were analyzed in R studio (R Team 2020).These subsets allow us to calculate possible ranges of nutrient endmembers and SGD loads, and offer probability ranges rather than absolute values.Each term in the equations below was treated as a stochastic variable in Monte Carlo simulations that were randomly calculated 1,000,000 times following theoretical distributions (Supporting Information Tables S4, S5) fitted to the empirical distributions.All distributions and loads are shown in the Supporting Information including subsets related to ocean basins, lithology, and water fluxes.
We first calculate nutrient loading via total SGD: where F SGD is the total SGD nutrient load, [C] gw is the groundwater endmember nutrient concentration, and Y sgd is the total SGD water flux.The estimated nutrient load results from the entire distribution of groundwater nutrient concentrations.
We also calculated net export via saline SGD: where F sSGD is the net saline SGD nutrient load, [C]s gw is the groundwater endmember nutrient concentration in samples with saliny > 10, and [C] sw is the seawater nutrient concentration (Garcia et al. 2013).This approach assumes that saline SGD (or seawater recirculation in aquifers) dominates SGD and estimates the net load after seawater infiltrates and exfiltrates sediments.Nutrient loads via fresh SGD only were calculated as: where F fSGD is the fresh SGD nutrient load, [C] fgw is the fresh groundwater endmember concentration, and Y fsgd is the fresh SGD water flux.Since fresh groundwater will flow through the subterranean estuary, we applied a correction term to the fresh endmember concentration (M) to account for transformations along the flow path.M was calculated for each site as: where M represents net transformations in the subterranean estuary (e.g., production, reduction), [C] conserv.represents the sum of conservative nutrient concentrations, and [C] observ.is the sum of observed concentrations.M quantifies the disparity between theoretical and observed nutrient concentrations across the subterranean salinity gradient.When M < 0, transformations reduce nutrient concentrations, when M > 0, concentrations increase, and when M = 0 no change occurs.The M term is applied to groundwater nutrient concentrations via: The distribution of M was applied to nutrient concentrations, [C fgw ], resulting in a transformed endmember distribution that was used to re-estimate fresh SGD fluxes.Nutrient concentrations are reported in Tmol m À3 and water fluxes are reported in m 3 yr À1 resulting in global loads in Tmol yr À1 with associated uncertainties from Monte Carlo simulations.

Global distribution of nutrients in subterranean estuaries
Coastal groundwater nutrient concentrations ranged widely with global averages (AE standard error) of 28 AE 8.7, 18 AE 9.5, and 0.9 AE 15 μM for dissolved inorganic nitrogen (DIN), dissolved inorganic phosphorus (DIP), and dissolved organic nitrogen (DON), respectively.Most observations were from Asia (22%), Europe (28%), and North America (33%), with scarce data from Central and South America, Africa, and polar coastlines.Salinity explained 34-52% of groundwater nutrient concentrations (Supporting Information Fig. S3).Mean concentrations of DIN and DON were approximately two times greater in low-salinity groundwater than in saline groundwater (Supporting Information Table S1).Wide variation in DIP concentrations was observed at all salinities.Latitude had the 2 nd largest influence (4-25%) on groundwater nutrient concentrations.DIN and DIP increased with latitude, whereas DON was lower at higher latitudes, possibly due to population density or fertilizer usage (Supporting Information Fig. S4).
Crop cover in a 1-km radius surrounding sites had an influence of 3-14% for all groundwater nutrient concentrations.Bare land was associated with increased groundwater NO À 3 and decreased NH þ 4 (Supporting Information Figs.S5, S6).In contrast, tree cover was associated with lower NO À 3 and higher NH þ 4 .Hence, anthropogenic activities, including agriculture and urbanization, may affect groundwater nutrient concentrations and speciation.Decreased groundwater DON concentrations were observed concurrently with increased built-up land cover.No significant relationships between groundwater nutrients and rainfall, evaporation, baseflow (Rodell et al. 2004), or local fresh SGD (Luijendijk et al. 2020) were observed at the global scale, though they may be important locally (Supporting Information Figs.S7, S8).

Nutrient speciation and ratios in subterranean estuaries
The speciation of N drives its fate and the biogeochemical impact of SGD on coastal ecosystems (Santos et al. 2021).Nutrient sources, redox conditions, organic matter, and microbial communities may drive N speciation.To assess N speciation patterns in coastal groundwater, subterranean estuaries were classified into four lithologies: mixed, muddy, rocky, and sandy.Low-salinity groundwater in mixed lithology systems was dominated by DON and NO À 3 (Fig. 2).Muddy sites exhibited mainly reduced forms of N, consistent with anoxic soils of marshes and mangroves minimizing nitrification (Joye and Hollibaugh 1995).Rocky karstic and volcanic aquifers often associated with rapid groundwater flow (Tovar-S anchez et al. 2014) were dominated by DON and NO x , suggesting minimal NH þ 4 production.
DON represented a significant portion (32-44%) of the groundwater N pool.DON is often overlooked in SGD and DON data are available for only $ 18% of this dataset.Different N forms may support different phytoplankton assemblages (Taylor et al. 2006;Cira et al. 2016), and primary producers may assimilate DON (Bronk et al. 2006).To further resolve SGD's role in coastal ecosystems, future studies should include DON and dissolved organic phosphorus.
Inorganic N : P ratios can play a major role in the ecology and biogeochemistry of receiving coastal waters (Downing 1997;Jickells 1998).All lithological site types had average N :P ratios above the Redfield Ratio (16 : 1) ranging from 17 AE 1.2 in rocky sites to 59 AE 1.2 in mixed sites (Fig. 2).Domination of N over P is determined by sources and redox conditions within subterranean estuaries (Slomp and van Cappellen 2004).N accumulation often results from widespread wastewater discharge, septic systems, fertilizer use, and manure leachate (Schlesinger 2009;Rahman et al. 2021).P has limited availability due to immobilization through adsorption and co-precipitation processes (Spiteri et al. 2007).Groundwater N : P ratios will rise as aquifers are contaminated with N. Since the coastal ocean is often N-limited, excess N from SGD may result in regime shifts toward P limitation, enhancing eutrophication risks (Slomp and van Cappellen 2004).

Nutrient transformations within subterranean estuaries
We assessed how subterranean estuaries modify the concentration of nutrients exported to the ocean via SGD.Patterns of nutrient behavior in subterranean estuaries were classified as conservative, consumptive, productive, or undefined (Fig. 3).Half of the sites exhibited either production (increased concentration) or consumption (reduced concentration) of nutrients supplied by fresh groundwater.
This analysis demonstrates the diverse behavior of subterranean estuaries.Many sites (35-53%) exhibited undefined behavior potentially due to spatial heterogeneity, temporal variability, and/or additional sources.Overall, this analysis revealed that subterranean estuaries on average produce DIN and DIP, but consume DON from fresh groundwater on the global scale (Supporting Information Tables S6-S8).

Nutrient loading to the global ocean
Calculating SGD-derived nutrient loads requires defining a nutrient endmember concentration and the SGD rate.Local-   et al. 2018;Sanial et al. 2021), and do not separate fresh from saline SGD fluxes.We first consider total SGD by simply multiplying the groundwater nutrient concentration by the total SGD flux as commonly done in local-scale studies (Cho et al. 2018;Santos et al. 2021).Here, Monte Carlo simulations with random sampling described the full distribution of groundwater DIN, DON, and DIP (Supporting Information Tables S4, S5).The total SGD rate was derived from a global radium mass balance model (120 AE 30 Â 10 12 m 3 yr À1 ) (Kwon et al. 2014).Total global SGD loads of DIN ranged from 2.0 to 15 Tmol yr À1 (25% and 75% quartiles, median 5.4 Tmol yr À1 ) and DIP loads ranged from 0.07 to 0.5 Tmol yr À1 (median 0.18 Tmol yr À1 ; Fig. 4).Our median total SGD DIN load is on the same order of magnitude as that previously estimated (2.3 Tmol yr À1 ) based on averages from 966 groundwater samples (Cho et al. 2018).The DIP flux is three times higher than earlier estimates (Cho et al. 2018).Despite its importance to coastal budgets and ecology, no previous global estimates of SGD-derived DON loading have been reported (Berman and Bronk 2003).DON represents $ 32% (1.0-6.6 Tmol yr À1 , median 2.6 Tmol yr À1 ) of the total SGD-N export globally.
Many SGD endmember and flux scenarios were considered (Supporting Information) to explore potential nutrient load ranges.A key scenario assumes that saline SGD dominates total SGD.In this case, nutrient fluxes via total SGD were calculated using concentrations from samples with salinity > 10.Since total SGD is largely ($ 99%) composed of saline SGD on a global scale, this approach likely produces the most reasonable estimate of global nutrient loads but is a conservative estimate where fresh SGD is relevant (Luijendijk et al. 2020).Median total SGD loads estimated with the saline groundwater nutrient distribution were 3.3, 1.3, and 0.15 Tmol yr À1 for DIN, DON, and DIP, respectively.These loads, calculated with saline groundwater endmembers, are 16-50% lower than those calculated with the full distribution of coastal groundwater samples.
The global average modified fresh SGD-associated nutrient loads imply that in general there is a net production of DIN and DIP, but net removal for DON in subterranean estuaries.Biogeochemical transformations increased fresh SGD DIN and DIP loads by 130% and 5.0%, respectively.In contrast, subterranean estuary reactions decreased DON loads by 9% (Supporting Information Tables S6-S8).A global average encompasses sites where nutrients are added as well as others where nutrients are removed; thus, the magnitude of modification on a global scale is smaller than what is observed in local-scale investigations.However, ignoring transformations within subterranean estuaries misrepresents fresh SGD loads at local and global scales.Our analysis cannot resolve individual biogeochemical processes at each site but demonstrates the net effect of transformations on nutrient loads.

Implications
The global coastal groundwater dataset compiled here was used to quantify nutrient loads and biogeochemical transformations in subterranean estuaries.Both fresh and saline groundwater nutrient concentrations are higher than those in seawater.Total SGD releases nutrients to the global ocean at rates that seem to outweigh riverine inputs (Fig. 4).In addition, SGD fluxes of DIN are of the same order of magnitude as global ocean nitrogen fixation rates (6.4-15Tmol yr À1 ) (Gruber and Sarmiento 1997;Luo et al. 2012;Tang et al. 2019).
Although loads transported via SGD and rivers are not completely analogous (e.g., rivers are a point source, SGD is diffuse along most shorelines and lithology controlled), the comparison gives perspective.In contrast to surface estuaries that filter out $ 10-20% of the global riverine DIN and DIP to the ocean (Sharples et al. 2017), subterranean estuaries enhance DIN and DIP loads from fresh groundwater.Fresh SGD accounts for only 0.6% of freshwater inputs to the ocean (Luijendijk et al. 2020).Even after accounting for transformations in the subterranean estuary, fresh SGD-derived nutrients represent < 1% of the global riverine nutrient input.Subterranean estuaries function analogously to surface estuaries, cycling terrestrially derived nutrients as they flow to the ocean via SGD.At the global scale, subterranean estuaries modify fresh groundwater nutrient concentrations and ratios (Fig. 4).
The median total (fresh and saline) SGD nutrient loads to the ocean are 3-to 4-fold greater than riverine exports.The distribution of riverine nutrient loads was simulated using reported standard errors rather than the raw observations (Seitzinger et al. 2005(Seitzinger et al. , 2010;;Letscher et al. 2013).Therefore, this distribution likely does not represent the true uncertainty of riverine nutrient loads.Our stochastic approach for calculating global SGD loads helps to resolve uncertainties and minimize biases from SGD hotspot areas.Total SGD delivers a mixture of new nutrients, derived from fresh and saline groundwater, and recycled nutrients that are flushed out from sediments and aquifers (Taniguchi et al. 2019).Although we cannot resolve their relative contribution, both new and recycled nutrients contribute to coastal primary production.With the potential for SGD to export nutrients at rates comparable to or greater than rivers, nutrient budgets that ignore SGD are incomplete.
Salinity and land-use (which also correlates to latitude) were major factors determining nutrient concentrations in SGD.Increasing human population and climate change will alter land uses, increase global temperatures, and drive groundwater salinization (Ferguson and Gleeson 2012;Taylor et al. 2013), likely enhancing SGD-derived nutrient loads to the ocean (Bowen et al. 2007;McDonough et al. 2020).SGD may impact events of coastal hypoxia, harmful algal blooms, and cause water quality degradation (Kwon et al. 2017) with potentially large economic and societal consequences.To develop innovative management practices and protect water quality, SGD must be considered a major driver of coastal biogeochemistry at local, regional, and global scales.

Fig. 4 .
Fig. 4. Contrasting the distribution of global river, total SGD, and fresh SGD fluxes.(a) Conceptual diagram of nutrient loads to the ocean based on endmembers from the Monte Carlo-derived median and 25-75% quartiles.Density plots relying on Monte Carlo analysis comparing most likely fluxes of (b) DIN, (c) DON, and (d) DIP to the global ocean via rivers (Seitzinger et al. 2005, 2010; Letscher et al. 2013), SGD (total SGD), transformed fresh SGD (tFSGD, overlaps with FSGD for DON and DIP), and fresh SGD (FSGD without accounting for transformation in the subterranean estuary).The total SGD flux was derived from a global scale radium mass balance model that quantified fresh and saline SGD (Cho et al. 2018).

Table 1 .
Behavior of groundwater nutrients during transit in subterranean estuaries shown as the percentage of sites exhibiting specific distribution profiles and the total number of sites assessed for each nutrient type.