Significant stream chemistry response to temperature variations in a high-elevation mountain watershed

High-elevation mountain regions, central to global freshwater supply, are experiencing more rapid warming than low-elevation locations. High-elevation streams are therefore potentially critical indicators for earth system and water chemistry response to warming. Here we present concerted hydroclimatic and biogeochemical data from Coal Creek, Colorado in the central Rocky Mountains at elevations of 2700 to 3700 m, where air temperatures have increased by about 2 °C since 1980. We analyzed water chemistry every other day from 2016 to 2019. Water chemistry data indicate distinct responses of different solutes to inter-annual hydroclimatic variations. Specifically, the concentrations of solutes from rock weathering are stable inter-annually. Solutes that are active in soils, including dissolved organic carbon, vary dramatically, with double to triple peak concentrations occurring during snowmelt and in warm years. We advocate for consistent and persistent monitoring of high elevation streams to record early glimpse of earth surface response to warming. Dissolved organic carbon in high-mountain streams respond strongly to temperature variability, peaking at snowmelt and increasing by up to three times in warm years, according to hydrological, meteorological and geochemical data from Coal Creek, Colorado.

H igh-elevation mountains (HEM) are "water towers for humanity" 1 . More than one-sixth of the Earth's population relies on glaciers and seasonal snowpacks for water supply 2 . In the western United States, mountain snowpack is estimated to provide nearly three-quarters of the freshwater supply for >60 million people 3 . HEMs are experiencing faster warming than low-elevation places [4][5][6] . Warming shifts precipitation from snow to rain; shrinking snowpack induces early snowmelt and prolonged summer droughts [7][8][9][10][11] .
These hydrological alterations can have profound impacts on ecosystems and water quality 12,13 . Droughts can alter biogeochemical transformation and water chemistry 14,15 , whereas large snow-melt events often flush out disproportionally large pulses of "stored" solutes 16,17 . It however has remained poorly understood how and by how much soil biogeochemical processes and weathering respond to warming to ultimately modify water quality in HEMs 3 . Alpine lakes and watersheds have observed paralleled increase in sulfate and cations [18][19][20] . Crawford et al. 21 suggested that this increase may arise from a lowered groundwater water table that exposes sulfide-containing minerals (e.g., pyrite) to oxic conditions and accelerates sulfate mobilization. In Alaska, Yukon River Basin underlain by discontinuous permafrost has seen elevated concentrations of sulfate and divalent cations but has seen no changes in Dissolved Organic Carbon (DOC) 22 . DOC however has been observed to gradually increase in Europe, North America, and Asia, which has been attributed to changing climate or recovery from acid rain [23][24][25] . These observations suggest that different solutes and biogeochemical reactions respond to warming and environmental drivers distinctively. By and large, however, we do not understand the impacts of warming on stream chemistry and how they will evolve as the pace of warming accelerates.
This work focuses on a HEM in Colorado and asks the following questions: (1) How do different soil biogeochemical processes and therefore solutes respond to temperature and hydrological variation distinctively? (2) How do solute variations in recent years compare to their long-term record? We present a rare, concurrent set of meteorological, hydrological, and geochemical data from Coal Creek, Colorado, a HEM watershed (elevation 2700-3700 m) in the central Rocky Mountains (Supplementary Fig. 1). The air temperature has increased by about 2.0°C since 1980, and the annual mean minimum temperature has just surpassed zero in 2018. We compare 4 years of highfrequency measurements (every other day) from 2016 to 2019 for >40 solutes spanning the entire periodic table of elements that have dissolved forms in water. The recent measurements show that solutes derived from distinct biogeochemical processes and source zones respond differently to inter-annual temperature variations. Solutes active in shallow soils have experienced concentration increases by a factor of more than three in recent years, whereas geogenic solutes derived from rock weathering in deeper subsurface have retained relatively similar inter-annual patterns. Juxtaposing recent frequent measurements with long-term bimonthly (once every 2 months) data, however, cannot reach conclusive insights because of the inconsistent data frequency. We advocate that HEMs experiencing rapid warming should be established as "sentinel watersheds" that record early glimpses of Earth surface response to warming. They should be monitored persistently and consistently so as to document long-term responses to warming.

