Trends and legacy of freshwater salinization: untangling over 50 years of stream chloride monitoring

Excessive use of road salts to maintain safe winter travel conditions leads to increasing chloride (Cl) concentrations in streams, damaging the structure and function of freshwater ecosystems. Long-term increasing stream Cl trends are generally attributed to increases in urban land cover, however recent research shows that even relatively rural streams can retain Cl and exceed water quality guidelines in summer after road salting has stopped. Untangling the relative influences of long-term changes in streamflow and urban growth on Cl trends is critical for making informed decisions about road salt management. The portion of Cl trends not explained by changes in streamflow or urban growth could be due to changes in road salt application rates and/or legacy Cl in groundwater that is slowly making its way to streams. This study assessed seasonal, long-term stream Cl trends across the Province of Ontario, Canada, where urbanization accelerated and road salt management plans started to develop since early 2000s. We compared stream Cl trends over salting and non-salting seasons with urban growth estimates from two independent time periods, 1965–1995 and 2002–2018. For a subset of sites with sufficient flow data in the periods analyzed, we parsed the seasonal trends into flow and management trend components. We found that most of the variance in the management trend component in the winter salting season could be explained by urbanization, while about half of it could be explained in the summer non-salting season. We further analyzed Cl estimates in low-flow conditions to explore the extent of subsurface contributions to Cl trends, and concluded with a summary of challenges and recommendations for future studies on road salt legacy in streams.


