Volcanic Arc Weathering Rates in the Humid Tropics Controlled by the Interplay Between Physical Erosion and Precipitation

Volcanic arcs are chemical weathering hotspots that may contribute disproportionately to global CO2 consumption through silicate weathering. Accurately modeling the impact of volcanic‐arc landscapes on the Earth's long‐term carbon cycle requires understanding how climate and physical erosion control weathering fluxes from arc landscapes. We evaluate these controls by examining the covariation of stream solutes, sediment geochemistry, and long‐term physical erosion fluxes inferred from cosmogenic 36Cl in magnetite in volcanic watersheds in Puerto Rico that span a ca. 15‐fold gradient in specific discharge. Analysis of this data using power‐law relationships demonstrates that CO2 consumption from arc‐rock weathering in the humid tropics is more strongly limited by physical erosion and the supply of primary minerals to the weathering zone than by temperature or the flux of fresh, chemically reactive waters through the critical zone. However, a positive correlation between long‐term physical erosion fluxes and specific discharge is also observed. This indicates that fresh mineral supply in arc environments may ultimately depend on precipitation rates, which may maintain a coupling between arc‐rock weathering fluxes and climate under principally supply limited weathering conditions.


Introduction
Volcanic arcs may play an important role in the evolution of the Earth's carbon cycle and climate.Weathering of volcanic-arc rocks (e.g., andesite or basalt) produces large fluxes of Ca and Mg to the ocean and may contribute disproportionately to global CO 2 consumption through silicate weathering coupled to marine carbonate precipitation.Consumption of CO 2 through arc-rock weathering may stabilize the Earth's climate through a negative feedback mechanism if weathering fluxes depend on temperature (Walker et al., 1981), either directly through the temperature-dependence of reaction kinetics or indirectly through the role of temperature in setting the vigor of the hydrologic cycle.In this scenario, increases (decreases) in global temperature with rising (falling) pCO 2 increase (decrease) weathering fluxes and CO 2 consumption, buffering atmospheric pCO 2 and global temperature.Conversely, arc weathering may destabilize Earth's climate if weathering fluxes are not closely linked to temperature.Recently, increases in arc-rock weathering have been hypothesized to have forced major periods of cooling through the Phanerozoic (Gernon et al., 2021;Macdonald et al., 2019;Swanson-Hysell and Macdonald, 2017), including the onset of modern ice-house conditions in the late-Cenozoic (Molnar and Cronin, 2015;Park et al., 2020).To assess whether volcanic arcs are a net stabilizing or destabilizing component of the Earth's climate system, a better understanding is needed of the relative sensitivity of arc-rock weathering fluxes to climate versus erosion.
Chemical weathering rates depend on temperature, runoff, and the supply of fresh minerals by physical erosion, and each of these factors may be rate-limiting (Brantley et al., 2023;Bufe et al., 2021;Gabet & Mudd, 2009;Maher & Chamberlain, 2014;West et al., 2005).These are termed "kinetic," "thermodynamic," and "supply" limits.Under a kinetic limit, weathering rates are restricted by the temperature-dependent rate constants of the underlying chemical reactions such that only an increase in temperature can increase weathering fluxes.At a thermodynamic limit, solute concentrations in stream waters approach chemical equilibrium and an increase in the flux of fresh water through the critical zone is needed to increase weathering fluxes.A thermodynamic limit is sometimes called an equilibrium, runoff, or fluid-transport limit (e.g., Brantley et al., 2023;Bufe et al., 2021).Under a supply limit, regolith is nearly completely depleted in reactive-mineral surfaces, such that only increases in the rate of the weathering front's advance into fresh rock can increase the weathering flux.Where regolith thickness is in steady state, the advance of the weathering front depends on the surface erosion rate.Thus, a supply limit is also termed an erosive-transport limit (e.g., Brantley et al., 2023;Stallard and Edmond, 1983).
Evidence for kinetic or thermodynamic limits to andesite and basalt weathering comes from studies of riverine solute export that show that dissolved fluxes increase exponentially with mean annual temperature (MAT) (Dessert et al., 2003;Li et al., 2016) and that dissolved load decreases by less than expected from simple dilution across runoff gradients (Bluth & Kump, 1994;Gaillardet et al., 2011).The strong control of temperature and runoff proposed for weathering of mafic rocks contrasts with wet, felsic landscapes, where weathering fluxes are thought to be set primarily by the supply of fresh minerals to the weathering zone through erosion, even in rapidly eroding mountainous settings, and thus independent of climate (Larsen et al., 2014;Riebe et al., 2001).Support for supply limited weathering in these environments comes from the observation that weathering intensity (the fraction of denudation that occurs through chemical weathering) inferred from the chemical depletion of regolith is nearly invariant across large gradients in long-term erosion rate (Ferrier et al., 2016;Riebe et al., 2001).However, the control of physical erosion on weathering fluxes has not been systematically investigated for volcanic arc rocks, despite their more rapid weathering kinetics, which suggest that a supply limit might be reached over shorter erosional timescales.
Here, we examine the limits to volcanic arc-rock weathering in mountains in the humid tropics, where temperatures are high, precipitation rates exceed potential evapotranspiration, and where weathering fluxes are greatest (Dupre et al., 2003), by characterizing physical erosion fluxes, weathering fluxes, and weathering intensity in 27 watersheds in Puerto Rico.These are underlain by an exhumed Cretaceous to Paleogene volcanic arc sequence and span a ca.15-fold gradient in specific discharge but a <5°C range in MAT (Table 1).This enables us to examine the relative importance of supply versus thermodynamic limits, while largely controlling for temperature, and to work in actively eroding watersheds in an ancient arc setting without the complicating factors of hydrothermal activity, volcanic acidity, and deep flow paths that are often present in young volcanic landscapes (Gaillardet et al., 2011;Li et al., 2016).
We determined long-term, watershed-averaged erosion fluxes using in situ-produced 36 Cl in magnetite (Moore and Granger, 2019a), a weathering-resistant mineral that persists in regolith under high weathering intensity.This method captures millennial-scale erosion fluxes that are unbiased by land-use, which may perturb modern sediment fluxes in the tropics (Brown et al., 1998;Hewawasam et al., 2003).We infer weathering fluxes using stream solute concentrations and estimates of specific discharge (e.g., Dessert et al., 2003;Gaillardet et al., 1999Gaillardet et al., , 2011;;Stallard, 2012a) and weathering intensity from stream sediment geochemistry.We use this data to parameterize power-law relationships and interpret the degree of supply versus thermodynamic or kinetic limitation to arc-rock weathering fluxes from the power-law slopes.

