# Antecedent Hydrologic Conditions Reflected in Stream Lithium Isotope Ratios During Storms

Jon K. Golla, Julien Bouchez, Jennifer L. Druhan

► To cite this version:

Jon K. Golla, Julien Bouchez, Jennifer L. Druhan. Antecedent Hydrologic Conditions Reflected in Stream Lithium Isotope Ratios During Storms. Geophysical Research Letters, 2024, 51, 10.1029/2024GL109624. insu-05081565

HAL Id: insu-05081565

https://insu.hal.science/insu-05081565v1

Submitted on 24 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Geophysical Research Letters®

## RESEARCH LETTER

10.1029/2024GL109624

**Key Points:**

* Stream lithium stable isotope ratios (δ⁷Li) recorded at high frequency over storm events are sensitive to antecedent conditions

* A reactive transport model cannot produce observed shifts in stream chemistry through variations in flow rate alone

* Flux-weighting of model fluid outputs based on time-varying fluid transit time distributions describes stream δ⁷Li over storm hydrographs

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

J. K. Golla,
golla1@llnl.gov

**Citation:**

Golla, J. K., Bouchez, J., & Druhan, J. L. (2024). Antecedent hydrologic conditions reflected in stream lithium isotope ratios during storms. *Geophysical Research Letters*, *51*, e2024GL109624. https://doi.org/10.1029/2024GL109624

Received 5 APR 2024
Accepted 16 AUG 2024

**Author Contributions:**

**Conceptualization:** Jon K. Golla,
Julien Bouchez, Jennifer L. Druhan

**Data curation:** Julien Bouchez

**Formal analysis:** Jon K. Golla,
Julien Bouchez, Jennifer L. Druhan

**Funding acquisition:** Jon K. Golla,
Julien Bouchez, Jennifer L. Druhan

**Investigation:** Jon K. Golla

**Methodology:** Jon K. Golla,
Julien Bouchez, Jennifer L. Druhan

**Resources:** Julien Bouchez, Jennifer
L. Druhan

**Software:** Jon K. Golla, Jennifer
L. Druhan

**Supervision:** Julien Bouchez, Jennifer
L. Druhan

**Visualization:** Jon K. Golla

© 2024. The Author(s).
This is an open access article under the
terms of the Creative Commons
Attribution License, which permits use,
distribution and reproduction in any
medium, provided the original work is
properly cited.

<figure><img src="image_2.png" /></figure>

## Antecedent Hydrologic Conditions Reflected in Stream Lithium Isotope Ratios During Storms

Jon K. Golla<fnref n="1" />,<fnref n="2" /> ID, Julien Bouchez<fnref n="3" /> ID, and Jennifer L. Druhan<fnref n="1" />,<fnref n="3" /> ID

**Abstract** Antecedent hydrological conditions are recorded through the evolution of dissolved lithium isotope signatures (δ⁷Li) by juxtaposing two storm events in an upland watershed subject to a Mediterranean climate. Discharge and δ⁷Li are negatively correlated in both events, but mean δ⁷Li ratios and associated ranges of variation are distinct between them. We apply a previously developed reactive transport model (RTM) for the site to these event-scale flow perturbations, but observed shifts in stream δ⁷Li are not reproduced. To reconcile the stability of the subsurface solute weathering profile with our observations of dynamic stream δ⁷Li signatures, we couple the RTM to a distribution of fluid transit times that evolve based on storm hydrographs. The approach guides appropriate flux-weighting of fluid from the RTM over a range of flow path lengths, or equivalently fluid residence times. This flux-weighted RTM approach accurately reproduces dynamic storm δ⁷Li-discharge patterns distinguished by the antecedent conditions of the watershed.

**Plain Language Summary** Storm events often cause characteristic shifts in stream solute chemistry. Interpreting these signals offers insight into the water-rock interactions occurring within watersheds. Here, we use lithium stable isotopes and reactive transport modeling to relate how long water spends in a catchment, or how deep water infiltrates through a catchment, to the extent of chemical weathering. We show that the first significant storm after a dry season exports more chemically evolved water, while a wet season storm releases less evolved, shallower, and younger water. Our results indicate that stream flow δ⁷Li in small watersheds offers a sensitive record of hydrological conditions prior to the storm, reflecting subtle shifts in the efficiency of the Critical Zone to generate, transport, and ultimately export solutes.

## 1. Introduction

Upland landscapes are commonly characterized by small drainage areas, short water residence times, and seasonal climatic variability (e.g., Anderson et al., 1997; McCormick et al., 2021; Sullivan et al., 2016; Whittaker et al., 1979). Individual storm concentration-discharge (C-Q) dynamics in these low-order systems exhibit pronounced responses in stream chemistry (Knapp et al., 2020; Rose et al., 2018), particularly based on the stable isotope ratios (δ) of weathering-derived solutes (Fernandez et al., 2022; Golla et al., 2022). These factors all suggest the capacity to interpret stream C-Q and δ-Q patterns as an integrated signature of transient environmental perturbations to upland landscapes (Brantley et al., 2023; Druhan & Benettin, 2023; Fernandez et al., 2022; Golla et al., 2022; Knapp et al., 2022; Li et al., 2022; Rose et al., 2018).

A recent steady-state reactive transport model (RTM) study (Winnick et al., 2022) indicated characteristically negative C-Q patterns in water dissolved lithium (Li) concentrations, corresponding to a weaker enrichment of isotopically heavy Li (lower δ⁷Li) as discharge increases. This model prediction is empirically supported by a survey of rivers and streams showing consistent δ⁷Li variations over seasonal time scales, with lower values during wet seasons (i.e., shorter water residence times) than in dry seasons (Zhang et al., 2022). Altogether, these studies illustrate broad-scale differences in stream δ⁷Li-Q dynamics as a lens into the generation and routing of rock-derived solutes across diverse environments. Here, we seek to extend the utility of this dynamic signal to event timescales using high-frequency data sets coupled to models capable of functioning beyond steady state assumptions.

In what follows, we utilize a small, upland catchment subject to a montane Mediterranean climate in southern France. This setting is ideally suited to produce dynamics in weathering-derived solute exports and associated δ⁷Li-Q patterns during storm events and across wet and dry antecedent conditions. We employ a previously

<footnote n="1">¹Department of Earth Science & Environmental Change, University of Illinois at Urbana-Champaign, Urbana, IL, USA,</footnote>

<footnote n="2">²Now at Physical and Life Sciences Directorate, Lawrence Livermore National Laboratory, Livermore, CA, USA,</footnote>

<footnote n="3">³Université Paris-Cité, Institut de Physique du Globe de Paris, CNRS, Paris, France</footnote>
**Writing – original draft:** Jon K. Golla,
Jennifer L. Druhan

**Writing – review & editing:** Jon K. Golla,
Julien Bouchez, Jennifer L. Druhan

