Neodymium Isotope Geochemistry of a Subterranean Estuary

Rare earth elements (REE) and Nd isotope compositions of surface and groundwaters from the Indian River Lagoon in Florida were measured to investigate the influence of submarine groundwater discharge (SGD) on these parameters in coastal waters. The Nd flux of the terrestrial component of SGD is around 0.7±0.03 μmol Nd/day per m of shoreline across the nearshore seepage face of the subterranean estuary. This translates to a terrestrial SGD Nd flux of 4±0.2 mmol/day for the entire 5,880 m long shoreline of the studied portion of the lagoon. The Nd flux from bioirrigation across the nearshore seepage face is 1±0.05 μmol Nd/day per m of shoreline, or 6±0.3 mmol/day for the entire shoreline. The combination of these two SGD fluxes is the same as the local, effective river water flux of Nd to the lagoon of 12.7±5.3 mmol/day. Using a similar approach, the marine-sourced SGD flux of Nd is 31.4±1.6 μmol Nd/day per m of shoreline, or 184±9.2 mmol/day for the investigated portion of the lagoon, which is 45 times higher than the terrestrial SGD Nd flux. Terrestrial-sourced SGD has an εNd(0) value of −5±0.42, which is similar to carbonate rocks (i.e., Ocala Limestone) from the Upper Floridan Aquifer (−5.6), but more radiogenic than the recirculated marine SGD, for which εNd(0) is −7±0.24. Marine SGD has a Nd isotope composition that is identical to the εNd(0) of Fe(III) oxide/oxyhydroxide coated sands of the surficial aquifer (−7.15±0.24 and −6.98±0.36). These secondary Fe(III) oxides/oxyhydroxides formed during subaerial weathering when sea level was substantially lower during the last glacial maximum. Subsequent flooding of these surficial sands by rising sea level followed by reductive dissolution of the Fe(III) oxide/oxyhydroxide coatings can explain the Nd isotope composition of the marine SGD component. Surficial waters of the Indian River Lagoon have an εNd(0) of −6.47±0.32, and are a mixture of terrestrial and marine SGD components, as well as the local rivers (−8.63 and −8.14). Nonetheless, the chief Nd source is marine SGD that has reacted with Fe(III) oxide/oxyhydroxide coatings on the surficial aquifer sands of the subterranean estuary.


INTRODUCTION
Previously, we hypothesized that submarine groundwater discharge (SGD) is an important, but poorly accounted for, source of rare earth elements (REE) and Nd isotopes to the ocean, and further that the SGD-flux of REEs may be sufficient to close the ocean REE budget and help resolve the Nd-paradox (Johannesson and Burdige, 2007). Several other studies have also investigated SGD fluxes of REEs to coastal waters (e.g., Duncan and Shaw, 2003;Prouty et al., 2009;Johannesson et al., 2011Johannesson et al., , 2017Kim and Kim, 2011, 2014Chevis et al., 2015a,b;Jiang et al., 2018;Paffrath et al., 2020). All these investigations have either inferred or demonstrated that SGD fluxes can be an important source of REEs to coastal waters. Hence, a growing consensus is that SGD is a fundamental component of the REE budget of the coastal ocean.
Submarine groundwater discharge includes all water flow across the sea floor and into the overlying surface waters along continental margins regardless of its composition or forcing mechanism (Church, 1996;Burnett et al., 2003;Moore, 2010). Thus, SGD can consist of terrestrial-sourced groundwaters that originate by recharging meteoric waters on the continents and/or ocean islands, as well as marine-sourced groundwater that cycles through coastal sediments and/or rocks owing to tidal pumping, wave set-up, bioirrigation, density gradients, and geothermal gradients (Taniguchi et al., 2002;Burnett et al., 2003). Numerous investigations have shown that SGD is an important flux of nutrients, carbon, and metals to the coastal ocean, and furthermore that the steep redox and salinity gradients that characterize subterranean estuaries are important "biogeochemical reactors" that can both remove solutes from solution or enhance fluxes of solutes to the ocean (Moore, 1999;Slomp and Van Cappellen, 2004;Charette and Sholkovitz, 2006;Dorsett et al., 2011;Johannesson et al., 2011;Gleeson et al., 2013;Suryaputra et al., 2015;Wang et al., 2015;Telfeyan et al., 2017Telfeyan et al., , 2018Santos et al., 2021).
The REE concentration of groundwaters in subterranean estuaries, and hence SGD fluxes of REEs to the coastal ocean, are largely controlled by biogeochemical reactions that occur when fresh, and commonly anoxic, terrestrial groundwaters of meteoric origin mix with recirculating, saline and oxic waters of marine origin, and in the process of mixing, react with solid phases that make up the subterranean estuary (e.g., minerals, sedimentary organic matter; Duncan and Shaw, 2003;Johannesson et al., 2011Johannesson et al., , 2017Kim, 2011, 2014;Chevis et al., 2015a,b;Jiang et al., 2018;Paffrath et al., 2020). For example, REE concentrations and input-normalized fractionation patterns of groundwater discharging to the Pettaquamscutt River estuary in Rhode Island appear to reflect weathering of accessory minerals like apatite that are common in local bedrock and related glacial deposits, followed by precipitation of secondary phosphate minerals during mixing of groundwaters and local stream water with incoming marine waters from Rhode Island Sound (Chevis et al., 2015a;Adebayo et al., 2020). Similarly, dissolution of basaltic glass in fractured basalt flows, as well as release of REEs from suspended particles or colloids during mixing of terrestrial-source groundwater with marine waters in subterranean estuaries appears to be an important process that mobilizes REEs into groundwaters that subsequently discharge to the coastal ocean along the arid Kona Coast of Hawaii (Johannesson et al., 2017). Other studies have also demonstrated the importance of SGD from volcanic, oceanic islands as a major source of REEs to the coastal ocean (Kim and Kim, 2011;Jiang et al., 2018;Molina-Kescher et al., 2018).
By comparison, REE cycling in the subterranean estuary beneath the Indian River Lagoon on Florida's Atlantic coast appears to be coupled to iron diagenesis that is driven, in part, by the influx of labile organic matter of marine origin (Roy et al., 2010Dorsett et al., 2011;Johannesson et al., 2011;Chevis et al., 2015b). More specifically, reductive dissolution of Fe(III) oxides/oxyhydroxides within the surficial aquifer beneath the Indian River Lagoon mobilizes Fe(II) and adsorbed REEs into the groundwater that subsequently discharges into the surface waters of the lagoon (Roy et al., 2010Johannesson et al., 2011;Chevis et al., 2015b). Submarine groundwater discharge fluxes of REEs along a nearshore transect are of the same magnitude as the effective riverine fluxes of REEs into the lagoon, whereas total SGD fluxes, which largely consist of lagoon waters that recirculate through the sediments of the subterranean estuary are roughly 15-fold larger than the local river inputs (Johannesson et al., 2011;Chevis et al., 2015b). These observations indicate that biogeochemical reactions occurring within the subterranean estuary beneath the Indian River Lagoon contribute substantially to the REE mass balance of the lagoon waters (Johannesson et al., 2011;Chevis et al., 2015b).
Despite the growing interest in SGD as a source of REEs to the ocean, to the best of our knowledge there have been no investigations of the Nd isotope geochemistry of SGD or any studies focused on the processes that control the geochemistry of Nd isotopes in subterranean estuaries. In this contribution we build on our previous investigations of SGD fluxes of REEs within the Indian River Lagoon system by presenting new Nd isotope data for local groundwaters, surface waters, and sediments from the subterranean and surface estuary. The Nd isotope data are employed to investigate possible sources of REEs to the subterranean and surface estuary (i.e., lagoon waters), and better constrain the importance of SGD to the REE mass balance within the Indian River Lagoon. Following convention, isotope ratios are expressed in epsilon notation, ε Nd (0), such that in which ( 143 Nd/ 144 Nd) Measured is the isotope ratio measured in our samples and ( 143 Nd/ 144 Nd) CHUR is the present day Nd isotope ratio (0.512638) of CHUR (i.e., Chondritic Uniform Reservoir; Jacobsen and Wasserburg, 1980).

