Climate Controls on River Chemistry

How does climate control river chemistry? Existing literature has examined extensively the response of river chemistry to short‐term weather conditions from event to seasonal scales. Patterns and drivers of long‐term, baseline river chemistry have remained poorly understood. Here we compile and analyze chemistry data from 506 minimally impacted rivers (412,801 data points) in the contiguous United States (CAMELS‐Chem) to identify patterns and drivers of river chemistry. Despite distinct sources and diverse reaction characteristics, a universal pattern emerges for 16 major solutes at the continental scale. Their long‐term mean concentrations (Cm) decrease with mean discharge (Qm), with elevated concentrations in arid climates and lower concentrations in humid climates, indicating overwhelming regulation by climate compared to local Critical Zone characteristics such as lithology and topography. To understand the CmQm pattern, a parsimonious watershed reactor model was solved by bringing together hydrology (storage–discharge relationship) and biogeochemical reaction theories from traditionally separate disciplines. The derivation of long‐term, steady state solutions lead to a power law form of CmQm relationships. The model illuminates two competing processes that determine mean solute concentrations: solute production by subsurface biogeochemical and chemical weathering reactions, and solute export (or removal) by mean discharge, the water flushing capacity dictated by climate and vegetation. In other words, watersheds function primarily as reactors that produce and accumulate solutes in arid climates, and as transporters that export solutes in humid climates. With space‐for‐time substitution, these results indicate that in places where river discharge dwindles in a warming climate, solute concentrations will elevate even without human perturbation, threatening water quality and aquatic ecosystems. Water quality deterioration therefore should be considered in the global calculation of future climate risks.

Despite its importance, the risks of changing river chemistry and degrading water quality under a warming climate have received marginal attention. An example is the 12-chapter, 1,300-page Intergovernmental Panel on Climate Change Assessment Report 6 (IPCC, 2021). The report discusses changing water chemistry in oceans, including ocean acidification and deoxygenation. The risks of changing chemistry of inland waters that sustain all aquatic and terrestrial lives, including human, have barely been discussed. This is in stark contrast to multiple chapters in IPCC report and extensive literature on the impacts of changing climate on water cycles, climate disasters, and hydrological extremes (AghaKouchak et al., 2021;Bales et al., 2006;Higuera & Abatzoglou, 2021;IPCC, 2021;Paschalis et al., 2020).
Such insubstantial attention on water quality risks may arise from limited understanding on how river chemistry evolves under a changing climate. The common perception is that river chemistry is regulated more by local material abundance instead of climate. For example, concentrations of geogenic solutes (e.g., elements released from rock weathering) are conceived as depending on lithology (Bluth & Kump, 1994;Gaillardet et al., 1999;Ibarra et al., 2016). In contrast, concentrations of biogenic solutes involved in biological processes (e.g., nutrients and carbon), are thought to hinge on vegetation, land cover, and anthropogenic input (Kim et al., 2020;Van Meter et al., 2018;Worrall & Burt, 2007). Water originates from precipitation and travels via subsurface flow paths and river corridors before exiting at river outlets. Along its flow paths, water mobilizes solutes by interacting with roots, microbe, soil, and rocks. The dynamics of solute concentrations in the stored water (V w ) are therefore regulated primarily by two competing processes: the addition of solutes by external input (I) and production of solutes by reactions (R) in soils, rocks, and streams, and the export of solutes by discharge (QC) (Equation 2). External input can happen but are often insignificant in natural, minimally-impacted watersheds. Under warming climate, some places will become more arid, leading to lower river discharge (Q = P -ET) and solute export (QC); some places will become more humid, resulting in lower solute concentrations but higher solute export fluxes to rivers. Watersheds function primarily as solute-producing reactors in arid, warm climates, and as solute-exporting transporters in humid, cold climates.

