Radium isotopes as submarine groundwater discharge (SGD) tracers: Review and recommendations

Submarine groundwater discharge (SGD) is now recognized as an important process of the hydrological cycle worldwide and plays a major role as a conveyor of dissolved compounds to the ocean. Naturally occurring radium isotopes ( 223 Ra, 224 Ra, 226 Ra and 228 Ra) are widely employed geochemical tracers in marine environments. Whilst Ra isotopes were initially predominantly applied to study open ocean processes and fluxes across the continental margins, their most common application in the marine environment has undoubtedly become the identification and quantification of SGD. This review focuses on the application of Ra isotopes as tracers of SGD and associated inputs of water and solutes to the coastal ocean. In addition, we review i) the processes controlling Ra enrichment and depletion in coastal groundwater and seawater; ii) the systematics applied to estimate SGD using Ra isotopes and iii) we summarize additional applications of Ra isotopes in groundwater and marine studies. We also provide some considerations that will help refine SGD estimates and identify the critical knowledge gaps and research needs related to the current use of Ra isotopes as SGD tracers.


Introduction
Radium (Ra) is the element number 88 in the Periodic Table and it belongs to Group IIA, the alkaline earth metals. There are four naturally occurring radium isotopes ( 224 Ra, 223 Ra, 228 Ra and 226 Ra), which are continuously produced by the decay of their Th-isotope parents of the U and Th decay series ( 228 Th,227 Th,232 Th and 230 Th, respectively). The most abundant isotopes are 226 Ra from the 238 U decay chain, an alphagamma emitter with a half-life of 1600 y, and 228 Ra from the 232 Th decay chain, a beta emitter with a half-life of 5.8 y. There are two shortlived isotopes that also occur naturally, 223 Ra and 224 Ra, which are alpha-gamma emitters from the 235 U and 232 Th decay chains with halflives of 11.4 d and 3.6 d, respectively. Uranium and thorium are widely distributed in nature, mainly in soils, sediments and rocks, and thus the four Ra isotopes are continuously produced in the environment at a rate that depends on the U and Th content and the half-life of each Ra isotope.
Whilst initially applied to study open ocean processes and to estimate different exchange processes across the continental margins, the most common application of Ra isotopes in the marine environment has undoubtedly been the quantification of water and solute fluxes from land to the coastal sea, driven by Submarine Groundwater Discharge (SGD) (Fig. 1) (Ma and Zhang, 2020). Radium isotopes have been used to quantify SGD in a wide range of marine environments, including sandy beaches (e.g., Bokuniewicz et al., 2015;Evans and Wilson, 2017;Rodellas et al., 2014), bays (e.g., Beck et al., 2008;Hwang et al., 2005;Lecher et al., 2016;Zhang et al., 2017), estuaries (e.g., Luek and Beck, 2014;Wang and Du, 2016;Young et al., 2008), coastal lagoons (e.g., Gattacceca et al., 2011;Rapaglia et al., 2010;Tamborski et al., 2019), salt marshes (e.g., Charette et al., 2003), large contintental shelves (e.g., Liu et al., 2014;Moore, 1996;Wang et al., 2014), ocean basins (e.g., Rodellas et al., 2015a) and the global ocean (e.g., Kwon et al., 2014). SGD is now recognized as an important process worldwide and plays a major role as a conveyor of dissolved compounds to the ocean (e.g., nutrients, metals, pollutants) (Cho et al., 2018;Moore, 2010;Rodellas et al., 2015a). Thus, SGD has significant implications for coastal biogeochemical cycles (e.g., Adolf et al., 2019;Andrisoa et al., 2019;Garces et al., 2011;Garcia-Orellana et al., 2016;Ruiz-González et al., 2021;Sugimoto et al., 2017). This review is focused on the application of Ra isotopes as tracers of SGD-derived inputs of water and solutes to the coastal ocean but also aims to provide some considerations that will help refine SGD estimates and identify critical knowledge gaps and research needs related to the current use of Ra isotopes as SGD tracers.  Radium was discovered from pitchblende ore in 1889 by Marie Sklodowksa-Curie and Pierre Curie, representing one of the first elements discovered by means of its radioactive properties (Porcelli et al., 2014). This historic discovery initiated a remarkable use of Ra for medical, industrial and scientific purposes, including its use as an ocean tracer. Early oceanic investigations revealed elevated 226 Ra activities in deep-sea sediments (Joly, 1908), but the first measurement of 226 Ra activities in seawater was performed in 1911 to characterize the penetrating radiation observed in the ionization chambers of ships at sea (Simpson et al., 1911). However, it was not until the publication of the first profile of Ra activities in seawater and sediments off the coast of California by Evans and Kip (1938), which suggested the influence of sediments on the increase of seawater 226 Ra activities with depth, that the use of Ra isotopes to study oceanographic processes was considered. The first true application of Ra as an oceanic tracer was by Koczy (1958), who showed that the primary oceanic source of 226 Ra was 230 Th decay in marine sediments, and who initiated the use of one-dimensional vertical diffusion models to describe the upward flux and characteristic concentration profiles of 226 Ra. Shortly thereafter, Koczy (1963) estimated that approximately 1-5% of the actual 226 Ra present in seawater was supplied by rivers. This was followed by a study of water column mixing and ventilation rates in the different oceans (Broecker et al., 1967). A global-scale picture of the distribution of long-lived Ra isotopes ( 226 Ra and 228 Ra) was first provided by the Geochemical Ocean Sections Study (GEOSECS -1972(GEOSECS - -1978e.g., Chan et al., 1977;Trier et al., 1972).

Early days: radium isotopes as tracers of marine processes
In the mid 1960s, Blanchard and Oakes (1965) reported elevated activities of 226 Ra in coastal waters relative to open ocean waters. The detection of significant activities of 228 Ra in both surface and bottom waters led to the understanding that the distribution of Ra depended on sediments but also on other sources that introduced Ra into the oceans (Kaufman et al., 1973;Moore, 1969aMoore, , 1969b. Few years later, Li et al. (1977) observed higher concentrations of 226 Ra in the Hudson River estuary compared with those observed either in the river itself or in the adjacent surface ocean water, proposing that 226 Ra was released by estuarine and continental shelf sediments, which represents an important source of 226 Ra to the ocean. This was further corroborated by measurements in the Pee Dee River-Winyah Bay estuary, in South Carolina, USA, by Elsinger and Moore (1980) who concluded that Ra desorption from sediments could quantitatively explain the increase of 226 Ra in brackish water. Several studies followed these initial investigations, highlighting the importance of coastal sediments in the release of long-lived Ra isotopes to coastal waters (e.g., Li and Chan, 1979;Moore, 1987).
Measurement of 226 Ra and 228 Ra in seawater samples was relatively complex prior to the wide-scale availability of high-purity germanium (HPGe) gamma detectors and the extraction of Ra from seawater using manganese-impregnated acrylic fibers. While 226 Ra was measured in a few liters of seawater (5-10 L) by the emanation of 222 Rn (Mathieu et al., 1988), the measurement of 228 Ra required greater volumes (often >> 100 L). The quantification of 228 Ra proceeded either via the extraction of 228 Ac and the measurement of its beta activity or via the partial ingrowth of 228 Th over a period of at least 4-12 months and the measurement via the alpha recoil of 228 Th or via its daughter 224 Ra by a Fig. 1. Record of published research articles that used the keywords 1) 'radium' and 'submarine groundwater discharge' (n = 511) and 2) 'radium' and 'ocean' (n = 477), as indexed by Web of Science. Articles related to NORM, radiation protection or sediment concentrations have been removed from the record of 'radium' and 'ocean'. proportional counter (Moore, 1981;Moore, 1969aMoore, , 1969b. The introduction of manganese-impregnated acrylic fibers for sampling Ra isotopes in the ocean was also a key advance in the use of Ra as a tracer in marine environments (Moore, 1976;Moore and Reid, 1973). Before the development of this approach, large volumes of seawater and long and laborious procedures were required to concentrate Ra isotopes via BaSO 4 precipitation (Kaufman et al., 1973;Moore, 1969aMoore, , 1969bMoore, 1972;Trier et al., 1972). Since then, Ra isotopes in marine waters are concentrated in situ with minimum time and effort by manganese-impregnated acrylic fibers. More recently, mass spectrometry has been used for 226 Ra determination in seawater, such as thermal ionization mass spectrometry (TIMS) (e.g., Ollivier et al., 2008), multicollector inductively coupled plasma mass spectrometry (MC-ICP-MS) (e.g., Te Hsieh and Henderson, 2011) or single-collector sector field (ICP-MS) (e.g., Vieira et al., 2021). The application of the short-lived 223 Ra and 224 Ra as tracers of marine processes lagged the long-lived Ra isotopes by several years due to analytical limitations and because their potential use in oceanic studies was limited due to their shorter half-lives (11.4 d and 3.6 d for 223 Ra and 224 Ra, respectively). Short-lived Ra isotopes were first applied to estuarine mixing studies, which have time-scales of 1-10 days. Moore (1983a, 1983b) published one of the first studies on the distribution of 224 Ra, 226 Ra and 228 Ra, which was focused on the mixing zone of the Pee Dee River-Winyah Bay and Delaware Bay estuaries (USA). They showed that the main source of Ra isotopes were desorption and diffusion from suspended and bottom sediments, which together contributed to the non-conservative increase of the three isotopes in the river-sea mixing zone. Shortly thereafter, Bollinger and Moore (1984) suggested that 224 Ra and 228 Ra fluxes in a salt marsh from the US eastern coast could be supported by bioirrigation and bioturbation (i.e., sediment reworking and pumping driven by benthic fauna, respectively). Levy and Moore (1985) concluded that the potential use of 224 Ra as a tracer required a much better understanding of its input functions in coastal zones. The authors classified Ra sources into two main groups: (1) primary sources that include desorption from estuaries and salt marsh particles and; (2) secondary sources such as dissolved 228 Th present in the water column, longshore currents that may transport 224 Ra from other areas and in situ production from 228 Th decay adsorbed on suspended particles or bottom sediments. (1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000) Whilst the concept of SGD had been introduced in several early hydrological-based studies (e.g., Bokuniewicz, 1992;Bokuniewicz, 1980;Capone and Bautista, 1985;Cooper, 1959;Freeze and Cherry, 1979;Glover, 1959;Johannes, 1980;Kohout, 1966;Toth, 1963), the role of groundwater as a major conveyor of Ra isotopes to the coastal ocean was not established until the benchmark papers of Burnett et al. (1990), Veeh et al. (1995) and Rama and Moore (1996). Burnett et al. (1990) suggested that the high 226 Ra activities in water from the Suwannee estuary (USA) and offshore were most likely supplied by submarine springs or seeps. Veeh et al. (1995) demonstrated that the "traditional" Ra sources to coastal waters (e.g., rivers, sediments) were insufficient to support the 226 Ra activities in waters from the Spencer Gulf (South Australia), suggesting an external source for the excess 226 Ra, such as the "submarine discharge of groundwater from granitic basement rocks". The importance of groundwater as a source of Ra to the coastal sea led Rama and Moore (1996) to apply for the first time the four naturally occurring Ra isotopes (which they defined as the radium quartet) as tracers for quantifying groundwater flows and water exchange (North Inlet salt marsh, South Carolina, USA). One key finding of this study was that extensive mixing between fresh groundwater and saline marsh porewater occurred in the subsurface marsh sediment before discharging to the coastal ocean, resulting in a groundwater discharge that was not entirely fresh but included a component of seawater circulating through the coastal aquifer. This concept was previously demonstrated in a hydrology-based study in Great South Bay (NY, USA) conducted by Bokuniewicz (1992), and it decisively contributed towards defining the term SGD.

Development period: radium isotopes as SGD tracers
Awareness of the volumetric and chemical importance of SGD was significantly increased by Moore (1996), who linked 226 Ra enrichments in coastal shelf waters of the South Atlantic Bight to large amounts of direct groundwater discharge. Moore suggested that the volume of SGD over this several hundred-kilometer coastline was comparable to the observed discharge from rivers, although the SGD probably included brackish and saline groundwater. In a commentary on this article, Younger (1996) questioned the inclusion of salty groundwater as a component of SGD, and its comparison to freshwater discharge to the coastal zone via rivers. The rebuttal by Church (1996) argued that SGD, regardless of its salinity, is chemically distinct from seawater and therefore can exert a significant control on coastal ocean biogeochemical cycles.
In the same year, Moore and Arnold (1996) published the design of the Radium Delayed Coincidence Counter (RaDeCC), an alpha detector system that allowed a simple and reliable determination of the shortlived Ra isotopes ( 223 Ra and 224 Ra) based on an original design of Giffin et al. (1963). Since its introduction, the RaDeCC system has become the reference technique for quantifying short-lived radium isotopes in water samples, as it allows a robust and rapid 223 Ra and 224 Ra determination with a simple setup at a relatively low cost. The full potential of short-lived Ra isotopes was described shortly after, when Moore published two studies providing conceptual models to derive ages of coastal waters and estimate coastal mixing rates using 223 Ra and 224 Ra (Moore, 2000a,b).
In the period of 1990-2000, these pioneering researchers opened a new era in Ra isotope application to SGD quantification, including the development of most of the conceptual models, approaches and techniques that are currently in use.

Expansion period: the widespread application of Ra isotopes as SGD tracers
The pioneering work on Ra isotopes of the 1990s, the increased understanding of the importance of SGD in coastal biogeochemical cycles, the technical improvements of Ra determination via HPGe and the commercialization of the first RaDeCC systems led to a rapid increase in the number of studies using Ra isotopes as tracers of SGD (Fig. 1). Hundreds of scientific articles have been published since then. During the 2000s, while some studies only used long-lived Ra isotopes as SGD tracers (e.g., Charette and Buesseler, 2004;Kim et al., 2005;Krest et al., 2000), most studies combined measurements of short-lived Ra isotopes to estimate water residence times (or mixing rates) and long-lived Ra isotopes to quantify SGD fluxes (e.g., Charette et al., 2003;Charette et al., 2001;Kelly and Moran, 2002;Moore, 2003). However, a few studies quantified SGD fluxes by using only short-lived Ra isotopes (Fig. 2;Boehm et al., 2004;Krest and Harvey, 2003;Paytan et al., 2004). One of the first studies where SGD fluxes were derived exclusively from short-lived Ra isotopes ( 223 Ra and 224 Ra) was conducted at Huntington Beach (CA, USA) by Boehm et al. (2004). Their flux estimates based on 223 Ra and 224 Ra were supported by a later analysis of 226 Ra at this location and by the results of a hydrological numerical model (Boehm et al., 2006). At the same period, Hwang et al. (2005) and Paytan et al. (2006) used only short-lived Ra isotopes to estimate SGD fluxes of water and associated nutrients. Over the last decade there has been a growing volume of studies using only the short-lived Ra isotopes as tracers of SGD (e.g., Baudron et al., 2015;Ferrarin et al., 2008;Garcia-Orellana et al., 2010;Krall et al., 2017;Shellenbarger et al., 2006;Tamborski et al., 2015;Trezzi et al., 2016). This relative abundance of publications, in which only short-lived Ra isotopes are used, is in contrast with those conducted in previous decades, in which long-lived Ra isotopes were the most applied tracers (Fig. 2).
The methodological advances in the measurement of short-lived Ra isotopes have also contributed to their widespread application as tracers of other coastal processes, such as water and solute transfer across the sediment-water interface (also referred to as porewater exchange (Cai et al., 2014), groundwater age as well as flow rates in coastal aquifers (Kiro et al., 2013;Kiro et al., 2012;Diego-Feliu et al., 2021), secondary permeability in animal burrows (Stieglitz et al., 2013), transfer of sediment-derived material inputs (Burt et al., 2013b;Sanial et al., 2015;Vieira et al., 2020;Vieira et al., 2019), and horizontal and vertical coastal water mixing rates Charette et al., 2007;Colbert and Hammond, 2007). Numerous new uses of the RaDeCC system have occurred in recent years, including the development of approaches to measure 223 Ra and 224 Ra in sediments (Cai et al., 2012;Cai et al., 2014) or 228 Ra and 226 Ra activities in water (Diego-Feliu et al., 2020;Geibert et al., 2013;Peterson et al., 2009;Waska et al., 2008), as well as quantification advancements and guidelines (Diego-Feliu et al., 2020;Moore, 2008;Selzer et al., 2021). There is still a widespread application of long-lived Ra isotopes, particularly to quantify SGD occurring over long flow paths (Moore, 1996;Rodellas et al., 2017;Tamborski et al., 2017a), SGD into entire ocean basins or the global ocean (Cho et al., 2018;Kwon et al., 2014;Rodellas et al., 2015a) and as a shelf flux gauge (Charette et al., 2016). Some investigations have also considered how to combine short-and long-lived Ra isotopes to discriminate between different SGD flow paths and sources (e.g., Charette et al., 2008;Moore, 2003;Rodellas et al., 2017;Tamborski et al., 2017a).

Submarine groundwater discharge: terminology
Initially, the concept of groundwater discharge into the oceans had been investigated largely by hydrologists, who considered the discharge of meteoric groundwater as a component of the hydrological cycle (Manheim, 1967). Thereafter, with the involvement of oceanographers interested in the influence of SGD on the chemistry of the ocean, the definition evolved, with SGD now encompassing both meteoric groundwater and circulated seawater (Bokuniewicz, 1992;Church, 1996;Moore, 1996;Rama and Moore, 1996). The interest in SGD has increasingly grown within the scientific community as it is now recognized as an important pathway for the transport of chemical compounds between land and ocean, which can strongly affect marine biogeochemical cycles at local, regional and global scales (e.g., Cho et al., 2018;Luijendijk et al., 2020;Rahman et al., 2019;Rodellas et al., 2015a;. SGD is also recognized to provide a wide range of ecosystem goods and services (Erostate et al., 2020;Alorda-Kleinglass et al., 2021) The definition of SGD has been discussed in several review papers in the marine geosciences field (e.g., Burnett et al., 2002;Burnett et al., 2001;Burnett et al., 2003;Church, 1996;Moore, 2010;Santos et al., 2012;Taniguchi et al., 2019;Taniguchi et al., 2002). Here, we will adopt the most inclusive definition, where SGD represents "the flow of water through continental and insular margins from the seabed to the coastal ocean, regardless of fluid composition or driving force" (Burnett et al., 2003;Taniguchi et al., 2019). We are thus including those centimeter-scale processess that can transport water and associated solutes across the sediment-water interface, which are commonly referred to as porewater exchange (PEX) or benthic fluxes, and are sometimes excluded from the SGD definition because of its short length scale (<1 m) (e.g., Moore, 2010;Santos et al., 2012). We are also implicitly considering that groundwater (defined as "any water in the ground") is synonymous with porewater (Burnett et al., 2003) and that the term "coastal aquifer" includes permeable marine sediments. The term "subterranean estuary" (STE), which is widely used in the SGD literature, is also used in this review to refer to the part of the coastal aquifer that dynamically interacts with the ocean Moore, 1999), determined in the hydrological literature as the freshsaline water interface or the mixing zone.
Within this broad definition of SGD we incorporate disparate water flow processes, some involving the discharge of fresh groundwater and others encompassing the circulation of seawater through the subterranean estuary, or a mixture of both. These processes can be grouped into five different SGD pathways according to the characteristics of the processes (George et al., 2020;Michael et al., 2011;Robinson et al., 2018;Santos et al., 2012) (Fig. 3): 1) Terrestrial groundwater discharge (usually fresh groundwater), driven by the hydraulic gradient between land and the sea; 2) Density-driven seawater circulation, caused by either density gradients along the freshwater-saltwater interface, or thermohaline gradients in permeable sediments; 3) Seasonal exchange of seawater, driven by the movement of the freshwater-saltwater interface due to temporal variations in aquifer recharge or sea level fluctuations; 4) Shoreface circulation of seawater, including intertidal circulation driven by tidal inundation (at beach faces, salt marshes or mangroves) and wave set-up; and 5) cmscale porewater exchange (PEX), driven by disparate mechanisms such as current-bedform interactions, bioirrigation, tidal and wave pumping, shear flow, ripple migration, etc. Notice that whereas all of these processes force water flow through the sediment-water interface, the discharge of terrestrial groundwater (Pathway 1) and, to a lesser extent, density-driven seawater circulation (Pathway 2), which also contains a fraction of freshwater, are the only mechanisms that represent a net source of water to the coastal ocean. Pathways 3-5 can be broadly classified as "saline SGD"; the summation of all five pathways is generally considered "total SGD".