Study Area
We studied small (1-18 km 2 ) montane catchments in northeastern Puerto Rico.These are divided into an eastern group around the Sierra Luquillo (n = 13, series A), and a western group in the foothills of the Cordillera Central (n = 14, series B) (Figure 1).Both areas have a tropical rainforest climate (Köppen-Geiger class Af) and similar MAT (21-25°C) (Daly et al., 2003).Mean annual precipitation (MAP) increases from the western group (ca.140-200 cm yr 1 ) to the highest elevation watersheds in the eastern group (>400 cm yr 1 ), where high precipitation rates are driven by adiabatic cooling of wet trade wind air rising over the Sierra Luquillo (Murphy and Stallard, 2012).Long-term records of catchment-scale water, solute, and sediment fluxes are available in several study watersheds in the Sierra Luquillo, which provides a unique opportunity for cross-comparison of results (McDowell et al., 2021;Stallard, 2012a;Wymore et al., 2017).
The study watersheds are underlain by late Cretaceous to Paleogene volcanic arc rocks with a compositional mode of basaltic andesite (Jolly et al., 1998;Krushensky, 2000;Moore et al., 2023).The original volcanic arc topography of Puerto Rico was eroded to base level in the middle Cenozoic.The island was uplifted to near present elevations in the Pliocene (ten Brink, 2005;van Gestel et al., 1998) and uplift since the late Pleistocene has been minimal, as evidenced by reef terraces at elevations of less than 10 m above modern sea level from marine isotope stage 5e (ca.125 ka), when eustatic sea level was 5.5-9 m higher than present (Dutton & Lambeck, 2012;Hubbard et al., 2008).The present geomorphology of Puerto Rico reflects ongoing fluvial adjustment to post-Pliocene base level and features both graded streams and rivers with prominent knickzones.Low relief landscapes above knickzones in the Cordillera Central were identified as relict surfaces by early workers (Lobeck, 1922;Meyerhoff, 1927) and this interpretation has since been extended to the Sierra Luquillo (Brocard et al., 2016).A transition from supply limited to kinetically or thermodynamically limited weathering with increasing erosion rate below knickpoints on quartz diorite in the Sierra Luquillo was inferred from base cation concentrations in shallow soils (Porder et al., 2015), however this relationship has yet to be examined on volcanic rocks elsewhere in Puerto Rico.
Study watersheds in the Cordillera Central were selected from both low gradient surfaces above knickzones and from high gradient streams coupled to modern base level to maximize the range of erosion rates captured, whereas watersheds in the Sierra Luquillo were chosen primarily for accessibility.Many study watersheds in the Cordillera Central are developed on epiclastic basalt that contains minor limestone (Jolly et al., 1998;Krushensky, 2000).Carbonates are less prevalent in the Sierra Luquillo.Soil maps indicate that the study watersheds are mantled by deep, clay-rich weathering profiles, which are classified predominantly as oxisols and ultisols (Moore et al., 2023;Staff, NRCS, 2023).Oxisols are more abundant in the Sierra Luquillo than in the Cordillera Central.

Assessing Weathering Limits Using Power Laws
Relationships between weathering, erosion, and runoff are conventionally characterized using power laws (e.g., Ferrier et al., 2016;Godsey et al., 2009;West et al., 2005).In watersheds with steady-state regolith thickness, the rate of advance of the weathering front is coupled to the surface erosion rate and the relationship between the longterm weathering and physical erosion fluxes (mass area 1 time 1 ) may be described by a power-law of the form: where W is the weathering flux, which in a watertight catchment equals the product of specific discharge (Q) (length time 1 ) and the average concentration of silicate-weathering derived solutes in stream water (C) (mass volume 1 ), a 1 is a scaling coefficient, E is the physical erosion flux, and b 1 is the power-law slope.A similar power-law relationship between weathering flux and specific discharge is given by: The power-law slopes in Equations 1 and 2 (b 1 and b 2 ) may be interpreted as reflecting the type and degree of ratelimitation.Under a strict supply or thermodynamic limit, weathering fluxes increase linearly with erosion or specific discharge, respectively (i.e., b 1 and b 2 equal unity).Conversely, when weathering fluxes are independent of erosion or specific discharge, the power-law slopes are equal to zero.Intermediate states are represented by 0 < b 1,2 < 1 (Figure S.1 in Supporting Information S1).
Equations 1 and 2 describe the weathering-erosion and weathering-runoff relationships separately.Yet, erosion and runoff may vary simultaneously.To analyze erosion and discharge in tandem it is convenient to divide the weathering and erosion fluxes in Equation 1 by the specific discharge and parameterize C as a power-law function of E/Q:

AGU Advances
10.1029/2023AV001066 where E/Q (mass volume 1 ) represents the long-term average concentration of sediment in stream water leaving a catchment.Equation 3 allows the relationship between weathering, erosion, and runoff to be visualized in two dimensions.Under a thermodynamic limit, solute concentrations in stream water approach saturation and are invariant with the sediment concentration (i.e., b 3 = 0), whereas under a strict supply limit, solute concentrations increase linearly with the sediment concentration (i.e., b 3 = 1).Intermediate states are represented by 0 < b 3 < 1 (Figure S.1 in Supporting Information S1).If weathering is principally kinetically limited, E/Q should be a poor predictor of weathering fluxes and model residuals (i.e., the difference between the observed and modeled values) should scale with temperature.
Where regolith thickness is in steady state, the relationship between weathering and erosion fluxes is reflected in the degree of depletion of the regolith in soluble elements.Ferrier et al. (2016) proposed a method for testing for supply limitation using a power-law relationship between weathering intensity (quantified from regolith geochemistry) and the long-term denudation rate from cosmogenic nuclides.Here, we modify this approach to compare ε, a metric derived from weathering intensity, with E/Q.We define ε as the constant of proportionality between the weathering and erosion fluxes: 4.7 ± 0.2 0.70 ± 0.05 9.1 ± 0.8 34.9 ± 1.1 24.9 ± 2.1 95.9 ± 3.1 0.048 ± 0.020 0.14 ± 0.14 129 ± 11 47 ± Figure 1.Precipitation, topography, and study watershed locations within northeastern Puerto Rico.Blue-yellow colors show mean annual precipitation (cm yr 1 ) from the PRISM model (Daly et al., 2003).Precipitation is greatest in the northeast at high elevations in the Sierra Luquillo, which displays a strong orographic precipitation gradient, and lowest in the southwest.Watersheds are colored in shades of orange by erosion flux determined from 36 Cl in magnetite.