10.1029/2021EF002603
3 of 17 Existing literature has extensively documented the response of river chemistry to weather conditions in streams and rivers at event to seasonal time scales, often encapsulated in concentration-discharge relationships (Dupas et al., 2016;Ebeling et al., 2021;Godsey et al., 2019;Godsey et al., 2009;Hooper et al., 1990;Knapp et al., 2020;Moatar et al., 2017;Pinder & Jones, 1969). Geogenic solutes often exhibit dilution patterns where concentrations decrease with discharge or chemostatic patterns where concentrations vary negligibly compared to discharge (Godsey et al., 2009;Ibarra et al., 2016;Torres & Baronas, 2021). On the other hand, biogenic solutes such as DOC often exhibit flushing patterns where concentrations increase with discharge (Boyer et al., 1997;Herndon et al., 2015;Zarnetske et al., 2018). The response of nitrate concentrations to discharge varies with land uses: flushing patterns emerge in most agriculture sites, whereas diverse patterns have been observed in natural rivers with minimal human influence (Botter et al., 2020;Ebeling et al., 2021;Moatar et al., 2020;. Mounting evidence has shown that these contrasting patterns arise from the influence of changing weather, flow paths, and distinct source water chemistry in shallow soils and in deeper groundwater Zhi et al., 2019;Zimmer & McGlynn, 2018).
Long-term river chemistry response to different climate conditions however has been scarcely studied. Riverine mean concentrations of geogenic solutes have been observed to decrease from arid to humid climates (Maher & Chamberlain, 2014;White & Blum, 1995).These concentration variations have been explained by the thermodynamic limits of chemical weathering and the approach to equilibrium at long water travel times under arid climates (Ibarra et al., 2016;Maher & Chamberlain, 2014). More recently, a wide range of solutes has been shown to similarly exhibit dilution patterns across climate gradients, with high concentrations under arid climates (Godsey et al., 2019). This has been postulated to depend on the time scales of soil buffering. In general, however, we lack conceptual framework and quantitative approaches that can mechanistically explain and project baseline concentrations for solutes of diverse origins.
How does climate control river chemistry? With space-for-time substitution, answers to this question can help project river chemistry in the future climate. Our aim is to identify patterns across climate gradients and understand predominant control of river chemistry at the continental scale. We compile and analyze chemistry data from 506 minimally-impacted rivers (412,801 data points) in the contiguous United States. In doing so we focus on the effects of climate without the entangled effects of human perturbation. We calculate long-term, baseline concentrations of pervasive solutes of diverse origins for each river. We further derive solutions to a parsimonious watershed hydro-biogeochemical reactor model to explain the observed patterns and quantitatively depict the relationship between long-term mean concentrations and mean specific discharge, a key measure of climate conditions.

The CAMELS-Chem Database
This work compiled a new dataset CAMELS-Chem. The new dataset augments the existing CAMELS (Catchment Attributes and Meteorology for Large-sample Studies) dataset (Addor et al., 2017), and comprises the U.S. Geological Survey chemistry data (412,801 data points, 1980 to 2014) from 506 minimally-impacted rivers (Sterle et al., 2022). It pairs atmospheric deposition and water chemistry data with the existing CAMELS (Catchment Attributes and Meteorology for Large-sample Studies) dataset, therefore filling the need of a data set with complementary watershed attributes, river water flow, and river chemistry data. The dataset includes 18 common river chemistry constituents: Al, Ca, Cl, DOC, Total Organic Carbon (TOC), HCO 3 , K, Mg, Na, Total Dissolved Nitrogen [nitrate + nitrite + ammonia + organic-N], total organic nitrogen (TON), nitrate (NO 3 ), dissolved oxygen (DO), pH (field), pH_l (lab), Si, SO 4 , and water temperature. This work used data of all solutes except DO, a solute that depends largely on river temperature, light, and flow dynamics (Bernhardt et al., 2018;Zhi et al., 2021), and pH_l (lab) that is similar to pH (field).
The solutes are loosely categorized based on their origins into biogenic (from soil biogeochemical reactions) and geogenic groups (from chemical weathering). The solute HCO 3 can be biogenic from soil respiration (that produces soil CO 2 ) and geogenic from carbonate weathering. Its concentration limits however are controlled by the thermodynamics of carbonate minerals and is therefore grouped into geogenic solutes.
The mean concentrations (C m ) of each solute were calculated as the arithmetic mean of available concentration data at each site in CAMELS-Chem. We purposely choose to use arithmetic mean, instead of the commonly used 10.1029/2021EF002603 4 of 17 flux-normalized mean concentrations (Godsey et al., 2019;Ibarra et al., 2016;White & Blum, 1995), because arithmetic mean gives equal weight to every data point. Flux-normalized mean concentrations are weighted by fluxes (discharge) that often vary by orders of magnitude, and therefore give more weight to concentrations at high discharge. This could lead to more biased mean concentrations, especially when the timing and frequency of sampling are inconsistent across different discharge regimes at different sites.

