# Global subterranean estuaries modify groundwater nutrient loading to the ocean

Stephanie J Wilson, Amy Moody, Tristan Mckenzie, M Bayani Cardenas, Elco Luijendijk, Audrey H Sawyer, Alicia Wilson, Holly A Michael, Bochao Xu, Karen L Knee, et al.

<figure><img src="image_1.png" /></figure>

► To cite this version:

Stephanie J Wilson, Amy Moody, Tristan Mckenzie, M Bayani Cardenas, Elco Luijendijk, et al. Global subterranean estuaries modify groundwater nutrient loading to the ocean. Limnology and Oceanography Letters, In press, 10.1002/lol2.10390. hal-04610022

HAL Id: hal-04610022

https://hal.science/hal-04610022

Submitted on 12 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Global subterranean estuaries modify groundwater nutrient loading to the ocean

Stephanie J Wilson, Amy Moody, Tristan Mckenzie, M Bayani Cardenas, Elco Luijendijk, Audrey H Sawyer, Alicia Wilson, Holly A Michael, Bochao Xu, Karen L Knee, et al.

► To cite this version:

Stephanie J Wilson, Amy Moody, Tristan Mckenzie, M Bayani Cardenas, Elco Luijendijk, et al. Global subterranean estuaries modify groundwater nutrient loading to the ocean. Limnology and Oceanography Letters, In press, 10.1002/lol2.10390. hal-04610022

HAL Id: hal-04610022

https://hal.science/hal-04610022v1

Submitted on 12 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
LETTER

Global subterranean estuaries modify groundwater nutrient loading to the ocean

Stephanie J. Wilson ID,¹,²* Amy Moody ID,³,⁴ Tristan McKenzie ID,⁵ M. Bayani Cardenas ID,⁶ Elco Luijendijk ID,⁷ Audrey H. Sawyer ID,⁸
Alicia Wilson ID,⁹ Holly A. Michael,¹⁰ Bochao Xu ID,¹¹ Karen L. Knee ID,¹² Hyung-Mi Cho ID,¹³ Yishai Weinstein ID,¹⁴ Adina Paytan ID,¹⁵
Nils Moosdorf ID,¹⁶ Chen-Tung Aurthur Chen ID,¹⁷ Melanie Beck ID,¹⁸ Cody Lopez,¹⁹ Dorina Murgulet ID,¹⁹ Guebuem Kim ID,²⁰
Mathew A. Charette ID,²¹ Hannelore Waska ID,¹⁸ J. Severino P. Ibánhez ID,²² Gwénaëlle Chaillou ID,²³ Till Oehler ID,²⁴
Shin-ichi Onodera ID,²⁵ Mitsuyo Saito ID,²⁵ Valenti Rodellas ID,²⁶ Natasha Dimova ID,²⁷ Daniel Montiel ID,²⁷ Henrietta Dulai ID,²⁸
Christina Richardson ID,¹⁵ Jinzhou Du ID,²⁹ Eric Petermann ID,³⁰ Xiaogang Chen ID,³¹ Kay L. Davis ID,³² Sebastien Lamontagne ID,³³
Ryo Sugimoto ID,³⁴ Guizhi Wang ID,³⁵ Hailong Li ID,³⁶ Américo I. Torres ID,³⁷ Cansu Demir ID,⁶ Emily Bristol ID,⁶ Craig T. Connolly ID,⁶
James W. McClelland ID,⁶,³⁸ Brenno J. Silva ID,³⁹ Douglas Tait ID,⁴⁰ BSK Kumar ID,⁴¹ R. Viswanadham ID,⁴¹ VVSS Sarma ID,⁴¹
Emmanoel Silva-Filho ID,⁴² Alan Shiller ID,⁴ Alanna Lecher ID,⁴³ Joseph Tamborski ID,⁴⁴ Henry Bokuniewicz ID,⁴⁵ Carlos Rocha,⁴⁶
Anja Reckhardt ID,¹⁸ Michael Ernst Böttcher ID,⁴⁴,⁴⁷,⁴⁸ Shan Jiang ID,³⁰ Thomas Stieglitz,⁴⁹ Houégnon Géraud Vinel Gbewezoun ID,⁵⁰
Céline Charbonnier,⁵¹ Pierre Anschutz ID,⁵¹ Laura M. Hernández-Terrones,⁵² Suresh Babu,⁵³ Beata Szymczycha ID,⁵⁴
Mahmood Sadat-Noori ID,⁵⁵ Felipe Niencheski ID,⁵⁶ Kimberly Null,⁵⁷ Craig Tobias ID,⁵⁸ Bongkeun Song ID,² Iris C. Anderson ID,²
Isaac R. Santos ID,⁵,⁴⁰

¹Smithsonian Environmental Research Center, Edgewater, Maryland, USA; ²Virginia Institute of Marine Science, Gloucester Point, Virginia, USA; ³Environnemental Protection Agency, Gulf of Mexico Division, Washington, DC, USA; ⁴University of Southern Mississippi, Hattiesburg, Mississippi, USA; ⁵University of Gothenburg, Gothenburg, Sweden; ⁶University of Texas at Austin, Austin, Texas, USA; ⁷University of Bergen, Bergen, Norway; ⁸Ohio State University, Columbus, Ohio, USA;

⁹University of South Carolina, Columbia, South Carolina, USA; ¹⁰University of Delaware, Newark, Delaware, USA;

¹¹Frontiers Science Center for Deep Ocean Multispheres and Earth System, Key Laboratory of Marine Chemistry Theory and Technology, Ministry of Education, Ocean University of China, Qingdao, China; ¹²American University, Washington, DC, USA; ¹³Inha University, Incheon, Korea; ¹⁴Bar-Ilan University, Ramat Gan, Israel; ¹⁵Department of Earth and Planetary Sciences, University of California at Santa Cruz, Santa Cruz, California, USA; ¹⁶Leibniz Centre for Tropical Marine Research (ZMT), Bremen, Germany; ¹⁷National Sun Yat-sen University, Kaohsiung, Taiwan; ¹⁸Institute for Chemistry and Biology of the Marine Environment (ICBM), Carl von Ossietzky Universität, Oldenburg, Germany; ¹⁹Center for Water Supply Studies, Texas A&M University—Corpus Christi, Corpus Christi, Texas, USA; ²⁰Seoul National University, Seoul, Korea; ²¹Woods Hole Oceanographic Institute, Woods Hole, Massachusetts, USA; ²²Instituto de Investigaciones Marinas, Spanish Council for Scientific Research (IIM-CSIC), Vigo, Spain; ²³University of Quebec at Rimouski, Rimouski, Quebec, Canada; ²⁴Leibniz

<footnote>*Correspondence: wilsonsj@si.edu</footnote>

<footnote>Associate editor: Dongyan Liu</footnote>

