Machine Learning Deciphers CO2 Sequestration and Subsurface Flowpaths from Stream Chemistry

Abstract. Endmember mixing analysis (EMMA) is often used by hydrogeochemists to interpret the sources of 10 stream solutes, but variations in stream concentrations and discharges remain difficult to explain. We discovered that machine learning can be used to reveal patterns in stream chemistry that pertain to information about weathering sources of solutes and also about subsurface groundwater flowpaths. The investigation has implications, in turn, for the balance of CO2 in the atmosphere. For example, CO2-driven weathering of silicate minerals removes carbon from the atmosphere over ~106-yr timescales. Weathering of another common mineral, pyrite, releases sulfuric acid that in 15 turn causes dissolution of carbonates. In that process, however, CO2 is released instead of sequestered from the atmosphere. Thus, to understand long-term global CO2 sequestration by weathering requires quantification of CO2versus H2SO4-driven reactions. Most researchers estimate such weathering fluxes from stream chemistry but interpreting the reactant minerals and acids dissolved in streams has been fraught with difficulty. We use a new machine learning technique in three watersheds to determine the minerals dissolved by each acid. The results show 20 that the watersheds continuously or intermittently sequester CO2 but the extent of CO2 drawdown is diminished in areas heavily affected by acid rain. Sulfide oxidation contributes ~23% to 62% of sulfate fluxes. Without the new algorithm to deconvolve the mineral weathering, CO2 drawdown was always overestimated. The new technique, which also elucidated the importance of different subsurface flowpaths and long-timescale changes in the watersheds, should have great utility as a new EMMA for investigating water resources worldwide. 25

apportioning. The machine-learning approach we describe here de-convolves sources of stream chemistry without pre-defining the endmembers. We demonstrate this with data from three well-studied watersheds with different characteristics. The new method discovers the endmember chemistries and, as a result, documents new findings of importance previously undiscovered with the other methods. 80 In the following, we focus first on Shale Hills, an acid rain-impacted shale watershed in Pennsylvania USA with extensive data for water/rock chemistry (Jin et al., 2010;Brantley et al., 2013a;Sullivan et al., 2016). This watershed allows the most complete understanding of solute sources. Although we do not show this here, if we use either of the two previously used models for source attribution, stream chemistry data for Shale Hills either does not separate acid rain and pyrite as a sulfate source (if we use the model of Torres et al., 2016) or yields a proportion for 85 acid rain which is larger than 100% (if we use the model of Burke et al., 2018). As shown below, the NMF model easily defines endmembers and proportions.
We then show the utility of the machine learning method for watersheds where less geological information has been published: we investigate sulfide oxidation in East River and Hubbard Brook catchments. Like Shale Hills, East River is shale-hosted, but it receives little acid rain (Winnick et al., 2017). In contrast, Hubbard Brook has been 90 extensively impacted by acid rain but is underlain by schist and glacial till (Likens et al., 2002). In both cases, NMF successfully determines endmembers and source proportions.

Study Sites
Where previous deconvolutions of endmembers were based on assumptions of the chemistry of dissolving minerals alone, data for watersheds shows that not only does mineral chemistry affect stream chemistry, but the flowpath of 100 the water affects this chemistry (e.g. Brantley et al 2017). We demonstrate this with data from three well-studied watersheds with different characteristics. We focus first on Shale Hills, a small (0.08 km 2 ), acid-rain impacted forested watershed underlain by Rose Hill Shale located in central Pennsylvania, USA (Brantley et al., 2018). The Rose Hill Formation shale contains ~0.14 wt% S as pyrite (FeS2) (Gu et al., 2020).
We then show utility of the method in East River (shale-hosted but it receives little acid rain) and Hubbard 105 Brook (extensively impacted by acid rain but is underlain by schist and glacial till) catchments. Specifically, East River is a large (85 km 2 ), mountainous watershed underlain by Mancos Shale that is located near Gothic, Colorado USA within the Gunnison River basin (Winnick et al., 2017). The Mancos is a black shale that contains ~1.6 wt% S as pyrite (Wan et al., 2019). Lastly, Hubbard Brook (Nezat et al., 2004), located in the White Mountains of New Hampshire USA, consists of a series of nine small (0.14-0.77 km 2 ), forested watersheds underlain by Rangeley

110
Formation metamorphosed shale (schist) and sandstone generally covered by glacial till derived mostly from the Kinsman granodiorite. The schist bedrock contains ~0.2-0.9 wt% S and till contains ~0.1-0.2 wt% S. Again, almost all S is present as iron sulfide (pyrite or pyrrhotite). Both bedrock and till are largely carbonate-free.