The Long-Term Energy and Water Balance of Watersheds
The amount of precipitated water (in both rainfall and snowfall forms) that eventually ends up in streams and rivers, or mean discharge, is influential in regulating water travel time and the extent of interactions among water, soil, and rocks, and therefore river chemistry (Ibarra et al., 2019;Keller, 2019;Li et al., 2021). Here we relate mean solute concentrations to climate conditions at each site using a widely used relationship, Budyko equation. Various forms of Budyko equations relate the long-term ratio of evaporative index to aridity index (Chen & Sivapalan, 2020;Gentine et al., 2012;Reaver et al., 2020). Here we used the original form (Budyko, 1974): where P is precipitation, ET is evapotranspiration, the summation of evaporation (from open water, bare soil, and vegetated surfaces), transpiration (from within plant leaves), and sublimation from ice and snow surfaces. Evapotranspiration is regulated by the land-surface interactions that determine the exchange of energy and water between land surface and atmosphere. The climatic potential ET (PET) measures the "drying power" regulated by both climate and land cover (Dingman, 2015). The dimensionless parameter β modifies the curvature of the evaporative index ET/P as a function of aridity index PET/P.
The CAMELS-Chem dataset includes PET, mean annual precipitation (P), and mean annual specific discharge (Q m, thereafter mean discharge) for each site. Actual evapotranspiration (ET) was calculated as the difference between P and Q m , assuming water-balanced watersheds.

The Watershed Hydro-Biogeochemical Reactor Model
Broadly, we can consider mean river chemistry as reflecting the response of long-term, baseline chemistry to climate, and their instantaneous variations from event to seasonal scales as reflecting its response to rapidly changing weather conditions. We can conceptualize watersheds as well-mixed hydro-biogeochemical reactors ( Figure 1). Solutes can come from external input such as atmospheric deposition and/or human input such as the addition of nutrients from agriculture or urban watersheds. As meteoric water infiltrates and travels through subsurface and river corridors, it mobilizes solutes by interacting with materials along its flow paths (at rates R), and ultimately exports solutes at the rates QC.

The governing equation and solutions.
The mass conservation of an arbitrary solute in a well-mixed watershed reactor can be written as follows: Here V w is water volume per drainage area (m 3 water/m 2 drainage area); C is aqueous concentration of the solute (g/m 3 water, equivalent to mg/L water); I is external input rate (g/m 2 drainage area/yr) (orange arrows across the land surface, Figure 1), which can be dry or wet deposition or both in a minimally-impacted watershed; R is net rate of reactions (g/m 2 drainage area/yr, red arrows) that produce or consume the solute; and Q is specific discharge (m 3 water/m 2 drainage area/yr, dark blue arrows in Figure 1). Individual solutes are often involved in multiple reactions. In the simple representation of Equation 2, the single R term may lump multiple reactions and can be considered as the net reaction rate of multiple reactions. For a water-balanced watershed, Q = P -ET, where P and ET are precipitation and evapotranspiration (m 3 water/m 2 drainage area/yr, light blue and green arrows in Figure 1), respectively.