<footnote>**Author Contribution Statement:** IRS conceived the project and initiated the data compilation. SJW completed the data compilation, performed most of the data analysis, made most of the figures, and drafted the manuscript with support from IRS and AM. AM aided in figure development, editing the manuscript, and created the interactive Shiny App. TM performed the random forest analysis. MBC assisted with the Monte Carlo analysis. EL provided land use data for database. AS provided site characteristic data. All authors provided data, contributed to planning and discussions, and edited the manuscript. All authors approve the content and submission.</footnote>

<footnote>**Data Availability Statement:** The data used for this meta-analysis, including all raw data used and characteristics determined for each site, have been submitted to the Pangaea data repository (https://doi.pangaea.de/10.1594/PANGAEA.955032). The data have also been published as an open-access Shiny App, which allows for interactive visualization and exploration of data by site or variable (https://marineresearch.shinyapps.io/Gobal_STE_Nutrients/).</footnote>

<footnote>Additional Supporting Information may be found in the online version of this article.</footnote>

<footnote>This is an open access article under the terms of the [Creative Commons Attribution License](https://creativecommons.org/licenses/by/4.0/), which permits use, distribution and reproduction in any medium, provided the original work is properly cited.</footnote>
## LETTER

# Global subterranean estuaries modify groundwater nutrient loading to the ocean

Stephanie J. Wilson,¹,²* Amy Moody,³,⁴ Tristan McKenzie,⁵ M. Bayani Cardenas,⁶ Elco Luijendijk,⁷ Audrey H. Sawyer,⁸ Alicia Wilson,⁹ Holly A. Michael,¹⁰ Bochao Xu,¹¹ Karen L. Knee,¹² Hyung-Mi Cho,¹³ Yishai Weinstein,¹⁴ Adina Paytan,¹⁵ Nils Moosdorf,¹⁶ Chen-Tung Aurthur Chen,¹⁷ Melanie Beck,¹⁸ Cody Lopez,¹⁹ Dorina Murgulet,²⁰ Guebuem Kim,²¹ Mathew A. Charette,²¹ Hannelore Waska,²² Gwénaëlle Chaillou,²³ Till Oehler,²⁴ Shin-ichi Onodera,²⁵ Mitsuyo Saito,²⁵ Valenti Rodellas,²⁶ Natasha Dimova,²⁷ Daniel Montiel,²⁷ Henrietta Dulai,²⁸ Christina Richardson,²⁹ Jinxin Du,³⁰ Eric Petermann,³¹ Xiaogang Chen,³¹ Kay L. Davis,³² Sebastien Lamontagne,³³ Ryo Sugimoto,³⁴ Guizhi Wang,³⁵ Hailong Li,³⁶ Américo I. Torres,³⁷ Cansu Demir,³⁸ Emily Bristol,³⁸ Craig T. Connolly,³⁹ James W. McClelland,⁴⁰ Brenno J. Silva,³⁹ Douglas Tait,⁴⁰ BSK Kumar,⁴¹ R. Viswanadham,⁴¹ VVSS Sarma,⁴¹ Emmanoel Silva-Filho,⁴² Alan Shiller,⁴³ Joseph Tamborski,⁴⁴ Henry Bokuniewicz,⁴⁵ Carlos Rocha,⁴⁶ Anja Reckhardt,⁴⁷ Michael Ernst Böttcher,⁴⁴,⁴⁷,⁴⁸ Shan Jiang,⁵⁰ Thomas Stieglitz,⁴⁹ Houégnon Géraud Vinel Gbewezoun,⁵⁰ Céline Charbonnier,⁵¹ Pierre Anschutz,⁵¹ Laura M. Hernández-Terrones,⁵² Suresh Babu,⁵³ Beata Szymczycha,⁵⁴ Mahmood Sadat-Noori,⁵⁵ Felipe Niencheski,⁵⁶ Kimberly Null,⁵⁷ Craig Tobias,⁵⁸ Bongkeun Song,² Iris C. Anderson,² Isaac R. Santos,⁵,⁴⁰

¹Smithsonian Environmental Research Center, Edgewater, Maryland, USA; ²Virginia Institute of Marine Science, Gloucester Point, Virginia, USA; ³Environnemental Protection Agency, Gulf of Mexico Division, Washington, DC, USA; ⁴University of Southern Mississippi, Hattiesburg, Mississippi, USA; ⁵University of Gothenburg, Gothenburg, Sweden; ⁶University of Texas at Austin, Austin, Texas, USA; ⁷University of Bergen, Bergen, Norway; ⁸Ohio State University, Columbus, Ohio, USA;

⁹University of South Carolina, Columbia, South Carolina, USA; ¹⁰University of Delaware, Newark, Delaware, USA;

¹¹Frontiers Science Center for Deep Ocean Multispheres and Earth System, Key Laboratory of Marine Chemistry Theory and Technology, Ministry of Education, Ocean University of China, Qingdao, China; ¹²American University, Washington, DC, USA; ¹³Inha University, Incheon, Korea; ¹⁴Bar-Ilan University, Ramat Gan, Israel; ¹⁵Department of Earth and Planetary Sciences, University of California at Santa Cruz, Santa Cruz, California, USA; ¹⁶Leibniz Centre for Tropical Marine Research (ZMT), Bremen, Germany; ¹⁷National Sun Yat-sen University, Kaohsiung, Taiwan; ¹⁸Institute for Chemistry and Biology of the Marine Environment (ICBM), Carl von Ossietzky Universität, Oldenburg, Germany; ¹⁹Center for Water Supply Studies, Texas A&M University—Corpus Christi, Corpus Christi, Texas, USA; ²⁰Seoul National University, Seoul, Korea; ²¹Woods Hole Oceanographic Institute, Woods Hole, Massachusetts, USA; ²²Instituto de Investigaciones Marinas, Spanish Council for Scientific Research (IIM-CSIC), Vigo, Spain; ²³University of Quebec at Rimouski, Rimouski, Quebec, Canada; ²⁴Leibniz

<footnote>*Correspondence: wilsonsj@si.edu</footnote>

<footnote>Associate editor: Dongyan Liu</footnote>

<footnote>**Author Contribution Statement:** IRS conceived the project and initiated the data compilation. SJW completed the data compilation, performed most of the data analysis, made most of the figures, and drafted the manuscript with support from IRS and AM. AM aided in figure development, editing the manuscript, and created the interactive Shiny App. TM performed the random forest analysis. MBC assisted with the Monte Carlo analysis. EL provided land use data for database. AS provided site characteristic data. All authors provided data, contributed to planning and discussions, and edited the manuscript. All authors approve the content and submission.</footnote>

<footnote>**Data Availability Statement:** The data used for this meta-analysis, including all raw data used and characteristics determined for each site, have been submitted to the Pangaea data repository (https://doi.pangaea.de/10.1594/PANGAEA.955032). The data have also been published as an open-access Shiny App, which allows for interactive visualization and exploration of data by site or variable (https://marineresearch.shinyapps.io/Gobal_STE_Nutrients/).</footnote>

<footnote>Additional Supporting Information may be found in the online version of this article.</footnote>

<footnote>This is an open access article under the terms of the [Creative Commons Attribution License](https://creativecommons.org/licenses/by/4.0/), which permits use, distribution and reproduction in any medium, provided the original work is properly cited.</footnote>
Centre for Tropical Marine Research, Bremen, Germany; ²⁵Hiroshima University, Higashihiroshima, Japan; ²⁶Universitat Autònoma de Barcelona, Bellaterra, Spain; ²⁷University of Alabama, Tuscaloosa, Alabama, USA; ²⁸University of Hawaii, Honolulu, Hawaii, USA; ²⁹State Key Laboratory of Estuarine and Coastal Research, East China Normal University, Shanghai, China; ³⁰Federal Office For Radiation Protection, Berlin, Germany; ³¹Key Laboratory of Coastal Environment and Resources of Zhejiang Province, School of Engineering, Westlake University, Hangzhou, China; ³²Australian Institute of Marine Science, Perth, Western Australia, Australia; ³³The Commonwealth Scientific and Industrial Research Organisation, Canberra, Australian Capital Territory, Australia; ³⁴Fukui Prefectural University, Eiheiji, Japan; ³⁵State Key Laboratory of Marine Environmental Science, College of Ocean and Earth Sciences, Xiamen University, Xiamen, China; ³⁶Southern University of Science and Technology, Shenzen, China; ³⁷Centro para el Estudio de Sistemas Marinos CCT-CONICET CENPAT, Chubut, Argentina; ³⁸Marine Biological Laboratory, Woods Hole, MA; ³⁹Federal University of Pernambuco, Recife, Brazil; ⁴⁰Southern Cross University, Lismore, New South Wales, Australia; ⁴¹National Institute of Oceanography, Council of Scientific and Industrial Research, Visakhapatnam, India; ⁴²Universidade Federal Fluminense, Niterói, Brazil; ⁴³Lynn University, Boca Raton, Florida, USA; ⁴⁴Old Dominion University, Norfolk, Virginia, USA; ⁴⁵Stony Brook University, Stony Brook, New York, USA; ⁴⁶Trinity College Dublin, Dublin, Ireland; ⁴⁷Marine Geochemistry, University of Greifswald, Greifswald, Germany; ⁴⁸University of Rostock, Rostock, Germany; ⁴⁹Aix Marseille Université, CNRS, IRD, INRAE, CEREGE UM 34, France National Centre for Earth Science Studies, Aix en Provence, France; ⁵⁰Laboratoire d'Hydrologie Appliquée, Institut National de l'Eau, Université d'Abomey-Calavi, Abomey-Calavi, Benin; ⁵¹University of Bordeaux, CNRS, EPOC, Arcachon, France; ⁵²Universidad del Caribe, Cancún, Mexico; ⁵³National Centre for Earth Science Studies, Ministry of Earth Sciences, Thiruvananthapuram, India; ⁵⁴Institute of Oceanology Polish Academy of Sciences, Sopot, Poland; ⁵⁵University of New South Wales Sydney, Sydney, New South Wales, Australia; ⁵⁶Universidade Federal do Rio Grande, Rio Grande, Brazil; ⁵⁷Moss Landing Marine Labs, Moss Landing, California, USA; ⁵⁸University of Connecticut, Storrs, Connecticut, USA

## 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, 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) fluxes globally were 5.4, 2.6, and 0.18 Tmol yr⁻¹ 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
Centre for Tropical Marine Research, Bremen, Germany; ²⁵Hiroshima University, Higashihiroshima, Japan; ²⁶Universitat Autònoma de Barcelona, Bellaterra, Spain; ²⁷University of Alabama, Tuscaloosa, Alabama, USA; ²⁸University of Hawaii, Honolulu, Hawaii, USA; ²⁹State Key Laboratory of Estuarine and Coastal Research, East China Normal University, Shanghai, China; ³⁰Federal Office For Radiation Protection, Berlin, Germany; ³¹Key Laboratory of Coastal Environment and Resources of Zhejiang Province, School of Engineering, Westlake University, Hangzhou, China; ³²Australian Institute of Marine Science, Perth, Western Australia, Australia; ³³The Commonwealth Scientific and Industrial Research Organisation, Canberra, Australian Capital Territory, Australia; ³⁴Fukui Prefectural University, Eiheiji, Japan; ³⁵State Key Laboratory of Marine Environmental Science, College of Ocean and Earth Sciences, Xiamen University, Xiamen, China; ³⁶Southern University of Science and Technology, Shenzen, China; ³⁷Centro para el Estudio de Sistemas Marinos CCT-CONICET CENPAT, Chubut, Argentina; ³⁸Marine Biological Laboratory, Woods Hole, MA; ³⁹Federal University of Pernambuco, Recife, Brazil; ⁴⁰Southern Cross University, Lismore, New South Wales, Australia; ⁴¹National Institute of Oceanography, Council of Scientific and Industrial Research, Visakhapatnam, India; ⁴²Universidade Federal Fluminense, Niterói, Brazil; ⁴³Lynn University, Boca Raton, Florida, USA; ⁴⁴Old Dominion University, Norfolk, Virginia, USA; ⁴⁵Stony Brook University, Stony Brook, New York, USA; ⁴⁶Trinity College Dublin, Dublin, Ireland; ⁴⁷Marine Geochemistry, University of Greifswald, Greifswald, Germany; ⁴⁸University of Rostock, Rostock, Germany; ⁴⁹Aix Marseille Université, CNRS, IRD, INRAE, CEREGE UM 34, France National Centre for Earth Science Studies, Aix en Provence, France; ⁵⁰Laboratoire d'Hydrologie Appliquée, Institut National de l'Eau, Université d'Abomey-Calavi, Abomey-Calavi, Benin; ⁵¹University of Bordeaux, CNRS, EPOC, Arcachon, France; ⁵²Universidad del Caribe, Cancún, Mexico; ⁵³National Centre for Earth Science Studies, Ministry of Earth Sciences, Thiruvananthapuram, India; ⁵⁴Institute of Oceanology Polish Academy of Sciences, Sopot, Poland; ⁵⁵University of New South Wales Sydney, Sydney, New South Wales, Australia; ⁵⁶Universidade Federal do Rio Grande, Rio Grande, Brazil; ⁵⁷Moss Landing Marine Labs, Moss Landing, California, USA; ⁵⁸University of Connecticut, Storrs, Connecticut, USA

## 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, 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) fluxes globally were 5.4, 2.6, and 0.18 Tmol yr⁻¹ 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
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ález et al. 2021). Phosphorus is attenuated by sorption or released by desorption from particles like metal oxides (Charette and Sholkovitz 2002). Subterranean estuaries thus determine the speciation and concentration of SGD-derived 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).

## Methods

### 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.955032). 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 ²²⁸Ra budgets (Cho et al. 2018) that can neither resolve the contribution of fresh and saline SGD nor the relative contribution
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ález et al. 2021). Phosphorus is attenuated by sorption or released by desorption from particles like metal oxides (Charette and Sholkovitz 2002). Subterranean estuaries thus determine the speciation and concentration of SGD-derived 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).

