Organic matter quantity and quality across salinity gradients in conduit‐ vs. diffuse flow‐dominated subterranean estuaries

Submarine groundwater discharge (SGD) is a source of water and bioreactive solutes to coastal zones but may be modified by organic matter (OM) remineralization dynamics within subterranean estuaries (STEs). We hypothesize that bioreactive solute fluxes should depend on water residence time in STEs, but links between OM transformations and residence time in STEs are poorly characterized. To test this hypothesis, we compare dissolved OM (DOM) quantity and quality in two hydrologically distinct STE systems: a reef lagoon on the east coast of the Yucatan Peninsula, Mexico, where semidiurnal mixing in submarine springs of a carbonate karst aquifer results in short residence times, and a barrier lagoon on the east coast of Florida, where slow flow through siliciclastic sediments results in long residence times. We measured dissolved organic carbon concentrations and characterized colored DOM (CDOM) with ultraviolet spectroscopy and fluorescence combined with Parallel Factor Analysis. Both sites exhibit similar shifts in OM quality with salinity and reflect a marine source of labile OM to the STEs. Nonconservative mixing and CDOM production occurs at all sites but the long water residence times in the siliciclastic STEs cause orders of magnitude greater production than the carbonate STE. Consistent CDOM production across sites with disparate characteristics indicates that STEs are common sources of CDOM to surface water. However, observed variation in the magnitudes of CDOM production indicates that estimating global, and even regional, solute fluxes associated with SGD will be complicated by hydrologic control on extents of OM remineralization.

Subterranean estuaries (STEs) occur where fresh (or brackish) groundwater and saline pore waters mix in coastal aquifers (Moore 1999). Similar to surface estuaries, they are zones of active biogeochemical transformation of ecologically relevant solutes such as nutrients, metals, and carbon (Santos et al. 2008;Roy et al. 2010;Gonneea et al. 2014). However, longer residence times and lower water-sediment ratios in STEs than surface estuaries allow for a greater range of biogeochemical reactions to occur (Moore 1999). Many reactions are driven by remineralization of organic matter (OM), which depletes terminal electron acceptors (TEAs) and lowers redox potential (Slomp and Van Cappellen 2004). These transformations alter fluxes of bioreactive terrestrial solutes delivered to coastal zones in submarine groundwater discharge (SGD) (Slomp and Van Cappellen 2004;Windom et al. 2006;Martin et al. 2007;Roy et al. 2013;Kwon et al. 2014), complicating evaluations of global fluxes from STEs (e.g., Windom et al. 2006;Roy et al. 2010). The extent of these reactions should rely on the availability and lability of OM to drive reactions.
Biogeochemical or physical (e.g., sorption and coagulation) processes within STEs may alter dissolved OM (DOM) concentrations as freshwater and saltwater end members mix, which will impact the availability of DOM to drive biogeochemical reactions. While such nonconservative behavior has been well-documented in surface estuaries (Guo et al. 2007;Spencer et al. 2007), OM processing in STEs has received far less attention (e.g., Kim et al. 2012;Suryaputra et al. 2015). The extent of nonconservative behavior of OM in STEs will also impact OM fluxes from groundwater to surface water: decreases in DOM fluxes may occur if OM mobility is reduced via processes such as absorption or adsorption (Linkhorst et al. 2017) and would also occur due to net remineralization of DOM. In contrast, increases in DOM fluxes may occur if DOM is produced from solid-phase OM through leaching or desorption or via microbial activity (Couturier et al. 2016).
OM originates from multiple sources in STEs (e.g., marine, terrestrial, sedimentary, and microbial), and each source likely has unique reactivity characteristics and contributions to biogeochemical reactions. Specifically, terrestrial OM is generally considered to be more recalcitrant than marine OM (Aller et al. 1996;Burdige 2005), whereas in situ OM sources, such as sedimentary organic carbon or microbial cell turnover, may have qualities distict from sources external to STEs. Production or consumption of DOM within STEs will impact DOM quality in SGD and thus processes within the water column. Therefore, characterization of OM across mixing zone of STEs is important to evaluate reactivity of OM within STEs and the quantity and quality of DOM transported to surface waters via SGD.
The fraction of DOM that absorbs and re-emits light in the ultraviolet-visible (UV-Vis) wavelength range is referred to as colored DOM (CDOM). The spectroscopic properties of CDOM may reflect CDOM sources (e.g., terrestrial, marine, and microbial), the extent of microbial processing or photodegradation, or molecular properties such as aromaticity and molecular weight (Table 1). While spectroscopic characterization of terrestrial compared to marine CDOM has been the focus of many studies in surface estuaries (e.g., Zepp et al. 2004;Guo et al. 2007;Kowalczuk et al. 2010;Cawley et al. 2012), to our knowledge, few studies have used spectroscopic techniques to evaluate CDOM transformations in STEs. We use these techniques to asses shifts in the sources of OM with salinity, as well as to make inferences about microbially driven biogeochemical processes and changes in molecular properies of CDOM in STEs. Our results provide direct assessments of quantity and quality of OM in SGD as well as indirect assessments of solutes, such as nutrients, metals, and greenhouse gases, that may be impacted by OM reactions.
Compared to surface estuaries, sharp and relatively immobile redox interfaces may exist in STEs and impose additional controls on biogeochemical dynamics. These interfaces result from gradients of TEA availability. For example, surface salt water contains reatively high concentrations of oxygen (Young et al. 2018) and sulfate, whereas fresh groundwater may contain high concentrations of nitrate, but little dissolved oxygen or sulfate (Kroeger et al. 2007;Kroeger and Charette 2008). OM remineralization mostly results from microbial processing; however, microbial metabolism of each TEA has a unique energy yield and leads to distinct reaction kinetics (Froelich et al. 1979). Therefore, the distribution of TEAs will also play a critical role in OM processing in addition to DOM reactivity.
Although OM reactivity and TEA availability control biogeochemical OM processing within STEs, the water residence times in STEs vary depending on hydraulic conductivity, head gradients, and thus groundwater flow rates. Although the extent of STE biogeochemical processes relative to residence times is poorly characterized, we hypothesize that long residence times will increase extents of biogeochemical reactions among water, solids, and microbes. We test this hypothesis through characterization of CDOM in STEs representing two