Machine Learning Model
To assign the proportion of sulfate in streams to sources, we first bootstrapped measurements to increase data volume 130 and then used a method of blind source separation (Alexandrov and Vesselinov, 2014;Vesselinov et al., 2018)  non-negative matrix factorization (NMF). NMF is unique from previously used methods in that it allows calculation of endmember compositions and mixing proportions simultaneously and does not rely on measurements or assumptions of endmembers a priori ( Fig. 2A; see SM section 1). Specifically, NMF decomposes the n x m matrix of stream sample chemistries, V, into two matrices W and H: Here, V is an n x m matrix whose cell entries are solute concentration ratios, [X]/[SO4 2-], for stream samples. Here n refers to the sampling date, m refers to different solutes X ( = Ca 2+ , Mg 2+ , Na + , K + , Cl -), and brackets refer to concentrations. W is the n x p matrix whose cell entries are proportions, a, for each endmember in each stream 140 sample. Again, n refers to sampling dates, but p is the number of sources of solutes (referred to as endmembers). The proportions refer to the fractions of sulfate in each sample that derive from an individual endmember, where the sum of proportions equals 1 ± 0.05 for each sample. As discussed later, we inferred after running the algorithm for the three watersheds that the endmembers represent different flowpaths in the subsurface. Therefore, these proportions of sulfate are referred to here as shallow, moderately shallow, and deep flowpaths, i.e. ashallow, amoderate, and adeep 145 respectively (only one of our target watersheds revealed all three flowpaths). H is the p x m matrix whose cell entries are the concentration ratios that define the chemical signature of each of the p endmembers. The key to NMF is that these concentration ratios are not determined prior to apportionment but rather are determined from the data itself. In addition, the chemical signatures of each endmember can vary temporally around central tendencies. Because the solution to eq. 1 is non-unique, we run the model 20,000 times, apply a filter to the models, and then calculate the 150 mean and standard deviation of the remaining models for trend and error analysis (see SM section 1; Eq. S1).
The only parameter that must be defined to run NMF a priori is the number of endmembers, p. We used principal component analysis (PCA) to determine the number of components needed to explain >90% of the variance in stream solute ratios, and trained NMF to the bootstrapped data while assuming that number of endmembers.
Machine learning determined the compositions defining the endmembers and the mixing proportions of each 155 endmember in each sample. After running NMF, we interpreted each endmember composition based on geological and watershed knowledge.
Based on the outputs of the NMF model, we calculated the weathering rates of sulfide, carbonate, and silicate minerals in the watersheds. Additionally, we calculated the relative contributions of sulfuric and carbonic acid driving those weathering reactions. For details on the weathering calculations see SM section 2.

Synthetic Dataset
NMF is an algorithm that has been used for many applications (e.g., spectral analysis, email surveillance, cluster analysis) but not, to our knowledge, for stream chemistry (Berry et al., 2007). To exemplify the validity of our modeling approach, we generated a dataset of synthetic stream chemistry versus time and ran it through our NMF ran NMF on the synthetic stream chemistry to determine the mixing proportions (a) and endmember compositions , for the five solutes.