## Methods

### 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.955032). 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 ²²⁸Ra budgets (Cho et al. 2018) that can neither resolve the contribution of fresh and saline SGD nor the relative contribution
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** (a) Global map of meta-dataset site locations (n = 216), including local-scale groundwater sample size (count) indicated by point size and lithology indicated by shape. Histograms of log-transformed groundwater (all samples) nutrient concentrations (µM) including (b) DIN (n = 5660), (c) DON (n = 1890), and (d) DIP (n = 4569). Scatter plots of groundwater (e) DIN, (f) DON, and (g) DIP concentrations (µM) vs. groundwater sample salinity (note the different Y-axes). Interactive sample map available at: https://marineresearch.shinyapps.io/Gabal_STE_Nutrients/.</figcaption></figure>

of different aquifer types and driving forces (Taniguchi et al. 2019). These ²²⁸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
<figure><img src="image_1.png" /><figcaption>**Fig. 1.** (a) Global map of meta-dataset site locations (n = 216), including local-scale groundwater sample size (count) indicated by point size and lithology indicated by shape. Histograms of log-transformed groundwater (all samples) nutrient concentrations (µM) including (b) DIN (n = 5660), (c) DON (n = 1890), and (d) DIP (n = 4569). Scatter plots of groundwater (e) DIN, (f) DON, and (g) DIP concentrations (µM) vs. groundwater sample salinity (note the different Y-axes). Interactive sample map available at: https://marineresearch.shinyapps.io/Gabal_STE_Nutrients/.</figcaption></figure>