Introduction
Chloride ('Cl') in surface water has long been considered the main indicator of freshwater salinization from road salts-a rising concern in seasonally frozen environments that require road salt application to maintain safe winter travel conditions (Mayer et al 1999, Oswald et al 2019. Elevated Cl levels in streams have damaging effects on the structure and function of many freshwater species, including the disruption of cell homeostasis, reproductive and developmental impairments, and death (Evans and Frick for more urban spaces. This poses a problem for urbanizing regions where stream Cl levels already regularly exceed federal water quality guidelines of 120 mg l −1 for chronic and 640 mg l −1 for acute toxicity in aquatic life, increasing even in summer when road salting has stopped (Canadian Council of Ministers of the Environment (CCME) 2011, Todd andKaltenecker 2012, Lake Simcoe Region Conservation Authority 2015). Increasing stream Cl trends in winter and summer seasons have led to concerns about the build-up of legacy Cl in soil and groundwater (Kelly et al 2008, Eyles et al 2013, Perera et al 2013. While many studies have established urbanization as the key driver of Cl increases in North American freshwaters since the 1980s (Chapra et al 2009, Corsi et al 2015, Kaushal et al 2018, Dugan et al 2020, few have examined the combined influence of urban growth, changes in streamflow and watershed management over time on stream Cl concentrations-especially at a local scale of small to medium ('mesoscale') catchments where these influences are more defined than in larger watersheds (Claessens et al 2006). Untangling these potential drivers of increasing stream Cl concentrations will help inform decisions to mitigate Cl pollution, along with an improved understanding of the portion of trends that can be attributed to legacy.
Previous studies show that watersheds can retain up to 90% of Cl applied as road salt (Kelly et al 2008, Novotny et al 2009, Gutchess et al 2018, Oswald et al 2019. While long-term mass balance approaches have proved essential for exploring legacy effects of other solutes such as nitrate (van Meter And Basu 2017), reliable data for Cl inputs to both public and private impervious surfaces are rarely available over large spatial and temporal scales (Dugan et al 2017, Oswald et al 2019. However, there is significant investment in long-term monitoring of Cl concentrations in streams draining mesoscale watersheds in many regions. While long-term data on Cl outputs will not provide a complete picture of Cl retention in watersheds, they may offer clues to the potential contribution of legacy Cl to long-term trends. This study focuses on the following questions: (a) what are the seasonal, long-term trends in stream Cl concentrations in mesoscale catchments across an urban-rural gradient? (b) What are the relative influences of changes in streamflow, watershed management, and urban growth on stream Cl trends? (c) What can we learn about Cl legacy from longterm summer low-flow (LF) trends? To address these questions, we utilized an extensive open dataset to quantify Cl concentration trends in mesoscale catchments across the Province of Ontario, Canada for two time periods (1965-1995 and 2002-2018). Many watersheds in this study area were undergoing a steady pace of urbanization in the earlier period, while some of the southern watersheds experienced a much higher rate of urban growth in the more recent period (Government of Ontario 2019). Changes in sampling protocols as well as the implementation of salt management programs in several areas from 2001 onwards also sets these two periods apart. This study extends several other works on seasonal water quality trends in Ontario (Todd andKaltenecker 2012, Stammler et al 2017), effects of streamflow variability (Hirsch et al 2010, Choquette et al 2019, Murphy and Sprague 2019) and legacy (Oswald et al 2019, Johnson andStets 2020) to parse the influences of urbanization, streamflow variability and watershed management as drivers of stream Cl trends. We begin by comparing observed Cl trends over winter (salting) and summer (non-salting) seasons with urban growth estimates in a large number of watersheds, then adjust the trends for non-stationarity in streamflow to get a flow and management trend component (MTC) in a small subset of sites. In the absence of reliable Cl input data, the portion of winter trends that are not explained by changes in streamflow or urban growth, can serve as an indicator of changes in winter maintenance practices and legacy contribution. As previous studies on Cl retention and summer trends have found that urbanization and road density outweigh other factors such as watershed-scale differences in structural, geological and agricultural inputs in these regions (Todd andKaltenecker 2012, Oswald et al 2019), the unexplained variance in summer trends includes the effects of legacy Cl from subsurface storage. We further examine this legacy contribution by estimating stream Cl trends during summer LF conditions, when changes in winter management is expected to have little to no effect.

Study area and data
Ontario is the second largest and most populated province in Canada, including parts of the watersheds for four of the five Great Lakes and the St. Lawrence River basin (Government of Ontario n.d.). With a latitude range of 42 • N-56 • N, Ontario includes a wide range of climate, land use, and geology. Climate studies in this region have shown significant variability in seasonal and annual streamflow patterns since the 1920s (Nalley et al 2016(Nalley et al , 2019. With rapid urbanization in Southern Ontario and addition of road salt in the Priority Substances List (Government of Canada 2001, Government of Ontario 2005), major urban centers now require salt management plans under the federal code of practice for environmental management of road salts (Environment Canada, E.P.S. 2004, Environment andClimate Change Canada 2020).
Long-term monitoring of Ontario's surface water quality started in 1964, with the provincial water quality monitoring network (PWQMN), though sampling locations and frequencies varied over time and space. More data are available in southern and southwestern Ontario, where there is a mix of  urban, rural and natural land uses (figure 1). No changes in analytical or sampling methodologies were recorded for Cl (Stammler et al 2017). However, an exploratory data analysis of the stream Cl data revealed disparities in sampling frequencies between cold and warm seasons. While more data are available for summer months, winter sampling declined after 1995, with a data gap until 2002 (figure 2). To avoid sampling bias and to observe any changes in seasonal trends after 2000, we divided the analysis into two independent time intervals: 1965-1995 and 2002-2018. To calculate seasonal trends, we divided monthly stream Cl data into winter salting season (November-April) and summer non-salting season (May-October). This seasonal division is similar to that observed in other water quality studies conducted in Ontario (Stammler et al 2017, Oswald et al 2019. However, unlike these studies which aggregated data into seasonal averages for determining trends, we preserved the monthly resolution of the dataset to account for inter-seasonal variability.

Chloride trend analysis 2.2.1. Seasonal Mann-Kendall trend analysis
Chloride data were filtered to contain no sites with a gap over one-third of the total record (Helsel and Hirsch 2002). Since summer samples were relatively abundant, we filtered all sites to have a minimum four months of data available out of the 6 month non-salting season (May-October), although this was not possible for winter where data were limited. After testing for normality using Shapiro-Wilk tests, a Seasonal Mann-Kendall test for nonparametric data was conducted on the monthly data for each season and time period separately. Trends were represented by Sen Slopes as rates of change in stream Cl (mg l −1 yr −1 ), which were then compared with urban growth in mesoscale catchments with drainage area of 10-1100 km 2 , following a similar criteria for watershed size used by Stammler et al (2017). Watersheds were delineated using the Ontario flow assessment tool (OFAT) (Ontario Ministry of Natural Resources and Forestry 2013). Nested catchments were retained in this analysis, since urbanization tends to occur in clusters around existing urban areas, affecting stream outlets from smaller urbanized catchments more heavily than the larger watershed area. While this may introduce some interdependence in the data, it is not expected to affect the results as all calculations were normalized by the upstream catchment area of the site (Lintern et al 2018). All statistical analyses were conducted using R version 4.0 (R Core Team 2017).

Urban growth analysis
To ensure consistency in land use categories (Penman et al 2003), urban areas were extracted from Agriculture and Agri-Food Canada land use maps, which were available for the years 1990, 2000, and 2010 (Government of Canada n.d.). The earliest land use for our time periods was available from a historical map of built areas by Canada Land Inventory (CLI), originally created from 1950 to 1966 data to estimate land capability for agriculture (Ministry of Natural Resources n.d., Department of Regional Economic Expansion 1970). Urban growth in the 1965-1995 period was estimated by subtracting urban areas (as % of watershed) in 1966 from that in 1990. The same estimation was carried out for 2000 and 2010 urban areas for the 2002-2018 period. This provided a sound, albeit conservative, approximation of the overall urban growth in these catchments (The Neptis Foundation n.d.). It should be noted that our urban growth estimates for 1965-1995 period do not account for changes in road area. As the historical CLI maps comprised areas of large urban 'settlements' (including open space and parks) with no separate data for roads, we used the similar settlements category from 1990 for comparison. Although the latest available land use maps were from 2010, the 2002-2018 urban growth estimates are presumably more precise as we were able to account for road areas as well as settlements. All geographical analyses were carried out in ArcGIS Pro 2.6.

Weighted regressions on time, discharge and season (WRTDS) analysis
For a subset of sites with co-located long-term daily flow data from the water survey of Canada (WSC), we normalized stream Cl concentrations for streamflow using the WRTDS method (Hirsch et al 2010). This model assumes that both the concentration-discharge relationship and discharge over time can vary for each site, and uses multiple sets of fitted coefficients to best describe water quality with different flow conditions in different seasons and time periods (Corsi et al 2015, Johnson andStets 2020). This is presented by the following equation: where C is concentration, β 0−4 are fitted coefficients, t is time, Q is discharge, and ε is the residual variability.
The WRTDS model was recently updated to control for both abrupt and long-term changes in flow over time (Choquette et al 2019). Murphy and Sprague (2019) used it to parse flow-normalized (FN) water quality trends from the United States Geological Survey (USGS) dataset into a flow trend component (QTC) and MTC, where MTC is calculated as the water quality trend assuming streamflow remained stationary over the period of analysis, and QTC is the difference between MTC and the trend allowing for non-stationarity in flow over time (i.e. FN trend − MTC = QTC). The QTC therefore represents influences of streamflow variability, while MTC represents all other influences from watershed management, including urban growth. While Murphy and Sprague defined trends as the difference in concentration estimates between start and end years of their time period, we represented trends as the rates of increase/decrease in stream Cl. We used the R packages 'tidyhydat' to extract WSC flow data and 'EGRET' (Exploration and Graphics of RivEr Trends) to run the model with the WSC and PWQMN data (Albers 2017, Hirsch et al 2018. Streamflow sites for the model were carefully selected to meet the USGS criteria recommended for WRTDS by Oelsner et al (2017). These flow gauges were matched with PWQMN sites by proximity and branching: in the same stream order within 2 km, with little to no tributaries in the middle (i.e. at least three orders smaller and producing <10% difference in drainage area). Nested catchments were removed from the subsets for this analysis. Model periods and seasons were set to be the same as the Mann-Kendall analysis for comparison.

Stream chloride trends in low-flow conditions
A Cl mass balance study of multiple watersheds in Southern Ontario found that watersheds with lower and more disconnected urban coverage had higher Cl retention, regardless of structural differences in watershed characteristics (Oswald et al 2019). This indicates that elevated summer Cl concentrations in less urbanized watersheds may, in part, be driven by the retention and eventual release of legacy Cl. To assess this legacy contribution, we calculated Cl trends during Low-Flow (LF) conditions in the driest months of summer (July-September), when groundwater is the dominant source of streamflow (Coulibaly and Burn 2005). This method was recently used by Johnson and Stets (2020), who defined LF as the 30th percentile of flow and extracted the corresponding concentrations from WRTDS as indicators of legacy nitrates in streams. While the authors used winter concentrations to observe nitrate legacy in times of lowered biotic activity, we chose the late summer months for Cl when chances of capturing any delayed spring-melt effects from the preceding winter salting season is low. We calculated LF as the 30th percentile of daily flow during July-September of each watershed for the two time periods, and extracted the Cl concentrations falling under these LF values assuming stationary flow (low-flow management trend component or LFMTC), to examine the management component in summer trends that could be from legacy. Since LF conditions do not occur every year, resulting in some missing data, we used the same approach as the Mann-Kendall analysis to estimate trends in LFMTC.

Seasonal Mann-Kendall analysis and urban growth
There were 82 mesoscale watersheds in the 1965-1995 period and 196 watersheds in 2002-2018 with sufficient long-term data for both seasons. Over 83% of sites in the 1965-1995 period resulted in statistically significant (p < 0.05) seasonal stream Cl trends (table 1; figure 3). In 2002-2018, about 42% of sites showed significant trends in the summer non-salting season, while only 15% of sites showed significant trends in the winter salting season. This was expected given the lack of regular monthly winter samples in this period, which limited the rank order correlation in the Mann-Kendall analysis. Even with the reduced portion of statistical significance in 2002-2018, significant trends from both periods for the winter salting season showed strong linear relationships with urbanization: about half of the variance in trends could be explained by variance in urban growth (table 1). In contrast, while there was little to no correlation between summer trends and urban growth in 1965-1995, the more recent trends from 2002 to 2018 show urban growth can account for about 45% of the variance in stream Cl trends in the summer nonsalting season. The linear Pearson correlations of significant trends with urban growth in all seasons and time periods were very similar to non-linear Spearman rank correlations, the rest of this paper therefore uses these linear relationships as it allows us to compare the shared variance with urban growth (de Winter et al 2016).
Overall, the Mann-Kendall results indicate that winter stream Cl trends are following the pace of urban growth over the time periods (figure 3). While associations of urban growth with summer Cl trends in 1965-1995 may be more difficult to interpret, rapid urbanization in the 2002-2018 period is portraying a pronounced impact on the recent summer trends (table 1). This agrees with previous studies that found that urban land cover is a primary driver of summer increases in stream Cl in Ontario (Todd and Kaltenecker 2012). It should be noted here that the lack of correlation of summer trends with urban growth in 1965-1995 could be due to the inability to include road density in our urban growth estimates for that period. Also, there may be a lag time between increases in urban land cover that is salted and the breakthrough of associated Cl in streams. Increasing winter Cl trends were ubiquitous across the study region and periods despite having less data coverage than summer.

Influence of streamflow on chloride trends
We selected ten sites from 1965 to 1995 and 15 from 2002 to 2018, that met the rigorous WRTDS model assumptions for flow and Cl data for both seasons over these periods (table 2). These sites span a range of drainage areas and urban growth representing the whole dataset. The ten sites from 1965 to 1995 had statistically significant Mann-Kendall trends (increasing or decreasing) in both seasons, while all 15 sites from 2002 to 2018 had statistically significant summer trends but only four were also significant in winter. Despite the poor sampling frequency in winter, these sites had sufficient long-term data available for WRTDS to model Cl estimates from daily flow data. Given the length of the time periods and lack of data for both seasons simultaneously, the same sites could not be used for both 1965-1995and 2002. The linear relationships of Mann-Kendall trends with urban growth in this subset of sites showed strong similarity with those in the larger sample (table 1), indicating this subset to be a good representation of the overall trends. Similar to the broader dataset, the 1965-1995 period had a higher proportion of sites with increasing trends (100% in winter and 90% in     (  in some sites (e.g. Canagagigue Creek, East Humber River, Jackson Creek and Rouge River in figure 5(a)). Not surprisingly, most of these changes in trend direction were observed in the 1970s or 1980s, when improvements in Cl trends were attributed to industrial controls and increases were attributed to rapid urbanization (Chapra et al 2009, Dugan et al 2017. In Jackson Creek, where both time periods were modeled, winter trends in the earlier period showed multiple step trends of increase and decrease in Cl. Further investigation revealed local concerns over flooding and year-to-year flow variation from sporadic weather patterns in this site (City of Peterborough 2010). These results indicate that a linear measure of change, such as concentration units per year, is not the best indicator of real change in such watersheds, although we resort to it for the scope of this study.

Influence of urban growth and road salt management on chloride trends
The WRTDS results show that the linear relationship of winter stream Cl trends with urban growth remained intact regardless of flow-normalization ( figure 4). In contrast, the difference in summer correlations of FN trends with urban growth in the two periods was not as pronounced as observed in the overall Mann-Kendall results (table 1). However, the WRTDS results were less scedastic in 2002-2018 compared to 1965-1995, with the exception of one potential outlier (figures 4(b) and (d)). Removing the suspect data point (Dodd Creek) from the group did not bring any change to the correlation, though this suggests FN trends may show stronger linear relationships with urban growth if we can achieve a larger sample size. It should also be noted that QTC was the highest in watersheds with very little urban growth (table 3). Further comparison using the MTC, i.e. the portion of trend not affected by changes in streamflow, showed strong linear correlations with urban growth in both seasons and time periods. These results show that the MTC is generally the major component of stream Cl trends and urban growth can account for most of its variance in winter salting seasons (r 2 of 80%-88%) and about half of it in summer non-salting seasons (r 2 of 45%-56%). The unexplained variance in the MTC-urban growth relationship could be due to changes in winter maintenance practices and/or legacy Cl. In summer, winter management practices do not directly influence stream Cl concentrations, hence we expect legacy Cl dominates the unexplained variance in this season. We recognize that changes in winter maintenance practices will ultimately impact the build-up of legacy Cl, but we make the assumption here that this impact is delayed enough for the two drivers to be conceptually separated. In both periods, agricultural areas were negatively correlated with summer MTC (although only significant in 2002-2018), confirming that these summer trends are not from agricultural sources.

Legacy contribution in stream Cl trends
While summer FN trends correlate fairly well with winter trends of the same sites (r 2 of 0.47 and 0.76 for 1965-1995 and 2002-2018, respectively), much of the variance remains unexplained. It is not clear how much of the summer trends are driven by the preceding salting season (e.g. delayed spring melt, shallow subsurface transport of Cl to stream). To further explore the potential for legacy Cl to be driving summer trends, we estimated LFMTC rates in the dry and late summer months of July-September to avoid any impacts from winter salt inputs of the same year. Both the FN trend and MTC shared the same values during LF conditions, confirming no influence from streamflow variability (QTC). Therefore, the observed LFMTC rates or changes in late-summer stream Cl suggest contributions from the subsurface. This may include legacy Cl in soil and groundwater from past road salting or other anthropogenic and natural sources (e.g. bedrock geology).
Our results showed that much of the variance in summer Cl MTC could be explained by the LFMTC rates (r 2 of 0.68 and 0.89 for 1965-1995 and 2002-2018, respectively). Higher rates of increasing LFMTC were not always associated with rapidly urbanizing watersheds (table 4)    Studies that focus on assessing the variability in road salt application rates to private areas (including parking lots), and examine spatial and temporal distribution of other Cl sources.
Long-term data on Cl inputs and outputs A greater challenge lies in acquiring the above input datasets for longer than a few years, whereas Cl retention and legacy can occur over several decades (Kelly et al 2019, Oswald et al 2019.
Likewise, consistent long-term water quality monitoring that is co-located with flow monitoring is not widely available.
Greater priority on co-located streamflow and water quality stations. Further exploration of the usefulness of longterm Cl, streamflow, and urban growth data for inferring the influence of legacy Cl on freshwater ecosystems.
Subsurface Cl storage Despite evidence of the importance of groundwater contributions to streams, there is a paucity of studies examining the role of groundwater-surface water interactions on stream water quality (Staudinger et al 2019, Condon et al 2020.
Process-based studies on the role of groundwater-surface water interactions in driving stream Cl dynamics. Water quality models that account for both overland and subsurface and estimate the timing and magnitude of stream Cl response to changes in road salting. Urbanizing areas or mixed land use Currently there are several available water quality models designed and widely used for nutrient modeling in agricultural watersheds. However, models that can be implemented for predicting watershed-scale Cl transport and account for urban land use and infrastructure are rare (Wellen et al 2015). Urban water quality models that can simulate urban drainage structures tend to overlook or homogenize natural watershed processes, making it difficult to model non-point sources of contamination and legacy contributions (Bach et al 2014, Fu et al 2020. Improved water quality models in areas with mixed land use that can represent natural hydrological flow paths (e.g. overland flow, soil, and groundwater flow) and artificial drainage (e.g. flow through storm sewers, stormwater control measures). retention is often significantly high in rural and other permeable watersheds where paved surfaces are located further from the outlets. In highly urban watersheds, most of the salt-laden meltwater from winter and spring is flushed to streams via sewers, with fewer opportunities for transport to the subsurface. Regionally, summer is the baseflow season and many of the watersheds with high LFMTC are in areas with high BFIs (Coulibaly and Burn 2005, Neff et al 2005, Sharpe et al 2014. Although specific values of BFI could not be obtained for these study periods (table 2), Moira, Sydenham, Ganaraska, Jock Rivers, and Jackson Creek (all of which have LFMTC ⩾ summer MTC) were previously identified as watersheds with the highest BFI in a study of 115 Ontario watersheds (Rudra et al 2015).
Interestingly, while urban growth showed very little correlation with the 1965-1995 LFMTC trends (r 2 of 0.10), this improved for the 2002-2018 sites (r 2 of 0.62). This aligns with our previous suggestions that the greater uncertainty in historic urban growth estimates and failure to include road density changes in the earlier period could be a limitation to inferring results. It is also possible that the recent rates of urban growth and sudden increase in salts may be contributing to Cl peaks beyond the winter (when interflow and infiltration are slow through frozen ground) and spring months in the continuously urbanizing areas. In these areas there still remain sufficient permeable soils between urban centers not directly connected to each other, to provide additional Cl from both groundwater storage and delayed lateral flow through the shallow subsurface, resulting in higher concentrations in summer than winter (Eyles et al 2013, Corsi et al 2015). Murphy and Sprague (2019) introduced MTC as the 'catch-all' term for all controllable changes in a watershed. Untangling the LF water quality trends to examine effects of any lagged sources can help understand how much of the watershed management can still be out of our immediate control.

Recommendations for long-term monitoring of Cl to support legacy studies
Despite progress in using long-term stream water quality data to infer legacy effects of nutrients in watersheds (Johnson and Stets 2020), challenges remain with inferring legacy Cl contributions to increasing stream Cl trends. Although this study provides an improved understanding of legacy contribution in some watersheds, quantifying the timing and magnitude of Cl retention and release is impossible without sufficient data on salt inputs. In the absence of high quality Cl input data, trends in Cl outputs should be tracked across urban-rural gradients to identify significant changes in Cl concentrations, especially as they relate to efforts to mitigate Cl pollution. However, given the rigorous data requirements for models such as WRTDS, we still find the current availability of co-located discharge and water quality data limits the number of sites that can be used in such studies. In table 5, we provide a list of challenges that prevent us from fully understanding and responding to short-and long-term Cl pollution. We also provide a set of recommendations, aimed at multiple stakeholders (e.g. water quality scientists, watershed managers, municipalities), to overcome these challenges.

Conclusion
Studies exploring the contribution of legacy Cl to long-term trends in stream Cl are rare and usually focused on individual watersheds (Kelly et al 2008, Novotny et al 2009, Perera et al 2013. A few studies have explored warm season Cl trends and annual Cl retention rates in Ontario streams (Todd andKaltenecker 2012, Oswald et al 2019), however they do not address the contribution of legacy Cl.
This study leverages open water quality and streamflow data to untangle potential drivers of stream Cl trends in Ontario, Canada since the 1960s: changes in streamflow and watershed management, which includes controllable changes in the watershed such as urban growth and road salt application rates. By estimating the portion of trends explained by these factors, we infer the potential contribution of legacy Cl. Our findings from a subset of sites indicate that, while long-term changes in streamflow can have a wide range of influence on stream Cl concentrations, they cannot explain seasonal trends in the study region. Urban growth is the primary driver of winter trends, and can explain up to half of the summer trends. Further analysis using stream Cl estimates during dry, LF conditions in the late summer months indicate that subsurface contributions to streams may be driving the summer trends in many of these watersheds, especially those that are relatively non-urban, suggesting the possible extent of effects by legacy Cl and any other sources not within the immediate control of watershed management.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: https:// data.ontario.ca/dataset/provincial-stream-waterquality-monitoring-network.