AGU Advances
10.1029/2023AV001066 where ε may vary between 0 and ∞.This term is related to, but distinct from, weathering intensity, which is the fraction of denudation that occurs through chemical weathering (i.e., W W+E ) and ranges between 0 and 1.In steadystate, ε may be derived from a mass-balance of weathering immobile elements (Section 3.5).Equation 3 can then be recast using Equation 4as: Dividing both sides of Equation 5 by E/Q gives a relationship between ε and E/Q that is analogous in form to Equation 3: In Equation 6, the power-law slope (b 3 -1) will equal 0 when weathering is completely supply limited (i.e., ε will be invariant with E/Q) because depletion of regolith in soluble elements progresses as far as possible.Conversely, when weathering is strictly thermodynamically limited ε will vary inversely with E/Q (i.e., b 3 -1 = 1).If weathering is kinetically limited, ε will vary inversely with E, which controls regolith residence time, but be invariant with Q.Thus, ε will be negatively correlated with E/Q only to the extent that variation in Q is small, else E/Q will be a poor predictor of ε and model residuals will scale with both erosion flux and temperature.
We use Equation 3 to quantify the limits to arc-rock weathering by fitting the power-law slope using solute concentration, erosion, and discharge data from the study catchments.We examine the relationship for both TDS sil (here defined as the sum of bedrock-weathering derived dissolved Na 2 O, K 2 O, CaO, MgO, and SiO 2 in stream water), as well as dissolved CaO and MgO only.These two elements are especially interesting because their supply to the ocean largely controls long-term CO 2 sequestration potential through marine carbonate precipitation (Berner et al., 1983).We express solute concentrations as oxides because this allows for direct comparison with physical erosion of oxide minerals (Stallard, 2012a).We then cross-check these results against weathering intensity metrics from sediment geochemistry using Equation 6.All geochemical data necessary for parameterizing Equations 3 and 6 are determined from measurements of stream sediment and water samples, while discharge (Q) is inferred from regional rainfall-runoff relationships.

Sampling
Water and sediment samples were collected between November 2020 and January 2022.Water from each stream was sampled in the wet (June-November) and dry (December-May) seasons to assess intra-annual variability in solute loads.Water samples were field-filtered through 0.2-μm polytetrafluorethylene (PTFE) syringe filters into HNO 3 -washed low-density polyethylene (LDPE) bottles and acidified to ∼pH 1 with trace metal grade HNO 3 .Instantaneous discharge was determined by measuring the channel depth at 10 evenly spaced intervals and determining flow velocity at each interval using a FLOWATCH™ flowmeter.Bulk sediment was collected from stream beds, channel bars, or channel marginal deposits during the dry season and captured a wide distribution of grain sizes, ranging from clay to cobbles.To obtain magnetite for determining erosion fluxes, magnetic sediment was collected from stream beds using magnetic rods (Moore and Granger, 2019b).

Specific Discharge
To estimate long-term mean annual specific discharge, linear rainfall-runoff relationships were constructed from regional stream gauge records and catchment-averaged MAP: where Q is specific discharge (cm yr 1 ), m is the slope, P is MAP (cm yr 1 ), and b is an intercept (cm yr 1 ).To fit the slope and intercept, 61 USGS stream gauges in Puerto Rico with complete, multidecadal runoff data were used (U.S. Geological Survey, 2022).These exclude stations downstream of reservoirs and in the karst region of northern Puerto Rico to avoid biases from water diversion and subsurface drainage.Record lengths range from 18 to 68 years (Moore et al., 2023).Catchment-averaged MAP was extracted from a 1963-1995 normals raster for Puerto Rico interpolated from rain-gauge data using the PRISM model (Daly et al., 2003).For watersheds in the rain-shadow of the Sierra Luquillo, MAP was estimated using polynomial elevation-precipitation equations from a reanalysis of regional rain gauge data (Murphy et al., 2017).Catchment areas were delineated in QGIS (version 3.16.9)using a 1/3rd arcsecond DEM from the USGS National Map.The slope and intercept in Equation 7were fit using ordinary least-squares regression.Separate relationships were constructed for the watersheds in the Sierra Luquillo region and those elsewhere in Puerto Rico to accommodate differences in hydrology; the Sierra Luquillo is more densely forested than the rest of the island and displays strong elevation gradients in evapotranspiration because of increasingly persistent cloudiness with increasing elevation (Murphy and Stallard, 2012) (Figure 2).

Solute Concentrations
Dissolved SiO 2 , Ca 2+ , Mg 2+ , K + , Na + , PO 4 3 , SO 4 2 , and Cl in the acidified water samples were measured using a Horiba Ultima-Expert™ inductively coupled-plasma optical emission spectrometer (ICP-OES) equipped with extended ultra-violet capabilities enabling measurement of Cl.Concentrations were determined using a linear calibration curve constructed from NIST traceable standard solutions.Measured P and S intensities are assumed to be entirely attributable to the PO 4 3 and SO 4 2 species.For the major cations, two emission lines were measured to cross-check concentrations and identify interferences.Instrument drift was minimized through frequent re-standardization and reproducibility was monitored through repeated measurements of an in-house reference material with a similar Cl concentration and solute load to Puerto Rican stream waters.The standard deviation of repeat measurements was ca.10% for Cl and less than 5% for all other ions.Bicarbonate was estimated from the charge balance.
Interpreting silicate weathering fluxes from solute loads requires correcting for atmospheric inputs and dissolution of non-silicate minerals (e.g., carbonates) (Gaillardet et al., 1999).Major ions were corrected for atmospheric deposition using measured Cl concentrations and sea salt ion/Cl ratios from Stallard (2012b): The lower intercept of the Sierra Luquillo regression may reflect higher evapotranspiration in this densely forested region than on average across the island and the higher slope a decline in evapotranspiration with increasing elevation as cloud cover becomes increasingly persistent (Murphy and Stallard, 2012).
where [X*] water is the concentration of ion X in water derived from rock dissolution, [X] water is the ion concentration measured in the water, [Cl] water is the concentration of Cl measured in water, and X Cl ) rain is the ratio of ion X to Cl in precipitation.Next, dissolved Ca 2+ and Mg 2+ concentrations were corrected for carbonate dissolution using element ratios in bedrock and water: where [X * sil ] water is the concentration of Ca 2+ or Mg 2+ derived from silicate weathering, Na X ) * water is the ratio of cyclic-salt corrected Na + to Ca 2+ or Mg 2+ in water, and Na X ) rock is the Na to Ca or Mg ratio in the silicate component of the bedrock.Silicate Na/Ca and Na/Mg ratios were determined from geochemical analyses of bedload clasts devoid of carbonate minerals (Section 3.5).This approach assumes that the carbonate endmember has Na/Ca and Na/Mg ratios near zero.