of different aquifer types and driving forces (Taniguchi et al. 2019). These ²²⁸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:

$$F_{SGD} = [C]_{gw} \times Y_{sgd} \qquad (1a)$$

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:

$$F_{SGD} = ([C]_{sgw} - [C]_{sw}) \times Y_{sgd} \qquad (1b)$$

where $F_{SGD}$ is the net saline SGD nutrient load, $[C]_{sgw}$ is the groundwater endmember nutrient concentration in samples with salinity > 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:

$$F_{SGD} = [C]_{fgw} \times Y_{fsgd} \qquad (2)$$

where $F_{SGD}$ 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:

$$M = \frac{[C]_{\text{observ.}} - [C]_{\text{conserv.}}}{[C]_{\text{conserv.}}} \qquad (3)$$

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:

$$[C]_M = (M \times [C_{fgw}]) + [C_{fgw}] \qquad (4)$$

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⁻³ and water fluxes are reported in m³ yr⁻¹ resulting in global loads in

Tmol yr⁻¹ with associated uncertainties from Monte Carlo simulations.

## Results and discussion

### Global distribution of nutrients in subterranean estuaries

Coastal groundwater nutrient concentrations ranged widely with global averages (± standard error) of 28 ± 8.7, 18 ± 9.5, and 0.9 ± 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 2nd 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ánchez et al. 2014) were dominated by DON and $NO_x$, suggesting minimal $NH_4^+$ production.
shown in the Supporting Information including subsets related to ocean basins, lithology, and water fluxes.

We first calculate nutrient loading via total SGD:

$$F_{SGD} = [C]_{gw} \times Y_{sgd} \quad (1a)$$

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:

$$F_{SGD} = ([C]_{sgw} - [C]_{sw}) \times Y_{sgd} \quad (1b)$$

where $F_{SGD}$ is the net saline SGD nutrient load, $[C]_{sgw}$ is the groundwater endmember nutrient concentration in samples with salinity > 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:

$$F_{SGD} = [C]_{fgw} \times Y_{fsgd} \quad (2)$$

where $F_{SGD}$ 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:

$$M = \frac{[C]_{\text{observ.}} - [C]_{\text{conserv.}}}{[C]_{\text{conserv.}}} \quad (3)$$

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:

$$[C]_M = (M \times [C_{fgw}]) + [C_{fgw}] \quad (4)$$

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⁻³ and water fluxes are reported in m³ yr⁻¹ resulting in global loads in

Tmol yr⁻¹ with associated uncertainties from Monte Carlo simulations.

## Results and discussion

### Global distribution of nutrients in subterranean estuaries

Coastal groundwater nutrient concentrations ranged widely with global averages (± standard error) of 28 ± 8.7, 18 ± 9.5, and 0.9 ± 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 2nd 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ánchez et al. 2014) were dominated by DON and $NO_x$, suggesting minimal $NH_4^+$ production.
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Ternary plots of N speciation in groundwater samples with salinity, indicated by color, with salinity gradations including 0–5 (n = 1861, across all lithology types), 5–10 (n = 826), 10–20 (n = 1106), 20–30 (n = 1121), and 30–40 (n = 524). Plots are presented by lithological type for (a) mixed, (b) muddy, (c) rocky, and (d) sandy sites. DIN: DIP ratios for (e) mixed, (f) muddy, (g) rocky, and (h) sandy sites are shown in the 2nd row with the Redfield ratio (16 : 1) line shown in black.</figcaption></figure>

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 ± 1.2 in rocky sites to 59 ± 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 con-
centration of nutrients exported to the ocean via SGD. Pat-
terns of nutrient behavior in subterranean estuaries were
classified as conservative, consumptive, productive, or
undefined (Fig. 3). Half of the sites exhibited either produc-
tion (increased concentration) or consumption (reduced con-
centration) of nutrients supplied by fresh groundwater.

Of the subterranean estuaries assessed for NO₃⁻, 22% exhibited consumption and 22% production (Table 1), indicating active microbial processing or additional pollution sources (Loveless and Oldham 2010; Robinson et al. 2018). Of the sites assessed for NH₄⁺, 28% were consumptive, and 26% were productive. Consumption of NH₄⁺ results from nitrification, assimilation, or sorption to sediments, whereas production results from desorption, remineralization, or NO₃⁻ reduction (Rigaud et al. 2013; Rodellas et al. 2018). Of sites assessed for DON behavior, 29% consumed DON, likely due to remineralization of organic matter. Remineralization and
<figure><img src="image_2.png" /><figcaption>**Fig. 2.** Ternary plots of N speciation in groundwater samples with salinity, indicated by color, with salinity gradations including 0–5 (n = 1861, across all lithology types), 5–10 (n = 826), 10–20 (n = 1106), 20–30 (n = 1121), and 30–40 (n = 524). Plots are presented by lithological type for (a) mixed, (b) muddy, (c) rocky, and (d) sandy sites. DIN: DIP ratios for (e) mixed, (f) muddy, (g) rocky, and (h) sandy sites are shown in the 2nd row with the Redfield ratio (16 : 1) line shown in black.</figcaption></figure>

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 ± 1.2 in rocky
sites to 59 ± 1.2 in mixed sites (Fig. 2). Domination of N over P
is determined by sources and redox conditions within subterra-
nean estuaries (Slomp and van Cappellen 2004). N accumulation
often results from widespread wastewater discharge, septic sys-
tems, fertilizer use, and manure leachate (Schlesinger 2009;
Rahman et al. 2021). P has limited availability due to immobili-
zation through adsorption and co-precipitation processes (Spiteri
et al. 2007). Groundwater N : P ratios will rise as aquifers are con-
taminated 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 con-
centration of nutrients exported to the ocean via SGD. Pat-
terns of nutrient behavior in subterranean estuaries were
classified as conservative, consumptive, productive, or
undefined (Fig. 3). Half of the sites exhibited either produc-
tion (increased concentration) or consumption (reduced con-
centration) of nutrients supplied by fresh groundwater.