Mechanisms controlling Ra in aquifers and SGD
One of the main characteristics that makes Ra isotopes useful tracers of SGD is that coastal groundwater is often greatly enriched with Ra isotopes relative to coastal seawater (Burnett et al., 2001). The enrichment of groundwater with Ra isotopes is originated from the interaction of groundwater with rocks, soils or minerals that comprise the geological matrix of the coastal aquifer. Some studies have summarized the various sources and sinks of Ra isotopes in coastal groundwater (Kiro et al., 2012;Krishnaswami et al., 1982;Luo et al., 2018;Porcelli et al., 2014;Tricca et al., 2001). The most common natural processes that regulate the activity of Ra isotopes in groundwater are: 1) radioactive Fig. 2. Number of research articles published in 5-year periods in which Ra isotopes have been used as tracers of SGD. The articles were classified according to the Ra isotopes used: only short-lived ( 223 Ra and/or 224 Ra), only long-lived ( 226 Ra and/or 228 Ra) or a combination of both. The search was performed using the words "SGD or Submarine Groundwater Discharge" AND "Ra or Radium" as keywords in "Web of Science" (n = 477). Reviews, method articles and studies not applying Ra isotopes as tracer of SGD-related processes were excluded from this classification. A total of 286 articles were considered in this classification analysis.
J. Garcia-Orellana et al. Before 2000200520102015 production from Th isotopes and decay; 2) adsorption and desorption from the aquifer solids and 3) weathering and precipitation. The transit time of groundwater in the coastal aquifer can also regulate the Ra activity in groundwater (Fig. 4).

Radioactive production and decay
Ra activity in groundwater is controlled, in part, by the production and decay rates of each radionuclide in both the groundwater and the geological matrix (Fig. 4). Radium isotope production depends on the continuous decay of Th isotopes ( 230 Th, 232 Th, 228 Th, and 227 Th) from the U and Th decay chains, which are predominantly contained in aquifer solids. Therefore, the U and Th content of the geological matrix regulates the production rate. Concentrations of U and Th can vary significantly depending on the host material (igneous rocks: 1.2-75 Bq⋅kg − 1 of U and 0.8-134 Bq⋅kg − 1 of Th; metamorphic rocks: 25-60 Bq⋅kg − 1 of U and 20-110 Bq⋅kg − 1 of Th; and sedimentary rocks: <1.2-40 Bq⋅kg − 1 of U and < 4 and 40 Bq⋅kg − 1 of Th) (Ivanovich and Harmon, 1992). For instance, carbonate minerals are enriched in U relative to Th, while clays are enriched in Th. As a consequence, groundwater flowing through formations with different lithologies have different activity ratios of Ra isotopes (e.g., 228 Ra/ 226 Ra AR, where AR means activity ratio), which can be used to identify and distinguish groundwater inflowing from different hydrogeologic units (Charette and Buesseler, 2004;Moore, 2006a;Swarzenski et al., 2007a) (see Section 6). On the other hand, radioactive decay depends only on the activity of Ra isotopes in groundwater and their specific decay constants.
Not all of the Ra produced in the aquifer is directly transferred to groundwater. The fraction of available Ra (i.e., the exchangeable Ra pool) includes the Th dissolved in groundwater (usually a negligible fraction) and mainly the Ra produced in the effective alpha recoil zone (i.e., by the Th bound in the outer mineral lattice and by the surfacebound, i.e., adsorbed Th) (Fig. 4). The decay of Th in the effective alpha recoil zone mobilizes part of the generated Ra from this zone to the adjacent pore fluid due to the alpha-decay recoil energy (Sun and Semkow, 1998;Swarzenski, 2007;Fig. 4). The extent of the effective alpha recoil zone is specific to each type of solid (i.e., mineralogy) and is a function of the size and of the surface characteristics of the aquifer solid grains Sun and Semkow, 1998;Swarzenski, 2007;Diego-Feliu et al., 2021). The greater the specific surface area of solids in the aquifer, the greater the fraction of Th in the effective recoil zone and thus, the pool of Ra available for solid-solution exchange (Copenhaver et al., 1993;Porcelli and Swarzenski, 2003). The influence of alpha recoil can produce deviations in the groundwater Ra isotopic ratios in relation to that expected from host rock ratios, since each Ra isotope is generated after a different number of decay events in each of the decay chains. For instance, when Ra isotopes in groundwaters are in equilibrium with aquifer solids, the alpha recoil process alone may produce equilibrium 226 Ra/ 228 Ra ratios up to 1.75 that of the host rock 238 U/ 232 Th ratio and 224 Ra/ 228 Ra equilibrium ratios ranging from 1 to 2.2 (e.g., Krishnaswami et al., 1982;Davidson and Dickson, 1986;Swarzenski, 2007;Diego-Feliu et al., 2021).