Weathering Intensity Metrics From Sediment Geochemistry
The chemical depletion factor (CDF) is a quantitative metric for weathering intensity that is calculated from geochemical measurements of bedrock and regolith (Riebe et al., 2001(Riebe et al., , 2003)).In steady state, the flux of rock into the weathering zone within a catchment must equal the flux of solutes and sediment out of the catchment.In this case, a watershed-scale CDF can be constructed from the geochemistry of bedrock and the long-term average sediment exported from the catchment.Here, we take the suspended sediment (defined as the <63 μm fraction of the bulk sediment sample) composition to represent the average sediment.This is a reasonable approximation because suspended load constitutes the majority of sediment export from volcaniclastic watersheds in Puerto Rico (Larsen, 2012).Similarly, we estimate watershed-averaged bedrock composition from amalgamated bedload clasts.To examine weathering intensity for each element measured in dissolved load, element specific CDFs (CDF x ) are used: where [X] sed is the concentration of element x in suspended sediment, [X ] rock is the concentration in bedrock, [I ] rock is the concentration of an insoluble index element in rock, and [I ] sed is the concentration in suspended sediment.A total CDF (CDF total ), integrating over all elements measured in dissolved load, can then be computed from the sum of element specific CDFs, weighted by each element's concentration in bedrock: An index that relates the erosion and weathering flux, ε (Equation 4), can be derived from CDF X .Consider that the weathering flux for element X (W x ) is equal to the product of the weathering intensity, represented by CDF X , the concentration of element X, and the regolith production rate (D) (mass area 1 time 1 ): In steady state, the supply of an insoluble element to the regolith must equal the export rate through physical erosion: Dividing Equation 13through by [I] rock , substituting this result and Equation 10 into Equation 12, and rearranging produces Equation 14: Thus, ε for element x (ε x ) is given by: Similar to Equation 11, an ε total may be defined as: Here, ε total is calculated by summing over all major rock-weathering derived elements measured in stream water (Na 2 O, K 2 O, CaO, MgO, and SiO 2 ) and ε CaO,MgO from CaO and MgO only.
To determine CDF and ε, suspended sediment and bedrock clast geochemistry were measured by ICP-OES.First, the homogenized bulk sediment samples were dried and split into >4 mm and <63 μm fractions.A bedrock fraction was selected by hand-picking the >4 mm fraction to remove saprock and crushed in a ring-and-puck mill.Suspended load and bedrock samples were fused in an equal mixture of lithium metaborate/lithium tetraborate in Pt crucibles at 1000°C and dissolved in 10% HNO 3 and sufficient HF to ensure complete solubilization of silica.The solutions were analyzed for major elements (Si, Ca, Mg, Na, K, Ti, Fe, and Al) and Zr.Accuracy and precision were monitored by repeat measurements of legacy USGS reference basalt (BCR-1) and andesite (AGV-1).Error was estimated from the standard deviation of the repeat measurements of the reference material.

Analytical Methods
Magnetite was separated from the magnetic sediment samples and prepared for 36 Cl measurement using methods modified from Moore andGranger (2019a, 2019b).The magnetic sediment was crushed in a ring-and-puck mill and non-magnetic fines were washed off with deionized water while using a magnet to retain the magnetic fraction.Magnetite was further concentrated through repeated wet-grinding on a wrist-action shaker with zirconia grinding media followed by magnetic separation.Secondary Fe-oxides were extracted in dithionite-citrate-bicarbonate and grain surfaces cleaned of adsorbed Cl using HNO 3 and NH 4 (OH) (Moore andGranger, 2019a, 2019b).Samples were spiked with a 35 Cl-enriched carrier and digested in oxalic acid.After dissolution, Fe and Ti precipitates were removed by centrifugation and Cl precipitated as AgCl.The precipitate was recovered by centrifugation, dissolved in NH 4 (OH), and then re-precipitated with HNO 3 to suppress the 36 S isobar.The product was washed with deionized water, dried, and loaded into AgBr filled Cu cathodes for 36 Cl/Cl and 35 Cl/ 37 Cl measurement by accelerator mass spectrometry (AMS) at PRIME Lab, Purdue University.Measurements were normalized to the standard described by Sharma et al. (1990).Magnetite target chemistry was measured by ICP-OES on aliquots dissolved separately in HF/HNO 3 .

Calculating Erosion Fluxes
We interpret 36 Cl concentrations as reflecting physical erosion fluxes.The study watersheds are mantled by thick (up to ca. 15+ meters on ridge crests), highly weathered regolith (Buss et al., 2017;Murphy et al., 2012).In this setting, chemical mass loss likely occurs primarily at a weathering front below the penetration depth of most cosmic radiation (Bazilevskaya et al., 2013) and the accumulation of cosmogenic nuclides in a weatheringresistant mineral, such as magnetite, is sensitive principally to physical erosion (Campbell et al., 2022;Riebe and Granger, 2013).Erosion fluxes were calculated from measured 36 Cl concentrations and target chemistries following an approach modified from Marrero, Phillips, Borchers, et al. ( 2016), Marrero, Phillips, Caffee, & Gosse (2016) and Moore and Granger (2019a).Error was propagated using Monte Carlo methods to capture covariance between stable Cl and 36 Cl.A detailed description of the erosion calculations is given in the Supporting Information S1 and the codes used to perform the calculations are provided in the data repository (Moore et al., 2023).