Of the subterranean estuaries assessed for NO₃⁻, 22% exhibited consumption and 22% production (Table 1), indicating active microbial processing or additional pollution sources (Loveless and Oldham 2010; Robinson et al. 2018). Of the sites assessed for NH₄⁺, 28% were consumptive, and 26% were productive. Consumption of NH₄⁺ results from nitrification, assimilation, or sorption to sediments, whereas production results from desorption, remineralization, or NO₃⁻ reduction (Rigaud et al. 2013; Rodellas et al. 2018). Of sites assessed for DON behavior, 29% consumed DON, likely due to remineralization of organic matter. Remineralization and
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Conceptual diagrams (upper row) of (a) conservative, (b) productive/pollution, (c) consumptive, and (d) undefined behavior in subterranean estuaries. The other rows show illustrations of individual sites with each classification for $NO_3^-$ (2nd row) and DIP (bottom row) in individual subterranean estuaries from the global meta-dataset. Similar graphs for all the 216 sites can be seen on https://marineresearch.shinyapps.io/Gabal_STE_Nutrients/. Examples include (e) Calanques of Marseille-Cassis (France), (f) Kasitsna Bay (USA), (g) Waquoit Bay (USA), (h) Killary (Ireland), (i) Rarotonga (Cook Islands), (j) Monterey Bay (Seabright, USA), (k) Gloucester Point (USA), and (l) Kona Coast (USA). Note the different axis scales in each figure.</figcaption></figure>

<table><caption>**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.</caption><thead><tr><td>Analyte</td><td>Conservative (%)</td><td>Productive (%)</td><td>Consumptive (%)</td><td>Undefined (%)</td><td>Sites</td></tr></thead><tbody><tr><td>NO<sub>3</sub><sup>-</sup></td><td>11.1</td><td>22.2</td><td>22.2</td><td>44.4</td><td>72</td></tr><tr><td>NH<sub>4</sub><sup>+</sup></td><td>10.8</td><td>26.2</td><td>27.7</td><td>35.4</td><td>65</td></tr><tr><td>DON</td><td>21.4</td><td>3.6</td><td>28.6</td><td>46.4</td><td>28</td></tr><tr><td>DIP</td><td>9.4</td><td>14.1</td><td>23.4</td><td>53.1</td><td>64</td></tr></tbody></table>

remobilization from particles (e.g., organic matter or Fe-oxides) may explain $PO_4^{3-}$ production (14%) while $PO_4^{3-}$ removal (23%) is likely due to sorption to manganese, iron, or calcium minerals (Charette and Sholkovitz 2002; Spiteri et al. 2008).

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-
<figure><img src="image_3.png" /><figcaption>**Fig. 3.** Conceptual diagrams (upper row) of (a) conservative, (b) productive/pollution, (c) consumptive, and (d) undefined behavior in subterranean estuaries. The other rows show illustrations of individual sites with each classification for $NO_3^-$ (2nd row) and DIP (bottom row) in individual subterranean estuaries from the global meta-dataset. Similar graphs for all the 216 sites can be seen on https://marineresearch.shinyapps.io/Gabal_STE_Nutrients/. Examples include (e) Calanques of Marseille-Cassis (France), (f) Kasitsna Bay (USA), (g) Waquoit Bay (USA), (h) Killary (Ireland), (i) Rarotonga (Cook Islands), (j) Monterey Bay (Seabright, USA), (k) Gloucester Point (USA), and (l) Kona Coast (USA). Note the different axis scales in each figure.</figcaption></figure>

<table><caption>**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.</caption><thead><tr><td>Analyte</td><td>Conservative (%)</td><td>Productive (%)</td><td>Consumptive (%)</td><td>Undefined (%)</td><td>Sites</td></tr></thead><tbody><tr><td>NO<sub>3</sub><sup>-</sup></td><td>11.1</td><td>22.2</td><td>22.2</td><td>44.4</td><td>72</td></tr><tr><td>NH<sub>4</sub><sup>+</sup></td><td>10.8</td><td>26.2</td><td>27.7</td><td>35.4</td><td>65</td></tr><tr><td>DON</td><td>21.4</td><td>3.6</td><td>28.6</td><td>46.4</td><td>28</td></tr><tr><td>DIP</td><td>9.4</td><td>14.1</td><td>23.4</td><td>53.1</td><td>64</td></tr></tbody></table>

remobilization from particles (e.g., organic matter or Fe-oxides) may explain $PO_4^{3-}$ production (14%) while $PO_4^{3-}$ removal (23%) is likely due to sorption to manganese, iron, or calcium minerals (Charette and Sholkovitz 2002; Spiteri et al. 2008).

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-
scale investigations usually take an average of nutrient con-
centrations in coastal aquifers (Dulaiova et al. 2010; Liu
et al. 2018; Sanial et al. 2021), and do not separate fresh from
saline SGD fluxes. We first consider total SGD by simply mul-
tiplying 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 ± 30 × 10¹² m³ yr⁻¹) (Kwon
et al. 2014). Total global SGD loads of DIN ranged from 2.0 to
15 Tmol yr⁻¹ (25% and 75% quartiles, median 5.4 Tmol yr⁻¹)
and DIP loads ranged from 0.07 to 0.5 Tmol yr⁻¹ (median
0.18 Tmol yr⁻¹; Fig. 4). Our median total SGD DIN load is on
the same order of magnitude as that previously estimated (2.3
Tmol yr⁻¹) based on averages from 966 groundwater samples
(Cho et al. 2018). The DIP flux is three times higher than ear-
lier 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⁻¹, median 2.6 Tmol yr⁻¹) 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⁻¹ 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.

We next consider fresh SGD with loads estimated at
0.0020-0.020 Tmol yr⁻¹ (25% and 75% quartiles) for DIN,
0.0010-0.010 Tmol yr⁻¹ for DON, and 0.000050-0.0050 Tmol
yr⁻¹ for DIP (Supporting Information Tables S6-S8). Fresh