Synthetic Data Model
After generating the synthetic dataset of stream samples, we utilized NMF to determine the mixing proportions and endmember compositions. We then filtered out the poor fitting models (see SM Eq. 4). As described more fully in the 180 SM, this left an average number of valid models of 62 (range: 42-77). The average variance between valid models was <10%. Without any prior information about the system, NMF accurately determined the correct mixing proportions (RMSE = 0.04; R 2 = 0.98; p < 0.001; Fig. 2B) and endmember compositions (RMSE = 0.21; R 2 = 0.99; p < 0.001; Fig. 2C). In effect, the model was able to use patterns in the data to accurately unmix the samples.  CO3)2), and pyrite (FeS2) dissolve in regional groundwaters that flow to the stream (Brantley et al., 2013a;Gu et al., 2020). These groundwaters thus contribute DIC, Ca 2+ , Mg 2+ , and SO4 2into the stream.
Like many catchments, water also flows to the stream in Shale Hills along a much shallower near-surface flowpath, referred to here as interflow (Fig. 3). Interflow occurs along a transiently perched water  (Sullivan et al., 2016;Li et al., 2017). Only one mineral, chlorite ((Fe 2+ 0.40Mg0.15Al0.35)6(Si0.76Al0.24)4O10(OH)8), is observed to begin to weather in the deep groundwater and continue weathering all the way to the surface ( Fig. 3; Gu et al., 2020). Chlorite thus weathers to release Mg 2+ to both interflow and deep groundwater. While most water entering the catchment leaves as interflow without entering deep groundwater, the wide reaction zone observed for chlorite is consistent with vertical infiltration of water to the deeper zone (Brantley et al., 2017). PCA for stream chemistry (2008)(2009)(2010) at Shale Hills revealed two sources of sulfate, and this was used to set up NMF, i.e., p = 2 (Table S1). By comparing the compositions from matrix H determined by NMF to our knowledge of the subsurface (Fig. 3), we interpreted the two endmembers as deep and shallow weathering along the 215 two flowpaths, respectively (Jin et al., 2014;Sullivan et al., 2016). Many lines of evidence back up these endmember attributions. The sulfate in the shallow endmember derives from flow well above the pyrite oxidation front through pyrite-depleted rock and is thus attributed to acid rain, while the sulfate in the deep endmember is attributed mostly to pyrite oxidation. Some sulfate from acid rain may make it to the regional groundwater, but the fraction is small. At Shale Hills, acid rain always contains Cland pyrite 225 oxidation always preferentially dissolves carbonate minerals, giving each flowpath endmember a unique signature.
To test the NMF deconvolution, we compared these attributions to isotopic data. Depletion in isotopically heavy S is observed to correlate with increasing concentrations of pyrite-derived sulfate (Fig. 4A), consistent with depleted d 34 S signatures in pyrite (e.g., -20‰; Killingsworth et al., 2018). In contrast, acid rain shows d 34 S values around +3-5‰ (Bailey et al., 2004), and low sulfate concentrations in stream samples are characterized by d 34 S 230 values within this range. Also, as pyrite oxidizes, the concentration of sulfate increases and the d 34 S values decrease to reflect the inferred composition of pyrite, -9.5‰ to -7.2‰ (Fig. 4A). Finally, Gu et al. (2020) showed that pyrite oxidation drives the calcite dissolution at Shale Hills, and this is corroborated by NMF results showing that samples near calcite equilibrium have the highest pyrite-derived sulfate concentrations (Fig. 4B).
However, the annual flux of acid rain-derived sulfate in the shallow endmember determined from NMF at 235 Shale Hills (Table 1) far exceeds the wet deposition of sulfate during the sampling period (Fig. 4C). But such inconsistencies have been noted elsewhere and attributed to travel-time delays over decades between acid rain input and stream output (Prechtel et al., 2001;Mörth et al., 2005). Fig. 4C thus allows us to calculate a ~19-31-year lag time between input and export of sulfate from the temporally changing acid rain input (see SM section 2.4).
Weathering profiles at Shale Hills, the chemistry of the composition (H) matrix, sulfur isotopes, calcite 240 saturation, and lag in acid rain export all support our interpretation that the two components in the NMF model are shallow and deep flowpaths and that sulfate largely derives from acid rain and pyrite respectively.

255
With these calculations we can use NMF results to elucidate the effect of sequestration or release of CO2 at Shale Hills over millennial timescales. CO2-driven weathering of the silicate minerals chlorite and illite removes carbon from the atmosphere and carries it as DIC in rivers to the ocean where it is buried as carbonate minerals (Figure 1, Table S2). In contrast, calcite and ankerite weathering coupled to pyrite oxidation instead releases CO2 to the atmosphere over those timescales. Additionally, acid rain can interact with silicate minerals competing with CO2 for 260 silicate dissolution. To determine the effect on CO2, we calculate a new parameter which we call the stream CO2 sequestration coefficient, kstream. This coefficient, calculated for a single stream sample with units of meq CO2/meq cation, reveals the moles of CO2 sequestered or released during weathering per cation equivalent: often have non-acid sources of sulfate (e.g., fertilizer, gypsum); thus, zstream is calculated by multiplying the fraction of sulfate associated with sulfuric acid (e.g., adeep + ashallow) by the total sulfate concentration and dividing that by the total base cation equivalents in the sample. This calculation shows that seasonally, Shale Hills switches between net production and net consumption of CO2 (Fig. 5D). Using the weathering reactions described in SM section 2.2, we calculate the CO2 flux and find that annually CO2 dynamics are net-neutral at Shale Hills (Table 1).

