Radium Mass Balance Sensitivity Analysis for Submarine Groundwater Discharge Estimation in Semi-Enclosed Basins: The Case Study of Long Island Sound

Estimation of submarine groundwater discharge (SGD) to semi-enclosed basins by Ra isotope mass balance is herein assessed. We evaluate 224Ra, 226Ra, and 228Ra distributions in surface and bottom waters of Long Island Sound (CT-NY, United States) collected during spring 2009 and summer 2010. Surface water and bottom water Ra activities display an apparent seasonality, with greater activities during the summer. Long-lived Ra isotope mass balances are highly sensitive to boundary fluxes (water flux and Ra activity). Variation (50%) in the 224Ra, 226Ra, and 228Ra offshore seawater activity results in a 63–74% change in the basin-wide 226Ra SGD flux and a 58–60% change in the 228Ra SGD flux, but only a 4–9% change in the 224Ra SGD flux. This highlights the need to accurately constrain long-lived Ra activities in the inflowing and outflowing water, as well as water fluxes across boundaries. Short-lived Ra isotope mass balances are sensitive to internal Ra fluxes, including desorption from resuspended particles and inputs from sediment diffusion and bioturbation. A 50% increase in the sediment diffusive flux of 224Ra, 226Ra, and 228Ra results in a ∼30% decrease in the 224Ra SGD flux, but only a ∼6–10% decrease in the 226Ra and 228Ra SGD flux. When boundary mixing is uncertain, 224Ra is the preferred tracer of SGD if sediment contributions are adequately constrained. When boundary mixing is well-constrained, 226Ra and 228Ra are the preferred tracers of SGD, as sediment contributions become less important. A three-dimensional numerical model is used to constrain boundary mixing in Long Island Sound (LIS), with mean SGD fluxes of 1.2 ± 0.9 × 1013 L y–1 during spring 2009 and 3.3 ± 0.7 × 1013 L y–1 during summer 2010. The SGD flux to LIS during summer 2010 was one order of magnitude greater than the freshwater inflow from the Connecticut River. The maximum marine SGD-driven N flux is 14 ± 11 × 108 mol N y–1 and rivals the N load of the Connecticut River.


INTRODUCTION
Submarine groundwater discharge (SGD) is a component of the hydrologic cycle and can act as an important vector for the transport of nutrients, carbon, trace elements, and pollutants to the coastal ocean (Moore, 2010;Knee and Paytan, 2011). SGD includes both terrestrial, meteorically-derived groundwater driven by a positive onshore hydraulic gradient and marine (i.e., saline) groundwater, driven by a variety of physical forcing mechanisms including density, tide, and wave driven flow (Santos et al., 2012). Naturally occurring radium isotopes are powerful tracers of SGD, as brackish groundwaters are typically enriched in dissolved Ra isotopes by several orders of magnitude over seawater (Swarzenski, 2007;Charette et al., 2008). The Ra quartet spans a wide range of half-lives ( 223 Ra = 11.4 d, 224 Ra = 3.66 d, 226 Ra = 1600 y, and 228 Ra = 5.75 y) and has thus been applied to trace and quantify inputs of SGD to the ocean on a variety of scales (Moore, 2010). For a given area, evaluation of the Ra source terms (e.g., rivers, particle desorption, diffusion, and bioirrigation), and Ra sinks (e.g., mixing, radioactive decay) can be used to quantify SGD. A Ra flux supplied by SGD is typically invoked to explain any imbalance between Ra sink and Ra source fluxes. Consequently, a SGD-driven Ra flux can be converted to a volumetric water flow with a proper characterization of the SGD endmember Ra activity. Multiple Ra isotopes are often used to quantify SGD; however, differences between shortlived and long-lived Ra isotope mass balances are often poorly constrained or not fully understood (Moore et al., 2006;Beck et al., 2007Beck et al., , 2008Garcia-Solsona et al., 2008;Knee et al., 2016;Tamborski et al., 2017b).
Each Ra isotope is sensitive to different source and sink terms, reflecting the time-scale of a particular process with respect to the Ra isotope half-life. For example, short-lived 223 Ra and 224 Ra are more rapidly lost via radioactive decay, while longlived 226 Ra and 228 Ra are not. Thus, the water column inventory of 223 Ra and 224 Ra must be well constrained to evaluate these short half-life tracers. Furthermore, flow paths of varying timescales may have unique short-lived and long-lived Ra activities, and thus these isotopes may trace different SGD and porewater exchange flow paths (Rodellas et al., 2017). In addition, different geological matrices can have unique ratios of uranium ( 223 Ra, 226 Ra), and thorium ( 224 Ra, 228 Ra) series isotopes, thus enabling the identification of different geologic sources Swarzenski, 2007).
This article synthesizes sediment, surface water, bottom water and groundwater Ra isotope data that has been previously collected in the semi-enclosed tidal estuary of Long Island Sound (LIS; Krishnaswami et al., 1982;Copenhaver et al., 1993;Turekian et al., 1996;Garcia-Orellana et al., 2014;Bokuniewicz et al., 2015;Tamborski et al., 2017a,b). In addition, we present new data on long-lived 226 Ra and 228 Ra, previously collected during 2009 and 2010. Here, we use short-lived 224 Ra and long-lived 226,228 Ra to quantify total SGD to LIS by mass balance. The main objective of this article is to evaluate the sensitivity of SGD estimated from Ra isotope mass balances. We evaluate Ra source and sink terms, and provide general recommendations on best-practices for Ra mass balances in semi-enclosed basins for future studies. We conclude with a revised estimate of SGD-driven NO 3 − loadings to LIS.

