A State-Of-The-Art Perspective on the Characterization of Subterranean Estuaries at the Regional Scale

Subterranean estuaries the, subsurface mixing zones of terrestrial groundwater and seawater, substantially influence solute fluxes to the oceans. Solutes brought by groundwater from land and solutes brought from the sea can undergo biogeochemical reactions. These are often mediated by microbes and controlled by reactions with coastal sediments, and determine the composition of fluids discharging from STEs (i.e., submarine groundwater discharge), which may have consequences showing in coastal ecosystems. While at the local scale (meters), processes have been intensively studied, the impact of subterranean estuary processes on solute fluxes to the coastal ocean remains poorly constrained at the regional scale (kilometers). In the present communication, we review the processes that occur in STEs, focusing mainly on fluid flow and biogeochemical transformations of nitrogen, phosphorus, carbon, sulfur and trace metals. We highlight the spatio-temporal dynamics and measurable manifestations of those processes. The objective of this contribution is to provide a perspective on how tracer studies, geophysical methods, remote sensing and hydrogeological modeling could exploit such manifestations to estimate the regional-scale impact of processes in STEs on solute fluxes to the coastal ocean.


INTRODUCTION
Along global coastlines, meteoric groundwater discharges into the ocean (Church, 1996;Taniguchi et al., 2002;Zhou et al., 2019;Luijendijk et al., 2020). The subsurface zone where meteoric groundwater mixes with saltwater was termed the subterranean estuary (STE) by Moore (1999). This term was defined as "a coastal aquifer where groundwater derived from land drainage measurably dilutes saltwater that has invaded the aquifer through a free connection to the sea" (Moore, 1999). The relevance of STEs for matter cycling and coastal ecology is increasingly recognized by the scientific community (Rocha et al., 2021).
The term estuary originates from hydrology and is defined as the area where fresh river water and saline seawater mix (Prandle, 2009). Riverine estuaries are biogeochemical reactors that, amongst other processes, filter out about 20% of the dissolved silicon through reverse weathering (Tréguer and De La Rocha, 2013) or degas substantial amounts of riverine transported organic carbon as CO 2 (Laruelle et al., 2010). By now, it is accepted that riverine solute fluxes to the coastal realm cannot be meaningfully quantified without considering transformations in the mixing zone of river estuaries (e.g., Kipp et al., 2020). For submarine groundwater discharge (SGD), the same should hold for subterranean estuaries .
STEs are characterized by waters with geochemical signatures controlled by steep (physico-)chemical gradients of mixing terrestrial groundwater and seawater (Moore, 1999). The groundwater-seawater mixing zone embraces thematically different scientific fields, sometimes leading to confusion in terminologies. The term STE, for instance, is partly overlapping with the term 'coastal aquifers', which describes a groundwater system at the interface of land and sea in hydrogeology (Duque et al., 2020;Jiao and Post, 2019). Mixing of fresh and saline groundwater occurs differently in porous sedimentary aquifers than in karst or volcanic aquifers since the former have longer residence times and a higher reaction area between solid and fluid phase. STEs of the latter type are termed 'anchialine' (Bishop et al., 2015). Other studies use the term STE for karstic environments (e.g., . In these situations, hydrology plays a major role in biogeochemical reactions (Brankovits et al., 2018), which is amplified by the karst's hydrogeological conditions. However, if karstic aquifers discharge in submarine springs with sufficient discharge rates, this can lead to freshwater mixing with seawater in the ocean itself in the form of plume-like structures (e.g., Fleury et al., 2007), whose nutrient fluxes can trigger algal blooms . Karst areas cover large parts of the global coastline (Goldscheider et al., 2020) and are essential groundwater-ocean interaction areas. Nevertheless, since processes in these conduit-systems differ sharply from those in porous sediments, we focus our assessment on STEs composed of porous sediments ( Figure 1). Thus, we here use the term STE to describe the zone where fresh and saline groundwater mix in porous coastal aquifers in the subsurface, which is in line with the marine scientific literature Duque et al., 2020).
While STEs share the name estuary and the freshwatersaltwater mixing zone, they differ sharply from their riverine counterparts. The term "estuary" itself is questionable for STEs, since some definitions of estuaries also encompass the shape of the water body, e.g., as "semi-enclosed", basically referring to rivers (Wolanski, 2007). One crucial difference is residence time, which is in the order of days in riverine estuaries (Rasmussen and Josefson, 2002) and can be decades in STEs (e.g., Grünenbaum et al., 2020). The different residence time leads to different mixing processes: while in rivers, mixing can be in turbulent flow driven by wind, flow in STEs is generally linear and advection-dispersion driven. While in river estuaries, water-solid interaction is mostly limited to suspended or surface sediment as solid phase, STEs provide a wide range of minerals and solids to interact due to a lower water/rock ratio. In some STEs, anoxic conditions develop along subsurface flow-paths or be promoted by the inflow of O 2free groundwater, which changes the biogeochemical reactions compared to STEs where oxic conditions prevail (Slomp and Van Cappellen, 2004) or to oxic surface waters flowing through riverine estuaries. Thus, despite river estuaries and STEs sharing parts of their name, the processes involved can be very different. Also, while rivers can be treated as point sources of solutes and particulates to the coastal ocean, STEs can occur along long stretches of coast, and their geochemistry can vary substantially at the meter scale (e.g., Beck et al., 2016;Ehlert et al., 2016;Beck et al., 2017;Waska et al., 2019b) and even more at the regional scale.
In general, biogeochemical reactions alter the composition of terrestrial groundwater and saltwater that flow through STEs. Thus, the biogeochemical composition of the resulting SGD will not represent a conservative mixing between the fresh and saline waters entering the STE. In many STEs, the contribution of marine water circulating through the sediment ("marine SGD") to the total water fluxes exceeds that of terrestrial groundwater (e.g., Lopez et al., 2020). Groundwater of marine origin in the seafloor sediment and STEs will here be called "saline groundwater." Flow and transport in STEs are highly dynamic and locally variable. While discharge from unconfined shallow aquifers occurs near the shoreline, discharge from confined aquifers has been observed several kilometers offshore, depending on the hydrogeological conditions Gustafson et al., 2019). Mixing of fluids triggers biogeochemical reactions that can change the discharge composition compared to the original end-members (i.e., fresh groundwater and saltwater). Understanding the coupling between the physical and biogeochemical processes is the prerequisite to understand the role of STEs in controlling the geochemical composition of SGD . Processes in STEs can substantially influence the water quality of the receiving water bodies and the biological processes in coastal waters (e.g., productivity, food web structure) and the full spectrum of the species inhabiting them, from bacteria (Adyasari et al., 2019a), algae, through grazers, detritivores, to invertebrates (Miller and Ullman, 2004), predatory fish (Pisternick et al., 2020), and birds (Kotwicki et al., 2014;Lecher and Mackey, 2018). The processes controlling the solute fluxes through STEs have been named among the main currently unsolved problems in hydrology (Blöschl et al., 2019).
Groundwater salinity distribution in STEs is typically characterized by a wedge of dense, saline groundwater that thins from the ocean inland underneath a fresh water body at the bottom of an aquifer ( Figure 2). Along these water bodies' interface, the fresh and saline groundwater mix and a zone of intermediate salinities forms. Another saline recirculation cell, where saltwater infiltrates into the sediment, the so-called 'upper saline plume' (USP), can occur in settings with a sloping land surface where tides and waves run up along the beach (e.g., Lebbe, 1981;Robinson et al., 2007). Between the USP and the saltwater wedge, a 'freshwater discharge tube' evolves, focusing fresh SGD along the low water line under an active SGD seepage zone. These distinct water bodies have been described to shift, expand and contract as a consequence of seasonal groundwater recharge variations (Michael et al., 2005), spring-neap tidal cycles (Abarca et al., 2013;Heiss and Michael, 2014) or intensified wave conditions . The shape of the STE is controlled by a large number of factors, such as geological heterogeneity , terrestrial groundwater discharge (Heiss and Michael, 2014), tidal amplitude (Abarca et al., 2013), sediment topography (Robinson et al., 2006;Waska et al., 2019b;Grünenbaum et al., 2020), or even seawater temperature (Kim et al., 2020).
STEs have been investigated from local, regional (Beck and Brumsack, 2012;Jurasinski et al., 2018) to global perspective Rahman et al., 2019), and groundwater dynamics and driving forces have been reviewed Robinson et al., 2018). For example, in Waquoit Bay, a comprehensive suite of STE studies related to biogeochemistry (Charette et al., 2005;Charette and Sholkovitz, 2006;Saenz et al., 2012;, hydrology (Michael et al., 2005;Spiteri et al., 2008b), and tracer application  has been conducted. The coastline of the Gulf of Mexico has provided a diverse geological background for local STE studies: from sandy beaches and lagoon systems in Florida (Santos et al., 2008;Roy et al., 2011;Pain et al., 2019), via an organic-rich STE in the eastern part of Mobile Bay, Alabama (Montiel et al., 2019), to developed karst system in Yucatan, Mexico Brankovits et al., 2017;Brankovits et al., 2018). Biogeochemical STE studies in North America have also been conducted in the eastern (Hays and Ullman, 2007;O'Connor et al., 2018;Tamborski et al., 2017) and western (Santoro et al., 2008;Boehm et al., 2014;Brown and Boehm, 2016) coastal areas of the contiguous United States, as well as boreal parts of Alaska (Lecher et al., 2016a) and Canada (Couturier et al., 2016;Sirois et al., 2018).
Many relevant field studies have been conducted outside North America. The development of an upper saline plume and the freshwater tube underneath was first demonstrated in Belgium by Lebbe (1981). He developed the conceptual model for the salinity distribution in unconfined aquifers influenced by tides based on field measurements and mathematical modeling. The Wadden Sea, located in the southern North Sea, is a wellstudied system where various geochemical Linkhorst et al., 2017;Reckhardt et al., 2017;Rullkötter, 2009;Seidel et al., 2015;Waska et al., 2019b), hydrological (Moore et al., 2011;Seibert et al., 2019), as well as microbiological and biogeochemical (Musat et al., 2006;Al-Raei et al., 2009) studies have been conducted. Here, STEs have been shown to function as a coastal filter for material exchange and transformation. They are characterized by a high input of organic matter, intensive nutrient recycling, and organic matter remineralization rates (Billerbeck et al., 2006;Al-Raei et al., 2009). Other studied STEs in Europe are the Baltic Sea coast (Szymczycha et al., 2012;Donis et al., 2017;Jurasinski et al., 2018;Virtasalo et al., 2019) and the coast of France (Anschutz et al., 2009;Charbonnier et al., 2013;Oehler et al., 2017). In Asia-Pacific, STE studies were conducted in South Korea (Kim et al., 2012;Lee et al., 2017) with particular attention paid to the volcanic Jeju Island Kim et al., 2013), China (Liu et al., 2012;Wang et al., 2015;Yang et al., 2015;Jiang et al., 2020), Japan (Uchiyama et al., 2000;Nakada et al., 2011), Indonesia (Adyasari et al., 2019b), Cook Island , and Australia (Robinson et al., 2006;Robinson et al., 2007;Sanders et al., 2012).
The delivery of nutrients by terrestrial groundwater has been plentifully addressed also at continental to global scales (e.g., Beusen et al., 2013;Sawyer et al., 2016;Luijendijk et al., 2020), typically by multiplying the estimated water flux by the nutrient concentration of the fresh groundwater to obtain the nutrient mass flux. Such estimates do not consider the biogeochemical reactions in the STE that modulate the actual inputs to the coastal ocean. Local-scale studies have revealed the critical influence of biogeochemical reactions in the STE on land-ocean solute fluxes and coastal ecology, even though STEs only cover a narrow stretch of the coast. This influence happens most certainly on a local and regional scale and is likely to have a significant effect on global scales. However, what remains unclear is how the knowledge and process descriptions at the local scale can be translated to the regional scale so that the impact of terrestrial SGD on the marine environment can be better quantified. We define the local scale as coastline lengths between 10 0 -10 2 m and the regional scale at 10 3 -10 5 m. Here, we summarize local scale knowledge about processes in STEs and outline perspectives toward the representation and extrapolation of these processes at the regional scale.