Adsorption and desorption
In fresh groundwater, Ra is usually and preferentially bound to aquifer solids, although a small fraction is found in solution as Ra 2+ . The adsorption of Ra on solid surfaces depends on the cation exchange capacity (CEC) of the aquifer solids, but also on the chemical composition of groundwater. The higher the CEC, the higher the content of Ra that may be potentially adsorbed on the grain surface of aquifer solids (e.g., Beneš et al., 1985;Kiro et al., 2012;Nathwani and Phillips, 1979;Vengosh et al., 2009). The ionic strength of the solution, governed mainly by the groundwater salinity, has been recognized as the most relevant factor controlling the exchange of Ra between solid and groundwater Gonneea et al., 2013;Kiro et al., 2012;Webster et al., 1995). High ionic strength (i.e., high salinity) hampers the adsorption of Ra 2+ due to competition with other cations dissolved in groundwater (e.g., Na + , K + , Ca 2+ ) and promotes the desorption of surface-bound Ra due to cationic exchange. As a consequence, Ra is typically more enriched in brackish to saline groundwater than in fresh groundwater. Exceptions may include carbonate aquifers where mineral dissolution is more important than desorption in driving groundwater Ra activities. Thus, Ra activities in coastal groundwater may vary substantially, depending on the subsurface salinity distribution, which is dynamic due to the interaction between the inland groundwater table elevation and marine driving forces (e.g., tides, waves and storms, mean sea-level). Other physico-chemical properties of groundwater, such as temperature and pH may also control the solidsolution partitioning (adsorption/desorption) of Ra in coastal aquifers.  reviewed the role of temperature and pH, concluding that while there is a Ra adsorption with increasing pH, there is no clear effect of temperature on the adsorption of Ra onto aquifer solids.
In order to estimate the relative distribution of Ra between solid and solution, a distribution coefficient, K D [m 3 ⋅kg − 1 ], is commonly used. This coefficient is defined as the ratio of the cocentration of Ra on the solid surface per mass of solid [Bq⋅kg − 1 ] to the amount of Ra remaining in mass of the solution at equilibrium [Bq⋅m − 3 ], and considers the Ra chemical equilibrium processes of adsortion-desorption (K D = Ra adsorbed /Ra desorbed ), but not the processes of weathering. The distribution coefficient is one of the most relevant parameters for understanding the Ra distribution in coastal aquifers and for applying transport models of Ra in groundwater. Radium distribution coefficient values span several orders of magnitude depending on the composition of groundwater and aquifer solids, ranging widely from 10 -2 to 10 2 m 3 ⋅kg -1 (Kumar et al., 2020;. Many authors analyzed the influence of different solid and water compositions within different experimental settings. The distribution coefficient is commonly determined by batch experiments (e.g., Colbert and Hammond, 2008;Gonneea et al., 2008;Rama and Moore, 1996;Tachi et al., 2001;Tamborski et al., 2019;Willett and Bond, 1995). The distribution coefficient of Ra has also been determined by other methods such as adsorption-desorption modelling (e.g., Copenhaver et al., 1993;Webster et al., 1995); chromatographic columns tests (e.g., Meier et al., 2015;Fig. 4. Conceptual model of the main processes determining the abundance of Ra isotopes in groundwater from a coastal aquifer. Blue and brown colour represents the liquid and solid phases, respectively. The light brown color represents the effective alpha recoil zone for Ra isotopes, which depends on the aquifer solids and the energy of the Ra alpha particle (the recoil distance usually ranges, on average, between 300 and 400 Å (Sun and Semkow, 1998)). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Relyea, 1982); and chemical equilibrium calculations (Puigdomenech and Bergstrom, 1995).

Weathering and precipitation
Radium isotopes are also supplied to groundwater by weathering processes, including dissolution and breakdown of rocks or minerals containing Ra that may occur during the flow of groundwater through the coastal aquifer. The intensity of weathering processes mainly depends on the physicochemical properties of groundwater, such as temperature, pH, redox potential or ionic strength, and on the characteristics of the geological matrix (e.g., mineralogy, specific surface area of aquifer solids) (Chabaux et al., 2003). Since the interaction between fresh and saline groundwater in the subterranean estuary is also often accompanied by a redox gradient (Charette and Sholkovitz, 2006;McAllister et al., 2015), the dissolution of hydrous oxides under reducing conditions may increase the activity of Ra isotopes in groundwater (Beneš et al., 1984;Gonneea et al., 2008). Conversely, mineral precipitation processes can remove Ra from groundwater. Due to the very low molar activities of Ra in groundwaters, the removal of Ra from groundwater usually occurs by co-precipitation with other phases. Ra 2+ may co-precipitate with Mn and Fe hydrous oxides under oxidizing and high pH (>7) conditions Porcelli et al., 2014), as well as with sulfates [e.g., (Ba, Ra, Sr)SO 4 ] or carbonate [(Ca, Ra) CO 3 ] (Kiro et al., 2013;Kiro et al., 2012;Porcelli et al., 2014). However, due to the relatively long time scales of weathering and precipitation processes, these vectors are often considered to exert negligible controls on the activities of 224 Ra, 223 Ra and 228 Ra in coastal groundwaters, although they might be relevant for 226 Ra (Porcelli and Swarzenski, 2003).

Groundwater transit time
The last fundamental factor that regulates the activities of Ra isotopes in aquifers is the time of groundwater to travel along a certain flow path, commonly known as groundwater transit time, which also represents the time that the groundwater is in contact with the solids in the aquifer Rajaomahefasoa et al., 2019;Vengosh et al., 2009). In SGD studies the groundwater transit time is instrumental to understand the degree of Ra isotope enrichment in groundwater since it entered the aquifer through one of its boundaries (e.g., aquifer recharge of freshwater or seawater infiltration through permeable sediments). Groundwater transit time largely depends on the hydraulic conductivity of the system and the physical mechanisms that are forcing the groundwater advective flow. The hydraulic conductivity describes the transmissive properties of a porous medium (for a given fluid) and thus clearly influences the groundwater flow velocity (Freeze and Cherry, 1979). For instance, hydraulic conductivity in unconsolidated silty-clay aquifers is typically on the order of 1 -100 cm⋅d − 1 , while in unconsolidated sandy aquifers it is significantly higher (10 2 -10 4 cm⋅d -1 ; Zhang and Schaap, 2019). An extreme example is represented by karstic or volcanic coastal aquifers, where localized hydraulic conductivity can reach values of ~10 7 cm⋅d − 1 (Li et al., 2020). On the other hand, physical mechanisms driving groundwater flow also control the groundwater velocity (as well as the spatial scale of the process), and therefore the groundwater transit times. In coastal aquifers, these driving forces include the terrestrial hydraulic gradient and their seasonal oscillations, shoreface circulation and tidal pumping, wave setup and wave pumping, bioirrigation, flow-and topography-induced advection, among others (Santos et al., 2012). These forces can be linked to the five SGD pathways outlined in Section 3. Obviously, groundwater transit times are not only determined by the driving mechanism itself, but also by the intensity and frequency of the physical forces (e.g., wave and tidal frequency and amplitude, magnitude and seasonality of the hydraulic gradient, recurrence and intensity of strong episodic wave events, etc.) Sawyer et al., 2013).
The common way to evaluate the effect of transit time on the degree of enrichment of Ra isotopes in groundwater is usually performed using 1-D advective-dispersive solute transport models. These models commonly assume: i) steady state aquifer conditions (i.e., activities do not vary with time; dA/dt = 0); ii) that the effects of dispersive and diffusive transport in relation to advective transport are negligible, and iii) that the exchange between the solid surfaces and the solution principally occurs via ion exchange, neglecting the processes of weathering or precipitation (Kiro et al., 2012;Michael et al., 2011;Tamborski et al., 2017a;Diego-Feliu et al., 2021). The underlying assumptions of these simple models may not be valid for all groundwater systems, although they allow understanding of Ra distribution within aquifers. Based on these models, the activities of Ra isotopes in groundwater increase along a flow path towards an equilibrium with Th activities in the aquifer solids (i.e., Th activities in the alpha recoil zone). This equilibrium, defined as bulk radioactive equilibrium by Diego-Feliu et al. (2021), is reached when the activities of Ra do not vary along the groundwater flowpath (dA/dτ = 0). The characteristic groundwater transit time needed to reach equilibrium between the isotopes of Ra and Th depends on the distribution coefficient of Ra (K D ) as well as on the decay constant of each Ra isotope (production and decay) Michael et al., 2011). The Ra loss due to advective transport when Ra is predominantly desorbed (low K D ; e.g., ~10 − 2 m 3 ⋅kg − 1 ) is higher than when Ra is mostly adsorbed onto grain surfaces (high K D ; e.g., ~10 2 m 3 ⋅kg − 1 ). Consequently, the characteristic groundwater time to reach the equilibrium with Th isotopes decreases as the distribution coefficient of Ra increases (Fig. 5). If Ra was completely desorbed (e.g., K D ~ 0 m 3 ⋅kg − 1 ), the time required to produce 50% of the equilibrium activity of a given Ra isotope would be equal to its half-life. However, in natural environments (K D ~ 10 − 3 to 10 3 m 3 ⋅kg − 1 ; , this activity would be reached within a significantly shorter time, ranging from minutes to hours for the short-lived Ra isotopes, from hours to days for 228 Ra and from days to years for 226 Ra (Fig. 5). Therefore, the enrichment rate of Ra isotopes strongly depends on the characteristics of the coastal aquifer (e.g., the higher the salinity, the longer the time to reach equilibrium concentrations). The relative difference between the characteristic transit time of each isotope is equivalent to the ratio of their decay constants. For example, the transit time needed for 228 Ra to Ko (m 3 kg-1 ) reach equilibrium with 232 Th is ~570 times higher (λ Ra− 228 /λ Ra− 224 ) than that for 224 Ra to reach the equilibrium with 228 Th. Characteristic transit times can be converted to equilibrium distances by assuming a constant velocity of groundwater through the aquifer. Considering a velocity on the order of 1 cm⋅d − 1 , all of the Ra isotopes would reach equilibrium within a few centimeters along the flow path in fresh groundwaters, while lengths would be on the order of tens of cm for 223 Ra and 224 Ra, ~1 m for 228 Ra and ~ 500 m for 226 Ra in saline groundwater. Notice that these lengths are specific for the assumed velocity and the distribution coefficients used. For instance, Tamborski et al. (2019) showed that short-lived Ra isotopes reached secular equilibrium after a flowpath of several meters in a sandy barrier beach from a hypersaline coastal lagoon, when groundwater velocities were on the order 10-100 cm⋅d − 1 . Given that physical mechanisms control the groundwater transit times in the subterranean estuary, the five pathways of SGD described in Section 3 are likely to be differently enriched in Ra isotopes, as they have different driving forces and spatio-temporal scales. The enrichment rate of Ra isotopes in each SGD pathway can be evaluated by comparing the ingrowth rates of different Ra isotopes with the most common spatiotemporal scales and groundwater salinities of the different pathways ( Fig. 6). Assuming a relatively homogeneous geological matrix, groundwater salinities are likely to control the K D of the different SGD pathways. As an orientative threshold, here we consider that a given Ra isotope becomes significantly enriched along a specific SGD pathway when 50% of its equilibrium activity has been reached ( Fig. 6). This 50% enrichment is an arbitrary (though reasonable) threshold because the ability of tracing Ra inputs from a given pathway depends on several site-specific factors, such as equilibrium activities, U/Th concentration in the coastal geological matrix, magnitude of SGD flow and interferences from other Ra sources (e.g., sediments, rivers). However, this qualitative comparison allows the identification of those SGD pathways that might be significantly enriched in the different Ra isotopes (Fig. 6). As illustrated in this assessment, short-lived Ra isotopes may become enriched in all of the SGD pathways. On the other hand, the transit time within the coastal aquifer of short-scale processes (e.g., cm-scale porewater exchange and shoreface seawater circulation) is not sufficient to produce measurable long-lived Ra isotopes activities (King, 2012;Michael et al., 2011;Moore, 2010;Rodellas et al., 2017).

Sources and sinks of Ra isotopes in the water column
The application of Ra isotopes as tracers of SGD in the marine environment (at local, regional and global scales) requires a comprehensive understanding of their sources and sinks in the water column. As summarized in several studies (e.g., Charette et al., 2008;Moore, 2010, Moore, 1999Rama and Moore, 1996;Swarzenski, 2007), the most common sources of Ra isotopes to the ocean can be classified into six groups: 1) Atmospheric input from wet or dry deposition (F atm ); 2) Discharge of surface water, such as rivers or streams (F river ); 3) Diffusive fluxes from underlying permeable sediments (F sed ); 4) Input from offshore waters (F in-ocean ); 5) Production of Ra in the water column from the decay of Th parents (F prod ) and 6) inputs from SGD (F SGD ). Major sinks of Ra in coastal systems are 1) Internal coastal cycling, which includes the biological or chemical removal of Ra isotopes through coprecipitation with minerals such as Ba sulfates or Fe (hydr)oxides and Ra uptake (F cycling ); 2) Radioactive decay of Ra (F decay ) and 3) Export of Ra offshore (F out-ocean ). A conceptual generalized box model summarizing all the pathways of removal or enrichment of Ra isotopes in the coastal environment is presented in Fig. 7.