of 17
As detailed in Supporting Information S1, the time-dependent analytical solution to Equation 2 is Here C m is the long-term, steady state, baseline mean concentration in individual sites. The second term [C 0 − C m ] − ∕ is time dependent with an exponential form. The characteristic time τ w = V w /Q m quantifies the mean water travel time from its entry on land to its exit point (river) in a watershed (Rinaldo et al., 2015;Sprenger et al., 2019;Tetzlaff et al., 2009). If time is sufficiently long (large t), the exponential term approaches zero such that the concentration becomes C m = (I m + R m ) / Q m , the steady state solution. For example, when t = τ w , − ∕ = 0.37; when t = 3τ w , − ∕ = 0.05. That is, when the time is longer than a few times of τ w , the second term becomes negligible.
Almost all solutes have varied sources. Solutes such as Cl, SO 4 , and NO 3 may have comparably larger atmospheric input (I) than other solutes. In urban and agricultural sites, the human input of nutrient-containing solutes can be significant. In minimally-impacted watersheds here, the external input however tends to be low compared to productions from reactions in the subsurface (Berner & Berner, 2012). We therefore focus primarily on the reaction rate R as the major source. In that case, the mean concentration has the following form with mean reaction rate R m and mean discharge Q m : Equation 4 reveals that long-term, steady-state mean concentrations are regulated by the relative magnitude of two processes. One is the reactions in the subsurface or in river corridors that produce a solute (at rate R m ), the other is the export process that carries a solute out of the watershed at the river exit (at rate Q m ). High R m and low Q m result in high concentrations, and low R m and high Q m lead to low concentrations in subsurface and by extension in rivers. This equation can be used broadly as an organizing framework to understand competing processes that regulate baseline solute concentrations.
The C m Q m relationships: the response of river chemistry to climate. The values of R m at the watershed scale are typically unknown. We do know the functional dependence of rates on measurable variables such as temperature, soil moisture, material abundance, intrinsic rates, and aqueous geochemistry from reaction kinetic theory. We can combine theories from hydrology and reaction kinetics and thermodynamics from biogeochemistry for different groups of solutes based on their origin and reaction types (details below and in Supporting Information S1) to derive steady state relationships between C m and Q m in the form of Equations 5 and 6.
For biogenic solutes, the derivation of C m Q m relationship involved the use of water storage-discharge relationships from hydrology and reaction kinetics with dependence on temperature and water content (or soil moisture) from biogeochemistry. The reaction rate R m depends on temperature (T) and soil moisture (or water content S w ) (Davidson et al., 1998;Mahecha et al., 2010). The rate dependence on water content takes the general form where k 0 is the maximum reaction rate constant, S A is the surface area that quantifies the abundance of organic materials, f(T) describes the rate dependence on temperature, which can take the form of Q 10 formula (e.g., f(T) = Q 10 (T−20)/10 , Lloyd & Taylor, 1994), and f(S w ) describes the rate dependence on soil moisture. Reactions such as soil respiration and carbon decomposition often reach maxima at intermediate water content (e.g., soil moisture of 0.5-0.7 in many places) (Yan et al., 2018), because rates tend to be low under extremely dry and wet conditions. The mean soil moisture in particular locations and at the watershed scale however rarely become higher than 0.5, even in wet places and at wet times (e.g., Lin et al., 2006). We therefore use a simple form f(S w ) = (S w ) n , where the exponent n quantifies the sensitivity of reaction rates to water content. Water content (or storage) and discharge can be linked by storage-discharge relationships in various forms (Kirchner, 2009). Here we use a simple form = 0 0 (Wittenberg, 1999). Large β 0 values indicate rapidly changing discharge responding to changes in S w , whereas small β 0 values indicate relatively stable discharge that varies insignificantly with S w .
With these reaction rate laws and functional dependence, the following C m Q m relationship can be derived for biogenic solutes (details in Supporting Information S1): The parameter A is a bulk measure of reaction rates and combines several rate-related parameters, including reaction rate constant k 0 , surface area S A , rate dependence on soil moisture n and temperature f(T), and parameters in S w -Q relationships. The parameter B is determined by the ratio of the rate dependence on water content (n) and discharge dependence on water content (β 0 ). The B value quantifies concentration sensitivity to variations in mean discharge. Here the symbols A and B are used to differentiate from the a and b in the instaneous CQ relationships C = aQ b .
For geogenic solutes, weathering reaction kinetics follow the transition state theory with a thermodynamic limit term R m = R 0 (1 − C/C eq ), with the last term indicating the distance from reaction equilibrium (Lasaga, 1998;Lasaga et al., 1994). This rate law leads to the following equation (details in Supporting Information S1): where Da = τ w /τ eq , the dimensionless Damköhler number that quantifies the relative magnitude of water travel time (τ w ) and the time to reach equilibrium τ eq = C eq V w /R 0 . This equation builds upon a rich history of reactive transport formulation in literature for geogenic solutes with thermodynamic limits (Ibarra et al., 2016;Maher, 2011;Maher & Chamberlain, 2014;Salehikhoo et al., 2013;Steefel, 2007;Torres & Baronas, 2021;Wen & Li, 2017;Wen & Li, 2018;Wymore et al., 2017).
Under the conditions that τ w is small (humid climate) or τ eq is large (e.g., for silicate weathering), Da << 1, reaction is far from equilibrium, Equation This equation has the same form as Equation 4 for biogenic solutes and can be further rearranged into the same form of Equation 5 with parameters specific for weathering reactions. For example, weathering of rocks such as silicates, shale, and evaporites releases cations (Na, K, Si, and Al) and anions (Cl, SO 4 , HCO 3 ). These reactions lead to the precipitation of clay, effectively reducing concentrations of solutes from silicates such that they hardly reach reaction equilibrium.
For geogenic solutes involved in carbonate weathering (Ca, Mg, and HCO 3 , and pH), they can rapidly reach equilibrium such that their C m Q m relationships follow Equation 6. For these solutes, when τ w is small (humid climate), Da << 1, C m decreases with increasing Q m and depends on R 0 . Arid climates entails that τ w >> τ eq and Da >> 1, such that Equation 6 becomes C m = C eq .

