Salinization of stream water and groundwater at daily to decadal scales in a temperate climate

Elevated salt concentrations in streams draining developed watersheds are well documented, but the effects of hydrologic variability and the role of groundwater in surface water salinization are poorly understood. To characterize these effects, we use long‐term data (12–19 yr) and high‐frequency specific conductance (SPC) data collected from 13 streams across New Hampshire, USA. Concentration–discharge (C–Q) relationships for chloride (Cl−) derived from high‐frequency SPC showed distinct seasonal variability. Diluting behavior was common, but flushing behavior occurred in autumn and winter, suggesting that both groundwater and surface runoff contribute salts to streams. Long‐term data show that although extreme flood events initially reduced salt concentrations in groundwater and rural streams, concentrations recovered to preflood conditions in about a decade. Chronic Cl− exceedances occurred in urban streams during all seasons. This research suggests that variation in stream flow, extreme events and application of deicing agents play a role in freshwater salinization.

Stream water chemistry provides an integrated signal of watershed processes, reflecting the balance between inputs and retention or transformation that drives solute export. Quantifying rates of retention or transformation can be difficult for solutes such as nitrate that are controlled largely by biotic processes (Duncan et al. 2015) but should be easier for conservative solutes such as chloride (Cl À ), which shows clear and rapid responses to watershed inputs (Kaushal et al. 2005). Cl À has direct implications for the health and integrity of freshwater ecosystems ) and the US Environmental Protection Agency (USEPA 1988) has established both chronic (230 mg L À1 ) and acute (860 mg L À1 ) Cl À toxicity criteria to protect aquatic life. In some highly populated and cold climate regions, the main input of Cl À into watersheds is road salt (NaCl) applied as a deicing agent (Trowbridge et al. 2010;Moore et al. 2020).
Numerous studies have documented elevated Cl À in streams draining developed watersheds (Daley et al. 2009;Trowbridge et al. 2010;Moore et al. 2020;Balerna et al. 2021) as well as long-term increases in Cl À (Daley et al. 2009;Kaushal et al. 2018;Moore et al. 2020) due to road salt application. These trends have resulted in the widely recognized freshwater salinization syndrome Kaushal et al. 2018). Road salt can be rapidly flushed to streams during winter storms leading to increased Cl À concentration, especially in urban settings where stormwater infrastructure connects roads and parking lots to stream networks (Niedrist et al. 2021). Nevertheless, Cl À concentration is highest during baseflow periods and diluted with increased flow in many rural and urban watersheds (Daley et al. 2009;Trowbridge et al. 2010;Coble et al. 2018). These differing responses of Cl À to hydrologic variability raise questions about watershed flow paths and the contributions of pavement and surface runoff vs. groundwater sources to the salinization of freshwater ecosystems.
The use of high-frequency sensors is providing unique insights into salinity dynamics in streams (Koenig et al. 2017;Moore et al. 2020;Galella et al. 2021). To our knowledge, no study has coupled sensor technology with long-term sampling of both stream water and groundwater to better understand watershedscale controls and patterns of salinity at daily to decadal scales. Here we use data from a network of high-frequency specific conductance (SPC) sensors as a proxy for salt concentrations along with long-term records (12-19 yr) of stream water and groundwater chemistry to investigate (1) variability in Cl À concentrationdischarge (C-Q) relationships at inter-(seasonal) and intra-annual timescales using SPC as a proxy for Cl À , (2) temporal changes in surface water and groundwater Cl À concentrations in response to extreme flood events, and (3) the frequency of Cl À exceedances above EPA aquatic life criteria in rural and urban streams.