Pain et al.
Organic matter processes in STEs hydrogeologic end members with disparate residence times: one with slow, widely distributed flow through siliciclastic sediments (Indian River Lagoon, Florida; e.g., Martin et al. 2007) and another with rapid, point-discharge through a carbonate karst aquifer (the Yucatan Peninsula, Mexico; Swarzenski et al. 2001;Null et al. 2014;Young et al. 2018).
We first assess changes in the quantity and quality of OM along salinity gradients and use salinity-based conservative mixing models (CMMs) to evaluate CDOM dynamics in the two STE end members. In addition, we sample three distict locations within the slow distributed flow STEs at the Florida site, each of which shows unique OM processes regardless of similarities in their hydrologic characteristics. We use these findings to discuss implications for controls on the magnitude of biogeochemical processes in STEs as well as fluxes of DOM to coastal zones from STEs of differing residence times.

Study locations
Our study sites are located in tropical (Yucatan) and subtropical (Indian River Lagoon) lagoonal settings (Fig. 1a). The Indian River Lagoon site is located on the east coast of Florida in a barrier island lagoon complex, which includes Indian River, Banana River, and Mosquito lagoons (Fig. 1b). The Yucatan site is located on the Yucatan Peninsula in a reef lagoon offshore of Puerto Morelos, Quintana Roo, Mexico (Fig. 1c,d). The two sites are separated by only a few hundred kilometers and thus have similar climates, including cool dry winters with warm humid summers and periodic impacts from tropical storms. Both sites are microtidal and have low wave energy as a result of barriers separating them from the ocean.
Indian River Lagoon spans approximately 250 km of coastline ( Fig. 1a,b). We collected samples at three different seepage faces discharging to the lagoon: Eau Gallie North (EGN), Riverwalk Park (RWP), and Banana River Lagoon (BRL; Fig. 1b). Indian River Lagoon STEs contain siliciclastic sediments that range from fine sand to clays. This fine grain size limits discharge rates to 0.02 to 0.9 m 3 d −1 m −1 of shoreline at EGN. (Table 2; Martin et al. 2007). Flow rates at RWP and BRL have not been measured but are likely of the same order of magnitude. This slow seepage makes STE salinity gradients static over time scales of days to weeks in contrast to the rapid exchange of the Yucatan aquifer. However, seasonal variation in lagoon water salinity and fresh groundwater head causes fluctuations in seepage face width (Roy et al. 2013), and storm-driven saltwater intrusion events can alter seepage face salinity for several months (Smith et al. 2008).
The Yucatan Peninsula is a karstic carbonate platform of Triassic to Holocene age. The platform contains extensive secondary porosity, which increases aquifer hydraulic conductivity (Table 2). Aquifer hydrological properties vary substantially among matrix material, fractures, and conduits. While 97.1% of storage occurs in the matrix, 99.8% of flow occurs in conduits and fractures due to the orders of magnitude higher hydraulic conductivity (Worthington et al. 2000; Table 2). Elevated hydraulic conductivity limits surface water to cenotes (water-filled sinkholes) as precipitation rapidly infiltrates the thin soil layers overlying the carbonate matrix to recharge the aquifer. Consequently, all freshwater drainage to the coast is submarine. Submarine springs provide~78.5% of the discharge at the Puerto Morelos lagoon while diffuse seepage from the beach shore face provides the remaining discharge (Beddows et al. 2007;Null et al. 2014). Submarine springs can reverse flow during extreme high tides, storm surges, wind setup or weakened Yucatan current, allowing surface seawater to intrude into conduits (Parra et al. 2015;Young et al. 2018). These events may occur daily and persist for several hours if tidally driven (Young et al. 2018) or may persist for several days if due to storm surges or wind setup (Parra et al. 2015). Intrusion events create a brackish mixing zone between freshwater and salt water that extends approximately 1-4 km inland (Beddows et al. 2007). Four springs (Hol Kokol, Gorgos, Laja, and Pargos) were sampled for this study. Although distributed along~5 km of shoreline, they are considered part of a single STE (Null et al. 2014;Fig. 1c,d).
The distinct hydrogeological characteristics between our two sampled regions (Table 2) lead to large estimated differences in residence times and the distance within the STE over which mixing occurs. Slow flow through Indian River Lagoon STEs (Martin et al. 2004) results in freshwater-salt water mixing zones that range between < 10 cm to > 1 m. These slow flow rates would result in residence times on the order of 1.5-28 d ( Table 2). However, shorter residence times of 0.33 d may occur in the upper sediment layers of Indian River Lagoon (< 70 cm) because of bioirrigation (Martin et al. 2006; Table 2). In contrast, rapid groundwater flow in Yucatan conduits generates turbulent mixing zones that may extend across 100s of meters within the STE (Parra et al. 2015;Young et al. 2018; Table 2). Despite the wide mixing zone compared to Indian River Lagoon, groundwater flow rates at the outlets of submarine springs have been measured at 0.1-0.5 m s −1 (Parra et al. 2015), implying that water residence times may range from < 1 to 20 min within the freshwater-salt water mixing zone.