Study Site
Long Island Sound is a tidal estuary bound by New York City at its western end, the southern shore of Connecticut at its northern boundary and the northern shore of Long Island (NY) at its southern boundary (Figure 1). At the beginning of the 20th century, urbanization, and pollution from the Metropolitan New York area led LIS to be nicknamed the "urban sea" (Koppelman et al., 1976;Latimer et al., 2013). Bottomwater hypoxia and eutrophication have been linked to excess nitrogen inputs from wastewater, sewage effluent, river inputs, and groundwater (NYSDEC and CTDEP, 2000). Indeed, Long Island's coastal embayments are vulnerable to excess nitrogen loading from submarine groundwater inputs (Bokuniewicz, 1980;Capone and Bautista, 1985). However, the volume of total SGD to LIS is not well known, despite increasing recognition of this important pollutant pathway. The volume of terrestrial SGD to LIS from both Connecticut and New York shorelines is generally well constrained from hydrogeologic models (Buxton and Modica, 1992;Scorca and Monti, 2001;Suffolk County, 2015). Ra isotope mass balances in LIS' smaller embayment, Smithtown Bay (Bokuniewicz et al., 2015;Tamborski et al., 2017b), and for the entire LIS (Garcia-Orellana et al., 2014) all suggest that total SGD is dominated by marine groundwater inputs. Tamborski et al. (2017b) found that SGD within the first 200 m of the shoreline, determined from short-lived radionuclide mass balances, was up to ∼55% greater than SGD estimates determined from long-lived radionuclides and physical seepage meter measurements. To explain this difference, they suggested that wave and tidal circulation SGD flow paths captured short-lived radionuclide fluxes due to their faster regeneration rates within sediments, while these flow paths would not capture 226 Ra and 228 Ra. Interestingly, the total SGD flux to Smithtown Bay determined by Bokuniewicz et al. (2015) from 224 Ra, which includes inputs to the entire bay, was one to three orders of magnitude greater than the total SGD flux within the first 200 m of the shoreline. This led Tamborski et al. (2017b) to hypothesize that there may be additional, deeper SGD flow paths farther offshore in Smithtown Bay and LIS, from either the deeper Magothy aquifer or seawater circulation through offshore permeable sediments driven by density-forcing mechanisms, in addition to short-scale circulation fluxes driven by bioturbation and wave pumping.

Analytical Methods
Surface and deep-water samples from LIS were collected aboard the R/V Seawolf during 24-30 April 2009, 29 July-04 August 2009, and 03-12 August 2010. These samples, along with an analysis of the short-lived Ra isotopes, are described in Garcia-Orellana et al. (2014). Separately, groundwaters were collected from intertidal cluster wells from a coastal bluff and a barrier beach subterranean estuary, during spring and FIGURE 1 | The study site of Long Island Sound, bound between Long Island (NY), and Connecticut. Surface and bottom sampling stations are indicated by white circles. The mouth of the Connecticut River is indicated by a black star; Orient Point is indicated by an orange star. Long Island Sound is subdivided into the western (light blue), central (green), and eastern (red) basins, after Garcia-Orellana et al. (2014). Purple numbers correspond to the cross-section numbers in Table 1. Mean water exchange transports through cross-sections 18 113, 180, and 257 (N-S oriented dark blue lines) are derived from the three-dimensional model of Crowley (2005); the two-dimensional grid used in the model calculations is shown for reference and overlaps the land-ocean boundary.
summer 2014 and 2015. Additional coastal stations were sampled via drive-point piezometer at the low tide shoreline. These groundwater samples, including an analysis of the Ra quartet, are described in Tamborski et al. (2017a).
Briefly, water samples (seawater = 20-60 L; groundwater = 1-4 L) were filtered through MnO 2 impregnated acrylic fibers at a flow rate of <1 L min −1 to quantitatively extract dissolved Ra from solution onto the fiber. Short-lived 223 Ra and 224 Ra were measured using a Radium Delayed Coincidence Counter (RaDeCC; Moore and Arnold, 1996) immediately after sample collection. An additional count was performed 1 month after sample collection to measure 228 Th, to determine the amount of unsupported (excess) 224 Ra. Subsequently, the Mn-fibers were leached in a HCl and Hydroxylamine hydrochloride mixture and Ra was co-precipitated with BaSO 4 . Precipitates were removed from the leachate, stored in glass vials and sealed for >3 weeks. Samples were counted on a Canberra Intrinsic Ge well detector, where the activity of long-lived 226 Ra and 228 Ra was determined from the 352 keV ( 214 Pb) and 911 keV ( 228 Ac) photopeaks, respectively. Gamma counting efficiencies were determined from NBS Standard Reference Material 4350B.