STUDY SITE
The study site is located along the central Atlantic coast of Florida, and within the Indian River Lagoon ∼4.5 km north of the City of Melbourne, Florida, and 0.5 km south of Eau Gallie, Florida (Figure 1). The Indian River Lagoon system extends for roughly 250 km between Daytona Beach in the north and West Palm Beach in the south and includes the Mosquito and Banana River Lagoons. The Indian River Lagoon is part of a transgressive barrier island-lagoonal coastal system (e.g., Kraft and Chrzastowski, 1985) surrounded by Pleistocene beach ridges to the west on mainland Florida, as well as to the east on the bounding barrier island (Figure 1). The lagoon is microtidal with tidal amplitudes generally <10 cm and has a mean depth of 1.5 m (Smith, 1993). Consequently, the shallow depth and wind-driven waves typically ensure that the Indian River Lagoon water column is well-mixed (Smith, 1987(Smith, , 1993Martin et al., 2007). Exchange of Indian River Lagoon waters with the Atlantic Ocean occurs via three inlets (i.e., Sebastian, Ft. Pierce, and St. Lucie Inlets), all of which are located at the southern end of the lagoon. The closest of these inlets (i.e., Sebastian Inlet) is about 33 km southeast of the study site. The residence time of lagoon water near the study site is estimated to be 18 days, on average, based on computer modeling using waterlevel variation measurements; however, residence times can be as high as a year in the northern reaches of the lagoon (Smith, 1993). Inputs of freshwater to the Indian River Lagoon within the study site include two small rivers (i.e., Eau Gallie River and Crane Creek; Figure 1), urban storm runoff, direct rainfall, and SGD . Previous studies evaluated SGD fluxes at the study site using a combination of techniques including Lee-type seepage meters (Lee, 1977), chemical and thermal tracers, and models of chemical profiles within the subterranean estuary (excess 222 Rn, Ra isotopes, Cl − measurements; Cable et al., 2004Cable et al., , 2006Martin et al., 2004Martin et al., , 2006Martin et al., , 2007Smith et al., 2008a,b). Lee-type seepage meters were installed at each of the nearshore EGN-x sites (Eau Gallie North, where x is distance offshore in m), except at the EGN-12.5 location (Figure 1C), as well as the offshore, CIRL sites (central Indian River Lagoon; Figure 1B) to estimate SGD fluxes within the Indian River Lagoon . Radium isotopes, excess 222 Rn, and Cl − concentrations were employed along with steady-state 1-D advective dispersive flow models to estimate the terrestrial SGD and bioirrigation-driven SGD fluxes of water across the nearshore seepage face (i.e., EGN sites; Figure 1C) and offshore FIGURE 1 | Location of the field site along the Atlantic coast of Florida. (A) Shows the general location of the portion of the Indian River Lagoon studied as well as the location of the Canova Beach sample. (B) Documents the studied portion (blue dashed polygon) of the Indian River Lagoon and the offshore "CIRL-x" transect (yellow circles), and (C) illustrates the nearshore "EGN-x" transect in which the colored, numbered circles depict the location of each multi-level piezometers [multisamplers of Martin et al. (2003)] where groundwaters and surface lagoon waters were sampled and analyzed for REEs and Nd isotope compositions. Note, the blue circles on (C) depict samples for which only REE concentrations were determined, whereas the red circles indicate locations where REE concentrations and Nd isotope compositions were measured. (B) Also shows location where the Eau Gallie River and Crane Creek samples were collected (see Johannesson et al., 2011), the location where the Indian River Lagoon surface water Nd isotope sample was collected (i.e., CIRL 39.5), and the location of the multi-level piezometer CIRL-39 at 39, and the Lee-type seepage meters at 39, 40, 41, and 42 installed across the width of the Indian River Lagoon. Sample locations shown as red circles on (A,B) also indicate locations where we collected samples for ε Nd (0) analysis. transect (CIRL sites; Martin et al., 2007;Smith et al., 2008a;Chevis et al., 2015b). The offshore marine SGD component largely consists of lagoon waters that circulate into and out of the sediments of the subterranean estuary . Our previous studies demonstrated that the magnitude of the terrestrial SGD flux decreases with increasing distance offshore, becoming negligible at the freshwater-saltwater seepage face at about 22.5 m offshore (Figure 2; Martin et al., 2007). Exchange of porewater across the seepage face by recirculating lagoon water (i.e., bioirrigation SGD) occurs across the entire width of the lagoon and becomes an increasingly larger fraction of SGD as terrestrial SGD diminishes with distances offshore Smith et al., 2008a;Chevis et al., 2015b).
The combination of the terrestrial SGD flux and SGD resulting from recirculating lagoon waters across the lagoon sediments has been referred to as the total SGD within the Indian River Lagoon subterranean estuary Smith et al., 2008a;Johannesson et al., 2011;Chevis et al., 2015b). Specifically, terrestrial SGD mixes with recirculated lagoon waters at the nearshore seepage face and bioirrigation drives substantial SGD fluxes into the overlying lagoon across the width of the Indian River Lagoon (i.e., offshore transect; Martin et al., 2007). The total SGD flux across the Indian River Lagoon was previously estimated to be ca. 117 m 3 /day per meter of shoreline with only 0.02-0.9 m 3 /day per meter of shoreline of the total flux consisting of terrestrial SGD . Therefore, recirculated lagoon water accounts for most of the groundwater discharged to the Indian River Lagoon. Again, because of low wave and tidal action, the chief mechanism driving the recirculation of lagoon waters beyond the nearshore transect appears to be bioirrigation Martin et al., 2006). Indeed, bioirrigation depths increase from <10 cm at the shoreline to around 60-70 cm depth seaward of the nearshore seepage face (i.e., located at ca. 22.5 m from the shoreline; Smith et al., 2008a).
The Surficial Aquifer that discharges terrestrial-sourced groundwater to the Indian River Lagoon consists of the Holocene Anastasia Formation, which is composed of interbedded layers of quartz sand and coquina Martin et al., 2007;Roy et al., 2013). The aquifer has a porosity that ranges between 37 and 45% and hydraulic conductivities that range from 10 −4 to 10 −10 m s −1 in the upper 300 cm but is more uniform (10 −4 to 10 −8 m s −1 ) in the top 70 cm Martin et al., 2004;Hartl, 2006;Smith et al., 2008a,b). The Anastasia Formation is underlain by the Miocene Hawthorn Group, which consists of sand, marl, clay, carbonate, and phosphorite-rich layers (Miller, 1986(Miller, , 1997Williams and Kuniansky, 2016). The Hawthorne Group acts as the upper confining layer for the underlying Oligocene and Eocene limestones and dolostones of the Upper Floridan Aquifer System (Miller, 1986(Miller, , 1997Randazzo, 1997). Consequently, exchange of water between the Upper Floridan and Surficial Aquifer is thought to be negligible at the study site Martin et al., 2007;Roy et al., 2013).
The sediments that make up the Indian River Lagoon subterranean estuary (i.e., < 2 m thick) are composed of three separate lithostratigraphic units (Hartl, 2006;Roy et al., 2010Roy et al., , 2011Roy et al., , 2013Chevis et al., 2015b). These include black, organic matter-and Fe sulfide mineral-rich sands at the sedimentwater interface, gray colored sands with low Fe sulfide and Fe(III) oxide/oxyhydroxide contents at intermediate depths, and orange, Fe(III) oxide/oxyhydroxide-rich sediments at depth (Roy et al., 2010(Roy et al., , 2013. The thickness of the black sulfiderich sediment increases offshore, ranging from ca. 17 cm at the shoreline to ca. 68 cm at 30 m offshore. Moreover, the thickness of the black sediment layer in a sediment core collected 250 m offshore and proximal to CIRL-39 ranged from the base of the core at 230 cm below seafloor (hereafter, cmbsf) up to 45 cmbsf, where it shifted to gray colored sediments, which extended up to the sediment-water interface (e.g., see Figure 2 from Roy et al., 2011). Orange colored sands rich in Fe(III) oxides/oxyhydroxides occur at this site at depths below 250 cmbsf (Hartl, 2006;Roy et al., 2010). The increase in the thickness of the black sulfide sediments with distance offshore is thought to reflect deposition of organic matter-rich sediment with the gradual inundation of seawater into the lagoon owing to sea-level rise (Roy et al., 2010.

Sample Collection
Groundwater and surface water samples for Nd isotopic analysis were collected in April 2012 from the Indian River Lagoon subterranean estuary using a peristaltic pump with acid-cleaned Teflon R tubing attached to the tip of a drivepoint piezometer (Chevis et al., 2015a,b). Sampling was conducted in April 2012 to capture similar seasonal conditions to those experienced during our previous sampling of the subterranean estuary for REE analyses, which were conducted in April 2007 and late March/early April 2009 (Johannesson et al., 2011;Chevis et al., 2015b). The previous sampling involved collection of 64 groundwater samples, 9 surface lagoon waters samples, and two inflow streams (Cane Creek and Eau Gallie River) for major solute, Fe, Mn, S(-II), and REE analysis as described previously Johannesson et al., 2011;Roy et al., 2011;Chevis et al., 2015b).
Groundwater samples for Nd isotope analysis were collected at 5 m, 12.5 m, 17.5 m, and 30 m offshore along the nearshore transect ( Figure 1C). These samples correspond to samples obtained from multilevel piezometers (hereafter, multisamplers; Martin et al., 2003) EGN-5, EGN-10/EGN-15, EGN-17.5, and EGN-30, respectively, that we previously sampled and analyzed for REE concentrations (Chevis et al., 2015b). Because many of the multisamplers that were originally installed in 2002 (Martin et al., 2003) were no longer accessible for sampling during April 2012 owing to either having been buried or presumably destroyed by natural processes, it was necessary to use the drive-point piezometer. Because the multisamplers originally located 10 and 15 m offshore (EGN-10 and EGN-15) could not be located, we collected a groundwater sample from halfway between where these original multisamplers were located (i.e., 12.5 m from shore), which is hereafter identified as EGN-12.5. Surface waters samples for Nd isotope analysis were taken from the Eau Gallie FIGURE 2 | Highly schematic cross-section through the nearshore transect (i.e., EGN -x) portion of the Indian River Lagoon subterranean estuary. (A) Shows hypothetical flow paths for SGD in the subterranean estuary underlying the Indian River Lagoon. Note, "sed-water interface" is the sediment-water interface. (B) Illustrates the approximate locations of the multisamplers along the nearshore transect where groundwater was collected for REE analysis (light blue circles; Chevis et al., 2015b) and Nd isotope analysis (red circles, this study). Dark blue circles show the location of surface water samples analyzed for REEs (Chevis et al., 2015b). Orange squares, gray square, and black squares show location of sediment samples analyzed for Nd isotope compositions color coded orange for the Fe production zone, gray where Fe is re-adsorbed in the Fe sink zone, and black where Fe-sulfide minerals (e.g., mackinawite) precipitate in the Fe sink zone (see Roy et al., 2010Roy et al., , 2011. Numbers adjacent to each symbol for the sediments represent measured ε Nd (0) values. Cross-section is modified after Martin et al. (2007), Cable et al. (2007), and Roy et al. (2010).

River, Crane
Creek, the Indian River Lagoon (ca. 400 m offshore, CIRL-39.5; Figure 1B), and Canova Beach on the seaward side of the barrier island (Figure 1). For all Nd isotope samples, 2 L of water was filtered through a 0.45 µm in-line filter cartridge (Gelman Science, polyether sulfone membrane) and collected in acid cleaned HDPE bottles following methods described previously (Johannesson et al., , 2011(Johannesson et al., , 2017Chevis et al., 2015a,b;Adebayo et al., 2018). The samples were then acidified to pH < 2 with ultra-pure HNO 3 (Optima TM grade) and stored cold (4 • C) until analysis.

Sample Analysis
Methods employed to quantify major solutes, field parameters (e.g., pH, electrical conductivity, temperature), Fe, Mn, and REE concentrations in Indian River Lagoon waters and sediments were described previously, and the details will not be repeated here Roy et al., 2010;Johannesson et al., 2011Johannesson et al., , 2017Chevis et al., 2015b). Briefly, the analytical precision of the REE analyses was always better than 5% relative standard deviation (RSD), and generally better than 2% RSD (Chevis et al., 2015b). To check accuracy of our analyses, we processed SLEW-3 using the identical procedures employed for the Indian River Lagoon water samples (e.g., Johannesson et al., 2017;Adebayo et al., 2020). Accuracy of the analyses for SLEW-3 were typically within ± 5% of the values reported by Lawrence and Kamber (2006). Additional figures of merit are reported in the Supplementary Material.
Surface water and groundwater samples collected from the Indian River Lagoon system for Nd isotope analysis were first preconcentrated by ferric iron coprecipitation using a FeCl 3 solution (∼10 mg Fe mL −1 ) prepared with ultra-pure Fe powder and ultra-pure HCl (Optima TM grade; Jeandel, 1993;Shannon and Wood, 2005). Approximately, 2 mL of the iron solution was added to each 2 L water sample (∼10 mg of Fe per liter of sample; Jeandel, 1993) and subsequently shaken to ensure thorough mixing. Then, 50 mL of ultra-pure ammonium hydroxide (Seastar Chemicals) was added to raise the pH to 8 -9, at which point the Fe(III) flocculated and settled to the bottom of the bottles (Shannon and Wood, 2005). The majority of the supernatant was then poured off, and the Fe(III) flocculant was collected into acid-cleaned 50 mL polypropylene centrifuge tubes and centrifuged at 3,000 rpm for 5 min to separate the Fe(III) flocculant/precipitate from the remaining supernatant. The precipitate was then rinsed and centrifuged three times with Milli-Q (18.2 MΩ cm) water, after which the separated Fe(III) precipitate was dissolved in 2 M ultra-pure HCl.
Approximately 0.2 g of each sediment sample from the Indian River Lagoon was placed into a Teflon R digestion tube and mixed with 5 mL of ultra-pure HF (Optima TM grade; Tang et al., 2004). The digestion tube with a watch glass covering the top was placed in a DigiPREP MS R digestion system (SCP Science; Champlain, NY) and allowed to reflux at 90 • C for ∼2 h and then evaporated to dryness. Then, 10 mL of ultra-pure HNO 3 (Optima TM grade) was added to the digestion tube and the mixture was again evaporated to dryness. Finally, another 5 mL of ultra-pure HF was added, and the mixture was again allowed to evaporate to dryness. The remaining residue was taken up in 2 M ultra-pure HCl.
The 2 M HCl solutions for the sediments and water samples were then loaded onto Bio-Rad Econo-Pac columns packed with 8 mL of Bio-Rad R AG 50 W-X8 (200-400 mesh, hydrogen form) cation-exchange resin to separate the REEs from any remaining major salts as well as Fe and Ba. Two acid rinses of 2 M ultra-pure HCl (35 mL) and 2 M ultra-pure HNO 3 (20 mL) were performed to elute Fe and Ba, respectively (e.g., Greaves et al., 1989). Then, 15 mL of 8 M ultra-pure HNO 3 was used to elute the REEs from the column, the eluent collected into a 30 mL Savillex R beaker, and subsequently evaporated to dryness.
A second column procedure employing a polypropylene column packed with 2 mL of Eichrom R Ln-spec resin was used to separate Nd from the rest of the REEs using 0.25 M ultrapure HCl as the eluent as described previously (Pin and Zalduegui, 1997;Dry et al., 2005). The sample residue from the first REE separation step was taken up in 1 mL of 0.05 M ultrapure HNO 3 . Then, 7 mL of 0.25 M HCl was passed through the column to remove La, Ce, and some Pr from the column, followed by 10 mL of 0.25 M HCl to remove Nd with the remaining Pr. The eluent containing Nd and Pr was collected into an acid cleaned 30 mL Savillex R beaker and then evaporated to dryness. A drop of ultrapure 16 M HNO 3 was then added to convert the samples to the nitrate form and subsequently evaporated. The residue was then redissolved in 0.23 M HNO 3 in preparation for analysis.
The 143 Nd/ 144 Nd ratio was measured for the subset of selected waters and the digested sediments using a Nu Instruments Nu Plasma II multi-collector ICP-MS at Stony Brook University. Several standard runs were performed prior to sample analysis and standard runs also bracketed each sample (JNdi standard; 143 Nd/ 144 Nd ratio = 0.512115; Tanaka et al., 2000). All samples and standards were run at a concentration of 20 ppb Nd. An inrun mass bias correction was applied for all samples and standard analyses using the in-run measured 146 Nd/ 144 Nd ratio, which is corrected to the natural abundance ratio of 0.7219 using a power law relation (Wasserburg et al., 1981). A second correction was made for drift in the mass bias corrected 143 Nd/ 144 Nd observed in the standard runs. The separation methods employed in this study resulted in low yield of Nd; therefore, replicate analyses were only performed for those samples containing enough Nd.

Computing SGD Fluxes of the REEs and Nd Isotopes
We followed a similar approach to that described by Martin et al. (2007), Johannesson et al. (2011), andChevis et al. (2015b) to estimate the submarine groundwater discharge (SGD) fluxes of the REEs from the nearshore transect ( Figure 1C) and the offshore transect ( Figure 1B). Previously, we underestimated the dimensions of the study site, which necessitated a re-evaluation of these fluxes. Specifically, Johannesson et al. (2011) reported 5 km as the distance between the Eau Gallie River and Crane Creek, and width of 1.8 km across the Indian River Lagoon from the mainland to the barrier island. However, close inspection of the field site indicates that these distances are 5.88 and 3 km, respectively. Details of our approach for estimating the SGD fluxes of the REEs in the Indian River Lagoon subterranean estuary and well as information on mass balance calculations are provided in the Supplementary Material.

Neodymium Isotope Ratios
Neodymium isotope compositions of surface waters and groundwaters from the Indian River Lagoon are presented in Table 1. The ε Nd (0) values for surface waters range from −8.62 ± 0.44 in Crane Creek to −5.11 ± 0.33 at Canova Beach, whereas for groundwaters, ε Nd (0) varies from −8.23 ± 0.39 at a depth of 48 cmbsf at the EGN-12.5 site to −5.01 ± 0.42 at 50 cmbsf at the EGN-5 multisampler (Figures 1-3). Samples collected closest to shore (i.e., groundwater from EGN-5 and seawater from the surf zone at Canova Beach) exhibit the most radiogenic Nd isotope ratios and waters from the two inflow streams (i.e., Eau Gallie River and Crane Creek) are the least radiogenic ( Table 1). The surface lagoon water sample (i.e., CIRL-39.5) collected 400 m offshore has an ε Nd (0) value of−6.47 ± 0.32 (Table 1). On average, groundwaters from the Indian River Lagoon subterranean estuary are slightly more radiogenic (−6.8) than the local surface waters (−7.1). There is a general trend toward more radiogenic ε Nd (0) values with increasing Nd concentrations for all of the Indian River Lagoon water samples (Figure 3), however, the relationship is not statistically significant (e.g., R 2 ∼ 0.02). Although the Nd isotope and concentration data for the Indian River Lagoon plot along a continuum with waters from the southern Sargasso Sea and Florida Straits (Figure 3), more data are required, especially for surface waters from the coastal Atlantic Ocean, to evaluate the meaning, if any, of this feature of the data.
The Nd isotope data for sediment samples from the Indian River Lagoon subterranean estuary are given in Table 2. The ε Nd (0) values of the sediments vary from −9.22 ± 0.49 for sediments from 2 to 4 cmbsf at the EGN-22.5 multisampler to −6.98 ± 0.36 in sediments collected 205-207 cmbsf at the EGN-0 multisampler (Table 2; Figure 1C). In general, the deeper, orange-colored sediments are the most radiogenic sediments from the Indian River Lagoon subterranean estuary exhibiting an average ε Nd (0) of −7.07, whereas the shallow, black colored sediments are the least radiogenic with an average ε Nd (0) of −9.63 ( Table 2). The gray sediment sample from mid depths in the subterranean estuary has an ε Nd (0) value that is intermediate (i.e., −7.55) between the deep, orange colored sediments and the shallow, black colored sediments ( Table 2). Johannesson et al. (2011) and Chevis et al. (2015b) discussed in detail the REE concentrations and their chemical speciation in groundwaters and surface waters of the Indian River Lagoon. Here, we only present the mean concentrations of groundwaters and surface waters from the Indian River Lagoon study site normalized to the Post-Archean Australian Shale composite (PAAS; Figure 4). The mean ± 1σ concentrations in pmol kg −1 for the groundwaters from the nearshore transect (seepage face) as well as pore waters from the multisampler installed at CIRL-39 are given in Supplementary Table 5. The entire REE concentration dataset can be found in Chevis et al. (2015b). The PAAS normalized REE patterns of the sediment samples collected from the Indian River Lagoon subterranean estuary are shown in Figure 5. A general feature of all the water and sediment samples collected from the Indian River Lagoon system is their enrichments in the heavy REEs (HREE) compared to the light REEs (LREE) and middle REEs (MREE) when normalized to shale composites like PAAS (Figures 4, 5).

SGD Fluxes of REEs
Estimated SGD fluxes of the REEs across the nearshore transect ( Figure 1C) are presented in Table 3, whereas the estimated SGD fluxes of each REE for the offshore (i.e., Marine SGD) transect ( Figure 1B) are given in Table 4. In each case, fluxes are computed for the portion of the Indian River Lagoon between the Eau Gallie River and Crane Creek by multiplying the fluxes (i.e., µmol m −1 day −1 ) across the 30 m nearshore transect, and the 3,000 m offshore transect, by the 5,880 m distance between these two streams (see Supplementary Material). This approach allows us to directly compare our new estimates with those of Chevis et al. (2015b) and Johannesson et al. (2011). Because the diffusive fluxes of the REE across the nearshore seepage face are more than 1000-fold lower, on average, than any of the advective SGD fluxes, their contribution to the REE mass balance of the Indian River Lagoon is negligible, and thus, not considered further in this contribution (see also Chevis et al., 2015b).
Our best estimate for the total SGD flux of Nd across the nearshore transect seepage face is 10.1 mmol day −1 ( Table 3). This flux is the sum of the terrestrial-sourced SGD and the flux that results from bioirrigation processes in the nearshore region Chevis et al., 2015b). Based on our bootstrap analysis of the error associated with our integrated approach, we estimate that the error on this flux is ± 5% or ± 0.2 mmol day −1 (see Supplementary Material). Chevis et al. (2015b) estimated a total SGD flux of Nd of 9.4 ± 1 mmol day −1 across the nearshore seepage face, which is identical to our new estimate (Table 3). Although Chevis et al. (2015b) did not consider the offshore transect in their investigation, our previous estimate for the recirculated lagoon water component of SGD to the entire portion of the Indian River Lagoon between the Eau Gallie River, Crane Creek, and the barrier island (Figure 1) of 184 ± 33.9 mmol day −1 (Johannesson et al., 2011) is also identical to our new estimate of 184 mmol day −1 for the offshore transect (Table 4). Again, we estimate an error of ± 5% or 9.2 mmol day −1 for the computed offshore marine SGD flux of Nd (see Supplementary Material). Taken together, the total SGD across the nearshore and offshore transects amounts to around 194 ± 9.2 mmol day −1 of Nd being delivered by SGD to the studied portion of the Indian River Lagoon. We note that this value is about 70% of our previous estimate of 287 ± 64.4 mmol day −1 (Johannesson et al., 2011). As discussed below,  Chevis et al., 2015b). b Mean ± 1σ Nd and Sm concentrations of groundwater collected from the 35 and 55 cmbsf samples from multisamplers EGN-10 and EGN-15 (Chevis et al., 2015b). c Nd and Sm concentrations of groundwater from the 75 cmbsf sample from multisampler EGN-17.5 (Chevis et al., 2015b). d Nd and Sm concentrations of groundwater from the 50 cmbsf sample from multisampler EGN-30 (Chevis et al., 2015b). (0) values for Indian River Lagoon (IRL) surface and groundwaters samples plotted against the inverse Nd concentration of each sample (see Table 1). Also included for comparison are Nd isotope and concentration data for the South Sargasso Sea (Piepgras and Wasserburg, 1987) and the Florida Straits (Osborne et al., 2014).
the higher previous estimate of the total SGD Nd flux reflects the much higher Nd concentration we previously employed, which was based on a limited number of analyses from the EGN-22.5 multisampler.

REE Cycling in the Subterranean Estuary
Based on the porewater Fe and S(-II) concentrations (Supplementary Figure 6), dissolved oxygen concentrations (Supplementary Figure 7), sediment color, and the distribution of Fe(III) oxides/oxyhydroxides and total sulfur content of sediments from the Indian River Lagoon subterranean estuary, Roy et al. (2011) hypothesized an "Fe production zone" in which Fe(II) is added to porewaters by reductive dissolution of Fe(III) oxides/oxyhydroxides that coat the Anastasia Formation aquifer sands, and an "Fe sink zone" where Fe(II) is removed from porewater by re-adsorption onto aquifer mineral surfaces and by precipitation of Fe-sulfide minerals near the sediment-water interface (Roy et al., 2010 (Hartl, 2006;Roy et al., 2010). b Collected from black, Fe sulfide-rich sediments (Hartl, 2006;Roy et al., 2010). c Collected from gray/white sediments (Hartl, 2006;Roy et al., 2010).
FIGURE 4 | PAAS-normalized REE concentrations of (A) groundwaters from the nearshore transect, (B) surface waters collected along the nearshore transect, (C) local surface waters from Eau Gallie River, Crane Creek, and Canova Beach, and (D) pore water from multisampler CIRL-39 and the overlying water column at the same location [data from Chevis et al. (2015b)]. The nearshore transect groundwaters shown in (A) and the CIRL-39 pore water are the mean values of the REEs measured at each sampling port from each multisampler (see Supplementary Table 5). Mean seawater is an average open ocean seawater computed using data from Piepgras and Jacobsen (1992), Westerlund and Öhman (1992), Sholkovitz et al. (1994), German et al. (1995), Nozaki and Zhang (1995), Zhang and Nozaki (1996), and Nozaki and Alibo (2003), and is included for comparison. PAAS is the Post-Archean Australian Shale composite that is given in McLennan (1989).

7).
The iron production zone extends downward from the porewater Fe(II) concentration maxima in each multi-sampler, whereas the Fe sink zone extends upward from the porewater Fe(II) maxima to the sediment-water interface, and consists of a deeper intermediate zone where re-adsorption appears to dominate Fe(II) removal from solution, and a shallow zone   Figure 1B). The total advective SGD flux across the nearshore transect (i.e., NS J totSGD ) represents the sum of the diffusive flux ( NS J diff ), the terrestrial advective SGD flux ( NS J TSGD ), and the advective flux from recirculated Indian River Lagoon water driven by bioirrigation ( NS J bioSGD ). A bootstrap analysis of the integrated fluxes indicates errors for each REE of ± 5% for these fluxes.
directly beneath the sediment-water interface where Fe(II) is chiefly removed from solution by precipitation of Fe-sulfide minerals (Roy et al., , 2012. The Fe production and sink zones of the Indian River Lagoon subterranean estuary were subsequently described quantitatively by a one-dimensional reactive transport model that assumed upward, vertical flow of terrestrial-sourced SGD, which was able to simulate the measured porewater Fe(II) concentrations across the Indian River Lagoon subterranean estuary relatively well [see Figure 4 in Roy et al. (2011)]. The reactive transport model  Figure 1B (i.e., multisampler and seepage meter sites CIRL-39,  employed pseudo-first order rate constants to simulate Fe(III) oxide/oxyhydroxide dissolution and a generalized rate constant to account for both Fe(II) adsorption and removal from solution by Fe-sulfide mineral precipitation . More specifically, the model supports reaction of anoxic, terrestrialsourced SGD with Fe(III) oxides/oxyhydroxides that coat the deeper, orange colored sediments of the subterranean estuary (i.e., Anastasia Formation), leading to Fe(II) mobilization by microbial-driven, reductive dissolution . The mobilization of Fe(II) in the Fe production zone by microbial respiration that employs Fe(III) as the electron acceptor is further supported by the isotopically light Fe(II) that characterizes porewaters from this zone when compared to the isotopically heavier solid-phase Fe(III) oxides/oxyhydroxides (Roy et al., 2012). The mobilized Fe(II) is subsequently transported upward with advecting SGD where some is removed from solution by re-adsorption onto aquifer mineral surfaces at intermediate depths (i.e., gray colored sediments of the subterranean estuary). The majority of the mobilized Fe(II), however, is removed from solution by precipitation of Fe-sulfide minerals (e.g., mackinawite) in the shallow, black sediments directly below the sediment-water interface . It is well-known that REEs are strongly adsorbed onto Fe(III) oxides/oxyhydroxides and that adsorbed REE can be released back to solution upon reductive dissolution of these oxidized mineral phases (Koeppenkastrop and De Carlo, 1993;Ohta and Kawabe, 2001;Quinn et al., 2006;Schijf and Marshall, 2011;Schijf et al., 2015;Liu et al., 2017). Statistical analysis of REE and Fe concentration data (Johannesson et al., 2011;Chevis et al., 2015b), in conjunction with the geochemical modeling of Roy et al. (2011), are consistent with the notion that REE cycling in the Indian River Lagoon subterranean estuary is coupled to the iron cycle. Specifically, Chevis et al. (2015b) used the Bonferroni correction method to compute Spearman's rank correlation coefficients of 10,000 sample bootstraps (i.e., Dunn, 1961) for each pore water sample from all the nearshore multisamplers (64 samples; Supplementary Figure 6). The bootstrap approach is expected to be more reliable than standard approaches because it does not require any assumptions about the theoretical sample distribution (e.g., normal distributions). Despite an apparent lack of correlation between dissolved Fe and REE concentrations for the entire data set from Chevis et al. (2015b), the bootstrap method indicated that Nd is moderately correlated with Fe in porewaters from multi-samplers EGN-15, EGN-22.5, and CIRL-39, and each correlation was significant at greater than the 80% confidence level (Chevis et al., 2015b). Furthermore, correlations between Fe and Nd in porewaters from multi-samplers EGN-15 and EGN-22.5 were significant at the 95% confidence level (Chevis et al., 2015b). The bootstrapping approach also identified correlations between Fe and Yb that were significant at the 80% confidence level for porewaters from multi-samplers EGN-10 and EGN-17.5, and at the 90% confidence level for porewaters from multi-sampler EGN-30 (Chevis et al., 2015b).
These moderate but statistically significant correlations between REEs and Fe reported by Chevis et al. (2015b) for the nearshore transect likely reflects the combination of close coupling of these trace elements during reductive dissolution of Fe(III) oxides/oxyhydroxides in the Fe production zone followed by decoupled behavior in the Fe sink zone. More specifically, pore water REE concentrations are positively associated with dissolved Fe concentrations in the "Fe production zone" [i.e., zone 3 of Roy et al. (2011)] as reductive dissolution of Fe(III) oxides/oxyhydroxides mobilizes Fe(II) along with sorbed and/or co-precipitated REEs into the advecting SGD. However, in the "Fe sink zone", where Fe-sulfide minerals like mackinawite precipitate (zone 1 of Roy et al., 2011), Fe and REEs are not correlated because the REEs do not adsorb onto, or co-reprecipitate with, precipitating Fe sulfide minerals (Supplementary Figure 6). This notion is further supported by close examination of dissolved Fe, Nd, and Yb concentrations in porewaters collected above the dissolved Fe(II) maxima in multisamplers EGN-10 (i.e., <35 cmbsf), EGN-17.5 (<95 cmbsf), and EGN-22.5 (<106 cmbsf), which reveal weak relationships, none of which is statistically significant. In contrast, REEs are moderately correlated with Fe (and Mn) in porewaters collected from intervals deeper than the Fe maxima in the multisamplers, and many of these relationships are statistically significant (Chevis et al., 2015b). For example, the correlation coefficient between Nd and Fe in the Fe production zone (r = 0.55) is statistically significant at the 90% confidence level, and the correlation between Nd and Mn (r = 0.63) is significant at ca. the 95% confidence level. By comparison, the correlation between Yb and Fe (r = 0.37) is not statistically significant, but the correlation between Yb and Mn (r = 0.46) is significant at >80% confidence level.
These relationships are consistent with mobilization of REEs, and especially the LREEs and MREEs from the orange, Fe(III) oxide/oxyhydroxide coated sand during microbial facilitated reduction of these oxides/oxyhydroxides, followed by some readsorption of mobilized REEs onto the gray aquifer sediments in the deeper portion of the Fe sink zone, but no change in dissolved REE concentrations near the sediment-water interface where Fe(II) is chiefly removed by precipitation of Fe sulfide minerals (Johannesson et al., 2011;Chevis et al., 2015b). That Fe sulfide mineral precipitation plays no role in removing REEs from solution, and thus does not inhibit advective SGD transport of REEs from the subterranean estuary to the overlying surface waters of the Indian River Lagoon, is also supported by the lack of statistically significant correlations between REEs and S(-II) of the porewaters [e.g., r = 0.14 for Nd and S(-II) and r = 0.18 for Yb and S(-II); Chevis et al., 2015b]. Taken together, these observations and relationships suggest that reductive dissolution of Fe(III) oxides/oxyhydroxides in the "Fe production zone" of the nearshore transect of the Indian River Lagoon subterranean estuary is the chief source of Nd discharged by SGD to the surface lagoon waters. The presence of Fe(III) oxide/oxyhydroxide coated sands beneath sulfide rich sediments along the offshore transect and rapid, bioirrigationdriven exchange rates of porewater to depths of ∼40 cmbsf, support the notion that the same processes drive the marine SGD fluxes of Nd and the other REEs within the Indian River Lagoon.

Nd Isotopes
A remarkable feature of the Nd isotope composition of waters and sediments from the studied portion of the Indian River Lagoon is the wide range of ε Nd (0) values (Tables 1, 2; Figure 6). The wide range of ε Nd (0) values measured in Indian River Lagoon samples likely reflects a combination of sources and processes such as input from relatively unradiogenic Saharan dust and other aeolian dusts from the Atlantic region, weathering of radiogenic Eocene and Oligocene carbonate rocks deposited prior to the closure of the Central American Seaway roughly 3.5 million years ago, weathering of phosphorite-rich Miocene rocks of the Hawthorn Group, biogeochemical reactions with Fe(III) oxides/oxyhydroxides within the Surficial aquifer (i.e., Anastasia Formation), and probable anthropogenic inputs from fertilizers, localized road and construction dusts, and surface runoff of domestic water that is sourced from the Upper Floridan aquifer (Goldstein et al., 1984;Burton et al., 1997;Grousset et al., 1998;Martin and McCulloch, 1999;Aubert et al., 2001Aubert et al., , 2002Kamenov et al., 2009;Newkirk and Martin, 2009;Douglas et al., 2012;G. D. Kamenov, 2021, pers. comm.).
Stream waters discharging into the Indian River Lagoon (i.e., the Eau Gallie River and Crane Creek) generally exhibit the least radiogenic Nd isotope compositions of waters from the study site (Table 1; Figures 2, 6). For example, the mean ε Nd (0) value of these two streams is −8.39 (Table 1). Apart from groundwater collected from the EGN-12.5 site, the other groundwater and surface water samples from the Indian River Lagoon have more radiogenic Nd isotope compositions (e.g., mean = −6.4). Hence, the source of Nd to the Eau Gallie River and Crane Creek differs from that of the groundwaters and the Indian River Lagoon water column. Kamenov et al. (2009) analyzed a sample from the Anastasia Formation that has an ε Nd (0) value of −7.2; which demonstrates that weathering of the Anastasia Formation cannot solely explain the Nd isotopic composition of the Eau Gallie River and Crane Creek (Figure 6). It is important to note that the Nd isotope composition of the Anastasia Formation chiefly reflects the 143 Nd/ 144 Nd values of the ubiquitous Fe(III) oxides/oxyhydroxides that coat the quartz grains and coquina lenses that characterizes this formation (G. D. Kamenov, 2021, pers. comm.). Nevertheless, an initial assumption that the Anastasia Formation is an important source of the Nd to local surface waters and groundwaters is reasonable, especially when considered along with other potential sources. For example, Kamenov et al. (2009) determined that the ε Nd (0) value of peat-rich sediments from the Blue Cypress Marsh Conservation Area located ∼45 km south of our study site likely reflects eolian-transported Saharan dust in addition to chemical weathering of the Anastasia Formation. Consequently, we employ a simple, two end-member mixing model using the Anastasia Formation (ε Nd (0) = −7.2; Kamenov et al., 2009) and Saharan dust with an ε Nd (0) of −15 to −11 (Goldstein et al., 1984;Grousset et al., 1988Grousset et al., , 1992Abouchami et al., 1999) as end members to evaluate the possible fraction of each respective endmember that would be necessary to produce an ε Nd (0) value equal to the average Nd isotope composition of the Eau Gallie River and Crane Creek (i.e., −8.39). The model indicates that the average Nd isotopic composition of the Eau Gallie River and Crane Creek is consistent with a source that is comprised of 70-85% Anastasia Formation and 15-30% Saharan dust. It is important to note that this simple mixing model only reflects the percentage of Nd that may be contributed by these hypothesized end-member sources based on their ε Nd (0) values, and does not account for their Nd contents nor how effectively Nd is extracted from these sources by local waters. Nonetheless, the mixing fractions are identical to those calculated for the Nd isotopic composition of the peatrich sediments from the Blue Cypress Marsh (Kamenov et al., 2009). Therefore, we suggest that organic matter-rich sediments within the watersheds of both the Eau Gallie River and Crane Creek may represent an important source of Nd, and hence the REEs in general, to these streams, and further, that the Nd isotope composition of such sediments reflects chemical weathering reactions of the Anastasia Formation and atmospheric deposition of Saharan dust.
The radiogenic ε Nd (0) values of the terrestrial SGD collected from the EGN-5 multisampler and the surf zone seawater sample from Canova Beach (i.e., average of −5.04; Figure 6) is similar to the ε Nd (0) value (i.e., −5.6) of the upper Eocene Ocala Limestone reported by Kamenov et al. (2009). Moreover, these water samples also have Sm/Nd ratios that are essentially identical to the Ocala Limestone (Figure 6). The Ocala Limestone is the most permeable unit of the Upper Floridan aquifer, which also consists of the middle Eocene Avon Park Formation and the lower Oligocene Suwannee Limestone, and consequently is a major source of drinking water to large regions of Florida (Miller, 1986(Miller, , 1997Plummer and Sprinkle, 2001). Because the Ocala Limestone was deposited during the upper Eocene, and thus prior to the closing of the Central American Seaway that began in the Miocene, its Nd isotope composition reflects mixtures of radiogenic, Pacific seawater and less radiogenic Atlantic Ocean seawater (e.g., Stille et al., 1994;Burton et al., 1997;Newkirk and FIGURE 6 | ε Nd (0) values vs. Sm/Nd for Indian river Lagoon (IRL) surface waters, groundwaters, and sediments collected from the IRL subterranean estuary. Sediment samples from the subterranean estuary are color coded to match their observed colors during sample collection [e.g., see Figure 2 from Roy et al. (2011)]. The large ± 1σ range for the Sm/Nd ratio shown for groundwater sample ENG-12.5 reflects the fact that the Nd and Sm values used to compute the ratio for this sample are based on the averages of four groundwater samples from the two closest multisamplers (see Table 1). Neodymium isotope composition of Indian River Lagoon sediments analyzed in this contribution are labeled "Orange IRL Sediments" from the Fe production zone (205)(206)(207)(185)(186)(187), "Gray IRL Sediment" from the portion of the Fe sink zone where Fe is re-adsorbed to the sediments (EGN-22.5, 50-52 cmbsf), and "Black IRL Sediments" where Fe sulfide minerals precipitate within the Fe sink zone (EGN-22.5, 2-4 cmbsf, and CIRL-39, 5-7 cmbsf; Table 2). Neodymium isotope data and REE contents for the Ocala Limestone, coquina from the Anastasia Formation, and peat from Blue Cypress Marsh in Florida are from Kamenov et al. (2009). Note, the ε Nd (0) values and Sm/Nd ratio for the Anastasia Formation coquina reflects the Fe(III) oxide/oxyhydroxide coatings on these materials (G. D. Kamenov, 2021, pers. comm.). Fertilizer data are from Martin and McCulloch (1999), Aubert et al. (2002), and Douglas et al. (2012), whereas the Mn nodule (MN-8) data are from the Blake Plateau to the east of the Florida-Hatteras Slope about 285 km due east of Sapelo Island, Georgia (Piepgras et al., 1979). The data for the "Atlantic" aeolian dust and Sahara dust are from Goldstein et al. (1984) and Grousset et al. (1998), respectively. The South Sargasso Sea and Florida Straits data are the same as shown in Figure 2. Martin, 2009;Osborne et al., 2014;Montes et al., 2015;Kirillova et al., 2019). It is reasonable to expect that groundwater from the Upper Floridan aquifer would have an ε Nd (0) value similar to the Ocala Limestone. Consequently, one explanation for the radiogenic terrestrial SGD and the Canova Beach water samples is upward seepage of Upper Floridan aquifer groundwater, followed by mixing with shallow groundwaters within the Surficial aquifer (e.g., Swarzenski et al., 2001). Alternatively, because the study area is densely settled, and owing to the fact that the EGN-5 and Canova Beach samples were collected close to the shoreline (Figure 1), it is possible that their radiogenic ε Nd (0) values may reflect runoff of water from domestic uses (e.g., lawn watering, golf courses) that originated within the Upper Florida aquifer (G. D. Kamenov, 2021, pers. comm.). Finally, we cannot rule out the possibility that road dust or dust from construction projects that used stone from, or cement made with, pre-Miocene carbonate rocks such as those of the Upper Floridan aquifer or phosphorites from the Hawthorn Group may have also played a role in generating the radiogenic ε Nd (0) values of these water samples (Figure 6).
The marine SGD component and surface waters from the Indian River Lagoon have ε Nd (0) values and Sm/Nd ratios that are nearly identical to the Fe(III) oxide/oxyhydroxide coatings that characterize the Anastasia Formation (Figure 6). In addition, the orange colored sediments from the subterranean estuary have ε Nd (0) values and Sm/Nd ratios that are similar to the marine SGD and Indian River Lagoon water column samples (Figure 6). These observations suggest that reductive dissolution of Fe(III) oxides/oxyhydroxides in the Fe production zone of the subterranean estuary (e.g., Roy et al., 2011) transfers the Nd isotope signature of the Fe(III) oxide/oxyhydroxide coatings on the sediments to the marine SGD component, which subsequently is discharged to the overlying Indian River Lagoon surface waters. Again, our previous investigations, along with the SGD Nd flux estimates presented herein indicate that the lagoon water that circulates through the subterranean estuary is the chief source of SGD and Nd to the overlying lagoon waters Johannesson et al., 2011;Chevis et al., 2015b). The Nd isotope compositions reveal that reductive dissolution of Fe(III) oxides/oxyhydroxides coating the Anastasia Formation sands and coquina is the principal source of Nd that is discharged by marine SGD to the overlying lagoon waters.
Local fertilizer application may also be responsible, in part, for the radiogenic Nd isotopic composition of some Indian River Lagoon groundwaters relative to the surface waters and sediments (Tables 1, 2; Figure 6). In highly populated areas such as the Atlantic coast of Florida, including the region immediately surrounding the Indian River Lagoon, the use of fertilizers is expected to be common. Aubert et al. (2002) argued that the Nd isotopic composition of the upper soil layer and associated porewaters from the Strengbach Catchment within the Vosges Mountains (France) could not be explained solely by the preferential weathering of local minerals like apatite; therefore, an atmospheric source that included fertilizer dust was invoked (Aubert et al., 2002). More specifically, two fertilizer samples analyzed by Aubert et al. (2002) exhibited more radiogenic ε Nd (0) values [i.e., ε Nd (0) = −5.81, −6.4] than either apatite [ε Nd (0) = −6.98] or plagioclase [ε Nd (0)= −11.5], two readily weathered minerals common in the Vosges Mountains. Fertilizer samples used on fields in Australia (Martin and McCulloch, 1999;Douglas et al., 2012) and those studied by Aubert et al. (2002) from France exhibit a mean ε Nd (0) value of −6.1 (n = 14), which is more radiogenic than all the surface and groundwaters from the Indian River Lagoon except the terrestrial SGD sample and waters from the surf zone at Canova Beach (Figure 6). Nevertheless, the range of ε Nd (0) values for these fertilizer samples overlap with the ± 2σ uncertainty values for the majority of the Indian River Lagoon surface and groundwaters. The only exceptions are the less radiogenic SGD from the EGN-12.5 site and the two inflow streams, the Eau Gallie River and Crane Creek (Figure 6).
Because the extensive Miocene phosphorite deposits that occur within the Hawthorn Group in west central Florida are chiefly used to manufacture fertilizer for agriculture in the USA, it is likely that fertilizers applied to crops and lawns in Florida as well as other parts of the USA have ε Nd (0) values similar to these phosphorite deposits (e.g., Riggs, 1979a,b;Compton, 1997). The mean ± 1σ ε Nd (0) value of Florida phosphorite samples reported by Kamenov et al. (2009) is −7.22 ± 0.89, which is identical to the ε Nd (0) values of Fe(III) oxide/oxyhydroxide coatings on quartz sand and coquina from the surficial aquifer as well as the recirculated, marine component of SGD (−7.16±0.33; Figure 6). Miocene phosphorite deposits along the North Carolina shelf exhibit similar, albeit, slightly more radiogenic ε Nd (0) values (mean ± 1σ = −6.34 ± 0.46; Stille et al., 1994Stille et al., , 1996. Hence, we cannot rule out the possibility that fertilizer application within the vicinity of the Indian River Lagoon impacts the Nd isotope composition of local surface and groundwaters, although it is not clear whether the Nd isotopic and REE composition of fertilizers used in Australia and Europe are applicable to Florida. Future research should examine the REE contents and Nd isotope ratios of fertilizers used locally.

SGD Fluxes and Nd Isotopes
Our updated best estimates of the Nd fluxes to the studied portion of the Indian River Lagoon summarized from Tables 3, 4 are presented along with measured and estimated ε Nd (0) values (± 1σ) in Table 5. Approximately 95% of the total SGD flux of Nd  Johannesson et al. (2011) and Chevis et al. (2015b). g Nd flux-weighted εNd(0) value for the Eau Gallie River and Crane Creek.
reflects recirculating lagoon water that is thought to be largely driven by biorrigation Martin et al., 2006Martin et al., , 2007Smith et al., 2008a;Chevis et al., 2015b). Although our new estimate for the marine SGD flux of Nd (184 ± 9.2 mmol Nd day −1 ) is identical to our previous estimate (i.e., 184 ± 33.9 mmol Nd day −1 ; Johannesson et al., 2011), the new estimate of the total SGD flux of Nd is 1.5 times lower ( Table 5). The difference reflects the much higher Nd concentration assumed by Johannesson et al. (2011) for the total SGD flux (i.e., 485 pmol kg −1 ), which was computed from a limited number of analyses from the EGN-22.5 multisampler, which included two samples with Nd concentrations exceeding 700 pmol kg −1 (Johannesson et al., 2011). Here we used the mean Nd concentration (174 pmol kg −1 ; Chevis et al., 2015b) of porewater samples from multisamplers EGN-30 and CIRL-39, both of which are located seaward of the freshwater-saltwater interface (i.e., nearshore seepage face) that occurs ∼22.5 m offshore along the nearshore transect (Figure 2). Porewater samples from both EGN-30 and especially CIRL-39 are expected to be more representative of the marine SGD flux that discharges across the offshore transect.
The terrestrial SGD flux of Nd to the studied portion of the lagoon (i.e., 4.07 ± 0.2 mmol Nd day −1 ) is twice as high as our initial estimate (2.01 ± 0.34 mmol Nd day −1 ; Johannesson et al., 2011) but more akin to the results of Chevis et al. (2015b; 6.2 ± 0.7 mmol Nd day −1 ). In contrast, we estimate that the Nd flux owing to bioirrigation across the nearshore transect (6 ± 0.3 mmol Nd day −1 ) is nearly 2-fold larger than the estimate of Chevis et al. (2015b), which was computed assuming non-local mass transfer following the approach of Boudreau (1984) and applied by Smith et al. (2008a) to the nearshore transect. Again, our estimate of the bioirrigation flux is based on the difference between the total specific discharge and the specific discharge attributed to terrestrial SGD across the nearshore transect as determined with Lee-type seepage meters and porewater Cl − concentrations  ; Supplementary Figure 2). Although seepage meters can lead to overestimates of SGD fluxes (e.g., Shinn et al., 2002), detailed investigations of the impact of physical processes like variable wind velocity, wave height, and current speeds do not appear to affect seepage rates in the Indian River Lagoon . The combination of the terrestrial SGD Nd flux and the bioirrigation Nd flux across the nearshore transect computed here (i.e., 10.1 ± 0.2 mmol Nd day −1 ) is identical to the estimate of 9.4 ± 1 mmol Nd day −1 from Chevis et al. (2015b).
The SGD flux of Nd across the nearshore transect is thus equivalent to the effective riverine flux of Nd to the studied portion of the Indian River Lagoon ( Table 5). The effective riverine flux of Nd represents that amount delivered to the lagoon after removal of around 70% of the river borne Nd by saltinduced flocculation of colloids that occurs during mixing of Eau Gallie River and Crane Creek waters with the saline waters of the Indian River Lagoon (e.g., Goldstein and Jacobsen, 1987;Sholkovitz, 1993Sholkovitz, , 1995Sholkovitz and Szymczak, 2000;Lawrence and Kamber, 2006;Johannesson et al., 2011Johannesson et al., , 2017Rousseau et al., 2015;Adebayo et al., 2018;Costa et al., 2021). However, these Nd fluxes are all dwarfed by the marine SGD flux of Nd, which dominates REE inputs to the Indian River Lagoon (Table 5). Johannesson et al. (2011) argued that the subterranean estuary beneath the Indian River Lagoon was a net source of LREEs and MREEs to the lagoon surface waters and a sink for the HREEs. Similar observations were recently reported for a subterranean estuary on a barrier island in northern Germany (Paffrath et al., 2020). Statistical relationships between dissolved Fe and REE concentrations along with the Nd isotope values further support that biogeochemical reactions occurring within the Indian River Lagoon subterranean estuary, and in particular, microbial-driven reductive dissolution of Fe(III) oxides/oxyhydroxides, release sorbed and/or co-precipitated REEs into solution, which are then transported to the overlying lagoon water by advecting groundwater (Johannesson et al., 2011;Chevis et al., 2015b). Based on this conceptual model, we can estimate the ε Nd (0) value of the total SGD using the expression: For the subterranean estuary, Equation 2 can be solved for ε Nd (0) TotSGD giving: (3) Using the appropriate values from Table 5 for each term in Equation 3 returns an ε Nd (0) value for the total SGD component of −6.49 ± 0.46 (± 1σ). This value is identical, within the uncertainty, to the Nd isotope value of surface water collected at CIRL-39.5, which we assume represents the Indian River Lagoon water column. It is also identical to the ε Nd (0) values of porewaters from EGN-17.5 and EGN-30, which represents SGD at the seepage face, the Nd isotope composition of the deep, orange-colored Fe(III) oxide/oxyhydroxide coated sands from the EGN-0 and EGN-22.5 sites ( Table 2), as well as the Fe(III) oxide/oxyhydroxide coated sand from the Anastasia Formation (Figure 6) reported by Kamenov et al. (2009). The Indian River Lagoon surface waters, and by implication the marine SGD flux, as well as the total SGD flux are too radiogenic to reflect leaching from either the intermediate depth gray sediments or shallow depth Fe-sulfide rich black sediments near the sediment-lagoon water interface (Tables 1, 2, 5; Figures 6, 7). Consequently, our data suggest that the ε Nd (0) values of the marine SGD, and hence the total SGD within the Indian River Lagoon, reflect a source from Fe(III) oxide/oxyhydroxide coated sands within the Fe production zone (e.g., Roy et al., 2011) as Nd is mobilized during reductive dissolution to marine SGD (Roy et al., 2010Johannesson et al., 2011;Chevis et al., 2015b) and ultimately transferred to the Indian River Lagoon surface waters. This process is consistent with our previous studies that proposed this subterranean estuary is a net source of Nd to the lagoon waters (Johannesson et al., 2011;Chevis et al., 2015b) and accords well with a recent investigation of REEs in a subterranean estuary along the North Sea coast of Germany (Paffrath et al., 2020). Neodymium mass balance within the surface waters of the Indian River Lagoon can be written most simply as J Nd in = J Nd out , assuming steady state conditions. Upon expansion we have: (4) in which J Nd riv−eff , J Nd seawater , J Nd TotSGD , and J Nd MarineSGD are the effective riverine Nd flux, the flux of Nd into the lagoon from "seawater exchange", the total SGD Nd flux, and the marine SGD Nd flux due to recirculating lagoon water, respectively, [Nd] lagoon is the average Nd concentration of the Indian River Lagoon surface waters, and Q out represents the volumetric outflow of water from the Indian River Lagoon to the coastal ocean (see Supplementary Figure 8 for more details). Neodymium inputs to the Indian River Lagoon from the effective river flux, total SGD, and exchange with seawater are estimated to be 217 ± 10.7 mmol Nd day −1 , whereas outflows from the lagoon owing to recirculating lagoon waters (i.e., marine SGD) and export to the coastal ocean total is 404 ± 118 mmol Nd day −1 (Table 5). Consequently, the surface water flux of Nd from the Indian River Lagoon to the coastal ocean is estimated to be 220 ± 118 mmol Nd day −1 , which has a Nd isotope composition of −6.47 ± 0.33 (Tables 1, 5). Because the Indian River Lagoon is already saline, this relatively large Nd flux to the coastal ocean is not expected to FIGURE 7 | Schematic cross section and conceptual model through the nearshore transect and a portion of the offshore transect (i.e., > 40 m offshore) of the Indian River Lagoon subterranean estuary showing the approximate distance offshore (in meters) and depth (in cm below the seafloor, cmbsf) of water samples (filled circles) and sediment samples (filled squares) along with their corresponding ε Nd (0) ± 2σ values. Note, the horizontal scale is in meters, whereas the vertical scale is in cm. The dashed red line refers to the 300 mmol kg −1 Cl − concentration that represents the seepage face as described by Roy et al. (2010). The orange dashed line shows the approximate upper boundary of the Fe production zone, which is defined by the dissolved Fe concentration maxima across the nearshore transect, whereas the black dashed line represents the approximate upper boundary of the Fe sink zone defined by the bottom of the Fe sulfide-rich, black sediment layer within the subterranean estuary (see Roy et al., 2010Roy et al., , 2011. Curved arrows show schematic flow paths for terrestrial SGD, the ε Nd (0) value of which (-5.01 ± 0.42) is characterized by groundwater from multisampler EGN-5, and recirculated lagoon water of largely marine origin (i.e., Marine SGD) characterized by groundwaters from multisamplers EGN-17.5 (ε Nd = −6.81 ± 0.32) and EGN-30 (ε Nd = −7.16 ± 0.33). Bioirrigation is schematically depicted for the nearshore transect by the double-headed arrow.
be diminished by salt-induced flocculation and removal as occurs when fresh river waters mix with seawater. We submit that the large Nd export to the coastal ocean suggested by the simple box model (Supplementary Figure 8) reflects biogeochemical reactions occurring in the subterranean estuary, and specifically, reductive dissolution of Fe(III) oxides/oxyhydroxides and release of REEs to advecting groundwater that over time has enriched total SGD, lagoon waters, and the recirculating marine SGD in the REEs.
The closest surface Atlantic seawater sample for which a Nd isotope composition has been reported is the OCE 63 sample from Piepgras and Wasserburg (1987), which was collected more than 600 km east by southeast from the Indian River Lagoon. This water sample is from a depth of 50 m and has an ε Nd (0) value of −9.6 ± 0.9. Surface waters (i.e., 15 m) from the Bermuda Atlantic Time Study (BATS) site reported in Lambelet et al. (2016) have ε Nd (0) values that range between −9.2 and 9.13 (mean ± 2σ = −9.17 ± 0.04). These authors also report ε Nd (0) values for surface waters (≤ 75 m) of −10.34 and −9.13 (mean ± 2σ = −9.43 ± 0.86) for their two closest stations to the Florida peninsula (i.e., St. 25 located about 1,400 km east southeast, and St. 21 ∼1,640 km northeast, from the Indian River lagoon, respectively; Lambelet et al., 2016). Two shallow waters samples from the Florida Straits reported by Osborne et al. (2014) have ε Nd (0) values of −10 ± 0.75 (40 m depth) and −8.05 ± 1.13 (80 m depth), and exhibit a mean ± 2σ of −9.03 ± 1.36. Therefore, it seems likely that modern Atlantic surface waters off the coast of Florida likely have ε Nd (0) values between −11 and −8. Based on Equations 2, 4, we estimate the ε Nd (0) value of seawater that exchanges with Indian River Lagoon using the simple mass balance model shown in Supplementary Figure 8 and the data in Table 5 as: which returns a value (± 2σ) of −6.49 ± 4.5 (Table 5). Although the ± 2σ uncertainty of our estimate is large, the computed value does capture the range of Atlantic surface waters discussed above. Nevertheless, more coastal surface water Nd isotope data are required to better refine this admittedly crude estimate, and more importantly ascertain the impacts of SGD fluxes on the coastal ocean. Finally, our re-evaluated terrestrial SGD Nd fluxes are of the same magnitude of those we reported previously for the Indian River Lagoon (Johannesson et al., 2011;Chevis et al., 2015b) as well as our investigations of the Pettaquamscutt River estuary in Rhode Island (26 mmol Nd day −1 ; Chevis et al., 2015a) and the arid Kona coast of Hawaii (1-3 mmol Nd day −1 ; Johannesson et al., 2017). Furthermore, our terrestrial SGD fluxes are nearly identical to those reported for a subterranean estuary on a barrier island in northern Germany (5.03 mmol Nd day −1 ; Paffrath et al., 2020). Our new estimates of the marine SGD and total SGD fluxes of Nd to the studied portion of the Indian River Lagoon, as well as those of Paffrath et al. (2020), are orders of magnitude lower than computed SGD fluxes of Nd to the coastal ocean surrounding Jeju Island (i.e., 120 ± 60 mol Nd day −1 ), as well as to Gamak Bay (57.5 ± 11 mol Nd day −1 ) and Hampyeong Bay (219 ± 71 mol Nd day −1 ), all of which are located in South Korea Kim, 2011, 2014). Paffrath et al. (2020) noted that the substantially higher SGD flux of Nd reported by Kim (2011, 2014) may reflect interaction with easily weathered basaltic rocks [although see Johannesson et al. (2017) and Molina-Kescher et al. (2018)] compared to quartz-rich sands, the substantially larger volumetric SGD fluxes determined for Jeju Island, the much higher Nd concentrations in its coastal groundwater, and the substantially longer shorelines (up to 333 km compared to 5.88 km for the Indian River Lagoon) considered. Normalizing the SGD flux of Nd to the coastal ocean surrounding Jeju Island by the length of the island's coastline (i.e., 333.19 km; Lv et al., 2016), returns an SGD flux of 360 ± 180 µmol Nd day −1 per m of shoreline. This estimate is 11fold higher than the shoreline length normalized SGD Nd flux for the studied portion of the Indian River Lagoon (i.e., 33 ± 1.6 µmol Nd day −1 per m shoreline). Similarly, normalizing the SGD flux of Nd for Gamak Bay and Hampyeong by the length of each bay's shoreline (i.e., 230 and 130 km, respectively), returns an SGD flux of 2,500 ± 476 µmol Nd day −1 per m of shoreline for Gamak Bay, and 1,685 ± 548 µmol Nd day −1 per m of shoreline for Hampyeong Bay. These SGD Nd fluxes are 76and 51-fold higher than for the studied portion of the Indian River Lagoon. Consequently, the higher SGD Nd fluxes reported by Kim (2011, 2014) cannot solely be ascribed to the length of shoreline investigated by these researchers. Instead, the exceptionally high REE concentrations measured in coastal groundwaters from Korea also plays a major role in generating the high SGD Nd fluxes to the local coastal ocean.

CONCLUSION
The chief source of Nd to the Indian River Lagoon appears to be the result of biogeochemical reactions occurring in the underlying subterranean estuary whereby microbially facilitated reductive dissolution of ubiquitous Fe(III) oxides/oxyhydroxides that coat the Anastasia Formation sands and coquina release sorbed and/or co-precipitated REEs to the advecting submarine groundwater discharge (SGD). This conclusion is supported by the statistically significant, positive relationships between dissolved Fe and LREE and MREE concentrations groundwater samples from the subterranean estuary, as well as the identical Nd isotope compositions of the Fe(III) oxide/oxyhydroxide coatings on the subterranean estuary sands, marine SGD samples, and the overlying lagoon waters. Radiogenic terrestrial sourced SGD and waters from the surf zone at Canova Beach on the Atlantic side of the barrier island suggest either upward seepage of groundwaters from the underlying Upper Floridan aquifer that is composed of middle Eocene to lower Oligocene limestones and dolostones, and/or the phosphorite-bearing rocks of the Miocene Hawthorn, which confines the Upper Floridan aquifer. All of these rocks were deposited prior to the closing of the Central American Seaway, and consequently have radiogenic ε Nd (0) values indicative of Pacific Ocean waters. Alternatively, these radiogenic waters may reflect runoff of domestic water that is sourced from the Upper Floridan aquifer or local use of fertilizer manufactured from the Hawthorn Group phosphorites. The nonradiogenic ε Nd (0) values of waters from the two local rivers appear to reflect weathering of the Anastasia Formation and atmospheric deposition of dust sourced from the Saharan Desert.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
DC conducted REE and Nd isotope analyses and contributed to writing the manuscript. TM and NY collected groundwater and surface water samples from the Indian River Lagoon for Nd isotope analysis. JC designed and implemented the study of submarine groundwater discharge at the field site, installed seepage meters, collected and analyzed water samples for major ion composition and radium isotope and radon, determined the radium and radon contents of groundwaters, and also assisted with collection of the Nd isotope samples. ER and SH assisted with and oversaw the Nd isotope analyses of the samples and contributed to the writing the manuscript. DB helped design the study for REE and Nd isotope sampling and contributed to the mass balance modeling. JM designed and implemented the study of submarine groundwater discharge at the field site, installed seepage meters, collected and analyzed water samples for major ion composition, and contributed to writing the manuscript. CW conducted statistical computations. KJ conceptualized the original research design and wrote the bulk of the manuscript. All authors contributed to the article and approved the submitted version.

FUNDING
This study was supported by US NSF awards OCE-0825920 and OCE-1850768 to KJ, OCE-0825895 to DB, EAR-0403461 to JC and JM, as well as a St. John's River Water Management District grant SG458AA to JM and JC.