Sample collection
Indian River Lagoon samples were collected four times starting in fall 2014 and ending spring 2016, at all three STE sites. The fall and spring sampling occurred during September-October and May, respectively, to evaluate seasonal variations in OM processing. STE widths varied between sites and extended about 20 m offshore at EGN, > 35 m offshore at RWP, and > 45 m offshore at BRL. Permanent multilevel piezometers (multisamplers; Martin et al. 2003) were installed perpendicular to the shoreline across the STE. Piezometers are constructed with multiple (4-8) well screenings at depths ranging from 7 cm to 2.5 m below the sediment-water interface (Fig. 2a). Tubing leads from the screened intervals to the surface and is sampled by pumping water using a peristaltic pump. Multisamplers were installed in 2004 at EGN and between May 2014-September 2015 at RWP and BRL. Wells are located 0, 10, 20, and 22.5 m offshore at EGN, 10, 20, and 35 m offshore at RWP, and 1, 11, 21, and 45 m offshore at BRL.
At the Yucatan site, water was sampled over a 2-week period in September 2014. During the sampling period, surface lagoon water periodically intruded into the springs. These intrusions caused salinity of the discharge to vary, which we use to trace salt water-freshwater mixing (Young et al. 2018) Table 2. Hydrologic properties of carbonate and siliciclastic STE in this study. Total SGD represents combined terrestrial and marine SGD. Estimates of STE residence time are based on literature reported values of specific discharge and mixing zone width. Residence times are calculated by multiplying the reported range in specific discharge to the reported range in mixing zone widths.

Indian River Lagoon Yucatan
Type Siliciclastic Carbonate karst Hydraulic conductivity (m s −1 ) Mixing zone width (m) 0.1-1 m 10-100* , † Residence time~1.5 to 28 d~20 s to 17 min *Zimmermann et al. (1985). and evaluate varations in OM processing related to mixing. Samples were collected from four springs ( Fig. 1d) during discharge periods only. Samples were collected by pumping water to the surface through 0.5 cm diameter flexible poly(vinyl chloride) tubing. At Indian River Lagoon, sample tubing was connected to multisampler piezometer ports (Fig. 2a). In the Yucatan, tubing was installed by SCUBA diving inside the conduit, < 1 m from the spring openings ( Fig. 2b). At both locations, a YSI Pro-Plus sensor was installed in an overflow cup in-line with the tubing to measure salinity, temperature, pH, dissolved oxygen, and oxidation-reduction potential (ORP) while pumping water. Once these parameters were stable, water samples were filtered through 0.45 μm trace-metal grade Geotech medium capacity disposable canister filters and collected and preserved in the field using methods appropriate for each solute. Samples for dissolved organic carbon (DOC) concentrations and CDOM analysis were collected in amber borosilicate vials that were combusted at 550 C prior to use. We restrict the analysis of OC to the < 0.45 μm size fraction. Samples therefore include dissolved and colloidal size fractions, although we refer to the < 0.45 μm size fraction as DOM. Samples for DOC concentration measurements were acidified with hydrochloric acid to pH < 2. CDOM samples were not acidified and kept frozen until analysis within 1 month of collection.