Radium sources
The relative magnitude of the different Ra sources vary according to the characteristics of the study site (e.g., presence of rivers or streams, characteristics of sediments, water column depth), as well as the half-life of the specific Ra isotope used in the mass balance (Fig. 8). The different source terms for Ra isotopes into the coastal ocean and their relative importance are described below (excluding the Ra inputs from offshore waters (F in− ocean ), which will be considered in the net Ra exchange with offshore waters (see section 5.2)).
The input of Ra through atmospheric deposition onto oceans (F atm ) 10-a 10-10 10-4 10-2 10° 10 2 104 Ko (m3 kg-1) + Salinity includes dissolved Ra in precipitation and, mainly, desorption from atmospheric dust or ash entering the water column. The highest relative contribution of Ra from atmospheric sources will be expected to occur in large basins (i.e., the higher the ratio of ocean surface area to coastline length, the higher the potential importance of atmospheric inputs is) and/or in areas affected by large atmospheric inputs (e.g., areas affected by the deposition of large amounts of dust or ash). However, even in large areas affected by Saharan dust inputs (e.g., the Atlantic Ocean or the Mediterranean Sea), the percentage of the atmospheric contribution is less than 1% of the total Ra inputs Rodellas et al., 2015a) (Fig. 8). Atmospheric input is expected to be smaller in areas such as bays or coves, and therefore it is commonly neglected in Ra mass balances . Surface water inputs (F river ) includes both the fraction of Ra dissolved in surface waters and the desorption of Ra from river-borne particles as these particles encounter salty water. Surface water input includes both rivers and streams discharge and inputs from freshwater and salt marshes, stormwater runoff and anthropogenic sources, such as discharge from wastewater treatment facilities. There are several studies that have quantified the contribution of Ra from rivers and streams to coastal and open ocean areas (Krest et al., 1999;Moore, 1997;Moore and Shaw, 2008;Ollivier et al., 2008;Rapaglia et al., 2010). Most of the studies conclude that the fraction of dissolved Ra in river waters is a minor contribution compared with the Ra desorbed from surface waterborne particles (Krest et al., 1999;Moore et al., 1995;Moore and Shaw, 2008). The importance of surface water as a Ra source obviously depends on the presence (and significance) of rivers and streams in the investigated area. The relative contribution of surface water might be significant in estuaries or other areas influenced by river discharge (e.g., Beck et al., 2008;Key et al., 1985;Luek and Beck, 2014;Moore et al., 1995;Moore and Shaw, 2008) (Fig. 8). Marshes are also commonly considered as a major source of Ra due to both erosion of the marsh and subsequent desorption of Ra from particles and porewater exchange Bollinger and Moore, 1993;Bollinger and Moore, 1984;Charette, 2007;Tamborski et al., 2017c). Runoff or ephemeral streams following a major rain event (Moore et al., 1998) is an additional input. Anthropogenic activities like channels, oil and gas installations, outfalls and waste water treatment plants may also supply Ra to coastal waters, although in most cases they prove to be a minor Ra source (e.g., Beck et al., 2007;Eriksen et al., 2009;Rodellas et al., 2017;Tamborski et al., 2020aTamborski et al., , 2020b. Sediments are a ubiquitous source of Ra isotopes (F sed ), but their relative importance as a Ra source is highly dependent on the characteristics of the sediment (e.g., sediment grain size, sediment mineral composition, porosity), the specific Ra isotope considered, the mechanism that releases Ra from the sediments and the spatial scale of the study area (Fig. 8). Mechanisms that may release Ra from sediments, excluding groundwater flow, are molecular diffusion, erosion, bioturbation, or sediment resuspension. High Ra in the porewater results in a molecular diffusion of Ra from the sediments to the water column (e. g., Beck et al., 2008;Garcia-Orellana et al., 2014;Garcia-Solsona et al., 2008). Other processes such as sediment resuspension, diagenesis and bioturbation may also enhance the exchange of Ra between sediment and the water column (e.g., Burt et al., 2014;Garcia-Orellana et al., 2014;Moore, 2007;Rodellas et al., 2015b;Tamborski et al., 2017c). For instance, bioturbation has been suggested to increase the 228 Ra flux by a factor of two over the flux due to molecular diffusion only (Hancock et al., 2000). In sediments of the Yangtze estuary, bioirrigation was found to be more important for the 224 Ra flux from the sediments than molecular diffusion and sediment bioturbation (Cai et al., 2014). The significance of Ra input from sediments largely depends on the production time of each Ra isotope relative to the Ra-releasing mechanism. Inputs from sediments can typically be ignored for long-lived Ra isotopes in small-scale studies due to their long production times (e.g., Alorda-Kleinglass et al., 2019;Beck et al., 2008;Beck et al., 2007;Garcia-Solsona et al., 2008). However, in basin-wide or global-scale studies, the large seafloor area results in long-lived Ra sediment input that might be comparable to SGD input Rodellas et al., 2015a) (Fig. 8). The fast production of short-lived Ra isotopes in sediments, which is set by their decay constants, results in a near-  continuous availability of 223 Ra and 224 Ra in sediments, which can result in comparatively larger fluxes of short-lived Ra isotopes to the water column relative to the long-lived ones. Particularly in areas with coarsegrained sediments with low Ra availability for sediment-water exchange, fluxes of longer-lived Ra isotopes from seafloor sediments usually account for a minor fraction of the total Ra inputs , Beck et al., 2007Garcia-Solsona et al., 2008). However, sediments can turn into a major source of Ra isotopes to the water column in shallow water bodies, in systems covered by fine-grained sediments (substrate with a high specific surface area), sediments with a high content of U and Th-series radionuclides, and/or in areas affected by processes that favor the Ra exchange between sediments and overlying waters, such as bioturbation (Cai et al., 2014), sediment resuspension (Burt et al., 2014;Rodellas et al., 2015b) or seasonal hypoxia .
The production of Ra isotopes from their dissolved Th parents (F prod ) is commonly avoided in SGD studies by reporting Ra activities as "excess" activities in relation to their respective progenitors. This is usually a minor Ra source because Th is a particle reactive element and is rapidly scavenged by particles sinking through the water column. The amount of 232 Th, which is introduced into the water column by dissolution from particles supplied by rivers, runoff, atmospheric deposition or resuspension, is very low and therefore water column dissolved 232 Th activities are usually orders of magnitude lower than its daughter 228 Ra. Dissolved activities of 227 Th (and 227 Ac), 228 Th and 230 Th in the water column are also low and thus the production of 223 Ra, 224 Ra and 226 Ra from their decay generally represents a minor source term.
Submarine Groundwater Discharge (F SGD ) is often a primary source of Ra isotopes to the ocean (Fig. 8) and it is the target flux in SGD investigations. This term includes Ra inputs from any water flow across the sediment-water interface, which can be supplied through the five different SGD pathways outlined in Section 3 (Fig. 3). Different SGD pathways occur at different locations (e.g., nearshore, beachface, offshore) and have different groundwater compositions (e.g., fresh, brackish or saline groundwater) and characteristic groundwater transit times within the subterranean estuary (e.g., hours-days for cm-scale porewater exchange and months-decades for terrestrial groundwater discharge). All of them are, however, potential sources of Ra isotopes to the ocean and thus need to be taken into account when Ra isotopes are used as SGD tracer.

Radium sinks
The two main Ra sinks for most of the coastal systems are the decay of Ra isotopes in the water column (F decay ) and net Ra exchanged with offshore waters (F out− ocean -F in− ocean ). As with the Ra sources, the relevance of the Ra sink terms is largely dependent on the characteristics of the water body (mainly water residence time) and on the Ra isotope used (Fig. 9). Other sinks are grouped together as internal cycling (F cycling ). These processes include Ra co-precipitation with salts (e.g. barium sulphate) or Fe-Mn (hydr)oxides that occur in estuaries, coastal lagoons or polluted areas (e.g., Alorda-Kleinglass et al., 2019;Kronfeld et al., 1991;Neff and Sauer, 1995;Snavely, 1989), the scavenging of Ra with sinking particles Moore and Dymond, 1991;van Beek et al., 2007), uptake by biota, including the incorporation into calcium carbonate, barium sulphate or calcium phosphate lattice of shells and fish bones (Iyengar and Rao, 1990;Szabo, 1967), and the adsorption of Ra to the outer surface of algae's or to their internal non-living tissue components (Neff, 2002). These various internal cycling processes are generally a negligible Ra sink compared to radioactive decay and exchange with offshore waters. However, when using long-lived Ra isotopes and conducting basin-scale and global ocean budgets (i.e., low decay and low exchange), F cycling should be taken into account Rodellas et al., 2015a) (Fig. 9).
The decay of Ra isotopes (F decay ) is characteristic output term of mass balances using radioactive isotopes as tracers. The loss of Ra due to decay is usually a term that is relatively easy to constrain because it only depends on the Ra inventory of the study site and the decay constant of the isotope used. Thus, when decay is a primary sink of Ra in the system studied, Ra removal can be accurately contrained provided that the Ra water-column inventory has been adequately determined . The importance of the decay term will depend on the halflife of the isotope used and the water residence time of the system studied, as illustrated in Figs. 9 and 10. The loss of Ra due to radioactive decay is negligible for long-lived Ra isotopes in coastal areas with relatively short residence times (< 100 days) (Fig. 9). In regional and global-scale studies, the decay of 228 Ra needs to be considered and usually represents the main sink of this isotope (e.g., Charette et al., 2015;Kwon et al., 2014;Liu et al., 2018;Rodellas et al., 2015a) (Fig. 9). On the contrary, 226 Ra can be considered as a stable element in SGD studies because decay is negligible on the timescale of processes occurring in coastal or regional systems. For the short-lived Ra isotopes, the removal of Ra due to decay must be The relative contribution of the exchange of Ra isotopes due to the mixing between coastal and open waters (offshore exchange -the difference between F out-ocean and F in-ocean ) depends on the study site. Nearshore waters usually have higher Ra concentrations than offshore waters, and thus there is usually a net export of Ra offshore. The loss of Ra due to the mixing between coastal and open waters is directly linked to the flushing time of Ra in the coastal system. In open coastal systems or sites with relatively short flushing times (e.g., systems with short water residence times and/or high dispersive mixing with offshore waters), the export of Ra isotopes offshore is frequently the primary removal term and thus, it is commonly one of the most critical parameters to be determined in Ra mass balances, particularly for the longlived Ra isotopes (Fig. 10) (Tamborski et al., 2020a(Tamborski et al., , 2020b. In semienclosed or enclosed coastal environments (e.g., bays, coastal lagoons, coves), flushing times are usually long, and Ra losses due to mixing with offshore waters is less important for the short-lived Ra isotopes. However, for the long-lived Ra isotopes this process is the main sink even in these semi-enclosed water bodies (Fig. 9), requiring an appropriate characterization of the mixing term to obtain an accurate quantification of SGD. Different Ra isotopes are commonly combined to constrain offshore exchange (e.g., Moore, 2000aMoore, , 2000bMoore et al., 2006) (see section 8.3), although this output term can also be estimated using other approaches such as hydrodynamic and numerical models (e.g., Chen et al., 2003;Lin and Liu, 2019;Warner et al., 2010), tidal prism (e.g., Dyer, 1973;Petermann et al., 2018;Sheldon and Alber, 2008) or direct current/flow measurements (e.g., Rodellas et al., 2012;Shellenbarger et al., 2006).

Ra-based approaches to quantify SGD
Ra isotopes are suitable SGD tracers mainly because i) activities in groundwater are typically 1-2 orders of magnitude higher than in coastal seawater; ii) Ra isotopes in coastal areas are usually primarily sourced from SGD; iii) they behave conservatively in seawater and iv) they have different half-lives (ranging from 3.7 days to 1600 years), therefore allowing the tracing of coastal processes on a variety of timescales. The common approach to quantify SGD is to quantify first the Ra flux supplied by SGD (F SGD ; Bq d − 1 ), regardless of the SGD pathway considered, and subsequently convert it into a volumetric water flow (SGD; m 3 d − 1 ) by characterizing the Ra activity in the discharging groundwater (i.e., the SGD endmember, C Ra-SGD ; Bq m − 3 ) (Eq. 1).

