Depth of Solute Generation Is a Dominant Control on Concentration‐Discharge Relations

Solutes in rivers often come from multiple sources, notably precipitation (above) and generation from the subsurface (below). The question of which source is more influential in shaping the dynamics of solute concentration cannot be easily addressed due to the general lack of input data. An analysis of solute concentrations and their dependence on discharge across 585 catchments in nine countries leads us to hypothesize that both the timing and the vertical distribution of the solute generation are important drivers of solute export dynamics at the catchment scale. We test this hypothesis running synthetic experiments with a tracer‐aided distributed hydrological model. The results reveal that the depth of solute generation is the most important control of the concentration‐discharge (C‐Q) relation for a number of solutes. Such relation shows that C‐Q patterns of solute export vary from dilution (Ca2+, Mg2+, K+, Na+, and Cl−) to weakly enriching (dissolved organic carbon). The timing of the input imposes a signature on temporal dynamics, most evident for nutrients, and adds uncertainty in the exponent of the C‐Q relation.


Introduction
Understanding how concentrations of different elements change in rivers is an open question for catchment hydrology, which has direct implications on water resources protection and management. In this respect, long-term water quality monitoring programs have been essential in collecting data of stream discharge and solute concentrations of major ions (Ca 2+ , Mg 2+ , Na + , Cl − , and K + ), nutrients (nitrate NO 3 -N and dissolved reactive phosphorus or DRP-P, called for brevity only NO 3 and DRP in the following) and dissolved organic carbon, DOC-C, called DOC. In the last decade, several studies investigated solute export in different catchments to identify commonalities based on concentration magnitude and temporal dynamics (e.g., Basu et al., 2010;Diamond & Cohen, 2018;Dupas et al., 2019;Kirchner & Neal, 2013). Most of these studies focused on regional contexts where climatic and morphologic variability is limited. The variability of solute concentrations in space and time has been studied almost exclusively in relation to the variability of the hydrological fluxes (e.g., Godsey et al., 2019;Herndon et al., 2015;Moatar et al., 2017;Wen et al., 2020) and hardly ever as an independent variable, which means investigating the order of magnitude of the concentration of each solute and its variability across different rivers.
Solute concentration magnitude is typically linked to discharge (Herndon et al., 2015) through a widely applied relation, which is the concentration-discharge (C-Q) relation. This is a fundamental descriptor of the spatiotemporal interaction between hydrological, geochemical, and biogeochemical processes in the catchment (Godsey et al., 2009;Herndon et al., 2015;Jawitz & Mitchell, 2011). The empirical C-Q relation C ¼ aQ b approximates a linear function in the log(C)-log(Q) space with the value of the b-exponent, that is, the slope of the linear regression of the C-Q relation in log-log space. This b-exponent provides insights on solute behavior in a given catchment. When the solute concentration is largely independent from discharge magnitude, that is, b ≅ 0, the behavior is called chemostatic. If the concentration increases with increasing discharge (i.e., b > 0), the solute shows an enriching behavior. If the concentration decreases with increasing discharge (i.e., b < 0), it shows a diluting behavior (Godsey et al., 2009). Recent studies have explored possible drivers of the C-Q relations, including hydroclimatic forcing (Aguilera & Melack, 2018;Blaen et al., 2017;Bouchez et al., 2017;Duncan et al., 2017;Godsey et al., 2019;Ibarra et al., 2016;Musolff et al., 2017) and catchment characteristics such as land-cover, catchment size, and lithology (Bouchez et al., 2017;Diamond & Cohen, 2018;Fields & Dethier, 2019;Ford et al., 2017;Hoagland et al., 2017;Minaudo et al., 2019). The tendency is toward identifying generalizations on the solute behavior in the C-Q space with some emerging pattern, such as the prevalent diluting behavior of geogenic solutes Godsey et al., 2009;Kim et al., 2017;Musolff et al., 2017;Thompson et al., 2011). However, observed behaviors of several other solutes, such as NO 3 , DRP, and DOC, are not easily amenable to generalizations due to contrasting observations across different catchments (Aguilera & Melack, 2018;Botter et al., 2019;Butturini et al., 2008), even though some attempts of generalization exist (Creed et al., 2015). These studies are often hampered by the limited sample of analyzed catchments and by the limited heterogeneity in climate.
Solutes observed in rivers originate from various sources including atmospheric deposition, bedrock weathering, and anthropogenic addition. For example, nutrients or other anthropogenic generated solutes can originate from either point sources (e.g., wastewater treatment plants) or nonpoint sources (e.g., agricultural fertilization). Additionally, once in the catchment, solutes also undergo different transformations such as removal through biogeochemical processes (e.g., denitrification) or sorption to soil particles. These processes, together with the diversity of possible input sources, enhance the spatial heterogeneity of solutes (Basu et al., 2010b;Musolff et al., 2015).
While annual estimates of fertilizer input in a catchment or a region is sometimes available, the spatially explicit and time-varying description of solute sources is typically insufficient for modeling applications (Costa et al., 2017;Ilampooranan et al., 2019). Consequently, recent studies investigated the spatiotemporal variabilities of solute inputs and their impact on surface water quality (Böhlke et al., 2002;Chanat & Yang, 2018;Garnier et al., 2018;Kennedy et al., 2018;McCrackin et al., 2017;McInerney et al., 2018;Zimmer et al., 2019). In some cases, inverse modeling approaches for solute input estimation are suggested (Goyette et al., 2019;Luscz et al., 2017;Miller et al., 2017).
Within the catchment, the physical and chemical contrasts between the shallow soil layer and the deep groundwater layer can impart geochemical signature on solute export patterns. Physically, shallow soils at the top 1-2 m of the subsurface often have much higher porosity and orders of magnitude higher permeability compared to the deeper, less weathered bedrocks (Heidari et al., 2017;Welch & Allen, 2014). This entails that infiltrated water often primarily flows through the very thin veneer of shallow soils with water ages mostly less than a few months old (Jasechko et al., 2016). In contrast, aquifers store much older water with ages from tens of years to millennia (Jasechko et al., 2017). Chemically, shallow soils are weathered materials that are often depleted in geogenic solutes such as Ca 2+ and Mg 2+ (Brantley et al., 2013;Jin et al., 2010); they are, however, enriched with organic materials from vegetation litter and microbial turnover (Billings et al., 2018;Lehmann & Kleber, 2015;Li, 2019). In contrast, deeper subsurface often harbor unweathered rocks that can release geogenic solutes but with low content of organic materials (Jobbágy & Jackson, 2000). Different experimental studies highlighted the link between the vertical distribution of the solute in the soil layers and the solute export at the catchment outlet (Bishop et al., 2004;Böhlke et al., 2002;Jobbágy & Jackson, 2001;Tardy et al., 2004), but only few recent studies attempted to explain it mechanistically through models (Musolff et al., 2017;Seibert et al., 2009;Wen et al., 2020;Zhi et al., 2019) linking hydrological and biogeochemical processes.
Here we address some of the above-mentioned limitations comparing solute concentrations at the outlet of 585 catchments around the world obtained as a combination of multiple data sets. We compute the coefficient of variation of the concentrations across stations in space (CV space ) and for each station in time (CV time ) as metrics to investigate the spatiotemporal variability of the concentration of each solute in absolute terms and in comparison with the others. We also compute the solute behavior, expressed through the b-exponents of the C-Q relation and compare these exponents across catchments and among solutes. Based on clues provided by the data analysis, we hypothesize that the solute input dynamics play a crucial role for the solute behavior and we test this hypothesis by means of synthetic experiments using a tracer-aided fully distributed hydrological model (Remondi et al., 2018(Remondi et al., , 2019. Specifically, we test two hypotheses: (1) the temporal dynamics of solute production determine solute behavior: continuous solute formation or injection lead to chemostatic behavior whereas sporadic solute formation or injections result in diluting or enriching behavior and (2) the depth of solute generation rather than the temporal dynamic is more influential in shaping solute behavior. The use of modeling tools and synthetic experiments to analyze the role of the 10.1029/2019WR026695

Water Resources Research
BOTTER ET AL. solute input on solute dynamics is critical, because of the consistent lack of solute input data. In summary, the objective is to unveil the role of the solute input timing and vertical distribution on the temporal variability of the solute exported from the catchment and its relation with discharge (i.e., C-Q relation) by combining extensive water quality data sets and numerical modeling.

Observations
We assembled a quasi-global water quality data set retrieving publicly available and published data worldwide. We looked for a combination of in-stream concentration data of major ions (Ca 2+ , Mg 2+ , K + , Na + , and Cl − ), nutrients (NO 3 and DRP), and DOC with discharge measurements, collected at the same gauge station. The data are from 585 stations in nine different countries (Canada, United States, Brazil, Spain, Switzerland, Germany, Sweden, Cameroon, and Australia), the detailed information of which is reported in Table S1 in the supporting information, while further information about the data sources are available in Text S1. The stations are reported in Figure 1, classified by the area-weighted mean annual discharge into hydrologically dry catchments (i.e., Q ≤ 200 mm/year), intermediate catchments (i.e., 200 < Q ≤ 500 mm/year) and wet catchments (i.e., Q > 500 mm/year).
Most of the water quality databases included discharge data colocated and synchronized with the concentration data. For Sweden and the region of Queensland in Australia, instantaneous concentration and hourly discharge data recorded in two different databases were combined.
To guarantee robustness to the results, we selected stations with at least 20 couples of concentrationdischarge observations for each analyzed solute. This is an arbitrary and pragmatic choice in order to maintain in the database as much station-solute combinations as possible but also preserving a certain representativeness of variability. Such a choice is consistent with other studies on C-Q relations (e.g., Gwenzi et al., 2017). We use data regardless of the frequency of acquisition, sampling continuity in time, or length of the monitoring period, maximizing the available information (Table S1). Although we are aware that some C-Q relations change in time Dupas et al., 2016), these changes are expected to Figure 1. Location of the water quality and discharge measurement stations representing the basin outlet. Data from 585 stations in Canada, United States, Spain, Switzerland, Germany, Sweden, Cameroon, and Australia (Queensland and Victoria) are shown. The stations are classified based on the mean annual discharge normalized by the catchment area, whenever available (549 stations). Light blue, blue, and dark blue stations are characterized by mean annual discharge <200, between 200-500, and >500 mm/year, respectively. White dots refer to stations that could not be classified due to the lack of either catchment area or mean annual discharge data. exert a minor role in comparison to across stations and across solutes variability. The assumption of temporal stationarity is consistent with the purposes of the study, which aims at characterizing general patterns in the long-term solute magnitude, variability, and C-Q relations. Consequently, for each station, the number of observations differs across solutes.

Data Analysis
We explore the magnitude of the observed in-stream concentrations for seven solutes including calcium (Ca 2+ ), magnesium (Mg 2+ ), potassium (K + ), sodium (Na + ), chloride (Cl − ), nitrate (NO 3 ), DOC, and DRP by means of the basic statistics (i.e., median, 75 th and 25 th percentiles). These solutes have different dominant sourcing processes, so that it is worth adopting definitions (albeit simplified) to cluster them in categories. We will refer to solutes mainly originating from weathering (Ca 2+ , Mg 2+ , and K + ) as weathering products and to those solutes mostly linked to anthropogenic sources, such as road salt spreading and fertilizers (Rose et al., 2018), as exogenous solutes (Na + , Cl − , NO 3 , and DRP). DOC is treated separately and not included in the two categories.
The magnitude analysis allows comparing catchments, which are clustered in three classes following their area-weighted mean annual discharge. This clearly generates significantly different discharge means across classes (two sample t test; Cressie & Whitford, 1986). We summarize the variability of the magnitude of the concentrations in space and in time with the indexes CV time and CV space , respectively, defined as follows: where CV all stations is the vector of temporal CVs of concentration for all stations, while C all stations is the mean concentration of each solute for each catchment. CV time represents the central estimate of the temporal variability of concentration in individual catchments, whereas CV space reflects cross-catchment concentration variations. Given the lack of specific knowledge of anthropic activities occurring in a given catchment, we also subdivide the catchments into disturbed and undisturbed catchments, based only on the mean concentration of NO 3 . The mean value of NO 3 should be correlated with nutrients introduced in catchments as inorganic fertilizers, even though some denitrification occur during subsurface retention and transport and it can be different across catchments (Basu et al., 2010;Musolff et al., 2015). We considered catchments with a mean nitrate concentration above 1 mg nitrate/L (i.e., 0.2259 mg NO 3 -N/L) as disturbed, since 1 mg nitrate/L is a background concentration value rarely exceeded in natural conditions Zobrist, 2010).
Subsequently, we compute the C-Q linear regression in a log(C)-log(Q) scale and classify the solute behaviors on the basis of the b-exponent. Given the large sample of analyzed stations, we are interested in assessing the magnitude and variability of the b-exponent and comparing it across different solutes and across catchments aggregating them in major classes corresponding to wetness conditions. We also take into consideration two alternative classifications, based on catchment area (three groups: area ≤1,000, between 1,000 and 10,000, and >10,000 km 2 ) and nitrate concentration, as a proxy of anthropic influence representing mostly undisturbed, impacted, and highly impacted catchments (three groups: NO 3 ≤ 0.23 mg/L, 0.23 < NO 3 ≤ 2.3 mg/L, and NO 3 > 2.3 mg/L). The combination of the impacted and highly impacted catchments is summarized in the class of disturbed catchments in terms of CV analysis as defined above. Differently from the previous analysis, which includes all data, we preprocessed data deleting outliers, defined as those with values that are more than three standard deviations away from the median. This data filtering is needed since linear regression is very sensitive to outliers. We also apply a Student's t test on the statistical significance of having the b-exponent different from zero setting the significance threshold α at 0.05. We compare our results with results obtained with other definitions of chemostatic behavior, that is defined imposing a threshold on the absolute value of the b-exponent below which the behavior is defined chemostatic. We assume thresholds of 0.1 and 0.2 as these thresholds have been often used in literature (