BIOGEOCHEMICAL CYCLING IN SUBTERRANEAN ESTUARIES
Subterranean estuaries are biogeochemically active, characterized by microbially mediated element release, fixation and transformation reactions (Charette et al., 2005;Spiteri et al., 2008a;McGrath et al., 2012;Beck et al., 2017;Duque et al., 2020). Biogeochemical reactions within the STE and the export of reaction products into coastal waters can have wide-reaching consequences for ecologically detrimental coastal processes such as eutrophication and associated oxygen depletion of bottom waters (Slomp and Van Cappellen, 2004;Dybas, 2005), modification of ocean derived substances (e.g., detrital matter including POC and S species: Shum and Sundby, 1996), buffering or enhancing ocean acidification (Wang et al., 2014;Santos et al., 2015;Liu et al., 2017), and transport of contaminants from the land to the sea (Brovelli et al., 2007;Wang et al., 2012). Depending on redox conditions and water flux rates, this can be crucial in determining element speciation and balances in coastal waters with effects for phytoplankton productivity and potentially eutrophication (Taniguchi et al., 2019). Increased primary productivity can further intensify the transfer of nutrients and possible contaminants from lower to higher trophic levels.
Organic carbon forms a primary electron donor in the coastal subsurface, which interacts with various electron acceptors (e.g., O, N, S, Fe and Mn species) derived from the sea or transported within the aquifer. Organic carbon stimulates a cascade of biogeochemical reactions relevant for local to global element turnover rates and budgets at the land-sea interface (Goodridge and Melack, 2014;Couturier et al., 2017;Cho et al., 2018;Seibert et al., 2019). Organic carbon quality, i.e., the fitness of organic molecules for microbial use as an electron donor (and sometimes acceptor), is dependent on its molecular composition, which is usually a complex mixture of terrestrial and marine-derived chemical moieties (Seidel et al., 2015).
The steep physicochemical gradients and enhanced microbial activity in STEs can also lead to mineralization processes of organic matter and hydrocarbons (Akam et al., 2020) as well as reoxidation reactions (Roy et al., 2013). These processes alter the redox conditions and the proton activity and may lead to the formation of authigenic solid phases, like carbonate, sulfide, sulfate or phosphate minerals (O'Connor et al., 2018;Riechelmann et al., 2020;Roy et al., 2013). Mineral formation may limit the concentrations of dissolved trace elements like metals being co-precipitated during iron sulfide formation (e.g., Huerta-Diaz and Morse, 1992). However, mineral formation also adds to the STE function as a filtering barrier between the terrestrial groundwater system and the sea for harmful chemical compounds. A well-known example of this is the socalled 'iron curtain,' an iron oxide rich zone that forms when oxic seawater is cycled through the beach sediments and mixes with anoxic or sub-oxic terrestrial groundwater (Chambers and Odum, 1990;Charette and Sholkovitz, 2002;Spiteri et al., 2006;Linkhorst et al., 2017;Sirois et al., 2018). The large sorption capacity of these iron oxides can be responsible for fixating aqueous pollutants, such as P and As species, before they become entwined in coastal food webs (Chambers and Odum, 1990;Bone et al., 2006;Beck et al., 2010).
Biogeochemically active areas can be unevenly distributed within the STE. This effect is amplified by the geological heterogeneity of flow regimes . Typically, certain conditions lead to the formation of biogeochemical 'hot-spots', areas of enhanced microbiological and chemical Frontiers in Earth Science | www.frontiersin.org May 2021 | Volume 9 | Article 601293 activity that drive the modification and transformation of ions and molecules in the subsurface (Santos et al., 2008). Microbial activity and community changes with salinity (Adyasari et al., 2019a). While the processes causing the observed spatial distribution of hot-spots are still poorly understood, an appropriate substrate (chemical and physical) is needed, as is an active flux of reactive electron acceptors and donors to maintain an energy flux for the microbial metabolism. Thus, knowing the hydrogeology of a site (water origin, flow paths, residence times, mixing processes, and aquifer properties) is crucial for understanding the distribution of hot spots in STEs (Goodridge and Melack, 2014). Modeling suggested that densitydriven flow can enhance mixing in an STE when dense saline water from the upper saline plume or periodic flooding sinks into the underlying freshwater lens due to 'salt fingering' (Greskowiak, 2014). Other studies show the importance of heterogeneity of hydrogeological parameters in controlling an STE's shape (Weinstein et al., 2007) and the associated biogeochemical reactions .
The temporal equivalent of 'hot-spot' is often termed 'hot moment'. The recognition that high element turnover (reactivity) areas are temporally dynamic is important Waska et al., 2019b). Hot moments imply that enhanced biogeochemical activity areas are stationary neither in time nor location and will fluctuate depending on hydrological and chemical boundary conditions (Seidel et al., 2015). The incorporation of hot moments into the framework of biogeochemical activity in the subsurface has helped understand the previously ambiguous distributions of elements and molecules, e.g., in the framework of groundwater-river interaction (Vidon et al., 2010). At the coast, storms may, for instance, physically impact the freshwater input, the availability of substrates, or the seafloor morphology, which would change STE behavior. Due to storms, SGD rates can strongly change in a matter of days (Cho et al., 2021). More generally, hot moments in STE can be created through movements of the salinity gradient, when, e.g., phosphorus can desorb from aquifer sediment with increasing salinity and cause a spike of P fluxes to the ocean that is not directly connected to terrestrial or marine inputs (Flower et al., 2017). Hot moments can also result from shifts in groundwater flow paths due to changed wave activity, which have been hypothesized to trigger pulses of arsenic transport into coastal waters due to instability of iron oxyhydroxides under changing redox conditions (Rakhimbekova et al., 2018). Water residence time is an essential factor influencing chemical reaction kinetics in the STE. For instance, in Waquoit Bay, longer residence times during winter caused attenuation of the groundwater nutrient load within the STE, whereas during periods of shorter residence time, SGD was nutrient-enriched .
The characterization of mineral authigenesis and its proxy potential in the STE requires sediment core material from the coastal and hinterland zone, best vertically resolved along transects through the STE to gain regional information. Whereas the physicochemical characterization of the aqueous solutions may allow for the identification of the current situation, the solid phases, when combined with appropriate age control, provide information about past SGD (Böttcher and Dietzel, 2010) as well as the movement of the geochemical mixing zone within the STE in response to sea-level change Hong et al., 2018). The sea-level rise since the last glacial maximum caused a landward movement of STEs. Isotopic composition of authigenic minerals documents past changes in the land-sea hydrological conductivity and the evolution of the paleoenvironment.
In the following, the biogeochemical cycling of individual compounds in STEs will be discussed.