<figure><img src="image_5.png" /><figcaption>**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).</figcaption></figure>
scale investigations usually take an average of nutrient con-
centrations in coastal aquifers (Dulaiova et al. 2010; Liu
et al. 2018; Sanial et al. 2021), and do not separate fresh from
saline SGD fluxes. We first consider total SGD by simply mul-
tiplying 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 ± 30 × 10¹² m³ yr⁻¹) (Kwon
et al. 2014). Total global SGD loads of DIN ranged from 2.0 to
15 Tmol yr⁻¹ (25% and 75% quartiles, median 5.4 Tmol yr⁻¹)
and DIP loads ranged from 0.07 to 0.5 Tmol yr⁻¹ (median
0.18 Tmol yr⁻¹; Fig. 4). Our median total SGD DIN load is on
the same order of magnitude as that previously estimated (2.3
Tmol yr⁻¹) based on averages from 966 groundwater samples
(Cho et al. 2018). The DIP flux is three times higher than ear-
lier 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⁻¹, median 2.6 Tmol yr⁻¹) 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⁻¹ 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.

We next consider fresh SGD with loads estimated at
0.0020-0.020 Tmol yr⁻¹ (25% and 75% quartiles) for DIN,
0.0010-0.010 Tmol yr⁻¹ for DON, and 0.000050-0.0050 Tmol
yr⁻¹ for DIP (Supporting Information Tables S6-S8). Fresh

<figure><img src="image_4.png" /><figcaption>**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).</figcaption></figure>
SGD passes through subterranean estuaries, mixes with saline porewater, and undergoes biogeochemical transformations prior to discharge. Assuming our observations represent global coastlines, we apply the calculated modification term derived from mixing analysis, *M*, to represent net transformations in the subterranean estuary. *M* ranged from 0.6 ± 0.1 for DON to 1.9 ± 0.50 for $NH_4^+$ (Supporting Information Table S3). The modified fresh SGD nutrient loads range from 0.00027 to 0.036 Tmol yr$^{-1}$ DIN (median 0.0082 Tmol yr$^{-1}$) and DIP loads were 0.000033–0.00048 Tmol yr$^{-1}$ (median 0.00013 Tmol yr$^{-1}$). These nutrient loads via modified fresh SGD demonstrate the many possible net outcomes (e.g., production and consumption) of biogeochemical cycling in subterranean estuaries.

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–15 Tmol 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, 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.

## References

* Basu, N. B., and others. 2022. Managing nitrogen legacies to accelerate water quality improvement. Nat. Geosci. **15**: 97–105. doi:10.1038/s41561-021-00889-9

* Berman, T., and D. A. Bronk. 2003. Dissolved organic nitrogen: A dynamic participant in aquatic ecosystems. Aquat. Microb. Ecol. **31**: 279–305. doi:10.3354/ame031279

* Bowen, J. L., K. D. Kroeger, G. Tomasky, W. J. Pabich, M. L. Cole, R. H. Carmichael, and I. Valiela. 2007. A review of land–sea coupling by groundwater discharge of nitrogen to New England estuaries: Mechanisms and effects. Appl. Geochem. **22**: 175–191. doi:10.1016/j.apgeochem.2006.09.002

* Boyle, E., R. Collier, A. T. Dengler, J. M. Edmond, A. C. Ng, and R. F. Stallard. 1974. On the chemical mass-balance in estuaries. Geochim. Cosmochim. Acta **38**: 1719–1728. doi:10.1016/0016-7037(74)90188-4

* Bronk, D. A., J. H. See, P. Bradley, and L. Killberg. 2006. DON as a source of bioavailable nitrogen for phytoplankton. Biogeosci. Discuss. **3**: 1247–1277. doi:10.5194/bgd-3-1247-2006
SGD passes through subterranean estuaries, mixes with saline porewater, and undergoes biogeochemical transformations prior to discharge. Assuming our observations represent global coastlines, we apply the calculated modification term derived from mixing analysis, *M*, to represent net transformations in the subterranean estuary. *M* ranged from 0.6 ± 0.1 for DON to 1.9 ± 0.50 for $NH_4^+$ (Supporting Information Table S3). The modified fresh SGD nutrient loads range from 0.00027 to 0.036 Tmol yr$^{-1}$ DIN (median 0.0082 Tmol yr$^{-1}$) and DIP loads were 0.000033–0.00048 Tmol yr$^{-1}$ (median 0.00013 Tmol yr$^{-1}$). These nutrient loads via modified fresh SGD demonstrate the many possible net outcomes (e.g., production and consumption) of biogeochemical cycling in subterranean estuaries.

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–15 Tmol 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, 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.

## References

* Basu, N. B., and others. 2022. Managing nitrogen legacies to accelerate water quality improvement. Nat. Geosci. **15**: 97–105. doi:10.1038/s41561-021-00889-9

* Berman, T., and D. A. Bronk. 2003. Dissolved organic nitrogen: A dynamic participant in aquatic ecosystems. Aquat. Microb. Ecol. **31**: 279–305. doi:10.3354/ame031279

* Bowen, J. L., K. D. Kroeger, G. Tomasky, W. J. Pabich, M. L. Cole, R. H. Carmichael, and I. Valiela. 2007. A review of land–sea coupling by groundwater discharge of nitrogen to New England estuaries: Mechanisms and effects. Appl. Geochem. **22**: 175–191. doi:10.1016/j.apgeochem.2006.09.002

* Boyle, E., R. Collier, A. T. Dengler, J. M. Edmond, A. C. Ng, and R. F. Stallard. 1974. On the chemical mass-balance in estuaries. Geochim. Cosmochim. Acta **38**: 1719–1728. doi:10.1016/0016-7037(74)90188-4

* Bronk, D. A., J. H. See, P. Bradley, and L. Killberg. 2006. DON as a source of bioavailable nitrogen for phytoplankton. Biogeosci. Discuss. **3**: 1247–1277. doi:10.5194/bgd-3-1247-2006
ocean. Adv. Water Resour. **115**: 315-331. doi:10.1016/j.adwatres.2017.10.041

Rodell, M., and others. 2004. The global land data assimilation system. Bull. Am. Meteorol. Soc. **85**: 381-394. doi:10.1175/BAMS-85-3-381

Rodellas, V., T. C. Stieglitz, A. Andrisoa, P. G. Cook, P. Raimbault, J. J. Tamborski, P. van Beek, and O. Radakovitch. 2018. Groundwater-driven nutrient inputs to coastal lagoons: The relevance of lagoon water recirculation as a conveyor of dissolved nutrients. Sci. Total Environ. **642**: 764-780. doi:10.1016/j.scitotenv.2018.06.095

Ruiz-González, C., V. Rodellas, and J. Garcia-Orellana. 2021. The microbial dimension of submarine groundwater discharge: Current challenges and future directions. FEMS Microbiol. Rev. **45**: 1-25. doi:10.1093/femsre/fuab010