Results
Persistent warming from 1981 to 2019. Data from the local Snow Telemetry (SNOTEL) station from 1981 to 2019 indicate the annual average minimum and maximum daily temperature have increased by~1.6 and 2°C, respectively (Fig. 1a). The daily data, corrected for sensor bias 26 , is comparable to an increase of2°C in Colorado and <1°C in the U.S., and <1°C globally (National Oceanic and Atmospheric Administration, National Centers for Environmental information, https://www.ncdc.noaa. gov/cag/). The mean and maximum discharge were highly responsive to Snow Water Equivalent (SWE, Fig. 1c). Declines in SWE since 1981 have occurred (25 mm/decade, or 6% of average) but are not statistically significant (p = 0.20). Snow fraction (snowfall to total annual precipitation) appear to be declining on the order of~1% per year, but the statistical significance is also poor (p = 0.46). Milly and Dunne 11 found that without changing precipitation, Colorado River annual streamflow will drop 9% for 2°C temperature increase. From 2016 to 2019, the annual temperature in 2016 was near the average for the 1981-2019 period, 2017 and 2018 were warmer years, and 2019 was a cold year. In 2019, average daily maximum (T max ) and minimum (T min ) temperatures fell below the warming trend but near the average for the 1981-2019 period. Therefore, the temperature increase during 2015-2018, superimposed upon an already accelerated warming trend, likely exacerbated drought conditions in 2018. Maximum SWE was 34, 54, 26, and 49 cm in 2016, 2017, 2018, and 2019, respectively, compared to an average of 37 cm in the long-term record. Therefore, 2016 was an average water year, whereas 2017 and 2019 were wet and 2018 was dry compared to average. Indeed, 2018 SWE saw the 4th lowest value over the data record, after 1989, 1981, and 2012. For solutes primarily derived from chemical weathering, such as Ca, Mg, SO 4 ( Fig. 2f), Na, Si, and Sr ( Supplementary Fig. 3), their high-frequency concentrations had similar year-to-year patterns. The Ca and Mg concentrations from the recent highfrequency measurements were consistently lower than the longterm USGS low-frequency measurements. This is because of the difference in sampling locations. The recent high-frequency sampling point is at about 2/3 of the creek length toward the outlet (at Coal-11, Supplementary Fig. 1), whereas the long-term sampling point is at the watershed outlet (USGS water quality, Supplementary Fig. 1). In between these two sampling locations, a water treatment plant often discharges lime-treated water to the stream, often elevating Ca and Mg concentrations at the stream outlet. Analysis of water treatment discharge indicates negligible addition of other solutes (e.g., Fe, Al). The concentration of SO 4 has changed slightly from year to year, although without a consistent pattern.
Export patterns and loads. The concentration-discharge (C-Q) figures show how and how much solute concentrations varied as streamflow varies (Fig. 3). DOC concentrations increased consistently with discharge (flushing patterns, maximum concentrations at the peak of snowmelt). The DOC C-Q patterns varied significantly in recent years but had maintained the flushing pattern. For Fe, however, the concentrations increased significantly under low discharge conditions (between August and November) such that their C-Q patterns transitioned from a pronounced flushing pattern in the long-term record to an almost chemostatic pattern, especially in 2018. In fact, concentrations of DOC, Fe, and Al have all elevated substantially in the low discharge drought period in 2018.
The C-Q patterns are encapsulated in the slope b of the powerlaw relationship C = aQ b , where C and Q are stream concentration and discharge, respectively; and a and b are fitting parameters. The slope b is a measure of the sensitivity of solute concentrations to streamflow changes. High absolute b values (negative or positive) indicate high sensitivity to changes in flow regimes, whereas close to zero b values indicate low sensitivity to hydrological changes (chemostatic). The ratios of the coefficient of variation (CV) compare the variation in concentration (CV c ) and in discharge (CV Q ) ( The loads of solute export were largely controlled by the discharge pattern, in particular, the timing and magnitude of pronounced snowmelt peaks each year, as illustrated by the daily loads (Load = C*Q) (Fig. 4). The DOC loads were much higher in the wetter years 2017 and 2019 compared to the dry year 2018, whereas the peaks in 2016 and in long-term records were in between these values. Although the year 2019 had a higher discharge compared to 2017, DOC export in 2019 was nearly half of the 2017 export. For other solutes, however, the loads in 2016 were the highest, elevated by the higher concentrations that year. Flux analysis under different hydrological conditions indicates that during the 6-8 weeks of snowmelt, the watershed exported >70% of annual loads for most solutes. The annual loads depended significantly on the annual discharge for DOC and Al but the correlation was less significant for Fe and Pb.
Bimonthly averages indicate significant changes under postsnowmelt, dry conditions. Additional analysis for concentrations averaged in every 2 months shows that concentration peaks generally occurred during snowmelt, but the percentage increase maximized during post-snowmelt, the low flow period from August to November (Fig. 5). Pb concentrations increased in 2016 and decreased to levels lower than the historical records, especially during snowmelt in 2017 and 2018. Its pre-and postsnowmelt concentrations however were still higher than longterm concentrations.
Annual average concentrations. Annual flow-weighted average concentrations ( Fig. 6) showed that among all these solutes, DOC saw the most significant variations from 3.9 to 8.5 mg/L, with higher concentrations in the warmer years of 2017 and 2018 (Fig. 6a, c). These high concentrations however did not occur at the high Q year of 2019, indicating a stronger control of T compared to discharge. The annual average of DIC varied Notably, the concentrations of DOC, Fe, and Al at low discharge increased significantly in recent years, shifting the export patterns from flushing (concentration increase with discharge) to patterns with no apparent trend but with high variations. In contrast, Pb concentrations remained at relatively similar levels at low flow but decreased at high flow. The gray shaded region in (e) represents the range of b and CV C /CV Q that define chemostasis ( b j j < 0.1 and CV C /CV Q < 0.2) 76 . All four solutes exhibit flushing patterns in the long-term record however transitioned to lower b values toward the b = 0 line, as indicated by the arrow, due to the concentration increase in the summer-fall droughts. Concentration-discharge figures for other solutes are in Supplementary Fig. 4. e Annual load vs. annual discharge. The numbers in the legend are the annual exports ("MegaG" for "mega grams", 10 6 grams, "kg" for kilogram) calculated from the area under load (C*Q) curve. For the LTM data that was sparsely sampled (i.e., monthly or bimonthly), daily flux from different years was put together in the same year with their corresponding dates to maximize LTM flux points for load estimation.  The flow-weighted sulfate concentrations increased from about 8.0 to 9.5 mg/L. It however did not have a clear correlation with temperature nor annual discharge. Ca and Mg mostly comes from rock weathering and exhibited similarly dilution patterns as sulfate. The ratio (Ca + Mg)/DIC is often used as a measure of recovery from acid rain or the influence of other chemicals such as those used in agricultural lands [27][28][29] . The ratio has varied between 3.5 to 5.0 in the past 4 years, with the wet years having lower ratios. The (Ca + Mg)/DIC, however, did not mirror the increasing temporal trend of sulfate, indicating recovery from acid rain may not be the dominant driver.
Extrapolation for long-term trend: Juxtaposing long-term, bimonthly data with short-term, frequent data. The recent high concentrations of some solutes (e.g., DOC, Fe, Al) in Coal Creek were in stark contrast with the decades-long, low-frequency USGS data. To analyze the long-term trend with consistent measurement frequency, we combined a machine learning approach to infer the full record of discharge from 2000 to 2019, and then used LOADEST and USGS bimonthly data to develop concentration time-series models for DOC and Mg as two representative solutes (details in "Methods" section). The predicted concentrations from the model had very similar dynamics as those of long-term bimonthly USGS data (Fig. 7). However, this pattern contrasted the high-frequency data in 2016-2019 that exhibited high DOC concentrations at discharge peaks during snowmelt.

Discussion
Persistent warming in Coal Creek and high elevation mountains. Global climate projection suggests a continued temperature increase in the coming century, in parallel with declining snowpack and freshwater storage 2,3,6,7 . HEMs generally warm at faster rates than other regions 4 . Inter-comparison Project Phase 5 (CMIP5) predicts a global temperature increase of 1.5°C could mean 2.1°C in HEMs and an increase in warming rate by 0.2-0.4°C/km 30 . In the European Alps, the temperature has risen pervasively, mounting to a mean annual temperature increase of about 2°C, twice as much as that of the northern-hemispheric average in the past century 31 . The projected temperature increase in western US mountains is about 0.5°C per decade for the twenty-first century 7 . The temperature increase of 1.6°C (T min ) to 2.1°C (T max ) in Coal Creek over the past four decades is on par with other mountains across the globe. Warming regulates the timing and magnitude of snowmelt and summer droughts. Longer summer droughts have occurred across the US 32 and particularly in the drainage area of Colorado River 11 . Together with other environmental drivers such as acid deposition, changes in temperature and snow/water dynamics can have profound impacts on soil biogeochemical transformation, weathering, and solute export 25,33 . For example, Hirmas et al. 34 showed that the Rocky Mountain regions have the largest increase in soil effective porosity in the US.
Solutes active in soil zones: large variation driven by warming or other drivers? Data from 2016 to 2019 in Coal Creek exhibited marked variations for some solutes but not as much for other solutes. They indicate that significant changes may be happening at least in the shallow soils and that different solutes (and therefore biogeochemical processes) diverge in their responses to warming. Solutes such as DOC, Fe, and Al (and many others in Supplementary Fig. 3, including P, PO 4 , NO 3 ) are among those that have seen large variations. These solutes generally exhibited flushing C-Q patterns, indicating a more abundant presence in shallow soils [35][36][37][38] . Fe-containing minerals, including various forms of iron oxides, commonly exist in soils and play an essential role in biogeochemical cycles 39 . Aluminum is a common element in organic matter that can stimulate or inhibit plant growth 40 . It is often used as an indicator of soil response to acid rain, as soil buffers increased acid loads by releasing Al. Toxic levels of Al have been observed in areas impacted by acid rain [41][42][43][44][45] .
The lack of TN variation may be because N is tightly cycled in forests with a limited release of reactive nitrogen species 46 . The co-occurrence of increased concentration (of Fe and Al, among other solutes) and sampling in 2016 may suggest changes have already taken place before 2016. The bimonthly data, because of their low frequency, may have missed important spikes in the early years and therefore do not allow rigorous comparison with high-frequency data. Without the presence of carbonate minerals at the site, DIC primarily originated from soil respiration and should also reflect changes if soil processes did alter. From Fig. 2, it appears that DIC did not change as much as DOC. However, a closer examination of time-series and C-Q figures revealed the DIC concentrations under high flow conditions in fact escalated in 2017 and 2018 ( Fig. 2 and Supplementary Fig. 4). These higher DIC may indicate higher soil CO 2 in these 2 years during or right after snowmelt. The higher flow-weighted concentrations in 2018 Comparison of USGS bimonthly data (black+) and recent high-frequency data (red square) with the model prediction of DOC and Mg at the high-frequency sampling point (gray lines). The predictions were from concentration time-series models developed using LOADEST, extended discharge record produced from a machine learning approach, and bimonthly USGS data (see the section of "Extrapolation for long-term trend" in Methods for details). The model did not use high-frequency data to avoid bias from inconsistent data frequency (or the problem of temporal aliasing) The predictions from models trained using bimonthly USGS data indicate that DOC and Mg should show very similar patterns as in the long-term record. This differed from the high-frequency data. and 2019 may reflect the overall "delayed" increase of DIC in 2019 responding to higher OM decomposition in 2017 and 2018.
Although concentrations have increased in both snowmelt and post-snowmelt seasons, the warm and dry post-snowmelt seasons have experienced a higher percentage increase compared to pre-snowmelt or snow-melt periods. This highlights the significant role of summer-fall drought on producing DOC, nutrients, and other water quality measures 13,14,47 . The high DOC percentage increase in post-snowmelt potentially indicates that DOC was produced in the warm, dry post-snowmelt period, stored in soil via sorption, and flushed out in the next snowmelt, as suggested in Wen et al. 47 . The concerted increase in Fe and DOC, albeit in different years, may indicate a linkage between Fe and DOC mobilization. Destabilization or reductive dissolution of iron was suggested as potentially mobilizing sorbed solutes such as DOC and PO 4 in ten mountainous sites 48 . One may argue that the wet conditions in 2017 drove DOC production. However, snowpack and precipitation have been observed at similar levels in earlier years (for example, 2005 and 2008). These similarly wet years did not produce as much DOC as observed in 2017 and 2018. The annual average concentrations (Fig. 6), although with only four data points, were much higher in high T years, especially in 2018 with average annual daily T min surpassing 0°C. This alludes to the stronger influence of temperature for DOC increases compared to discharge. Whether temperature or precipitation is the major driver of organic matter decomposition is contentiously debated in literature 25,49 . The sensitivity of organic matter (OM) decomposition to temperature and hydrological conditions has been supported with extensive evidence of soil carbon loss and increased soil respiration in warming experiments [50][51][52] .
DOC increase has also been observed worldwide in recent decades, especially in Europe, North America, and Asia [23][24][25] . These places, however, typically observe a gradual increase, instead of the abrupt increase observed in Coal Creek. Multiple hypotheses have been put forth to explain the widespread DOC increase, including recovery from acid rain 23 , changes in land use, and climate change 25 . Coal Creek did not experience a land-use change in the past decades. One may also associate stream DOC increase with tree mortality. Various areas in Colorado have experienced the bark beetle epidemic under warm and drought conditions 53,54 . Their impacts on water quality however are typically significant in areas with >40% infested area, which is not the case in Coal Creek (US Forest Service; personal communication). Even in highly impacted areas, total organic carbon has exhibited a gradual increase over the decades, instead of abrupt escalation within a short period of time 54 .
Geogenic solutes remaining relatively constant: is there influence of acid deposition? Solutes derived from rock weathering typically exhibit dilution C-Q patterns, with high concentrations under low discharge condition reflecting mostly groundwater composition, and low concentrations under high discharge conditions signaling soil water composition 35 . The relatively similar inter-annual dynamics of geogenic solutes (e.g., Ca, Mg, Na) from rock weathering in recent years may indicate that weathering minerals are not as influenced by hydroclimatic conditions because reactive, weathering minerals residing in deeper subsurface (e.g., bedrock) are not as easily influenced by surface temperature and hydrology conditions 55,56 .
Sulfate can come from both rock weathering and acid deposition. Its dilution pattern in concentration-discharge relationship, however, suggests higher concentrations in groundwater become diluted by shallow soil water during snowmelt. This indicates rock weathering is a more dominant source compared to atmospheric deposition. The annual averages of sulfate exhibited an increasing temporal trend that is consistent with observations in a nearby Colorado mountain site 21 and other alpine waters [18][19][20] , pointing to the possibility of accelerated sulfate mobilization from pyrite oxidative dissolution as influenced by lower groundwater table under warming conditions 21 . Ca and Mg, however, did not exhibit a paralleled increase as sulfate, as observed in other mountain sites.
Some areas in Colorado have experienced changes in water chemistry as a result of recovery from acid deposition 57 . In Coal Creek, pH has remained within a relatively stable range between 7 and 8.5 (Supplementary Fig. 3). Recovery from acid precipitation however may not express as increasing pH but as decreasing ratios of (Ca + Mg) to alkalinity, as watersheds gradually lose sulfuric acid as an additional acid source for weathering 58 . The (Ca + Mg)/DIC ratios from 2016 to 2019 do not indicate a strong correlation with sulfate but decrease with discharge. This is possibly due to higher proportions of water from shallow soil zone with low Ca and Mg concentrations in wet years. Analysis for long-term trend line will require long-term data of sulfate, Ca, Mg, and alkalinity or DIC. In Coal Creek, the lack of these longterm data prevents such analysis such that we cannot fully rule out the possibility of acid rain impacts.
Quantifying the long-term trend: a call for establishing HEMs as sentinel sites for consistent and persistent monitoring. The contrasting variations across different solutes accentuate the complex response of Earth's surface processes to warming. Comparison analysis for long-term and short-term data does not lead to decisive conclusions about the long-term response of stream chemistry to environmental conditions in this work. As highlighted in Fig. 7, the LOADEST model developed using longterm data captured the dynamics and concentration range exhibited in these data, and projected the solute dynamics in very similar patterns as in long-term data. In other words, models can only be as good as data. It is important to keep in mind that the infrequent, bimonthly data typically recorded only one-or twopoints during snowmelt, often missing the snowmelt peaks where the highest concentrations of DOC and lowest concentrations of geogenic solutes occur. The increase in Fe and Al in 2016 potentially indicates that variations observed in recent years have already occurred in earlier years. But low-frequency measurements cannot capture these dynamics. As such we cannot draw unequivocal conclusions about the long-term trend.
The concerted escalation of biogenic solute concentrations in topsoils may collaboratively suggest changes occur fastest and earliest in the top thin veneer of the Earth surface with abundant stores of OM. These marked increases may be also an early sign of a tipping point of climate change, above which the rates of soil carbon decomposition may increase substantially. However, this message will remain equivocal with the lack of long-term data at consistent, sufficient frequency. This work underscores therefore the needs for long-term data with sufficient sampling frequency that capture the full dynamics of stream solute variation, in particular under temperature and flow conditions where rapid changes occur. Such data can rule out the possibility of temporal aliasing, differentiate multiple mechanisms, and enable hypothesis testing.
Many HEMs are similarly in remote areas with limited accessibility for temporally dense sampling as in Coal Creek, such that we cannot understand the full magnitude and impact of climatic shifts and environmental perturbations. The present time is often labeled as an era of "big data". The past decades have witnessed rapid advances in technology and an unprecedented generation of earth surface data 59 . Detailed, high-resolution earth surface characteristics data have now become available at unprecedented rates and quantities 60 . The availability of longterm, consistent observations that reflect the functioning of earth's surface (e.g., water chemistry), however, remains to be the bottleneck for understanding earth system response to environmental perturbations. Large collaborative research networks have galvanized to collect consistent data for cross-site comparison 61,62 . Most of the available data however are still in easy-to-access places. Places such as HEMs that experience rapid changes are in dire need of earth system response data to understand ongoing changes 63 . In places such as tropical glaciers in the Andes, the issue of data limitation is even more exacerbated.
Given the sensitivity of HEMs to warming, their water quality response to warming can offer early glimpses of climate change impacts. It is essential to establish watersheds such as Coal Creek as sentinel sites, defined as sites "charged to guard or keeps watch over an area and sounds an alarm if danger comes near." 64 . Extensive and multiple-perspective climate monitoring systems are urgently needed in HEMs in order to record long-term alteration of water quantity and quality in changing climate. Such monitoring has historically recorded many human-induced perturbations such as the discovery of ozone depletion and the response to acid rain 65 . As the pace of climate change accelerates, we will need these sentinel sites that similarly identify early signals of thresholds and tipping points of climate change.

Methods
Field site. Coal Creek (~53 km 2 ) is a representative high elevation watershed in the central Rocky Mountains of Colorado with an average slope of 16°(Supplementary Fig. 1). The mean annual minimum and maximum daily temperatures are −1.5 and 9.0°C, respectively. Annual rainfall and snowfall were~612 mm and 551 cm, respectively, in 2016. The watershed is snow-covered from approximately November to June. The watershed consists primarily of evergreen forest (65%) and herbaceous vegetation (20%), followed by deciduous forest (9%), barren land (3%), and woody wetland (3%). The lithology includes primarily sandstone (39%) and mudstone (15%) that belong to the Mesaverde, Tertiary Wasatch, and Ohio Creek Formations 66,67 . Plutonic rock (15%) originated during the Middle Tertiary. A small area (10%) in the valley is covered by Quaternary Glacial Drift, which consists of surficial deposits including unsorted glacial till and associated sand and gravel deposits. The site was mined for metals between 1874 and 1974 at three primary mines ( Supplementary Fig. 1c) 68 . Although mining has ceased, heavy metals continue to flush into Coal Creek. Coal Creek is part of the designated testbed of Lawrence Berkeley National Laboratory Watershed Function Science focus area (http://watershed.lbl.gov) that explores how mountainous watersheds store and release water and solutes across perturbation gradients, and represents an active collaboration between local, federal partners and university partners.
SNOTEL climate data. Precipitation, snow water (e.g., snowpack depth, snowwater-equivalent (SWE)), and temperature data were retrieved from the United States Department of Agriculture (USDA) Snow Telemetry (SNOTEL) database (CO #380, elevation 3100 m). Daily minimum and maximum temperature data were corrected for sensor bias following Oyler et al. 26 . Warming trends were developed for calendar years (Jan 1-Dec 31). Correlation statistics of temperature to streamflow were on a water year basis (Oct 1-Sept 30). Annual precipitation and SWE data are from 1981 to the present.
Water chemistry and discharge data. Stream chemistry has two sets of measurements. One is the long-term low-frequency measurements (LTM) by the United States Geological Survey (USGS) from November 2000 to November 2019, measured at the very outlet of the watershed (Supplementary Fig. 1). The other is the high-frequency, every-other-day measurements (HFM, a total of >40 solutes) from April to October and every week from November to March from 2016 to 2019 35 , measured at Coal-11, about 1/3 of the creek length into the watershed from the outlet (Supplementary Fig. 1). The sampling location was chosen to reflect watershed chemical conditions. It does not coincide with the USGS sampling point because the USGS sampling point (with long-term data) has a direct anthropogenic addition of Ca and Mg from the wastewater treatment plant discharge. The USGS daily streamflow data are for during ice-free periods (April to November) starting from October 2014 (site ID: 09111250). A regression was developed using a downstream USGS site (ID 385106106571000) to infill winter discharge and extrapolate the record back to Oct 1, 2006. The regression was done in two-parts to integrate low flow data from 2018 (flow ≤ 0.42 m 3 /s, r 2 = 0.87; flow >0.42 m 3 /s, r 2 = 0.92; RMSE for all flows = 5%). All water data were retrieved from the website (https://waterdata.usgs.gov/usa/nwis/).
Water chemistry analysis. Water samples for geochemical analysis were collected and filtered in the field using 0.45 μm Millipore filters. Anions samples were collected in no-headspace 2 mL polypropylene vials. DIC samples were collected in no-headspace 40 mL glass vials with polypropylene open-top caps and butyl rubber septa (VWR® TraceClean®). Cations samples were collected in high-density polyethylene 20 mL vials and acidified to 2% nitric acid with ultra-pure concentrated nitric acid. The samples were transported to the laboratory on ice and stored in 4 o C refrigerator until analysis. Anion concentrations were measured using a Dionex ICS-2100 Ion Chromatography (IC) system (ThermoScientific, USA), while cation concentrations were analyzed using an inductively coupled plasma mass spectrometry (ICP-MS) (Elan DRC II, PerkinElmer SCIEX, USA) 69 . DIC concentrations were measured using a TOC-VCPH analyzer (Shimadzu Corporation, Japan). Total dissolved nitrogen (TDN) was analyzed using a Shimadzu Total Nitrogen Module (TNM-1) combined with the TOC-VCSH analyzer (Shimadzu Corporation, Japan). TNM-1 is a non-specific measurement of TN. All nitrogen species in samples are combusted to nitrogen monoxide and nitrogen dioxide, then reacted with ozone to form an excited state of nitrogen dioxide. Upon returning to the ground state, light energy is emitted. Then, TN is measured using a chemiluminescence detector.
The method detection limits (MDLs) for ICP-MS and DIC and TDN are determined using the US EPA recommended method. Definition and procedures for the determination of the method detection limit, Revision 2 70 . The MDLs for anions are cited from Thermo Fisher Scientific Application Note 154: Determination of inorganic anions in environmental waters using a hydroxideselective column 71 . Generally, the RSD or uncertainty for ICP-MS is <3% based on 5 replicate measurements for concentrations higher than MDLs. The DIC and DOC have an RSD < 3% for concentrations higher than MDLs based on 3-5 measurements. For measurement of anions (Cl − , NO 3 − , and SO 4 2− ), it generally has an RSD < 5% for concentrations higher than MDLs.
Annual export calculation. Daily loads (i.e., mass/time) were calculated as the product of daily discharge and solute concentration (Fig. 4). For the recent 4 years (2016-2019), the annual export was calculated as the area under the load curve (from Jan to December) using the integral trapz function in Matlab. For the sparse LTM data, daily loads from different years were put together in the same year with their corresponding dates to maximize data points for LTM load estimation.
Extrapolation for long-term trend. As an attempt to examine the long-term changes with consistent sampling frequency, we used a combination of a machine learning approach and the USGS tool LOADEST for two representative solutes, DOC and Mg. The bimonthly data were used to avoid sampling frequency bias. The discharge and DOC data are from 2006 to 2019. To infer the long record to 2000, we used the machine learning statistical approach Gaussian Process Regression (GPR) 72 , an approach commonly used to infer discharge data [73][74][75] . We trained a discharge model using the odd years of climate and discharge data from 2006 to 2019, and tested the model using even year data. The GPR uses time-series of precipitation, radiation, temperature, pressure, wind speed, SNOTEL data, and discharge data as inputs to train the model (Matlab R2018a, Machine Learning Regression Toolbox, 5-fold cross-validation), and output missing discharge data in 2000-2006. The Nash-Sutcliffe Efficiency (NSE) values for the training and testing were 0.9 and 0.7, respectively, both are higher than the NSE satisfactory criteria of 0.5 and therefore are satisfactory.
With the full record of 2000-2019 discharge data and model output, the LOAD ESTimator (LOADEST, https://water.usgs.gov/software/loadest/) was used to develop temporal relationships between discharge and concentration for two representative species of DOC and Mg. LOADEST used decimal date, discharge, and concentration data to build linear regressions and automatically selected the best regression model. The training process used the data from 2006 to 2015 and then tested the model using 2016-2019 USGS bimonthly data with the same sampling frequency. The test models have a satisfactory performance (NSE > 0.8). The model was then used to infer concentrations in 2000-2005 with the provided decimal date and daily discharge. Continuous temporal predictions of DOC and Mg that were then inferred from the regression model from 2000 to 2019 (gray lines in Fig. 7).