Regression and Statistical Methods
The data in this study are subject to both analytical uncertainties, which are reasonably well-characterized, as well as unconstrained geological uncertainties from, for example, incomplete sediment mixing, dilution of solute concentrations at high discharge, bedrock compositional variability, or differences in averaging timescales between data types.Importantly, there are large uncertainties in the specific discharge estimates from the rainfallrunoff relationships that may be inadequately characterized by the prediction intervals of the regressions.Therefore, specific discharge error is also treated as unconstrained geological error.Both analytical and geological error must be considered when fitting power laws.
Power-law slopes are typically fit using linear regression on log-transformed data.Equations 1-3, and 6 require a regression method that can accommodate analytical and geological error in the x and y directions.We adopt a fitting algorithm based on the York regression, which allows individual data points to be weighted in both x and y (York et al., 2004).First, we calculate a York regression and the bivariate mean-square of weighted deviates (MSWD) of the fit (Mahon, 1996) weighting by analytical errors only.A MSWD greater than unity indicates scatter beyond analytical error because of geological uncertainties.In the absence of other information, it is reasonable to assume that these uncertainties are equal for all observations.Thus, to account for geological error, we iteratively add variance symmetrically to each log-transformed x and y observation and recalculate the regression until the MSWD reduces to unity.In effect, this algorithm estimates the geological variance from the excess scatter in the data.The final variance of each observation is equal to the sum of the analytical and geological variances.If analytical errors are negligible, this approach produces regression parameters that are similar to those from the well-known reduced major axis method.The quality of the fit is quantified by the uncertainties on the slope and intercept.

Specific Discharge and Solute Concentrations
There are strong spatial patterns in annual specific discharge and cyclic salt and carbonate corrected solute concentrations across the study area.Specific discharge estimated from the rainfall-runoff relationships ranges from 23 cm yr 1 to 355 cm yr 1 and increases from the Cordillera Central to the Sierra Luquillo.Cyclic salt corrected wet and dry season solute concentrations are within 2σ error of each other in 24 of 27 watersheds, despite a ca.40% difference in instantaneous discharge measured in the field, on average, between the two samplings.This indicates a high degree of intra-catchment chemostasis (Godsey et al., 2009), wherein solute concentrations do not vary with streamflow, and allows long-term average solute concentrations, to a first order, to be approximated from the mean of the wet and dry season samples without weighting by the annual hydrograph.
Measured solute concentrations and ratios are consistent with mixing between atmospheric, carbonate, and silicate sources.The dominant anions are HCO 3 and Cl , whereas cationic charge is divided sub-equally between Ca 2+ , Mg 2+ , and the alkali elements (Figure S.2 in Supporting Information S1).A plot of Na-normalized Ca and Mg concentrations shows that many water samples are offset from bedrock ratios, indicating that carbonate dissolution contributes significantly to the solute load (Figure S.3 in Supporting Information S1).Mean annual total dissolved silicate solids corrected for atmospheric inputs and carbonate dissolution and expressed as oxides (TDS sil ) range from 21.8 μg g 1 to 100 μg g 1 and increase with decreasing specific discharge.The sum of corrected CaO and MgO concentrations averages 33 ± 9% of TDS sil (mean and standard deviation) (Moore et al., 2023).

Sediment Geochemistry and Weathering Indices
Suspended sediment and bedload clasts are compositionally distinct.Insoluble element concentrations in bedload clasts are consistent with expected values for bedrock.For example, Zr and TiO 2 concentrations in clasts from A1 (74 ± 5 μg g 1 and 0.53 ± 0.03 wt.%) are broadly comparable to the concentrations in bedrock reported for a subcatchment of this watershed (82 ± 4 μg g 1 and 0.65 ± 0.02 wt.%) (Buss et al., 2017).Insoluble element concentrations in suspended sediment are generally higher than in the bedload clasts and concentrations of base cations are uniformly lower (Moore et al., 2023).
Calculating weathering indices requires choosing an index element that is immobile during chemical weathering.Two relatively immobile elements, Ti and Zr, were measured, but Ti might be preferentially leached relative to Zr in weathering profiles in Puerto Rico (Buss et al., 2017).Weathering indices do not show strong spatial patterns across the study area.Values of CDF CaO and CDF MgO are higher on average than CDF total (with a mean and standard deviation of 0.66 ± 0.17 and 0.51 ± 0.22 vs. 0.25 ± 0.11).This indicates less complete loss of Si than of base cations from regolith.In 26 of 27 watersheds, ε total is less than unity, implying that weathering fluxes are lower than erosion fluxes.The value of ε CaO,MgO averages 28 ± 14% of ε total (mean and standard deviation).This represents the fraction of weathering attributable to CaO and MgO dissolution on average and is nearly identical to the 33 ± 9% inferred independently from solute concentrations (Section 4.1).

Erosion Fluxes
The data capture a ca.8-fold range in physical erosion fluxes (49-380 tons km 2 yr 1 ) (Table 1).Erosion fluxes increase from the Cordillera Central to the Sierra Luquillo (Figure 1) and are positively correlated with specific discharge (Figure 3).Although the study watersheds span a broad range in mean basin slope, measured erosion fluxes vary more strongly with precipitation and specific discharge than with slope (Table 1).The functional relationship between specific discharge and erosion and the underlying drivers of this correlation are beyond the scope of the present analysis; however, it is notable that similar correlations have been inferred in other likely nonequilibrium basaltic and andesitic landscapes using reconstructions of valley volumes and thermochronology (e.g., Ferrier, Huppert, & Perron, 2013;Reiners et al., 2003).This suggests that a runoff dependence of erosion rates may be a common feature of these landscapes.
The erosion flux estimated from 36 Cl in magnetite can be compared with independent constraints based on longterm sediment and solute gauging in watershed A1, the Rio Mameyes.The Rio Mameyes was sampled in Jan. 2019 during exploratory fieldwork and again in Jan. 2022.The fluxes inferred from both samples agree within 1σ, demonstrating reproducibility of the 36 Cl data (the inverse-error weighted mean of 2 measurements of the 2019 sample is 256 ± 32 tons km 2 yr 1 compared with 248 ± 57 tons km 2 yr 1 for the 2022 sample).Both measurements are consistent with the erosion flux determined from long-term sediment gauging (270 tons km 2 yr 1 ) and from an equilibrium model between solute fluxes and sediment geochemistry (265 tons km 2 yr 1 ) (Stallard, 2012a).This agreement supports the accuracy of the 36 Cl-based erosion fluxes in volcaniclastic watersheds in Puerto Rico and suggests that modern erosion fluxes in the Rio Mameyes are similar to average erosion fluxes over the last ca.6.5 ka, the averaging timescale of the cosmogenic nuclide signal in this watershed (Text S.1 in Supporting Information S1).