Site descriptions and data collection
We collected weekly to monthly water samples at a total of 13 streams and rivers in New Hampshire, of which 10 sites were instrumented with sensors. Six of the sites (five also included sensors) were located in the Lamprey River Hydrologic Observatory (LRHO; Wymore et al. 2021a; Fig. 1). Altogether these sites drain a gradient in watershed development (Table 1) and were categorized as rural (< 25% developed) or urban (≥ 25% developed). Eight of the sensor sites were loworder streams and two were mainstem rivers, Lamprey River (LMP; USGS 01073500) and Merrimack River (GOF; USGS 01092000). Site LMP drains the LRHO which received 100-yr scale floods in 2006 and 2007 (Coble et al. 2018) and is a rural watershed experiencing increases in development (3.4% of developed watershed area in 2006 and5.8% in 2016;NOAA 2006NOAA , 2016. Sensor sites were instrumented during 2013-2014 and data used in this study (previously published by Potter et al. 2020) were collected year-round (except in stream site BDC which freezes solid during winter) at 15-min intervals through 2017. An EXO2 multiparameter sonde (YSI, Xylem Inc.) was used to measure SPC. Discharge (Q) was measured by the USGS in the two mainstem rivers (LMP and GOF) and estimated in the eight low-order sensor sites using site-specific rating curves and stage height data loggers (Snyder et al. 2018).
Groundwater was collected from two forested riparian well fields in the LRHO (JF, six wells and WHB, five wells). Wells were installed in 2004 and are 5 cm diameter polyvinyl chloride wells with 53-76 cm of perforated well casing. Wells were hand augured to the depth of refusal (approximately 1-2.4 m below the ground surface). Wells were sampled monthly through May 2007, then generally sampled quarterly until May 2021 (Table 1).
Stream and groundwater samples were analyzed for Cl À using ion chromatography with suppressed conductivity detection (Dionex 1000 ICS with autosampler) and data are published by Shattuck et al. (2022). We used the relationship between Cl À and SPC at the 10 sensor sites to estimate Cl À concentrations at 15-min intervals (linear regression; r 2 = 0.74, n = 2213, p < 0.0001; Cl À = 0.1889 Â SPC À 2.709). Other studies have used the relationship between Cl À and SPC to understand temporal trends and evaluate Cl À exceedances (Daley et al. 2009;Trowbridge et al. 2010). We considered exceedance days as any day where estimated or measured Cl À exceeded the EPA chronic (230 mg L À1 ) or acute (860 mg L À1 ) criteria.

Moving window analyses of sensor data
To investigate the annual, seasonal, and intra-annual variability in SPC C-Q relationships, we used a moving window linear regression to generate C-Q slopes from daily mean SPC data. A 30-day window at a timestep of 1 d was used to produce regression coefficients that incorporate antecedent conditions of the previous 30-d period (Fazekas et al. 2020;Wymore et al. 2021b). For each 30-day interval, we obtained the slope b (c = aQ b ), where b is a unit-less exponent representing the slope of the log-transformed C-Q relationship, and a is the intercept in units of concentration. We classified C-Q relationships as chemodynamic when slopes were significantly different than 0. Slopes > 0 were classified as flushing and slopes < 0 as diluting (Moatar et al. 2017;Wymore et al. 2021b). Concentration and Q data were log-transformed to normalize data distribution prior to analysis. To assess the strength of the SPC C-Q relationship, we calculated Pearson's correlation coefficient (r) across each 30-d window. We defined seasons as winter (January-March), spring (April-June), summer (July-September), and autumn (October-December). We used a one-way ANOVA to test for differences in chemodynamic SPC C-Q slopes among seasons for each site (n = 10 for flushing; n = 10 for diluting) followed by Tukey's HSD, significance was determined with a Bonferroni correction.

Temporal changes in Cl À concentrations
We performed three separate analyses to assess temporal changes in surface water and groundwater Cl À concentrations at sites with more than 10 yr of data (n = 5 for surface water; n = 2 for groundwater). First, we used the weighted regression on time, discharge, and seasons (WRTDS) model implemented as the EGRETci package in R (Hirsch and De Cicco 2015) to analyze trends in surface water Cl À concentrations. This method uses weighted linear regression to estimate daily concentration, mean annual flow normalized concentration and requires daily Q throughout the period of concentration data. WRTDS was only suitable for sites LMP (daily Q from USGS station 01073500) and WHB. At WHB we measured daily Q directly from 2012 to 2021, using its site-specific rating curve and 15-min stage height. For earlier years when stage height was not available (2002-2012) we used relationships obtained from 2012 to 2021 between WHB daily Q and the two closest USGS stations (LMP and Oyster River USGS station 01073000; 4.8 and 4.3 km from site WHB, respectively). Our approach used forward multiple linear regression and lagged flows of up to 2 d to account for differences in watershed area (R 2 = 0.74, n = 3054, p < 0.001; Table S1). The second approach used to assess temporal trends was piecewise linear regression (PLR) using the segmented package in R (Muggeo 2017) on surface water sites (n = 5; Table 1) and groundwater sites (n = 2) to identify significant changes in Cl À concentrations through time and any inflection points within the trends. The third approach taken, if no inflection point was detected, was simple linear regression analysis of Cl À concentrations through time.  (2016). Note LMP represents the location for long-term data at LMP73 and sensors located 1 km upstream (LMP72). CB represents site CB03 and PB represents site PB02. WHB represents the stream water site with sensors and the groundwater site.   At sites with Cl À exceedances and more than 10 yr of data (n = 3), we also used simple linear regression analysis to assess temporal changes in annual exceedances. We tested for differences in the number of days of exceedances among streams (n = 4) using a two-way ANOVA with site, season and their interaction. Analyses were conducted in R version 4.0.1 (R Core Team 2021) with significance level set at 0.05.

Results
The C-Q behavior for SPC was typically chemodynamic, displaying a diluting pattern on an intra-annual timescale (82.5% of the time on average), with a range across sites of 70.5 to 92.8% (Figs. 2, S1, S2). However, SPC flushing occurred on average 17.5% of the time (range across sites of 7.2% to 29.5%). Many streams displayed periods of distinct seasonality in their flushing behavior (Fig. 2). All 10 sites exhibited significantly different diluting SPC C-Q slopes among seasons, but which seasons were different varied by site (Table S2). Winter and autumn had diluting behavior that was most commonly significantly different than the other seasons, but whether there was more or less dilution varied by site. Flushing behavior was seldom observed during spring and summer among any of the 10 sites but was often observed in autumn and winter across both years and sites. Winter of 2015 exhibited a sharp increase in C-Q slope at eight of the nine sites, including in one site that normally did not display winter flushing behavior (Fig. 2J) and resulted in an order of magnitude increase in C-Q slope in one stream (Fig. 2G). Several streams displayed a steadily increasing C-Q slope from diluting to flushing starting in spring, through summer and into autumn (e.g., Fig. 2d,J).
Long-term sampling showed Cl À concentrations declined during the beginning of our sampling period with a singular inflection point followed by increasing Cl À concentrations at surface water sites LMP and WHB (Figs. 3, S3, S4; Tables S3, S4) and groundwater sites (JF and WHB; Figs. 3, S4; Table S4). Inflection points occurred 0-9 yr after the last 100-yr flood. Inflection points identified via WRTDS and PLR modeling both occurred in water year (WY; starts 1 October) 2010 at site LMP, but differed at site WHB (WY 2011 for WRTDS; WY 2007 for PLR). The inflection point at groundwater sites occurred in WY 2010 at JF and WY 2016 at WHB. Groundwater and LMP Cl À concentrations exceeded preflood concentrations by 2020. Surface water Cl À concentrations declined faster and exhibited a slower recovery to preflood levels than groundwater at site WHB (Table S4). Although surface water Cl À concentrations at WHB have increased to 2005 levels, they had not reached the highest preflood concentrations by 2021. Cl À concentrations increased over time at urban sites CB03 and MLB (Fig. S5) with no inflection point identified. There was no statistically significant temporal change in Cl À concentrations at the one urban site with Cl À data prior to the 100-yr floods (PB02; Fig. S5) or in annual Cl À exceedances at any site.
Exceedances of the EPA chronic Cl À criterion were observed in the four urban sites where mean Cl À concentrations ranged from 125 to 281 mg L À1 . In the urban stream where sensors were used to estimate Cl À from SPC (McQ), an average of 14 d per year exceeded the chronic Cl À criterion, but exceedances did not occur during summer (Figs. 4, S6). Among the three urban sites sampled monthly, average chronic Cl À exceedances ranged from 3 to 7 d per year with exceedances observed in all seasons (Fig. 4). There was no significant difference between winter or summer exceedances except at site McQ. Winter and summer exceedances were significantly less than annual exceedances at two sites (CB03 and MLB for winter; CB03 and McQ for summer; Table S5).
The acute Cl À criterion was exceeded at the urban sensor site (McQ) 2 d per year on average (range 0-5 d per year; Fig. S6), and only once in 2 of the urban sites (Fig. S5a,b) where monthly samples were collected. Had we relied solely on sample collection to quantify exceedances at the urban sensor site (McQ), we would have documented only one chronic exceedance among the 59 samples collected over the almost 3-yr period, and no acute exceedances. No exceedances were documented at any rural site, where mean Cl À concentrations reached up to 55 mg L À1 .

Discussion
The moving window approach provides a more complete assessment of temporal structure in C-Q behavior than is possible to obtain through the analysis of long-term trends in samples collected (e.g., Daley et al. 2009) or analysis of a series of individual storms with high-frequency sensors (e.g., Koenig et al. 2017). An analysis of more than 700 individual storm events from this same data set also documented a general dilution response of SPC, as well as anticlockwise hysteresis (Koenig et al. 2017). Other regional studies also report dilution of Cl À at high flows (Trowbridge et al. 2010;Coble et al. 2018;Balerna et al. 2021) and periodic chemostatic (Balerna et al. 2021) C-Q behavior. Despite this prevalence of dilution at high flows, our analysis documents that flushing behavior also occurs with a distinct seasonality that tends to be overlooked in more conventional analyses. These seasonal patterns are generalizable across the entire LRHO watershed, and their quantification and characterization may be critical for the development of regional watershed management protocols aimed at protecting aquatic life.
Multidecade time series are needed to understand how linkages between surface waters and groundwaters are altered by extreme flooding events. The coupled response of Cl À in surface water at site LMP and groundwater at site JF to two 100-yr scale flood events in the LRHO suggests that groundwater is an important driver of surface water Cl À in this watershed. The decoupling of temporal trends in surface water and groundwater Cl À concentrations at site WHB suggests that surface water concentrations are less tightly coupled to groundwater in this small tributary watershed, and more responsive to local watershed conditions. We infer from the decadal-scale responses in both surface and groundwaters that extreme events can alter groundwater and surface water salinity for many years, as suggested previously by studies of soil and groundwater residence times (Gutchess et al. 2016;Robinson and Hasenmueller 2017;Kelly et al. 2019). Management actions to reduce surface water salinity thus need to take groundwater residence time for Cl À into account.
The combination of sensor data to capture seasonal patterns in peak Cl À concentrations and decadal sampling to document the long-term effects of extreme events provides important insights into management strategies for documenting and mitigating the impacts of freshwater salinization. Flushing of salt in the short-term, such as we have demonstrated with the moving window approach to C-Q behavior, provides insights into when periods of high Cl À and possible exceedance of aquatic life criteria are likely to occur. Timing and total amount of road salt application are likely to affect the periods when highest Cl À concentrations are likely to occur, and both might be altered to reduce exceedances of aquatic life criteria. Reduction in groundwater Cl À concentrations can also occur during extreme flood events that deliver relatively low-Cl À rainwater to groundwaters in sufficient volume to dilute Cl À levels. These long-term records provide useful evidence of groundwater transit time that can be valuable in devising realistic management strategies that reduce overall Cl À in impaired waters.
Our results are consistent with earlier results that established watershed development thresholds for Cl À impairment, and underscore the importance of high-frequency sensors for determining the extent of Cl À exceedances. Throughout Mid-Atlantic and New England watersheds Cl À concentrations exceeded the chronic criterion beginning at 10% impervious cover/25% urban cover (Corsi et al. 2015) or at median Cl À concentrations of 30-80 mg L À1 (Moore et al. 2020). Another New Hampshire study found Cl À exceedances when mean concentrations were greater than 102 mg L À1 (Trowbridge et al. 2010). Without the use of high-frequency sensors, we would have missed 98% of the chronic exceedance days at the most developed site (McQ) and 100% of the acute exceedance days. Although a monthly sampling regime can approximate seasonal patterns of minimum and maximum SPC (Timpano et al. 2018), our results suggests that Cl À exceedances can be severely underestimated when relying solely on lower-frequency manual sample collection or when samples are collected only during one season. The most developed site (89% developed watershed area; McQ) did not experience any summertime exceedances, suggesting that road salt is flushed through surficial flow paths during the winter-spring transition and does not accumulate in the groundwater. Other studies of urban watersheds have also documented peak Cl À concentrations in the winter (Corsi et al. 2015;Bird et al. 2018;Niedrist et al. 2021), likely due to the high degree of impervious cover and urban infrastructure that quickly routes surface runoff to streams through stormwater piping (Niedrist et al. 2021). Roadside ditches and swales are common among our other sites and road salt may accumulate in these areas (Lazarcik and Dibb 2017). Road salt in these near-road areas may be either quickly washed to streams during storms or infiltrate into the groundwater system that provides summertime baseflow, leading to exceedances during winter when road salt is applied and also in summer, when salts are concentrated in baseflow.
Although another study of 19 streams across the northern US has documented increased frequency of Cl À exceedances from 1990 to 2011 (Corsi et al. 2015), we have not yet observed a similar change. We have, however, observed an increase in Cl À concentrations at two of the urban sites with Cl À exceedances (Fig. S5) and using SPC as a proxy for Cl À , salinity increased approximately three times at LMP since records began in 1953 ( Fig. S7 and historical data). Corsi et al. (2015) also show that Cl À concentrations have approximately doubled over two decades, outpacing the rate of urbanization. Watersheds in the mid-Atlantic region of the United States have exhibited increasing salinity despite steady-state impervious surface cover (Bird et al. 2018). Both studies attribute these Cl À increases to the fact that deicing salt sales outpaced urban landcover growth by 40% from 1987 to 2021 (Corsi et al. 2015). While small changes in development over our study period have occurred, we believe that the increased Cl À concentrations at our sites are likely due to an increase in road salt application rates.
Continued population growth, land use change and climate change could affect road salt application rates in cold climate regions and the salinization of freshwater systems through several confounding mechanisms. Many surface waters are predicted to become toxic to freshwater life within the next century if current salinization rates continue due to changes in impervious surface coverage (Kaushal et al. 2005). The Northeastern US is projected to experience increases in total precipitation during winter and spring with a greater percentage of precipitation falling as rain (Dupigny-Giroux et al. 2018). Although less snowfall would result in lower road salt application rates, an increased frequency of freezing rain may require more use of deicing materials (Minnesota 2022).
With earlier snowmelt-related runoff (Dupigny-Giroux et al. 2018) and more frequent winter high flow pulses associated with rainstorms rather than snowstorms, we expect increased winter flushing of Cl À . Warmer summers with lower minimum flows are also predicted for the Northeast (Dupigny-Giroux et al. 2018) and these conditions may concentrate Cl À in groundwater and baseflow. The latter climate scenario may result in more streams reaching the chronic Cl À toxicity level unless more frequent extreme precipitation events are able to dilute groundwater reservoirs of Cl À , resetting surface water to lower salinity levels at base flow. Coupling these climate change scenarios with future increases in impervious surfaces, it is likely that negative impacts to aquatic life will continue until effective and ecologically sound road salt alternatives are implemented or until road salt use is optimized.