Model Assumption, Implications, and Utilities
The model conceptualizes watersheds as well-mixed reactors that drain to rivers. It is meant to broadly capture controls of long-term, baseline mean concentrations. It is not meant to represent detailed short-term temporal dynamics and spatial structure such as flow paths and source water chemistry at different subsurface depths.
In-stream processes may also influence concentrations, especially under low flow conditions (Dodds, 2006;Perdrial et al., 2014). In-stream reaction processes are implicitly counted for as in the governing Equation 2, as one can consider in-stream reactions as contributing partially to R m . The parsimonious representation entails a slim number of parameters and salient relationship that reveal first-order mechanisms. As we will illustrate in Section 3, these relationships capture concentration variations across climate gradients and reveal competing mechanisms that ultimately drive river chemistry.
For solutes that are not thermodynamically limited by equilibrium, the derivation here accidently arrives at a power law form (C m = AQ m B ). This C m Q m power law relationship however differs from the short-term, instantaneous CQ patterns (C = aQ b ) at event to seasonal time scales that have been extensively studied in literature, as we will elaborate more in Section 4. It depicts the dependence of long-term, baseline concentration on mean discharge at annual time scales and beyond.
Although this work focuses on minimally-impacted watersheds, the application of the model is not limited to such watersheds. For example, the input term I in Equation 2 is loosely defined: it can be atmospheric input and/ or human input such as the addition of nutrients in agriculture or urban watersheds. The long-term accumulation of nutrients on agriculture lands can also change some of reaction parameters, for example, the surface area S A that quantifies the abundance of nitrogen-and phosphorus-containing materials. The model can be used broadly as an organizing framework to understand competing processes that regulate mean solute concentrations. More complex reactive transport models at the watershed scale can be used to complement such simple models and to examine detailed processes, variables, and spatial heterogeneities that reflect idiosyncrasies of specific sites (Li, 2019;Wen et al., 2021;Wen et al., 2020).

Large Concentration Variations at the Continental Scale
The compiled data indicate that mean discharge varies by three orders of magnitude and is higher in the eastern US and along the coasts, and lowest in the Great Plains from North Dakota down to Texas (Figure 2). Mean concentrations also vary by orders of magnitude across sites, depending on specific solutes. They are highest in the Great Plains for many solutes, including Cl, SO 4 , and cations (e.g., Na, K, Ca, and Mg). Concentrations are typically higher in arid climates and lower in humid climates, although considerable variations occur. Although eastern US has historically received more acid deposition (Berner & Berner, 2012), rivers in the western US generally have higher SO 4 concentrations, possibly due to gypsum from evaporites and lack of dilution due to low discharge. NO 3 concentrations are highest in mid-west Corn Belt area, possibly elevated by dust and dry deposition from nearby agricultural areas.

Biogenic Solutes: Higher Concentrations at Lower Mean Discharge
Biogenic solutes include TOC, DOC, TON and dissolved organic nitrogen (DON), total nitrogen (TN), and nitrate (NO 3 ) (Figure 3). They are primarily involved in vegetation-and microbe-mediated reactions including soil respiration, nitrification, and denitrification. Some solutes, including NO 3 , also come from atmospheric deposition and fertilizer application (Berner & Berner, 2012). The climate conditions of these rivers cover a wide range, as manifested in evaporative (ET/P) and aridity index (PET/P) values (red lines in left figures of each pair, Figure 3). Concentrations decrease as mean discharge increases and are typically higher under arid climates with high ET/P and PET/P (Figure 3, right). Red lines in the right figures of Figure  Values of B vary across solutes, indicating distinct sensitivity to mean discharge variations. For TOC and DOC, B values are close to −0.9; for TON and DON, they are around −0.5 to −0.6. Organic carbon is typically more abundant than nitrogen-containing solutes in natural systems. Carbon-and nitrogen-containing solutes also experience distinct reactions with different rate sensitivity to temperature and water content. These differences could lead to dissimilar B values and sensitivity to mean discharge. At similar discharge, concentrations vary considerably and hinge upon mean reaction rates, as indicated by the red lines with different A values. The model inferred that rates of TOC and DOC vary from 0.1 to 5 g/m 2 /yr, whereas the rates of nitrogen-containing solutes are typically one to two orders of magnitude lower. This is in fact consistent with the observed range of C:N ratios between 1 and 100 (Cleveland & Liptzin, 2007), as carbon is much more abundant than nitrogen in natural soils.