Weathering-Erosion and Weathering-Discharge Relationships
Total silicate weathering fluxes were determined from the product of specific discharge and TDS sil , and CaO and MgO weathering fluxes from specific discharge and [CaO + MgO].Total weathering fluxes span one order of magnitude (19 tons km 2 yr 1 to 187 tons km 2 yr 1 ) and are higher on average in the Sierra Luquillo than in the Cordillera Central (Table 1).Total weathering fluxes increase with erosion flux (Figure 4).The best-fit power-law slopes between weathering and erosion fluxes (Equation 1) are 0.90 ± 0.17 for total weathering and, almost identically, 0.95 ± 0.19 for CaO and MgO weathering (Figure 4).When interpreted in isolation, these relationships support a near linear dependence of weathering on erosion.However, weathering fluxes also increase with specific discharge (Figure 4).The power-law slope between the total weathering flux and specific discharge (Equation 2) is 0.63 ± 0.06, whereas the corresponding value for CaO and MgO is 0.56 ± 0.10.This implies that solute concentrations are not simply diluted with increasing specific discharge on the landscape scale (Figure S.4 in Supporting Information S1), consistent with expectations for a largely thermodynamically limited landscape.Disentangling the relative strength of supply versus thermodynamic controls on weathering fluxes in Puerto Rico requires addressing the strong covariation observed between erosion and specific discharge.

Differentiating Between Supply and Thermodynamic Limits
To differentiate between supply and thermodynamic limits, we analyze solute concentrations versus E/Q (Equation 3).The data set captures a 20-fold range in E/Q (34 μg cm 3 to 701 μg cm 3 ), which increases with decreasing specific discharge from the Sierra Luquillo to the Cordillera Central.The best-fit power-law slope is 0.59 ± 0.09 for all silicate-weathering derived oxides.A power law slope equal to 1 implies a purely supply limit, a slope of 0 a purely thermodynamic limit, and 0.5 equivalency between the two.Thus, this result indicates a slightly stronger supply than thermodynamic limit to total weathering.Likewise, the best-fit slope is 1.01 ± 0.18 for dissolved CaO and MgO, suggesting a near complete supply limit for Ca and Mg release and long-term CO 2 sequestration potential (Figure 5).The higher slope for CaO and MgO than for all silicateweathering derived oxides reflects a stronger thermodynamic limit for Si weathering than for the base cations.Plots for each of the five silicate-weathering derived constituent measured in water are presented in Supporting Information S1 (Figure S.5).The power-law slope for [SiO 2 ] versus E/Q is 0.31 ± 0.07.This indicates that Si weathering fluxes are mostly thermodynamically limited, consistent with the formation of low-solubility aluminosilicate clay minerals during weathering and with observations of Si chemostasis across global stream water data sets (Godsey et al., 2019).Significantly, analysis of solute concentrations versus E/Q indicates a much weaker dependence of weathering on specific discharge than obtained from the weathering-discharge relationship directly.
Although African dust deposition is important in Puerto Rico (McClintock et al., 2019;Pett-Ridge et al., 2009), dust does not challenge the conclusion of a significant supply limit.Modern dust deposition rates inferred from open-fall collection in the Sierra Luquillo vary by less than a factor of two and average 8.4 tons km 2 yr 1 (mineral + dissolved dust) (McClintock et al., 2019).Assuming complete weathering of dust, subtracting this flux from the measured weathering fluxes, and refitting parameters in Equation 3 reduces the power-law slope to 0.46 ± 0.09 (for TDS sil ).This likely represents a maximum effect because dissolution of quartz and Fe and Al oxyhydroxides in dust may be minimal and because dust deposition rates may scale with precipitation and be lower in the drier watersheds in the Cordillera Central, which have higher E/Q.
The largest weathering fluxes are found in watersheds at high elevations in the Sierra Luquillo.Many of these are located within the contact aureole of the Rio Blanco stock and may contain sulfide minerals.Sulfide oxidation
Although sulfide oxidation has been shown to be locally important to weathering front propagation in the Sierra Luquillo (Moore et al., 2019), sulfate fluxes are not correlated with total silicate weathering fluxes in the study watersheds (r 2 = 0.01, Figure S.6 in Supporting Information S1) and anionic load is dominated by HCO 3 and Cl (Figure S.2 in Supporting Information S1).This indicates that carbonic acid is the dominant source of acidity driving weathering and that the weathering-erosion-discharge relationship is not an artifact of spatially variable acid-sulfate weathering.
Although the study watersheds span a 4.2°C range in MAT, analysis of residuals in the regression model indicates that temperature is not a strong control on the observed weathering fluxes.Model residuals are poorly correlated with temperature, with r 2 = 0.13 for TDS sil and r 2 = 0.01 for [CaO + MgO] (Figure S.7 in Supporting Information S1).This suggests that chemical weathering fluxes from volcanic rocks in Puerto Rico are not subject to a kinetic limit.

Support From Sediment Geochemistry
The existence of a stronger supply than thermodynamic limit can be cross-checked by analysis of ε versus E/Q (Equation 6).The best-fit power-law slopes (b 3 -1) between ε values derived from sediment geochemistry and E/Q are 0.33 ± 0.15 for ε total and 0.23 ± 0.12 for ε CaO,MgO , which are equivalent to b 3 values in Equation 6 of 0.67 ± 0.15 and 0.77 ± 0.12 (Figure 6).The ε data set is subject to substantial geological uncertainties from, for example, non-representative sampling of watershed-averaged bedrock and regolith using stream sediment, geochemical mobility of Ti and Zr, or dust deposition.Yet, if the sum of all geological uncertainties affecting ε leads to random scatter, rather than a unidirectional bias, then the best-fit slopes should be robust.The best-fit slopes are consistent with those indicated by analysis of stream solutes versus E/Q within 2σ and support a stronger supply than thermodynamic limit to weathering fluxes.Residuals are poorly correlated with both temperature and erosion (Figures S.8 and S.9 in Supporting Information S1, r 2 = 0.09 in both instances for ε total and r 2 = 0.29 for temperature and 0.31 for erosion for ε CaO+MgO ), supporting the absence of a kinetic limit.Furthermore, because a water sample captures the instantaneous solute concentration, whereas regolith geochemistry develops across the thousand-year timescales of regolith residence on hillslopes, the agreement in the power-law slopes from Equations 3 and 6 indicates that long-term and short-term weathering fluxes are, on average, nearly in equilibrium.
The modest decrease in ε with increasing E/Q inferred from sediment geochemistry is further supported by the spatial distribution of soil orders across the landscape.The dominant soil order shifts from oxisols to ultisols from the Sierra Luquillo (series A) to the Cordillera Central (series B) as E/Q increases.Inceptisols are the most abundant soil order in the two watersheds with the highest E/Q (Moore et al., 2023;Staff, NRCS, 2023).This qualitatively confirms that regolith is highly weathered across the landscape, but that there is a small decrease in weathering intensity with increasing E/Q.