Carbon
Subterranean estuaries receive carbon (C) from both terrestrial and marine sources. Terrestrial dissolved inorganic carbon (DIC) is principally derived from biological processes taking place in the soil, and further influenced by biogeochemical processes in the groundwater (Bogli, 1980;Clark and Fritz, 1997;Deines et al., 1974). Terrestrial dissolved organic carbon (DOC) in STEs is either transported from the soil zone of the surrounding watersheds to the coast via groundwater or leached from vascular-plant material buried in coastal aquifers, such as paleosols (Sirois et al., 2018) or peat layers. Marine DOC is supplied to STEs by infiltrating seawater and produced locally in the STEs through remineralization of particulate organic matter (POM), for example, from coastal phytoplankton and macrobenthic beach wrack (Kim et al., 2019). DOC from both sources is transformed in STEs by microbial activity into CO 2 and CH 4 .
Generally, terrestrial groundwater contributions will enhance the regional potential of coastal seawater to degas CO 2 to the atmosphere and contribute to primary productivity. However, DIC from terrestrial groundwater may also induce intense carbonate mineral dissolution (if terrestrial groundwater is oversaturated in CO 2 ) or lead to the precipitation of secondary carbonates (Deines et al., 1974;Wigley et al., 1978). The dissolution of carbonates in beach sands could be a substantial pH-buffer preventing acidic groundwater from changing coastal marine environments. Mixing between different water masses may impact the saturation states of carbonate minerals in an aqueous solution (Bogli, 1980;Hanshaw and Back, 1979), which amplifies the STE effect on the carbon cycle. Carbonate dissolution in the groundwater-seawater mixing zone can form caves and caverneous structures (Mylroie and Carew, 1990). At the regional scale, the total amount of carbonate dissolution through saturation effects of mixing groundwater and seawater was estimated on the Yucatan Peninsula (Hanshaw and Back, 1980). Therefore, carbonate dissolution structures may potentially be used to derive information about STE activity at regional scales or trace back STEs over geological time scales. Precipitation or dissolution of authigenic carbonate minerals could develop characteristic C and O isotope as well as trace metal signatures that can be used to deduce the composition of past solute gradients in an STE (Böttcher and Dietzel, 2010).
SGD can also be a source of methane (CH 4 ) to coastal waters (Lecher et al., 2016b;O'Reilly et al., 2015), and CH 4 gradients have been used to identify sites of SGD (e.g., Cable et al., 1996;Böttcher et al., 2021). In the Baltic Sea, seepage of terrestrial Frontiers in Earth Science | www.frontiersin.org May 2021 | Volume 9 | Article 601293 groundwater associated with CH 4 releases was observed (Whiticar and Werner, 1981;Whiticar, 2002). STEs, both in karstic and porous settings, can be hotspots for CH 4 oxidation that removes CH 4 from groundwater before discharging into surface waters (Schutte et al., 2016;Brankovits and Pohlman, 2020). High variability in temporal (seasonal) and spatial (among geographically similar STEs) CO 2 and CH 4 concentrations and fluxes have been observed (Pain et al., 2019;Pain et al., 2020), complicating extrapolations to the regional scale. Mineralization of DOM and CH 4 in STEs changes ratios of total alkalinity (TA) over DIC depending on the involved electron acceptors and superimposing biogeochemical processes (Akam et al., 2020). The mineralization can result in characteristic carbon isotope signatures (Deines et al., 1974;Meister et al., 2019). Sheltered, organic-rich STEs with sulfate reduction as a dominant remineralization pathway may become enriched with DOC and TA (e.g., Sippo et al., 2016). The mineralization of organic carbon in STEs is often linked to several overlapping mechanisms, which can vary over seasonal cycles, such as the advective flow of terrestrial groundwater and saltwater and the availability of organic carbon Liu et al., 2017;Kim et al., 2019). Marine-derived organic matter is assumed to be more abundant and more labile than terrestrially derived organic matter (Seidel et al., 2015) and is accompanied by an injection of oxygen, which accelerates mineralization rates. This "priming" by labile organic carbon and supply of oxygen may impact remineralization of groundwater-imported terrestrial DOC, but how and to which extent is weakly known. Since seawater circulation volumes exceed terrestrial groundwater throughputs in many STEs, marine DOC inputs could primarily determine the activity of the STE microbial reactor. Hence, most DOC released back into the coastal water column could be recycled/transformed. However, the relative contributions of terrestrial and marine organic carbon, as well as the function of the STE as a net source or net sink of organic carbon, are highly variable in space and time (Webb et al., 2019).
Carbon cycling in the STE on a local scale is best understood based on parallel profiling on land and at the seaside of the groundwater-seawater mixing zone since both terrestrial-and marine processes can control DIC and DOC fluxes. Besides direct tracking and benthic measurements of SGD leaving an STE (Donis et al., 2017), stable C isotope compositions of DIC, DOC and CH 4 are useful variables in the characterization and even quantification of internal processes and modulations in STEs (e.g., Winde et al., 2014;Donis et al., 2017;Meister et al., 2019;Pain et al., 2019;Pain et al., 2020). Also, DOM molecular traits can elucidate organic carbon cycling in STEs (Seidel et al., 2015). At the regional scale, the C cycle is particularly complex to represent because of the considerable heterogeneity of unconfined aquifers and the formation of metabolically diverse 'hot spots'. Therefore, an abundance of information is necessary both from the land-and marine sides of the STE.

Nitrogen
Subterranean estuaries may act as either a source or sink of nitrogen (N). Microorganisms catalyze different reactions of reactive N, which include nitrate (NO − 3 ), nitrite (NO − 2 ), and ammonium (NH + 4 ). The core nitrogen cycle involves four reductions [N fixation , assimilatory nitrate reduction (NO − 3 → NH + 4 ), and dissimilatory nitrate reduction to ammonium (DNRA, NO − 3 → NH + 4 ) as well as two oxidation pathways [nitrification (NH + 4 → NO − 3 )] and anaerobic ammonium oxidation (anammox, NH + 4 → N 2 ) (Kanehisa and Goto, 2000). The nitrogen cycle further includes organic nitrogen mineralization, or ammonification, which is the breakdown of organic N (N org ), leading to the release of NH + 4 . Nitrogen usually enters the system via N fixation (Fulweiler et al., 2007;Rao and Charette, 2012), organic input from terrestrial and marine sources, or terrestrial groundwater in the form of nitrate (Weinstein et al., 2011). STEs were classified as nutrient-contaminated STEs when they receive elevated amounts of nitrate from terrestrial groundwater and uncontaminated STEs when they produce nitrate internally from the mineralization of organic matter (Santos et al., 2009b;Loveless and Oldham, 2010;Robinson et al., 2018). In the latter case, tidal pumping provides N org to the STE, which later is remineralized into NH + 4 or nitrified to create NO − 2 and NO − 3 (Ullman et al., 2003;Santos et al., 2008;Charbonnier et al., 2013). The nitrifying microorganisms are aerobic chemoautotrophs; thus, nitrification generally occurs at the oxic surface layer of STEs where the ammonification rate is high, and N org input and oxygen concentration are abundant (Santoro et al., 2008).
It has been shown that nitrate may be removed in permeable intertidal sediments from circulating seawater (Marchant et al., 2016), in some locations by up to 70% (Wong et al., 2020). Permanent biological removal of reactive N (NO − 3 and NH + 4 ) from the STE may be achieved through biological assimilation, denitrification and anammox, provided that contact time between the aqueous and solid phase of an STE is sufficient. Denitrification is reportedly the primary biogeochemical process responsible for N loss in coastal and marine systems (Canfield et al., 2010) and performed by a wide range of bacteria and archaea, mostly heterotrophic microorganisms. Particulate organic carbon (POC) in the coastal aquifer was reported as a major control on denitrification in STE (Kim et al., 2020). Under conditions where NO − 3 inputs exceed the availability of carbon substrate for denitrification, DNRA may occur (Tiedje, 1988;Gardner et al., 2006), which, however, preserves N in the system as NH + 4 . Anammox, another process of biological N removal (Jetten et al., 1998), may occur in an anaerobic, NH + 4 -abundant, low organic matter environment (Sáenz et al., 2012). The cooccurrences of N transformation processes, such as simultaneous nitrification-denitrification, has been found in the STE due to the oxygen and organic matter stratification (Hays and Ullman, 2007;Erler et al., 2014) or rapid mixing of different water masses, creating both oxic and anoxic microzones (Uchiyama et al., 2000;Kroeger and Charette, 2008).
The N transformation processes are controlled mainly by groundwater residence time (Santos et al., 2008;, redox condition (Slomp and Van Cappellen, 2004), availability of electron donors (Santoro, 2010), and mixing rate of freshwater and saltwater (Kroeger and Charette, 2008). They have been detected and measured by geochemistry or molecular biology approaches. Stable isotopes (δ 15 N) and concentrations of reactive N help understand N behavior in STEs (Kroeger and Charette, 2008) and can also be manipulated as a tracer . Denitrification rates can also be calculated from STE sediment cores (DeSimone and Howes, 1996;Nowicki et al., 1999). Microbial N transformation can be determined by molecular detection, which involves extracting the microbial DNA from water or sediment samples (Santoro et al., 2006;Rogers and Casciotti, 2010;Hong et al., 2018;Adyasari et al., 2020).
Reactions in the nitrogen cycle in STEs are controlled primarily by microbiology. While at the local scale geochemistry or molecular biology assessment methods exist, regional scale tracers for these processes have yet to be identified. They could be indirectly inferred from STE physical properties (e.g., organic matter content and quality in host sediment, residence time, redox conditions) that control reaction rates and types. One example is the use of hyperspectral imaging to identify phytoplankton and CDOM concentration or total suspended solids (e.g., Brando and Dekker, 2003) and benthic information (e.g., Vahtmäe et al., 2020) in coastal waters if the signal is strong and clear enough given possible depth and turbidity. Lastly, if substantial quantities of gaseous N 2 are produced by denitrification, it could be detected by geophysical methods under favorable conditions.

Phosphorus
Phosphorous (P) inputs into terrestrial groundwater are mainly derived from fertilizers, waste and sewage (Slomp and Van Cappellen, 2004), to a minor extent from the mineralization of organic matter (Froelich et al., 1979) or released from minerals as geogenic sources (Kazmierczak et al., 2020;Tao et al., 2020). Geogenic P sources are mostly of local importance and depend on the aquifer mineralogy. In saline groundwater, P is usually derived from the mineralization of organic matter (Froelich et al., 1979). In karstic STEs, where groundwater can reach surface waters via rapid conduit flow pulse, inputs of P are associated with high-intensity rainfall events. In such a system in southern Java (Indonesia), elevated P concentrations were linked to events with high groundwater discharge rates, leading to exceptionally high P inputs into coastal waters (Oehler et al., 2018).
Dissolved P has a high affinity to adsorb at mineral surfaces of carbonates (Gaudette and Lyons, 1980;Burton and Walter, 1990) and iron oxyhydroxides (Einsele, 1936;van der Grift et al., 2014). These minerals can precipitate in STEs, which thereby would reduce the amount of P that is transported via SGD (Pain et al., 2020). In Waquoit Bay, P concentrations are 5-7 times higher in iron oxide rich sands than in the overlying surface sands, indicating how effectively processes in STEs can bind P (Charette and Sholkovitz, 2002).
Phosphorus concentrations along the salinity gradient of STEs usually behave non-conservatively. Elevated P concentrations are often associated with elevated groundwater salinity (Gaudette and Lyons, 1980), e.g., due to desorption from particles (Suzumura et al., 2000), Fe reduction associated with organic matter mineralization, direct P liberation (Froelich et al., 1979), or the release of colloid-bound nutrients (Prouty et al., 2017a). Sporadic P release can also follow occasional saltwater intrusions. Thus, to understand the temporal P transport behavior in an STE, the mobility of the groundwater-seawater mixing zone needs to be known. Depending on the distance to the sediment surface and the SGD flow regime, particulate and dissolved P may be transported to coastal waters (e.g., Lipka et al., 2018). Therefore, the mobility and final release of P to surface waters indirectly depend on redox conditions and the specific composition of the aquifer (Gaudette and Lyons, 1980;Lewandowski et al., 2015).
Phosphorous cycling in the STE is on a local scale best investigated based on the analysis of pore waters and sediments along transects across the terrestrial groundwater end-member through the mixing zone. Porewater analyses should be combined with continuous monitoring of the composition in the pelagic system and complemented by element budgeting. Care has to be taken in separating P in porewater derived from the groundwater against that derived from processes within an STE (Suzumura et al., 2000;Price et al., 2010;Prouty et al., 2017b). Redox conditions in the terrestrial groundwater, as well as the abundance of iron oxyhydroxides and carbonate minerals in the STE sediment need to be known or approximated to estimate the regional scale effect of STEs on P transport and release to the coastal ecosystem.

Sulfur
Subterranean estuaries gain most sulfur (S) from seawater, of which sulfate (SO 2− 4 ) is a major constituent. In brackish-marine sediments, microbial dissimilatory SO 2− 4 reduction is the primary anaerobic process responsible for the mineralization of organic matter (Jørgensen, 1982), leading to the formation of dissolved sulfide. Sulfate reduction is also responsible for the oxidative conversion of methane into dissolved carbon dioxide (Boetius et al., 2000). Nevertheless, STEs also receive S from terrestrial groundwater, where dissolved SO 2− 4 may originate from the dissolution of aquifer minerals like gypsum or oxidation of iron sulfides, like pyrite (Zhang et al., 2012). In young water bodies, anthropogenic sources, like acid rain, acid mine drainage and fertilizers, may further enhance SO 2− 4 loads (Clark and Fritz, 1997;Alorda-Kleinglass et al., 2019).
In STEs containing sufficient electron donors, SO 2− 4 reduction might be enhanced and can outcompete methanogenesis (Slomp and Van Cappellen, 2004). These systems can differ substantially in their biogeochemical processes and thereby in their release of climate-relevant gas emissions (Böttcher et al., 2021). Hydrogen sulfide may furthermore be re-oxidized by solid and aqueous species or be precipitated as iron sulfides (Luther et al., 1991;Rickard, 1997), thus acting as a sink for other dissolved metals (Huerta-Diaz and Morse, 1992). Overall, the dissolved SO 2− 4 availability controls the biogeochemical element cycling in the groundwater-seawater mixing zone and specifically the coupled sulfur-carbon-metal cycle. In particular, the coupled stable sulfur and oxygen isotope composition of dissolved SO 2− 4 in the STE pore waters provides information about sources, sinks and cycling of sulfur (Fritz et al., 1989;Zhang et al., 2012) that may even be traced in an impacted shallow water column. Using sulfur isotopes as a tracer will be difficult at the regional scale due to the many involved controls and high sampling effort necessary.

Metals
Dissolved metal concentrations in pore waters of STEs (e.g., Ba, Cd, Cu, Fe, Mn, Mo, Pb, Zn) depend strongly on the mineralogy of the aquifer and may be enriched by anthropogenic contamination (Bone et al., 2007;Knee and Paytan, 2011). Many of the trace metals have particular ecological relevance as (micro)nutrients or toxins (Salt et al., 1995;Beck et al., 2010), and their signatures may be recognizable in near-shore STE sediments (Knee and Paytan, 2011). Submarine groundwater discharge is relevant for marine alkaline earth metal composition (Mayfield et al., 2021), a significant source of rare earth elements to the ocean and could be the missing link in the global distribution of Neodynium (Nd), a key proxy for oceanic water-mass mixing (Johannesson and Burdige, 2007;Chevis et al., 2015;Paffrath et al., 2020). Formation of complexes and colloids, sorption on the surface of particles, ion exchange, and changes in speciation are some of the relevant processes that control metal mobility, reactivity and toxicity (Charette and Sholkovitz, 2006;Waska et al., 2019a). The interactions between solutes and aquifer particles exert a firm control on trace metal transport in STEs and are strongly influenced by redox and pH conditions, salinity, and ligand availability (Knee and Paytan, 2011).
One of the quantitatively most abundant metals in STEs is iron (Fe). Its mobility is highly sensitive to changes in redox milieu and pH (Spiteri et al., 2006). Upon precipitation, iron-oxyhydroxides (FeOOH) and iron sulfides (FeS) provide a substrate for sorption and incorporation of metals, Si, As, DOM, and especially P (Huerta-Diaz and Morse, 1992; van der Grift et al., 2014). Furthermore, a substantial impact of FeOOH precipitation on the fractionation of different DOM fractions was described (Linkhorst et al., 2017). Iron mineral formation in the STE may remove selected trace metals from the aqueous solution (Charette and Sholkovitz, 2002;Charette and Sholkovitz, 2006;Böttcher and Dietzel, 2010) and release them later, when the reactive front may have shifted. Karst systems, in particular, can be efficient pathways for the exfiltration of metal-enriched waters through an STE into the coastal zone due to short aquifer residence times, enhanced flow rates, and a low mineral surfaces to groundwater volume ratio (Knee and Paytan, 2011;Pain et al., 2020).
The use of trace metals to identify SGD in the coastal waters depends on their reactivity in redox and pH gradients. Mn, for instance, remains mobile for some time even under oxic conditions (Kowalski et al., 2012;Winde et al., 2014), whereas Fe is efficiently fixed in the STE. This stability difference leads to fractionation of these metals (Balzer, 1982) and the trace metals that sorb on them. On a regional scale, selected trace metals like Mn or Ba in the water column may help to detect and even quantify SGD in tidal areas (Kowalski et al., 2012;Winde et al., 2014).

Geochemical Tracers
Geochemical tracers can help to identify locations of STEs as well as processes within them. In general, tracers are characterized by distinct differences in concentrations or isotopic compositions between groundwater and seawater. Salinity or conductivity is a ubiquitous tracer since freshening of coastal ocean waters without a nearby surficial terrestrial source is an indicator of STE influence. However, traces of freshening diminish quickly due to mixing with seawater, and interpretations can be ambiguous due to other potential freshwater sources, such as rain.
Thus, other geochemical tracers with higher sensitivity, i.e., large concentration differences between groundwater and seawater, are applied. These are e.g. the natural radionuclides of radon ( 222 Rn; t ½ 3.8 days; Rn in following) and radium ( 224 Ra; t 1/2 3.7 days; 223 Ra; t 1/2 11.4 days; 228 Ra; t 1/2 5.75 years; 226 Ra; t 1/2 1,600 years, Ra in the following), which are generally enriched in STE groundwater compared to surface waters. They can provide qualitative and quantitative information on sources and types of water sources to the STE, estimates of seawater residence times in the STE, as well as an overall quantification of the water flux out of the STE (i.e., SGD) into the coastal sea and open ocean (Taniguchi et al., 2019), while the results often include large uncertainties (Rodellas et al., 2021).
In addition to imports by sea-and groundwater, a fraction of the Ra and Rn in STEs pore waters is produced locally in sediment mineral grains by the decay of the radioactive parents Thorium (Th) and 226 Ra, respectively. In pore waters with low chloride content Ra is typically immobilized by sorption to clay and Fe and Mn oxide surfaces whereas, at salinities > ∼3 practical salinity units (PSU), it behaves conservatively and can be used to trace water flow paths and fluxes (Webster et al., 1994). As an inert noble gas, Rn is not affected by the chemical composition of pore waters and geochemical reactions within the STE. The Ra and Rn pore water concentration will depend on the thorium/radium content of aquifer sediments, the recoiled efficiency, the residence time (i.e., the time elapsed since the water entered the STE) as well as Mn and Fe redox cycling. As seawater Ra and Rn concentrations usually are low, enrichment of these isotopes in coastal waters can be used to locate water fluxes out of the STE and quantify them in a mass balance approach assuming a known groundwater endmember (e.g., Garcia-Solsona et al., 2008).
A sampling of radionuclides can be done by boat at regional scales (Moore, 2000;Schubert et al., 2019;Stieglitz et al., 2010) (Figure 3). Time series measurements provide information on the temporal component of Ra and Rn fluxes and thus can be used to understand the dynamics of SGD in response to changing boundary conditions such as tides or storm events (Burnett and Dulaiova, 2003;Santos et al., 2009a). However, identifying STEs and associated terrestrial groundwater flux remains very challenging due to overlapping groundwater circulation processes and pathways, which are difficult to differentiate based on Rn and Ra only.
The quantification of SGD using Rn and Ra is based on solving a mass balance of these radionuclides, considering all their major sources and sinks. It is usually done on sub-regional scales of up to 100s of meters, with notable exceptions for the Atlantic Ocean (Moore et al., 2008), the Mediterranean Sea  and the global oceans (Kwon et al., 2014). The quantification can be done either assuming a steady-state or using a non-steadystate approach when influxes from the STE are transient Moore, 2003;Hwang et al., 2005;Burnett et al., 2006). The steady-state assumption needs to be critically appraised when analyzing field data used in the mass-balance calculation, especially when using Rn. Turbulent conditions caused by wind and wave action can lead to excessive degassing before and during field campaigns so that the assumption of steady state between source and loss terms is not met. Indeed, it can be difficult to measure any Rn in surface waters following rough conditions. The 'memory' effect of variable degassing for the Rn balance tends to be around 5 days, but depends on factors such as water depth. In general, a non-steady state mass-balance together with continuous Rn measurements is better suited to dynamic conditions in the coastal zone (e.g., Santos et al., 2009c;Gilfedder et al., 2015).
The concentrations of tracers in SGD, i.e., the end-member concentrations of Ra and Rn, must be known to translate a mass flux into a water flux. Constraining these end-member concentrations is one of the most critical and uncertain steps of mass balance approaches. Groundwater concentrations may have substantial temporal and spatial variations caused by the heterogeneity of parent isotopes, different groundwater circulation pathways, hydrological and marine forcing, changes of geochemical processes within the STE as well as variability in terrestrial sources Cho and Kim, 2016;Rocha et al., 2016;Cerda-Domenech et al., 2017). To overcome this problem, averaging a large number of end-member concentration measurements or taking the maximum and minimum end-member concentration to provide an envelope of possible SGD flux rates are common approaches (Moore, 1996;Beck et al., 2007). Nevertheless, the variability in the end-member concentrations remains the primary challenge for reliable quantification of SGD since shortcomings in accounting for the heterogeneity of end-members at local-and regional scales remain unsolved.
On a local scale, seepage meter measurements, as well as water balances, allow independent verifications of radionuclide-based FIGURE 3 | Example of a regional scale Rn survey in Queensland, Australia (Stieglitz et al., 2010). Permission for reprint was provided by the copyright-holder.
Frontiers in Earth Science | www.frontiersin.org May 2021 | Volume 9 | Article 601293 SGD rates (Povinec et al., 2012). However, on regional scales, verifications are still lacking or show systematic differences between tracer and model approaches (Prieto and Destouni, 2011). Nevertheless, Ra and Rn isotopes can be applied to map spatial and temporal changes in SGD, which under certain conditions would be a first indicator for the potential presence of an STE. Ra and Rn isotopes were used in regional scale SGD estimates. Another isotopic tracer that can be applied for STE studies is the stable isotope composition of water ( 2 H/ 1 H, 18 O/ 16 O, here shortly called "water isotopes"), which is a well-established tracer in hydro (geo)logical studies for the quantification of mixing processes and the deduction of water mass sources. Terrestrial groundwater is substantially more depleted in the heavier isotopes than seawater (Gat, 1996;Hoefs, 2018) and contains further information about the integrated meteorological conditions in the recharge area. Therefore, water isotopes are widely used in studies of SGD or STE (e.g., Burnett et al., 2006;Povinec et al., 2011;Rocha et al., 2016;Oehler et al., 2017;Duque et al., 2019). Time series can reveal the variability of terrestrial groundwater inputs on different time scales, particularly in dynamic karst STEs. Although still in its infancy, the extension of stable isotope characterization toward a triple isotope approach considering 17 O (resp. 17 O excess Sharp et al., 2018) may allow freshwater sourcing on the scale of regional STEs. Minerals formed in the mixing zone may record the parent solution's isotopic signature, providing opportunities for studies of past STE conditions.
Besides the stable isotopes of water, groundwater formed after the early 1960s may contain tritium ( 3 H), which has been used for dating purposes (Begemann and Libby, 1957;Bethke and Johnson, 2008). Meanwhile, the anthropogenic component, the atmospheric contamination of the water cycle by surface nuclear tests, has reached a natural background level. However, combined with its daughter isotope 3 He and Ne (Tolstikhin and Kamenskii, 1969;Sültenfuss et al., 2009), tritium still can be used for regional age estimation of the freshwater component in STEs (Röper et al., 2012;Post et al., 2019;Grünenbaum et al., 2020). The accumulation of 4 He is furthermore used in older groundwater dating (Bethke and Johnson, 2008). Other noble gas isotopes useful in the hydrogeology of STEs are 39 Ar and 85 Kr (e.g., Sánchez-Úbeda et al., 2018).
Local-scale geochemical investigation of STEs ideally includes vertical transects through the STE, a detailed physicochemical analysis, and solid-phase profiles from core material for a further microanalysis to identify precipitated authigenic solid phases. Measuring transects along the beach face enables mapping the spatial distribution of the biogeochemical environment, reactions and residence times over distances from tens to hundreds of meters. Larger scale mapping is challenging due to laborintensive field sampling and increasing heterogeneity of the subsurface at larger scales meaning that extrapolation between transects is difficult.
Different circulation pathways in the STE can also be examined by exploiting the different half-lives of Ra isotopes, ranging from tidal to seasonal time scales. Assuming a steadystate hydrological situation, residence times can be estimated applying a general advective transport model in which the concentration of Ra and Rn in pore waters is a function of the supply of the radionuclides and their loss by decay and advection Tamborski et al., 2017). One of the main uncertainties in this approach is the lithological composition (mineral content, grain size, porosity) of the STE, which is assumed to be homogenous in space and in time throughout the STE. In reality, however, the heterogeneity in the STE causes a wide range of radionuclide supply rates . Furthermore, temporally and spatially variable geochemical cycling of Fe and Mn and associated oxides profoundly influence the Ra pore water release . Since Ra is conservative in saline waters only, Ra-based residence times are only meaningful when seawater dominates the circulation through the STE. In contrast, Rn can be applied in both the fresh and brackish parts of the STE but is sensitive to 226 Ra distributions in the sediments as the parent isotope. Short term hydrological forcing like tides, wave set-up and induced circulation may also cause non-steadystate situations resulting in biases in residence times estimates. Finally, mixing of waters with different ages can lead to errors in the residence time estimates (which is similar to the real mean age of a water parcel) (Post et al., 2013;Gilfedder et al., 2019). As mixing is a linear process while radioactive decay and radioactive ingrowth are exponential (Bethke and Johnson, 2002), a systematic underestimation of the real mean age can occur if waters of different ages mix. One approach to circumvent this problem is the use of isotope ratios like, e.g., 224 Ra/ 223 Ra. This method is very accurate if one water mixes with another one having zero activity (e.g., negligible Ra in freshwater) but does not work if the age difference between the water masses is large.
Residence time measurements tend to be point measurements at specific times. Verification of residence times based on radioisotopes is complicated to obtain. Comparing Ra-and Rn-based residence times may be one approach (Tamborski et al., 2017). Alternatively, pore water residence times may be compared to those derived by the ratio of water volume in the STE and the water flux out of the STE (Colbert et al., 2008).
Ra and Rn isotopes can be combined with other tracers, such as temperature (Cranswick et al., 2014), salinity (Example in Figure 3), or dissolved silicon (DSi) (Waska and Kim, 2011;Oehler et al., 2019b). Since DSi behaves rather conservatively during transport through the STE, it can be used to trace processes that change other solutes in groundwater in STE, such as N or P (Oehler et al., 2019b). Another example is the combination of CO 2 measurements and Rn measurements (e.g., Cyronak et al., 2014;Santos et al., 2015;Macklin et al., 2019). While this ties CO 2 concentrations to SGD intensity, it does not yet differentiate if the CO 2 is a product of STE processes or if the coastal groundwater just has elevated CO 2 concentrations that originated from the terrestrial groundwater. Similarly, CH 4 was combined with Rn as a tracer of SGD (Cable et al., 1996;Dulaiova et al., 2010) but could also be used for tracing processes in STE that generate CH 4 . Using C stable isotopes as a tracer (Winde et al., 2014;Donis et al., 2017) and considering elemental mass balances (Deines et al., 1974) may be a successful approach.
The measurement of chemical tracers of STE processes at larger scales is associated with the problem to observe and measure the parameters with sufficient temporal and spatial resolution to resolve the underlying processes. Because tracer concentrations are modified upon emerging from the subsurface, sensors near the sea bottom are best suited to detect STE processes. The metabolites of typical benthic biogeochemical anaerobic OM degradation processes, like Mn 2+ and NH + 4 , may, after benthic-pelagic coupling, survive the oxidation within the oxic water column for some time (Kowalski et al., 2012;Winde et al., 2014). Since some liberated trace elements will be reoxidized and sorbed to suspended particles, regional investigations will have to consider both the dissolved and the solid phases (Kowalski et al., 2012). Whereas the dissolved substances may be detected using discrete or continuous sampling techniques (Petersen et al., 2011;Kowalski et al., 2012), solid-phase sampling for geochemical analysis is limited to discrete sampling. Measurements can also be performed by autonomous underwater vehicles, which could automatically trace signals, like temperature, conductivity or Rn-isotopes, of groundwater discharge and scout for STEs (Tholen et al., 2019).
New combinations of sensors would help to detect STEs. Optical spectrometers in the UV-visible range of the light spectrum can be used to determine the concentration of nitrate, nitrite, HS − , humic acids or DOC rapidly and without the use of reagents. It is to be expected that the optical absorption spectra will allow identifying water influenced by STEs. Sensors for the oxidation-reduction-potential (or, when referenced to the potential of a Standard Hydrogen Electrode, E h ) have been used for the estimation of concentrations of redoxsensitive elements in the redoxcline of the Baltic Sea (Meyer et al., 2014) or detecting hydrothermal plumes in the deep ocean (e.g., Baker et al., 2005). Combining several sensor types should allow the detection of STE signals in waters close to the bottom.
Geochemical tracers can help to regionally assess the amounts of SGD and point to the presence of STEs. They can also be applied at the local scale to draw quantitative conclusions about residence times in STEs and thus infer information about their reaction kinetics. At the regional scale, discrete and continuous measurements can be done. Still, these measurements will consist of discrete sampling points, and the feasibility of sampling will limit their resolution and coverage, which needs to be scaled according to the variability of the processes and tracers considered.

