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 (Bögli, 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₂ and CH₄.

Generally, terrestrial groundwater contributions will enhance the regional potential of coastal seawater to degas CO₂ 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₂) 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 (Bögli, 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₄) to coastal waters (Lecher et al., 2016b; O'Reilly et al., 2015), and CH₄ 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 groundwater associated with CH₄ releases was observed (Whiticar and Werner, 1981; Whiticar, 2002). STEs, both in karstic and porous settings, can be hotspots for CH₄ oxidation that removes CH₄ 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₂ and CH₄ 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₄ 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 (Kim and Kim, 2017; 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₄ 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.

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 ( N 2 → NH 4 + ) , denitrification ( NO 3 − → N 2 ) ], 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 co-occurrences 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; Gonneea and Charette, 2014), 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 (δ¹⁵N) and concentrations of reactive N help understand N behavior in STEs (Kroeger and Charette, 2008) and can also be manipulated as a tracer (Erler et al., 2014). 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₂ are produced by denitrification, it could be detected by geophysical methods under favorable conditions.

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.

Subterranean estuaries gain most sulfur (S) from seawater, of which sulfate ( SO 4 2− ) is a major constituent. In brackish-marine sediments, microbial dissimilatory SO 4 2− 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 4 2− 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 4 2− loads (Clark and Fritz, 1997; Alorda-Kleinglass et al., 2019).

In STEs containing sufficient electron donors, SO 4 2− 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 4 2− 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 4 2− 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.

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 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 (²²²Rn; t_½ = 3.8 days; Rn in following) and radium (²²⁴Ra; t_1/2 = 3.7 days; ²²³Ra; t_1/2 = 11.4 days; ²²⁸Ra; t_1/2 = 5.75 years; ²²⁶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 ²²⁶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 (Rodellas et al., 2015) and the global oceans (Kwon et al., 2014). The quantification can be done either assuming a steady-state or using a non-steady-state approach when influxes from the STE are transient (Burnett et al., 2003; 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 (Gonneea et al., 2013; 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 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 (²H/¹H, ¹⁸O/¹⁶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 ¹⁷O (resp. ¹⁷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 (³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 ³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 ⁴He is furthermore used in older groundwater dating (Bethke and Johnson, 2008). Other noble gas isotopes useful in the hydrogeology of STEs are ³⁹Ar and ⁸⁵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 labor-intensive 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 steady-state 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 (Bokuniewicz et al., 2015; 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 (Beck and Cochran, 2013). Furthermore, temporally and spatially variable geochemical cycling of Fe and Mn and associated oxides profoundly influence the Ra pore water release (Gonneea et al., 2008). 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 ²²⁶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-steady-state 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., ²²⁴Ra/²²³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₂ measurements and Rn measurements (e.g., Cyronak et al., 2014; Santos et al., 2015; Macklin et al., 2019). While this ties CO₂ concentrations to SGD intensity, it does not yet differentiate if the CO₂ is a product of STE processes or if the coastal groundwater just has elevated CO₂ concentrations that originated from the terrestrial groundwater. Similarly, CH₄ 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₄. 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²⁺ 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 redox-sensitive 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.

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₂ 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 (Sugimoto et al., 2017; Grzelak et al., 2018). Simultaneous measurements of ²²²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 δ¹⁵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.

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). Frequency-domain (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² 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 (Siemon et al., 2009; 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, semi-airborne 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 magnetic-resonance 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₄, N₂) 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, gas-related 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 provides spatially continuous information on a scale from 10¹ to 10⁹ m² (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⁰ to 10² 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 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.

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; Kim et al., 2017). 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²⁺ 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 (Kim et al., 2017). 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₂ and N₂ 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 (Geng et al., 2020) 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.