Estimates of Water Transport
LIS is made up of three basins (western, central, and eastern); water flux between basins and within LIS's boundaries were estimated to evaluate Ra mixing. Crowley (2005) determined that the transport through cross-sections to the west of Orient Point had a long term mean of 2.05 × 10 13 L y −1 throughout the LIS basin (section 257; Figure 1), as determined from three-dimensional numerical simulations of circulation with meteorological and boundary sea-level forcing. The crosssectionally-averaged transport exhibited significant temporal variability which was spatially coherent, with residence times on the order of several weeks to 2 months (Crowley, 2005). The differences in cross-sectionally-averaged transport can be attributed to localized effects of river discharge and to storage of water between sections. Water exchange fluxes were derived from the model of Crowley (2005) at cross-sections that approximately separate LIS into its three unique basins (western, central, and eastern; Figure 1) over a period between 2009 and 2010. Long term mean exchange transports through cross-sections exhibit an east to west net transport ( Table 1). Note that these mean exchange transports have a standard error on the order of 10% because of the large temporal variability in the exchange. The mean surface freshwater inflow to the entire basin is on the order 1.34 × 10 13 L y −1 , of which the Connecticut River at the far end of the Sound contributes approximately 1.19 × 10 13 L y −1 , depending upon the seasonal conditions. The exchange transports at section 257 are corrected to accommodate the surface freshwater inflow between sections 18 and 257 (Figure 1). It should be noted that section 257 is located approximately at the longitudinal position of the Connecticut River, and that recent dye simulations emphasize that much of the Connecticut River waters exit to the east (Jia and Whitney, 2019).

Ra Distribution
Long-lived 226 Ra and 228 Ra activities generally increased from east to west in LIS ( Table 2). Compared to surface waters, deep-water Ra activities were greater on average for the western basin and similar to surface waters for the central and eastern basins ( Table 2). Surface and deep-water Ra activities generally Cross-section numbers correspond to the sections depicted on Figure 1.

Ra Mass Balance
Long-lived 226 Ra and 228 Ra mass balances for LIS were constructed after the short-lived 224 Ra mass balance developed by Garcia-Orellana et al. (2014) for spring 2009 and summer 2010, for the entire LIS basin. We have further developed Ra mass balances for each individual basin (western, central, and eastern; Figure 1). The Ra mass balances have been updated to reflect new terms, as described below. Briefly, a surface and deep box Ra inventory (dpm) is calculated for the western, central, and eastern basins of LIS ( Table 4). The Ra inventory is calculated as the product of the water volume within each basin and the mean Ra activity of the basin. The Ra mass balance is written as: Where the left-hand side represents Ra sinks and the right-hand side represents Ra sources. The Ra mass balance includes loss from mixing (J out ) and radioactive decay (J decay ). Ra sources include mixing (J in ), riverine input (J river ), the desorption of Ra from resuspended particles (J desorp ), molecular diffusion and bioturbation (J diffusion ), and SGD (J SGD ). Each term in Eq. 1 is assessed for 224 Ra, 226 Ra, and 228 Ra for the western, central, eastern, and total basins during spring 2009 and summer 2010, thus resulting in 24 unique Ra mass balances. All Ra mass balances are summarized in the Supplementary Material (Supplementary Table S1). Each term in Eq. 1 is described in further detail below, including a sensitivity analysis on the final estimated Ra flux supplied by SGD (section "Mass balance sensitivity").

Boundary Fluxes
Boundary fluxes determined from a three-dimensional numerical model (section "Estimates of water transport") are used to quantify the exchange of water and Ra isotopes between New York Harbor and the East River with the western basin of LIS (section 18), between the western basin and the central basin of LIS (section 113), between the central basin and the eastern basin of LIS (section 180), and between the eastern basin of LIS and Block Island Sound (section 257; Table 1 and Figure 1). We use the same boundary flux for spring 2009 and summer 2010, although we acknowledge that there may be minor differences in seasonal water volume and thus exchange. Ra mixing fluxes (J in and J out ) are calculated from the sectional mean water exchange transports and the associated endmember Ra activity. Because there is a net east to west water transport in LIS (Table 1), we use the mean deep-water Ra activities ( Table 2) to represent Ra mixing from east to west. Mean Ra surface water activities ( Table 2) are used to represent Ra mixing from west to east, for each respective basin under consideration. Boundary fluxes are summarized in Figure 4.
The 224 Ra, 226 Ra, and 228 Ra activities of the East River are, respectively, taken as 9 dpm 100 L −1 (Garcia-Orellana et al., 2014), 10.7 dpm 100 L −1 (Li et al., 1977), and 64 dpm 100 L −1 (Turekian et al., 1996) for both seasons. Block Island Sound seawater was previously sampled for long-lived Ra at 5 m and 25 m depth during summer 1991 (St101; Turekian et al., 1996). This same station was sampled during spring 2009 in this study. We use a 224 Ra, 226 Ra, and 228 Ra activity of 3.6, 5.3, and 29 dpm 100 L −1 , respectively, for spring 2009 and activities of 3.6, 8.5, and 47 dpm 100 L −1 , respectively, for summer 2010. It is important to note that only one sample is used to characterize the Ra endmember for each of these mixing input terms.
The Connecticut River is the largest river entering LIS. Connecticut River discharge was estimated as the mean discharge over the two-week period preceding each seasonal LIS survey, based on the average residence time of the eastern basin (Crowley, 2005). The Connecticut River discharge is taken as 2.52 × 10 13 L y −1 during spring 2009 and 0.68 × 10 13 L y −1 during summer 2010 from the USGS gaging station at Middle Haddam (ID 01193050). Note that the spring 2009 and summer 2010 Connecticut River discharges were, respectively, above and below the mean annual river discharge (∼1.2 × 10 13 L y −1 ). The Connecticut River is only considered in the eastern basin and total basin Ra mass balances. Several other rivers enter LIS from both Long Island (e.g., Nissequogue River) and Connecticut (e.g., Housatonic River, Thames River) shorelines, which together account for a cumulative discharge on the order of ∼0.4 × 10 13 L y −1 (Koppelman et al., 1976); here we assume that one third of this water flux enters each of the three LIS basins. The riverine Ra flux to LIS (J river ) is calculated as the product of the river discharge to each basin and the dissolved Ra Frontiers in Environmental Science | www.frontiersin.org activity of the Connecticut River mouth (Figure 1). Surface water salinity at the station sampled near the mouth of the Connecticut River was ∼28; therefore, the Ra flux estimated here intrinsically includes the desorption of Ra from suspended riverine particles. This calculation assumes that the Ra endmember at the mouth of the Connecticut River is representative of all rivers entering LIS, as done by Turekian et al. (1996); riverine fluxes are summarized in Figure 4.

