Lithium isotope ratios (δ7Li) of rivers are increasingly serving as a diagnostic of the balance between chemical and physical weathering contributions to overall landscape denudation rates. Here, we show that intermediate weathering intensities and highly enriched stream δ7Li values typically associated with lowland floodplains can also describe small upland watersheds subject to cool, wet climates. This behavior is revealed by stream δ7Li between +22.4 and +23.5‰ within a Critical Zone observatory located in the Cévennes region of southern France, where dilute stream solute concentrations and significant atmospheric deposition otherwise mask evidence of incongruence. The water‐rock reaction pathways underlying this behavior are quantified through a multicomponent, isotope‐enabled reactive transport model. Using geochemical characterization of soil profiles, bedrock, and long‐term stream samples as constraints, we evolve the simulation from an initially unweathered granite to a steady state weathering profile which reflects the balance between (a) fluid infiltration and drainage and (b) bedrock uplift and soil erosion. Enriched stream δ7Li occurs because Li is strongly incorporated into actively precipitating secondary clay phases beyond what prior laboratory experiments have suggested. Chemical weathering incongruence is maintained despite relatively slow reaction rates and moderate clay accumulation due to a combination of two factors. First, reactive primary mineral phases persist across the weathering profile and effectively “shield” the secondary clays from resolubilization due to their greater solubility. Second, the clays accumulating in the near‐surface profile are relatively mature weathering byproducts. These factors promote characteristically low total dissolved solute export from the catchment despite significant input of exogenous dust.
]]>Upland landscapes reach a geomorphic steady state when the mass flux of bedrock exhumation is balanced by a combination of physical erosion and chemical weathering (e.g., Anderson et al., 2021;Riebe et al., 2004;Waldbauer & Chamberlain, 2005). These Critical Zone import and export fluxes operate across disparate timescales. Bedrock uplift and exposure through erosion at the land surface typically occurs over geological timescales (10 3 -10 6 years; Schaller & Ehlers, 2022), whereas infiltration and drainage of meteoric fluid occurs over timescales of water storage and transit through catchments (hours to 10 3 years; Sprenger et al., 2019). Together, these distinct solid and fluid fluxes facilitate a suite of water-rock-life interactions that yield the major weathering products of the Critical Zone: soils and the suspended sediments and solutes carried by streams. The dissolved elemental composition of rivers reflect a flux-weighted average of fluid ages closely tied to the drainage rates and groundwater storage characteristics of the watershed, while soils commonly record a much longer history of weathering profile development. The complementary nature of these shortand long-timescale integrated signals offers a compelling means of characterizing geochemical reactivity of the near-surface, ultimately affording improved quantitative understanding of the coupling between the hydrologic and geomorphic functioning of natural landscapes.
Reactive transport models (RTMs) provide a way to couple such disparate fluid and mineral transit timescales to produce a mass-balance based quantitative description of weathering profiles (Maher & Navarre-Sitchler, 2019). However, the coupled multi-component balance between primary mineral solubilization and secondary mineral formation embedded in these models often leads to a high degree of equifinality when only element abundances are considered. In this regard, the signatures of stable isotope fractionation associated with weathering-derived solutes offers a powerful means of parsing within and among these reactive pathways. RTMs that take advantage of the sensitivity of metal(loid) stable isotope ratios to discern weathering reactions are now emerging (e.g., Fernandez et al., 2019;Golla et al., 2021Golla et al., , 2022;;Guinoiseau et al., 2021;Winnick et al., 2022). However, modern RTMs that account for geomorphic transport of the solid phase tied to a stoichiometrically balanced, multicomponent description of fluid-rock reactivity remain rare (Giambalvo et al., 2002;Guinoiseau et al., 2021;Maher et al., 2009;Navarre-Sitchler et al., 2011). Among these studies, only one has leveraged a metal(loid) stable isotope system to interrogate the development of a weathering profile, revealing a sequence of past weathering stages driven by distinct hydrologic conditions over million-year timescales (Guinoiseau et al., 2021).
Leveraging the appropriate metal(loid) isotope ratio(s) within RTM frameworks presents an opportunity to improve quantitative understanding of weathering processes and their relation to climate in natural systems. Lithium (Li) is particularly well suited to this application, given that its isotopes ( 7 Li and 6 Li) fractionate during the formation of secondary mineral byproducts (Huh et al., 2001;Kısakurek et al., 2004Kısakurek et al., , 2005;;Pistiner & Henderson, 2003;Rudnick et al., 2004;Tomascak, 2004) and are generally considered to be relatively insensitive to plant uptake (Lemarchand et al., 2010;Pogge von Strandmann et al., 2016) with the potential exception of tropical environments (W. Li et al., 2020). This behavior leads to a characteristic relationship between riverine lithium isotope ratios (δ 7 Li) and the weathering intensity (W/D) of a given landscape, which is defined as the proportion of total landscape denudation (D, which is the sum of chemical weathering flux (W) and physical erosive flux (E)) that is attributable to W (Bouchez et al., 2014;Dellinger et al., 2015). At low W/D, highly erosive upland systems covered by thin soils reflect low weathering intensity and are generally limited to congruent dissolution of bedrock (Dellinger et al., 2015;Pogge von Strandmann & Henderson, 2015). In contrast, high W/D conditions promote the development of a thick layer of weathered material reflecting a high weathering intensity, where soluble cations tend to be entirely leached from soils. Both of these regimes yield similarly low δ 7 Li values in rivers (i.e., close to that of bedrock), but for very different Li export fluxes (Bouchez et al., 2013;Dellinger et al., 2015). Intermediate weathering intensities (W ≈ E) that fall between these two extremes produce much more enriched δ 7 Li values in rivers.
To date, correspondence between intermediate weathering intensities and enriched fluid δ 7 Li has largely been documented for large rivers with extensive alluvial plains, where sediments are allowed to interact with water over relatively long time scales (Dellinger et al., 2015). However, a modeling study of Li isotopes at the scale of the Amazon Basin suggests that such enriched river δ 7 Li can be produced within Andean hillslopes rather than in floodplains, depending on the nature of the rock subjected to weathering (Maffre et al., 2020). In addition, high dissolved δ 7 Li export is not exclusively controlled by the residence time of solids within a catchment. A recent isotope-enabled RTM suggested that fluid travel times exert an influential control over the extent of secondary mineral formation within a given system and thus the resulting fluid δ 7 Li composition (Winnick et al., 2022). This inference is empirically supported by hydrologic sensitivity in δ 7 Li over a global assessment of seasonal river Li signatures, illustrating a consistent negative correlation between annual runoff and fluid δ 7 Li across a range of spatial and temporal scales (F. Zhang et al., 2022). These studies imply a variety of intermediate weathering intensity environments and the capacity for RTMs embedding the lithium isotope system to offer quantitative and predictive frameworks for these settings.
Here, we consider tectonically active landscapes where fresh bedrock is relatively shallow (e.g., Anderson et al., 2021;Waldbauer & Chamberlain, 2005), but erosion rates are lower than the high mountain endmember in Dellinger et al. (2015). In these upland systems weathering rates can be strongly regulated by climate (i.e., temperature and runoff; Brantley et al., 2023;Maher & Chamberlain, 2014;Velbel, 1993;West, 2012;West et al., 2005), particularly where precipitation rates are high and thus stream solute concentrations are low. Such dilute solute export regimes can deviate from the typically chemostatic behavior (i.e., negligible to weak dependence of concentration on discharge) of weathering-derived solutes (e.g., Diamond & Cohen, 2018;Godsey et al., 2009;Herndon et al., 2015;Koger et al., 2018;Sullivan et al., 2019;Wlostowski et al., 2018). Under these conditions, even slight alterations in inputs (e.g., precipitation rate, dust deposition), subsurface fluid routing (e.g., variable groundwater table elevation during wet and dry seasons) and climate forcing are capable of influencing stream chemical composition (Brantley et al., 2023). Hence, these environments offer natural laboratories to assess the functioning of landscapes over time and subject to environmental perturbations.
Identifying the processes that underlie incongruent weathering conditions in these dynamic environments from the dilute stream chemistry alone remains difficult and presents an opportunity to leverage the sensitivity associated with δ 7 Li signatures. In addition, such a characterization requires a mechanistic approach that can link the disparate timescales of mineral and fluid residence in the weathering profile. The present study reports a multiyear geochemical monitoring effort over the period 2013-2018 and develops an δ 7 Li-enabled RTM for a granitic hillslope catchment mantled by thin soils in the Cévennes region in southern France that is typified by a cool, wet mountainous Mediterranean climate. We use the model to quantify the causes of apparently strong weathering incongruence manifested in elevated stream δ 7 Li ratios that are difficult to resolve from the dilute concentrations of major solutes alone. In doing so, we expand recent characterization of δ 7 Li dynamics across actively eroding, upland watersheds (Golla et al., 2021(Golla et al., , 2022) ) and offer a forward and predictive modeling framework for the subsurface weathering reaction network in these environments.
Sapine Creek sits 1,150-1,450 m above sea level and drains a 0.54-km 2 catchment with an average slope of 18°l ocated on the southern flank of Mt. Lozère within the Cévennes National Park in Southern France (Figure 1). The site is only ∼80 km north of the Mediterranean Sea, yet the creek feeds into the larger Alignon stream and ultimately the Tarn and Garonne rivers, which discharge to the Atlantic Ocean (Figure 1). Sapine has been continuously monitored since the 1980s and is now part of the Observatoire Hydrométérologique Méditerranéen Cévennes-Vivarais (OHM-CV; Boudevillain et al., 2011) but has recently been integrated into the French Critical Zone Observatory network OZCAR (Observatoires de la Zone Critique: Applications et Recherche; Gaillardet et al., 2018). Approximately 80% of the catchment is covered by a mature beech forest (Fagus sylvatica), which has remained relatively undisturbed for more than 60 years (Cognard-Plancq et al., 2001;Durand et al., 1991). Sapine has been used as a "natural" reference point in prior studies focusing on the hydrologic impacts of anthropogenic encroachment in watersheds (Cognard-Plancq et al., 2001;Didon-Lescot & Martin, 2000;Durand et al., 1992;Lelong et al., 1990;Marc et al., 2001;Martin et al., 2003). Most recently, the stream was included in a international cross-site study of silicon stable isotope signatures recording the biogeochemical dynamics of watersheds during storm events (Fernandez et al., 2022).
Mt. Lozère is underlain by porphyritic granodiorite from the Pont-de-Montvert-Borne complex that is intruded by aplite dykes as a result of hydrothermal alteration upon pluton emplacement (Chauvet et al., 2012;Talbot et al., 2004). At Sapine, the granodiorite is widely observed in outcrops and is composed of large orthoclase phenocrysts surrounded by finer-grained albite and quartz with biotite, chlorite, and enstatite as accessory minerals (Text S1 in Supporting Information S1; Talbot et al., 2004). Soil is 0.6 m thick on average and described as acidic and humiferous with a coarse texture (50%-60% corse sand, 15% fine sand, 20% silt, and 5% clay) (Cognard-Plancq et al., 2001;Marc et al., 2001). Seismic surveys show that the weathering profile extends as deep as 10 m below the land surface (Dangeard et al., 2019). Regional incision rates over the last 4 Ma based on cosmogenic isotopes are 83 +17/-5 m/Ma (Malcles et al., 2020), which is consistent with vertical velocities of <0.5 mm/yr estimated from Global Positioning System stations in the surrounding area (Masson et al., 2019).
The field site is characterized by a mountainous Mediterranean climate with monthly mean temperatures varying between -5°C in January and 17°C in July (Martin et al., 2003). The site receives a mean annual rainfall of 2,000 mm. Previous estimates of evapotranspiration suggest a mean annual value of ∼600 mm (Martin et al., 2003). More recent estimates over the interval 2016-2019 based on a Penman-Monteith relationship yield ∼670 mm (Text S2 in Supporting Information S1). Precipitation is largely delivered by storms which occur between autumn and spring. Discharge in Sapine responds rapidly to these events, yet subsurface flow is the dominant source to the stream even during the storms (Marc et al., 2001;Martin et al., 2003).
Discharge at Sapine is measured by a V-notch weir at the stream gauge (Figure 1) based on an established rating curve. Rain gauges are installed in clearings close to the gauging stations at Sapine and at the neighboring Cloutasse watershed ∼3 km due east within Mt. Lozère. These high-frequency data are publicly available and distributed by the OHM-CV (https://ohmcv.osug.fr/). Stream water samples were collected monthly from 2013 to 2018. Rain water was sampled in 1-month aggregates both in Sapine and the neighboring Cloutasse catchment in Mt. Lozère. Samples were filtered using 0.22-μm cellulose acetate membranes and stored in acid-washed polyethylene bottles. Aliquots for major and trace element analyses were acidified to pH ∼2 with purified 16 M HNO 3 . An 80-cm soil core was obtained in 2014 (∼41 m away from the stream) using a manual auger. In addition, four porphyritic granite samples were collected from nearby outcrops.
Elemental chemistry was analyzed using the Plateau d'Analyses Haute Résolution Analytical (PARI) Platform at the Institut de Physique du Globe de Paris (IPGP). Major cation and anion concentrations were measured on a Dionex using 3,000 and 120 columns, respectively. Trace element concentrations of water samples and all element concentrations for digested soil and rock samples were analyzed on an Agilent 7900 quadrupole ICP-MS. The river water reference material SLRS-5 was repeatedly measured to assess the accuracy of water sample measurements. Replicate measurements of this external standard yielded typical accuracy of ≤2% for major elements and ≤4% for trace elements. The accuracy of soil and bedrock measurements (≤3%) was based on repeated measurements of external reference materials GS-N (Vosges granite; Centre de Recherches Pétrographiques et Géochimiques), BHVO-2 (Hawaiian basalt; United States Geological Survey), JB-2 (Oshima volcano basalt; Geological Survey of Japan), and NIST2709a (San Joaquin soil; National Institute of Standards and Technology). Further characterization of bedrock via X-Ray diffraction is described in Text S1 of Supporting Information S1.
Chemical purification of samples for Li isotopic analyses were carried out at IPGP with a robotic pipetting system "ChemCobOne" (CleanChemLab; https://cleanchemlab.com/) and the protocol developed by Kuessner et al. (2020). The sample preparation and separation protocol is an adaptation of the method established by James and Palmer (2000). In order to eliminate matrix effects associated with organic matter, samples were pre-treated with HF and H 2 O 2 . The column was cleaned with alternate rinsings of 3 M HCl and ultrapure H 2 O and conditioned with 0.2 M HCl. Samples contained in a 0.2-M HCl solution were introduced to the column filled with 4 mL of BioRad AG50-X12 (200-400 mesh) resin for separation and elution of Li.
Analytical work for lithium isotopes was performed at IPGP using a Thermo Fisher Scientific MC-ICP-MS. The purified yield from chromatographic separation was evaporated to near-dryness and taken up in 0.5-M HNO 3 to obtain a 20 ng/mL Li solution and then introduced into an Elemental Scientific APEX-HydroFluoric desolvating nebulizer. Analyses were performed with operating conditions similar to those reported by Millot et al. (2004). Results of sample measurements are reported in standard delta per-mille notation
where std-I and std-II refer to the L-SVEC Li carbonate standard (Flesch et al., 1973) evaluated directly before and after the sample. The analytical uncertainty for each sample is reported as the 95% confidence interval
where σ is the standard deviation of n replicate measurements of a sample and t n-1 is the critical value of a Student's t-distribution over a 95% confidence level. Long-term precision of Li isotope analysis is based on repeated measurements (x ± 2σ) of NASS-6 seawater (δ 7 Li = 30.99 ± 0.48‰; n = 45). We also report analyses of reference materials for basalt (BHVO-2: 4.50 ± 0.15‰, n = 6; JB-2: 4.57 ± 0.01‰, n = 6), granite (GSN; -0.57 ± 0.25‰, n = 15), and soil (NIST-2709a: -0.40 ± 0.60‰, n = 9). Our reference material analyses are in agreement with measurements reported in previous studies (Clergue et al., 2015;Lin et al., 2016;Magna et al., 2004;Rosner et al., 2007;Ryu et al., 2011;Weynell et al., 2017).
We use mass transfer coefficients or "tau values" (Anderson et al., 2002;Brimhall & Dietrich, 1987) as a metric for the relative loss or accumulation of a mobile element relative to the abundance of an immobile element in the weathering profile
where j is any (mobile) element, i is the reference immobile element, w is the soil or regolith, and p is the parent bedrock. Titanium (Ti) is used as the reference immobile element, consistent with prior studies in granite lithology (e.g., Ackerer et al., 2016;Hayes et al., 2019;White et al., 1998). When τ j = 0, there has been no loss of element j between the parent material and the regolith, whereas τ j = -1 illustrates total depletion compared to the parent material, and τ j > 0 indicates accumulation.
A numerical multicomponent RTM is developed using the open-source, isotope-enabled CrunchTope software (Druhan et al., 2013;Steefel et al., 2015). The global implicit scheme (Steefel & Lasaga, 1994) is used for this application, which allows simultaneous solution of transport and the multicomponent reaction network. For a given primary solute C i in mol/kg-H 2 O the governing equation is given as:
where ϕ is the porosity of the medium (m 3 void/m 3 porous medium), D i is the hydrodynamic dispersion coefficient (m 2 /s), and u is the Darcy flux (m 3 H 2 O/m 2 porous medium/s). The two summations in the equation correspond to kinetically controlled homogeneous (aqueous, subscript r) and heterogeneous (mineral, subscript m) reactions where v is the number of moles of the solute involved in a given reaction, and R r (mol/kg H 2 O/s) and R m (mol/m 3 /s) are the overall rates of a given aqueous and mineral reaction, respectively. The aqueous solution is composed of a set of components (primary species) that are coupled to a larger set of secondary aqueous species and gas phases through their respective equilibria (Table S1 in Supporting Information S1).
The mass-balance basis of the model means that reactions alter both the fluid solute chemistry and the mineral abundances through space and time. Mineral volumes are updated at the end of each time step (Lichtner, 1996;Steefel & Lasaga, 1994):
where ϕ m is the volume fraction of a mineral m (m 3 mineral/m 3 porous medium), V s is the mineral molar volume (m 3 mineral/mol mineral), R m is the mineral reaction rate (mol mineral/m 3 porous medium/s), and Δt is the size of the current time step (s). This means that the subscript t -Δt corresponds to the mineral volume from the previous time step. Initial mineral-specific reactive surface areas are updated as (Lichtner, 1996;Steefel & Lichtner, 1998):
where A m,o , ϕ o , and ϕ m,o are the respective starting surface area (m 2 /m 3 porous medium), porosity (m 3 void/m 3 porous medium), and mineral volume fraction (m 3 mineral/m 3 porous medium) specified as initial conditions and ϕ is the current porosity (Lichtner, 1985(Lichtner, , 1988;;Steefel & Lasaga, 1994;Steefel & Lichtner, 1998):
Uplift and erosion are represented in the model through a prescribed steady-state volumetric uplift rate. Solid phase from the bottom of the domain is advected upwards while material at the top of the domain is removed at the same rate, offering a simplified representation of erosion.
The model is constructed to produce the modern-day contemporaneous geochemical composition of both the solid weathering profile and dissolved solutes, based on a steady state balance between fluid infiltration and solid phase uplift. We accomplish this by constructing a representative 1D vertical column initially composed of fresh, unaltered granite subjected to constant rates of (a) meteoric fluid infiltration and (b) uplift/erosion. The domain is discretized into 200 0.1-m long grid cells. This 20-m length is twice as long as the inferred maximum vertical depth of the weathered regolith based on shallow seismic surveys (Section 2.1; Dangeard et al., 2019) and is designed to approximately represent a flow path length from vertical infiltration and drainage to lateral flow and discharge into the stream, ultimately supplying baseflow to Sapine Creek. Fluid flow is specified at a constant rate of 1.4 m/year, which is based on the local mean annual precipitation (2,000 mm; Section 2.1) minus the mean annual evapotranspiration (600 mm; Section 2.1). The rate of uplift and erosion (0.0001 m/year) is consistent with what has been reported for the surrounding region (Section 2.1; Masson et al., 2019;Malcles et al., 2020). From this starting point, the simulation requires 10 5 year timescales to reach steady state (Section 3.2), achieving a relatively simplified representation of watershed structure including average infiltration, drainage, uplift and erosion rates.
A Dirichlet upper boundary condition is prescribed with an infiltrating fluid solute (Figure S1 in Supporting Information S1) and isotopic composition (δ 7 Li = +8.6‰) based on rainwater samples collected at the site (discussed further in Section 4.3). A flux boundary condition is set at the base of the model domain. The initial condition of fluid within the domain is equilibrated with the granite following Maher et al. (2009) and Winnick et al. (2022) but this is quickly replaced with infiltrating meteoric water when timestepping begins. The initial condition of the solid phase is based on the elemental composition of the bedrock (Section 2.5.2) and the average measured δ 7 Li of individual porphyritic granite samples (-0.85‰). The simulations are run at the mean annual temperature of the field site, 7.0°C (Fernandez et al., 2022). Parameterization of the reaction network is summarized in Table 1 (fluid-mineral reactions) and S1 (aqueous equilibria) and details of development, including incorporation of Li stable isotopes, are outlined below.
The bulk elemental composition of the unaltered granite is represented in the model by an assemblage of primary minerals (quartz, K-feldspar, plagioclase; Text S1 and Figure S2 in Supporting Information S1) consistent with the granodiorite classification scheme for igneous rocks (Le Bas & Streckeisen, 1991). Biotite is used as a representative of the possible accessory minerals (Section 2.1) owing to its common occurrence and rich Li content (discussed further below). The volume fraction of each mineral is specified such that the distribution of elements in the model matches the bulk elemental composition of the mineral assemblage measured in our porphyritic granite samples (discussed further in Section 3.2).
Smectite, illite, and kaolinite are initially allowed to form as weathering byproducts. These minerals commonly occur in weathered granitic regolith profiles (e.g., Banfield, 1985;White et al., 2001). However, illite and smectite are ultimately omitted in the final simulations as they are not favored to form and remain thermodynamically undersaturated across the entirety of the steady state model domain. Dissolution of primary minerals composing the unaltered granite and precipitation of secondary clays occur simultaneously and are inherently coupled Note. Equilibrium and kinetic rate constants are for a temperature of 298 K. a Adjusted values to achieve best fit between observed and modeled solid-phase elemental depletion profiles. b Adjusted values to achieve best fit between observed and modeled bedrock elemental composition. c Palandri and Kharaka (2004). d EQ3/6 database (Wolery et al., 1990). e Arnorsson and Stefansson (1999). f Malmstroem et al. (1995). g x varies between 3E-5 and 1.455E-2 moles (see Section 3.3). h Initial bulk volume fraction for kaolinite (0.01%) partitioned into 7 Li and 6 Li endmembers. i Offset in rate constant to implement kinetic isotope fractionation (i.e., α = 7 k/ 6 k). (Maher et al., 2006). Each fluid-mineral reaction is described using transition state theory reversible kinetic rate laws (Lasaga et al., 1994)
where A m is the reactive mineral surface area (m 2 mineral/m 3 porous medium), k m is the kinetic rate constant (mol mineral/m 2 mineral/s), Q m is the ion activity product, which evolves both spatially and temporally through the simulations, and K eq is the temperature-dependent equilibrium constant. Parameter values for Equation 8 are given in Table 1 for each mineral in the model. Equilibrium constants are largely taken from the EQ3/6 thermodynamic database (Wolery et al., 1990) while kinetic rate constants are taken from the compilation of Palandri and Kharaka (2004). The equilibrium constant for oligoclase is sourced from Arnorsson and Stefansson (1999) and that for biotite is taken from Malmstroem et al. (1995).
A sequence of steps are taken in development and validation of the RTM. After running the simulation to steady state, the resulting model output is first challenged to reproduce the major element distribution we measured in the soil core relative to bedrock outcrop samples across the solid phase regolith profile. The rates produced for each mineral weathering reaction using Equation 8 are constrained to the specific elemental distribution of the Sapine weathering profile by reasonable adjustments to the starting reactive surface areas (discussed further in Section 3.2). Such adjustments are common to all RTM models applied to natural weathering profiles given the high degree of uncertainty associated with field-scale mineral surface area constraints. All values remain within the ranges suggested by prior studies for these mineral phases (e.g., Maher et al., 2006;Navarre-Sitchler & Brantley, 2007;White & Brantley, 2003).
We only introduce Li into the simulations after the major element chemistry accurately reproduces our observations. This sequence of steps first constrains the major element distributions across the system, and then introduces the trace phase, such that we are able to verify that incorporation of Li does not influence the overall bulk mineralogy or solute composition of the model. Rather, introducing Li into a functioning reaction network offers an additional means of testing the fidelity of the simulations by requiring accurate representation of a trace phase without the capacity to adjust any of the parameters associated with the major element distributions in the model.
The original stoichiometry of the primary minerals is updated to include Li consistent with the typical distributions of this trace element (Maneta & Baker, 2019). Lithium is not incorporated into quartz due to its relatively low solubility and resulting lack of reactivity in the simulations. Lithium substitution in the remaining primary phases is based on our own solid phase analyses (Text S1 in Supporting Information S1) and previously published data of mineral separates in a felsic lithology, such that correlation between a given element and Li was used to guide substitution (Table 1; Maneta & Baker, 2019). The total Li content of a given mineral is divided into 6 Li and 7
Li isotopologues and the relative amount of these "species" is set such that the δ 7 Li of the unaltered granite in the model is equivalent to the observed -0.85‰ for this bedrock.
Incorporating Li into the secondary mineral stoichiometry is inherently more difficult given that these phases are virtually absent in the initial mineral assemblage. Further, the partitioning of Li between the solid and the fluid from which the clay forms is dynamic and strongly dependent on temperature and solution chemistry (Barshad, 1957;Eberl et al., 1984). Synthesizing clays in laboratory settings remains an outstanding challenge, which impedes direct constraint of empirical relationships for the incorporation of Li into these phases. To date, only a few studies report relevant clay synthesis in a laboratory setting (Berger et al., 1988;Decarreau et al., 2012;Hindshaw et al., 2019;Vigier et al., 2008;Williams & Hervig, 2005), none of which include kaolinite and only two of which are conducted at low temperatures relevant to Earth surface conditions (Hindshaw et al., 2019;Vigier et al., 2008). Given these constraints, we rely upon batch experiments (W. Li & Liu, 2020) that expose fluid to pre-existing kaolinite to obtain an initial set of values to parameterize Li partitioning between these two phases (Text S3 and Figure S3 in Supporting Information S1). This approach is in agreement with our prior models (Golla et al., 2021(Golla et al., , 2022)). In the present reaction network, Li substitutes for Al in the stoichiometry of kaolinite as they share a similar ionic radius and prior evidence suggests this is a common behavior (Horstman, 1957;Odom et al., 1984;Starkey, 1982). In order to maintain charge balance, the stoichiometric content of Al substituted by Li is offset by the difference in their ionic charge (Table 1; Golla et al., 2021;Winnick et al., 2022).
The isotope composition of kaolinite is initially specified to match that of the unaltered granite bedrock (δ 7 Li = -0.85‰). Isotopic fractionation occurs as this clay forms from solution within the model domain. This mineral is treated as a solid solution with 7 Li and 6 Li endmembers that react in parallel (following Druhan et al. (2013)). A difference between the kinetic rate constants of the two endmembers reflects the appropriate fractionation factor for these reactions (i.e., α clay-diss = 7 k m / 6 k m ). Where clay (re-)dissolves, fractionation is omitted within the CrunchTope software (Golla et al., 2021;Winnick et al., 2022) as is consistent with expectation (Pistiner & Henderson, 2003;Wimpenny et al., 2010). The fractionation factor used in the model is initially based on a laboratory-derived value (α kaolinite-diss = 0.992; W. Li & Liu, 2020) and is later updated based on our data (Section 3.3). Clay formation is the single reaction pathway through which fractionation occurs in the model and is used to encapsulate what is likely a suite of fractionating mechanisms for Li. This is consistent with the modeling approaches used in Golla et al. (2021) and Winnick et al. (2022).
All geochemical measurements and accompanying data are provided in an open-access data repository (Golla, Kuessner, et al., 2024).
Mass transfer coefficients across the 80-cm soil core generally show τ < 0 for the major elements, indicating depletion relative to the granite protolith (Figure S4 in Supporting Information S1). The one exception is Ca, which maintains a τ > 0. This feature is considered further in Section 4.3. Across the major elements and Li, there are no clear monotonic trends with depth, and for the purpose of this study, we use depth-average values over the upper 80 cm for each element in further analyses (presented in Section 3.2). This averaging is justified given that the core is <5% of the 20-m length scale of the model domain. The average δ 7 Li of the soil (-2.98‰; 1σ = 0.62‰; n = 4) is significantly lighter than that of the porphyritic granite (-0.85‰; 1σ = 0.97‰; n = 5).
Stream water samples were collected over a wide range of discharge rates (0.5-90 L s -1 ) at coarse (approximately monthly) sampling frequency between 2013 and 2018. These values are taken as representative of average background conditions across multiple seasons, years, and flow rates. Solute and δ 7 Li values are not corrected for atmospheric contributions as they were in a prior geochemical study of this site (Fernandez et al., 2022). Rather, the use of a multi-component RTM allows us to explicitly account for the actual geochemical composition of the rainwater and embed this input into the upper boundary condition of the simulations. Importantly, the observed rainwater chemistry (dissolved Li concentration Li = 0.04 μM, Li/ Na = 1.2 μM/mM, and δ 7 Li = +8.3‰; Golla, Kuessner, et al., 2024) does not resemble a marine signature (Li/ Na = 5 × 10 -2 μM/mM and δ 7 Li = +31‰; Millero, 2016;Misra & Froelich, 2012) as is typically assumed for coastal sites such as our study area. The significance of this boundary condition is further explored in Section 4.3. We note that the model does not include direct dilution of solute concentrations in the stream by rainwater. This effect is further considered in Section 3.2, but we largely emphasize comparison between model and measured elemental and isotopic ratios to avoid such complications (e.g., Gaillardet, Dupré, Louvat, & Allègre, 1999;Meybeck, 1987;Stallard & Edmond, 1983). For this purpose, the solute concentrations of the stream are normalized to Na, which is generally not retained in the byproducts of silicate weathering (Sawhney, 1972). This Na-normalization eliminates the effects of dilution and/or evaporation which could otherwise impact absolute solute concentrations (e.g., Dessert et al., 2003;Gaillardet, Dupré, Louvat, & Allègre, 1999).
The major element chemistry of stream water at Sapine is consistent with the range of molar ratios observed in rivers draining primarily aluminosilicate lithologies (Figure 2; Gaillardet, Dupré, Louvat, & Allègre, 1999) and aligns with other small upland granitic watersheds reported by Oliva et al. (2003). This indicates that the porphyritic granite observed in outcrops at the site is reasonably representative of the bedrock composition. The relative abundance of lithium in solution (given as Li/Na) ranges from 1.3 to 2.4 μM/mM, while the corresponding δ 7 Li values vary from +21.7 to +24.2‰. This places the δ 7 Li -weathering intensity relationship of Sapine near the apex of the downward parabolic trend (Section 1) reported by Dellinger et al. (2015), and within the range of δ 7 Li measured in rivers draining larger catchments with floodplains displaying intermediate weathering intensity (Figure 2b). To date, we are aware of only one comparably small, upland watershed subject to similar climate which we can add to Figure 2b: Elder Creek in the Eel River watershed of northern California (Golla et al., 2021(Golla et al., , 2022)). In tandem, Sapine and Elder Creek reveal a characteristic class of small, wet and cool upland landscapes that sustain enriched riverine δ 7 Li at intermediate weathering intensity.
All RTM simulation results compared against the data from our site are first run to steady state, which is operationally defined as the point in time at which elemental mass influx and efflux are balanced. Total mass influx is the sum of infiltrating rainwater solutes and the elemental composition of uplifting granite at the base of the domain. Total mass efflux is the sum of solutes draining from the base of the domain and erosive removal of weathered solid phase at the surface. For the present model parameterization, this balance is achieved after approximately 200 kyr (Figure 3). This duration is reasonable assuming that the watershed has remained relatively undisturbed over a comparable timescale. The extent to which the Mont Lozère region surrounding Sapine was covered in ice during the Quaternary is still under debate. Indirect evidence from moraines in nearby valleys suggest the possibility of glaciation at altitudes >1,300 m (Etlicher & De Goër De Hervé, 1988;Veyret, 1978), which rules out the drainage area in Sapine (Figure 1). Therefore, we assume that there is minimal to negligible influence of potential glacial erosion and/or deposition over the 10 5 -year timescale simulated by the model.
The initial condition is representative of an unaltered granite, which is prescribed in the model as 30% quartz, 29% oligoclase, 23% K-feldspar, and 8% biotite with 10% porosity. This primary mineral assemblage is consistent with the elemental composition of the granite samples taken from our site to within ±5%. Furthermore, inserting the model bedrock mineral stoichiometries into Equation S1 (Text S1 in Supporting Information S1) returns elemental ratios that are in agreement with observations to within ±10% (Figure S5 in Supporting Information S1).
Once the model reaches steady state, the simulated solid-phase geochemistry of the upper boundary is constrained against our measured soil mass transfer coefficients, while the simulated solute composition at the base of the x" markers). The range of values for rivers draining a specific lithology (silicate, carbonate, and evaporite) are illustrated based on Gaillardet, Dupré, Louvat, and Allègre (1999). The molar ratios are calculated from major solute chemistry that is not corrected for atmospheric inputs. Panel (b) is the parabolic relationship between dissolved riverine δ 7 Li and weathering intensity originally presented by Dellinger et al. (2015), which now includes Sapine (this study) and Elder Creek (Golla et al., 2021). Weathering intensity is defined as the proportion of total landscape denudation (D is the sum of chemical weathering flux (W) and physical erosive flux (E)) that is attributable to W (Bouchez et al., 2014;Dellinger et al., 2015). For Sapine, the estimated value of E is used to prescribe steady-state uplift/denudation in the reactive transport model (0.0001 m/ yr or ∼270 t/km 2 /yr) and W is calculated from the long-term major ion chemistry data (Golla, Kuessner, et al., 2024). For the similarly small, mid-altitude Elder Creek catchment in northern California, USA, E is reported in Golla et al. (2021) and W is calculated from the long-term solute chemistry record afforded by the USGS Hydrologic network benchmark program (Clark et al., 2000). For both Sapine and Elder, W values were calculated from stream solute concentrations without correction for atmospheric inputs. At this scale, the uncertainty (taken as the greater value between the sample-specific and long-term analytical confidence interval) associated with the δ 7 Li of Sapine (±0.48‰) and Elder Creek (±0.76‰; Golla et al., 2021Golla et al., , 2022) ) is smaller than the size of the symbols shown in subpanel (b).
Journal of Geophysical Research: Earth Surface 10.1029/2023JF007359
domain is constrained against our measured stream solute concentrations. The modeled bulk solid phase is characterized by a spatial gradient from zero τ values at the bottom of the domain to elemental depletion at the soil-atmosphere interface (Figure 4). This compositional gradient is a consequence of simultaneous dissolution of primary minerals and formation of kaolinite (Figure 5), and the model accurately reproduces all major element τ values with the single exception of Ca. The corresponding steady state fluid major element concentrations exiting the base of the domain are slightly higher than the range measured in the stream, and the aqueous solutes of Sapine generally follow a dilution trend between the simulated effluent fluid and the rainwater boundary condition (Figure S7 in Supporting Information S1), suggesting a range of water flow paths through the watershed which are generally flux-weighted toward lengths shorter than the 20-m model domain. In total, we take these results as evidence of model fidelity in that we are simultaneously reproducing the majority of contemporaneous solid and fluid major element chemistry governed by the fluid-rock reactivity that underlies the development of this weathering profile over geologic time.
The model is next tasked with reproducing observed modern day Li distributions in Sapine by incorporating Li into the initial condition and spin-up to steady state. As described in Section 2.5.2, this sequential addition of Li to the constrained major element simulation is intentional, as the incorporation of a trace element phase should not alter the behavior reported in Section 3.2 and discussed further in Section 4.1.
Infiltrating fluid in the model is assigned an Li composition reflecting that of measured rainwater (Li = 4.77 × 10 -3 μM; Li/Na = 1.2 μM/mM, n = 17; δ 7 Li = +8.62‰, n = 1) and Li is added to the appropriate minerals composing fresh bedrock (i.e., K-feldspar, oligoclase, and biotite). As solutes concentrate in solution due to the dissolution of these primary phases, the clays that begin to form are allowed to incorporate Li. This incorporation can be approximated to first-order with the following expression
Figure 3. The reactive transport model is run from an initially unweathered granite to steady state, which is achieved when the total elemental mass influx (rainwater entering the top of the model domain plus bedrock uplifting from the base) and efflux (fluid draining from the bottom of the domain plus solid phase eroded at the soil surface) are balanced. These mass fluxes are calculated based on the major element concentrations (Na, K, Ca, Mg, and Si) of the modeled aqueous and solid phases. The schematic representing the model domain at the start and at steady state is modified from Winnick et al. (2022). The time-series illustrates approach of the reactive transport simulation to steady state at ∼200 kyrs. The influx is time-invariant as these are prescribed in the boundary conditions of the model. The vast majority of mass influx is derived from the uplifting bedrock. Mass efflux is also strongly dominated by solid-phase erosion. However, the contribution of fluid to this summation of effluxes is more significant than rainwater is to the summation of influxes due to drainage of geochemically evolved subsurface fluid.
where Li kaolinite is the unit mass of Li present in kaolinite, Li diss is the concentration of Li in solution, and K Li d is a coefficient that describes the degree of Li partitioning between kaolinite and solution. At present, the RTM software does not allow for a simultaneous treatment of a dynamic Li kaolinite stoichiometry and a dynamic δ 7 Li of this mineral. Following the approach of recent Li-isotope enabled RTMs (Golla et al., 2021(Golla et al., , 2022;;Winnick et al., 2022), we specify a Li content of kaolinite corresponding to a range of conditions supporting Li partitioning (Text S3 in Supporting Information S1). As a result, the quantity Li kaolinite is built into the stoichiometry of kaolinite in the RTM (Table 1) and therefore represents the moles of Li taken from solution per mole of clay formation (Winnick et al., 2022). To constrain this value, we require a K Li d estimate. As a starting point, we employ the experimental adsorption data set from W. Li and Liu (2020) to estimate the extent to which Li is taken up by kaolinite. However, as expected, these data show a relationship between the magnitude of Li partitioning and the amount of Li present in solution, which requires translation into a fixed parameter value for the RTM (Winnick et al., 2022). To arrive at this estimate, we utilize a supporting batch-reactor simulation that allows us to define the appropriate range of dissolved Li concentrations that are produced when meteoric fluid is allowed to interact with the mineral assemblage composing our porphyritic granite (Text S3 in Supporting Information S1). Using these concentrations, an appropriate partition coefficient is calculated from the empirical relationship derived from the data of W. Li and Liu (2020). The extent of Li isotopic fractionation during incorporation into kaolinite is similarly based on a laboratory-derived, intrinsic value of α kaolinite-diss = 0.992 (Section 2.5.2).
When we run our simulations using this literature-derived parameterization, the Li/Na of fluid exiting the base of the model domain is ∼10x higher than what is observed in the stream water of Sapine (Figure S8 in Supporting Information S1), while the corresponding simulated solid-phase Li/Na at the top of the domain is 4 μM/mM lower than what is seen in our soils. The δ 7 Li of the modeled water at the outlet (-0.52‰) is far too low compared to that measured in the stream (+23.2‰), while the modeled solid at the top of the domain (-0.86‰) is too enriched relative to the average of our measured soils (-2.98‰). The model produces very little evolution in the simulated water or bulk solid phase Li/Na and δ 7 Li from their respective boundary conditions (i.e., precipitation and bedrock) over the depth profile, all suggesting that the model is not producing sufficient incorporation of Li into the secondary clay in comparison to our field observations. These clear discrepancies indicate that the natural weathering profile developed at Sapine, which this model accurately simulates for the major elements (Section 3.2), must be strongly retaining Li, leading to diminished export by the stream relative to what the laboratory-based parameterization of Li-clay incorporation suggests. We take this as evidence that the combination of pathways through which solubilized Li is reincorporated into the weathering profile (i.e., secondary clay formation and sorption) are enhanced in this watershed relative to what has been previously constrained under isolated experimental conditions (e.g., Decarreau et al., 2012;W. Li & Liu, 2020) and in other field locations based on simpler, less process-based modeling approaches (e.g., Bagard et al., 2015;Dellinger et al., 2015).
The vast majority of parameters specified in our RTM cannot be altered without undermining the agreement between model output and the major ion geochemistry of Sapine (e.g., reaction rate constants, mineralogy; Figures 4 and5). Accordingly, the only remaining avenue is to propose an increase in the amount of Li taken up by kaolinite (from 0.03 to 14.5 mmoles Li/mol mineral). This is clearly a significant increase, but similar adjustments have been necessary in recent multi-component RTM simulations of lithium (e.g., Golla et al., 2021Golla et al., , 2022;;Winnick et al., 2022). With this adjustment, the modeled Li/Na values of water draining from the base of the 3) calculated from the average of measured soil compositions over the 80 cm core (colored dots) and simulated solid phase major element distributions across the model domain (colored dashed lines). As the model simulation directly tracks changes in mineral volume fraction (Equation 5), model outputs can be readily converted into τ estimations. The modeled Ca profile (which perfectly overlaps and is plotted on top of that of Na) is the one instance in which simulated results are inconsistent with measured soil values. This mismatch is likely associated with an exogenous dust input to the system, which is discussed further in Section 4.3. The shaded region represents model sensitivity derived from shifting all mineral surface areas by ±10% around the best-fit value, which has minimal effect on Li behavior (Figure S6 in Supporting Information S1). The goodness of fit of the model is quantified as the RMSE between the modeled solid at the top of the domain and the soil observations for Na, K, and Mg (7.5 × 10 -3 for the bestfit profile, n = 3; 2.38 × 10 -2 for the +10% profile, n = 3; 3.95 × 10 -2 for the -10% profile, n = 3). With this adjustment, simulated δ 7 Li ratios show a greater degree of fractionation due to the influence of clay formation (Figure S8 in Supporting Information S1). However, there is still a ∼18‰ difference between modeled fluid δ 7 Li and the stream and a ∼1.8‰ difference between the modeled solidphase at the top of the domain and measured soil. Once again, the model offers evidence for discrepancy between field-and laboratory-constrained fractionation factors, now illustrating that the Li isotope fractionation associated with clay formation must be larger in Sapine than prior literature-reported values (α kaolinite-diss = 0.992; W. Li & Liu, 2020), as there are no other parameters remaining which can be used to mediate this disparity without losing accurate representation of the major elements and lithium across the domain. Accordingly, the fractionation factor for kaolinite is adjusted from α kaolinite-diss = 0.992 to α kaolinite-diss = 0.972 ± 0.001, which leads to a greater degree of isotopic fractionation and agreement between the model and our data (Figure 6 and Figure S8 in Supporting Information S1). Such a large fractionation factor for kaolinite has been reported previously both in field settings (α kaolinite-diss = 0.976; L. Zhang et al., 1998) and in the laboratory (α kaolinite-diss = 0.976; X. Y. Zhang et al., 2021), although the latter study is most relevant to marine environments. In addition, this adjusted value is comparable with bulk fractionation factors previously used in other Li-based In subpanel (c), the uncertainty (taken as the greater value between the sample-specific and long-term analytical confidence interval) associated with the δ 7 Li of Sapine (±0.48‰) is incorporated into the shaded gray area representing the stream observations and is smaller than the size of the symbols representing the solid-phase observations. The shaded region surrounding the best-fit model profile in panels (a) and (b) represents the solution space generated from adjusting the degree of Li partitioning (14.5 ± 1 mmoles Li/mol mineral) to reproduce the maximum and minimum observed stream Li/Na shown as a gray shaded area in panel (a). Based on this Li/Na constraint, the corresponding model results for Li (μM) overlaps with the observed range of stream concentrations and reaches higher values only across the deepest sections of the representative flowpath. The solution space in panel (c) is similarly generated by adjusting the isotopic fractionation factor (α kaolinite-diss = 0.972 ± 0.001) to reproduce the maximum and minimum observed stream δ 7 Li. The goodness of fit of the model is quantified as the RMSE between the modeled solid at the top of the domain and the soil observations and between the modeled effluent and median value of stream observations (δ 7 Li = 1.7‰, n = 2; Li/Na = 2.1 μM/mM, n = 2). RTM studies (α clay-diss = 0.978, Golla et al., 2021Golla et al., , 2022; α clay-diss = 0.950-0.975, Wanner et al., 2017; α clay-diss = 0.982, Winnick et al., 2022).
With these adjustments, guided by the sequential construction of the model beginning with major elements then including Li and ultimately Li isotopes, we produce contemporaneous simulated solid and aqueous phase steady state spatial profiles of Li signatures across the model domain (Figure 6). Over the first ∼5 m, the fluid Li/Na values slightly increase from the initial rainwater signature (1.2 μM/mM) but show minimal change with further distance along the flow path (Figure 6a) due to the simultaneous dissolution of Li-bearing primary minerals and formation of kaolinite (Figure S9 in Supporting Information S1). Over the same spatial interval, the Li isotope ratios are strongly enriched, beginning from a value (+10.6‰) close to that of the imposed rainwater boundary condition (+8.6‰) and quickly reaching a maximum of +27‰ (Figure 6c). This behavior highlights the sensitivity of δ 7 Li to secondary mineral formation within the relatively steady Li/Na signal (Figure 5; Figure S9 in Supporting Information S1). With increasing distance, the rate of solubilization of primary minerals that contain the most Li (i.e., biotite and oligoclase; Figure S9 in Supporting Information S1) increases, leading to a more pronounced influence of dissolution to the water (i.e., imparting solid-phase Li/Na and δ 7 Li values). The result is a slight increase in the dissolved Li/Na ratio from 1.2 μM/mM at 5 m to 1.6 μM/mM at the outlet. Corresponding dissolved δ 7 Li values gradually decrease from +27‰ to +24‰ at the terminus of the flow path.
The capacity of the RTM to embed both water and solid phase transit times allows us to study the contemporaneous changes in mineralogy across the simulated flow path and attendant Li in the solid phase. Starting from bedrock Li/Na (14 μM/mM) and δ 7 Li (-0.85‰) values at the base of the domain, bulk solid Li/Na increases toward the surface, while corresponding δ 7 Li values slightly decrease (Figure 6), both in agreement with measured soil values. It is important to recognize that these bulk values reflect the mass-weighted sum of all mineral phases that contain lithium, each with their own δ 7 Li signature. These phase-specific δ 7 Li values remain -0.85‰ for the primary minerals (K-feldspar, oligoclase, and biotite), but the composition of kaolinite changes as a function of depth and evolves through time to steady state. This is particularly notable near the upper boundary condition where the accumulation of clay is sufficient to drive the bulk solid phase composition to values as low as -2.8‰ (Figure 6c).
In this study, an RTM approach is required to describe the dynamics of Li transport and reactivity across a small upland watershed. Here, the difference in major element and Li concentrations as well as Li/Na ratios between precipitation and stream water is extremely low, and yet there is a large increase in the δ 7 Li values in the creek relative to atmospheric sources (Figure 6). Such behavior could not be reproduced with simpler Rayleigh or steady state open flow through models. Similar limitations were demonstrated for Li behavior in natural systems by Bohlin and Bickle (2019) and Golla et al. (2022). At Sapine, as in these previous studies, a process-based and multi-component RTM is necessary to contextualize the internal reactivity that underlies the evolution of dissolved Li between infiltration and discharge, with a boundary condition strongly influenced by atmospheric deposition. Notably, the RTM used in this study expands upon previous applications by simulating the establishment of a steady state weathering profile over geologic timescales.
Using a sequential set of steps in our data-model comparison from major element chemistry to Li and ultimately δ 7 Li, we find it is necessary to facilitate greater Li uptake by clays and larger isotopic fractionation than what would be produced if we base our parameters on prior experimental data. Alternative adjustments to parameterization could achieve the same behavior in Li, but they are ruled out given that they would also impact the major elements. For example, changing the adjusted values for mineral surface areas by ±10% (Figure 4) significantly alters the bulk composition of the resulting weathering profile developed from the unaltered granite (Figure 5) and causes the model results to deviate from our observed mass transfer coefficients. Similarly, if one or more of these mass transfer coefficients were changed such that mineral surface area values(s) required revision, the corresponding behavior of Li and δ 7 Li would be strongly impacted. Hence, constraint of the Li system prior to or independent from the major element chemistry fails to contextualize the trace element in the overarching geochemical behavior of the system, and may lead to apparently reasonable modeled Li behavior that would yield correspondingly erroneous major element mass transfer coefficients. This demonstrates the Journal of Geophysical Research: Earth Surface 10.1029/2023JF007359 advancement offered by our stepwise multi-component RTM approach in which we initially constrain the overall behavior of the system such that Li dynamics are embedded in and determined by the geochemical reactivity of the weathering profile.
In our approach, we use the available empirical constraints as a starting point. From here, we employ a logical sequence of model development and observation-based constraint to ultimately achieve an Li partition coefficient and fractionation factor appropriate to these environmental conditions. The necessity to adjust these parameters in the model suggests behavior in natural systems that is not presently captured by available laboratory experiments.
Sensitivity of the values of these parameters to a variation of ±50% in both fluid drainage and erosion/uplift rates does not alter the magnitude of these adjustments (Table S2 in Supporting Information S1), reflecting a strong affinity for Li by the clays in this system attended by large mass-dependent isotopic fractionation (Text S4 in Supporting Information S1). Such strong uptake of Li by clays in our model is consistent with other reactive transport frameworks used to characterize Li behavior in field settings. For the Rivendell hillslope draining into Elder Creek in Northern California, an isotope-enabled RTM required increased illite formation rates relative to laboratory values in order to reproduce the extent of Li uptake and isotopic fractionation in fluids draining unsaturated bedrock (Golla et al., 2021). The Li-enabled RTM developed by Winnick et al. (2022) estimated an appropriate partition coefficient for Li using observations from the Amazon River basin (Dellinger et al., 2015), suggesting ∼90% of solubilized Li is incorporated into newly forming clays. In the Alaknanda river basin that serves as headwaters of the Ganges, Bohlin and Bickle (2019) Such elevated incorporation of Li into clays at the field scale also appears to be associated with relatively large isotopic fractionation factors (e.g., Golla et al., 2021;Lemarchand et al., 2010). These observations are seemingly inconsistent with the general expectation of muted fractionation as multiple flow paths are flux-averaged across heterogeneous systems (Druhan & Maher, 2017). Silicon isotopes measured in Sapine stream water appeared to exhibit this muting behavior, albeit over the course of a storm event and likely subject to the compounding effects of isotopic fractionation associated with biological cycling (Fernandez et al., 2022). Stream δ 7 Li during a storm event have been reported for Elder Creek in northern California, USA (Golla et al., 2022), yet the observed isotope ratios (+26 to +29‰) still reflect extensive isotopic fractionation. It is possible that muting of Li isotopic fractionation does occur at larger length scales (e.g., Bagard et al., 2015). Discerning such a scale-dependence on the muting effect would require further expansion of the data sets now available for riverine δ 7 Li records in nested catchments. Given our present observations, it appears there is emerging evidence of large Li fractionation factors in natural watersheds which have not been fully captured by laboratory environments. Such extensive enrichment of riverine δ 7 Li in natural systems appears to be inconsistent with contemporaneous δ 30 Si, but this observation is based on sparingly little data. Experiments that can account for multiple reaction pathways and the complexity associated with natural weathering reaction assemblages (e.g., Pogge von Strandmann et al., 2022) may provide a means of uncovering the factors driving these disparities.
Given the apparently significant role of secondary phases in setting the magnitude and composition of chemical weathering in the Sapine system, we now leverage the forward RTM framework to study this behavior in further detail. Kaolinite is the principal solid weathering byproduct formed in the present model at steady state, but in general, we often observe a sequence of secondary mineral formation across a compositional gradient (Goldich, 1938;Langmuir, 1997;Railsback, 2003). In the early stages of weathering, precipitation of clay minerals like smectite and illite are driven by partial dissolution of primary aluminosilicates such as feldspars and micas. However, these clays are only metastable under near-surface conditions due to relatively weak bonding of elements with low ionic potential, like K + and Mg 2+ , which are more readily removed from the crystalline structure. As these cations are leached by infiltrating dilute, weakly acidic fluid, such metastable phases are transformed into a suite of more mature weathering byproducts typified by kaolinite, oxides, and oxyhydroxides. The maturity of these phases is evident in the relative enrichment of lithophilic elements that less readily form soluble cations (e.g., Al 3+ , Fe 3+ ). Hence, the formation and predominance of kaolinite in the weathering profile of Sapine, which is consistent with observations in other granitic weathering profiles (e.g., Frings, Schubring, et al., 2021;Grant, 1962;Jeong, 2000;Kretzschmar et al., 1997;Sequeira Braga et al., 2002), suggests conditions that favor removal of more soluble cations from the system by high infiltration rates and large volumes of fluid drainage through the system. Clay-forming elements (Si and Al) are additionally extracted from solution by kaolinite formation. Given the inherent coupling between primary mineral dissolution and clay formation (L. Li, 2019;Maher & Navarre-Sitchler, 2019), removal of solutes by a combination of high infiltration rates and formation of mature secondary phases sustains undersaturation and promotes further dissolution of primary minerals. Critically, this behavior does not imply a high chemical weathering flux, rather, that high infiltration rates combined with moderate erosion rates such as those observed at Sapine (Martin et al., 2003) facilitates both maturation of secondary clays and moderate export of more soluble cations. Apparently, Li remains retained in the kaolinite that forms from this process, a behavior which deserves more targeted studies.
The RTM allows us to further appraise the stability of secondary mineral phases in Sapine through a sensitivity test of flow rate. Adjusting the flow rate by ±10% relative to the current value of 1.4 m/yr does not appreciably influence the approach to steady state or the time at which this point is reached (Figure S10 in Supporting Information S1). Furthermore, kaolinite accumulation in the shallow layers of the model is only barely altered. A flow rate of 1.56 m/yr produces a 0.1% increase in kaolinite volume fraction (Figure S11 in Supporting Information S1). This slight increase in accumulation is coincident with a deepening of the reaction front of K-feldspar, causing a greater degree of depletion from the initial unaltered granite across the domain. In comparison, a lower flow rate yields less accumulation of clay due to attenuation of the K-feldspar dissolution rate across the spatial profile. These relationships demonstrate the tight coupling between primary mineral dissolution and clay formation. Under much wetter conditions, kaolinite may no longer be stable in favor of more mature oxide or oxyhydroxide weathering byproducts (Brantley et al., 2023;Eberl et al., 1984). In a much drier climate, smectite and illite would likely predominate due to a lack of intermittent flushing and accumulation of soluble cations in the weathering profile (Deepthy & Balakrishnan, 2005;Moravec et al., 2020;Tardy et al., 1973). Nonetheless, kaolinite accumulation at Sapine is relatively unaffected by small (±10%) differences in the infiltration rate used to develop the model to steady state, suggesting fidelity in our modeling framework. Furthermore, the stability of kaolinite in this system is predicated on the capacity for weatherable primary minerals (e.g., K-feldspar) to reach the shallow soil layers, offering a means of protection against further resolubilization of the clay. This shielding of kaolinite is based on the greater solubility of these primary minerals that are favored to weather first (Figure S9 in Supporting Information S1). Therefore, the capacity of the clay formed in the Sapine weathering profile to strongly retain Li appears to be resilient under environmental conditions typified by this location, suggesting a potential prevalence of such behavior in similar cool, wet climates with low dissolved solute exports (Figure 7a).
Atmospheric deposition impacts the isotope composition of Li in weathering profiles and streams in environments drained by dilute waters (Chapela Lara et al., 2022;Clergue et al., 2015;Fries et al., 2019). The RTM we use to account for weathering processes at Sapine incorporates an upper boundary condition for water Li concentration of 0.04 μM, with Li/Na = 1.2 μM/mM, which is the average of 16 precipitation samples. These samples individually range in Li/Na ratios from as low as 0.29 μM/mM to as high as 5.2 μM/mM. The average values are essentially indistinguishable from a flux-weighted average based on accumulated precipitation in each sample, implying that there is no clear relationship between the intensity or duration of a precipitation event and the chemistry of the rainwater. Only one δ 7 Li measurement is available at this time for precipitation (+8.3‰, Figure 6), owing to the large sample volume necessary for analysis, but the Li/Na ratio of this sample is indistinguishable from the mean 1.2 μM/mM value used in the model.
Clearly, the geochemical composition of this precipitation is distinct from the common assumption of marine Li signatures in precipitation (Li/Na = 5 × 10 -2 μM/mM and δ 7 Li = +31‰; Millero, 2016;Misra & Froelich, 2012). The possibility of anthropogenic contributions can be ruled out given that the δ 7 Li signatures of fertilizers and soil amendment are very enriched in 7 Li (Millot et al., 2010;Négrel et al., 2010) and any such input would overwhelm the δ 7 Li signature of such a pristine environment (Cognard-Plancq et al., 2001;Durand et al., 1991). Rather, our rainwater data must reflect mixing between a marine signal and an input with relatively high Li/Na and low δ 7 Li ratio. Furthermore, the elevated Ca mass transfer coefficient (Figure 4 and Figure S4 in Supporting Information S1) is far higher than that produced by any reasonable modeled mineral assemblage. A plausible explanation for these observations is an exogenous dust input to the watershed, likely derived from sections of the Saharan desert (Goudie & Middleton, 2001;Israelevich et al., 2012), which commonly contains Journal of Geophysical Research: Earth Surface 10.1029/2023JF007359 carbonate (Krueger et al., 2004) and has been shown to influence weathering profiles across the European continent (e.g., Angelisi & Gaudichet, 1991;Castorina & Masi, 2015;Menéndez et al., 2007). Taking an average Li composition for Saharan dust (Li/Na = 65 μM/mM and δ 7 Li = +0.7‰; Clergue et al., 2015) and a typical marine aerosol composition, the precipitation signature measured at Sapine suggests the observed rainwater signal (Li/Na = 1.2 μM/mM) is composed of 96% dust input based on a mixing expression for Li/Na (Golla et al., 2022, their equations S3 and S4). A comparable calculation based on δ 7 Li yields a dust contribution of 74%, though again we caution the use of only one measurement. Based on our Li/Na data (n = 16), the rainwater composition is closely bounded between marine aerosol and dust end members and strongly influenced by the Figure 7. A comparison of Sapine streamflow signatures to global trends. Subpanel (a) compares the major element composition (as the sum of the concentrations of major dissolved cations, Ca 2+ , Mg 2+ , K + , and Na + , and Si) of Sapine to that of other granitic catchments compiled by Oliva et al. (2003). Subpanel (b) overlays Sapine and Elder Creek (Golla et al., 2021(Golla et al., , 2022) ) on the δ 7 Li-Li/Na trend classified by different geomorphic environments from F. Zhang et al. (2022). Subpanels (c) and (d) show the dissolved Si and Li fluxes of Sapine relative to those observed in other upland watersheds in Oliva et al. (2003) and the range of geomorphic environments featured in F. Zhang et al. (2022) and the similarly small, midaltitude Elder Creek catchment (Golla et al., 2021(Golla et al., , 2022)), respectively. The major solute chemistry used to calculate (a) the sum of cation and Si concentrations, (c) Si yield, and (d) Li yield are not corrected for atmospheric inputs. At this scale, the uncertainty (taken as the greater value between the sample-specific and long-term analytical confidence interval) associated with the δ 7 Li of Sapine (±0.48‰) and Elder Creek (±0.76‰; Golla et al., 2021Golla et al., , 2022) ) is smaller than the size of the symbols shown in subpanel (b).
Journal of Geophysical Research: Earth Surface 10.1029/2023JF007359 latter (Figure S12 in Supporting Information S1). This suggests the capacity for even lower δ 7 Li values in precipitation, as well as a potentially significant range of variability.
Such consideration of the origin(s) of atmospheric input(s) is necessary given that the generally low Li content of rainwater derived from marine aerosols can be easily overwhelmed by exogenous sources. This implies a potentially significant contribution of rock material to streams that is not sourced from the local watershed. RTM sensitivity testing, now based on the δ 7 Li composition of rainwater, is again used to evaluate the importance of this input to the development of Li signatures across the weathering profile (Figure S13 in Supporting Information S1). Changing the δ 7 Li of the upper boundary condition fluid by ±10‰ causes shifts in the δ 7 Li of both the fluid and solid phases down to 5 m into the domain (Figure S13 in Supporting Information S1). Below this depth, even such a large range of variation in precipitation δ 7 Li is attenuated, and the composition of the discharging fluid at the terminus of the 1D domain is unaffected. However, this δ 7 Li atmospheric input variability is influential over length scales <5 m, which can be thought of as reflecting regolith depths where shorter flowpaths to the stream are hosted. Therefore, any enhanced contribution of these shorter flow paths to the streamflow (e.g., during storm events) could make such variability in atmospheric δ 7 Li values influential on the composition of the Li export.
Variation in the δ 7 Li of rainwater is also propagated into the δ 7 Li of new solid phases forming in the shallow section of the weathering profile, reflecting isotopic fractionation in an open system. Parsing the clay signal (Figure S13 in Supporting Information S1) from the bulk solid composition (Figure 6c) clearly shows the dependence of the secondary mineral δ 7 Li on the composition of the infiltrating fluid. A more 7 Li-enriched infiltrating fluid produces a more positive δ 7 Li in the clays and vice-versa. This mechanism is consistent with observations from the Luquillo Critical Zone Observatory in Puerto Rico (Chapela Lara et al., 2022). In this tropical, intensely weathered site, extremely negative δ 7 Li values in the deep regolith are attributed to formation of clay from fluid that had attained a 7 Li-depleted signature as a result of solubilization of highly weathered material at shallower depths (Chapela Lara et al., 2022). Although the degree of weathering is much less advanced in Sapine (Figure 7a), inheritance of the infiltrating fluid signature by clays forming in the shallow subsurface demonstrates the capacity for soil δ 7 Li to reflect sensitivity to the source of atmospheric input(s) in upland environments described by intermediate weathering intensities.
In field-based studies, weathering-derived solute signatures of streams and rivers are commonly corrected for rainwater contributions (Dellinger et al., 2015;Gaillardet, Dupré, & Allègre, 1999;Négrel & Roy, 1998), often assuming a marine-like signature of rain in catchments proximal to or within coastal regions or in cases where there is no direct data. Though Sapine is only 80 km from the Mediterranean Sea, application of such an assumption to the stream chemistry would erroneously reduce the measured stream dissolved Li content by as much as 85%. Here, the strong influence of exogenous dust inputs negate such corrections even using the observed rainwater composition, which would suggest >100% of dissolved Li in the stream is derived from atmospheric sources. Clearly the Sapine stream δ 7 Li values rule this out. In such a situation, stream major solute compositions could be easily misinterpreted to reflect low or negligible weathering rates.
Our use of a multi-component RTM allows us to account for the role of a dust-influenced rainwater on a dilute effluent exiting the base of the domain through a coupled set of water-rock reactions in an open system, and sets the stage for further treatment of aeolian inputs in RTM frameworks. For example, at present we have no direct samples of dust deposition at the field site, much less a temporal record of these events sufficient to consider the timing of dry deposition events and their transmission into existing long-term stream chemistry. However, this work clearly points to the plausibility of both wet and dry deposition events as important contributions to chemical weathering fluxes from such pristine upland landscapes, and future studies should consider the relative timing and interplay of these inputs to upland systems. Further developments to this model allowing for direct aeolian deposition onto the soil surface could offer this capacity and improve model fidelity (e.g., Ca behavior; Figure 4). Here, we show that the typical marine-based rainwater correction used to approximate atmospheric input is neither appropriate nor sufficient to describe our data. This is likely to hold true across similar upland watersheds subject to exogenous aeolian deposition.
There is consistency between the low dissolved load of Sapine and that of other granitic catchments subject to similarly wet, cool climates (Oliva et al., 2003). Through the lens of major element chemistry, Sapine does not Journal of Geophysical Research: Earth Surface 10.1029/2023JF007359 diverge from expected patterns (Figures 7a and7c; Figure S14 in Supporting Information S1). In a similar manner, the streamflow Li signatures of Sapine fall within range of what has been observed from a diversity of geomorphic environments defined by tectonically stable lowlands and highly erosive, high mountain catchments (Figure 7b; F. Zhang et al., 2022). Yet the relatively low Li/Na ratios and high δ 7 Li values of the Sapine stream are more consistent with those of lowland sites that are typified by more developed weathering profiles and greater clay content. The only high mountain site that appears close to Sapine in this parameter space is the Taramakau River, which drains the western side of the New Zealand Alps. The rivers in this region are characterized by high dissolution rates (Lyons et al., 2005;Pogge von Strandmann & Henderson, 2015) and correspondingly high chemical erosion fluxes (Robinson et al., 2004). Such high dissolution rates may naturally encourage greater secondary mineral formation, and hence δ 7 Li enrichment, but ample supply of unweathered bedrock is maintained in these sites given that physical erosion overwhelms the total yield (Lyons et al., 2005). The result is a very different system from that of Sapine, capable of producing a high dissolved Li yield (Figure 7d). In contrast, strong retention of lithium in the solid phases forming at Sapine (Figure 6) limits dissolved Li export fluxes.
Both Li/Na and Si/Na stream water ratios at Sapine exhibit depleted values relative to that of fresh bedrock (Figure 6; Fernandez et al., 2022) reflective of retention of both elements in the weathering profile. However, stream δ 30 Si ratios reported for Sapine (+0.32‰, n = 17; Fernandez et al., 2022) are relatively low within the range of values reported for rivers (-0.2 to +4.7‰; Frings et al., 2016) and fairly close to the composition of the granitic bedrock (δ 30 Si = -0.23‰; Fernandez et al., 2022). In contrast, stream δ 7 Li ratios are high (+23.2‰; n = 15) and strongly indicate the preservation of secondary mineral formation signatures. This preservation occurs due to both (a) slow primary mineral dissolution over the whole weathering profile and (b) limited redissolution of secondary mineral weathering byproducts due to the persistence of primary minerals in the upper layers of the weathering profile, both of which would contribute a low δ 7 Li signal to the dissolved load resembling either bedrock or an even more 7 Li-depleted clay. Hence, for this system, the Li isotope tracer offers a more direct accounting of secondary mineral accumulation, whereas the δ 30 Si signal suggests a competition between the formation of secondary minerals and other fractionating pathways, such as the effect of ecosystem nutrient cycling suggested by Fernandez et al. (2022). Juxtaposing these two elements demonstrates a potential means of parsing the individual effects of secondary mineralization and plant uptake when coupling isotope systems with differing (bio)geochemical sensitivities (Frings, Oelze, et al., 2021;Pogge von Strandmann et al., 2022;Xu et al., 2021).
These comparisons also motivate further exploration of mid-altitude upland catchments defined by intermediate weathering intensity (Figure 2b). The Li stream signatures that describe Sapine are consistent with a similarly small (17 km 2 ), cool (11°C), wet (mean annual runoff of 1,744 mm) and steep (∼32°) catchment composed of shale lithology in the Northern California coastal ranges, referred to as Elder Creek (Kim et al., 2014(Kim et al., , 2017;;Rempe & Dietrich, 2018;Salve et al., 2012). As in Sapine, Elder Creek Li/Na ratios are relatively low with elevated δ 7 Li (Figure 7b). However, a similar RTM as well as direct sampling of vadose zone fluid draining to Elder Creek show that Li concentrations increase by ∼0.9 μM in 12 m despite a ∼2x faster fluid drainage rate (Golla et al., 2021). In comparison, Li accumulation in clays at Sapine is much greater and not solely restricted to the shallowest portions of the weathering profile, such that dissolved Li/Na ratios remain <12% that of fresh bedrock (Figure 6a). The result is that the Li yield for Elder Creek (210 mol/km 2 /yr) is almost twice that of Sapine (113 mol/km 2 /yr) even though the shale protolith of Elder has a lower Li content (61 ppm; Golla et al., 2021) than the granite of Sapine (91 ppm; Golla, Kuessner, et al., 2024). This more efficient export of Li in Elder Creek is due to greater, albeit incomplete, re-dissolution of secondary weathering byproducts in the shallow portions of the deep regolith as well as dissolution of clay phases composing the shale protolith (Golla et al., 2021(Golla et al., , 2022)).
The dissolved signatures that characterize Sapine and Elder Creek are a result of specific silicate weathering patterns and expand the range of environments that show δ 7 Li evidence of intermediate weathering intensities (Figure 2b). In these cool, wet upland sites, pronounced weathering incongruence coupled with a low total dissolved solute export suggest slow dissolution rates characteristic of kinetic limitation (Riebe et al., 2004;West et al., 2005). At Sapine, the flux of clay formation is moderate, but the clay that does accumulate is not significantly re-dissolved, such that the Li taken from solution is effectively "locked" into the solid phase over the development of the weathering profile. This behavior underscores the signatures of weathering incongruence through clay formation carried in a dilute stream and reveals the subtle nuances of geochemical sensitivity to shifts in climate and hydrologic forcing.
In a cool, wet upland catchment, we use a novel reactive transport approach to simulate the evolution of a granitic weathering profile to steady state based on fluid infiltration and bedrock uplift rates. We show that Sapine and other mid-altitude upland catchments appear to be expanding the range of environments that represent intermediate weathering intensity conditions (Figure 2b; Dellinger et al., 2015). In our model, the development of the weathering profile is directly connected to the solute signatures of the stream draining the catchment, which features enriched Li isotope signatures within low solute mass fluxes. In order to reconcile these observations, it is necessary to allow strong scavenging of Li from the dissolved load, effectively locking Li into the regolith. Such strong retention of Li is facilitated by the stability of clays promoted by the simultaneous maturation of weathering byproducts and the presence of more soluble primary phases even in the shallow subsurface. Altogether, this apparent shielding of clays (Section 4.2) and sequestration of Li therein (Section 4.1) illustrate a mechanism by which the dissolved weathering export from Sapine remains very low (Figure 7) despite substantial atmospheric deposition that we attribute to exogenous dust. The extent to which this behavior is common to such upland environments drained by dilute streams remains to be explored, but increasingly intense periods of drought and desertification associated with climate change are likely to exacerbate dust production and aeolian deposition rates (Avila & Peñuelas, 1999). Additionally, RTMs capable of explicitly treating such exogenous inputs could be used to identify the geochemical characteristics of soil paleorecords indicative of drier periods in Earth history.
21699011, 2024, 3, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JF007359 by Cochrane France, Wiley Online Library on [23/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License