SGD
There are three basic strategies to quantify total SGD fluxes (F SGD ) using Ra isotopes: i) mass balances; ii) endmember mixing models and iii) offshore flux determination from horizontal eddy diffusive mixing. The most comprehensive and widely applied approach is the Ra mass balance, where the flux of Ra supplied by SGD is usually quantified by a "flux by difference approach", which considers all the potential Ra sources and sinks identified in Fig. 7: where A is the average Ra activity in the study area, V is the volume affected by SGD and t is time (i.e., ∂AV ∂t is the change of Ra activity in the study area over time). This approach is often simplified in coastal areas by neglecting the commonly minor Ra sources and sinks (atmospheric inputs, production and internal cycling) and assuming that the system is in steady state (i.e., ∂AV ∂t = 0), and thus all quantifiable Ra input fluxes are subtracted from the total output with the residual being attributed to SGD: where A ocn is the Ra activity of the open ocean water that exchanges with the study area, t f is the flushing time of Ra in the system due to mixing and λ is the radioactive decay constant of the specific Ra isotope used. Notice that the first and second terms on the right side of the equation describe offshore exchange (i.e., F out-ocean -F in-ocean ) and radioactive decay, respectively. Flushing time is an integrative time parameter used to describe the solute (Ra isotopes) transport processes in a surface water body (due to both advection and dispersion) and it relates the mass of a tracer and its renewal rate due to mixing (Monsen et al., 2002). Notice also that we refer to flushing time of Ra due to mixing only, as if it were a conservative and stable solute. Some authors use the concept of water residence time (i.e., the time a water parcel remains in a waterbody before exiting through one of the boundaries; t r ) to refer to t f , but this term is not appropriate for systems influenced by dispersive mixing (such as many of coastal sites) or evaporation where tracers and water are transported due to different mechanisms. This approach is commonly applied in environments where the distribution of Ra isotopes in the study site is influenced by multiple Ra sources and sinks (Fig. 7). The accuracy of the estimated Ra flux supplied by SGD will thus largely depend on the accuracy with which the most relevant sources and sinks of Ra are determined . Indeed, given that SGD is quantified by a difference of fluxes, large uncertainties are often associated with the SGD estimates in systems where SGD fluxes only represent a minor fraction of total inputs. Most studies quantify total SGD with a mass balance by using a single Ra isotope and use other short-lived isotopes to estimate the flushing time (t f ) (see section 8.3) that is required in SGD estimations (Eq. 3) (e.g., Gu et al., 2012;Kim et al., 2008;Krall et al., 2017). Ra isotopes can also be combined through concurrent mass balances to simultaneously quantify SGD derived from different aquifers or pathways . The second approach estimates SGD fluxes using mixing models between different endmembers (e.g., Charette and Buesseler, 2004;Charette, 2007;Moore, 2003;Young et al., 2008). This approach is essentially a simplification of the mass balance approach where all Ra inputs are attributed to water flows (e.g., groundwater, rivers) and where the mixing offshore is assumed to represent the main Ra sink (Moore, 2003). This approach is thus recommended for environments where mixing between SGD and the ocean controls the distribution of Ra in the study site and where other Ra sources (e.g., sediment inputs) and sinks (e.g., decay) can be neglected. If various Ra isotopes are used, this method can also be used to distinguish the relative contribution of SGD and other water sources (e.g., rivers; Dulaiova et al., 2006aDulaiova et al., , 2006b or to separate different SGD components (e.g., SGD fluxes from confined and unconfined aquifers (Charette and Buesseler, 2004;Moore, 2003), terrestrial SGD and marsh porewater (Charette, 2007). For example, one can consider a system with two Ra sources (e.g. SGD1 and SGD2), aside from offshore exchange. SGD inputs can be estimated by defining the major contributors to the Ra isotope budgets of two isotopes (for example, 226 Ra and 228 Ra) and solve a series of simultaneous equations (Moore, 2003).
where 228 Ra and 226 Ra represent the measured Ra activity for a given time period, f is the water mass fraction contributed by coastal ocean mixing (f ocn ), the SGD source 1 (f SGD1 ) and the SGD source 2 (f SGD2 ), and the Ra isotope endmembers are identified by the same series of subscripts. With this approach, the f SGD can be determined for each location where there is a Ra isotope measurement ( 22X Ra). Short-lived Ra isotopes can also be used in the equations by incorporating their decay term. This approach is particularly useful to resolve changes in mixing between different Ra sources over tidal time-scales, e.g., within marsh creeks (Charette, 2007). The fractional SGD contributions derived from this mixing model can be converted into volumetric fluxes by considering the time scales of water mass transport (t r ) in the system under study, as follows: If the water outflow of the system can be directly measured (e.g., using an Acoustic Doppler Current Profiler -ADCP), then the endmember fraction can be directly multiplied by the measured flow (e.g., Rodellas et al., 2012;Tamborski et al., 2021).
The third approach is based on using offshore transects of short-lived Ra isotopes to estimate coastal mixing rates via the offshore coefficient of solute dispersivity, K h [m 2 ⋅s − 1 ], which can be used in conjunction with the offshore gradient of long-lived Ra isotopes to estimate the offshore export of 228 Ra or 226 Ra (Moore, 2015;Moore, 2000b). Several assumptions are required to apply this model: i) there is no additional input of Ra beyond the nearshore source, ii) the system is in steady state on the timescale of the isotope used (see section 7.3), iii) advection is negligible in any direction; iv) the open ocean Ra isotope activities are negligible and v) a flat seabed or the presence of a stratified water column of a constant thickness Moore, 2015;Moore, 2000b). Considering these assumptions, the measured log-linear decrease in the activity of 223 Ra or 224 Ra offshore can be used to determine K h from the following simplified advection-diffusion equation (Eq. 8).
where A 0 and A x are the Ra activities at the coast and at a distance x from the coast, respectively. A detailed discussion on the model and its assumptions is included in Moore (2015Moore ( , 2000b. The mixing coefficient K h derived from the logarithmic-linear scale plot of the short-lived Ra isotopes activities is then combined with the offshore gradient of longlived Ra isotopes (linear-linear decrease) to estimate the export offshore of 228 Ra or 226 Ra. This long-lived Ra flux offshore needs to be balanced by Ra inputs from all the potential sources, and thus estimating SGD from this approach requires subtracting the contribution of all the sources aside SGD from the estimated flux offshore (similar to the mass balance approach) (Moore, 2015, Moore, 2000b. This approach is recommended for open coastal systems where the export of long-lived Ra isotopes offshore due to diffusive mixing is frequently the primary removal term for 226 Ra and 228 Ra (e.g., Bejannin et al., 2020;Boehm et al., 2006;Dulaiova et al., 2006aDulaiova et al., , 2006bRodellas et al., 2014). One of the advantages of using this approach is that it allows combining shortlived Ra estimates of K h with solute gradients offshore to obtain the rate of solute transport offshore (e.g., Charette et al., 2007). For these reasons, this approach has been widely applied in oceanographic studies focused on the transport of solutes supplied by a combination of boundary exchange processes and is generally not applied to discriminate different Ra sources (e.g., SGD, rivers, sediments) (e.g., Charette et al., 2016;Jeandel, 2016;Vieira et al., 2019Vieira et al., , 2020.

Determination of the Ra concentration in the SGD endmember
The above-mentioned approaches of SGD quantification rely on an appropriate characterization of the Ra activity in the discharging groundwater (i.e., the SGD endmember) (Eq. 1). This characterization and selection of the Ra activity in the discharging groundwater often represents the main source of uncertainty for the final SGD estimates, since Ra activities in the coastal aquifer can vary significantly (1-3 orders of magnitude) over space and time depending on the geochemical and hydrogeological characteristics at each study site (e.g., Cerdà-Domènech et al., 2017;Cho and Kim, 2016;Cook et al., 2018;Duque et al., 2019;Gonneea et al., 2008Gonneea et al., , 2013Michael et al., 2011). Several studies highlighted the inherent difficulties and uncertainties related to the selection of representative SGD endmembers (Cerdà-Domènech et al., 2017;Cho and Kim, 2016;Cook et al., 2018;Gonneea et al., 2013;Michael et al., 2011). There is no general framework to characterize SGD endmembers valid for all study sites. The SGD endmember is commonly determined from near-shore piezometers or inland wells (e.g., Charette et al., 2013;Garcia-Solsona et al., 2010a;Tovar-Sánchez et al., 2014), porewaters collected at the seafloor, usually within the first meter below the sediment-water interface Tamborski et al., 2018), direct measurements of SGD actually discharging to the system through seepage meters or springs Montiel et al., 2018;Weinstein et al., 2007) or laboratory experiments to estimate the Ra activity in equilibrium with sediments Garcia-Solsona et al., 2008). Once groundwater samples are collected at different locations and/or times, in most published studies Ra values are averaged to obtain a 'representative average' SGD endmember (e.g., Beck et al., 2007;Kwon et al., 2014;Lee et al., 2012;Rapaglia et al., 2010). An important consideration when determining the SGD endmember is that the selected Ra concentration needs to be representative of the discharging groundwater. In coastal systems with a dominant pathway, the Ra endmember concentration should thus be representative of this pathway. Whilst this recommendation is self-evident, many studies have overlooked or they assumed that the sampled endmember is representative for the SGD pathways. Given the general and inherent difficulties of sampling SGD, a common strategy to obtain SGD endmembers is measuring Ra concentrations in existing wells or piezometers, which are generally installed in the freshwater (or brackish) zone of the coastal aquifer. However, these wells may not span the range of Ra activities along the salinity gradients in a subterranean estuary and thus they are not necessarily representative of the SGD endmember in sites dominated by saline SGD (Pathways 2 to 5) (Michael et al., 2011). Cho and Kim (2016) highlighted and discussed this issue at a global scale and showed that the endmembers frequently used to estimate total SGD are often low salinity groundwater samples (i.e., with relatively low Ra concentrations), even if SGD fluxes are composed mainly of circulatated seawater (Pathways 2 to 5). These authors re-evaluated previous SGD estimates for the Atlantic Ocean and for the global ocean by using only endmembers with salinities ~ > 10 and showed these fluxes were likely overestimated two-to three-fold only because non-representative lowsalinity groundwaters were included in the determination of the SGD endmember (Cho and Kim, 2016). On the other hand, if SGD mainly consists of seawater circulation with short transit times within the subterranean estuary (e.g., shoreface seawater circulation or cm-scale porewater exchange), the long-lived isotope activities (e.g., 228 Ra) of circulated seawater (Pathways 2 to 5) may be much lower than terrestrial water activities (even in fresh water of Pathway 1). This bias towards on-shore borehole/well data (long transit times) may actually result in an underestimation of SGD, in particular in basin-scale studies.
In more complex settings where several SGD pathways coexist, the SGD endmember should account for the different Ra enrichments in the SGD pathways that occur in the study area and should reflect the relative contributions of these different pathways. Obtaining the representative Ra endmember concentration thus requires, in addition to determining the Ra activity in each pathway, constraining the relative contribution of each SGD pathway. This is difficult and requires previous knowledge of the fluxes from different pathways, which is usually not available Michael et al., 2011). In addition, when seawater and groundwater mix, there are salinity gradients and redox interfaces, which promote geochemical reactions Moore, 1999;Rocha et al., 2021), that can further complicate differentiation between the different groundwater discharge processes. To simplify this approach, previous studies have argued that some of these pathways are negligible either because the driving forces are not significant (e.g., no tidal-driven recirculation in microtidal environments; Alorda-Kleinglass et al., 2019;Krall et al., 2017;Rodellas et al., 2017;Trezzi et al., 2016) or because low permeabilities of ocean sediments restrict inputs from some SGD pathways (e.g., Beck et al., 2007;Michael et al., 2011;Michael et al., 2005). Some authors also identified groundwater pathways, which were weakly enriched in the specific Ra isotope used, and then argued that the system was not significantly influenced by Ra inputs from this pathway (e.g., Beck et al., 2007;Garcia-Solsona et al., 2008;Moore et al., 2019). It should be additionally noted that different SGD pathways often mix before discharge and do not occur as isolated processes. When two pathways mix before discharge, this mixed SGD can be treated as a single pathway, provided that the Ra concentration in this mixed flow is appropriately characterized .

Improving the current use of Ra isotopes as SGD tracers
Ra isotopes have been widely applied to quantify SGD in all kind of locations worldwide and in a large range of physical environments (e.g., lagoons, wetlands, estuaries, bays or ocean basins). The estimates of SGD derived from Ra isotopes are based on conceptual assumptions that are often not properly justified or validated, and may involve large uncertainties that are usually not quantificable and could lead to inaccuracies or unrealistic of the SGD estimates. In this section, we discuss some key considerations with the aim of improving future applications of Ra isotopes as SGD tracers.