Ecological Tracers
Biodiversity in terrestrial groundwater recently received increased scientific attention (e.g., Hancock et al., 2005;Humphreys, 2008). Environmental gradients affect marine organisms in coastal waters. N and P are growth-limiting nutrients for most phytoplankton and macrophytes species that form the primary producers of the marine food web. Externally added nutrients alter the cycles of energy-flow between the pelagic and benthic zones, changing the community structure and population dynamics of both pelagic and benthic systems (Johannes, 1980;Sugimoto et al., 2017;Grzelak et al., 2018). A similar adaption to an increase in the concentration of nutrients can also be observed in the case of changes in other physical or chemical parameters of the environment: salinity, amount of light in the water column, amount of C species, heavy metals, O 2 concentration, the presence of hydrogen sulfide, or changes in pH. For example, observed changes in the chemical composition of bioavailable C in the Yucatan coastal waters result in lower coral cover, smaller size and reduced species richness (Crook et al., 2012). Benthic meio-and macrofaunal organisms may react to changes in the environment. Differences in benthic microbiology were observed in STE compared to purely saline beach waters, but little is known about microbial communities in STE (Santoro et al., 2006;Santoro et al., 2008;Santoro, 2010;Adyasari et al., 2019a;Adyasari et al., 2020). Thus, depending on location or environmental factors, STE may change species richness and diversity of meiofauna (Kotwicki et al., 2014;Encarnação et al., 2015;Welti et al., 2015;Grzelak et al., 2018). Similarly, macrofaunal communities can respond; their biodiversity can increase due to food supply (Waska and Kim, 2010;Pisternick et al., 2020) or may be reduced due to salinity stress or changes in pH (Zipperle and Reise, 2005;Migné et al., 2011;Utsunomiya et al., 2017). However, while the named ecological responses can highlight an STE presence, they will be hard to interpret as tracers of processes within that STE.
Since the impact of microbially mediated biogeochemical processes in the physicochemical gradients of STEs can be diverse (e.g., Beck et al., 2011), it is useful to apply the established methods to study the communities of organisms and associated biogeochemical processes occurring in the entire ecosystem. Methods of assessing abundance, activity, biomass and diversity include taxonomic phytoplankton analyses in water samples, fauna and microbial activity analyses in sediment cores, or fish assessment in the water. For ecological studies on the impacts of processes in STEs, these methods are combined with thermal infrared cameras or Ra/Rn tracing Grzelak et al., 2018). Simultaneous measurements of 222 Rn with other biological parameters allow assessing the reaction of biota activity on products of STEs (Taniguchi et al., 2019). Stable isotopes of N and C can also be used to trace STE products in the food web. Several studies have traced the behavior and fate of groundwater DIN using the isotope δ 15 N of seagrass and macroalgae (Winde et al., 2017;Andrisoa et al., 2019) and gradients in the stable carbon isotope ratio between fresh and marine end-members (Gramling et al., 2003;Winde et al., 2017).
Consistent integration of hydrographic and biogeochemical research with bio-monitoring, elementary approaches and stable isotope analyses while considering higher trophic levels like fish and trophic flows in coastal ecosystems would enable the estimation of the direct impact of processes in STEs on ecosystems.
While not yet done, specific species may have the potential to indicate STE processes in specific environments. Research in this direction is emerging (Lecher and Mackey, 2018), but additional knowledge is required before ecology can be used as a tracer for STE processes.