Internal Fluxes
The desorption of long-lived 226 Ra and 228 Ra from resuspended particles (J desorp ) throughout LIS is negligible due to the slow regeneration time of 226 Ra and 228 Ra in sediments. The desorption of 224 Ra is estimated from suspended particle concentrations for spring 2009 (1.5 ± 1.3 mg L −1 ) and summer 2010 (4.0 ± 2.5 mg L −1 ), and a surface-exchangeable 224 Ra activity of 0.75 dpm g −1 . Desorption of 224 Ra from tidally resuspended particles is 1.6 ± 1.4 × 10 4 dpm m −2 y −1 for spring 2009 and 4.4 ± 2.7 × 10 4 dpm m −2 y −1 for summer 2010 (Garcia-Orellana et al., 2014; Figure 4). Sediment diffusion and bioturbation have been shown to be an important source of short-lived Ra isotopes to LIS, due to their relatively rapid regeneration rates within sediments (Garcia-Orellana et al., 2014). This source term should be relatively less important for the long-lived Ra isotopes, due to their slower regeneration rates in sediments. Tamborski et al. (2017b) conducted sediment core incubation experiments (at 20 • C and 3 • C; oxic, and hypoxic conditions) to determine the sediment flux of long-lived Ra isotopes for a sandy (intertidal) core and a silty (offshore) core from LIS. 228 Ra sediment fluxes in LIS were previously determined from the disequilibrium between solid-phase 228 Ra and its parent 232 Th (Cochran, 1979;Turekian et al., 1996), which integrates over several half-lives of 228 Ra and thus represents a mean-annual 228 Ra flux. Each of these methods includes a Ra flux from both molecular diffusion and bioirrigation. These various flux estimates were averaged together for coarse-grained sediments and fine-grained sediments. Fluxes were further separated by temperature (warm = summer; cold = early spring) and oxygen content (hypoxic = summer; oxic = early spring) to differentiate between spring and summer conditions ( Table 5). Sediment-mediated 228 Ra and 226 Ra inputs were calculated using a LIS bottom surface area of 0.90 × 10 9 m 2 , 1.08 × 10 9 m 2 , and 0.96 × 10 9 m 2 for the western, central and eastern basins. Further, we assume a fine-grained and coarsegrained surficial sediment distribution for the western (70% fine, 30% coarse), central (50% fine, 50% coarse), and eastern (0% fine, 100% coarse) basins following Poppe et al. (2000). Sediment Ra fluxes are summarized in Figure 4.
Assuming steady-state, the radioactive decay of 224 Ra and 228 Ra is calculated as the product of the respective Ra inventory for each basin and the Ra isotope decay constant (0.189 d −1 for 224 Ra and 0.121 y −1 for 228 Ra). Decay of 226 Ra is negligible due to its long half-life (Figure 4). The excess Ra inventory in each basin ( Ra sinks -Ra sources; Eq. 1) is presumed to be balanced by SGD (J SGD ). A negative SGD flux implies either that Ra sources approximately balance Ra sinks (within the propagated uncertainties), or that one (or more) known Ra flux is improperly characterized. SGD-derived Ra fluxes are converted into volumetric water flows by dividing the Ra flux by the Ra activity of the SGD endmember, sampled from intertidal wells and shoreline piezometer stations along both NY and CT shorelines (Garcia-Orellana et al., 2014;Tamborski et al., 2017a). Selection of the SGD Ra endmember is discussed below (section "SGD volumetric water flow to Long Island Sound").