established, isotope-enabled RTM (Golla, Bouchez, Kuessner, & Druhan, 2024) built for the site to explore the capacity to produce stream solute and isotope ratio variations at these event to seasonal timescales. We equate depth in the 1-D vertical weathering profile in the model to characteristic water transit time, which allows us to expand our RTM into a range of fluid flow paths, reproducing observed dynamics in stream chemistry during storms across a range of seasons and antecedent conditions. To our knowledge, this study is the first demonstration of the capacity for a multi-component RTM coupled to a transient fluid residence time distribution to reproduce observed event scale riverine C-δ-Q dynamics. Our approach offers a new bridge between disparate hydrological and (bio)geochemical frameworks commonly used to describe solute transport in catchments (Benettin et al., 2022; Li et al., 2021).

## 2. Methods

### 2.1. Site Description

The field site is the small (0.54 km²), upland (1,150–1,450 m.a.s.l.; 18° slope), granitic catchment of Sapine, covered by European beech (Fagus sylvatica) located on the southern slope of Mont Lozère in the Cévennes region of the French Massif Central. The catchment is part of the French Critical Zone Observatory network OZCAR (Observatoires de la Zone Critique: Applications et Recherche; Gaillardet et al., 2018) through the long-term monitoring program Observatoire HydroMétérologique Cévennes-Vivarais (OHM-CV; Boudevillain et al., 2011) since the 1980s (https://ohmcv.osug.fr/). Composition of the porphyritic granodiorite bedrock and of overlying soils are provided in Golla, Kuessner, Reina, et al. (2024). The proximity of the catchment to the Mediterranean Sea strongly regulates local climate. The annual average temperature is 7°C with mean annual precipitation of 2,000 mm (Martin et al., 2003), which is partitioned into characteristically dry (typically April–August) and wet (typically September–March) seasons. Most of this precipitation is concentrated during storm events in the fall season, typifying the “épisodes cévenols” of the Cévennes region (Marc et al., 2001; Martin et al., 2003) and delivering heavy rainfall to the site over short intervals of time. For example, 358 mm of precipitation fell over a span of 30 min during an event in August 1999 (Martin et al., 2003). Response in stream discharge is immediate and despite such intense and flashy precipitation, runoff is largely routed through subsurface flow paths and the occurrence of overland flow is minimal (Durand et al., 1992; Marc et al., 2001; Martin et al., 2003).

### 2.2. Storm Discharge Sampling

We focus on two storms in this study, which occurred between November 2–4, 2015 (hereafter referred to as the “wet season storm”) and October 12–14, 2016 (hereafter referred to as the “dry season storm”) (Figure 1a). During these events, the creek was sampled at much higher frequency than the supporting long-term geochemical monitoring effort. The wet season storm was characterized by two rain pulses. Sampling for this event began during the recession of the first rain pulse and continued across the second storm hydrograph. Samples were taken hourly through the second pulse, and every 4 hours over the “recovery” period following peak discharge. For the dry season storm, samples were taken every 4 hours leading up to the start of the storm, every hour from the start to the height of the storm hydrograph, and every 4 hours over the recover period. Samples were collected using a Hach Sigma SD900 water sampler deployed a few meters upstream from the gauging station. Water samples were passed through 0.22-µ filters into acid-washed polyethylene bottles within a few hours after collection. An aliquot was treated with distilled 16-M HNO₃ to achieve a pH of ~2 for subsequent analysis of cations and isotope ratios.

### 2.3. Analytical Measurements

Measurements of solute chemistry and lithium isotope ratios were performed using the PARI Analytical Platform of the Institut de Physique du Globe de Paris (IPGP) following the methods reported in Golla, Bouchez, Kuessner, and Druhan (2024). Analytical accuracy for solute concentrations are ≤2% and ≤4% for major and trace elements, respectively, and the uncertainty for δ⁷Li measurements is expressed as the 95% confidence interval (Text S1 in Supporting Information S1).

### 2.4. Reactive Transport Model Development

A full description of the conceptualization, parameterization, and validation of our isotope-enabled RTM for the Sapine catchment using the CrunchTope software can be found in Golla, Bouchez, Kuessner, and Druhan (2024).
<figure><img src="image_3.png" /><figcaption>**Figure 1.** (a) Time series of rainfall (upper light gray bars) and Sapine stream discharge (lower black line) from July 2015 to July 2017, encompassing two full water years. Vertical shading indicates the two storms sampled at high frequency reported in the present study. The wet season (b, green) and dry season (c, orange) storms are individually illustrated with accompanying observations of lithium isotope ratios (circles) and Li/Na ratios (triangles). Both the 95% confidence interval calculated from sample replicate measurements and the 2σ associated with repeat analyses of NASS-6 (Text S1 in Supporting Information S1) are shown as error bars for δ<sup>7</sup>Li.</figcaption></figure>

Briefly, the model was originally developed to describe the present-day geochemical signatures of the local weathering profile and the Sapine streamflow solute chemistry under ambient (i.e., multi-year average) hydrological conditions. A 20-m domain is used to encompass the transit of water from vertical infiltration through shallow soil and the deeper vadose zone to lateral drainage and discharge to the small headwater stream driven by constant meteoric fluid infiltration (Darcy flux of 1.4 m/yr) and solid-phase uplift and erosion (0.0001 m/yr). The model reaches steady state in ~200,000 years, at which point simulation results show agreement with measured soil elemental mass depletion, stream solute compositions, Li/Na ratios and δ<sup>7</sup>Li (Golla, Bouchez, Kuessner, & Druhan, 2024).

The current study makes an important update to this calibrated and data-validated RTM. We subject the prior model to transient hydrologic forcing occurring over the much shorter (days - weeks) duration characteristic of our observed storm events. To allow flow rates to vary smoothly in accordance with measured storm precipitation and drainage rates, we implement an update to the CrunchTope source code which produces smooth variations in flow rate with time. Each storm event is run over a duration of 14 days starting from the steady state described in Golla, Bouchez, Kuessner, and Druhan (2024). This encapsulates 7 days of antecedent conditions, which are unique to each storm, as well as 7 days of the storm event and recovery period. For the wet season storm (November 2015), the antecedent conditions include the first storm pulse which occurred prior to sampling (Figure 1), whereas the antecedent condition for the dry season storm (October 2016) corresponds to a decline in flow rate. The November 2015 storm occurred during the wet season and multiple peaks are observed over the hydrograph (Figure 1). These are mimicked through a sequence of two lognormal curves. The October 2016 storm was the first after a prolonged dry period, and we mimic this antecedent condition using an exponential decline in flow rate. Following this decay in streamflow, the storm is simulated by a single lognormal curve. All storms are simulated using the same shape parameters for the lognormal distribution, while the magnitude is varied to
produce changes in flow through time resembling a short and steep rising limb and a relatively long and gradually declining falling limb, consistent with the Sapine hydrograph (Figure 1). Further details of the storm parameterizations can be found in Text S2 in Supporting Information S1.

## 3. Results

### 3.1. Storm Characteristics and Discharge Patterns

Despite the different antecedent conditions of these two storm events (Figure 1), the duration and magnitude of precipitation are fairly similar. 97 mm of precipitation were delivered over 38 hours during the November 2015 wet-season storm (Figure 1b), whereas a total of 126 mm of rain fell on the catchment over a span of 42 hours in the dry-season storm (Figure 1c). Resulting event-averaged rainfall rates are 2.6 and 3.0 mm/hr for the wet-season and dry-season events, respectively.

Although the precipitation characteristics of the two storms are similar, the resulting stream discharge patterns are significantly different (Figure 1). In the wet season, pre-storm discharge rates were approximately $3.8 \times 10^{-2} m^3/s$ (Figure 1b), whereas in the dry season pre-storm values were much lower, averaging approximately $5.0 \times 10^{-4} m^3/s$ (Figure 1c). This 100-fold difference in pre-storm discharge rates highlights the bimodal nature of the montane Mediterranean climate which characterizes the Cévennes region. Such distinct antecedent conditions led to discharge rates in the dry season storm roughly an order of magnitude lower than in the wet season, although the dry season storm featured higher rainfall rates. Across both storm events, there is minimal to no lag time between maximum rainfall and peak discharge. This relatively instantaneous response is consistent with observations at other small, “flashy” catchments (e.g., Anderson et al., 1997; Hooper et al., 1990; McDowell & Likens, 1988).

### 3.2. Lithium Signatures

For our current purpose, correction for atmospheric inputs is unnecessary given that our RTM directly incorporates the geochemical composition of rain water in the upper boundary condition and thus factors the contribution of this input into the modeled subsurface fluid signatures that ultimately produce streamflow chemistry. Furthermore, typical rainwater corrections (i.e., assuming a marine signature) to stream chemistry are not applicable at this site given the strong influence of wet and dry dust deposition (Golla, Bouchez, Kuessner, & Druhan, 2024).

Sapine stream dissolved Li concentrations for the November 2015 (0.078–0.9 µM) and October 2016 (0.066–0.094 µM) events (Figure S1 in Supporting Information S1) are consistent with background Li concentrations in Sapine based on monthly and inter-annual sampling frequencies (0.07–0.14 µM; Golla, Bouchez, Kuessner, & Druhan, 2024). However, across the storm hydrograph, Li concentrations decrease during the wet season storm (Figure S1a in Supporting Information S1) and increase during the dry season storm (Figure S1b in Supporting Information S1). Normalization to Na concentration facilitates comparison of observations to model results given that the model does not explicitly account for the influence of direct rain input to stream solute concentrations. The largely monotonic trends in Li concentrations are replaced by a more nuanced pattern of Li/Na ratios that are closely associated with the storm hydrograph (Figures 1b and 1c). $\delta^7$Li values recorded during both the wet season (+21.1‰ to +22.5‰) and dry season (+22.7‰ to +25.2‰) storms are generally consistent with values observed over long-term sampling campaigns (+22.4‰ to +24.2‰; Golla, Bouchez, Kuessner, & Druhan, 2024), but the maximum and minimum storm $\delta^7$Li ratios expand this long-term range. Both storm events produce an inverse relationship between $\delta^7$Li and discharge (Figures 1b and 1c). The highest $\delta^7$Li ratio during the wet season storm is observed in between the two discharge peaks (Figure 1b) whereas the maximum values during the dry season event occur in the very early stages of the storm when discharge is still quite low (Figure 1c). The net decrease in $\delta^7$Li from baseflow to peak discharge is larger during the dry season storm (~2.5‰) than the difference recorded during the wet season storm (~1.4‰).

### 3.3. Model Results

The initial geochemical conditions of our RTM simulation are provided by the steady state spatial profiles of the dissolved load across the 1-D, 20-m flow path (Figure 2; Golla, Bouchez, Kuessner, & Druhan, 2024). The modeled Li/Na starts from the rainwater boundary condition (1.2 µM/mM) and remains relatively steady over the first ~5 m of the model domain. In tandem across the same 5-m interval, the simulated $\delta^7$Li begins from a value
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Model results for the wet season (a) and dry season (b) storm simulations. Each set of results consists of a time series of the model Darcy flux (black dashed line; Equations S3 and S4 in Supporting Information S1) above three vertical spatial profiles of fluid Li/Na, lithium isotope ratios, and flux weighting (Equation S9 in Supporting Information S1). The first 7 days of the model run correspond to the simulated antecedent hydrological conditions associated with each storm event while the latter 7 days in the simulation mimic the measured storm event itself (Section 2.4). The set of colored circles in the Darcy flux time series plot denotes different points in time during the simulation. For the wet season storm, the coloring highlights the increase (from dark purple to yellow) in flow rates during a given storm pulse. For the dry season, the coloring illustrates the recession of flow rates (from dark purple to yellow) prior to the storm. Vertical spatial profiles of Li/Na and $\delta^7$Li are colored using the same gradient, but these are essentially indiscernible given the lack of variation in the 1-D RTM solute chemistry. Gray vertical shading of these spatial profiles illustrate the range of stream observations across both storm events. Note that the simulated spatial profiles during the second storm pulse in the modeled wet season event are identical to those generated and highlighted from the first storm pulse.</figcaption></figure>

(+10.4‰) resembling that of the prescribed rainwater boundary condition (+8.6‰) and rapidly reaches the highest value (~+26‰) across the entire model domain. The Li behavior simulated across this initial section of the model domain is a result of the balance between dissolution of minerals supplying Li to the solution and the cumulative formation of kaolinite clay that takes up Li and fractionates its isotopes. Further into the model domain, Li/Na increases to ~1.7 µM/mM and $\delta^7$Li decreases to ~+23‰. These trends are associated with the progressively decreasing extent of secondary clay accumulation and thus a greater effect of primary mineral dissolution. Subjecting the RTM to variations in water flow associated with both the wet season (Figure 2a) and dry season (Figure 2b) storms has no discernible effect on the modeled Li/Na or $\delta^7$Li fluid ratios across the 1-D depth profile.

# 4. Discussion

## 4.1. Attenuation of Events Along a Flowpath

In our RTM framework, all water takes the same amount of time to transit the domain and hence, for a given flow rate, discharge is of a single age. The baseflow Darcy flux (1.4 m/yr) and porosity (~10%) allows fluid to transit the 20-m domain in 1.43 years. At the peak Darcy flux modeled during the wet season storm (9.1 m/yr), this transit time is still ~80 days. Even if this peak flow rate were maintained for the entire 7-day duration of a storm event, the fluid infiltrating the domain could only traverse ~1% of the entire spatial profile (i.e., the first two grid cells). The short duration of these events, combined with the dilute solute composition of the Sapine catchment and the typical timescales of reactivity between water and silicate minerals collectively yield model simulations in which pore water chemistry is insensitive to such short-term perturbations in flow rate. The result is that the geochemical variations measured in Sapine creek during these two storms (Figure 1) are not a result of changes in the solute composition of near-surface water due to the storm events. This disconnect between observations and model results necessitates an alternative model representation of the dynamic storm solute signatures produced in this upland watershed.
Piston-like behavior, that is, all fluid moving at a single velocity along a common flowpath, does not characterize the complexity of many natural watersheds (McGlynn et al., 2003; McGuire & McDonnell, 2006; Sukhija et al., 2003). In particular, during storm events, small upland catchments are thought to increasingly deviate from piston flow due to their subsurface heterogeneity (Harman et al., 2010). Rapid attenuation of storm pulses in the simulated Li depth profiles of our 1-D RTM (Figure 2) further support this inference and suggest that the Li dynamics observed in the stream must result from event-scale shifts in the way subsurface water is released to streamflow due to these storms. This may result from alteration in runoff generation pathways (e.g., preferential flow through fractures and shallow porous layers), an increase in catchment wetness or hydrological connectivity, and/or a change in the relative contribution of a range of flow paths through the system. As a basic representation of these complex and coupled flow routing mechanisms, in the following we consider the capacity for (a) a variety of fluid transit times to contribute to streamflow and (b) dynamics in these fluid transit time distributions (TTDs) over the course of the storms.

## 4.2. Stream Chemistry as a Mixture of Flowpaths

TTDs appropriate to upland watersheds are well established in the literature (Benettin et al., 2022; Heidbüchel et al., 2012; Hrachowitz et al., 2016; McGuire & McDonnell, 2006; Sprenger et al., 2019) and their utility in describing event-scale C-Q and δ-Q patterns is becoming increasingly evident (e.g., Benettin et al., 2017; Fernandez et al., 2022; Torres & Baronas, 2021). Our novel contribution is to connect such fluid TTDs to the geochemical rigor offered by our multi-component RTM simulations. To do so, we make a simple assumption in which depth across the 1-D RTM model is converted into corresponding travel times for a given rate of fluid draining through the profile (Equation S5 in Supporting Information S1). With this translation, we may think of each of the domain grid cells as the end of an individual flow path with its own characteristic fluid transit time. In isolation, these flow paths act as independent 1-D domains operating in parallel. Using an appropriate TTD as a guide, we may then decide how to sample across these flow paths to produce a flux-averaged ensemble of fluid travel times which contribute to streamflow. From this starting point, the parameters of the TTD may be reasonably varied such that the behavior of the storm events are described. Doing so allows transient shifts in the way water is subsampled out of the 1-D RTM domain, guided by the geochemical dynamics observed in Sapine creek during these storm events. While the dissolved Li signatures within the model domain do not change over the timescale of an event (Figure 2), the way in which we mix these waters to produce stream chemistry does. We strongly emphasize that this approach is predicated on the underlying spatial and temporal stability of subsurface solute and fluid isotope ratio profiles. Where perturbations to the system are of sufficient duration and/or intensity to alter pore water chemistry, then more advanced RTM approaches are needed (discussed further Supplementary Text S3 and Figure S3 in Supporting Information S1).

Empirical evidence across multiple catchments suggests TTDs are commonly biased toward young water (i.e., short flow paths) but also heavily skewed by a long tail (Godsey et al., 2010; Heidbüchel et al., 2013; Hrachowitz et al., 2010; Jasechko et al., 2016; Kirchner et al., 2000; McGuire & McDonnell, 2010; Segura et al., 2012). In keeping with these observations, we apply a simple exponential TTD inherently predisposed to a positively skewed distribution (Equation S6 in Supporting Information S1; Druhan & Maher, 2017; Fernandez et al., 2022; Maher, 2011). The functional form of the TTD is prescribed, but the mean transit time, τ(t), is dynamic following the behavior of the storm events (Figure 2). The value of τbaseflow is ~8.5 months, which already produces a TTD that is relatively skewed toward the shortest flow paths. In the simulated wet season storm, the changes in TTD profile with time are identical across the two modeled storm pulses (Figure 2a). As flow increases over the course of a simulated storm pulse, τ(t) decreases and the TTD profile skews even further toward short flow paths. At peak flow, τpeak = ~6 days. In contrast, the initial TTD profile of the simulated dry season storm first increases from the same τbaseflow starting point, reaching a value sufficient to cause the TTD profile to approach nearly uniform flux-weighting of flow paths across the model domain (Figure 2b). The absolute value of this pre-storm τdry (~7,143 years) is essentially arbitrary given that any τ(t) > 30 years creates uniform flux-weighting along the depth profile and hence no further variation in modeled stream solute chemistry or isotope ratios (Figure S5 in Supporting Information S1). The onset of the subsequent storm causes τ(t) to decrease following TTD behavior similar to the simulated wet season storm. However, the preference for water from short flow paths is weaker in this dry season condition as a result of the lower peak discharge rate. This occurs even though approximately the same total rainfall was received over comparable time intervals for the two storms (Section 3.1; Figure 1).
<figure><img src="image_5.png" /></figure>

**Figure 3.** Model-data comparison of fluid lithium isotope ratios (a) and Li/Na ratios (b) with flow rate. Measured stream chemistry (circles and triangles) are plotted against the corresponding discharge of Sapine Creek normalized by the maximum value of the storm events. Both the 95% confidence interval calculated from sample replicate measurements and the $2\sigma$ associated with repeat analyses of NASS-6 (Text S1 in Supporting Information S1) are shown as error bars for $\delta^7$Li. Model results correspond to the flux-weighted fluid lithium compositions (orange line for wet season simulation and green line for dry season simulation; Equation S10 in Supporting Information S1) and the model Darcy flux is again normalized to the maximum value. Such normalization to maximum flow rate is to facilitate juxtaposition of model results and observations given that the dynamic subsurface flow rate in the simulated storm events are parameterized to capture the same order-of-magnitude variations observed in stream discharge. The shaded regions surrounding these modeled $\delta^7$Li and Li/Na profiles represent the solution space that captures the range of stream observations considering a plausible range of variation in the Li isotope fractionation factor and degree of Li partitioning in clay (Golla, Bouchez, Kuessner, & Druhan, 2024). The black arrow in panel (b) illustrates an inferred trend associated with enhanced dust dissolution during high-discharge conditions (Section 4.2).

The range of RTM depths is subsampled using the storm-driven variations in mean transit time to produce a flux-weighted solute signature that varies over the storm events (Figure 2). In doing so, dynamics in stream $\delta^7$Li are immediately observable (Figure S2 in Supporting Information S1). In general, this combined TTD-RTM model produces the characteristically inverse relationship between $\delta^7$Li and discharge previously documented in streams and rivers (Golla et al., 2022; Zhang et al., 2022). This systematic trend is also in agreement with our observations at Sapine (Figure S2 in Supporting Information S1). The starting $\delta^7$Li in both wet and dry season simulations (+23.6‰) is slightly higher than the piston-flow $\delta^7$Li value at the base of the 1-D RTM profile (+23.2‰; Figure 2) due to preferential weighting of somewhat shorter flow paths, where the RTM $\delta^7$Li ratio is generally $\ge$+25‰, over the spatial interval between 2.5 and 18 m. When this flux-weighted sampling is subjected to the behavior of the wet season storm, peak discharge yields a $\delta^7$Li of +19.6‰ (Figure S2 in Supporting Information S1). In the dry season storm, the pre-event low-flow period holds water transit times to longer values even during the storm, such that $\delta^7$Li is generally more stable and only produces a minimum of +21.6‰ at peak discharge. Corresponding shifts in flux-weighted Li/Na ratios during the modeled storm simulations are generally bounded between values of 1.2 and 1.5 µM/mM (Figure S2 in Supporting Information S1). As a result, any changes in these values that occur during periods of high flow rates are extremely small. This is a result of the narrow range of Li/Na values produced in the RTM across the spatial domain, which negates changes in the flux-weighting of water during high flow conditions (Figure 2).

We do not attempt to fit the detailed characteristics of the individual storm $\delta^7$Li-Q or Li/Na-Q patterns beyond what is achieved by simply applying an appropriate scaling relationship between $\tau(t)$ and flow rate and mapping the resulting time-varying TTDs onto the 1-D RTM profile. However, direct comparison of the combined model output to our observations in Sapine Creek (Figure 3) are generally quite robust. For this catchment, we find that $\tau \propto \frac{1}{q(t)}(1 + \log(\frac{q_{\text{dunec}}}{q(t)}))$ offers an accurate linkage between variations in discharge and the emergent flux-weighting of our solute profiles such that the measured data are reproduced with minimum parameterization (see Supplementary Text S4, Equations S7 and S8, and Figure S4 in Supporting Information S1 for detailed explanation). Both the absolute value and the shift in relation to discharge for Li/Na and $\delta^7$Li depend on the background hydrologic conditions preceding a given event. When the catchment is dry, flow paths and fluid transit times are
relatively long. As flow rates increase from these initial hydrologic conditions, the coupled TTD-RTM accesses fluid from a more isotopically enriched portion of the spatial profile, resulting in a small (~0.2‰) increase in the observed δ⁷Li (Figure 3a). However, when the catchment is already wet, flux-weighted δ⁷Li ratios are more biased toward relatively enriched fluid from these shorter flow paths. Pushing the system even further toward increasingly short fluid travel times guides a preference toward fluid from the shortest fluid flow paths, where δ⁷Li ratios are strongly influenced by the upper boundary condition with characteristically low values (Golla, Bouchez, Kuessner, & Druhan, 2024). This produces a sharp decrease in streamflow δ⁷Li of ~4‰, which is mirrored by measured values in the stream during the wet season event (Figure 3a).

In comparison to δ⁷Li, the agreement between observations and modeled Li/Na (Figure 3b) is generally good at low flow conditions, but the behavior during the storms is not as well captured. Sapine streamflow Li/Na shows an increase with discharge of ~1 µM/mM whereas the simulations project a small decrease of ~0.2 µM/mM over the same range of flow rates. This misfit is likely associated with variability in Saharan dust inputs that influence local rainwater geochemical composition (Golla, Bouchez, Kuessner, & Druhan, 2024). Variation in the prescribed rainwater upper boundary condition of our RTM can propagate as deep as ≤5 m (Golla, Bouchez, Kuessner, & Druhan, 2024), and these short flow paths are heavily weighted in the flux-averaged streamflow chemistry during intermediate to high discharge events (Figure 2). Even a small enhancement in dust dissolution (Li/Na = 65 µM/mM; Clergue et al., 2015) during these storm events could create a significant increase in Li/Na. Such an effect is also consistent with the associated decrease in observed δ⁷Li given the crust-like signature of Saharan dust (δ⁷Li = +0.7‰; Clergue et al., 2015). Significant solubilization of dust leading to higher solute concentrations in rainwater has been suggested to occur over brief periods of time as short as 2 hours (Desboeufs et al., 1999). This is well within the timescales of the observed storms in Sapine (Figure 1). Presently, we lack the temporal resolution in precipitation samples necessary to account for such boundary condition dynamics in the model, but the propagation of a more dust-influenced rainwater signature over the course of a storm simulation would not affect the general trend of decreasing flux-averaged model δ⁷Li with increasing flow rates (Figure 3a) since fluid from the ≤5-m flow paths affected by the boundary condition are already associated with the lowest δ⁷Li values in the model domain (Figure 2).

### 4.3. Implications for Fluid Routing and Solute Export During Storms

Our coupled TTD-RTM model offers a mechanistic underpinning for the observed variability in streamflow δ⁷Li (Figure 3) recorded in Sapine creek across both wet and dry season storms. These results demonstrate the importance of antecedent conditions in subsurface fluid routing during storm events. Under dry conditions, baseflow is characteristically maintained by relatively long or deep flowpaths as represented in our model by a more even flux-weighting across various depths in the 1-D RTM profile. This would suggest that dry season streamflow is composed of fluid that has generally taken longer to drain through the watershed. In comparison, relatively wet periods of the year route more water into storage, enhance hydrologic connectivity, and generally sustain higher flow rates. Under such conditions, a higher water table combines with enhanced connectivity to route more fluid to the stream, and this supply is generally skewed toward shorter flow paths (Figure 2a). From either starting point, imposing a storm event on the system increases drainage rates and shortens transit times over a timescale of days to weeks. The corresponding variation in δ⁷Li during the wet season storm was fairly small (difference between maximum and minimum values ~1.5‰; Figure 1b) and around distinctly lower values (average 21.8‰) than those of the dry season (average 24.0‰), because the characteristically longer flow paths of the drier season sustain higher baseflow δ⁷Li. Furthermore, although the peak discharge of the dry season storm was ~50% smaller than that of the wet season storm, the corresponding variation in δ⁷Li across the dry season event was much larger (difference between maximum and minimum values of ~2.4‰; Figure 1c).

These non-uniform and dynamic responses are a result of the sensitivity of Li isotopes to characteristic weathering reactions and therefore highlight their utility across variations in discharge and seasonal antecedent conditions (Golla et al., 2022; Zhang et al., 2022). Therefore, dissolved Li isotope ratios are markers of distinct compartments across a weathering profile, such that the resulting composition of streamflow could be used to infer the flow path(s) predominantly contributing to drainage (Golla et al., 2021, 2022). The hydrogeochemistry of Sapine produces a sharp contrast between low-δ⁷Li water in the shallowest, most highly weathered portions of the subsurface and higher-δ⁷Li water which evolves over longer flow paths across less weathered regolith (Figure 2).
Although less pronounced, distinctions between intermediate flow paths (2.5 m ≤ x ≤ 12.5 m; δ⁷Li = 25–26‰) and longer length scales (>12.5 m; δ⁷Li = 24–25‰) in the modeled spatial profile are also diagnostic and discernible within the stream δ⁷Li signatures. Such flow paths characterize sections of the weathering profile where secondary clay precipitation is increasingly balanced by primary mineral dissolution (Golla, Bouchez, Kuessner, & Druhan, 2024).

Variations in stream δ⁷Li as demonstrated here have already been documented in similar upland sites (Golla et al., 2021, 2022; Zhang et al., 2022). Furthermore, such dynamics have been illustrated in other metal(loid) isotope systems which are not exclusively fractionated by secondary clay formation. A systematic trend between Si isotope ratios (δ³⁰Si) and discharge was shown in Sapine Creek for the same dry season storm reported herein, which should occur as an effect of preferential routing of water from shorter flow paths, as well as the combination of fractionating pathways associated with clay formation and biological uptake (Fernandez et al., 2022). Such consistency between isotope systems encourages expansion and use of δ-Q observations alongside more common hydrological and biogeochemical (i.e., C-Q) data to better understand the response of the Critical Zone to transient perturbations such as storms.

## 5. Conclusions

This study addresses a disconnect between modern multi-component RTMs and the transient hydrological routing of water hosted in structurally complex catchments during short-term infiltration events. Widely available RTM software affords a critical capability to represent geochemical structure in both solute and mineral compositions along a given flowpath. This spatial resolution is necessary to offer a mechanistic underpinning for element-specific geochemically distinct stores of solutes within the near surface environment. However, in these frameworks, periodic short-term pulses or declines in long-term average flow rates are rapidly attenuated with depth into the profile, which prohibits characteristically dynamic storm C-Q and δ-Q patterns. Our mapping of a transient fluid TTD onto structured RTM profiles resolves this issue through a simple space-for-time substitution, enabling representation of a network of discrete flow paths. Merging of these two frameworks offers a more appropriate means of representing fluid flow and geochemical reactivity, as called for in recent synthesis papers (Benettin et al., 2022; Li et al., 2021). Here, this integration highlights the potential of streamflow Li isotope ratios as a tracer of hydrological routing through catchments during storm events.

## Data Availability Statement

This work builds upon previously published data (Golla, Kuessner, et al., 2024) and CrunchTope reactive transport models developed by Golla, Bouchez, Kuessner, and Druhan (2024). The associated model input and output files can be accessed on Zenodo (Golla, Bouchez, & Druhan, 2024).

## References

* Anderson, S. P., Dietrich, W. E., Torres, R., Montgomery, D. R., & Loague, K. (1997). Concentration-discharge relationships in runoff from a steep, unchanneled catchment. *Water Resources Research*, 33(1), 211–225. https://doi.org/10.1029/96WR02715

* Benettin, P., Bailey, S. W., Rinaldo, A., Likens, G. E., McGuire, K. J., & Botter, G. (2017). Young runoff fractions control streamwater age and solute concentration dynamics. *Hydrological Processes*, 31(16), 2982–2986. https://doi.org/10.1002/hyp.11243

* Benettin, P., Rodriguez, N. B., Sprenger, M., Kim, M., Klaus, J., Harman, C. J., et al. (2022). Transit time estimation in catchments: Recent developments and future directions. *Water Resources Research*, 58(11), e2022WR033096. https://doi.org/10.1029/2022WR033096

* Boudevillain, B., Delrieu, G., Galabertier, B., Bonifait, L., Bouilloud, L., Kirstetter, P.-E., & Mosini, M.-L. (2011). The Cévennes-Vivarais mediterranean hydrometeorological observatory database. *Water Resources Research*, 47(7), W07701. https://doi.org/10.1029/2010WR010353

* Brantley, S. L., Shaughnessy, A., Lebedeva, M. I., & Balashov, V. N. (2023). How temperature-dependent silicate weathering acts as Earth's geological thermostat. *Science*, 379(6630), 382–389. https://doi.org/10.1126/science.add2922

* Clergue, C., Dellinger, M., Buss, H. L., Gaillardet, J., Benedetti, M. F., & Dessert, C. (2015). Influence of atmospheric deposits and secondary minerals on Li isotopes budget in a highly weathered catchment, Guadeloupe (Lesser Antilles). *Chemical Geology*, 414, 28–41. https://doi.org/10.1016/j.chemgeo.2015.08.015

* Desboeufs, K. V., Losno, R., Vimeux, F., & Cholbi, S. (1999). The pH-dependent dissolution of wind-transported Saharan dust. *Journal of Geophysical Research*, 104(D17), 21287–21299. https://doi.org/10.1029/1999JD900236

* Druhan, J. L., & Benettin, P. (2023). Isotope ratio – Discharge relationships of solutes derived from weathering reactions. *American Journal of Science*, 323, 5. https://doi.org/10.2475/001c.84469

* Druhan, J. L., & Maher, K. (2017). The influence of mixing on stable isotope ratios in porous media: A revised Rayleigh model. *Water Resources Research*, 53(2), 1101–1124. https://doi.org/10.1002/2016WR019666

* Durand, P., Robson, A., & Neal, C. (1992). Modelling the hydrology of submediterranean montane catchments (Mont-Lozère, France) using TOPMODEL: Initial results. *Journal of Hydrology*, 139(1), 1–14. https://doi.org/10.1016/0022-1694(92)90191-W

## Acknowledgments

We are especially grateful to Marie Kuessner, Pierre-Alain Ayral, Jean-François Didon-Lescot, Jean-Marc Domergue, Nadine Grard, Didier Josselin, Philippe Martin (UMR ESPACE), and Yannick Manche (PNC), for support during field work and insightful discussions. The "Parc National des Cévennes" (PNC) is acknowledged for allowing us to conduct scientific work within the park. We thank Pierre Maffre, and a second anonymous reviewer for their constructive feedback, as well as GRL editor Hari Rajaram for his handling of the paper. Funding support is provided by NSF-EAR-2047318 awarded to J.L.D. Sample analysis was supported by the People Programme (Marie Curie Actions) of the European Union's Seventh Framework Programme FP7/2007–2013/ under REA agreement [608069] (ITN "IsoNose"), the IPGP multidisciplinary program PARI (with analytical help from Pascale Louvat, Julien Moureau, Jessica Dallas, Thibaud Sondag, Caroline Gorge, Pierre Burckel, and Laëticia Faure), and the Region Île-de-France SESAME Grant 12015908. J.K.G. acknowledges support from the NSF Graduate Research Fellowship Program. Part of this work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344 (LLNL-JRNL-866541).
Fernandez, N. M., Bouchez, J., Derry, L. A., Chorover, J., Gaillardet, J., Giesbrecht, I., et al. (2022). Resiliency of silica export signatures when low order streams are subject to storm events. *Journal of Geophysical Research: Biogeosciences*, 127(5), e2021JG006660. https://doi.org/10.1029/2021JG006660

Gaillardet, J., Braud, I., Hankard, F., Anquetin, S., Bour, O., Dorfliger, N., et al. (2018). OZCAR: The French network of critical zone observatories. *Vadose Zone Journal*, 17(1), 1–24. https://doi.org/10.2136/vzj2018.04.0067

Godsey, S. E., Aas, W., Clair, T. A., de Wit, H. A., Fernandez, I. J., Kahl, J. S., et al. (2010). Generality of fractal 1/f scaling in catchment tracer time series, and its implications for catchment travel time distributions. *Hydrological Processes*, 24(12), 1660–1671. https://doi.org/10.1002/hyp.7677

Golla, J. K., Bouchez, J., Kuessner, M. L., Rempe, D. M., & Druhan, J. L. (2022). Subsurface weathering signatures in stream chemistry during an intense storm. *Earth and Planetary Science Letters*, 595, 117773. https://doi.org/10.1016/j.epsl.2022.117773

Golla, J. K., Kuessner, M. L., Henehan, M. J., Bouchez, J., Rempe, D. M., & Druhan, J. L. (2021). The evolution of lithium isotope signatures in fluids draining actively weathering hillslopes. *Earth and Planetary Science Letters*, 567, 116988. https://doi.org/10.1016/j.epsl.2021.116988

Golla, J. K., Bouchez, J., & Druhan, J. L. (2024). Accompanying model files to Golla et al. (2024, GRL) "Antecedent hydrologic conditions reflected in stream lithium isotope ratios during storms". *Zenodo*. https://doi.org/10.5281/zenodo.13371813

Golla, J. K., Bouchez, J., Kuessner, M. L., & Druhan, J. L. (2024). Weathering incongruence in mountainous mediterranean climates recorded by stream lithium isotope ratios. *Journal of Geophysical Research: Earth Surface*, 129(3), e2023JF007359. https://doi.org/10.1029/2023JF007359

Golla, J. K., Kuessner, M. L., Reina, C. A., Grard, N., Domergue, J.-M., Ayral, P.-A., et al. (2024). Geochemistry of a weathering profile and natural waters in the upland Mediterranean catchment of Sapine, Mt-Lozère, southern France. *Zenodo*. https://zenodo.org/doi/10.5281/zenodo.10794566

Harman, C. J., Reeves, D. M., Baeumer, B., & Sivapalan, M. (2010). A subordinated kinematic wave equation for heavy-tailed flow responses from heterogeneous hillslopes. *Journal of Geophysical Research*, 115(F1), F00A08. https://doi.org/10.1029/2009JF001273

Heidbüchel, I., Troch, P. A., & Lyon, S. W. (2013). Separating physical and meteorological controls of variable transit times in zero-order catchments. *Water Resources Research*, 49(11), 7644–7657. https://doi.org/10.1002/2012WR013149

Heidbüchel, I., Troch, P. A., Lyon, S. W., & Weiler, M. (2012). The master transit time distribution of variable flow systems. *Water Resources Research*, 48(6), W06520. https://doi.org/10.1029/2011WR011293

Hooper, R. P., Christophersen, N., & Peters, N. E. (1990). Modelling streamwater chemistry as a mixture of soilwater end-members — An application to the Panola Mountain catchment, Georgia, U.S.A. *Journal of Hydrology*, 116(1), 321–343. https://doi.org/10.1016/0022-1694(90)90131-G

Hrachowitz, M., Benettin, P., van Breukelen, B. M., Fovet, O., Howden, N. J., Ruiz, L., et al. (2016). Transit times—The link between hydrology and water quality at the catchment scale. *WIREs Water*, 3(5), 629–657. https://doi.org/10.1002/wat2.1155

Hrachowitz, M., Soulsby, C., Tetzlaff, D., & Speed, M. (2010). Catchment transit times and landscape controls—Does scale matter? *Hydrological Processes*, 24(1), 117–125. https://doi.org/10.1002/hyp.7510

Jasechko, S., Kirchner, J. W., Welker, J. M., & McDonnell, J. J. (2016). Substantial proportion of global streamflow less than three months old. *Nature Geoscience*, 9(2), 126–129. https://doi.org/10.1038/ngeo2636

Kirchner, J. W., Feng, X., & Neal, C. (2000). Fractal stream chemistry and its implications for contaminant transport in catchments. *Nature*, 403(6769), 524–527. https://doi.org/10.1038/35000537

Knapp, J. L. A., Li, L., & Musolff, A. (2022). Hydrologic connectivity and source heterogeneity control concentration–discharge relationships. *Hydrological Processes*, 36(9), e14683. https://doi.org/10.1002/hyp.14683

Knapp, J. L. A., von Freyberg, J., Studer, B., Kiewiet, L., & Kirchner, J. W. (2020). Concentration–discharge relationships vary among hydrological events, reflecting differences in event characteristics. *Hydrology and Earth System Sciences*, 24(5), 2561–2576. https://doi.org/10.5194/hess-24-2561-2020

Li, L., Stewart, B., Zhi, W., Sadayappan, K., Ramesh, S., Kerins, D., et al. (2022). Climate Controls on River Chemistry. *Earth's Future*, 10(6), e2021EF002603. https://doi.org/10.1029/2021EF002603

Li, L., Sullivan, P. L., Benettin, P., Cirkpa, O. A., Bishop, K., Brantley, S. L., et al. (2021). Toward catchment hydro-biogeochemical theories. *WIREs Water*, 8(1), e1495. https://doi.org/10.1002/wat2.1495

Maher, K. (2011). The role of fluid residence time and topographic scales in determining chemical fluxes from landscapes. *Earth and Planetary Science Letters*, 312(1), 48–58. https://doi.org/10.1016/j.epsl.2011.09.040

Marc, V., Didon-Lescot, J.-F., & Michael, C. (2001). Investigation of the hydrological processes using chemical and isotopic tracers in a small Mediterranean forested catchment during autumn recharge. *Journal of Hydrology*, 247(3), 215–229. https://doi.org/10.1016/S0022-1694(01)00386-9

Martin, C., Didon-Lescot, J.-F., & Cosandey, C. (2003). Le fonctionnement hydrologique des petits bassins versants granitiques du Mont-Lozère: Influence du couvert végétal sur les crues et les étages. *Etudes de géographie physique*, 30, 3–25.

McCormick, E. L., Dralle, D. N., Hahm, W. J., Tune, A. K., Schmidt, L. M., Chadwick, K. D., & Rempe, D. M. (2021). Widespread woody plant use of water stored in bedrock. *Nature*, 597(7875), 225–229. https://doi.org/10.1038/s41586-021-03761-3

McDowell, W. H., & Likens, G. E. (1988). Origin, Composition, and Flux of Dissolved Organic Carbon in the Hubbard Brook Valley. *Ecological Monographs*, 58(3), 177–195. https://doi.org/10.2307/2937024

McGlynn, B., McDonnell, J., Stewart, M., & Seibert, J. (2003). On the relationships between catchment scale and streamwater mean residence time. *Hydrological Processes*, 17(1), 175–181. https://doi.org/10.1002/hyp.5085

McGuire, K. J., & McDonnell, J. J. (2006). A review and evaluation of catchment transit time modeling. *Journal of Hydrology*, 330(3), 543–563. https://doi.org/10.1016/j.jhydrol.2006.04.020

McGuire, K. J., & McDonnell, J. J. (2010). Hydrological connectivity of hillslopes and streams: Characteristic time scales and nonlinearities. *Water Resources Research*, 46(10), W10543. https://doi.org/10.1029/2010WR009341

Rose, L. A., Karwan, D. L., & Godsey, S. E. (2018). Concentration–discharge relationships describe solute and sediment mobilization, reaction, and transport at event and longer timescales. *Hydrological Processes*, 32(18), 2829–2844. https://doi.org/10.1002/hyp.13235

Segura, C., James, A. L., Lazzati, D., & Roulet, N. T. (2012). Scaling relationships for event water contributions and transit times in small-forested catchments in Eastern Quebec. *Water Resources Research*, 48(7), W07502. https://doi.org/10.1029/2012WR011890

Sprenger, M., Stumpp, C., Weiler, M., Aeschbach, W., Allen, S. T., Benettin, P., et al. (2019). The Demographics of Water: A Review of Water Ages in the Critical Zone. *Reviews of Geophysics*, 57(3), 800–834. https://doi.org/10.1029/2018RG000633

Sukhija, B. S., Reddy, D. V., Nagabhushanam, P., & Hussain, S. (2003). Recharge processes: Piston flow vs preferential flow in semi-arid aquifers of India. *Hydrogeology Journal*, 11(3), 387–395. https://doi.org/10.1007/s10040-002-0243-3

Sullivan, P. L., Hynek, S. A., Gu, X., Singha, K., White, T., West, N., et al. (2016). Oxidative dissolution under the channel leads geomorphological evolution at the Shale Hills catchment. *American Journal of Science*, 316(10), 981–1026. https://doi.org/10.2475/10.2016.02
Torres, M. A., & Baronas, J. J. (2021). Modulation of Riverine Concentration-Discharge Relationships by Changes in the Shape of the Water Transit Time Distribution. *Global Biogeochemical Cycles*, **35**(1), e2020GB006694. https://doi.org/10.1029/2020GB006694

Whittaker, R. H., Likens, G. E., Bormann, F. H., Easton, J. S., & Siccama, T. G. (1979). The Hubbard Brook Ecosystem Study: Forest Nutrient Cycling and Element Behavior. *Ecology*, **60**(1), 203–220. https://doi.org/10.2307/1936481

Winnick, M. J., Druhan, J. L., & Maher, K. (2022). Weathering intensity and lithium isotopes: A reactive transport perspective. *American Journal of Science*, **322**(5), 647–682. https://doi.org/10.2475/05.2022.01

Zhang, F., Dellinger, M., Hilton, R. G., Yu, J., Allen, M. B., Densmore, A. L., et al. (2022). Hydrological control of river and seawater lithium isotopes. *Nature Communications*, **13**(1), 3359. https://doi.org/10.1038/s41467-022-31076-y

## References From the Supporting Information

Demars, B. O. L., Friberg, N., & Thornton, B. (2020). Pulse of dissolved organic matter alters reciprocal carbon subsidies between autotrophs and bacteria in stream food webs. *Ecological Monographs*, **90**(1), e01399. https://doi.org/10.1002/ecm.1399

Field, M. S. (2004). Forecasting versus predicting solute transport in solution conduits for estimating drinking-water risks. *Acta Carsologica*, **33**(2). https://doi.org/10.3986/ac.v33i2.295

Flesch, G., Anderson Jr, A., & Svec, H. (1973). A secondary isotopic standard for 6Li/7Li determinations. *International Journal of Mass Spectrometry and Ion Physics*, **12**(3), 265–272. https://doi.org/10.1016/0020-7381(73)80043-9

James, R. H., & Palmer, M. R. (2000). The lithium isotope composition of international rock standards. *Chemical Geology*, **166**(3), 319–326. https://doi.org/10.1016/S0009-2541(99)00217-X

Konrad, C. P., Schmadel, N. M., Harvey, J. W., Schwarz, G. E., Gomez-Velez, J., Boyer, E. W., & Scott, D. (2020). Accounting for temporal variability of streamflow in estimates of travel time. *Frontiers in Water*, **2**, 29. https://doi.org/10.3389/frwa.2020.00029

Kuessner, M. L., Gourgiotis, A., Manhès, G., Bouchez, J., Zhang, X., & Gaillardet, J. (2020). Automated Analyte Separation by Ion Chromatography Using a Cobot Applied to Geological Reference Materials for Li Isotope Composition. *Geostandards and Geoanalytical Research*, **44**(1), 57–67. https://doi.org/10.1111/ggr.12295

Lin, J., Liu, Y., Hu, Z., Yang, L., Chen, K., Chen, H., et al. (2016). Accurate determination of lithium isotope ratios by MC-ICP-MS without strict matrix-matching by using a novel washing method. *Journal of Analytical Atomic Spectrometry*, **31**(2), 390–397. https://doi.org/10.1039/C5JA00231A

Magna, T., Wiechert, U. H., & Halliday, A. N. (2004). Lowblank isotope ratio measurement of small samples of lithium using multiplecollector ICPMS. *International Journal of Mass Spectrometry*, **239**(1), 67–76. https://doi.org/10.1016/j.ijms.2004.09.008

Millot, R., Guerrot, C., & Vigier, N. (2004). Accurate and High-Precision Measurement of Lithium Isotopes in Two Reference Materials by MC-ICP-MS. *Geostandards and Geoanalytical Research*, **28**(1), 153–159. https://doi.org/10.1111/j.1751-908X.2004.tb01052.x

Rosner, M., Ball, L., Peucker-Ehrenbrink, B., Blusztajn, J., Bach, W., & Erzinger, J. (2007). A Simplified, Accurate and Fast Method for Lithium Isotope Analysis of Rocks and Fluids, and 7Li Values of Seawater and Rock Reference Materials. *Geostandards and Geoanalytical Research*, **31**(2), 77–88. https://doi.org/10.1111/j.1751-908X.2007.00843.x

Ryu, J.-S., Jacobson, A. D., Holmden, C., Lundstrom, C., & Zhang, Z. (2011). The major ion, 44/40Ca, 44/42Ca, and 26/24Mg geochemistry of granite weathering at pH=1 and T=25°C: Power-law processes and the relative reactivity of minerals. *Geochimica et Cosmochimica Acta*, **75**(20), 6004–6026. https://doi.org/10.1016/j.gca.2011.07.025

Weynell, M., Wiechert, U., & Schuessler, J. A. (2017). Lithium isotopes and implications on chemical weathering in the catchment of Lake Donggi Cona, northeastern Tibetan Plateau. *Geochimica et Cosmochimica Acta*, **213**, 155–177. https://doi.org/10.1016/j.gca.2017.06.026

Zhang, J.-W., Zhao, Z.-Q., Yan, Y.-N., Cui, L.-F., Wang, Q.-L., Meng, J.-L., et al. (2021). Lithium and its isotopes behavior during incipient weathering of granite in the eastern Tibetan Plateau, China. *Chemical Geology*, **559**, 119969. https://doi.org/10.1016/j.chemgeo.2020.119969