Geogenic Solutes Without Limits: Elevated Concentrations at Higher Aridity
Soil and rock weathering release cations (Na, K, Si, and Al) and anions (Cl, SO 4 , HCO 3 ). Solutes such as Cl, NO 3 , and SO 4 can also cone from rainfall, dust, and pollution sources (e.g., acid rain), although atmospheric deposition is considered low compared to weathering (Berner & Berner, 2012). Weathering of silicates, shale, and evaporites rarely reach equilibrium such that the C m Q m relationships of solutes from these rocks have the same form as The data and model solution indicate that the reaction rates (A values) vary by orders of magnitude, with more abundant solutes having higher rates. For anions, the values of A vary from 0.1 to 10 g/m 2 /yr for Cl and 0.2-25 g/m 2 /yr for SO 4 . Among cations, Na and Al have the highest (0.03-1) and lowest (0.0002-0.1 g/m 2 /yr) A values, respectively. This is expected, as Na is highly soluble and has the highest concentration range of all cations. In contrast, Al released from silicate dissolution is usually immobilized rapidly via clay precipitation. Equation 5 also predicts that B values of Na and K are around −0.6 and −0.5, whereas that of Si is around −0.2. The B values of Cl and SO 4 are close to −1. Cl can come from dissolution of evaporites. SO 4 can originate from atmospheric acid deposition but also from pyrite oxidation in shale and gypsum dissolution (Berner & Berner, 2012). These source rocks dissolve rapidly with high solubility. In particular, pyrite dissolution is known to hinge upon the rise and fall of water tables and the availability of oxygen (Crawford et al., 2019). These characteristics can lead to its high sensitivity to discharge and therefore more negative B values.

Geogenic Solutes From Carbonate Weathering: Higher Concentrations at Higher Aridity but With Limits
Geogenic solutes related to carbonate weathering (Ca, Mg, and HCO 3 , and pH) can reach equilibrium rapidly because of high reaction rates. Their C m Q m relationship (Equation 6) indicates a form of dependence on mean discharge that differs from biogenic and non-thermodynamic-limiting geogenic solutes discussed in the previous section. Instead of monotonic increase with decreasing discharge, values of C m reach maxima at a threshold discharge. Data in Figure 5 confirm the existence of such threshold: at Q m > ∼200 mm/yr, concentrations decrease sharply with discharge. Carbonate as a parent rock mostly occurs in Midwestern US (Moosdorf et al., 2010). However, pedogenic carbonate is known to prevail in dry soils in the western US and under arid climates around the globe (Monger et al., 2015;Zamanian et al., 2016). Concentrations of each solute vary considerably, depending on the reaction rate R 0 and equilibrium constant C eq . Pedogenic carbonates often have different solubility due to dust and impurities (Macpherson & Sullivan, 2019), and therefore can have variable C eq   , C m = C eq R 0 /(Q m C eq + R 0 ), with different equilibrium concentrations C eq (mg/L) and reaction rates R 0 (g/m 2 /yr). Mean concentrations increase with decreasing Q m until reaching maxima at a threshold Q m around 200 mm/ yr. values. Between solutes, HCO 3 has the highest R 0 and C eq , possibly because HCO 3 can also come from biological processes such as soil respiration in addition to carbonate weathering.