275
The switch in consumption and production is related to the dominant seasonal flowpath. In the dry season when water tables are low, the stream water is often dominated by deeper groundwater flow that interacts with the deep pyrite reaction front and has little contribution of acid rain. Even though the dry season at Shale Hills is characterized by higher contributions of pyrite-derived sulfate, the watershed acts predominantly as a sink of CO2 during this time of the year because the drawdown of CO2 from silicate weathering is larger than the efflux of CO2 280 from pyrite-driven weathering of carbonate. In the wet season when water tables are high, the stream is dominated by shallow interflow that does not interact with pyrite but has a large contribution of acid rain. Especially during this wet season, acid rain reduces the CO2 consumption by dissolving silicates that could otherwise have interacted with carbonic acid. This is similar to the ideas described in Kanzaki et al. (2020). The impacts of acid rain are greatest during the wet season, during which time the rain often causes the watershed to become a CO2 source by interacting 285 with the deep carbonate minerals and reducing the CO2 consumption by silicate minerals.
To test the accuracy of these inferences based on NMF, we compare to previous results for Shale Hills.
Based on soil pore-water chemistry and rain fluxes at Shale Hills, Jin et al. (2014) estimated the CO2 drawdown from silicate weathering to be 44 mmol m -2 yr -1 . We find that if we assume all silicate weathering is CO2-driven, then the silicate weathering drawdown is 38 mmol m -2 yr -1 , which is consistent with the estimate of Jin et al (2014). But this 290 value is an overestimate that ignored pyrite and acid rain. When we correct this value for H2SO4-weathering of carbonate and silicate derived from pyrite and acid rain, we find that Shale Hills is actually neutral when averaged over annual timeframes and does not sequester CO2 on net (Table 1).

East River
Shale Hills is unique in that it is a monolithologic catchment and the data volume to constrain endmember apportionment is large. But NMF also works well for watersheds in which the subsurface flow and reactions are less constrained partly due to the more complex subsurface geology. The weathering profile at East River (underlain by black shale) shows that pyrite and carbonate are depleted in upper layers but start dissolving at ~2-4mbls (Wan et al.,310 2019). PCA shows that the number of components is 2. The composition of the endmembers for East River are similar to Shale Hills (Table S1); however, the endmember composition indicates a higher proportion of carbonate dissolution by sulfuric acid (see SM section 2).
Based on NMF for East River, pyrite contributes 62% of the annual sulfate flux (Table 1). Sulfuric acid drives 29% to 69% of carbonate dissolution depending on the season, and this compares well with previous estimates 315 of 35-75% (Winnick et al., 2017). Unlike Shale Hills, pyrite oxidation at East River is the dominant source of sulfate because acid rain is less important, and the black shale is pyrite-rich (Fig. 5B).
Both Shale Hills and East River intermittently switch between acting as a source or sink of CO2 (Fig. 5), and both are net neutral within error with respect to CO2 on an annual basis. Interestingly, the seasonality of the switch between Shale Hills and East River is reversed. During baseflow (i.e., flow sustaining stream in between periods of 320 precipitation), Shale Hills is predominantly a sink of CO2, but sometimes it switches to a source of CO2 in the wet season because of the impacts of acid rain. Without the large acid rain influx, East River instead acts as a sink of CO2 during the wet season of snowmelt and then switches to a source during baseflow. This reflects that most of the https://doi.org/10.5194/hess-2020-537 Preprint. Discussion started: 27 October 2020 c Author(s) 2020. CC BY 4.0 License.
sulfate at East River is derived from the pyrite-rich shale (i.e., 65%) and that H2SO4-weathering of carbonates is a more dominant reaction at East River than at Shale Hills. Our results are consistent with previous interpretations 325 (Winnick et al., 2017) suggesting CO2 efflux rates are highest in baseflow-dominated and lowest in snowmeltdominated flow regimes. We agree with the interpretations of Winnick et al. (2017), but NMF adds the new insight that on an annual basis, the efflux of CO2 from H2SO4-weathering of carbonates is balanced by CO2-weathering of silicates making East River net neutral with respect to CO2.

Hubbard Brook
Monolithologic shale watersheds are not the only target chemistries that can be deconvoluted with NMF. Large variations in the d 34 S composition of the bedrock at Hubbard Brook (Bailey et al., 2004) mean that sulfur isotopes in stream water are difficult to unambiguously apportion sulfate sources. Weathering fluxes from pyrrhotite (FeS), the sulfide mineral known to be present, have therefore never been fully constrained (Mitchell et al., 2001).

335
At Hubbard Brook, PCA shows three endmember sources of sulfate. As described below, we attribute these to three inferred flow lines, two in till and one at depth: waters flowing through i) shallow soil developed from till, ii) moderately-deep, less-weathered till, and iii) weathering bedrock. A three-layered weathering profile has been observed in other till-covered areas of New Hampshire as well (Goldthwait and Kruger, 1938). We used these ideas to identify endmembers as described below.