Sanial, V., W. S. Moore, and A. M. Shiller. 2021. Does a bottom-up mechanism promote hypoxia in the Mississippi Bight? Mar. Chem. **235**: 104007. doi:10.1016/j.marchem.2021.104007

Santos, I. R., W. C. Burnett, T. Dittmar, I. G. N. A. Suryaputra, and J. Chanton. 2009. Tidal pumping drives nutrient and dissolved organic matter dynamics in a Gulf of Mexico subterranean estuary. Geochim. Cosmochim. Acta **73**: 1325-1339. doi:10.1016/j.gca.2008.11.029

Santos, I. R., and others. 2021. Submarine groundwater discharge impacts on coastal nutrient biogeochemistry. Nat. Rev. Earth Environ. **6**: 108-119. doi:10.1016/B978-0-12-409548-9.11482-4

Schlesinger, W. H. 2009. On the fate of anthropogenic nitrogen. Proc. Natl. Acad. Sci. USA **106**: 203-208. doi:10.1073/pnas.0810193105

Seitzinger, S. P., J. A. Harrison, E. Dumont, A. H. W. Beusen, and A. F. Bouwman. 2005. Sources and delivery of carbon, nitrogen, and phosphorus to the coastal zone: An overview of Global Nutrient Export from Watersheds (NEWS) models and their application. Global Biogeochem. Cycl. **19**: 1-11. doi:10.1029/2005GB002606

Seitzinger, S. P., and others. 2010. Global river nutrient export: A scenario analysis of past and future trends. Global Biogeochem. Cycl. **24**. doi:10.1029/2009GB003587

Shan, X., A. Zhu, D. Fu, Y. Song, Q. Li, Z. Huang, L. Pei, and H. Zhao. 2023. Modeling nutrient flows from land to rivers and seas—A review and synthesis. Mar. Environ. Res. **186**: 1-13. doi:10.1016/j.marenvres.2023.105928

Sharples, J., J. J. Middelburg, K. Fennel, and T. D. Jickells. 2017. What proportion of riverine nutrients reaches the open ocean? Global Biogeochem. Cycl. **31**: 39-58. doi:10.1111/1462-2920.13280

Slomp, C. P., and P. van Cappellen. 2004. Nutrient inputs to the coastal ocean through submarine groundwater discharge: Controls and potential impact. J. Hydrol. **295**: 64-86. doi:10.1016/j.jhydrol.2004.02.018

Spiteri, C., C. P. Slomp, P. Regnier, C. Meile, and P. van Cappellen. 2007. Modelling the geochemical fate and

transport of wastewater-derived phosphorus in contrasting groundwater systems. J. Contam. Hydrol. **92**: 87-108. doi: 10.1016/j.jconhyd.2007.01.002

Spiteri, C., P. van Cappellen, and P. Regnier. 2008. Surface complexation effects on phosphate adsorption to ferric iron oxyhydroxides along pH and salinity gradients in estuaries and coastal aquifers. Geochim. Cosmochim. Acta **72**: 3431-3445. doi:10.1016/j.gca.2008.05.003

Tang, W., Z. Li, and N. Cassar. 2019. Machine learning estimates of global marine nitrogen fixation. Eur. J. Vasc. Endovasc. Surg. **124**: 717-730. doi:10.1029/2018JG004828

Taniguchi, M., H. Dulai, K. M. Burnett, I. R. Santos, R. Sugimoto, T. Stieglitz, G. Kim, N. Moosdorf, and W. C. Burnett. 2019. Submarine groundwater discharge: Updates on its measurement techniques, geophysical drivers, magnitudes, and effects. Front. Environ. Sci. **7**: doi:10.3389/fenvs.2019.00141

Taylor, G. T., C. J. Gobler, and S. A. Sañudo-Wilhelmy. 2006. Speciation and concentrations of dissolved nitrogen as determinants of brown tide *Aureococcus anophagefferens* bloom initiation. Mar. Ecol. Prog. Ser. **312**: 67-83. doi:10.3354/meps312067

Taylor, R. G., and others. 2013. Ground water and climate change. Nat. Clim. Change **3**: 322-329. doi:10.1038/nclimate1744

Tovar-Sánchez, A., and others. 2014. Contribution of groundwater discharge to the coastal dissolved nutrients and trace metal concentrations in Majorca Island: Karstic vs detrital systems. Environ. Sci. Technol. **48**: 11819-11827. doi:10.1021/es502958t

Ullman, W. J., B. Chang, D. C. Miller, and J. A. Madsen. 2003. Groundwater mixing, nutrient diagenesis, and discharges across a sandy beachface Cape Henlopen, Delaware (USA). Estuar. Coast. Shelf Sci. **57**: 539-552. doi:10.1016/S0272-7714(02)00398-0

Zhou, Y. Q., A. H. Sawyer, C. H. David, and J. S. Famiglietti. 2019. Fresh submarine groundwater discharge to the near-global coast. Geophys. Res. Lett. **46**: 5855-5863. doi:10.1029/2019GL082749

## Acknowledgments

During the writing of this paper, our colleague Jordi Garcia-Orellana passed away too early. Jordi was a well-respected member of the SGD community and will be missed by the many co-authors of this paper who worked alongside Jordi for many years. We thank the numerous students and scientists that contributed to the work reported in this manuscript, which is the result of decades of SGD research around the world. The data compilation and some of the analysis was initiated with funding from the Australian Research Council (FT170100327) and completed with funding from Sweden's Knut and Alice Wallenberg Foundation (2022.0096). SJW and AM were funded by a US National Science Foundation and Association for Limnology and Oceanography LOREX fellowship to work in Sweden (IRES 1831075). Field and laboratory investigations to obtain the multiple datasets were funded by the US National Science Foundation (1658135 and 1737258), the Australian Research Council (DE140101733, DP120101645), the Swedish Research Council (2020-00457), the Fundacao de Amparo a
ocean. Adv. Water Resour. **115**: 315-331. doi:10.1016/j.adwatres.2017.10.041

Rodell, M., and others. 2004. The global land data assimilation system. Bull. Am. Meteorol. Soc. **85**: 381-394. doi:10.1175/BAMS-85-3-381

Rodellas, V., T. C. Stieglitz, A. Andrisoa, P. G. Cook, P. Raimbault, J. J. Tamborski, P. van Beek, and O. Radakovitch. 2018. Groundwater-driven nutrient inputs to coastal lagoons: The relevance of lagoon water recirculation as a conveyor of dissolved nutrients. Sci. Total Environ. **642**: 764-780. doi:10.1016/j.scitotenv.2018.06.095

Ruiz-González, C., V. Rodellas, and J. Garcia-Orellana. 2021. The microbial dimension of submarine groundwater discharge: Current challenges and future directions. FEMS Microbiol. Rev. **45**: 1-25. doi:10.1093/femsre/fuab010