Climate Regulation of River Chemistry: Watersheds as Reactors in Arid Climates and Transporters in Humid Climates
This work used a large data set to identify patterns of baseline river chemistry at the contiguous United States. A universal C m Q m pattern emerges: mean concentrations of all 16 solutes increase with aridity and decrease with increasing mean discharge. This pattern contrasts the general perception that river chemistry is primarily controlled by local material abundance in the Critical Zone, including lithology (Gaillardet et al., 1999), vegetation, and organic matter (Raymond et al., 2008;Van Meter et al., 2018;Worrall & Burt, 2007).
A watershed hydro-biogeochemical reactor model was solved analytically to understand the processes that drive the pattern. The derivation of the solution and C m Q m relationships (Equations 5 and 6) brings together theories from traditionally separate disciplines: storage-discharge relationship in hydrology and reaction kinetics in biogeochemistry and chemical weathering. The steady state solution illuminates two competing processes that drive the continental-scale pattern (Equation 4, Figure 1): solute addition via reactions (and possibly input), and solute export by discharge. Equation 4 indicates that in arid climates, reaction rates R m are higher compared to Q m , leading to high C m . In other words, watersheds function as reactors that primarily produce and accumulate solutes. In contrast, in humid climates, reaction rates R m are higher compared to Q m , leading to low C m . Watersheds therefore mostly function as transporters that export solutes.
Such solution has been developed before for geogenic solutes (chemical weathering solutes) (e.g. Maher & Chamberlain, 2014) but not for biogenic solutes such as carbon and nutrient solutes. The dilution patterns of geogenic solutes, where high concentrations under arid climates and lower, rapidly decreasing concentrations under humid climates, echo earlier observations in literature (Maher & Chamberlain, 2014;White & Blum, 1995). This work therefore brings a unified framework for solutes of both biogenic and geogenic origins.
The model underscores the first-order control of mean discharge and the secondary influence of material abundance, as expressed in reaction rate parameters (A or R 0 ). As shown in Table 1, the A values are lowest (0.0002-0.1 g/m 2 /yr) for Al, which happens to have the lowest concentration range. The A values are highest (9-80 g/m 2 /yr) for HCO 3 , which also has the highest concentration range. In other words, high A or R 0 indicate rapid solute production and high concentrations. This is hardly surprising, as these values encapsulate the intertwined influence of reactivity (k 0 ), mineral abundance (S A ), climate (f(T)), and water flow (α 0 ) (Equations 5 and 6). For biogenic solutes, values of A may also indicate the quantity and quality of organic matter (and therefore vegetation) and microbial communities. The lumped nature of A and R 0 however alludes to the challenges of differentiating the influence of individual factors.
The sensitivity of concentrations to discharge variation is quantified by the value of B (= n/β 0 -1), which depends on the reaction rate dependence on water content (n) and discharge dependence on water content (β 0 ). The former (n) likely hinges upon solute origin and chemical properties of materials (e.g., organic matter and minerals) that produce solutes. In particualr, Cl and SO 4 have B values of 01, potentially indicating they are primarily from external atmospheric input in many sites (see Equation S3 for mostly airborne solutes). The latter (β 0 ) depends on the physical structure that regulates water storage and flow generation. For example, steep topography and clay-rich soils respond rapidly to precipitation and have high β 0 values (Kirchner, 2009;Wlostowski et al., 2021;Xiao et al., 2019). When n > β 0 , or reaction rates are more sensitive to water content than discharge, B > 0;

Table 1
Summary of Parameters in Equations 5 and 6 for All Solutes when n < β 0 , or reaction rate is less sensitive to water content than discharge, B < 0. Values of B are negative for all solutes, indicating reaction rates are less sensitive to water content than discharge. This is expected, as discharge predominately depends on water content, whereas reaction rates depend on an array of factors other than water content, which may subdue the impacts of water content. Overall, the expression of B underscores the importance of the subsurface physical and chemical structure in regulating the sensitivity of solute concentrations to discharge variations.
It is important to note that in CAMELS-Chem, there are more sites in the eastern than the western US. Between different sites, measurement duration, frequency, and time periods vary for different solutes. Some biogenic solutes were measured using different approaches in different time periods. These data inconsistencies can influence the values of parameters including A, B, R 0 , and C eq . Despite these uncertainties, mean concentrations of all 16 solutes consistently decrease with increasing discharge, underscoring the first-order control of mean discharge. We therefore expect the general conclusion will remain the same when more data become available.