Numerical Experiments 2.3.1. The WATET Model
We designed numerical experiments to test if the timing and the vertical distribution of the solute input are important determinants of the solute behavior in the C-Q space. We used the WATET (Water Age and Tracer Efficient Tracking) model, which couples a spatially fully distributed hydrological model and a module that simulates conservative solute transport. The model was introduced with a full description by Remondi et al. (2018). Briefly, the hydrological component of WATET runs at 5-min time steps and keeps track of the water mass balance in each cell of the gridded domain accounting for main hydrological processes such as precipitation, evapotranspiration, infiltration, recharge to the aquifer, lateral flow from soil and groundwater, overland, and channel flow. Lateral transfer of soil water, overland, and channel flow are simulated using the kinematic wave approximation (Ciarapica & Todini, 2002), so that the stream water can come from surface runoff, lateral flow from the shallow soil profile, and groundwater flow from an aquifer storage ( Figure 2). The groundwater flow is represented with a nonlinear reservoir (Benettin et al., 2015) for each cell, meaning that the drainage of water from one groundwater storage cell to the adjacent cell follows a power function that is assumed to be the same for each cell of the catchment.
The module of transport tracks the solute concentrations in each storage compartment, that is, soil, aquifer, surface, and channel (for each cell) of the domain assuming that the solute is advected with the water and that water is well mixed in each cell, therefore neglecting dispersion and considering the solute as conservative. In thise numerical experiments, the solute input is injected through precipitation or directly in the soil or groundwater storage and is distributed uniformly in space.
For the numerical experiments, we used meteorological forcing for the year 2008 and land-use, soil, aquifer, and topographic features of the Hafren catchment (3.6 km 2 ) in middle Wales (UK), which is part of the Plynlimon experimental catchment managed by the Center for Ecology and Hydrology (CEH). The catchment has mildly heterogeneous soil type and vegetation cover (Supplementary Information in Remondi et al., 2018). The domain was discretized using a digital terrain model with 24.4 m × 24.4 m resolution; each cell is characterized by the own land-cover and soil properties. Hourly precipitation input in each cell was the result of the spatial interpolation of the measurements at two stations located in the proximity of the catchment (Carreg Wen and Tanllwyth). Potential evapotranspiration was also provided to each cell as input and was computed as function of vegetation cover and meteorological forcing (Remondi et al., 2018). Since 1980s, a considerable number of studies have been investigating the Plynlimon catchment (Brandt et al., 2004;Kirby et al., 1991Kirby et al., , 1997, with particular attention devoted to stream water chemistry (Durand et al., 1994;Kirchner et al., 2010;Knapp et al., 2019;Neal et al., 2003;Neal & Kirchner, 2000). This catchment was also used to develop and test the WATET model. Remondi et al. (2018) report details of meteorological inputs, parametrization, and validation results.

Numerical Experiments-Vertical Distribution of Solute Generation
To evaluate the influence of the timing and vertical distribution of the solute input on water chemistry, we kept the model setup unchanged and we forced the model with different solute input time series, given at different depths and with different frequencies. Since WATET does not solve explicitly hydrological Figure 2. Conceptualization of the lateral fluxes in the numerical model. In each cell of the simulated domain, the water and solute fluxes can move laterally from each compartment (i.e., surface, soil, and groundwater storage) to one or more neighboring hillslope or channel cells; water can also move vertically among compartments.

10.1029/2019WR026695
Water Resources Research dynamics in multiple soil layers, the vertical differentiation is obtained by injecting input in different model compartments (soil and aquifer) and using soil moisture thresholds for the activation of transport as a proxy of vertical processes. Under wetter conditions, the shallow paths in the soils are more likely to be activated. A similar consideration could be made for spatial activation of distant paths within a cell. Specifically, we provide a fixed quantity of solute input conserving the total mass of solute over the year, corresponding to the observed mass of chloride in precipitation measured at Carreg Wen meteorological station in the Plynlimon catchment over 2008.
A conceptual scheme of the numerical experiments is illustrated in Figure 3. The vertical distribution of solute input is produced injecting the input in different compartments ( Figure 3a): (i) on the surface through deposition, (ii) 100% in the soil compartment, (iii) 70% of the input in the soil compartment and 30% in the aquifer, and (iv) 70% in the aquifer and 30% in the soil. This group of experiments aims at representing a range of sourcing of solutes from deposition in the surface, to nutrient leaching due to soil organic matter turnover in shallow soils up to weathering of primary minerals that can happen at variable depths in the soil profile. We also include the case where the solute is mobilized (both on the horizontal and vertical directions) only when a soil moisture threshold (the same for vertical and horizontal) is exceeded (Figure 3b). Below this threshold, the hydrological flux leaving the soil cell carries a null concentration of solute in both directions. Such a case aims at representing a mechanism similar to transport of a solute present in the topsoil layer only, meaning that the horizontal and vertical solute transport capacity are impeded until a fixed threshold of soil moisture is reached. Given the saturated water content of 0.61 used in the simulations (Remondi et al., 2018) and aiming at the activation of solute export only in particularly wet conditions, we fix soil moisture thresholds of 0.50, 0.56, and 0.60 to mimic such a mechanism.
Additionally, we ran a series of experiments where solute is added by deposition but we included solute degradation dynamics in the groundwater (Figure 3c), expressed with a first order kinetic. The concentration change was computed with the degradation law where C ini is the concentration at the beginning of the time step, t is time, and K deg is the degradation constant that defines the rate of solute removal. We ran the experiments setting the K deg values at 10 −3 and 10 −4 h −1 . This simple kinetic is widely used to represent processes such as denitrification or solute degradation (Costa et al., 2017;Husic et al., 2019;van der Velde et al., 2010). This experiment is representative of the dynamics of those solutes such as nitrate, which is deposited on the surface (i.e., fertilization) and can undergo degradation once they leach into the groundwater compartment. In these experiments where the impacts of the soil moisture activation thresholds and, of the solute degradation, was tested, the solute was injected only through deposition.

Numerical Experiments-Timing of Solute Input
To test the impact of the timing of the input, we run for each set of experiments 10 simulations, each one characterized by a different number of input applications, spanning from a single application, where all Simulations were run injecting all the solute in the surface through deposition (I), generating all the solute in the soil layer (II), generating 70% of the total solute in the soil layer and 30% in the groundwater storage (III), and generating 70% of the total solute input in the groundwater storage and 30% in the soil layer (IV). (b) Soil moisture activation threshold. The solute was injected on the surface through deposition; solute export from the soil (in both vertical and horizontal directions) was activated only above a fixed threshold of soil moisture: 0.50 (I), 0.56 (II), and 0.60 (III). (c) Solute degradation in the groundwater. The solute was injected through deposition. Once in the groundwater, it undergoes degradation accordingly to a first order kinetic with different degradation constants: The total mass of the injected/geo-generated input is preserved across all experiments. For each experiment, the impact of the timing of the input was tested by running 10 simulations with different solute input frequencies. 10.1029/2019WR026695

Water Resources Research
the input amount is given at once, to a simulation that gives inputs of solute every day ( Figure 4). While certain combinations of temporal frequency and experiment might not be realistic, for instance, sporadic injections applied mostly in the groundwater, we simulate all the input frequencies for all the experiment sets for sake of completeness. The goal is to compare the variance induced by the different numbers of application in the different cases. In case of input entering the catchment per deposition, the maximum number of injections corresponds to 227, the number of rainy days (i.e., days with at least 1 mm of precipitation). In case of input generation in the soil or aquifer, the maximum number of injections corresponds to 366 (days in the year 2008). Except for the cases where the solute is applied in any (rainy) day, the choice of the days of application is arbitrary, therefore, we replicate each experiment 10 times consisting of multiple stochastic simulations, selecting the days of the applications randomly ( Figure 4). In this way, we address the uncertainty associated with the period of application. We removed the impact of the initial conditions running each simulation characterized by a different input for a spin-up period of 20 years.

Surface-Groundwater Concentration Ratio
We summarize the results of different experimental setups in terms of b-exponent and C sw /C gw ratio, a ratio that has been found to be crucial in determining C-Q relationships (Zhi et al., 2019). Here C sw is the mean solute concentration in the soil water compartment, and C gw is the mean solute concentration of the groundwater compartment. Both variables are averaged in space (i.e., across the domain) and in time (i.e., across the simulation period corresponding to the year 2008). The C sw /C gw ratio represents a measure of the homogeneity of the vertical distribution of the solute. That is, a C sw /C gw ratio closer to 1 means a more homogeneous vertical distribution of the solute concentration. Values lower than 1 mean higher concentration in the groundwater than in the soil water, while values higher than 1 imply higher concentrations in the shallow soil layers than in the groundwater.

Space-Time Variability of Observed Solute Concentrations
The boxplots in Figure 5a synthesize the observed mean concentration of each solute across the catchments grouped per mean annual discharge in mm year −1 . The concentrations of different solutes span different orders of magnitudes, with the exogenous Na + and Cl − and the geogenic Ca 2+ exhibiting particularly high concentrations and variability, with concentrations almost one order of magnitude higher than the geogenic solutes Mg 2+ and K + . DOC concentration is on average 1 order of magnitude higher than NO 3 and 2 orders of magnitude higher than DRP. The classification in classes of discharge highlights the expected tendency to obtain lower concentrations, that is, higher dilution, with higher average specific discharge for all solutes except for the nutrients NO 3 , DRP, and partially K + that exhibit higher concentrations in intermediate wetness catchments rather than dry or wet catchments. The t test confirms that data in class 2 of NO 3 have a significantly different mean than in the other classes (Table S2a). When classified by catchment size, mean concentrations are typically increasing with catchment area except for NO 3 and DRP and partially DOC ( Figure S1). However, catchment area represents a quite weak control on mean concentrations, as supported by the few statistically different means ( Table S2b). The classification of catchments based on the mean concentration of NO 3 ( Figure S2), one of the main components of inorganic fertilizers, shows higher values of concentrations for a number of solutes in the catchments with higher NO 3 . Except for DOC, the value NO 3 is a quite strong control in separating average solute concentrations ( Table S2c).
The geogenic solutes mostly show diluting behavior, that is, negative b-exponent of the C-Q relation, with the 25 th and 75 th percentiles of b included between −0.37 and −0.008 (Figure 5b). Na + and Cl − also exhibit diluting behavior, but spanning a wider range of negative b values, being −0.24 the median and reaching the lowest value of −0.93. The nutrients NO 3 and DRP span both positive and negative b values, while DOC shows a slightly positive b with limited variability across catchments. Generally, the exogenous solutes exhibit larger b spreads than geogenic solutes and DOC and show a slight tendency of decreasing b-exponent with increasing discharge. The behaviors of Ca 2+ , Mg 2+ , and K + are consistently diluting regardless of the different hydrological regimes. On the basis of these results, we summarize the behaviors of solutes in four main groups. The geogenic solutes (Ca 2+ , Mg 2+ , and K + ) show a consistent diluting behavior, solutes associated with salt (Na + and Cl − ) show a pronounced diluting behavior but with strong differences across catchments, meaning that the b-exponent assumes negative and highly variable values. The C-Q patterns of nutrients NO 3 and DRP vary from strong dilution to strong enrichment, whereas DOC shows smaller variability with a weak enriching behavior.
A conceptual scheme is illustrated in Figure 6a to guide the interpretation of the CV time /CV space and CV time / b-exponent relations (see also Musolff et al., 2015). The results show that for both undisturbed (left) and disturbed catchments (right), almost all points fall above the 1:1 line, indicating higher CV space than CV time (Figure 6b). This is the case although the difference is more pronounced in disturbed catchments where the spatial heterogeneity is enhanced. The higher CV space of disturbed catchment is not justified by higher discharge variability ( Figure S3).
In disturbed catchments, as expected, exogenous solutes exhibit higher CV time compared to Ca 2+ , Mg 2+ , K + , and DOC. In undisturbed catchments, nutrients (NO 3 and DRP) and DOC still have higher CV time than geogenic solutes and salt solutes.

Water Resources Research
In the b-exponent versus CV time figure, two clusters emerge especially in the undisturbed catchments ( Figure 6c). One group comprises the geogenic solutes (Ca 2+ , Mg 2+ , and K + ) and the salt solutes (Na + , and Cl − ) with relatively low temporal variability (CV time ) and pronounced dilution behavior (negative b) remarkably similar across solutes in undisturbed catchments. The other group includes nutrients (NO 3 and DRP) and DOC with larger and highly variable CV time and nearly zero or slightly positive b-exponents. For disturbed catchments, the difference between these two clusters is less pronounced because Na + and Cl − b-exponents are lower, temporal variability is larger, while CV time of DOC and NO 3 is reduced.

Chemostatic Behavior: How Often?
Some C-Q relation studies consider b-exponents non significantly different from zero when their absolute value is lower than a fixed threshold, for example, equal to 0.1 (Herndon et al., 2015) or even to 0.2 (Zhi et al., 2019). In Table S3, we compared the percentage of b-exponents computed from the observations significantly different from zero applying both a statistical test (the Student's t test based on an α ¼ 0.05 significance level) or the arbitrary (but commonly used) thresholds of 0.1 and 0.2. Applying a threshold for b of 0.1 and 0.2 generates, respectively, on average 42.6% and 67.1% of cases for all solutes indistinguishable from 0 (chemostatic behavior), while applying the statistical test shows that only 4.5% of the station/solute combinations have indeed a chemostatic behavior. This indicates that most of solute/stations are not geochemical  (c), the left plots refer to undisturbed catchments (defined as catchments with mean nitrate concentration lower than 1 mg/L) and the right plots refer to disturbed catchments (defined as catchments with mean nitrate concentration higher than 1 mg/L). The uncertainty bars span between the 25th and the 75th percentile.

10.1029/2019WR026695
Water Resources Research stationary and the use of fixed thresholds might overemphasize this behavior. This difference is especially pronounced for the geogenic species, for which data from observations fail the chemostatic test and rather show a consistent diluting behavior.

Numerical Experiments
The Results of the numerical experiments with different injection frequencies at randomly selected days of the year represent the effect of timing (and seasonality) of the input, while cases with different vertical distribution of the input, soil moisture activation thresholds, and degradation dynamics represent impacts of vertical distribution (Figure 7).
The three sets of experiments exhibit clearly different responses. The generation of solute at different depths leads to different vertical distribution of solutes in the model compartments, which we summarize by means of the C sw /C gw ratio. Injections per deposition or solute generation in the soil do not lead to a strong vertical stratification (C sw /C gw ≅ 1), whereas deeper solute generation leads to low C sw /C gw ratio, which tends toward zero in the case where 70% of the solute is generated in the groundwater (Figure 7, red axis).
Values of the C sw /C gw ratio higher than 1 are obtained when the solute is injected through deposition, but it is mobilized only in particularly wet conditions (Active threshold set of experiments in Figure 7, red axis). The higher the threshold on the soil moisture for solute mobilization, the higher the resulting C sw /C gw ratio. The case with residence time-dependent depletion of the solute in the groundwater also exhibits C sw /C gw ratio large than 1, higher with increasing degradation rate. The timing of the input adds variability to the observed C sw /C gw ratios, especially evident in the "Activation Threshold" experiments. For the red axis, the dashed horizontal red line is set at 1 and represents the reference homogeneous vertical distribution of the solute. The scale of the y axes is logarithmic. For the blue axis, the dashed horizontal blue line is set at 0, which represents the biogeochemical stationarity line, below which the behavior leans toward dilution and above which the behavior leans toward enrichment. In both cases, each boxplot refers to a different experimental setup and includes all simulations with different frequencies of solute injection at randomly selected days. The experiments "Deposition," "Soil," "70% soil 30% GW," and "30% soil 70% GW" (white boxes) refer to the different compartments where the solute input is injected. The "Activation threshold" experiments (gray boxes) consist in setting thresholds on soil moisture below which the solute cannot be mobilized. The selected thresholds are 0.50, 0.56, and 0.60. The "Degradation" experiments (black boxes) are characterized by deposition input and by the nonconservative behavior of the solute, that is, the solute undergoes degradation while residing in the groundwater.

Water Resources Research
The b-exponent decreases with increasing depth of solute generation, spanning between biogeochemical stationarity, in case of input through deposition, to strong dilution (i.e., b-exponent ¼ −0.5) when the solute is mainly generated in the groundwater compartment. The timing of the input adds uncertainty to the definition of solute behavior, resulting in higher variance of the b-exponent. This variance decreases with increasing depth because of the slower aquifer response. Imposing a soil moisture activation threshold for solute export generates an enriching behavior, which becomes stronger with the threshold approaching saturation, that is, with solute that is mobilized only during the wettest conditions. For this set of experiments, the variance due to the timing of the input application is higher than for the other two cases.
Although a first order degradation kinetic is a simplification of the real processes that are removing solutes from the catchment, the numerical experiments show that reactive properties can play a role on the solute behavior. The results show that degradation changes behavior from chemostatic to enrichment with increasing degradation rates. Given the same residence times dictated by the hydrological processes, which are identical in all simulations, a larger degradation constant is removing more solute from the groundwater resulting in low concentration during low-flows and thus positive b-exponents.
Combining the C sw /C gw ratio with the b-exponent, a clear pattern emerges synthesizing the different numerical experiments (Figure 8). When the solute is mainly generated deep in the groundwater, that is, the C sw /C gw ratio is close to zero due to the high concentration in the groundwater, resulting in strong diluting behavior with b-exponent values up to −0.5 (an extreme case based on observed b-exponents). When the solute is injected through deposition, if the solute is mobilized independently from the wetness condition and it does no undergo depletion in the groundwater, there is no evident vertical stratification in the catchments (i.e., C sw /C gw ≅ 1) and the behavior is nearly chemostatic (i.e., b ≅ 0). In contrast, when the solute is injected through deposition and the concentration in the groundwater is depleted by degradation, the C sw /C gw ratio is higher than 1 and the export dynamics exhibit enriching behavior. The enrichment is even more pronounced and the C sw /C gw ratio assumes higher values (i.e., up to 4) when the solute mobilization is activated only in wet conditions, for example, when surficial paths are likely to be activated. High threshold means that solutes remain high in the shallow zones without being flushed out until very wet conditions occur, therefore maintaining high C sw /C gw ratios.

Clues From Observations on Temporal and Spatial Variability
The mean concentrations across catchments classified with area-weighted discharge are expected to decrease with increasing mean annual discharge, that is, the solute should be more diluted when the water availability is more per unit area if solute input does not change systematically with wetness. This is somehow intuitive and indeed observed for most solutes, but K + , NO 3 , and DRP do not follow this pattern, exhibiting higher concentrations under intermediate wetness conditions. This might be due to the exogenous anthropic forcing, which changes the natural cycling of these solutes and increases variability. Nutrients are also prevalently located in the shallow soil and are correlated with biological activity. A certain level of wetness is needed to generate biomass and nutrient leaching and also to connect enough nutrient sources to the stream network. However, when wetness increases further, lower concentrations may prevail as a result of nutrient depletion.
The comparison of observations in undisturbed and disturbed catchments based on the mean concentration of NO 3 ( Figure S2) indicate that the anthropic forcing increases both spatial variability and time variability of each solute, especially the exogenous ones (Figure 6b). Nutrients mostly originate from soil biogeochemical reactions in undisturbed catchment (Cleveland et al., 2013;Fatichi et al., 2019;Schlesinger & Bernhardt, 2013) and from large but sporadic fertilization events in disturbed catchments. Nutrient b-exponents are characterized by strong uncertainty spanning from strongly diluting to strongly enriching, as similarly

Water Resources Research
shown in other studies (Aguilera & Melack, 2018;Botter et al., 2019;Butturini et al., 2008;Hunsaker & Johnson, 2017). The sporadic input application is likely the driver of the high CV time and of its high variability for the exogenous solutes (Figure 6a). This differs from DOC, the sources of which are mostly related to primary productivity and associated reactions of litter and soil organic matter decomposition and therefore it is produced in a relatively continuous way with notable exceptions (as too dry or frozen soils). As a result, the b-exponents of DOC are expected to be more associated to the vertical or spatial distribution and less affected by input timing. In disturbed catchments, DOC exhibits CV time comparable to geogenic solutes (Figure 6b), indicating low temporal variability. Its b-exponents, however, are similar to nutrients (Figure 6c), as they reflect the higher organic carbon content in shallow soils, in a way similar to NO 3 and DRP (Figure 6c).
Regardless of the criteria of catchment classification, the C-Q relation reveals very consistent diluting behavior for the geogenic solutes. This agrees with some literature (e.g., Botter et al., 2019;Diamond & Cohen, 2018;Moatar et al., 2017;Thompson et al., 2011;Winnick et al., 2017;Wymore et al., 2017) but contrasts the idea of a pervasive geochemical stationarity (Godsey et al., 2009;Herndon et al., 2015;Zhi et al., 2019). This is likely related to the use of fixed thresholds on the b-exponent rather than statistical tests, as illustrated by our analysis, or generally to the challenge of classifying behavior for b-exponent values close to zero (Thompson et al., 2011). Despite the high variability of mean concentrations for Na + and Cl − , they also show a highly diluting behavior. Geogenic solutes are mostly generated through chemical weathering, which is a continuous process, while Na + and Cl − are generated through weathering, but they may enter the catchments also through deposition and through external sourcing, such as in the case of the spread of deicing salt during winter. For geogenic solutes, including Na + and Cl − , the relation between CV time and the b-exponent in Figure 6c shows a cluster of solutes in undisturbed catchments, consistently with the findings of Jawitz and Mitchell (2011), who provide an analytical representation of the relation between temporal CV and the b-exponent. This observation, together with the evidence of high CV space to CV time ratio in both undisturbed and disturbed catchments for geogenic solutes, suggests that the timing of the input is not the main driver of the geogenic solutes export dynamics. This is consistent with the fact that unweathered rocks tend to be more abundant at depth compared to shallow soils with more weathered material and less Ca 2+ and Mg 2+ . In addition, slow flow at depth prolongs water-rock contact time and elevates concentrations to higher levels. Deeper subsurfaces are also less sensitive to temporal variability in surface weather. Instead, they are more likely to be shaped by lithological difference that can span from carbonates to shales to siliclastic rocks across the 585 sites. All factors contribute to an increase CV space and to downplay the role of temporal variability versus the one of vertical variability.
The larger intercatchment variability (i.e., higher CV space to CV time ratio) of nutrients, especially in disturbed catchments, could be related to geological substrates differences affecting nutrient transport and/or homogenizing role of biological activities or fertilization (Figure 6a), even though interpretations are speculative based on concentration data alone. The smaller b-exponent suggests that the depth of solute generation is important as discussed later (Figure 6c). This supports the recent findings of Zhi et al. (2019), which highlighted the important role of contrasts between soil water and groundwater as the main driver of the C-Q relation. They show that a much higher concentration in the soil water than in groundwater translates into an enriching behavior, whereas higher concentrations in groundwater lead to a diluting behavior. A similar concept but connected to riparian areas was reported by Seibert et al. (2009) who provided an analytical formulation connecting the hydrochemical vertical profile in the riparian zone to the observed stream solute loads. Musolff et al. (2017) also combined data analysis and modeling experiments and showed that the spatial heterogeneity structure of solute sources and solute losses along flow paths is the main driver of C-Q relations.
Overall, observations suggest that both the timing of input and the vertical distribution of solute source can affect the solute behavior with the latter prominent for geogenic solutes. A broader generalization and quantification of these two key players is provided by the numerical experiments with the tracer-aided model.

Interpretation Through Numerical Experiments-b-Exponents and C sw /C gw Ratios
Results of the numerical experiments show decreasing b-exponent with increasing depth of solute generation. According to data reported in Jobbágy and Jackson (2001), Na + and Cl − concentrations are higher at

10.1029/2019WR026695
Water Resources Research depth than K + , Ca 2+ , and Mg 2+ , due to the higher solubility of Na + and Cl − , which favors their leaching.
Observations show a moderately diluting behavior of Ca 2+ , Mg 2+ , and K + and a highly diluting behavior of Na + and Cl − . A highly diluting behavior is consistent with the results of the modeling experiments in which the solutes are injected at depth (Figure 7). Na + and Cl − dilution is more evident in disturbed catchments where concentrations are higher, most likely due to the spreading of deicing salt on the surface (Sansalone & Glenn, 2002), and where their dynamics are linked to the snow melting process, which can both favor water percolation in depth and also activate shallow flow paths deprived of solutes. The b-exponents from the numerical simulations spans the range from 0.2 to −0.5, which encompass the observed b-exponents. Deeper solute generation leads to more dilution behavior (more negative b values). This is because, when solutes are mostly generated deeper in the soil and transported in groundwater (e.g., Ca 2+ , Mg 2+ , and K + ), the stream water has relatively constant supply of solutes from the aquifer but negligible contribution from shallow soil and surface water such that these sources dilute the more enriched groundwater. The C sw /C gw ratio supports this explanation, since deeper solute generation leads to lower concentration ratio.
The activation thresholds in the numerical experiments conceptually represent the conditions under which export occur from otherwise inactive layers. At high thresholds, the solute was mobilized only when soil water content approached saturation. In reality, this can be seen as the shallowest soil layers that contribute only under very wet conditions. Low to high thresholds therefore conceptually represents increasing hydrological connectivity. Pronounced enriching behavior occurred only at high threshold conditions (Figure 7). In this case, the solute accumulated in shallow soil layers until solute export is activated under close-tosaturation conditions. The accumulated solutes are exported in both horizontal and vertical directions when the soil approaches saturation. In such saturated conditions, the horizontal export capacity is more efficient than the vertical one. The joined effect of these features leads to the observed enriching behavior. This is also the case with the highest C sw /C gw ratio.
Observations suggest that enriching behavior occurs often for NO 3 , DRP, and, moderately, also for DOC. Literature supports the hypothesis that enriching behavior for NO 3 and DOC is due to the increased hydrological connectivity to zones of enriched organic carbon under very wet conditions (Herndon et al., 2015;Laudon et al., 2011;Outram et al., 2016;Zarnetske et al., 2018;Zhi et al., 2019). The results from the numerical experiments suggest that the increased hydrological connectivity is likely to happen because of the activation of surficial water paths. In the case of DRP, the explanation can be different, as DRP dynamics are related to suspended sediments, so that the triggering factor in case of wet conditions might be soil erosion and solid transport rather than increased hydrological connectivity (Thompson et al., 2011).
An enrichment behavior, although of smaller magnitude, is also detected when a degradation kinetic is introduced in the groundwater compartment. When the solute is injected through deposition, we observe biogeochemical stationarity ( Figure 5b); however, assuming a degradation kinetic in the groundwater leads to solute depletion at low-flow (with generally longer residence time), which corresponds to a C sw /C gw ratio higher than 1 (Figure 7, blue axis). This effect is increasing with increasing degradation rate (Figure 8). The same results were obtained by Musolff et al. (2017). These mechanisms are not mutually exclusive: activation of surficial paths during wet periods and degradation of solute in groundwater can both contribute to the observed positive b slopes. This can be the case of NO 3 , which often exhibits a mixed behavior. Nitrates can accumulate in the shallow soil layers in agricultural catchments and export under wet conditions; they can also percolate into the groundwater, where they undergo denitrification before reaching the channel through slow groundwater release Benettin et al., 2020;Dupas et al., 2016). The results here remark and provide numerical support to the importance of vertical distribution of solutes (Bishop et al., 2004;Koven et al., 2013) and solute transport activation in specific areas or depths shown before (Basu et al., 2010;Musolff et al., 2015;Seibert et al., 2009;Wen et al., 2020;Zhi et al., 2019).
The CV time of the synthetic experiments does not span the range of CV time from observations, even when considering different temporal aggregations of the simulation outputs ( Figure S4). This can be due to the use of the Hafren catchment for the numerical experiments that experience a wet climate throughout the year and relatively low spatial heterogeneity in soil and vegetation properties (Brandt et al., 2004;Kirby et al., 1991) compared to the variety of climatic and land-cover conditions encompassed by the observations. Regardless of this limitation, we identified some correspondences between our experiments and the clusters in the CV time /b-exponent relation. The input frequency affects CV time and adds some variability in the C-Q relation. However, its effects on b-exponent are not as significant as the vertical distribution of the input. The increase in variability of the C-Q relation without a consistent change of the b-exponent was also detected numerically by Musolff et al. (2017) modifying the velocity of the immobile-mobile mass exchange. Independently from the modeling setup, the C sw /C gw ratio emerges as the dominant control on the solute export behavior (Figure 8) in strong agreement with the recent study of Zhi et al., 2019. The fact that b values from numerical experiments cover a much larger range than data, especially for low b values, may indicate that some assumptions in the numerical experiments may not necessarily represent natural catchments, as natural systems rarely have C sw /C gw ratio lower than 0.4 (Figure 8).

Limits of Interpretation
In order to properly evaluate the results, it is important to highlight the limitations of the study and especially of the numerical experiments. The hydrological model includes a number of simplifications, such as the non-linear representation of the groundwater reservoir with spatially homogeneous parameters, the well-mixing condition within each storage compartment or the fully advective transport that neglects any form of solute dispersion. Most importantly, we compared the data obtained from 585 catchments with the results of synthetic numerical experiments in a single catchment. Modeling results are influenced by the specific climatic forcing and topographic and land-use characteristics of the Hafren catchment. This limitation emerges, for example, in the results of CV time and of the b-exponent. The CV time from the numerical experiments spans a small range of values compared to observations, because they represent only one of all possible climatic conditions and catchment characteristics embedded in the quasi-global data sets. Other model simplifications can also influence the values of b-exponents, CV time, and C sw /C gw ratios, such the lack of representation of certain processes as stream solute retention and removal, or concurrent combination of different solute sources. The simplifications can be problematic if the model is used for prediction without calibration. In this work, however, the model has a diagnostic role and it serves the purpose of unraveling the mechanisms and identify key controls of patterns arising from observations. In fact, the agreement between modeling results, the interpretation of the observed patterns, and conclusions of previous studies using more complex models or other data does reinforce the main conclusions from the numerical experiments. If the simulations were run under more diverse conditions (e.g., multiple catchments), we do expect the overall patterns to be similar to results shown here, although the exact numbers may vary considerably (e.g., Figure 8).

Summary
The observations from 585 sites generally indicate a higher degree of intercatchment spatial variability than within-catchment temporal variability. The concentrations of geogenic solutes (Ca 2+ , Mg 2+ , and Na + ) generally decrease from semiarid to humid climatic conditions and respond negatively to increasing discharge (dilution behavior). In contrast, the concentrations of anthropogenic (exogenic solutes), including NO 3 and DRP, correlate more to the level of human disturbance rather than climate conditions. They also tend to be insensitive (chemostatic) or respond positively (enriching behavior) to increasing discharge.
Numerical experiments using a simplified solute generation and transport model complement the observations by discerning the impacts of different mechanisms. Results show that the vertical distribution of solute generation and activation thresholds of solute export are more influential in shaping the solute behavior (C-Q relation) than the timing of solute input. The depth at which the solute is supplied or generated plays the major role in defining the b-exponent for highly diluting (Na + and Cl − ), weakly diluting (Ca 2+ , Mg 2+ , and K + ) and enriching solutes (e.g., DOC). Sporadic input applications increase temporal variability and add uncertainty on the solute behavior especially for NO 3 and DRP, the response of which remains the most difficult to constrain. This study offers a more generalizable view on C-Q relation than currently available.

Data Availability Statement
The database used in this study can be downloaded at the web sites reported in Text S1. The code of the model WATET is available at the link (10.6084/m9.figshare.12618854).