Geophysical Methods
From a geophysical perspective, the changes of physical and geochemical conditions under the seafloor associated with the presence of STEs are detectable by measuring properties such as electrical conductivity (EC), acoustic impedance (AI), or temperature. STEs can affect 1) the seafloor itself, e.g., its morphology, 2) the sediment or rock below it and 3) the properties of the water column, e.g., turbidity or sea-surface temperature.
At the local scale, geophysical methods can delineate the shape of an STE in detail. The EC difference between the terrestrial groundwater and saltwater is the most prominent parameter for local scale geophysical STE observation. The EC of freshwater varies depending on the total amount of dissolved ions, and literature puts it below 300 μS/cm (Kirsch, 2006), 1,000 μS/cm (Jiao and Post, 2019), or 2000 μS/cm (Langguth and Voigt, 2004). The EC of surficial seawater is orders of magnitude higher, > 40,000 μS/cm, depending on the temperature and thus latitude (Tyler et al., 2017). However, mineral grains, clay aggregates and organic matter in the STE also affect EC and need to be considered. Particularly sediment porosity and clay content strongly influence the overall EC of the bulk (rocks plus groundwater) volume (Archie, 1942). The high EC of clay leads to the general challenge of distinguishing freshwater saturated clays from brackish water saturated sands that may have very similar EC values. Gamma radiation of potassium, thorium and uranium at the surface helps to distinguish between clayey and sandy sediments (Siemon et al., 2020), as does an acoustic impedance contrast that seismic methods can show.
Many approaches that use electromagnetic fields (EM) emitted and detected by coils are available to represent subsurface properties (frequency-domain EM, time-domain-EM and slingram/ground-conductivity meter: Kirsch, 2006). Frequencydomain (sinusoidal transmitter currents) and time-domain EM (on/off transmitter currents) are typically used to obtain depth resolution at a single measurement point or along transects (Auken et al., 2003). Moving ground-conductivity meters are used for fast lateral mapping of EC but with limited depth resolution (e.g., De Smedt et al., 2013). Locating freshwater and thus STEs, at the regional scale can be realized using airborne or ship-based EC measuring methods.
Airborne methods use helicopters or fixed-wing aircraft to carry the geophysical systems at about 20-100 m above ground, covering 100-200 km of profiles or typically 10-50 km 2 per hour (Figure 4). In airborne EM, frequency-domain systems focus on near-surface (1-100 m) investigation of the spatial bulk EC (Siemon et al., 2015;Siemon et al., 2019;Siemon et al., 2020), whereas time-domain systems enable some deeper (5-500 m) investigation Steuer et al., 2009), depending on the conductivity distribution of the subsurface. Apart from electrical conductivity or resistivity, airborne radiometry maps the gamma radiation of the upper few decimeters of the surface (Wilford et al., 1997;IAEA, 2003). This technique has not yet been applied toward Ra/Rn concentrations in seawater but might be a pathway forward to map those isotopes with broader spatial coverage. Besides, semiairborne methods using a ground-based transmitter and a helicopter-borne receiver have been developed recently to increase the penetration. While deep penetration is less important for STE investigation, UAV-based concepts can supplement the helicopter by a drone carrying the EM sensor in the future and hopefully enable cheaper investigations at the several kilometer scale. Finally, ship-based EM measurements have been applied to detect SGD (Müller, 2010;Müller et al., 2011).
Another EC method is electrical resistivity tomography (ERT), which, on land, uses steel electrodes pinned to the ground measuring the resistivity (the reciprocal of electrical conductivity) of the sediment by injecting a direct current (Stieglitz, 2005;Swarzenski et al., 2006). ERT has been used to observe the saltwater recirculation zone (Morrow et al., 2010) and its temporal changes (e.g., Johnson et al., 2015;Sutter and Ingham, 2017). ERT also can detect induced polarization effects caused by iron oxides (Mansoor and Slater, 2007). While this has not yet been applied to STEs, it could analyze Fe-cycling in the mixing zone. Measurements of magneticresonance also provide information on iron oxides and their concentrations (Keating and Knight, 2007;Costabel et al., 2018) and provide subsurface porosity and hydraulic conductivity (e.g., Müller-Petke and Yaramanci, 2015).
Marine ship-based measurements of ERT tow floating or submerged streamers behind a boat (Manheim et al., 2004;Day-Lewis et al., 2006;Hermans and Paepen, 2020). A profile of 5 km length can be measured in 1 h, but the conductive seawater limits penetration depth and therefore, streamer lengths of several hundred meters are needed to achieve sufficient penetration depth.
Seafloor morphological SGD proxies, such as pockmarks (e.g., Hoffmann et al., 2020), submarine terraces (Jakobsson et al., 2020), and depressions in carbonate rocks (Oehler et al., 2019a), are visible in the bathymetry. There are several techniques available to map the seafloor. Sidescan sonar data can reveal the location of pockmarks (e.g., Schlüter et al., 2004;Virtasalo et al., 2019), and echo sounders generally can also be applied for this task (Feldens et al., 2018;Papenmeier et al., 2020). Multibeam echo sounders provide high lateral coverage across-track on the seafloor, which is several times the water depth, enabling a fast mapping of larger scales.
Obtaining information on the geological structures by using high-resolution seismic profiling can be relevant for delineating the geometry of STEs (Mosher and Simpkin, 1999). The resulting image shows the reflection amplitudes and times (or depths) resulting from material interfaces with an impedance contrast, where the acoustic impedance is the product of bulk density and compressional velocity. Seismic profiles can be measured simultaneously with geoelectric streamers providing the same estimated profile length of 5 km per hour.
Geophysical methods also allow analyzing biogeochemical processes that occur within STEs, since some products can be identified. Shallow free gas accumulations (e.g., CH 4 , N 2 ) can be mapped on larger scales using multibeam echo sounders. Gas is visible as acoustic blanking, bright spots, acoustic turbidity, or gas chimneys in seismic data (Hovland and Judd, 1988;Judd and Hovland, 1992). Free gas can have multiple origins in the ocean, so SGD does not necessarily cause its occurrence. However, gasrelated features along the Belgian coast (Missiaen et al., 2002), in Brazil, Argentina and South Africa (Weschenfelder et al., 2016) and the North Yellow Sea (Wang et al., 2018) were linked to localized groundwater discharge.
Finding transfer functions linking discrete biogeochemical data and continuous geophysical methods is a promising methodology to extrapolate STE processes to the regional scale. In particular, methods going beyond the "classical" mapping of resistivity (Figure 4), which would allow to spot offshore freshening of pore water and thus STEs, have only rarely been applied to characterizing biogeochemical processes in STEs but show a substantial potential. A critical challenge is integrating local information into regional mapping by connecting different scales (point-like measurements, ground-based spatial information, airborne and satellite). Recent developments allow to include point-like or drilling information into ground-based geophysics using different concepts (e.g., Wunderlich et al., 2018) and coupling ground-based geophysics with airborne data (Dickson et al., 2014), providing perspectives for regional-scale assessments that are unviable based on ground-based measurements alone.