Sanial, V., W. S. Moore, and A. M. Shiller. 2021. Does a bottom-up mechanism promote hypoxia in the Mississippi Bight? Mar. Chem. **235**: 104007. doi:10.1016/j.marchem.2021.104007

Santos, I. R., W. C. Burnett, T. Dittmar, I. G. N. A. Suryaputra, and J. Chanton. 2009. Tidal pumping drives nutrient and dissolved organic matter dynamics in a Gulf of Mexico subterranean estuary. Geochim. Cosmochim. Acta **73**: 1325-1339. doi:10.1016/j.gca.2008.11.029

Santos, I. R., and others. 2021. Submarine groundwater discharge impacts on coastal nutrient biogeochemistry. Nat. Rev. Earth Environ. **6**: 108-119. doi:10.1016/B978-0-12-409548-9.11482-4

Schlesinger, W. H. 2009. On the fate of anthropogenic nitrogen. Proc. Natl. Acad. Sci. USA **106**: 203-208. doi:10.1073/pnas.0810193105

Seitzinger, S. P., J. A. Harrison, E. Dumont, A. H. W. Beusen, and A. F. Bouwman. 2005. Sources and delivery of carbon, nitrogen, and phosphorus to the coastal zone: An overview of Global Nutrient Export from Watersheds (NEWS) models and their application. Global Biogeochem. Cycl. **19**: 1-11. doi:10.1029/2005GB002606

Seitzinger, S. P., and others. 2010. Global river nutrient export: A scenario analysis of past and future trends. Global Biogeochem. Cycl. **24**. doi:10.1029/2009GB003587

Shan, X., A. Zhu, D. Fu, Y. Song, Q. Li, Z. Huang, L. Pei, and H. Zhao. 2023. Modeling nutrient flows from land to rivers and seas—A review and synthesis. Mar. Environ. Res. **186**: 1-13. doi:10.1016/j.marenvres.2023.105928

Sharples, J., J. J. Middelburg, K. Fennel, and T. D. Jickells. 2017. What proportion of riverine nutrients reaches the open ocean? Global Biogeochem. Cycl. **31**: 39-58. doi:10.1111/1462-2920.13280

Slomp, C. P., and P. van Cappellen. 2004. Nutrient inputs to the coastal ocean through submarine groundwater discharge: Controls and potential impact. J. Hydrol. **295**: 64-86. doi:10.1016/j.jhydrol.2004.02.018

Spiteri, C., C. P. Slomp, P. Regnier, C. Meile, and P. van Cappellen. 2007. Modelling the geochemical fate and

transport of wastewater-derived phosphorus in contrasting groundwater systems. J. Contam. Hydrol. **92**: 87-108. doi:10.1016/j.jconhyd.2007.01.002

Spiteri, C., P. van Cappellen, and P. Regnier. 2008. Surface complexation effects on phosphate adsorption to ferric iron oxyhydroxides along pH and salinity gradients in estuaries and coastal aquifers. Geochim. Cosmochim. Acta **72**: 3431-3445. doi:10.1016/j.gca.2008.05.003

Tang, W., Z. Li, and N. Cassar. 2019. Machine learning estimates of global marine nitrogen fixation. Eur. J. Vasc. Endovasc. Surg. **124**: 717-730. doi:10.1029/2018JG004828

Taniguchi, M., H. Dulai, K. M. Burnett, I. R. Santos, R. Sugimoto, T. Stieglitz, G. Kim, N. Moosdorf, and W. C. Burnett. 2019. Submarine groundwater discharge: Updates on its measurement techniques, geophysical drivers, magnitudes, and effects. Front. Environ. Sci. **7**: doi:10.3389/fenvs.2019.00141

Taylor, G. T., C. J. Gobler, and S. A. Sañudo-Wilhelmy. 2006. Speciation and concentrations of dissolved nitrogen as determinants of brown tide *Aureococcus anophagefferens* bloom initiation. Mar. Ecol. Prog. Ser. **312**: 67-83. doi:10.3354/meps312067

Taylor, R. G., and others. 2013. Ground water and climate change. Nat. Clim. Change **3**: 322-329. doi:10.1038/nclimate1744

Tovar-Sánchez, A., and others. 2014. Contribution of groundwater discharge to the coastal dissolved nutrients and trace metal concentrations in Majorca Island: Karstic vs detrital systems. Environ. Sci. Technol. **48**: 11819-11827. doi:10.1021/es502958t

Ullman, W. J., B. Chang, D. C. Miller, and J. A. Madsen. 2003. Groundwater mixing, nutrient diagenesis, and discharges across a sandy beachface Cape Henlopen, Delaware (USA). Estuar. Coast. Shelf Sci. **57**: 539-552. doi:10.1016/S0272-7714(02)00398-0

Zhou, Y. Q., A. H. Sawyer, C. H. David, and J. S. Famiglietti. 2019. Fresh submarine groundwater discharge to the near-global coast. Geophys. Res. Lett. **46**: 5855-5863. doi:10.1029/2019GL082749

## Acknowledgments

During the writing of this paper, our colleague Jordi Garcia-Orellana passed away too early. Jordi was a well-respected member of the SGD community and will be missed by the many co-authors of this paper who worked alongside Jordi for many years. We thank the numerous students and scientists that contributed to the work reported in this manuscript, which is the result of decades of SGD research around the world. The data compilation and some of the analysis was initiated with funding from the Australian Research Council (FT170100327) and completed with funding from Sweden's Knut and Alice Wallenberg Foundation (2022.0096). SJW and AM were funded by a US National Science Foundation and Association for Limnology and Oceanography LOREX fellowship to work in Sweden (IRES 1831075). Field and laboratory investigations to obtain the multiple datasets were funded by the US National Science Foundation (1658135 and 1737258), the Australian Research Council (DE140101733, DP120101645), the Swedish Research Council (2020-00457), the Fundacao de Amparo a
Pesquisa do Rio de Janeiro and an Institutional Grant to the Texas General Land Office, Coastal Management Program pursuant to National Oceanic and Atmospheric Administration (No. NA15NOS419012), and the Texas A&M University—Corpus Christi, Center for Water Supply Studies. SJW thanks the Virginia Institute of Marine Science at the College of William & Mary and the University of Gothenburg for support.

Submitted 09 September 2023

Revised 16 March 2024

Accepted 23 March 2024
Pesquisa do Rio de Janeiro and an Institutional Grant to the Texas General Land Office, Coastal Management Program pursuant to National Oceanic and Atmospheric Administration (No. NA15NOS419012), and the Texas A&M University—Corpus Christi, Center for Water Supply Studies. SJW thanks the Virginia Institute of Marine Science at the College of William & Mary and the University of Gothenburg for support.

Submitted 09 September 2023

Revised 16 March 2024

Accepted 23 March 2024