Validation of model assumptions and consideration of conceptual uncertainties
The magnitude of SGD can be quantified using a variety approaches, including field methods (e.g., seepage meters, natural tracers and piezometric level measurements) or computational approaches (e.g., water budget models and numerical solutions of the groundwater flow equation) (Taniguchi et al., 2019). Each of these approaches is prone to certain biases and is limited in both the spatial coverage and the SGD pathways it is capturing Zhou et al., 2018). Using different approaches can yield order-of-magnitude differences in SGD at individual locations and, thus, the selection of the method itself constitutes the first source of uncertainty in any SGD study Zhou et al., 2018).
Tracer approaches are based on conceptual models that represent an abstraction of a complex and dynamic natural system. Thus, a second additional source of uncertainty is linked to the conceptualization of the natural system, the so-called conceptual or structural uncertainty (Regan et al., 2002). In the case of Ra isotopes, there are three models (or concepts) commonly applied to quantify SGD (see section 6.1): "Mass balance model" and "mixing model", which are both based on Ra budgets using lumped parameter models, and a "offshore export model", which is based on a simplified analytical solution to describe Ra transport. Each one of these models is extremely sensitive to its specific simplifications, assumptions and boundary conditions used to approximate the real system. For instance, Ra "mass balances" are strongly influenced by the sources and sinks of Ra included in the model (e.g., some studies assume that all Ra inputs can be directly attributed to SGD) and, mainly, by the parameterization of the different input and output terms and their intrinsic assumptions. Indeed, the use of different methods or parameterizations to quantify the same term in the mass balance (e.g., different approaches to estimate offshore exchange rates, sediment diffusion or to produce an averaged Ra concentration in the study area) may lead to order-of-magnitude discrepancies in Ra-derived SGD estimates . The mass balance approach also assumes that the integrative parameters can appropriately represent the system (e.g., the model assumes that offshore exchange occurs with an average Ra concentration). Similarly, the "offshore export" model is essentially a simplification of the 1-D advection-dispersion equation, which is very sensitive to the boundary conditions and basic assumptions. These assumptions may not be correct in environments with along-shore and across-shore currents (advection cannot be neglected), with Ra inputs occurring beyond the nearshore source, with not uniform vertical profiles of Ra isotopes or non-negligible Ra concentrations offshore, or in systems with a sloped seabed (e.g., Hancock et al., 2006;Webster, 2019a, 2019b). Additionally, all the Ra models are often based on the assumption that mixing processes affect all the isotopes alike, either through the estimation of flushing times or water ages (t f or t w ) using pairs of isotopes (see section 8.3) or by using an offshore diffusivity coefficient (K h ) derived from one isotope to evaluate the transport of another one (Lamontagne and Webster, 2019b). Since Ra transport offshore is often controlled by dispersion, which integrates processes operating at different temporal and spatial scales, mixing parameters (e.g., t f or K h ) are scale dependent and should be independently evaluated for each tracer (Lamontagne and Webster, 2019b;Moore, 2015;Okubo, 1976). Finally, and most importantly, almost all studies assume that the system is in steady state and that all model parameters are constant on the timescale the tracer resides in the system, assumptions that are very often not metin particular for longlived isotopesand that should be carefully evaluated (see section 7.2).
The reliability of Ra approaches to quantify SGD may thus be limited and scientists should be aware of these limitations and try to validate their assumptions and constrain the uncertainties of the estimated flows. Overcoming these limitations is not straightforward and there is not a general framework for accounting for these conceptual uncertainties or validating assumptions. Some authors have suggested evaluating the most sensitive parameters using multiple approaches and obtaining final SGD fluxes as an ensemble of results (e.g., Rodellas et al., 2021), combining multiple independent methods to constrain the magnitude of SGD (e.g., Zhou et al., 2018) or developing more complex transport models (e.g., Lamontagne and Webster, 2019b). For example, combining 'terrestrial' hydrogeological, oceanographic or geophysical approaches with tracers tools may help to improve SGD estimates. Producing more robust estimates of SGD is clearly a major open research topic in SGD and Ra investigations (see section 9.2).

The assumption of steady-state needs to be validated
A key general assumption in all the Ra-based approaches to quantify SGD is that all Ra fluxes (e.g., SGD-driven Ra flux, Ra diffusion from sediments, radioactive decay, Ra export offshore) are constant with respect to the timescale of the Ra isotope used and that the system is in steady state (i.e., Ra inputs equal Ra outputs). Therefore, any appropriate SGD study should validate these assumptions considering the time the tracer resides in the system, i.e. the tracer residence time . The average residence time of Ra isotopes in the surface water system depends on the rate of removal of the isotope and thus on their different sinks. Given that the residence time of Ra isotopes (t Ra ) in most of the systems is mainly controlled by radioactive decay and export offshore (see Fig. 9), it can be estimated following Eq. 9 : where λ is the radioactive decay constant of the specific Ra isotope used and t f is the flushing time in the system due to mixing. This flushing time is equivalent to mean water residence time in systems where the influence of evaporation and dispersion are negligible. Characteristic Ra residence time for study sites with different mixing timescales (i.e., different flushing times) are shown in Fig. 11. In systems where the timescales of transport processes exceeds that of radioactive decay (1/t f >> λ), the Ra residence time in the system is comparable to the temporal scale of coastal mixing. Thus, the assumptions of constant fluxes and steady state need to be validated on timescales of transport processes (Moore, 2015;Rodellas et al., 2021). For systems with negligible mixing processes relative to the isotope half-life (1/t f << λ), the Ra temporal scale is controlled by the radioactive decay of each isotope. In this latter case, the assumption of steady state and constant fluxes need to be validated on a timescale comparable to the mean life of the specific isotope used. Validating this assumption might thus be particularly challenging for long-lived Ra isotopes in systems with long flushing times. In those dynamic systems that present a relevant temporal variability on the timescale of the tracer used (e.g., systems with variable surface water inputs, temporal changes of SGD or variable mixing offshore), the evolution of Ra concentration in sources, sinks and inventories in the water column need to be understood. Neglecting the transience of these parameters could lead to significant errors on the final SGD estimates. In such non-steady state systems, a dynamic modelling approach might provide a better representation of the studied system (Gilfedder et al., 2015).

Ra isotopes cannot be used indistinctly to trace SGD
Ra isotopes have been commonly used to estimate total SGD, often without consideration as to whether the Ra isotope used is actually capturing the entire discharge pathway(s) underlying these flows. This raises the question of whether or not the four Ra isotopes provide the same SGD estimates and if they therefore can be used indistinctly.
In an idealized case where there is a unique SGD pathway or source (e.g., a karstic or volcanic area with a main coastal spring that dominates over any diffuse, non-point source discharge), all Ra isotopes should yield the same SGD flux, provided that i) all the Ra sources (aside from SGD) and sinks are appropriately constrained, ii) groundwater is sufficiently enriched with the Ra isotope used as a tracer, and iii) the SGD endmember is properly characterized. However, most natural systems do not satisfy these idealized conditions because of the ubiquitous presence of multiple driving forces (e.g., hydraulic gradient, wave and tidal pumping, bioirrigation) that result in concurrent inputs of Ra from different SGD pathways (Fig. 3). As illustrated in Fig. 6, the different SGD pathways are likely to have different groundwater composition and characteristic groundwater transit times, resulting in isotope-specific Ra enrichments, which may differ by orders of magnitude among different SGD pathways (  Ra isotopes are thus incorporating different pathways in a characteristic proportion, and thus they are likely to produce different SGD estimates. For instance, whilst SGD estimates based on 223 Ra and 224 Ra are likely to incorporate inputs from all the pathways, long-lived Ra isotopes may not incorporate SGD induced by short-scale processes (cm-scale circulation or shoreface saline circulation) (King, 2012;Michael et al., 2011;Moore, 2010;Rodellas et al., 2017) (Fig. 6). Therefore, when SGD is discharging via multiple pathways, all four Ra isotopes cannot be used indistinctly to obtain total SGD. Most Ra-based SGD studies aim at obtaining a single estimate of the total magnitude of SGD. In these cases, SGD is often treated as a single source, regardless of the difference in Ra activities of the multiple pathways, and only a single endmember is determined to represent the 'average' Ra activity in groundwater discharging into the study site (i.e., C Ra-SGD in Eq. 1). However, this Ra endmember should not be obtained from simple averages of all the SGD samples collected, because of possible biases towards the most sampled pathway. Michael et al. (2011) illustrated this issue by comparing SGD estimates obtained by sampling different pathways and showed that this could produce differences in final SGD estimates by more than an order of magnitude. Therefore, unless there is a comprehensive understanding of the system that allows constraining the relative contribution of different pathways, SGD estimates derived from a single Ra isotope might not be appropriate at sites with multiple pathways.
The combination of multiple Ra isotopes can be instrumental in those systems where several pathways are hypothesized to significantly contribute to Ra budgets (see also section 7.5). Both total SGD and the relative contribution of different SGD pathways can be accurately obtained from the concurrent application of multiple Ra isotopes. The number of Ra isotopes to be used in the models will depend on the number of SGD pathways to be determined, as well as the amount of other potential unknowns (e.g., residence time, river contribution). This approach can be applied through concurrent mass balances (see section 6.1; e.g., Rodellas et al., 2017;Tamborski et al., 2017a) or mixing models (e.g., Charette, 2007), provided that the Ra isotopes used have distinct relative enrichments (i.e., distinct Ra isotopic ratios) in the different pathways. Ra isotopes can also be combined with other tracers (e.g., silica, radon, salinity, stable isotopes, etc.) (e.g., Burnett et al., 2006;Garcia-Orellana et al., 2010;Oehler et al., 2019;Schubert et al., 2015) or approaches (e.g., seepage meters, Darcy estimates, heat tracing, tidal prism, hydrogeological models) (e.g., Garcia-Solsona et al., 2010a,b;Povinec et al., 2012;Prieto and Destouni, 2005;Rosenberry et al., 2020), which can decisively contribute towards the quantification of the different pathways and the appropriate characterization of SGD.

The Ra isotope(s) used need to be carefully selected
The selection of the tracer used should ideally be based on the process to be studied and the characteristics of the study site, rather than on the methods available (e.g., salinity or radon are much better SGD tracers than radium in certain systems). Likewise, the best suited Ra isotope for any specific study also needs to be chosen based on the target SGD pathway of interest, the area studied and according to the sensitivity of the final estimates of Ra inputs and outputs, which is specific for each coastal water system (Tamborski et al., 2020a(Tamborski et al., , 2020b. In principle, the time scale of the half life of the chosen Ra isotope should scale with the size and time escales of the system, such that the short-lived Ra isotopes are most effective at nearshore and embayment scales (e.g., Boehm et al., 2004;Knee et al., 2010;Rapaglia et al., 2012), while longlived Ra isotopes are most effective at regional and global-scales (water residence times >100 days; e.g., Charette et al., 2015;Kwon et al., 2014;Liu et al., 2018;Rodellas et al., 2015a). In either scenario, radioactive decay can be easily estimated from the Ra inventory in the water column and this facilitates characterizing total output fluxes in systems where radioactive decay is the primary Ra sink. On the contrary, mixing losses of Ra isotopes can be highly uncertain, including both water exchange at the boundaries of the system and the Ra endmember concentration (Tamborski et al., 2020a(Tamborski et al., , 2020b. Given that long-lived Ra isotopes are highly sensitive to mixing in coastal areas (i.e., estuaries, coastal lagoons), they may not be adequate tracers of SGD in coastal environments when mixing is uncertain (Ku and Luo, 2008;Rutgers van der Loeff et al., 2018). Short-lived Ra isotopes are ideal for coastal areas where timescales are on the order of days; however, sediment fluxes must be adequately characterized because they often represent a major source of short-lived Ra isotopes (Burt et al., 2014;Rodellas et al., 2015b). Water column inventories of short-lived Ra isotopes integrate over time-scales similar to the Ra isotope halflife, and therefore multiple samplings are necessary if seasonal or annual variations in SGD want to be captured.

Ra isotopes can contribute identifying the pathways of SGD
A proper evaluation of SGD requires the understanding of the dominant SGD pathways. There are several hydrogeological, geophysical and geochemical characterization techniques (e.g., Folch et al., 2020;Gonneea et al., 2008;Swarzenski et al., 2007b;Swarzenski et al., 2006;Zarroca et al., 2014), as well as numerical modelling approaches (e.g., Amir et al., 2013;Anwar et al., 2014;Danielescu et al., 2009), that provide fundamental information about the origin and spatio-temporal scales of the different SGD pathways through the subterranean estuary (Fig. 6). The characterization of Ra isotopes and especially their activity ratios (e.g., 228 Ra/ 226 Ra and 224 Ra/ 228 Ra) in the subterranean estuary can also provide very useful information to constrain the presence of different SGD pathways in the study site. For instance, since 228 Ra and 226 Ra belong to different decay chains (U and Th), the activity ratio of 228 Ra/ 226 Ra may provide information about the Th/U ratio in the host material and thus be used to identify and distinguish groundwater inflowing from different geological systems (e.g., Charette and Buesseler, 2004;Moore, 2006aMoore, , 2006bSwarzenski et al., 2007aSwarzenski et al., , 2007b. Additionally, the 224 Ra -228 Ra parent-daughter relationship enables the use of the 224 Ra/ 228 Ra activity ratio for characterizing different SGD pathways. Whilst short-time-scale SGD pathways are likely to be more enriched in short-lived Ra isotopes relative to long-lived ones (Fig. 6) presenting thus high 224 Ra/ 228 Ra activity ratios (> 1), long-scale pathways are likely to present similar 224 Ra and 228 Ra activities (i.e., 224 Ra/ 228 Ra ~ 1) . Moreover, the activity ratios of Ra isotopes in coastal seawater can also be used to evaluate the predominant discharge pathways of SGD, provided that Ra is mainly supplied to the coastal ocean by SGD (i.e., absence of any significant source such as diffusion from sediments, seawater or rivers). In this case and if radioactive decay is accounted for, the activity ratios of coastal seawater must reflect the combined contribution of the different SGD pathways, being biased towards the more relevant ones (Boehm et al., 2006;Moore, 2006b;Rodellas et al., 2014). For instance, at a site dominated by short-scale pathways (e.g., porewater exchange), the 224 Ra/ 228 Ra activity ratio of coastal sea water is expected to be comparable to that of porewater (i.e., 224 Ra/ 228 Ra > 1).

Estimates of SGD-driven solute fluxes need to account for the variable composition of different pathways
The most commonly applied method to estimate fluxes of dissolved chemicals (e.g., nutrients, metals, contaminants) driven by SGD is by multiplying either the (Ra-derived) SGD water flow by the average concentration of the chemical solutes in the SGD endmember or the SGD-driven Ra flux by the average solute/Ra ratio in the SGD endmember (Charette et al., 2016;Hwang et al., 2016;Santos et al., 2011;Tovar-Sánchez et al., 2014). This straight-forward approach might be appropriate for systems with one dominant SGD pathway (provided that the endmember concentration is appropriately determined), but it should not be directly applied in systems with multiple SGD pathways or sources. This is mainly because the biogeochemical composition of discharging fluids from different SGD pathways may be considerably different depending on the origin of solutes, their transformations within the subterranean estuary and the temporal and spatial scales of the different pathways (e.g., Rodellas et al., 2018;Santos et al., 2012;Tamborski et al., 2017a;Weinstein et al., 2011). Therefore, as discussed for Ra isotopes, the solute concentration obtained from averaging groundwater samples (and thus the flux estimates) is likely to be biased towards the most sampled pathway, and is not representative of the relative SGD flux from that pathway. For instance, terrestrial groundwater discharge (Pathway 1) often contains high concentrations of nutrients and other compounds released from anthropogenic activities, whereas other SGD pathways may not be influenced by anthropogenic sources. Using only solute concentrations from wells and boreholes located within the low-salinity part of the subterranean estuary may result in average solute concentrations having a strong anthropogenic signal (e.g., high nutrient concentrations) that might not be representative of SGD in sites influenced by pathways other than terrestrial groundwater discharge (see section 6.2). An appropriate understanding of the magnitude of solute fluxes driven by SGD, requires thus a previous identification of the dominant SGD pathways, which needs to be based on a detailed knowledge of the studied system. Sampling efforts should then be taken to target all the relevant SGD pathways for the particular system, particularly over the full-salinity gradient of the subterranean estuary. By doing so, the conservative and/or non-conservative behavior of the solute of interest can be determined, since some solutes can be chemically modified in the subterranean estuary by varying the ratio of solute to Ra along the SGD flow path. This is an essential step towards a correct quantification of SGD-driven solute fluxes in sites with multiple pathways and thus, towards a comprehensive understanding of the role of SGD for coastal biogeochemical cycles.

Additional applications of Ra isotopes in groundwater and marine studies
Aside from the application of Ra isotopes to quantify SGD, Ra isotopes can also provide instrumental information that can help to constrain the different terms needed for SGD evaluations, as well as key information for characterizing hydrological and oceanic systems. In this section, we briefly summarize parallel applications of Ra isotopes in groundwater and marine studies, which include: i) assessing transit times in coastal aquifers; ii) estimating solute fluxes across the sediment-water interface; iii) estimating ages of coastal surface waters; and iv) quantifying shelf-scale solute fluxes.

Assessment of groundwater transit times in coastal aquifers
Constraining groundwater transit times in a coastal aquifer or in the subterranean estuary is a key question in characterizing coastal hydrogeology. Applications include deriving characteristic distances to the upgradient recharge points, evaluating the connectivity of the coastal aquifer with the sea, determining ages and velocities of groundwater, differentiating variable spatio-temporal scale processes and evaluating the potential transformation of solutes in the subterranean estuary Gonneea and Charette, 2014;Lerner and Harris, 2009;Werner et al., 2013). Since observational techniques are limited in the subsurface, hydrogeological tracers (e.g., H and O stable isotopes, 3 He, 36 Cl, CFCs) are commonly used to obtain key information on groundwater flow (Leibundgut and Seibert, 2011). In recent years, the activities of Ra isotopes have been used to evaluate groundwater dynamics in the subterranean estuary Kiro et al., 2013;Tamborski et al., 2019;Tamborski et al., 2017b), which can be a complementary method to other approaches based on atmosphericintroduced tracers, which mainly decay while traveling through an aquifer (Cook and Herczeg, 2012). To evaluate the distribution of Ra in aquifers, each of the inputs and outputs of Ra to groundwater needs to be constrained. These input and output terms, which are summarized in Section 5, include production, decay, desorption, adsorption, weathering and precipitation, as well as dispersion. The most common analysis of Ra distribution in groundwater often utilizes a one-dimensional transport model (Eq. 10),