Mass Balance Sensitivity
Analysis of Ra sources and sinks reveals that boundary fluxes dominate the long-lived Ra isotopes while sediment contributions (diffusion and desorption) and radioactive decay  Groundwater samples are classified by salinity: fresh groundwater = 0-1; brackish groundwater = 4-24; and marine groundwater = 25-28. Error bars indicate ± 1 standard deviation from the mean. Note that brackish groundwaters from Connecticut are not included in this analysis. *Groundwaters from Tamborski et al. (2017a). Groundwaters from Copenhaver et al. (1993) and Krishnaswami et al. (1982).ˆˆGroundwaters from Garcia-Orellana et al. (2014) and analyzed in this study.
dominate the short-lived Ra isotopes (Figure 4). Ra mass balance sensitivity to each of these terms is discussed below. Sediment diffusion, bioirrigation, and desorption support ∼50% of the 224 Ra inventory for LIS but only ∼6% of the 226 Ra inventory and ∼10-20% of the 228 Ra inventory (Figure 4). The 224 Ra sediment flux is constrained from sediment incubation experiments from five different locations (Garcia-Orellana et al., 2014); the 228 Ra sediment flux is constrained from four sediment incubations and two estimates from 228 Ra: 232 Th disequilibrium (Cochran, 1979;Turekian et al., 1996). Is this representative of the entirety of LIS? Few Ra-based SGD studies capture sedimentmediated Ra inputs from more than a few measurements (Beck et al., 2007Garcia-Solsona et al., 2008;Rodellas et al., 2012Rodellas et al., , 2015Cai et al., 2014;Tamborski et al., 2017b).
The molecular diffusive flux of Ra from shallow porewater to overlying seawater is governed by the Ra concentration gradient between porewater and seawater (Fick's first law). Bioturbation and bioirrigation can further facilitate Ra transport by enhancing the effective sediment surface-area, thus complicating diffusive flux estimates. Due to its short half-life, 224 Ra regenerates rapidly in sediments from the decay of its particle-reactive parent 228 Th. The activity of 230 Th and 232 Th (the parents of 226 Ra and 228 Ra) are similar in LIS sediments, resulting in similar production rates in muddy and sandy sediments (Cochran, 1979). Therefore, seasonal redox changes may affect the diffusive flux of short-lived and long-lived Ra isotopes differently, as the redox interface migrates between bottom waters and the shallow subsurface (Garcia-Orellana et al., 2014). This seasonal control impacts Ra diffusion and more importantly, the influence of the benthic fauna, which are presumed to be less active during colder periods. We have attempted to capture this variability by using different sediment incubations for spring and summer conditions ( Table 5).
Are the incubations used to quantify diffusion and bioirrigation from sediments accurate? 228 Ra fluxes for muddy LIS sediments, determined from the deficit of solidphase 228 Ra relative to its parent 232 Th (Cochran, 1979;Turekian et al., 1996), are in general agreement with 228 Ra fluxes determined from sediment chamber incubations (33-82 dpm m −2 d −1 vs 12-207 dpm m −2 d −1 ). The solid-phase 228 Ra: 232 Th approach integrates over the half-life of 228 Ra, therein representing a mean-annual 228 Ra flux; therefore, these estimates seem reasonable at the seasonal time-scale considered here. A recent study found the traditional sediment incubation approach for 224 Ra flux determination to be similar to 224 Ra fluxes determined from 224 Ra: 228 Th disequilibrium (Shi et al., 2018). Thus, it seems that the incubations used to quantify the Ra flux from diffusion and bioirrigation are accurate, although it remains to be seen how representative these several cores are for the entirety of LIS. A 50% increase in the diffusive flux of 224 Ra, 226 Ra, and 228 Ra, while keeping all other parameters equal (Eq. 1), results in a ∼30% decrease   in the 224 Ra SGD flux to LIS, but only a ∼6% and ∼10% decrease in the 226 Ra and 228 Ra SGD flux, respectively (total-basin mass balances). More realistic is a decrease in the 226 Ra and 228 Ra sediment diffusive flux. Many studies neglect diffusion of long-lived Ra isotopes altogether (Rama and Moore, 1996;Beck et al., 2007Beck et al., , 2008. Indeed, exclusion of this term would result in a 6-7% increase in the 226 Ra SGD flux and a 10-18% increase in the 228 Ra SGD flux (total-basin mass balances). Desorption of 226 Ra and 228 Ra from resuspended sediments is assumed negligible. Desorption of 224 Ra from resuspended sediments is estimated to account for 14-19% of the total 224 Ra inventory for the total-basin mass balance. Rodellas et al. (2015) note that it is difficult to accurately constrain the Ra flux from resuspended sediments, which can be a significant Ra source in embayments with fine-grained sediments. The surface-exchangeable 224 Ra estimated for LIS (0.75 dpm g −1 ) is based on shallow porewater (0-2 cm) 224 Ra activities from two sediment cores (Garcia-Orellana et al., 2014) and a mean K d for Ra of 50 L kg −1 (Cochran, 1979;Sun and Torgersen, 2001). Even more difficult to constrain is the rate at which sediments are resuspended to release 224 Ra into the water column. Here, sediments are assumed to be resuspended on tidal time-scales. This term represents the second largest source of uncertainty in the 224 Ra mass balance (Figure 4). A 50% increase in the 224 Ra desorption flux, while keeping all other parameters equal (Eq. 1), results in a ∼14-19% decrease in the 224 Ra SGD flux to LIS (total-basin mass balances).
With respect to the loss by radioactive decay, the basinwide 224 Ra inventory is more sensitive to this term than are the long-lived Ra inventories. Indeed, decay of 224 Ra dominates over mixing losses, regardless of the basin or season (Figure 4). This implies that a large number of water column samples are necessary to accurately capture the (near) instantaneous 224 Ra inventory. The 224 Ra inventory integrates over the half-life of 224 Ra, such that this inventory will accurately reflect the times of the year when the sampling was conducted, but this inventory may not be representative of a seasonal or annual 224 Ra flux. Importantly, this suggests that multiple sampling campaigns are necessary to capture any seasonality in the 224 Ra inventory for SGD flux determination. Fortunately, 224 Ra decay is the simplest flux term to constrain, as it is only dependent upon measuring 224 Ra concentration. This is the opposite case for long-lived 226 Ra and 228 Ra, where a lower number of water column samples may be adequate to capture the long-lived Ra water column inventory, while mixing terms must be wellconstrained. Seasonality in long-lived Ra sources and sinks may be difficult to resolve with these isotopes because they integrate over long temporal scales.
Multiple water column Ra samples are required at the boundaries (East River and Block Island Sound in this case) to accurately constrain the boundary mixing Ra flux (Figure 4). The boundary water flux may be difficult to quantify in environments where numerical models or instrument deployment (e.g., ADCP) are unavailable. The 224 Ra inventory will be less sensitive to boundary mixing when the time-scale of mixing is sufficiently large with respect to the half-life of 224 Ra (Figure 4). Western and eastern basin long-lived Ra isotope mass balances for spring 2009 reveal that Ra sources approximately balance Ra sinks, without the need to invoke SGD. However, the uncertainties on these mass balances are quite large, owing to the uncertainty of the mixing endmember activity (East River for the western basin and Block Island Sound for the eastern basin). The Ra flux into LIS from mixing with the East River (J river ) and Block Island Sound (J in ) were each determined from one sampling station (Li et al., 1977;Turekian et al., 1996). A 50% change in the East River 224 Ra, 226 Ra, and 228 Ra activity, while keeping all other parameters equal (Eq. 1), results in only a 1-2% change in the 224 Ra SGD flux, and a 6-10% change in both the 226 Ra and 228 Ra SGD fluxes (total-basin mass balances). This is minor in comparison to the Ra flux exchanged between Block Island Sound. A 50% change in the Block Island Sound 224 Ra, 226 Ra, and 228 Ra activity, while keeping all other parameters equal (Eq. 1), results in only a 4-9% change in the 224 Ra SGD flux, but a 63-74% change in the 226 Ra SGD flux and a 58-60% change in the 228 Ra SGD flux. This highlights the critical importance of accurately determining boundary water fluxes and endmember activities, especially for the long-lived Ra isotopes ( Table 1).
The western and eastern basin Ra activities are properly characterized (Table 2); thus, the central basin Ra mass balances are more adequately constrained for mixing compared to the total-basin mass balances. We emphasize that, by capturing mixing gains and losses, the central basin Ra mass balances adequately characterize the Ra-derived SGD flux, regardless of the Ra isotope used (Figure 4). We note that the real water exchange uncertainties may be higher than what is determined from the numerical model (Table 1); uncertainty in determining boundary mixing will produce large uncertainties in the longlived Ra SGD fluxes.

The Magnitude and Significance of Submarine Groundwater Discharge to Long Island Sound
Endmember Mixing Model: An Alternative Approach to Estimate SGD Results from the Ra mass balance in LIS reveal that SGD estimates derived from long-lived Ra isotopes are highly sensitive to boundary exchange processes (Figure 4). Potential uncertainties on the estimation of these boundary exchange fluxes might thus compromise the SGD estimates. An endmember mixing model is an alternative approach to estimate the importance of SGD as a source of 226 Ra and 228 Ra to LIS, as proposed by Moore (2003). In order to simplify the system, and considering that other sources such as rivers, desorption, diffusion and decay are less important (Figure 4), we assume that long-lived Ra isotopes in LIS are exclusively supplied by SGD and open ocean water through boundary exchange. This approach is therefore qualitative and serves to constrain the relative magnitude of SGD only as a comparison to the mass balance. Coastal groundwater endmembers from Long Island and Connecticut shorelines have relatively distinct 228 Ra/ 226 Ra ratios (Table 3 and Figure 3) and are thus treated as separate sources of Ra. We may thus approximate the relative contribution of Ra measured in LIS from Long Island groundwater (f LI ), Connecticut groundwater (f CT ), and seawater (f sea ) using a three-endmember mixing model (Moore, 2003): 226 Ra sea * f sea + 226 Ra LI * f LI + 226 Ra CT * f CT = 226 Ra measured (4) where Ra measured is the measured 226,228 Ra activity of LIS. The three-endmember mixing model is only applied to bottom water samples from the western and eastern basins, in order to reduce riverine contributions and uncertainty in selecting a proper seawater (i.e., boundary) endmember. For the eastern basin, the seawater endmember is taken as the seasonal deep-water Ra activity of Block Island Sound. For the western basin, the seawater endmember is taken as the seasonal deep-water Ra activity of the central basin ( Table 2). We note that the relatively large standard deviation of the endmember 228 Ra/ 226 Ra ratios ( Table 3) reflects natural variability in endmember composition and thus produces considerable uncertainties in the endmember mixing models (∼30%), such that these results should be interpreted as qualitative only.
On average, Long Island groundwater contributed to 5% (spring 2009), and 6% (summer 2010) of the measured eastern basin Ra activities, with negligible contributions from Connecticut groundwater. In contrast, Long Island and Connecticut groundwaters contributed approximately equal Ra proportions to the western basin during summer 2010 (4 and 5%, respectively), while Connecticut dominated western inputs during spring 2009 (10%) with negligible inputs from Long Island. Multiplying these percentages by the volume of water in LIS (5.20 × 10 13 L; Table 4) and dividing by a residence time of 60 days (Crowley, 2005) results in an SGD flux on the order of ∼1-3 × 10 13 L y −1 . While this is a clear simplification of the system and thus the results can only be used in a qualitative manner, this independent approach helps constrain the relative magnitude of SGD to LIS and further helps to constrain the groundwater Ra endmember. These qualitative results demonstrate that ≥90% of the long-lived Ra isotopes in LIS are derived from boundary exchange, highlighting the critical importance of properly characterizing boundary endmembers and water exchange transports. Therefore, when boundary exchange mixing is not well constrained, it is not recommended to use 226 Ra and 228 Ra to quantify SGD in semi-enclosed basins.

SGD Volumetric Water Flow to Long Island Sound
The Ra mass balances indicate a significant amount of 224 Ra, 226 Ra, and 228 Ra unaccounted for, that must be balanced by inputs from SGD (Figure 4). The SGD Ra flux to each basin is converted into a water flow by dividing by the SGD Ra endmember activity. Results from the three-endmember mixing analysis (section "Endmember mixing model: An alternative approach to estimate SGD") reveals SGD from Connecticut is insignificant in the eastern basin and is dominated by Long Island marine groundwater; thus, Long Island groundwater is used as the Ra endmember for the eastern basin ( Table 3). The western basin was impacted by Connecticut groundwater during spring 2009, and therefore we use the mean Connecticut marine groundwater Ra activity as an endmember during spring. In contrast, the western basin was impacted by approximately equal proportions of Long Island and Connecticut groundwater during summer 2010, and therefore an average of each endmember is used. For the central basin, we simply use an average of the Long Island and Connecticut marine groundwater Ra activities (salinity 25-28) for SGD flux determination ( Table 3). SGD determination is highly sensitive to the selection of the groundwater endmember (Cook et al., 2018); endmember sensitivity is not evaluated in this study. Ra activities in the SGD endmembers are assumed to be constant throughout the year and therefore we only consider seasonal differences in SGD flux (Luek and Beck, 2014). Marine groundwater 224 Ra activities from LIS are not seasonal (Tamborski et al., 2017a), as 224 Ra will be quickly regenerated within shallow marine sediments, integrating all SGD flow paths. Seasonality in marine SGD to LIS may be driven, in part, by seasonal differences in density-dependent dispersive mixing within permeable sediments (Tamborski et al., 2017a), and movement of the freshwater-saltwater interface in response to regional precipitation (Michael et al., 2005;Tamborski et al., 2017b).
Volumetric SGD flows are summarized in Figure 5. Averaging all three Ra isotope mass balances for the total LIS basin (±standard deviation) results in an SGD flux of 1.2 ± 0.9 × 10 13 L y −1 during spring 2009 and 3.8 ± 0.7 × 10 13 L y −1 during summer 2010. These estimates are within the range of estimates qualitatively determined from the three-endmember mixing model and slightly lower than the 224 Ra-derived SGD flux of 3.2-7.4 × 10 13 L y −1 , previously estimated by Garcia-Orellana et al. (2014). As noted above, the long-lived Ra isotope mass balances are significantly impacted by boundary mixing (section "Mass balance sensitivity"). The central basin Ra mass balances accurately constrain the 226 Ra and 228 Ra endmember activities of the western and eastern basins ( Table 2), which account for the 226,228 Ra mixing sources and sinks, assuming that the water flow estimations are accurate ( Table 1). The average SGD flux to the central basin of LIS is 1.1 ± 0.8 × 10 13 L y −1 during spring 2009 and 1.9 ± 0.9 × 10 13 L y −1 during summer 2010,

Nitrogen Loads
Excess nitrogen loading stimulates phytoplankton growth and can lead to adverse ecological conditions in LIS. Nitrogen loads from wastewater effluent and the Connecticut River are typically assumed to be the dominant sources of N to LIS (NYSDEC and CTDEP, 2000;Suffolk County, 2015). However, recent work suggests that SGD may rival the N load of wastewater effluent and rivers to LIS (Tamborski et al., 2017b). The mixing zone between groundwater and seawater in the coastal aquifer, i.e., the subterranean estuary, is critically important for controlling a variety of chemical reactions, which can add or remove chemical elements from the system (Moore, 1999). In the subterranean estuary, nitrogen in groundwater can be removed by biological processes including denitrification, or groundwater nitrogen may merely be diluted by mixing with seawater (Kroeger and Charette, 2008). Conversely, nitrate can be generated as a product of organic matter remineralization, when seawater rich in organic matter infiltrates into permeable sediments from waves and tidal forcing mechanisms. Below, we provide a revised estimate of the SGD-driven N load to LIS during spring and summer conditions, using our revised Ra isotope mass balances (Figure 5).
Marine SGD includes circulating seawater flow paths driven by physical forcing mechanisms, including density-driven flow and tidal pumping (Santos et al., 2012); importantly, marine SGD is separate from terrestrial (i.e., meteoric, fresh) groundwater in the ensuing analysis. Tamborski et al. (2017b) estimated a marine SGD NO 3 − endmember of 23 ± 13 µM during spring and 37 ± 29 µM during summer. This endmember is a non-conservative N enrichment, corrected for binary mixing between seawater and terrestrial groundwaters. Stable isotope analyses of 15 N-NO 3 − and 18 O-NO 3 − suggest that the marine SGD NO 3 − is derived from the remineralization of organic matter within the subterranean estuary, rather than from an atmospheric or anthropogenic source. The marine SGD flux to LIS is estimated as the total SGD flux (totalbasin = 1.2 ± 0.9 × 10 13 L y −1 during spring 2009 and 3.8 ± 0.7 × 10 13 L y −1 during summer 2010) corrected for fresh groundwater contributions, estimated from numerical models (Scorca and Monti, 2001). As a first-order approximation, the marine SGD-driven NO 3 − flux is 2.8 ± 2.7 × 10 8 mol N y −1 for spring 2009 and 14 ± 11 × 10 8 mol N y −1 for summer 2010. Just considering the central basin of LIS, where SGD estimates are the most accurate (Figure 5), the marine SGDdriven NO 3 − flux is 2.6 ± 2.3 × 10 8 mol N y −1 for spring 2009 and 7.1 ± 6.4 × 10 8 mol N y −1 for summer 2010. We note that more work is required to fully constrain the spatial and temporal variability of the marine SGD NO 3 − endmember for the entire LIS basin. This first-order marine SGD NO 3 − flux is approximately two orders of magnitude greater than the N flux determined within the first 200 m of the Smithtown Bay shoreline by Tamborski et al. (2017b). Importantly, this suggests that SGD supplies a N load nearly equivalent to that of the Connecticut River and wastewater effluent, such that N loss via burial or denitrification may be greater than currently estimated (Vlahos et al., 2020). The mean annual NO 3 − + NO 2 − load of the Connecticut River, measured at Middle Haddam (approximately 97% of the Connecticut River drainage area) is 4.5-4.7 × 10 8 mol N y −1 , and comprises 42-49% of the total N flux of the Connecticut River (Mullaney et al., 2018). The N load from wastewater treatment plants to LIS is approximately 5.6 × 10 8 mol N y −1 (NYSDEC and CTDEP, 2000), although this load is in decline due to improving wastewater treatment conditions (Suffolk County, 2015). Importantly, these lines of evidence suggest that far more attention should be paid to monitoring SGD-driven N loads to LIS.

CONCLUSION
Physical measurements and hydrologic models often fail to capture total SGD (Burnett et al., 2006). Ra isotopes integrate over a larger spatial area, making their use to quantify SGD a popular tool among scientists. In certain large-scale embayments like LIS, the monitoring of SGD and its associated chemical load to the sea is equally as important as monitoring riverine fluxes. However, long-term SGD monitoring is seldom performed, due to the unforeseen nature of SGD and its broad difficulty in quantification. A sensitivity analysis of Ra isotope mass balances to the semi-enclosed LIS basin reveals: 1. The selection and interpretation of the Ra isotope used will ultimately depend on the target process and flow path of interest. The different ingrowth rates of the Ra quartet enable tracing different time-scale processes. Short-lived Ra isotopes may trace short-scale length processes (Santos et al., 2012), such as wave-pumping, that are not fully captured by long-lived Ra isotopes (Rodellas et al., 2017). 2. Short-lived Ra mass balances are highly sensitive to sediment (diffusion and desorption) fluxes, but are less significantly impacted by boundary mixing (scaledependent). When mixing is uncertain, 224 Ra is the preferred tracer of SGD. Studies using short-lived Ra mass balances should direct their attention toward accurately constraining sediment Ra contributions. The advantage of short-lived Ra isotopes is that their major sink is radioactive decay (scale-dependent); thus, an adequate sampling strategy can reasonably constrain the shortlived Ra inventory with minor uncertainty. Short-lived Ra mass balances can provide a reasonable first-order approximation of SGD if the Ra inventory is wellconstrained and sediment contributions are minor. 3. A large number of water column samples are necessary to accurately capture the (near) instantaneous short-lived Ra inventory. The 224 Ra inventory integrates over the half-life of 224 Ra, such that this inventory will accurately reflect the times of the year when the sampling was conducted. The 224 Ra inventory may not be representative of a seasonal or annual 224 Ra flux associated with SGD, suggesting that multiple sampling campaigns are necessary to capture any seasonality in the 224 Ra inventory for SGD flux determination. A lower number of water column samples may be adequate to capture the long-lived Ra water column inventory; however, a greater number of samples are required to evaluate spatial variability.
4. Long-lived Ra mass balances are highly sensitive to fluxes represented by exchange at the boundaries of the system. Studies using long-lived Ra mass balances in semi-enclosed environments should direct their attention toward accurately constraining mixing gains and losses. This requires characterization of long-lived Ra activities in the inflowing and outflowing water, as well as flows of water across the boundaries. Water exchange across boundaries must be well-constrained in order to minimize the final uncertainty of the Ra SGD flux. 5. Long-term SGD monitoring using short-lived Ra isotopes will require frequent water column surveys to accurately constrain the short-lived Ra inventory. Long-term SGD monitoring using long-lived Ra isotopes should focus sampling efforts on accurately constraining mixing gains and losses.
SGD to the central basin of LIS is estimated as 1.1 ± 0.8 × 10 13 L y −1 during spring 2009 and 1.9 ± 0.9 × 10 13 L y −1 during summer 2010, equivalent to 45 and 280% of the freshwater inflow of the Connecticut River during the same time period. This SGD flux supplies a bioavailable N load of 2.6-7.1 × 10 8 mol N y −1 , similar to that of the Connecticut River (Mullaney et al., 2018). Long-term SGD monitoring is required to fully understand the magnitude and temporal variability of SGD-driven N loads to LIS.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
JC, HB, and JG-O developed the seawater sampling strategy. JG-O, VR, JC, and CH were responsible for seawater sample collection. JT developed the groundwater sampling strategy and was responsible for groundwater and sediment core sample collection. JT, JC, JG-O, and CH were responsible for radium analyses. RW produced water exchange estimates from the model of Crowley (2005). JT developed and wrote the manuscript, with the assistance of JC, HB, JG-O, VR, CH, and RW. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
We thank David Bowman and the crews of the R/V Seawolf for their assistance in sampling.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fenvs. 2020.00108/full#supplementary-material TABLE S1 | Summary of Ra isotope mass balances, arranged by basin and season.