Distinct Controls of Long-Term C m Q m Versus Short-Term CQ Patterns
The long-term, cross-site C m Q m patterns here differ from the short-term, within-site instantaneous CQ patterns (C = aQ b ) at event to monthly time scales. The latter reflects instantaneous river chemistry to short term hydrological changes whereas the former carries signature of baseline response to climate conditions. The instantaneous CQ patterns have been studied extensively whereas the cross-site C m Q m patterns have received scarce attention (Godsey et al., 2019;Maher & Chamberlain, 2014;White & Blum, 1995). The instantaneous concentrations of geogenic solutes have shown dilution patterns with concentrations decreasing with increasing discharge (Ibarra et al., 2016;Torres & Baronas, 2021;. Instantaneous concentrations of DOC however have seen increase with discharge in more than 80% of the US sites (Zarnetske et al., 2018), the opposite of the C m Q m pattern here. The instantaneous CQ patterns at individual sites reflect switching flow paths under different hydrological conditions and depth distribution of source water chemistry [Botter et al., 2020;Dwivedi et al., 2018;Ebeling et al., 2021;Knapp et al., 2020;. In other words, instantaneous CQ patterns arise from mechanisms that differ from long-term C m Q m patterns. Even if these patterns at different time scales are "apparently" similar, they should not be interpreted as originated from the same mechanisms.
The C m Q m patterns of all solutes exhibit dilution patterns, albeit with differences in the forms for solutes with and without thermodynamic limitation and parameter values. This is remarkable given the diverse origin and reactions across solutes, and the wide range of site idiosyncrasies in lithology, soil type, vegetation cover, among other characteristics. This indicates that climate is the most predominant control on baseline concentrations. In contrast, the short-term CQ relationships are mostly governed by switching flow paths and distinct source water chemistry at different depths, characteristics that are more regulated by subsurface physical and biogeochemical properties, or the subsurface structure of the Critical Zone (Sullivan et al., 2022). For example, the extent of concentration variation with instantaneous discharge, quantified as the power law slope b, has been shown to depend on the magnitude of chemical contrasts in shallow soil water and in deeper groundwater, which is in turn regulated by the abundance of materials such as organic matter and weathering materials at different depths . Solutes that are more abundant in shallow soil water tends to exhibit flushing CQ patterns; solutes that are more abundant in deeper groundwater tends to exhibit dilution CQ patterns (Botter et al., 2020;Seibert et al., 2009;Stewart et al., 2022;Zimmer & McGlynn, 2018). Under conditions where only one major flow path exist, or where the shallow and deep subsurface has similar water chemistry, CQ patterns tend to be chemostatic (Zhi et al., 2019).

The Future of River Chemistry and Water Quality
How does river chemistry change under a changing climate? Climate disasters, including flooding, droughts, and wildfire, can modify material fluxes and water quality (Robinne et al., 2021;Whitehead et al., 2009). Warmer climates are expected to change river discharge, as well as the timing and magnitude of precipitation events, although the extent and direction will vary by region. River discharge is also expected to decrease due to increasing water demand by a growing population (AghaKouchak et al., 2021). Intermittent streams have already seen longer dry periods, especially in the western US (Zipper et al., 2021). Groundwater aquifers are over extracted and are subject to surface contaminations (Hartmann et al., 2021;Jasechko & Perrone, 2021;Jasechko et al., 2017;Kumar et al., 2020).
A warmer climate can possibly accelerate soil biogeochemical reactions and chemical weathering, although their response to warming depends on a variety of complex and coupled ecosystem processes that generally remain poorly understood (Greaver et al., 2016;Perdrial et al., 2018;Stegen et al., 2011). For example, vegetation structure will shift in response to shifts in precipitation amount and timing, with ensuing effects on water demand and eventually discharge (Keller, 2019). This in turn will influence river chemistry.
Nonetheless, many places have already seen elevating solute concentrations, although it is often challenging to tease apart the convoluted influence of climate versus human drivers (Guo et al., 2019;Kaushal et al., 2013;Lintern et al., 2018;Noacco et al., 2017). Salinity and electrical conductance, a collective measure of geogenic solutes (cations and anions), have been observed to rise in many places across the United States (Kaushal et al., 2018). Riverine DOC has observed increasing trend since 1980s in Europe and North America (Monteith et al., 2007), often attributed to factors including changing climate, recovery from acid rain, and modified land uses (Adler et al., 2021). In places without acid rain and human activities, DOC exhibits elevated concentrations in warmer and drier years, indicating strong climatic influence . Dissolved inorganic carbon (DIC) and alkalinity has also seen continued increase (Drake et al., 2020;Kaushal et al., 2013;Najjar et al., 2020;Raymond & Cole, 2003;Zamanian et al., 2018). Low discharge also entails longer residence time (Benettin et al., 2020) and possibly more extensive transformation of nutrients and carbon into greenhouse gases (Campeau et al., 2019).
Using the space-for-time substitution, we can interpret the data and model here for the future climate. Because all sites are minimally impacted by human perturbations, the results particularly speaks to the effects of a changing climate without the convoluted effects of human perturbations. In places that will become wetter, higher discharge can lower solute concentrations but possibly increased fluxes. This in fact has been predicted for nutrient contamination (Sinha et al., 2017). In places that will become drier, mean concentrations are expected to escalate even without human perturbation. The magnitude of concentration change will vary among solutes and hinges upon the sensitivity of concentrations to discharge variations. For example, some solutes (e.g., TOC, DOC, Cl, SO 4 ) are highly sensitive (B ∼ −1) whereas other solutes (e.g., Si) will vary less (B ∼ −0.2). With B values close to −1, mean concentrations of TOC and DOC can double when mean discharge decreases by half. This is alarming as carbon solutes are highly influential in regulating biological activities and sustaining aquatic life in rivers and streams.
The Clean Water Act, the water quality regulation in the US, establishes total maximum daily loads and MCLs in Drinking Water Standards for many solutes discussed here, including cations, nitrate, and carbon solutes (often categorized as chemical or biological oxygen demand). Increasing solute concentrations can elevate costs and energy use for water treatment, and demand renovated or augmented infrastructure for growing treatment demands (Sanders & Webber, 2012). Water quality risks therefore are an important aspect of climate risks and should be accounted for when we consider future adaptation and mitigation strategies.