340
Concentrations of sulfate in acid rain have declined over time in northeastern USA (Lynch et al., 2000;Lehmann et al., 2007). Of the three NMF-determined endmembers at Hubbard Brook, two of them show declining sulfate concentrations with time. We therefore attributed the first and second endmembers to acid rain (Fig. S1).
Only one endmember showed little to no decline in sulfate concentration over time, and we therefore attributed that endmember to deep weathering in water interacting with the underlying bedrock. Moreover, the metamorphic conditions that produce pyrrhotite also produce biotite and chlorite, and those three 350 minerals tend to be located together in schist foliations (Carpenter, 1974). Pyrrhotite oxidation at Hubbard Brook apparently causes dissolution of biotite + chlorite because these are the most susceptible minerals in close proximity to the sulfide. Thus, several lines of evidence underlie our interpretation that component 3 is the deep weathering source of sulfate.
As summarized in Table 1, pyrrhotite can account for 29% of the total sulfate flux at Hubbard Brook. The 355 schist and till contain essentially no carbonate; therefore, weathering is always a net sink for CO2. In this watershed, however, the story is complicated by the sulfuric acid from pyrrhotite oxidation and acid rain that dissolves silicate minerals. If we had assumed all of the base cations detected in Hubbard Brook were caused by CO2-weathering, we would have overestimated the net drawdown of CO2 out of the atmosphere (Fig. 1). NMF deconvolutes this signal. into account the fraction of base cations weathered, the fraction of base cations from carbonates, and the capacity of the bedrock to produce H2SO4: Here, krock is the time-integrated CO2 sequestration coefficient determined from the solid phase weathering products

Conclusions
By not requiring a priori assignments of endmembers, a new machine learning technique not only successfully reproduced source apportionments made in more traditional endmember analysis for streams, but it also revealed new information about how watersheds work. At the same time, the method also solved some issues related to source apportionment for streams with time variations of large acid rain inputs. The approach documented that two 415 https://doi.org/10.5194/hess-2020-537 Preprint. Discussion started: 27 October 2020 c Author(s) 2020. CC BY 4.0 License.
carbonate-containing shale watersheds (Shale Hills, East River) act intermittently as sources or sinks of CO2 to the atmosphere but on net are neutral with respect to CO2. In contrast, because it has no carbonate minerals, Hubbard Brook is a constant sink for CO2 (Fig. 5). These observations were compared and confirmed by comparing stream chemistry to rock chemistry.
NMF also emphasized the importance of different water flowpaths in determining endmembers: the 420 endmembers were not strictly defined by mineralogy but by patterns of subsurface flow. These flowpaths lead to patterns in stream water chemistry that were easily deciphered by our newly developed machine learning-based mixing model. In particular, for three streams, signals in the chemical variations were observed to reveal dissolution of the most reactive mineral in proximity to sulfide oxidation. Many watersheds have flowpaths distinguished by geochemical signatures from mineral reactions (Brantley et al., 2017) but we do not know these paths a priori when 425 we investigate stream chemistry. Machine learning will be useful to model mineral reactions on broader spatial scales and will help constrain global weathering-related CO2 dynamics because it can delineate endmembers without a priori assumptions.
Beyond these attributes, the machine learning approach also revealed other new attributes of weathering. In Shale Hills, we discovered that sulfate inputs from acid rain are not exported completely for two decades, which 430 impacts mass balance and weathering-related CO2 dynamics. Although not discussed explicitly here, this decadal time-lag was also observed at Hubbard Brook. NMF also showed that Hubbard Brook, recovering from the impacts of acid rain, is only recently returning to its full potential as a CO2-sequestering rock system. In other words, prior to acid rain, Hubbard Brook sequestered more CO2 per mole of weathered bedrock than it does today. But acid rain dissolved some of the silicates with H2SO4, lowering the CO2 sequestration capability of the watershed. NMF led us 435 to discover this new attribute of acid rain, namely that it diminishes the capacity of a rock to sequester CO2 at millennial timescales ( Figure 1) by replacement of CO2 by H2SO4 as a weathering agent. Regardless of the net CO2 dynamic, we discovered that without considering sulfide oxidation or acid rain, the CO2 weathering sink considered over millennial timescales is always overestimated.

Data Availability
Data used in analysis for this work can be found at the Susquehanna Shale Hills Critical Zone website (https://criticalzone.org/shale-hills/data), the Hubbard Brook data catalog (https://hubbardbrook.org/d/hubbard-455 brook-data-catalog), and in Wan et al. (2019) and in Winnick et al. (2017). Codes used in the modeling are available upon request.

Author Contributions
SLB, TW, and ARS conceived the project. TW and ARS developed the codes for the model. SLB, XG, and ARS

460
interpreted the data and model outputs. All authors contributed in the preparation of the manuscript.