Laboratory methods
We characterize OM quality through the use of spectroscopic methods that provide both qualitative and quantitative information through UV-Vis and fluorescence spectroscopy (Table 1). While CDOM is a subset of total DOC, DOC and CDOM are commonly collinear in coastal systems (Del Castillo 2005). DOC concentrations were analyzed on a Shimadzu TOC-VCSN total organic carbon analyzer (Shimadzu Scientific Instruments, Inc.), and the coefficient of variance for check standards was less than 2%. Fluorescence measurements were collected on a Hitachi F-7000 Fluorescence Spectrophotometer (Hitachi, Ltd.) to generate three-dimensional excitation-emission matrices (EEMs). Scans were collected at 700 V and at excitation wavelengths ranging from 240 to 450 nm at 5-nm intervals and emission wavelengths ranging from 250 to 550 nm at 2-nm intervals. Instrumentspecific effects were corrected for differences in lamp intensity across the excitation-emission wavelength range. Inner filter effects due to organic carbon content were corrected with UV spectra according to methods outlined in Ohno (2002). Aliquots of samples collected for CDOM measurements were used to measure UV absorption at 1 nm intervals from 240 to 550 nm on a Shimadzu 1800 UV Spectrophotometer (Shimadzu Scientific Instruments, Inc.).

Statistical modeling
The EEMs were deconvolved into statistically resolvable components through the use of Parallel Factor Analysis (PARAFAC) (Murphy et al. 2013) using the drEEM toolbox in Matlab version 2015b (Murphy et al. 2013;MathWorks ). Preprocessing of EEMs included applying corrections for instrument-specific effects, innerfilter effects, masking to eliminate signals from firstand second-order Rayleigh scattering, and conversion to Raman Units (R.U.) (Ohno 2002;Zepp et al. 2004). To reduce collinearity between samples due to dilution effects and to better model low-concentration samples, EEMs were normalized to total sample fluorescence intensity before modeling with PARAFAC. A total of 322 samples were included in the PARAFAC model, which was run with nonnegativity constraints and split-half validated. Model results were reverse-normalized before being exported. The abundances of PARAFAC components are reported in R.U. that are normalized to the fluorescence intensity of water to remove variability due to instrument drift (Lawaetz and Stedmon 2009). While R.U. are quantitative, they cannot be converted to molar units as relationships are unknown between concentration and fluorescence intensity for each PARAFAC component. Details of these relationships would require an instrument-specific calibration for each component. The subterranean estuary, or mixing zone between freshwater and saltwater end members, occurs across the entire freshwater-salt water interface. Permanent piezometers (vertical gray lines) installed in sediment are used to sample across the seepage face. Black circles represent pore-water sampling locations on piezometers. Modified from Martin et al. (2007). (b) Yucatan model depicting a carbonate karst STE with point discharge through conduits. Mixing occurs within conduits that have much higher hydraulic conductivity than the carbonate matrix. Samples were collected < 1 m from the discharging spring vent.
Component abundance is discussed in terms of their relative concentrations as well as quantity in R.U. Relative abundances are proportional to total fluorescence (%C n ), where n represents the component number, calculated as component abundance (in R.U.) divided by total CDOM. Total CDOM is defined as the sum of all PARAFAC components (Total CDOM = P n 1 C n ) to allow direct comparison between the quantity of CDOM and the relative proportion of each PARAFAC component.

Mixing models
Changes in DOC and CDOM concentrations caused by reactions within the STEs are identified based on deviations from concentrations expected from conservative mixing between freshwater and saltwater end member values. The freshwater end member compositions for the Indian River Lagoon STEs are taken as the average of the freshest porewater sample collected at each STE site at each sampling time, and the saltwater end member composition is taken as the average of surface lagoon water collected at each STE site at each sampling time. Because only one sampling time is available for the Yucatan, Yucatan STE freshwater and saltwater end members are assumed to be the freshest discharging groundwater collected and surface seawater outside the influence of the spring plume (blue star in Fig. 1d), respectively. In most cases, the surface seawater sample has the highest salinity in the dataset. However, some pore water at EGN has higher salinity than surface water because surface and pore water exchange can be slower than salinity changes in surface water resulting from evaporation and precipitation (Martin et al. 2006).
Changes in DOC and CDOM abundances due to reactions (DOC Rx and CDOM Rx ) are expressed both as absolute deviation from CMM value (measured in mg L −1 for DOC and R.U. for CDOM), as well as percent deviation from CMM:

PARAFAC components and salinity relationships
PARAFAC modeling resulted in a five-component model that explained > 99% of the variability among the EEMs wavelengths. Colors scale from dark blue (low) to yellow (high) and indicate relative fluorescence intensity exhibited by components. Components C1, C2, and C4 are characterized as terrestrial humic like, C3 is microbial humic like, and C5 is protein like (Table 3). (Fig. 3). We compared our modeled components to those identified in the literature using OpenFluor (http://www. openfluor.org/), where matches are identified when the correlation between excitation-emission spectra of our modeled components to literature-reported components yields r 2 > 0.95 (Table 3). All identified components have statistically significant positive correlations (p < 0.001) with each other; however, correlations have variable strengths with r 2 values ranging from 0.13 to 0.95 (Table 4).

Organic carbon concentrations and CMM results
Total CDOM shows a positive correlation to total DOC for all sites but each site has a unique slope (Fig. 5). Both DOC and CDOM are correlated (p < 0.001) to salinity at the Yucatan, EGN, and RWP sites, but not at the BRL site (Fig. 6). Correlations with salinity are negative at the Yucatan for both DOC (r 2 = 0.92) and CDOM (r 2 = 0.96) and EGN (DOC r 2 = 0.47; CDOM r 2 = 0.13). Salinity is weakly negatively correlated with DOC (r 2 = 0.13, p < 0.001) and CDOM (r 2 = 0.10, p < 0.001) at RWP. Visual inspection of the distribution of DOC and CDOM with salinity of Indian River Lagoon samples indicated no clear seasonal differences. We therefore do not display Indian River Lagoon data per sampling season and instead examine changes in DOC and CDOM across the salinity gradient as a whole.
Results of CMMs for DOC and CDOM indicate smaller deviations between measured and CMM predictions at the Yucatan compared to Indian River Lagoon sites (Table 5; Fig. 7). DOC residuals are predominantly negative for the Yucatan site but CDOM residuals are positive (Fig. 7a). For Indian River Lagoon sites, DOC residuals are predominantly positive for BRL and RWP sites but are both positive and negative at EGN (Fig. 7b-d). CDOM residuals are predominantly positive for all Indian River Lagoon sites. DOC and CDOM exhibit maximum deviation from conservative mixing at a higher salinity (~20) at the Yucatan, and lower salinity (~5) at BRL. At EGN and RWP, the deviation from conservative mixing occurs across the range of salinity (0-25).

UV indices
The relationships between salinity and UV indices of CDOM vary between all STEs (Fig. 8). SUVA-254 values of saltwater end members are lower than freshwater end members for the Yucatan and RWP sites, but SUVA-254 is similar between the two end members at BRL and EGN. S R and S 275-295 are higher for seawater samples than freshwater   (Fig. 8). S 350-400 is greater in saltwater end members compared to freshwater end members for BRL and EGN, but values are approximately the same at RWP and the freshwater end member is greater than the saltwater end member for the Yucatan. Water samples with salinities intermediate to freshwater and saltwater end members display variable relationships between salinity and UV indices. SUVA-254 and S R values increase above freshwater and saltwater end member values (Fig. 8). Particularly for Yucatan, BRL, and EGN sites, SUVA-254 and S R reach maximum values at approximately the same salinity interval with maximum CDOM and DOC residuals from conservative mixing lines (Figs. 7-8), while S 275-295 and S 350-400 decrease in these salinity intervals.

Discussion
The following discussion describes differences in OM quantity and quality of the sampled STEs through analyses of the deviations expected from conservative mixing. We evaluate end member sources for OM in STEs, the OM qualities based on spectroscopic characteristics, and shifts in the quantity and quality of OM with salinity. We discuss these characteristics with respect to impact of residence time on varations in OM processing in STEs and implications for fluxes of OM in SGD.

Mixing model end members
Mixing models have been used to study chemical transformations in STEs (e.g., Santos et al. 2008) but require the definition of appropriate freshwater and saltwater end members. We assume a two end member model and consider that STE pore waters are initially derived from mixing of fresh groundwater and surface saltwater that enters the STE via processes such as bioirrigation, tidal pumping, or diffusion. The chemical composition of water samples is determined by the relative proportions of freshwater and saltwater end members but may be modified by reactions between dissolved species or by interaction with solid STE materials. However, because our mixing model end members describe the origin of pore waters rather than OM, we do not consider the addition of external OM sources as end members but as processes that contribute or remove OM and alter its quality in the STE.
Freshwater end member compositions may vary due to heterogeneity in aquifer compositions, groundwater flow paths, and perhaps most importantly from variable catchment land cover. In Indian River Lagoon, fresh groundwater exhibits considerable chemical variability among the sampled STEs, and catchment land cover appears to be the most important control. For instance, EGN is directly offshore of a heavily developed area with a high proportion of paved surfaces, which may limit groundwater infiltration, while RWP is directly offshore of the last remaining natural wetland on the Indian River Lagoon, which may enhance both groundwater infiltration and its organic carbon content. BRL is located at a shoreline colonized by mangroves, which may increase DOC concentrations and drive biogeochemical reactions from root respiration. These variations are reflected in OM content: DOC concentrations in the freshwater end member at EGN is 7 times lower than at RWP and 10 times lower than at BRL (Table 5). Saltwater end member composition varies less at the Indian River Lagoon STEs than fresh groundwater because of relatively fast flushing of the lagoon which ranges between approximately 6 and 43 d (Kim 2003). However, despite the rapid flushing, some salinity variations occur, possibly from evaporation (e.g., Martin et al. 2006) and/or from localized freshwater input from surface streams or SGD. Because of variable compositions of both saltwater and freshwater end member compositions, we use unique concentrations for each Indian River Lagoon STE.
In the Yucatan, higher hydraulic conductivity and more rapid groundwater flow than Indian River Lagoon (Table 2) leads to similar freshwater composition along the coastline, which has led to the recognition that offshore springs represent discharge from a single STE system (Null et al. 2014). This single system permits multiple submarine springs to be appropriately modeled by a common freshwater end member. Salt water in the Yucatan lagoon has an average residence time of 3 h but can be as short as 0.35 h following storm swells (Coronado et al. 2007). The saltwater end member composition for all sampled springs is thus best represented by seawater collected outside the direct influence of discharging submarine springs (Fig. 1).

Organic carbon quality across salinity gradients
Using the above-defined end members, salinity-based mixing models indicate nonconservative behavior in STEs: CDOM is produced at all STEs despite variability in end member salinity and OM concentrations (Table 5), while DOC is both produced (BRL and RWP at Indian River Lagoon) and consumed (EGN at Indian River Lagoon and Yucatan; Fig. 7). We examine changes in the spectroscopic qualities of CDOM between sites to better depict the shift in OM sources and quality with increasing salinity and to characterize the addition of CDOM within STEs that is reflected in nonconservative mixing. These results have implications for the reactivity of OM, which is a parameter that strongly impacts the rates and magnitudes of many other biogeochemical processes and fluxes of DOM due to SGD.

Freshwater OM quality
The majority of OM in STEs is derived from fresh groundwater, as indicated by negative correlations between salinity and DOC and CDOM at BRL, RWP, and the Yucatan sites (Fig. 6) and identification of three terrestrial humic components in PARAFAC modeling (Table 3). UV spectroscopic indices suggest that freshwater CDOM is more aromatic than marine CDOM, as evidenced by elevated SUVA-254 values (particularly for the Yucatan and RWP sites; Fig. 8a,d), and may result from CDOM produced through the degradation of aromatic terrestrial soil OM (e.g., Mcknight et al. 2001). Lower S 275-295 values also occur in freshwater compared to saltwater end members, which may be explained by increases in the  proportion of terrestrial OM, as was observed in Fichot and Benner (2012). The exception to this trend occurs at RWP, which has a freshwater S 275-295 value similar to the saltwater end member. At this site, however, S 275-295 has a positive correlation with salinity (p < 0.001), which also suggests the freshwater contains relatively more terrestrial OM than salt water (Fig. 8d). Decreases in S 275-295 have also been shown to reflect increasing molecular weight, which suggests that terrestrial OM is of greater molecular weight than marine OM, while higher SUVA-254 values suggest it is also relatively more aromatic. Both degree of aromaticity and increased molecular weight are associated with increased recalcitrance of OM (Arndt et al. 2013), suggesting that terrestrial OM delivered by freshwater may be relatively less labile than marine OM.

Saltwater OM quality
Compared to freshwater, salt water has a relatively greater abundance of the protein-like PARAFAC component (C5) at all STE sites (Fig. 4). Proteins tend to be more labile than other OM types such as humic or fulvic acids (Arndt et al. 2013) and OM reactivity correlates positively with protein-like OM, measured with fluorescence or direct quantification (Fellman et al. 2009;Cory and Kaplan, 2012). These associations suggest C5 may represent a more labile fraction of CDOM in STE settings. UV indices also suggest that marine OM may be of higher reactivity than terrestrial OM. For instance, higher S R values of marine OM may reflect higher degrees of photodegradation in surface salt water, which may lead to the production of labile molecules that are readily microbially degraded (Helms et al. 2008). These findings are consistent with increasing S R values with salinity in surface estuaries, which result from more photodegradation in seawater compared to turbid river water (Helms et al. 2008).
While salt water appears to be a source of reactive OM, its delivery to STEs requires mixing of salt water into the subsurface. Delivery of marine OM to the STE may occur by advection of marine surface water across the sediment-water interface through processes such as bioirrigation, tidal pumping, and wave setup (Moore 1999;Martin et al. 2006) The addition of reactive OM to all our sampled STEs occurs regardless of flow rates or residence time of water in the STEs. Consequently, the presence of marine OM in STEs appears to be a more critical factor to biogeochemical reactions than flow rates, although the location of those reactions reflects the delivery rate of the labile OM.

Shifts in OC quality with nonconservative mixing
Changes in OM spectroscopic properties occur where CDOM deviates from conservative mixing in the freshwater-salt water mixing zone, indicating that its quality is distinct from either freshwater or saltwater end members and may be derived from sedimentary OM or in situ microbial biomass. Where CDOM and DOC residuals are highest (Fig. 7), CDOM is more aromatic, as reflected by local maxima of SUVA-254 values, and may be more microbially processed, as reflected by local minima of S R , S 275-295 , and S 350-400 (Fig. 8). In situ production of CDOM coincided with increases in SUVA-254 and S R in other STE settings and was interpreted to reflect the production of aromatic, low-molecular-weight compounds (Couturier et al. 2016). Higher aromaticity and lower molecular weight suggest that produced CDOM may be less reactive than that contained in freshwater or saltwater end members (e.g., size-reactivity continuum model; Benner and Amon 2014). Decreases in both S 275-295 and S 350-400 that coincide with increased SUVA-254 and S R values, particularly at BRL and EGN (Fig. 8b,d), could reflect microbial processing or increased molecular weight (Table 1). Because higher S R values would indicate lower molecular weight of produced CDOM, decreases in S 275-295 and S 350-400 are more likely due to greater microbial processing of OM.

Link between hydrogeology and nonconservative mixing of OM
Mixing models of DOC and CDOM are distinct between Indian River Lagoon and Yucatan STEs: in the karstic carbonate aquifer of the Yucatan, DOC and CDOM concentrations are close to those predicted by conservative mixing, similar to the findings of Beck et al. (2007) and Kim et al. (2012), whereas Indian River Lagoon STEs have larger residuals from conservative mixing in both absolute (mg L −1 DOC and R.U. CDOM) and relative magnitudes (%) and residuals indicate CDOM production (Fig. 7). This behavior contrasts to surface estuaries, where CDOM is more typically removed by physicochemical processes such as coagulation and flocculation resulting from salinity changes (Cauwet 2002). These results suggest that CDOM production is a common trait of biogeochemical processing in STEs: for instance, CDOM production has also been observed in Gulf of Mexico STEs and interpreted to reflect contributions to the DOM pool by microbial degradation of particulate OM (Santos et al. 2009;Suryaputra et al. 2015). Our results additionally suggest that the magnitude of DOM production depends on hydrogeologic characteristics and flow rates with the STEs, which strongly differ between siliciclastic and carbonate karst aquifers.
Although the degree of nonconservative mixing strongly contrasts between Indian River Lagoon and Yucatan sites, considerable variability is observed among the Indian River Lagoon STEs and suggests that additional parameters besides residence time impact the magnitude of biogeochemical reactions that occur. These parameters may include heterogeneous sediment composition, including its OM content, as well as differences in OM quality of the inflowing fresh groundwater. Both of these variables may impact the relative proportion of organic carbon compared to TEA concentrations in fresh groundwater. This then impacts how OM is processed in STEs as fresh groundwater mixes with salt water that contains labile organic carbon and potentially more favorable TEAs such as oxygen and sulfate. Chemical variations in fresh groundwater composition could impact OM processing even in hydrogeologically similar STE settings in Indian River Lagoon and consequently, we explore below linkages between TEA delivery and OM processing in the Yucatan and Indian River Lagoon STEs.

Link between hydrogeology and TEA delivery
While CDOM is produced at all sites, the maximum deviation from conservative mixing occurs at distinct salinities for each STE. Specifically, CDOM production is greatest in fresher portions of BRL (Fig. 7b), at mid-salinity regions of the Yucatan (Fig. 7a), and in both fresh and saline portions of EGN and RWP (Fig. 7c,d). These salinity-CDOM production relationships may reflect availability of TEAs and the resulting reactions. Oxygen-and sulfate-rich surface salt water should inhibit methanogenesis, which is the terminal metabolic process in most reduced fresh groundwater (Froelich et al. 1979). The freshest groundwater in this study is reducing, with the exception of water with salinity < 1 at EGN (Fig. 6), and thus mixing with oxygen-and sulfate-rich saltwater mixing should enhance OM remineralization in the fresh groundwater. At EGN, where fresh groundwater has positive ORP values ( Fig. 6c), the largest residuals from conservative mixing occur at higher salinities, suggesting delivery of new TEAs such as sulfate or oxygen may not increase OM remineralization.
The rate of salt water delivery to the STEs may influence OM remineralization rates because redox reaction kinetics decrease from oxygen to less energetically favorable TEAs and TEA delivery depends on flow rate. Fastest delivery occurs at the Yucatan STE during seawater intrusion into conduits at high tide and during storm setup and results in rapid consumption of OM through biogeochemical reactions (Young et al. 2018). In contrast, slow flow rates in the Indian River Lagoon STEs result in complete consumption of oxygen within a few centimeters of the sediment-water interface, leading to negative ORP values of the saline portion of the STE (Fig. 6). Two of the three sites at Indian River Lagoon (BRL and RWP) have freshwater with low redox potential (and thus limited TEA availability) as well as greater DOC and CDOM concentrations in freshwater than salt water (Fig. 6), suggesting TEA limitation of OM remineralization. The introduction of sulfate-rich seawater into the freshwater portion of Indian River Lagoon STEs would thus enhance OM processing as reflected in maximum CDOM residuals at the freshwater end of the salinity gradient at BRL (Fig. 6b,d).
Alternatively, if organic carbon rather than TEA availability is limiting, enhanced OM remineralization may occur from delivery of reactive TEAs (likely sulfate) in salt water allowing more extensive processing of sedimentary OM reservoirs, as was observed in Suryaputra et al. (2015). OM limitation may occur at EGN, where freshwater contains lower DOC and CDOM concentrations than salt water and has positive ORP values, reflecting elevated TEA concentrations (Fig. 6c). Consequently, CDOM production may occur in the saline portions of the STE as salt water interacts with sedimentary OM reservoirs. Alternatively, CDOM production in saline portions of STEs may result from higher concentrations of labile organic substrates such as the protein like PARAFAC component C5. These results reflect the importance of hydrogeology in regulating delivery rates of TEAs and reactive marine OM to STEs.

Implications for SGD fluxes of OM
The nonconservative behavior of CDOM at our study sites suggests that its biogeochemical production is a common process in STEs regardless of flow rates and residence times. Additionally, produced CDOM quality is similar among STEs, and its spectroscopic characteristics reflect greater aromaticity and lower molecular weight than OM in freshwater or saltwater end members. This OM quality would be consistent with that produced from the microbial degradation of terrestrial solidphase OM in STE sediments or particulate terrestrial OM within karst conduits. The CDOM production would enhance fluxes of CDOM from STEs to surface water, as well as decrease OM storage in STE sediments as OM is transferred from sedimentary to dissolved and mobile forms. The freshwater-salt water mixing zone should thus reduce sedimentary OM storage as it is preferentially remobilized by microbial activity, but the magnitude of mobilization depends on groundwater flow rates, which govern STE residence times.
The production of CDOM in STEs due to freshwater-salt water mixing has important implications for the role of SGD in coastal OM budgets. Many studies estimate SGD fluxes by focusing on "new" solutes delivered to the system in the fresh component of SGD, with little consideration of solutes regenerated from marine sediments. These estimates often depend on the concentration of a solute of interest present in the fresh groundwater end member which is then multiplied by estimates of total fresh SGD that is assessed using seepage meters, chemical tracers, water balance approaches, piezometers, or numerical methods (Burnett et al. 2006;McCoy and Corbett 2009;Knee and Paytan 2011). Our results reflect the importance of freshwater-salt water mixing to production of CDOM in STEs and suggest that concentrations, and thus fluxes, of other associated solutes such as nutrients and metals may be impacted by the introduction of labile marine OM and TEAs.
While CDOM production occurs regardless of flow rate, greater CDOM production is observed in Indian River Lagoon due to relatively longer residence times. The magnitudes of CDOM fluxes from siliciclastic compared to carbonate karst STE types depend on both the concentrations of CDOM in SGD as well as SGD volumes. The combined fresh and saline SGD flux (commonly referred to as total SGD) at Indian River Lagoon SGD is estimated at approximately 43 × 10 6 m 3 yr −1 km −1 of shoreline (Martin et al. 2007) while at the Yucatan is 112 × 10 6 m 3 yr −1 km −1 of shoreline (Table 2; Null et al. 2014). The similarities in water fluxes suggest that solute fluxes will depend less on flow than concentrations, which vary by several orders of magnitude between sites (Fig. 6). Because CDOM production is greater at Indian River Lagoon than the Yucatan, OM fluxes due to SGD are likely to be greater there as a result of longer residence times and enhanced biogeochemical alteration of OM.
This study highlights that biogeochemical processing in STEs leads to CDOM production, but that the magnitude of production depends on STE residence times. However, as observed in Indian River Lagoon STEs, even similar hydrogeological conditions among STEs may lead to variable CDOM production because of OM quality in the inflowing groundwater or sediment composition that may alter TEA distributions within STEs. Despite differences in nonconservative mixing, CDOM quality exhibits similar shifts with salinity between STE types and suggest that biogeochemical processes impacting CDOM quality are similar between STEs regardless of residence time. These systematic differences between STEs with disparate residence times indicate that heterogeneity in hydrogeology and groundwater composition in STEs complicate upscaling of local results in estimates of global impacts of SGD on coastal OM cycling, as well as other solutes linked to biogeochemical reactions. However, ubiquitous production of CDOM at STE sites implies that OM fluxes via SGD may play an important role in carbon budgets in individual study settings.