Effects of Dilution at High Discharge
In this analysis stream water chemistry is assumed to reflect the chemistry of the long-term average pore water in the regolith.However, the wet and dry season water samples do not capture the concentration-discharge relationship across the entirety of the annual hydrograph.To what extent are solute concentrations diluted at high discharge?More rigorous evaluation of chemostasis at high discharge is possible in three watersheds (A1, A10, and A16) in the Sierra Luquillo for which long-term, high-resolution concentration-discharge data are available (Wymore et al., 2017).Flow-weighted concentrations of major rock-weathering derived constituents are lower than the average of the wet and dry season samples by factors that range between 1.0 and 2.1 and increase with specific discharge (Q) (Table S.1 in Supporting Information S1).
If all stream water infiltrates through the regolith, then this implies that long-term average pore water solute concentrations are lower than the concentrations we measured.If the deviation between measured and flowweighted concentrations increases with Q, then solute concentrations will be reduced by a larger factor at high Q than low Q.This will increase the concentration versus E/Q power law slope.Approximating this effect using the observed dilution in the flux-weighted concentrations in the Sierra Luquillo gives power-law slopes of 0.88 ± 0.16 for TDS sil and 1.25 ± 0.26 for [CaO + MgO], which remain consistent with a supply limit (Figure S.10 in Supporting Information S1).
It is also possible that as Q increases overland flow becomes an increasingly important contributor to total streamflow and stream and pore water chemistry diverge from one another.In this case, an inverse relationship between solute concentrations in stream water and Q across runoff gradients may develop, even when weathering is thermodynamically limited, and pore water concentrations are invariant.However, it is unlikely that our water samples capture overland flow.Infiltration rates in the Sierra Luquillo are higher than typical precipitation rates (McDowell et al., 1992) and the solute concentrations are essentially chemostatic over the range of sampled discharges.If the lower flow-weighted solute concentrations in the Sierra Luquillo are a result of increasing overland flow that does not interact with the subsurface, overland flow should be subtracted from Q before calibrating the C versus E/Q power law, however the effect would also be to increase the power law slope and the apparent degree of supply limitation.

Implications for the Global Carbon Cycle
The trends observed in solute concentrations and ε across an E/Q gradient indicate that arc-rock weathering in montane, tropical landscapes is subject to a stronger supply than thermodynamic limit up to the maximum erosion fluxes observed in Puerto Rico but does not preclude a transition to a stronger thermodynamic or kinetic limit at higher erosion fluxes.Are the fluxes observed in Puerto Rico representative of basaltic and andesitic landscapes globally?Although few long-term erosion flux data are available in tropical volcanic landscapes, estimates of denudation from topographic reconstructions of volcanic edifices are broadly similar to the rates measured in Puerto Rico, for example, ranging up to 335 tons km 2 yr 1 on basalt in Kaua'i (Ferrier, Perron, et al., 2013), or 1350 ± 550 tons km 2 yr 1 on young arc volcanoes in the Lesser Antilles (Ricci et al., 2015).By comparison, denudation fluxes measured in Puerto Rico range from 84 ± 7 tons km 2 yr 1 to 568 ± 56 tons km 2 yr 1 when considering the sum of the physical and chemical fluxes (Table 1).Therefore, a strong supply limit may be generally applicable for modeling weathering of volcanic rocks in the humid tropics.Importantly, this result does not exclude a thermodynamic or kinetic limit in drier or colder climates.
An expression for the weathering flux from arc rocks as a function of specific discharge and erosion flux can be obtained by multiplying both sides of Equation 3 by specific discharge (Q).Using the best-fit slope for all ions from Figure 5 this gives: A similar expression may be written for CaO and MgO fluxes, but where weathering fluxes and long-term CO 2 sequestration potential are nearly linearly dependent on erosion and independent of discharge.The stronger supply than thermodynamic limitation to arc-rock weathering in the humid tropics supports the importance of tectonic uplift and associated increases in erosion flux in driving weathering fluxes from arc-rocks, similar to felsic lithologies (Ferrier et al., 2016;Larsen et al., 2014;Riebe et al., 2001).However, the impact on global CO 2 drawdown of uplift and erosion may be larger for mafic arc rocks than felsic rocks because of their greater CaO and MgO concentrations.This result is consistent with the hypothesis that volcanic arc landscapes are destabilizing features in the Earth's climate system.
In addition to the poor correlation of model residuals with temperature, two lines of evidence indicate that temperature does not control arc-rock weathering in Puerto Rico.First, this data set captures a 10-fold range in chemical weathering fluxes across a ca.4.2°C range in MAT.An Arrhenius relationship with the highest activation energy reported for basalt weathering (89.6 kJ/mole) (Brantley et al., 2023) predicts that weathering rates should increase by only 68% across this temperature range.The existence of a much larger spread in rates indicates that temperature cannot be the primary limit to weathering.Second, the high weathering intensity observed for CaO and MgO (CDF CaO and CDF MgO ) (Moore et al., 2023) implies that there is only modest potential for increased weathering of these elements from regolith with increasing temperature because, on average, 66% of CaO and 51% of MgO have already been lost.Even in the absence of any thermodynamic limitation, CaO and MgO fluxes can be increased by, at most, less than a factor of two without a concomitant increase in the supply of fresh rock.Instead, reaction rates are adequate at present temperatures in Puerto Rico to ensure release of most CaO and MgO from rock over the observed regolith residence timescales.Thus, temperature is not the primary limit to arc-rock weathering at MAT greater than or equal to ca. 20°C.