Remote Sensing
Remote sensing provides spatially continuous information on a scale from 10 1 to 10 9 m 2 (Böttcher et al., 2021), depending on the platform (ground-based application, remotely operated vehicle, uncrewed erial vehicle (UAV), kite gyrocopter, airplane, satellite). The chosen scale constrains the ground resolution (ground sampling distance), which may vary between 10 0 to 10 2 m. Scale and ground resolution intrinsically determine which STE indicator may be observed (Böttcher et al., 2021). Diffuse discharging groundwater is unlikely to be observed using satellite-based remote sensing, but can be observed with close-range applications or UAVs. However, the effort to cover regional scales is unequally larger for UAVs than for airplane-or satellite-based applications. Given the negligible depth penetration of remote sensing earth observation techniques in most situations, remote sensing can only investigate STE processes indirectly through its surface expressions. Remote sensing can detect terrestrial SGD, map seafloor and coastline morphology to find potential STEs, indicate hydrogeological STE characteristics, and identify STE processes by, e.g., classifying submerged aquatic vegetation that may be influenced by STEs.
Available sensors cover a wide range of the electromagnetic spectrum and can detect/differentiate surface properties that may act as indicators for STE processes or characteristics. The most commonly applied sensor in STE investigations is the thermal infrared (TIR) sensor, which can exploit temperature differences between seawater and terrestrial groundwater at the sea-surface ( Figure 5). TIR detects groundwater inflow locations by identifying a thermal anomaly based on mono-temporal investigations (Fischer et al., 1964;Mejias et al., 2012;Wilson and Rocha, 2012;Kelly et al., 2013;Mallast et al., 2014;Xing et al., 2016), or multi-temporal investigations (Schubert et al., 2014;Oehler et al., 2018), and can be used to quantify freshwater fluxes given in-situ reference data of currents and bathymetry (Roseen, 2002;Johnson et al., 2008;Danielescu et al., 2009;Tamborski et al., 2015). The extent and shape of sea-surface temperature anomalies can indicate specific STE processes and characteristics (Chen, 1991;Jirka, 2004). Elongated anomalies oriented perpendicular to the coastline (Shaban et al., 2005) or extensive areas of thermal anomalies (Kelly et al., 2013) point to focused groundwater discharge at a high rate, suggesting a low residence time in STEs. Areas with smaller anomalies close to the coastline distributed over tens of meters alongshore (Tamborski et al., 2015) suggest diffuse discharge. The former is typical for karst or volcanic aquifers from which groundwater discharges along conduits, whereas the latter is commonly associated with sedimentary aquifers. Since higher temperature differences and distinct anomalies are easier to measure by TIR sensors, focused, high-volume discharge is more easily detected than slow, diffuse discharge. Besides thermally identifying terrestrial groundwater discharge, remote sensing can classify benthic microbiology (e.g., Kazemipour et al., 2012) and submerged aquatic vegetation, which has been found to allow inferences about STEs (e.g., Leitão et al., 2015;Welti et al., 2015). While these observations are connected to local scale in-situ investigations, optical remote sensing using multi-and hyperspectral sensors (Klemas, 2016) has successfully been applied in shallow waters but not yet linked to potential terrestrial groundwater influence.
The relationship between STE characteristics indicated by seabed morphology and even operationally provided seafloor maps over regional scales using air-and satellite-borne platforms is similarly unexploited (Siermann et al., 2014;Eugenio et al., 2015). Remotely sensed bathymetry would allow a first-order approximation of potential pockmark sites. Since one of the forces that can produce pockmarks is fluid flow (Hovland et al., 2002), which can be an STE defining process, pockmarks can be seen as hinting toward the potential existence of STEs. Coastline morphology can be shaped by the same processes and may represent a similar indicator (Johannes, 1980). A drawback of remote sensing methods is that it is challenging to link discrete, in-situ measurements to remote sensing data that integrate a multitude of effects other than those caused by an STE presence (e.g., ocean currents, chlorophyll presences). This weakness calls for a multi-scale approach (Lausch et al., 2013) in which an identified STE location is observed synchronously and over long periods with an appropriate in-situ sensor array and various airborne to satellite platforms covering different spatial scales. Such long-term sites could be associated with coastal observatories (Schofield et al., 2003) equipped with sensors to measure, e.g., turbidity, EC, pH, oxygen level, and thus be investigated through interdisciplinary approaches (Mollenhauer et al., 2018).
Second, since the temporal variability of the STE-associated processes is high, an interesting approach would be to have a TIR remote sensing approach that enables a multi-temporal observation possibility. The recent advance of UAVs do provide such a possibility, as shown in Mallast and Siebert (2019), but cannot cover a regional scale simultaneously. The revisiting times will be reduced for Landsat with the launch of Landsat-9 in spring 2021 and the Indian-French THRISHNA mission, whose launch is planned for 2024 (Lagouarde et al., 2018). Nano-to mini-satellites can even revisit individual FIGURE 5 | TIR map displaying an area of cool diffuse seepage of terrestrial SGD that indicates an STE (Tamborski et al., 2015). Permission for reprint was provided by the copyright-holder.
locations several times per day (Van Ryswyk, 2020). If these vehicles would also carry a TIR sensor and resemble a similarly very high ground resolution of <3 m, a regional to global scale localization of SGD sites and thus STE locations could become possible.