∂A Ra ∂t
where A Ra [Bq m − 3 ] is the activity of Ra in groundwater, t [s] is time, λ [s − 1 ] is the decay constant, and R Ra [− ] is the linear retardation factor of Ra. Retardation is defined as the retention of a solute due to interaction with solid phases relative to bulk solution in a dynamic system (McKinley and Russell Alexander, 1993). The first term on the right hand-side of the eq. (10) Kiro et al., 2013;Michael et al., 2011;Krest and Harvey, 2003;Schmidt et al., 2011;Tamborski et al., 2017b).
where A Ra, 0 is the activity of groundwater entering a streamline. This simple solution implicitly assumes that production rates are constant and thus it neglects the effect of spatial and/or temporal variations in the parent isotope activities (e.g., geological heterogeneity). Moreover, the determination of groundwater ages, velocities or transit times using this equation strongly relies on an accurate determination of the retardation factor of Ra (R Ra ), which may also be spatially and temporally variable. The retardation factor of Ra is commonly determined by using the distribution coefficient of Ra (Eq. 12; Michael et al., 2011).
where K D [m 3 kg − 1 ] is the distribution coefficient, ρ s [kg m − 3 ] is the dry bulk mass density, and ϕ is the porosity. If the production rate (P), the retardation factor of Ra (R Ra ) and the inflow boundary activities of Ra (A Ra, 0 ) are known, groundwater transit times can be easily derived from eq. 11 using a single Ra isotope. However, determinations of P, R Ra and A Ra,0 are usually difficult and induce large uncertainties that must be propagated in the calculation of groundwater transit times. In practice, this may be difficult due to spatial and temporal variability that affects the four Ra isotopes differently. Although the processes are integrated along a single flowpath, short-lived isotopes have a much shorter 'memory' than long-lived isotopes, so short-lived Ra activities are more sensitive to production rates and retardation factors more proximal to the sampling point. Alternatively, when the model parameters are unknown, groundwater transit times may be determined by using multiple Ra isotopes or activity ratios (e.g., 224 Ra/ 228 Ra) Tamborski et al., 2019). When the Ra production rate changes with depth, a depth-dependent 1D reactive transport model can be applied . Non-steady state conditions like variable hydrological forcing, wave set-up as well as tidal fluctuations may cause non-steady state conditions which may obscure groundwater transit times. Since this model assumes that there is no mixing between two different types of groundwater (with different ages and Ra content), it can cause problems in the interpretation of the transit times when mixing does occur, largely because mixing is a linear process whereas radioactive decay and ingrowth are exponential (Bethke and Johnson, 2002). Since the determination of groundwater transit times requires some assumptions in solving the transport equation (eq. 11) and involved terms, independent verifications of Raderived transit times are strongly recommended.