Reconciling a Supply Limit With Concentration-Discharge Relationships
The trends in ε and solute concentrations with E/Q support a stronger supply than thermodynamic limit to arc-rock weathering fluxes in the humid tropics and a nearly complete supply limit for CaO and MgO fluxes and long-term CO 2 sequestration potential.However, landscape-scale chemostasis is often cited in support of a strong thermodynamic limit to basalt and andesite weathering fluxes, wherein solute concentrations approach saturation and weathering is limited chiefly by the through-flux of water (Bluth & Kump, 1994;Dixon et al., 2016;Gaillardet et al., 2011;Maher & Chamberlain, 2014).Although in Puerto Rico solute concentrations decrease with increasing specific discharge by less than expected from simple dilution (indicative of partial landscape-scale chemostasis), this trend may be largely explained by the observed increase in erosion flux and fresh mineral supply with discharge for all elements except Si.
More generally, in landscapes where erosion fluxes increase in tandem with runoff, the associated increase in fresh mineral supply may support chemostatic concentration-discharge relationships across runoff gradients under supply limited weathering conditions.A dependence of erosion on precipitation similar to that found in Puerto Rico has been identified in other volcanic landscapes (e.g., Ferrier, Huppert, & Perron, 2013;Ferrier, Huppert, & Perron, 2013;Reiners et al., 2003).If weathering in these landscapes is typically supply limited, then an increase in erosion flux with discharge will generate an increase in weathering flux with discharge.Importantly, this would preserve a negative feedback control on climate under a supply limit, wherein increased intensity of the hydrologic cycle under a warming climate increases weathering fluxes and CO 2 consumption by increasing runoff, erosion, and fresh mineral supply.However, the mechanism for this feedback is hydrologic modulation of hillslope transport and channel incision processes, and thus fresh mineral supply, rather than a thermodynamic limit due to waters approaching equilibrium saturation.

Conclusions
This paper reports new measurements of solute concentrations, sediment geochemistry-derived weathering indices, and erosion fluxes across a runoff gradient on volcanic arc rocks in Puerto Rico.This data was used to evaluate the degree to which erosion and fresh mineral supply limit weathering fluxes on arc rocks in the humid tropics.Long-term erosion fluxes inferred from 36 Cl in magnetite illustrate a clear dependence of erosion flux on specific discharge.Analysis of erosion flux normalized by specific discharge (E/Q) versus solute concentrations and ε indicates that weathering fluxes from arc rocks in the humid tropics are more strongly limited by the ability of erosion to remove depleted regolith than by a thermodynamic limit due to waters approaching equilibrium saturation and that temperature is not a significant limit to weathering fluxes from arc rocks at typical erosion rates and temperatures in the humid tropics.We propose that landscape-scale chemostatic concentration-discharge relationships across groups of catchments may be reconciled with a strong supply limit where, as observed in Puerto Rico, erosion fluxes are positively correlated with specific discharge and that a coupling between erosion and specific discharge could maintain a negative feedback between arc-rock weathering and climate in supply limited weathering regimes.

Figure 2 .
Figure2.Specific discharge (Q) versus mean annual precipitation (MAP).The solid gray line shows the best-fit line for data outside of the Sierra Luquillo and the blue line for data in the Sierra Luquillo, which receives the highest precipitation in Puerto Rico.Dashed lines represent 2σ confidence intervals.The lower intercept of the Sierra Luquillo regression may reflect higher evapotranspiration in this densely forested region than on average across the island and the higher slope a decline in evapotranspiration with increasing elevation as cloud cover becomes increasingly persistent(Murphy and Stallard, 2012).
The measured[Zr]  rock /[Zr] sed and[Ti]  rock /[Ti] sed ratios differ by more than 2σ analytical errors in 14 of 27 watersheds.However, the offset is not clearly systematic; the mean and standard error of the ratios of [Zr] rock /[Zr] sed to[Ti]  rock /[Ti] sed is 0.91 ± 0.07.This indicates that any bias due to Ti mobility is likely small relative to the degree of stochastic variability.Therefore, to reduce noise, we calculate weathering indices using the average and standard error of the [Zr] rock / [Zr] sed and [Ti] rock /[Ti] sed ratios.

Figure 3 .
Figure 3. Erosion flux (E) versus specific discharge (Q).Erosion fluxes are positively correlated with Q and increase markedly above Q of approximately 200 cm yr 1 .Watersheds at low Q occupy different geomorphic positions, either lowrelief, relict surfaces or steep catchments that are graded to modern baselevel, which may account for the weaker correlation between E and Q at low Q.Data points are colored by TDS sil , which generally declines toward higher Q and lower E. Error bars show 1σ errors.

Figure 4 .
Figure 4. Weathering flux (W) versus erosion flux (E) and specific discharge (Q).Panels (a, b) show total weathering fluxes, whereas panels (c, d) show CaO and MgO weathering fluxes.Error bars show 1σ analytical errors.Uncertainties in Q are treated as unconstrained geological error when calculating linear regressions and are excluded from the error bars.The equations for the best-fit line and 1σ uncertainties are given at the top of each panel.Confidence intervals (2σ) for the best-fit line are shown by the dashed lines.Weathering fluxes are positively correlated with both E and Q.The correlation between W and Q in panels (b, d) is artificially elevated because W is defined as the product of Q and the solute concentration (Figure S.4 in Supporting Information S1 plots concentration vs. Q, which is not subject to a spurious correlation).

Figure 5 .
Figure 5. Dissolved silicate-weathering derived ion concentrations versus erosion flux divided by specific discharge (E/Q).Panel (a) shows TDS sil and panel (b) [MgO + CaO].Panel (c) is a schematic plot of the expected change in concentration (C) with E/Q under a supply and a thermodynamic limit.Solute concentrations are given in μg g 1 , which is equivalent to μg cm 3 for water of unit density.Best-fit linear equations with 1σ errors are given in each panel.The best-fit line is shown by the solid line and confidence intervals (2σ) by dashed lines.Error bars represent 1σ analytical errors.

Figure 6 .
Figure 6.ε versus erosion flux divided by specific discharge (E/Q).Panel (a) shows ε total and panel (b) ε CaO,MgO .Panel (c) is a schematic plot of the expected trend in ε versus E/Q under a supply and a thermodynamic limit.Equations for the best-fit lines are given at the top of each panel with 1σ errors.The best-fit line is shown by the solid line and confidence intervals (2σ) by dashed lines.Error bars represent 1σ analytical errors.