Hydrogeological Modeling
Another way to investigate STEs is through reactive transport models (RTMs), which remains underutilized in STE research (Robinson et al., 2017). Published modeling studies analyzed STEs as a high reactivity zones for, e.g., organic contaminants, nitrate and sulfate (Robinson et al., 2009;Anwar et al., 2014;. The extent to which reactions proceeded was dependent on the degree of dispersive mixing and residence time in the STE, which was controlled by hydrological boundary conditions, such as terrestrial groundwater flux, beach slope and tidal amplitude. Another modeling study also showed the importance of dispersive mixing in STE on organic contaminant removal from groundwater (Nick et al., 2013).
The RTMs mentioned above provided process understanding but lacked confirmation by field data. There are only two field sites where RTMs and observations were integrated. The studies at Waquoit Bay by Spiteri et al. (2008a), Spiteri et al. (2008b), and Spiteri et al. (2008c), formed the earliest applications of RTMs of the STE and showed that the pH change at the interface between terrestrial groundwater and the lower saltwater wedge is an essential factor that controls the oxidation of Fe 2+ and subsequent adsorption of land derived phosphorus (Spiteri et al., 2008a). Further, model results suggested that enhanced nutrient turnover occurs close to the SGD exit point due to high transverse dispersive mixing of oxic and anaerobic groundwater, resulting from converging groundwater flow toward the sea (Spiteri et al., 2008b;Spiteri et al., 2008c). The other site where model simulations were used to support field data interpretation is at Cape Henlopen, Delaware, USA . At this site, the field data showed that aerobic respiration led to a depletion of oxygen within the upper saline plume and that nitrate concentrations were lowest where ammonium and particulate organic carbon were highest. Even though their numerical model was able to simulate the general spatial O 2 and N 2 concentration trends, discrepancies remained evident, which the authors attributed to uncaptured heterogeneity and transient processes, as well as the omission of iron and sulfate reduction processes in the model. Thus, this is an example of the difficulties to model STEs due to their high spatial complexity. Modeling can also be applied in combination with isotopic data to elucidate reliable dating information on the STE water age (Bethke and Johnson, 2008;Post et al., 2019), which could provide information for reaction rates.
It is unlikely that numerical, reactive transport models of biogeochemical reactions in STEs can be applied with any reasonable expectation of success at the regional scale. The first reason is that the computational resources required will remain prohibitive in the foreseeable future. The second and perhaps foremost reason is that the input data requirements are unattainable at the regional scale due to the high spatial and temporal variability of STEs. Therefore, regional-scale models of STE nutrient transformations will have to rely on simplified process representations that can be parametrized with input data obtainable at a reasonable effort and cost. There are several essential questions to consider, including: What is an appropriate spatial level of heterogeneity for capturing regional scale processes that could be similar to the representative elementary volume (REV) concept for groundwater flow (cf. Freeze, 1975)? How can processes be integrated along the vertical dimension? How to account for local-scale variability? What is the temporal variability at the intra-and inter-annual scale? Can coastal landscape units be classified according to their geochemical reactivity?
Concerning the last question, one feasible approach to upscale modeling of STE processes might be to identify distinctive zones where controlling factors on the reactions are relatively uniform and then use small-scale model information to extrapolate to a regional-scale area. Working with distinctive uniform zones would imply a classification (e.g., Bokuniewicz et al., 2003). However, levels of homogeneity representative for larger areas are not yet adequately defined, as is the answer to whether and how small-scale temporal variations extrapolate into a larger scale. The terrestrial component of SGD was already modeled at the regional scale (Jarsjö et al., 2008;Befus et al., 2017;Hajati et al., 2019) and global scale (Zhou et al., 2019;Luijendijk et al., 2020) using similar approaches. These models produce spatially explicit estimates of terrestrial SGD that is a necessary prerequisite for STEs. Local scale RTMs highlight the importance of understanding mixing (hydrodynamic dispersion) and sediment organic carbon content, so it would seem reasonable to expect that these at least will have to be accounted for at the regional scale as well. Both are linked to sediment heterogeneity (regarding hydraulic properties and geochemistry) and the variability of hydraulic gradients so that the starting point would be a deterministic understanding of the longshore geological variability and a mechanistic description of the driving forces.
Data products of sufficient quality that represent controls of STE processes are required to model STE processes on the regional scale. The amount of data available at a global scale and high resolution is quickly increasing, but many aspects are not yet adequately covered. Foremost, since STEs occur along the coast, a high-resolution coastline is necessary and available (Sayre et al., 2018). STE sediment properties can be represented by recent datasets focusing on coastal sediment heterogeneity (de Graaf et al., 2017;Zamrsky et al., 2020). Other datasets focusing on the terrestrial aquifer (Gleeson et al., 2014) or marine sediment thickness (Straume et al., 2019) and its properties (Dutkiewicz et al., 2016) are hardly useable at the coast due to their global focus and low spatial resolution. Additional hydrogeological background knowledge would be, e.g., tectonic activity or stress maps, where areas of high tectonic activity would contain more faults and thus more potential pathways for groundwater.
Mixing intensity and processes in STEs will also be controlled by the amount of freshwater discharging (Luijendijk et al., 2020) or the amount of seawater recirculating (Mayfield et al., 2021). Water quality of the incoming freshwater is not known at a global scale, but, as a first crude guess, could be inferred by proxy, e.g., from land-cover (e.g., Arino et al., 2007) or soil (FAOIIASA, ISRIC et al., 2009) datasets. Another option would be to extrapolate from datasets of multiple local groundwater quality measurements (e.g., NAWQA from USGS). In any case, groundwater quality is far less constrained than aquifer properties at large scales. Synthesized data sets of geological and geochemical properties of the subsurface below the seafloor on a regional scale are virtually absent.
Any of the large-scale datasets available have problems resolving the spatial heterogeneity of coasts, so the resulting estimates will not be applicable at the local scale  but provide first-order estimates for larger scales. This problem is not unique to the coast (e.g., Moosdorf et al., 2010) but amplified by the steep gradients in that setting. Given that the uncertainty of any model-based estimate is going to be very large, rather than focusing on the accuracy of the absolute numbers, it would be preferable to develop approaches that reveal the lateral differences. Such an approach would allow the ranking of different areas in terms of their effect on the coastal systems, highlighting hotspots to be targeted by more detailed investigations or specific measures.

CONCLUSIONS AND PERSPECTIVE
This review of biogeochemical processes in subterranean estuaries highlighted the spatial and temporal variability of these systems. Their reactions and corresponding rates depend on the hydrology of the STEs and the chemical composition of sediments and source waters (terrestrial groundwater and seawater). These controls have to be constrained to budget the processes within STEs and estimate the regional scale effects of STEs.
Characterizing STEs at the regional scale can be achieved through their products ( Figure 6). The first step of estimating STE activity at a regional scale is to identify locations, amount, and the general extent of STEs. Locations can be identified through continuous measurements of geochemical tracers, e.g., electrical conductivity or Ra/Rn isotopes in the ocean.
Biological indicators can point to the existence of STEs. Geophysical methods can provide information about freshened pore water and gas occurrences below the seafloor. Hydrogeological modeling can also estimate the flow behavior, but the necessary data are usually not available to represent larger areas in numerical 3D models realistically. Remote sensing information, for example, TIR images, can additionally show thermal signatures of freshwater influx. Combining observational methods from different disciplinary backgrounds, paired with modeling and hydrogeological data, is most promising for a regional scale effort to locate STEs.
Products of biogeochemical processes can be used to trace processes occurring in STEs ( Figure 6). These products can be specific chemicals, like, e.g., CO 2 , N 2 gas or Mn 2+ , but also isotopic tracers like δ13C in DIC, δ18O and δ15 N in nitrate, as well as δ2D and δ18O in water. With geophysical methods, gaseous products of STEs can be identified, as well as iron oxides in the sediment, which are an essential control on the biogeochemistry of STEs. Remote sensing allows tracing changes in sea bottom that could indicate properties of the STEs in some cases. Nevertheless, since processes in STEs are spatially and temporally highly variable, discrete in-situ observations of specific STE products (e.g., CO 2 , N 2 , or Mn 2+ ) need to be combined with temporally and spatially continuous and temporally higher resolved data (e.g., remote sensing) using transfer functions.
Reactive transport models (RTM) can act as temporal and spatial link between discrete local scale in-situ measurements and spatially continuous regional scale data to extrapolate STE processes onto larger scales. To account for the high variability, quantitative local scale measurements will have to be combined with the concept of representative elementary volume (REV) based on regional-scale data sets such as from geophysical measurements or maps containing hydrogeological information. An RTM would represent each REV to quantify fluxes at scales relevant to the coastal zone and coastal zone management. On that basis, it would be even possible to sort STEs in classes when discussing their regional impact. Particularly, STEs with diffuse porous flow need to be distinguished from STEs in conduits, e.g., in karstic or volcanic settings.
FIGURE 6 | Consequences of hydro (geo)logical and biogeochemical processes in STEs that can be used for extrapolation to the regional scale.
Frontiers in Earth Science | www.frontiersin.org May 2021 | Volume 9 | Article 601293 A combination of methods from the presented disciplines promises a pathway toward regional scale estimates of processes in STEs and their impact on land-ocean matter fluxes. That impact may be substantial since local scale studies strongly hint toward their relevance. The processes are known; the products are known; the methods are available. Now is the time to combine methods across disciplines and understand these critical processes.

AUTHOR CONTRIBUTIONS
NM conceived the paper and wrote about upscaling; NM and UM coordinated the paper; UM, HW, and NM got funding for the project; MEB, BG, DA, A-KJ, TO, RP, JS, CE, and HW, contributed the biogeochemical part; LK contributed the biological part; EE, MM-P, BS, and TW, contributed the geophysical part; JG, GM, VP, and MW, contributed the hydrogeological modeling part; UM contributed the remote sensing part. All authors contributed to writing the manuscript and discussing the concept and formulating the overarching themes.

ACKNOWLEDGMENTS
This publication is a result of the DFG supported KiSNet project (MA7041/6-1). MEB and AKJ wish to thank DFG for financial support during research training group BALTIC TRANSCOAST (GRK 2000) and MEB and CvA the DAAD for a Ph.D. stipend for a stay of CvA at IOW. This is BALTIC TRANSCOAST publication No. GRK 2000/0045. HW was funded by the Niedersächsisches Ministerium für Wissenschaft und Kultur (MWK) in the scope of project "BIME" (ZN3184). JS acknowledges the support through the SEAMOUNT BONUS project (art. 185), which is funded jointly by the EU and the Federal Ministry of Education and Research of Germany (BMBF, Grant Nos. 03F0771B). The authors acknowledge the helpful review comments by D. Brankovits, G. Chaillou and A. Pain, and the editorial handling by F. Frappart. Finally, the authors would like to express their gratitude to Bill Burnett, who turns 75 this year. Bill inspired many of us to pursue research in the field of SGD and supported wherever he could. Thank you, Bill.