Solute flux across the sediment-water interface derived from Ra isotopes
The transfer of solutes across the sediment-water interface, through both diffusion and advection, can have a major effect on the chemical composition of the overlying waters and shallow sediments. Traditional approaches to quantify water exchange and solute fluxes using Ra isotopes rely upon one-dimensional reactive-transport models to reproduce porewater activities as a function of sediment production, bioturbation, molecular diffusion, dispersion and advection (Cochran and Krishnaswami, 1980;Krest and Harvey, 2003;Sun and Torgersen, 2001). Recent analytical developments have allowed the measurement of surfaceexchangeable 224 Ra and particle-bound 228 Th from sediment cores (Cai et al., 2012). Radioactive disequilibrium between the soluble 224 Ra daughter and its particle-bound parent 228 Th (T 1/2 = 1.9 y) may be produced in sediments and their interstitial pore fluids from the transport of 224 Ra as a result of molecular diffusion, bioturbation and irrigation and/or SGD. Assuming that the in-situ decay of 228 Th is the sole source of sediment 224 Ra, a steady-state mass balance can be written (Cai et al., 2014;Cai et al., 2012) as Eq. 13.
where F Ra is the flux of 224 Ra across the sediment-water interface (Bq⋅m − 2 d − 1 ), z is the depth in the sediment where disequilibrium occurs (m), λ Ra is the 224 Ra decay constant (0.189 d − 1 ) and A Th,B and A Ra,B are respectively the activities of bulk 228 Th and 224 Ra (Bq⋅cm − 3 ) of wet sediment, which is obtained from a mass unit considering sediment dry bulk density (ρ s = 2.65 g⋅cm − 3 ) and porosity (ø). A key advantage of the 224 Ra/ 228 Th disequilibrium method over traditional reactive-transport models (Boudreau, 1997a;Meysman et al., 2005) is that 224 Ra production is directly measured (via 228 Th) for each sediment section, and thus this approach accounts for depth-varying production rates. Derivation of a benthic 224 Ra flux from fine-grained (muddy) sediments can in turn be used to quantify the transfer rate (F i ) of a dissolved species across the sediment-water interface (Cai et al., 2014) using Eq. 14.
where D S i and D S Ra are the in-situ molecular diffusion coefficients for dissolved species i and 224 Ra, respectively, corrected for tortuosity and in-situ temperature (Boudreau, 1997b). The concentration gradients of dissolved species i (∂A i /∂z) and 224 Ra (∂A Ra /∂z) are taken as the concentration difference between overlying surface waters and shallow porewaters, typically measured at 1 cm depth, and are thus net solute fluxes (Cai et al., 2012(Cai et al., , 2014. Note that the flux obtained in the Eq. 14 is translated from an adjusted Fick's First Law (Eq. 15).
where ξ is an area enhancement factor representing the extended subsurface interface and advective influences on the diffusive flux, and Ф is porosity. Although Eq. 14 only includes molecular diffusive coefficients of dissolved species i and 224 Ra, it represents the sum of all processes that affect solute transfer across the sediment-water interface. The inherent assumption is that molecular diffusion is the ratelimiting step for solute transport. This approach is suitable to quantify fluxes across the sediment-water interface from shallow (< 20 cm) sediment cores in benthic and coastal (i.e., salt marsh) mud environments for oxygen (Cai et al., 2014;Dias et al., 2016), dissolved inorganic carbon and nutrients (Cai et al., 2015), rare earth and trace elements (Hong et al., 2018;Shi et al., 2018Shi et al., , 2019b. Note that solute fluxes derived from Eq. 14 integrate over a time-scale of ~10 days and may greatly exceed fluxes determined from porewater gradients (Cai et al., 2014(Cai et al., , 2015 and traditional sediment incubations, which alter the physical conditions of the sediments, and therefore may not capture small-scale (mm to cm) advective processes (Hong et al., 2018;Shi et al., 2019b). More complicated diagenetic models of 224 Ra have been developed for systems subject to significant bioturbation and particle reworking, and may thus require parallel sediment measurements of excess 234 Th to accurately constrain 224 Ra fluxes (Cai et al., 2014(Cai et al., , 2015. Concurrent sediment and water column 224 Ra mass balances may be used to separate short-scale advective-diffusive processes (i.e., PEX) from other SGD pathways (Hong et al., 2017). In deeper sediment systems, where seepage occurs laterally at depth (i.e., as SGD in coarsegrained sediments), Eq. 13 can be used for the horizontal (1-D) export flux of 224 Ra, and thus horizontal water exchange rates (Q; L⋅m − 2 ⋅d − 1 ) can be estimated (Shi et al., 2019a) using Eq. 16, where A PW and A sea represents the dissolved 224 Ra activity in porewater from the seepage layer and seawater, respectively. In turn, solute fluxes may be estimated by considering the concentration difference between porewater and seawater of a given solute (Shi et al., 2019a). Similarly, a two-dimensional advective cycling model can be used to estimate water exchange and solute fluxes in sandy seabeds subject to wave and tidal pumping (Cai et al., 2020). In conclusion, application of the 224 Ra/ 228 Th disequilibrium method has been significantly expanded upon since its first introduction (Cai et al., 2012), and likely will continue to be expanded upon as more researchers apply this technique to different coastal environments. Whereas the approaches outlined above describe fluxes across the sediment-water interface, diapycnal mixing i.e., the mixing away from deep-sea sediments into the interior of the water column is one further process which can be quantified using Ra isotopes. As bottom sediments are the source of Ra, the concentration gradient above the bottom sediments can be used to calculate vertical eddy diffusivity (see Eq. 8). So far very few studies applied this method (Huh and Ku, 1998;Koch-Larrouy et al., 2015). Nevertheless it may be an important approach to investigate the fate of benthic sourced solutes in the water column.

Estimation of timescales in coastal surface waters
The functioning and vulnerability of a coastal ecosystem (e.g., accumulation of contaminants, risk of eutrophication or harmful algal blooms) are often closely related to the retention of solutes (e.g., nutrients, contaminants, suspended biomass) and thus the transport mechanisms of solutes in the system studied. Therefore, understanding the temporal scales of water and solute fluxes is crucial for coastal oceanography. Obtaining a flushing time of Ra in a coastal study is also a key parameter needed to estimate SGD through the Ra mass balance or the mixing model approach (see Section 6.1). Ra isotopes can provide instrumental information in this regard because the different decay constants of each isotope can help infer mixing timescales in coastal systems. There are two basic approaches to estimate water ages using Ra isotopes (Moore, 2015): the "mummy" model, appropriate for systems with Ra inputs occurring near the shoreline (Moore, 2000a) and the "continuous input" model, recommended for systems where Ra is added over the entire study area, such as a shallow estuary . Both approaches are based on using activity ratios of a pair of Ra isotopes to determine water apparent ages, which is the time elapsed since the water sample became enriched in Ra and isolated from the source (Moore, 2000a). Notice that these approaches are used to estimate "apparent ages", defined as the time spent since a water parcel became isolated from the Ra source. "Apparent ages" are thus not equivalent with "flushing times", although they are often used indistinctly in SGD literature (Monsen et al., 2002;Moore et al., 2006). However, evaluations of water apparent ages and flushing times might yield similar time estimates if the system is under steady-state and the comparison of flushing times and apparent ages is made across an entire study area (Tomasky-Holmes et al., 2013).
Average water apparent ages for a system (t w ) can be estimated as follows, depending on whether the Ra inputs occur at the shoreline (Mummy model; Eq. 17) or continuously (Continuous input model; Eq. 18): where AR In and AR Sys are the activity ratios of the shorter-lived Ra isotope to the longer-lived one in the source (considering all the potential sources) and in the system, respectively, and λ s and λ l are the decay constants of the shorter-lived and longer-lived Ra isotope, respectively. The term λ l can be neglected when a long-lived Ra isotope is used ( 226 Ra and/or 228 Ra). Both models assume that i) Ra activities and ARs are highest in the source (AR in > AR Sys ), ii) the AR in is constant, iii) the only Ra losses are due to mixing and radioactive decay, and iv) the open ocean contains negligible activities of the Ra isotopes used (Charette et al., 2001;Moore, 2006b;Moore, 2000a). The activity ratio of any pair of Ra isotopes can be used to estimate water apparent ages, provided that the half-life of the shorter-lived Ra isotope is appropriate for the water mixing time-scales expected in the study area (Moore, 2015). Typical water apparent ages estimated using 223 Ra and 224 Ra as the shorter-lived isotope range from ~1 to ~50 days (e.g., Dulaiova et al., 2009;Hancock et al., 2006;Krall et al., 2017;Moore et al., 2006;Rengarajan and Sarma, 2015;Sanial et al., 2015). The accuracy of these AR-based apparent water ages will depend on the analytical uncertainties associated with the isotopes used (e.g., uncertainties of 50-100% when water apparent ages are shorter than 3-5 days . It should be noted that if water masses with different Ra content mix, the relative fractions of each water should be known for reliable apparent age estimates, otherwise the apparent age will tentatively be biased towards the younger water mass (Delhez et al., 2003;Hougham and Moran, 2007).

Ra-228 as tracer of offshore solute fluxes
Solute fluxes from the coastal ocean derived from SGD, shelves and other continental sources have the potential to contribute to biogeochemical cycles of the open ocean via shelf-ocean exchange processes. However, while solute inputs to the coastal zone are relatively straightforward to obtain (e.g., Homoky et al., 2016), calculating their net input to the ocean has been a greater challenge owing to solute removal or addition processes that take place in estuaries or over the shelf. To address this shortcoming, Charette et al. (2016) proposed the use of 228 Ra as a shelf solute flux gauge, which takes advantage of the global shelf 228 Ra flux model developed by Kwon et al. (2014). The method is most easily applied where shelf-ocean exchange is primarily driven by eddy diffusion, whereby the net cross-shelf solute (S) flux can be linearly scaled with the net cross-shelf 228 Ra flux as follows (see also 6.1.): where S shelf and 228 Ra shelf are the average concentrations of the solute (e. g., trace metals, nutrients, dissolved organic carbon) of interest and 228 Ra over the shelf water column (<200 m). Because the shelf water is exchanging with the open ocean via mixing, the shelf solute and 228 Ra concentrations must be corrected for their concentrations in the open ocean (S ocean and 228 Ra ocean ).
This technique only requires paired measurements of 228 Ra and the solute of interest along a shelf-ocean transect. Given that lateral inputs of solutes to the ocean have been shown to be important on a global basis, for example as is the case with dissolved iron (Tagliabue et al., 2014), this method is particularly valuable for oceanographic field studies where solute mass balance budgets are required. However, the application of this method requires validating the assumptions of negligible advection and steady-state conditions on the scale of the tracer used, which are not valid for all the systems.

Guidelines to conduct a SGD study
This article reviews the application of Ra isotopes as tracers of SGDderived inputs of water and solutes to the coastal ocean. In this final chapter, we provide a step-by-step protocol that should serve as simplified guidelines to perform a SGD study using Ra isotopes. This protocol is based on the 7 steps illustrated in Fig. 12.
Step 1 -Definition of the objective of the study: SGD studies often focus on determining the significance of SGD in the water cycle or on its implications for coastal biogeochemical cycles. This objective determines the spatial scale of the study (e.g., nearshore vs shelf scale; evaluate spatial variability vs produce integrated estimates; spatial resolution), as well as the temporal scale (e.g., snap-shot observations vs continuous observations; base conditions vs episodic events; temporal resolution). The method(s) used to quantify SGD (e.g., water budgets, groundwater flow models, seepage meter or natural tracers such as radium, radon, salinity or heat) should be selected according to the specific objective of the study and the characteristics of the system. In case of using Ra isotopes, a clear definition of the objective will facilitate progressing towards the following steps.
Step 2 -Characterization of the study site: Identifying the hydrological (e.g., magnitude of surface water inputs, groundwater table elevation, system geology) and oceanographic characteristics (e.g., circulation drivers, system timescales, water depths and stratification) characteristics, the potential Ra sources and sinks, as well as the potential SGD pathways discharging to the study site (see Fig. 3). This characterization usually involves a combination of literature review, potentially preliminary samplings and the use of complementary techniques (e.g., salinity profiles, electrical resistivity tomography (ERT), thermal infrared (TIR) images, Rn surveys, water level measurements) from different disciplines (e.g., hydrogeology, geophysics, oceanography).
Step 3 -Construction of a conceptual model: A conceptual model should consider the main SGD pathways, the predominant sources and sinks of Ra isotopes and the characteristics of the study site (e.g., water residence time). The conceptual model is fundamental to selecting the appropriate Ra isotope(s) to be used, and to determining the quantification approach, assumptions made and overall sampling approach. It is important that researchers clearly identify the assumptions of the approach they are applying and validate them. This step is crucial to characterize the limitations of the method used, to constrain the uncertainties of the final estimates and, thus, to produce reliable and justifiable SGD estimates (see section 7.1).
Step 4 -Selection of the Ra isotope(s): Ra isotopes applied in the study should be chosen according to the target process or pathway, the enrichment in Ra of discharging groundwaters, the potential sources and sinks in the study site and the residence time of the isotope in the study area (see Sections 7.3 and 7.4).
Step 5 -Sampling Ra in all the compartments: Samples for Ra isotopes should be collected both in seawater and in all the potential endmembers (SGD from different pathways, open ocean, surface water, surface discharges, etc.). The collection of sediment samples or cores might also be required to evaluate inputs of Ra from sediments. Aside from Ra isotopes, concurrent samples for other parameters (e.g., salinity, water composition, stable isotopes (δ 2 H, δ 18 O) of water, concentrations of nutrients, metals, contaminants) might also be collected depending on the objectives of the study.
Step 6 -Quantification of SGD using Ra: The quantification approach (mass balance, endmember mixing model or export offshore) must be chosen according to the characteristics of the study site and considering the validity of model assumptions (see section 6.1). The quantification of SGD requires constraining the sources and sinks of Ra isotopes in the study area (step 6.1) and characterizing the appropriate Ra concentration in the SGD endmembers (step 6.2). In complex settings, the contribution from the different SGD pathways needs to be accurately accounted for to provide reliable SGD estimates. We note that this may not be possible in some study sites, and so the use of other tracers should be considered. Ideally, the estimates derived from Ra isotopes should be compared with other methods (e.g., seepage meters, other tracers, hydrogeological models) to validate and better constrain the results obtained.
Step 6.1 -Constraining Ra sources and sinks: The accuracy of SGD estimates strongly rely on the appropriate quantification of the relevant sources and sinks of Ra at the study site. The significance of sources and sinks (and thus the accuracy with which these parameters need to be determined) depend on the characteristics of the study site ( Fig. 8 and 9) and the Ra isotope used (see 7.3).
Step 6.2 -Determination of the SGD endmember(s): The determination of the Ra concentration in the SGD endmember(s) is often the major source of uncertainty of SGD estimates. In coastal system with a dominant SGD pathway, the Ra endmember concentration should be representative of this pathway. However, in more complex settings, the respective Ra endmembers for each of the pathways as well as its proportion of SGD contribution should be determined when possible (see Sections 6.2 and 7.5). Ra-based SGD assessments should focus their effort on accurately determining these endmembers.
Step 7 -Quantification of SGD-driven solute fluxes: A representative concentration of solutes or solute/Ra ratios in the SGD endmember/s needs to be determined to estimate SGD-driven solute fluxes. Importantly, the solute composition of different SGD pathways may be considerably different and specific endmembers should be determined for the different components. This is especially important for studies that aim to quantify SGD-driven solute fluxes (nutrients, metals, contaminants, etc) or to evaluate the significance of SGD in coastal biogeochemical cycles (see Section 7.6).

Knowledge gaps and research needs
Despite the advances in the use of Ra isotopes as SGD tracer over the last decades reviewed in this article, several major research gaps remain open. We believe that the following issues should be addressed in future studies to improve the use of Ra isotopes as a tracer in SGD studies.
Discriminating SGD pathways As highlighted in this review, identifying the SGD pathways or mechanisms that are most significant in contributing Ra isotopes at a specific study site is crucial to obtaining realistic SGD estimates. There is thus a need to obtain a conceptual understanding of the systems prior to the application of Ra isotopes to trace SGD. The combination of different Ra isotopes can be instrumental to quantify not only the total magnitude of SGD, but also the respective water flows supplied by different pathways. Their discrimination will allow obtaining accurate estimations of the solute fluxes supplied by SGD, which will decisively contribute towards a better understanding of the role SGD plays in coastal biogeochemical cycles.

Multi-method studies
Comparisons of Ra-based SGD estimates with independent methods (e.g., other tracers, seepage meters, hydrological modelling) are a key step towards the validation and refinement of Ra-derived estimates. However, estimates derived from different methods are not always directly comparable because different approaches often capture different components of SGD or integrate over unique spatio-temporal scales. Using Ra isotopes to discriminate SGD pathways will thus facilitate multi-method comparisons and, at the same time, multi-method studies can contribute to distinguishing different SGD components.

Spatial and temporal Ra variability within the subterranean estuary
There are still few studies constraining the spatial and temporal variability of Ra isotopes in the subsurface and Ra behaviour in the subterranean estuary is not properly understood. Field measurements should also be combined with reactive transport models to understand the spatio-temporal variability of Ra isotopes within the subterranean estuary and the role of physical mechanisms driving groundwater.

Development of new analytical techniques for the measurement of Ra isotopes
Most of the methods applied to measure Ra isotopes require preconcentration of Ra from large sample volumes (usually 10-200 L) and measurements via radiometric techniques such as the RaDeCC system, gamma spectrometry, ICP-MS and TIMS systems. However, there is still the need for an accurate and easy methodology to measure the four Ra isotopes in small volumes of water (<1 L), allowing for the routine analysis of a large number of samples, and for determining Ra concentrations in compartments where only small volumes can be sampled (e. g., porewaters, deep ocean). These developments would facilitate the improvement of the spatial and temporal resolution of Ra samples collected both in the subterranean estuary and the ocean and would allow, for instance, to assess Ra concentration changes over time, obtaining detailed Ra profiles in the water column or high-resolution Ra distributions in the subterranean estuary. These methodological improvements are thus crucial to move towards a better understanding of the different SGD pathways and their spatio-temporal scales.
Uncertainties of Ra-derived SGD estimates SGD fluxes obtained from Ra isotopes have large uncertainties that are frequently overlooked. Propagated uncertainties associated with estimates should always and clearly be reported in SGD studies to facilitate understanding of the precision of tracer approaches. Importantly, an accurate assessment of uncertainties should not only consider the uncertainties linked to individual parameters (e.g., analytical errors of Ra measurements, the standard deviation of Ra averages), but also those errors linked to the conceptualization of the system (e.g., assumption of steady state, selection of the endmember, assumptions linked to mixing loss assessments).
Towards more robust Ra-derived SGD estimates All studies using Ra isotopes to quantify SGD processes are based on numerous assumptions that should be validated in order to produce accurate SGD estimates. In other words, most of the investigations using Ra isotopes produce SGD estimates, but the estimated fluxes are only meaningful and reliable if the assumptions and uncertainties of the model are properly understood, acknowledged, quantified, and accounted for. Authors working on SGD-related investigations should always be aware of the limitations of the tracers and approaches used. In some cases, traditional approaches (or the tracer itself) might not be appropriate for the system studied or the objective of the investigation, requiring the development of more complex models or the use of alternative/complementary methods.

Declaration of Competing Interest
None