Lithological Control of Stream Chemistry in the Luquillo Mountains, Puerto Rico

Meteoric waters move along pathways in the subsurface that differ as a function of lithology because of the effects of chemical and physical weathering. To explore how this affects stream chemistry, we investigated watersheds around an igneous intrusion in the Luquillo Mountains (Puerto Rico). We analyzed streams on 1) unmetamorphosed country rock (volcaniclastic sedimentary strata, VC) surrounding an igneous intrusion, 2) the quartz-diorite intrusion (QD), and 3) the metamorphosed aureole rock (hornfels-facies volcaniclastics, HF). These lithologies differ physically and chemically but weather under the same tropical rain forest conditions. The sedimentary VC lithology is pervasively fractured while the massive QD and HF lithologies are relatively unfractured. However, the QD fractures during weathering to produce spheroidally-weathered corestones surrounded by cm-thick rindlets of increasingly weathered rock. Meteoric waters flow pervasively through the network of already-fractured VC rock and the spheroidally weathered rindlets on the QD, but only access a limited fraction of the HF, explaining why streams draining HF are the most dilute in the mountains. This results in various thicknesses of regolith from thick (VC) to moderate (QD) to thin or nonexistent (HF). The pervasive fractures allow groundwater to flow deeply through the VC and then return to the mainstem river (Río Mameyes) at lower elevations. These “rock waters” drive concentrations of rock-derived solutes (silica, base cations, sulfate, phosphate) higher in the lower reaches of the stream. Water also flows through weathering-induced fractures on the QD at high elevations where rindletted corestones are present in stacks, and this water flux dissolves plagioclase and hornblende and oxidizes biotite. This “QD rock water” is not generated at lower elevations in the Río Icacos watershed, where stacks of corestones are absent, and contributions to stream solutes derive from weathering of feldspar- and hornblende-depleted saprolite. The stream chemistry in the QD-dominated watershed (Río Icacos) thus varies from concentrated QD-rock water at channel heads below steep ridgelines toward more diluted “saprolite water” downstream. These observations emphasize the importance of lithology and fracture patterns in dictating water flowpaths, stream chemistry, and regolith development in headwater catchments.

Introduction

A central tenet of critical zone science is to explore how knowledge of geological structure and its evolution can help predict water flowpaths in the near subsurface to understand groundwater and stream water hydrologic and chemical linkages (Banks et al., 2009; Riebe et al., 2016; Bailey et al., 2019; Hayes et al., 2020; Brantley and Lebedeva, 2021). In some cases, the most important control on development of flowpaths and weathering may be dictated by the permeability and fracture density of underlying bedrock (Rempe and Dietrich, 2014), including fractures that are induced by the interplay between tectonic stresses, topography, and exhumation (St. Clair et al., 2015). In other cases, subsurface flowpaths may be dictated by porosity generated or occluded by geochemical reactions (MacQuarrie and Mayer, 2005; Lebedeva and Brantley, 2013; Brantley et al., 2017; Li et al., 2017). One way to explore subsurface flowpaths indirectly is to investigate longitudinal profiles of stream chemistry (e.g., Lawrence and Driscoll, 1990). Following that strategy in this paper, we compare watersheds on three lithologies in the Luquillo Experimental Forest (LEF) of the Luquillo Mountains in eastern Puerto Rico to explore relationships between lithology and hydrogeochemistry. In this paper, our intent is to emphasize physical and chemical attributes rather than ecological aspects of the system, recognizing that more work is needed to integrate physical, chemical, and ecological understanding of the critical zone in this location.

We examine the importance of lithology by combining subsurface observations with spatial analysis of stream solutes in three adjacent drainages—Icacos, Mameyes, and Sonadora—in the LEF. The Luquillo Mountains, rising from sea level to nearly 1,100 m elevation on Puerto Rico, are weathering and eroding at very high rates (e.g., Brown et al., 1995; White et al., 1998; Riebe et al., 2004; Dixon and von Blanckenburg, 2012; Brocard et al., 2015). The rapid denudation of the Luquillo Mountains is related in part to their high mean rainfall (MAP ∼3.5–4.5 m/a) and temperature (19–24°C) (Murphy et al., 2017). In addition, landsliding is common throughout the LEF regardless of lithology (Larsen et al., 1999), contributing to both physical and chemical denudation (e.g., Bhatt and McDowell, 2007).

The comparison of regolith development, chemical weathering, and stream chemistry on different lithologies is possible because of the well exposed contact metamorphic aureole at the core of the Luquillo Mountains. Here we define regolith as any weathered material lying above the parent lithology. Cretaceous volcaniclastic strata (VC) of basaltic and basaltic andesite composition were intruded by the Río Blanco quartz diorite (QD) pluton, and metamorphosed to hornfels (HF) facies near the contact (Seiders, 1971) at ∼46–48 Ma (Cox et al., 1977; Smith et al., 1998). The intrusion resulted in juxtaposition of our three study lithologies. These rocks and their contacts are hypothesized to have been been exhumed during rapid uplift of Puerto Rico since ∼4 Ma (ten Brink, 2005). This uplift, and the corresponding wave of erosion (Brocard et al., 2015), removed the Cretaceous VC and younger marine sediments, exposing the HF and QD. Today, the HF is exposed mostly in ridges and peaks of the mountains, and the QD underlying the bowl-shaped Río Icacos watershed (Figure 1).

FIGURE 1

FIGURE 1. Map of the study watersheds draining the Luquillo Mountains, northeastern Puerto Rico and their underlying major lithologies (color). Study sites (Table 1) within each watershed are shown with symbols: diamonds (Mameyes), squares (Espiritu Santo), circles (Rio Blanco). Additional sites sampled within a day of one another (synoptic sites) are also shown (see Supplementary Table S6). Within the Río Espiritu Santo watershed, the study profile is along the Quebrada Sonadora and within the Río Blanco watershed, the primary study profile is along the Río Icacos. The study profile in the Río Mameyes is along a major tributary (Río La Mina) in the upper reaches and along the mainstem Mameyes in the lower reaches. Groundwater wells in the Mameyes and Blanco watersheds are also indicated, including US Forest Service (USFS) wells. Letters for groundwater wells are explained in Table 3. The major knickzone is shown with hachured lines. The inset map shows sites in the Icacos catchment.

We use new geochemical data from streams and groundwater to elucidate subsurface flowpaths and weathering reactions. We also summarize new data for the Sr isotopic composition of waters and for atmospheric tracers (CFCs, SF₆, ³H, noble gases) in groundwater (see Supplementary Material) to infer dominant weathering reactions and estimate residence time for a subset of groundwater samples.

Study Area

The Three Study Watersheds

The Luquillo Mountains lie in tropical latitudes under strong marine influence, resulting in a warm, wet climate (McDowell et al., 2021). Temperature and rainfall (Table 1) are correlated with elevation and aspect because of the predominantly northeasterly winds, which cause orographic effects resulting in higher precipitation rates in the LEF than coastal areas (Murphy et al., 2017).

TABLE 1

TABLE 1. Location and watershed characteristics of stream study sites.

Three gauged watersheds, each with >95% of their area within the LEF, were the focus of our research: the Río Mameyes near Sabana (USGS 50065500; also referred to herein as Mameyes Puente Roto- MPR), Río Icacos near Naguabo (USGS 50075000), and Quebrada Sonadora near the El Verde Field Station (USGS 50063440). Data were also used from additional sites, some with stream gauge records, in the downstream reaches of the watersheds. The three watersheds were chosen because they each typify the hydrology and geochemistry of one of the three main lithologies for which we hypothesize the importance of lithological controls on weathering and solute export via streams: VC underlies much of the Río Mameyes watershed (from Table 1, range in mean annual precipitation (MAP) for the stream sites studied here = 3.23–4.09 m); QD underlies the Río Icacos (MAP = 3.25–4.06 m) and most of its larger parent watershed, Río Blanco; and HF largely underlies the Quebrada Sonadora (MAP = 3.52–4.34 m) within the larger watershed of the Río Espiritu Santo (Figure 1). Temperatures vary little seasonally, and MAT is 23.7 at 352 masl (Bisley) and 20.1 at 1051 masl (East Peak) (McDowell et al., 2012). Regolith thickness and production rates both decrease in magnitude from the VC to the QD and then to the HF, as described below. The watersheds, long studied as part of the U.S. Long Term Ecological Research (LTER) and Critical Zone Observatory (CZO) networks (see summary in Table 1), are described in the next sections.

Río Mameyes Watershed: Predominantly VC

The Río Mameyes drains northward to the Atlantic Ocean and is dominated by the lithology (VC) that defines the lower elevations of the mountains. This is the sedimentary rock into which the QD intruded. Only the very uppermost reach of the Mameyes (subcatchment, Río La Mina) is partially underlain by QD and HF. Unlike the Icacos, no major knickpoint is observed along rivers on the VC (Porder et al., 2015) and the Mameyes watershed is significantly steeper on average than that of the Icacos (Larsen, 1997). Estimates of regolith formation on VC using U-series nuclides yield rates that are ∼5–10 higher than those on QD (Dosseto et al., 2012). Consistent with this, regolith on VC in the intermediate elevations of the watershed can be > 100 m thick (Fletcher and Brantley, 2010; Buss et al., 2013). The Mameyes watershed is thus a deeply incised, VC-dominated watershed with very thick regolith, and a steep longitudinal profile.

Río Icacos Watershed: Predominantly QD

The Río Icacos is the upper portion of the Río Blanco watershed that drains southward from the Luquillo Mountains to the Caribbean Sea. The watershed is underlain almost completely by the Río Blanco quartz diorite (QD) pluton (Figure 1). Significant erosion into the pluton created a steep-walled, bowl-shaped watershed defined by ridges, typically mantled with the contact metamorphosed HF, and QD at moderate elevations underlying the ridges and the rest of the watershed (Orlando et al., 2016). Within this bowl, the floodplain of the Río Icacos bottom is broad, gently sloping, and characterized by thick and highly weathered soils and saprolite with few exposures of QD corestones. Here, saprolite is defined as regolith formed isovolumetrically in place that often retains evidence of the original structure of parent material.

The “bowl” created by the Río Icacos is truncated at its lower elevation by a major knickzone that roughly occurs where the river has cut through the ring of HF on its southern limit (Figure 1): in the upper part of the knickzone the elevation drops from 600 to 25 masl over a distance of ∼3.5 km (Comas et al., 2019). Several other tributary streams join the Icacos in this knickzone; below this zone the river is known as the Río Blanco. Where the river has breached the ring of HF at about 450 masl, the knickzone steepens, dropping 350 m in ∼1.5 km. Above the knickzone, the denudation rates and regolith production rates on the QD estimated from ¹⁰Be and U-series disequilibrium equal 43 m/Myr (Brown et al., 1995) and 45 m/Myr (Chabaux et al., 2013) respectively. Landscapes in the knickzone are eroding approximately twice as fast as those above (Brocard et al., 2015).

Where seen at depth in outcrop, the unweathered QD is largely unfractured. Most of the exposure is along Route 191. At elevations in the watershed above the knickzone (see Figure 1), widely spaced lineaments cross the road as viewed in planform (Comas et al., 2019). These lineaments, observed in aerial photographs, are manifested in the field as gullies roughly perpendicular to the road and the Río Icacos. These are inferred to be deep fracture zones at least partly because ground penetrating radar (GPR) reveals them to be zones of deep GPR reflections (Orlando et al., 2016; Hynek et al., 2017; Comas et al., 2019). They are interpreted as fracture zones that permit and focus downward flow of meteoric water. Where we have drilled through them or observed them in outcrop, they reveal rounded corestones that have formed by spheroidal weathering one on top of the other in vertical stacks (Turner et al., 2003; Buss et al., 2008; Comas et al., 2019). In contrast, on the interfluves between the gullies, only single rindletted corestones (not stacks of corestones) are observed to overlie massive QD (e.g., Guaba Ridge, LGW1 core (Orlando et al., 2016)). In other words, weathering in the deep fracture zones results in stacks of corestones one on top of each other while weathering between the vertical fracture zones (in the interfluves) is characterized by a single corestone layer separated by rindlets from the underlying massive QD. Regardless of whether a single corestone is observed or a stack of multiple corestones, quartz-rich saprolite uniformly overlies corestones throughout the watershed.

Quebrada Sonadora Watershed: Predominantly HF

In the LEF, no major watershed is entirely underlain by the ridge-forming HF. To focus on a watershed largely on HF, we sampled a steep high-elevation tributary, the Quebrada Sonadora, which is underlain by 100% HF in its upper reaches; these upper sites also have the highest MAP (Table 1). Quebrada Sonadora is a sub-catchment of the Río Espiritu Santo watershed that drains the northwestern slopes of the Luquillo Mountains; the lower reaches of the Río Espiritu Santo predominantly drain VC. Very little is known about weathering and erosion rates on HF, but based on ¹⁰Be, the ridges and peaks comprised of HF erode at 4 m/Ma (unpub., Brocard, pers. comm.), at least an order of magnitude slower than either the QD or the VC.

Mineralogy of Bedrock and Regolith

Volcaniclastic

The VC is comprised of the Hato Puerco and Tabonuco Formations. These marine strata are clinopyroxene-bearing volcanic sandstones of basaltic to andesitic composition with mudstone, volcanic breccia or conglomerate interbeds (Seiders, 1971). The Tabonuco contains calcareous mudstones that are pyrite-bearing.

In decreasing abundance, the VC contains plagioclase > chlorite > quartz > pyroxene > epidote > alkali feldspar > amphibole ± trace prehnite, illite, calcite, and kaolinite. As a more mafic rock than the QD pluton, the VC generally has higher concentrations of Mg, Ca, and Fe. The VC contains roughly 35% plagioclase, 5.8% potassium feldspar, 24% chlorite, 9.4% pyroxene, 3.6% amphibole, 8% epidote, and 10% (but variable) quartz (Buss et al., 2013). Much of the volcaniclastic strata were hydrothermally altered, explaining why the compositions show a bimodal distribution of SiO₂ content ranging from 46 to 61 wt%: in general, the lower values represent the parent SiO₂ composition while the higher values represent parent composition with added hydrothermal silica. Alteration replaced volcanic glass, plagioclase, and biotite with chlorite, illite, and microcrystalline quartz. Alteration of pyroxenes (augite) to hornblende (actinolite) is also observed (Buss et al., 2017; Moore et al., 2019).

Much of what is known about weathering of the VC within the Mameyes was learned from studies in the subcatchment known as Bisley that extends from ∼50 to 400 masl and where the weathering profile on the VC is very thick (Fletcher and Brantley, 2010; Buss et al., 2013; Buss et al., 2017). For example, a geophysical survey documented regolith that is at least 60 m thick in places (Schellekens et al., 2004). In the VC-underlain landscapes, corestones are very angular as compared to the spherical corestones in the QD landscape and are observed everywhere on the land surface except at the highest elevations. Fletcher and Brantley (2010) used the observed distribution of angular corestones exhumed at the land surface at moderate elevations to calculate that the depth of regolith in Bisley is ∼135 m, reasonably consistent with observations from drilling (Buss et al., 2013). Evidence from Schellekens et al. (2004) and Buss et al. (2013) emphasizes that the regolith extends to great depths below the stream bed, documenting the importance of subsurface groundwater flow through the VC. Buss et al. (2017) reports that VC regolith in Bisley is dominated by quartz, hematite, goethite, kaolinite (up to almost 70 wt%), and illite (up to almost 20 wt%).

Quartz diorite

The SiO₂ content of the Río Blanco QD pluton sampled from outcrops ranges between 55 and 70 wt%, with an average of 63.1 ± 4.3 wt% (1σ, n = 12) (Turner et al., 2003). The composition varies from quartz diorite (typically 5–20% quartz) to very rare diorite (i.e., <5% quartz) (Orlando, 2014). The QD contains ∼60% plagioclase and 25% quartz, and thus is more felsic than the VC (Navarre-Sitchler et al., 2013). The plagioclase crystals have an average composition of An₄₈ and range from 200 to 1,000 µm in dimension (Buss et al., 2008). The abundance of minerals in the QD decrease roughly in the order, plagioclase > quartz > hornblende > biotite ± alkali feldspar. Accessory minerals include Fe-Ti oxides, apatite, zircon, titanite, tourmaline, epidote, chlorite, Fe sulfides, pyroxene, and calcite (Seiders, 1971; White et al., 1998; Turner et al., 2003; Buss et al., 2008). Veins and fractures containing pyrite have been pervasively oxidized from the surface down to ∼6 m depth (Brantley et al., 2011) and in some fractures to greater depths. For example, oxidized pyrite has been observed deeper than 35 m in the QD where it is capped with a thin veneer of HF (Orlando, 2014; Orlando et al., 2016).

The massive QD stock, where it is exposed and unweathered, exhibits only a few widely spaced fractures. These and other exhumation-related vertical fractures apparently allow spheroidal weathering to initiate, perhaps because biotite oxidation causes a slight expansion (Fletcher et al., 2006; Buss et al., 2008). Both micro- and macro-fractures are observed to form. Macro-fractures tend to wrap around spheroidal corestones like onion skins to sequentially form rindlets that are roughly 1–2 cm in thickness. This fracturing allows meteoric water to percolate along the newly exposed mineral surfaces, oxidizing more biotite and dissolving plagioclase and hornblende within and along the rindlets, releasing alkali cations Na and K (Buss et al., 2008; Eq. 4) and alkaline earth cations (Mg, Ca) (Buss et al., 2008; Eq. (5)) respectively. Fractures are large enough to also allow transport of partially weathered particles in the subsurface (Gu et al., 2021).

The release of base cations and dissolved Si as plagioclase and hornblende dissolve is limited mostly to rindlets. In the outermost rindlets, weathering leaves behind disaggregated saprolite that mostly contains quartz, oxidized biotite, recalcitrant accessory minerals, and secondary kaolinite (Murphy et al., 1998; Schulz and White, 1999). The predominant reaction in the saprolite is weathering of biotite and quartz to release K, Mg, and Si (White et al., 1998).

These observations document that at the scale of an individual corestone, the depth interval across which plagioclase dissolves (the reaction front for plagioclase) typically comprises about 20 rindlets of ∼2 cm each (Fletcher et al., 2006). However, if the entire profile from the massive igneous bedrock to plagioclase-free saprolite is considered, the reaction front is much thicker (Turner et al., 2003; Brantley et al., 2011; Comas et al., 2019): it can cross through inferred-fracture zones of focused weathering and highly altered rindlets surrounding multiple stacked corestones above bedrock.

Hornfels facies rock

The protolith for HF and VC is compositionally identical, but the HF has undergone high temperature hornfels-facies metamorphism because of emplacement of the Río Blanco stock into the VC. Orlando (2014) measured SiO₂ contents from 48 to 72 wt% in HF, and, like the VC, generally higher concentrations of Mg, Ca, and Fe compared to QD. The mineralogy of the HF includes both volcanic and metamorphic clinopyroxene, plagioclase, and chlorite, and locally it can include wollastonite, garnet, epidote, actinolite, blue-green hornblende, biotite and scapolite (Seiders, 1971). Microcrystalline quartz is also locally abundant along fractures and bedding planes.

The regolith thickness on HF is less well constrained, but on average is relatively thin, consistent with the observation of ridges with large outcrops of HF and long reaches of the stream channel that flow over bedrock. Angular corestones and clasts litter the land surface along the HF ridges (see discussion in Lebedeva and Brantley, 2017) and intact bedrock can be observed where small springs emit. Thus, observations of depth to HF bedrock on ridges range from 0 (rocky outcrops) to 5–10 m in zones of HF landsliding (augerable depth (Orlando, 2014; Orlando et al., 2016).

Drainage on the HF ridges also tends to be poor. For example, the only true bogs in the Luquillo Mountains are located on topographic highs underlain by HF (c.f. Ogle, 1970). One of these bogs was mapped near well site EP1 (site G on Figure 1), the only borehole drilled into HF. The hole is situated on a ridgeline in a saddle between two rocky HF peaks (peaks rise 35–55 m above the drill site). The upper 7.5 m of regolith at EP1 has been interpreted as HF-derived whereas the lower ∼35 m is QD-derived regolith (Orlando, 2014). Thus, the borehole traversed through the entire HF regolith into the underlying QD pluton.

Materials and Methods

Sample Collection

Streams

We report sampling of streams during 1–2 day campaigns in periods without large storms between 1997 and 2007 (full data presented in Supplementary Material). We refer to these as synoptic samplings because they were taken within a few days of one another on the same streams; thereby proving a spatial snapshot of the stream chemistry under the same hydrologic conditions. In the Río Icacos/Río Blanco, samples were taken at 13 sites ranging from the mountains at 640 masl to almost sea level (3 masl); in the Río Mameyes, samples were taken at seven sites from 754 to 4 masl; and in the Quebrada Sonadora/Río Espiritu Santo samples were taken at 12 sites from 992 to 14 masl (cf. Figure 1 and Table 1). Owing to differences in difficulty of access, ∼70 synoptic samplings were completed in the Mameyes and four to eight in the other watersheds. Major solutes were analyzed for all ∼700 samples (full data available in Supplementary Material).

Several additional sets of stream samples were taken for targeted geochemical measurements including a few for Sr isotopic analysis. First, baseflows from the three drainages were sampled between March 10 and 28, 2013 during a relatively dry period. Select locations within the Río Icacos watershed and at Río Mameyes Puente Roto (MPR) were also sampled beginning 29 March 2013 during a strong storm from the northwest which lasted ∼48 h. In February 2014, synoptic sampling of stream solutes in Río Icacos at near-median discharge values was paired with discharge measurements made by acoustic Doppler velocimeter (Sontek Flowtracker).

Groundwater

Springs are relatively rare on the VC and HF. On QD, groundwaters were sampled at seeps or springs on an occasional basis at several elevations along Route 191. Groundwater chemistry for four wells drilled in the Río Icacos watershed between 2008 and 2014 are also reported. These values are compared to previous groundwater compositions sampled from U.S. Forest Service potable-water wells drilled to >100 m depth along Route 191 (these are referred to throughout as the USFS wells) as reported by Scholl et al. (2015). A description of the wells is provided in Supplementary Material (see Supplementary Figure S3) and locations are shown on Figure 1 with key in Table 3.

Analytical Methods

Water chemistry

Stream water samples were filtered directly into acid washed polyethylene bottles in the field using pre-combusted glass microfiber filters with a 0.7 µm pore size (Whatman GF/F). Samples were frozen for analysis of major cations and anions by ion chromotagraphy as well as analysis of dissolved organic carbon (DOC) and total dissolved nitrogen (TDN) using high-temperature combustion. Total dissolved silica, soluble reactive phosphorus, and ammonium were measured using robotic colorimetric analysis (details given in McDowell et al., 2021). Stream samples were shipped to the NH Water Resources Research Center, University of New Hampshire for analysis. Analytical reproducibility for each analysis was typically <5% and the glass fiber filters had no measurable effects on Si concentrations (McDowell et al., 2019).

Groundwater (including seeps and wells) and surface water samples from 2012 to 2014 were filtered with 0.45 µm nylon syringe filters in the field into two 30 ml acid washed polyethylene bottles and stored refrigerated until analysis at the Laboratory for Isotopes and Metals in the Environment (LIME), Penn State University. One bottle was acidified in the field with several drops of ultrapure HNO₃, and used for analysis of cations, minor elements, and dissolved silica by inductively coupled plasma-atomic emission spectroscopy (ICP-AES). The other bottle was used for analysis of anions by ion chromatography. For major elements, the analytical error is on the order of a few percent.

A subset of water samples was analyzed for Sr isotopic composition using a Thermo-Finnigan Neptune Plus ICP-MS at the University of Utah using published methods (Chesson et al., 2012). Reported ⁸⁷Sr/⁸⁶Sr ratios were corrected for mass bias using an exponential law, normalizing to ⁸⁶Sr/⁸⁸Sr = 0.1194. During the course of analysis the ⁸⁷Sr/⁸⁶Sr ratio of SRM 987 was measured as 0.71030 ± 0.00002 (2σ SD, n = 16, MSWD = 2.36).

A small subset of waters were also measured for groundwater tracer concentrations and interpreted with respect to residence times in the subsurface as described in the Supplementary Material with data in Table 5 and in Supplementary Material (Supplementary Table S5).

Results

Stream chemistry

Here we discuss a few general trends in water chemistry for streams sampled in this study (Tables 1, 2; Supplementary Tables S1–S3, S6; Figure 1). At upper elevations, most of the anions (Cl⁻, NO₃⁻, SO₄^2-, PO₄^3-) and the ammonium cation (NH₄⁺) are found in all three streams at concentrations similar to or slightly above that of precipitation (Supplementary Figures S1, S2). For example, most sulfate concentrations in all three watersheds (18.2 ± 13.7 µM; 2σ) lie within the range for precipitation (8.8 ± 7.7 µM; 2σ) collected at El Verde field station (NADP site PR20). In all streams where solutes are present at higher concentrations than the chemistry of rainfall, cation concentrations were positively correlated with Si, as expected if they are indicative of silicate mineral-water reactions (Supplementary Figure S2).

TABLE 2

TABLE 2. Average stream chemistry values.

Overall, concentrations of most solutes (e.g., Si, phosphate, nitrate, sulfate, and dissolved inorganic carbon (DIC)) in the stream on HF (Quebrada Sonadora, Supplementary Table S2) did not change to any great extent as a function of downstream position (plotted in Figure 2 as watershed area drained at each mainstem position). Solute concentrations on Sonadora are also generally low relative to the other two watersheds (Figure 2; Table 2, Supplementary Table S2).

FIGURE 2

FIGURE 2. Stream chemistry plotted as a function of upstream watershed drainage area for each of the three watersheds shows different patterns depending upon lithology. Average values for dissolved Si and the sum of the seasalt-corrected base cations (Σbc*) are shown for sites in each watershed that are >95% within LEF. Groundwater chemistry from wells in the watersheds (see Table 3 and Supplementary Data) are shown within the appropriate watershed, with reference to approximate sampling depth in meters below the land surface (mbls). Wells near drainage divides are indicated for both watersheds.

In contrast, concentrations of these analytes generally increased downstream on the VC (Río La Mina/Mameyes) and decreased downstream on the QD (Figure 2, Supplementary Table S3). Specifically, high sulfate (∼35 µM), moderately high phosphate (∼0.4 µM), and moderately high Si (∼370 µM) concentrations were observed at low elevations in the Mameyes watershed while high nitrate (∼20 µM), phosphate (∼1 µM), and Si concentrations (∼440 µM) were observed high in the QD watershed (Río Icacos) as shown in Table 2 and Figure 2.

Atmospheric input-corrected water chemistry

To use the chemistry of weathering to trace fluid flow in the subsurface, we focus on cations derived from weathering of bedrock. To do this, we employed a standard correction for seasalt inputs from atmospheric deposition (Stallard, 2012). These data were corrected based on chloride concentrations in water samples and are denoted with an asterisk (*) in Table 2. The sum of the corrected base cations Na*, K*, Mg*, and Ca* is noted here as Σbc* and is used as a measure of the total cation concentration (µM) that is likely to be derived from weathering of rocks and dust. Consistent with the correction successfully isolating only solutes contributed from rock or dust weathering, Σbc* does not exceed dissolved Si (µM) in all three watersheds unless a portion of the drainage area is outside the LEF where anthropogenic activity generally increases, especially at the lowest elevations. The ratio of Σbc*/Si in most non-tropical watersheds can vary based on the weathering intensity, but in this tropical landscape, weathering is of very high intensity overall. In addition, although the Si and Σbc* concentrations in groundwater vary with depth (Table 3), they generally bracket the stream solute data (Figure 2).

TABLE 3

TABLE 3. Average groundwater chemistry values.

The seasalt-corrected concentrations as well as uncorrected concentrations demonstrate longitudinal patterns distinct to each lithology (Figure 2). For example, the concentrations of Si, Σbc*, phosphate, and dissolved inorganic carbon (DIC) remain roughly constant in the Sonadora (on HF) from upstream to downstream. In contrast, these analytes become more concentrated with increasing watershed area within the Mameyes watershed (on VC), and more diluted with increasing watershed area in the Icacos (on QD) (Figure 2). Thus, weathering derived solutes either stay constant (HF), increase (VC), or decrease (QD) downstream on the study rivers.

Groundwater chemistry and borehole observations

Table 3 shows that the Si concentrations in groundwater vary with depth, and Figure 2 shows that these concentrations generally bracket the stream solute data. The concentrations of Si, Σbc*, and specific conductance in groundwater were observed to be highest on the VC in the USFS wells, located at low elevations within the Mameyes watershed (Figure 1, 3). These groundwaters were also the deepest waters sampled and were likely to have therefore followed the longest flowpaths.

FIGURE 3

FIGURE 3. Summary of groundwater chemistry in the Luquillo Mountains (see Table 3 and Supplementary Data): (A) Wide ranges of pH values are observed at sample sites with generally minor differences in specific conductance. Strong differences exist between LGW2 wells B and C, which are less than 2 m apart. The deepest groundwater (USFS wells) show notably higher specific conductance. (B) Dissolved Si generally increases with the sum of base cations corrected for seasalt inputs (Σbc*) as expected for weathering-derived solutes.

On the QD, in contrast, the highest specific conductance and concentrations of Si in waters were sampled in EP1, LGW2B and LGW2C (as shown in Figure 1, Supplementary Table S4), located high in the watershed (Figure 4). As shown in Figure 4, the EP1 well was drilled through HF at the surface, but cut through a 40-m thick section that appears to be a highly-weathered stack of QD corestones that have developed in situ one on top of one another with intervening depth intervals of saprolite-like material. The deepest oxidized fracture was observed at about 37 m. The LGW2B and LGW2C wells were also drilled high in the watershed but directly in QD on route 191. Those wells were drilled through stacks of heavily weathered corestones characterized by many fractures (Figure 4). In some cases, voids were encountered that caused the drill bit to drop during drilling (Figure 4). The deepest oxidized fractures were observed at about 25 m in the LGW2 wells. These LGW2 stacks of corestones are thought to have formed in one of the deep vertical fracture zones that were identified with ground penetrating radar (GPR) (Orlando et al., 2016; Hynek et al., 2017; Comas et al., 2019). These zones, observed at the land surface as gullies or valleys, crisscross the watershed at wide spacing intervals high in the watershed but at close spacing near the large knickpoint on the river (Comas et al., 2019).

FIGURE 4

FIGURE 4. Logs of two sites almost entirely drilled through QD (LGW1, LGW2) and one drilled through HF at the surface and QD at depth (EP1). LGW2 hosts three boreholes (A,B,C). Sampling was only possible in B and C. Depth intervals of saprolite, fractures, veins, etc. are all summarized to show the structure of the regolith. Water table measurements are from field campaigns spanning 2012–2015. Measurements of in situ parameters as a function of depth throughout the full water column in well LGW2B are also shown for three dates as labelled by symbol shading. Locations of wells is summarized in Table 3.

In contrast to these deeper, high-elevation wells, the specific conductance and Si concentrations were lower in groundwaters sampled in LGW1, a borehole drilled at lower elevations on the QD (location D on Figure 1). Borehole LGW1 was drilled in an interfluve between the linear vertical fracture zones where the total depth of weathering is shallower than in EP1 or LGW2. Specifically, borehole LGW1 was drilled into a lower toeslope near Guaba located approximately halfway between the uppermost elevation of the watershed and the upper part of the knickzone where extensive GPR and geochemical analysis has been completed on one outcrop (Orlando et al., 2016; Comas et al., 2019). The entire zone from land surface to massive unaltered bedrock drilled in LGW1 (Supplementary Figure S3) was ∼5 m: the borehole crossed one corestone and one 50-cm thick set of rindlets across which plagioclase and hornblende weathered. The last ∼20 m of the borehole intersected unfractured, unweathered rock (Figure 4). Waters were sampled from the shallow depths near the rindlets (see description below). The deepest oxidized fracture was located just below the rindlet zone.

Waters in LGW1 were chemically similar to those sampled in seeps at moderate elevations along Route 191. At those seep sites, a few corestones were exposed (Orlando et al., 2016; Comas et al., 2019). Occasionally, pores between corestones were observed to be so large that they could be entered. It was common to see fast-flowing springs emit from such stacks of still-in-place corestones during wetter periods. The fast-flowing water also was observed in some places to transport relatively large particles in the subsurface (Gu et al., 2021). At lower elevations, one stream, the Guaba, was fed by such seeps near locations A and B on Figure 1. The specific conductance and Si content of the seeps and the Guaba were similar to waters in LGW1 (Figure 3).

Groundwaters were measured for dissolved oxygen (DO), pH, and specific conductance in LGW2B on three separate days over 3 years by depth-profiling (Figure 4). Strikingly large differences in all three of these parameters were observed both as a function of depth and over time. On one day, very little variation was observed with depth; however, on two other days the DO decreased downward while the specific conductance increased downward (Figure 4).

Strontium isotopes

During March 2013, we sampled surface water at baseflow as well as some waters on the rising limb of the hydrograph of a late March/early April storm to constrain weathering reactions in the different river systems (Table 4). Our average baseflow measurement of ⁸⁷Sr/⁸⁶Sr for the Río Icacos on QD (0.7052) is similar to previously reported data collected for the river (0.7055) (Pett-Ridge et al., 2009a; Pett-Ridge et al., 2009b). These values are lower than average values for streams draining HF (0.7067) and higher than values for the VC (0.7044) (Supplementary Table S1). In contrast, the Sr in whole-rock VC and QD are both reported to be 0.7041 (Porder et al., 2015); therefore, ⁸⁷Sr/⁸⁶Sr ratios in baseflow in the Mameyes (on VC) are similar to the bedrock. The VC has approximately twice the Sr content of HF whole rock.

TABLE 4

TABLE 4. Strontium isotope data for surface and groundwater samples.

The differences are also shown on endmember mixing diagrams in Figure 5A,B for the VC/HF and QD respectively. Figure 5A demonstrates that baseflow stream chemistry Sr on VC/HF can largely be explained as a mixture of cloudwater, throughfall, and plagioclase Sr (mineral data from published work in Pett-Ridge et al., 2009a; Pett-Ridge, 2009; Pett-Ridge et al., 2009b). Likewise, Figure 5B,D show that Río Icacos water (on QD) is a mixture of Sr from plagioclase and hornblende as well as throughfall. Mixing diagrams for storms are shown in Supplementary Figure S4.

FIGURE 5

FIGURE 5. Sr isotope geochemistry of streams in the Luquillo Mountains. Two endmember mixing models are consistent with the ⁸⁷Sr/⁸⁶Sr of the high concentration endmember and broadly constrain the low concentration endmember. Samples during the March 2013 period of baseflow and on the rising limb of the hydrograph during the early April 2013 storm demonstrate that: (A) baseflow on VC- and HF-dominated watersheds can be explained as a mixture between a high concentration endmember with low Sr isotope ratios and a low concentration endmember with ⁸⁷Sr/⁸⁶Sr values similar to cloudwater and throughfall precipitation, (B) the storm response of the QD streams also predicts a high Sr endmember with low Sr isotope ratios and a low concentration endmember dominated by rainfall, (C) baseflow ⁸⁷Sr/⁸⁶Sr in VC- and HF-dominated watersheds is a function of lithology and these relationships predict that Sr flux is dominated by a mineral(s) with ⁸⁷Sr/⁸⁶Sr values of ∼0.7040, presumably Ca-rich plagioclase with Sr isotope composition similar to that of the Río Blanco stock, and (D) streams on quartz diorite (QD) show a hydrochemical response to the storm that also predict a high concentration endmember with low ⁸⁷Sr/⁸⁶Sr (∼0.7044) that is similar to groundwater in the Río Icacos watershed (with exception of LGW2B). Precipitation inputs in the area have been reported to yield ⁸⁷Sr/⁸⁶Sr values of 0.7103; when additions of radiogenic Sr from Saharan dust (Pett-Ridge et al., 2009a,b) are included, the total atmospheric input is estimated at 0.7134. Other data for atmospheric inputs and minerals are from Pett-Ridge et al. (2009a, b). See Figure 1 and Table 3 for groundwater sample sites labelled in panel (d) (also, Supplementary Data).

The higher Sr concentrations in the VC than the HF are consistent with the mixing line for the Río Mameyes (VC) and Espiritu Santo (HF) plotted as a function of percent HF and VC (Figure 5C). Figure 5C shows that most of the Sr in rivers overlying VC could derive from a mineral with ⁸⁷Sr/⁸⁶Sr values near 0.7040, presumably plagioclase. Clearly, another mineral(s) with higher ⁸⁷Sr/⁸⁶Sr releases Sr to the river on HF where the Sr concentrations in the bedrock are lower. Only one groundwater sample (from groundwater well LGW2B high in the Icacos watershed as shown in Figure 1 symbol E) also showed this higher ⁸⁷Sr/⁸⁶Sr signature (Figure 5D).

Figure 5D emphasizes that even closely spaced wells on QD (LGW2B and LGW2C) can show very different ⁸⁷Sr/⁸⁶Sr values. While data for LGW2C groundwater are plotted on Figure 5D because the borehole is situated in the Icacos (mostly QD) watershed, the water shows higher ⁸⁷Sr/⁸⁶Sr values that are more similar to Sr released from HF and some HF was observed in the LGW2 borehole. Figure 5D also shows that ⁸⁷Sr/⁸⁶Sr varies in the Guaba and Icacos (on QD) during baseflow versus stormflow (see, also, Supplementary Figure S5). The baseflow-stormflow data series plot along a mixing line documenting a low endmember value at about 0.7044 (like Sr in plagioclase, and like waters in LGW2C, LGW1, and a Route 191 seep, shown as symbols F, D, and C on Figure 1) that is observed at base flow, and an endmember value slightly higher than Sr in hornblende or in LGW2B (symbol E) that is observed in storm waters.

Discussion

Groundwater residence times

Water flows quickly through most of this tropical landscape, with all groundwater samples recharged since 1960 and the majority of groundwater recharged within the last 8 years. Some differences in residence time were observed among the wells, which we attribute to differences in critical zone structure as has been observed in other locations (e.g. Braun et al., 2009). For example, using a piston-flow model for tritium, residence times were calculated around 5–8 years for water samples from the LGW1 well located at moderate elevations in the Río Icacos watershed on QD (Figure 6). These are the maximum residence times estimated from ³H concentrations and are in general agreement with other estimates for the residence time for waters in Luquillo (White et al., 1998). The LGW1 well is located in an interfluve between the deep vertical fracture zones: thus, the waters in the borehole travelled only through saprolite and a single rindletted corestone layer (see Figure 4; Supplementary Figure S3).

FIGURE 6

FIGURE 6. Tritium concentration and weathering-derived solutes for surface and groundwater in the Luquillo Mountains (see Supplementary Material). (A) The measured concentration in tritium units (TU = 1 tritium atom per 10¹⁸ hydrogen atoms) presented on the y-axis is used to estimate mean residence time of the water (x-axis) using the radioactive decay equation (half-life = 12.43 years) assuming piston flow. The age models are predicated on the assumptions that variability in precipitation is approximately equal to the analytical precision, that residence times are sufficient to attenuate differences in precipitation inputs, and that this input has been constant since the 1990s. Two different assumptions for the tritium content of precipitation in the Luquillo Mountains are shown as model 1 and model 2. (B) A comparison of the residence time of water to its weathering derived solutes (the sum of sea-salt corrected base cations + dissolved silicon) based on model 1 (1.1 TU). Waters were sampled on quartz diorite (QD) or on highly weathered QD (QD soil/sap), and on hornfels (HF) or highly weathered HF (HF soil/sap).

Residence times were probably longer in some of the other wells, especially those drilled into deep fracture zones in QD such as the LGW2 wells high in the Icacos watershed. However, inconsistencies were observed among the model ages calculated from different tracers in those wells. For example, the tracers at first yielded conflicting residence times for water in the LGW2B well at the top of the Rio Icacos watershed (see Supplementary Material). By applying a correction with the noble-gas-derived value for excess air in LGW2B (32 cm³/L), the SF₆ age (31 years) and the CFC derived ages (25–33 years) were brought into agreement. Likewise with reasonable assumptions about 1980s tritium concentrations in precipitation in the Luquillo Mountains, the tritium data for this well also yield an age >30 years. In that case, the data are all consistent with water from 15 to 17 m depth in LGW2B that was recharged more than 30 years ago. This long residence time emphasizes that even at high elevations in the watershed, flow paths that start very high on the ridge may require decades of flow time before emerging at the land surface.

Sr isotope data are consistent with this observation of a very long flowpath for water at the top of the Icacos watershed. Specifically, in Figure 5D the Sr isotopic composition of that same well water (0.7067) is very different from the adjacent LGW2C well (0.7049). The Sr isotopic composition of LGW2B is consistent with that of water sampled in HF drainages rather than on QD (Figure 5C,D). The most reasonable explanation of Sr isotopes and chemistry in LGW2B is that there is a contribution from weathering of Ca- (and Sr-) containing minerals in the HF–pointing toward reactions occurring higher on the ridge above the well. Two example minerals that are not reported in the VC but are reported in the HF, wollastonite and scapolite, could be releasing Sr with higher ⁸⁷Sr^/86Sr ratios from the HF. Both minerals are highly susceptible to weathering alteration. If this hypothesis is correct, then some of the HF contact aureole remains even at the location of LGW2B. Indeed, a piece of HF was recovered during drilling of the LGW2 wells (Orlando et al., 2016).

Major elements in weathering reactions

To understand the differences in stream chemistry from river to river, we partitioned water chemistries to specific mineral weathering reactions. Figure 7A,B shows analyte concentrations in the rivers delineated by the different rock types that emphasize the controls on water chemistry by the felsic (QD) versus mafic (HF, VC) mineralogies. The rivers on the more mafic lithologies are characterized by higher concentrations of Ca and Mg compared to felsic lithologies. Consistent with lithological control of chemistry, some of the samples from the upper Mameyes, where the Río La Mina tributary flows on QD, plot near the data for the Río Icacos on QD.

FIGURE 7

FIGURE 7. (A) (B) Stream water draining different lithologies define geochemical groups. The quartz diorite field includes only samples on >98% QD. The hornfels field is delineated by samples draining 100% HF. Several water samples draining some portion of volcaniclastic rocks plot in the HF field. The volcaniclastic field is delineated by the three lowest samples in the Río Mameyes where VC is the dominant lithology. Silicate weathering on HF and VC produces relatively more Ca and Mg than weathering of QD. Silicate weathering of HF produces the lowest proportion of Na and K. (C) (D) Theoretical lines showing the effects of dissolution of different phases (as labelled). Icacos waters are largely dominated by plagioclase dissolution while other phases (augite, chlorite, epidote) contribute to waters on VC and HF.

In the bottom two panels on Figure 7, lines are shown for hypothetical dissolution of each phase as labelled. On the three lithologies, solutes probably derive predominantly from calcic plagioclase, chlorite, clinopyroxene (augite), and to a lesser degree amphibole and epidote where those minerals are present. As shown in Figure 7C, the dominant signal on Mg- and Fe-poor QD is weathering of Río Blanco plutonic An₄₈ plagioclase with some Mg from hornblende or biotite. In contrast, water draining the HF and VC shows the greater influence of dissolution of the more mafic and more soluble An₆₀ plagioclase and augite that are present in those rocks (releasing Ca and Mg but very little Na and K). Thus, QD releases the highest concentration of alkali elements (Na, K) and the lowest concentration of alkaline earth elements (Mg, Ca) relative to dissolved Si. Furthermore, the Río Icacos samples (Figure 7D) demonstrate almost no variation in Ca*/Mg* ratios, again emphasizing the control of feldspar dissolution with only a minor contribution of those elements from hornblende and biotite.

The chlorite line on Figure 7D is consistent with the conclusion that chlorite is the dominant source of Mg across the Sonadora and Espiritu Santo drainages. This conclusion is similar to the conclusion made by other workers that baseflow δ²⁶Mg values at the Río Mameyes Puente Roto (MPR) site on VC reflects dissolution of Mg-rich chlorite (Chapela Lara et al., 2017). Some Mg is also derived from weathering of augitic pyroxene in the Quebrada Sonadora and Espiritu Santo.

Isotopic evidence for weathering reactions

Evidence for weathering reactions is also derived from Sr isotopic data. Consistent with a dominantly weathering-derived source for Sr, the Sr isotopic values measured for VC streams are almost identical to the average values (0.7041) reported for bulk VC bedrock for the two VC formations (Fajardo and Hato Puerco) and the less common basaltic bedrock in those areas (0.7038) (Jones and Kesler, 1980; Frost et al., 1998; Jolly et al., 1998; Chabaux et al., 2013; Porder et al., 2015). The differences between streamwater and bedrock are at least partly caused by atmospheric inputs. As shown in Figure 5A,B, baseflow samples from VC- and HF-dominated catchments can be explained by mixing bedrock- and regolith-derived Sr with throughfall-derived Sr (see also, Supplementary Material).

The data for QD baseflow during the same sampling campaign form a mixing line with a shallower slope than observed for VC/HF (i.e. they are less dominated by precipitation inputs) and yield a predicted high concentration endmember value of ⁸⁷Sr/⁸⁶Sr ≈ 0.704 (Figure 5D). This high concentration endmember which contributes to baseflow is probably a mixture of solutes derived from plagioclase and hornblende.

Precipitation and weathering as anion sources

The chemistry of precipitation in the LEF is strongly affected by marine aerosols, with lesser contributions from African dust and anthropogenic sources (e.g. McDowell et al., 1990; Reid et al., 2003; Stallard, 2012). Between 361 and 1,050 masl, a relationship between elevation and sulfate in precipitation is observed (Gioda et al., 2013.). For this reason, atmospheric inputs of SO₄ to surface systems are greater at the higher elevations of the catchments on QD (Icacos) and HF (Quebrada Sonadora) (McDowell and Asbury, 1994) as compared with the lower-elevation VC (Yi Balan et al., 2014). For example, waters from the QD landscape (from Guaba and the nearby LGW1 well) show SO₄ concentrations and S isotopic compositions that are similar to those of precipitation.

The other source of sulfate to ground and surface waters in the LEF, especially in the Mameyes stream waters on VC at lower elevations, is pyrite oxidation. For example, some sulfate concentrations in samples from the Mameyes (Puente Roto (MPR) site) exceed those of precipitation and have S isotopic compositions consistent with a bedrock source (Yi-Balan et al., 2014). Sulfate attributed to deep oxidation of pyrite is also observed at higher concentrations in the two USFS wells on the VC than in stream waters. The contribution of pyrite-derived sulfate to river waters on the QD is generally difficult to detect and some sulfate concentrations are negative after precipitation correction. However, we have observed pyrite-containing fractures below the water table with evidence of oxidation at depths approaching 40 m in wells drilled in the QD landscape (Figure 4).

Atmospheric deposition also contributes phosphate (Saharan dust) to the LEF (Pett-Ridge, 2009) but most of the trends in dissolved concentrations of PO₄ downstream are best explained as deriving from weathering. This is because these trends correlate with the trends in concentrations of base cations and Si in the watersheds on QD (downstream decrease) and VC (downstream increase) that are shown in Figure 2. Apatite begins to weather as deep as 18 m within the VC but a fraction of the P is retained over the entire regolith (Buss et al., 2017). This retention (and accompanying slow release over meters) may explain the relatively dilute phosphate concentrations in water on the VC as compared to the Icacos. In contrast, apatite weathers in rindlets and saprolite at <5 m depth in QD (Buss et al., 2010). This relatively shallow source of P in the QD may explain why the highest P concentrations on QD (at higher elevations) are approximately twice that of the highest concentrations on VC (at lower elevations), even though the P content of VC bedrock and soils are both ∼2X that of QD (Mage and Porder, 2013).

Identifying sources of solutes on QD

To explore what controls the high-elevation trends in QD water, we used concentration-ratio plots. As shown in Figure 8A, lines were calculated to show different hypothetical contributions of plagioclase and hornblende weathering to water samples. For example, water sampled on QD high in the watershed (in borehole LGW2C groundwater and Río Mameyes stream water at LM1) plots on the 99% plagioclase +1% hornblende dissolution line. Borehole LGW2C was drilled through multiple rindletted corestones within a vertical fracture zone (Orlando et al., 2016; Comas et al., 2019). Thus, the groundwaters sampled at 17.8 m depth in LGW2C (as well as water from ∼38 m depth in borehole EP1 drilled through HF and QD) have interacted with multiple corestones, rindlets, and fractures (Figure 3). EP1 and LGW2C samples also have a higher component of alkaline earth cations, as expected since both boreholes are thought to intersect QD rock with plagioclase and hornblende. This QD-specific rock water chemistry is labelled generically as “rock water”. The rock endmember was estimated by projecting a linear fit from the saprolite water through the Icacos water samples (r² = 0.77) to the 99:1 plagioclase:hornblende line.

FIGURE 8

FIGURE 8. Mixing of rock water and saprolite water in the Río Icacos drainage. (A) the geochemistry of stream water on QD (Río Icacos) forms a linear array on this element ratio plot, consistent with two-endmember mixing. Water dominated by rock- and saprolite-derived solutes comprise the two endmember water compositions that define QD stream water. Data for groundwater from wells are also shown (LGW1, EP1, LGW2B, LGW2C). Shallow groundwater from the LGW1 well is compositionally similar to saprolite water. The coordinates of the sample for LGW2B are offscale as indicated in parentheses. (B) Plot showing the decrease in fraction of rock water (F_rw) at stream sites with increasing drainage area for samples along the Icacos mainstem. See text for description of calculation of F_rw. See Table 1 for location of sample site RIS3 which plots off the trend.

The water from LGW1 is consistent with water that not only is interacting with plagioclase, but also contains some Mg from biotite—a reaction that occurs almost exclusively in the saprolite (White et al., 1998). In contrast to LGW2C which is drilled within one of the lineaments crossing the Icacos and thus intersects a vertical fracture zone with several corestones, borehole LGW1 is drilled on an interfluve and intersects saprolite and only one corestone. The water sampled from LGW1 has thus infiltrated through 2–8 m of saprolite (White et al., 1998) before reaching the saprolite/bedrock interface where it subsequently interacts with only one ∼30 cm thick rindlet (shown in Figure 4; see also Supplementary Figure S3). The chemistry in LGW1 is thus labelled “saprolite water” (and is based on Table 5 of White et al. (1998)) because it reflects a larger contribution from weathering of saprolite (quartz + biotite + kaolinite) than rindlet. The rest of the measured values of Río Icacos stream water (open symbols, data from Table 2) plot on a mixing line in Figure 8 between these endmembers of “saprolite water” and “rock water”.

TABLE 5

TABLE 5. Atmospheric tracer data for surface and groundwater samples.

Strontium isotope ratios for water in Río Icacos are also consistent with solute source contributions from saprolite and rock water sources (Pett-Ridge et al., 2009a; Pett-Ridge et al., 2009b; Chabaux et al., 2013). For example, using molar ratios in Figure 8A, the fraction of rock water (F_rw) in Río Icacos water was calculated to range from 33 to 77%, with F_rw decreasing downstream (Figure 8B). Downstream at the Río Icacos stream gauge (site RI in Table 2) we calculated F_rw in baseflow to be 0.42. Together, the major cation and Sr isotope data suggest that baseflow solute flux in streams draining QD is derived 20–77% from rock water along the study reach; downstream at the Río Icacos stream gauge the range is 33–50% (Figure 8B).

The decreasing fraction of rock water in the Icacos mainstem on QD was explored in 1 day of sampling along the stream reach on 25 February 2014 when discharge at the Río Icacos gauge was 244 L/s, slightly less than the median discharge over the period of record. These data were used to calculate watershed-scale estimates of solute fluxes. These results confirm that groundwaters draining from high ridges of the Icacos watershed–ridges where we have observed stacks of rindletted corestones—yield a disproportionately high solute flux per unit area because they have a higher fraction of rock water (Figure 8). Downstream, we attribute the rapid decrease in the elemental yield normalized by watershed area to the increasing contribution of saprolite water which begins to dominate as watershed size increases and much of the incoming rainwater passes only through the saprolite-covered Icacos floodplain (Figures 8B, 9). Ultimately, the elemental solute flux only increases again toward the Río Icacos stream gauge (see also, Shanley et al., 2011). This downstream increase near the gauge is attributed to rock water draining the high ridges and reaching the stream along deeper flowpaths that may follow the upper interface of unweathered, massive QD: at the Río Icacos stream gauge, these flow paths may account for 30–50% of baseflow discharge and the majority of weathering-derived solutes.

FIGURE 9

FIGURE 9. Mass fluxes plotted for the mainstem of the Icacos as a function of upstream drainage area for Río Icacos synoptic sampling on 25 February 2014. (A) Dissolved Si and base cation (Σbc*) fluxes from discharge measurements paired with stream samples, supporting data can be found in Supplemental Table S6. (B) Total weathering flux (Si + Σbc*) for the February 25, 2014 synoptic compared to indendent long term flux estimates at USGS gauges in the Luquillo Mountains, see Supplementary Table S7 for details of those estimates. The flux at the site highest in the watershed is equivalent or higher than the flux at the Río Icacos gauge. Fluxes per unit area are higher than the baseflow estimates from Shanley et al. (2011) assuming a baseflow runoff of 0.10–0.20 mm/h at Río Icacos stream gauge (shown in shaded bars on both panels).

In summary, the high weathering fluxes near the upper part of the QD stream (Icacos) identify the steep ridges surrounding the watershed as weathering hotspots, likely driven by meteoric water rapidly infiltrating to, reacting with, and draining from rindletted corestones on bedrock (Bhatt and McDowell 2007). The intermediate reaches of Río Icacos are diluted by shallow-flowing saprolite water which produces lower weathering-derived solute fluxes. Above the Río Icacos stream gauge, the estimated F_rw reaches its minimum along the study reach. Further downstream at the gauge, solute fluxes again increase and are likely the result of rock water derived from longer flowpaths which may infiltrate to greater depth along fracture zones and then return to the mainstem flowing along the upper interface of massive QD (Orlando et al., 2016; Hynek et al., 2017; Comas et al., 2019).

Weathering fluxes on the three lithologies

We can also pair the average observed stream chemistry with median discharge values derived from decades of streamflow data to explore watershed scale solute fluxes (Supplementary Table S7, Figure 9). This simple approach to calculate fluxes agrees well with previous approaches based on more in-depth estimates (e.g., Murphy and Stallard, 2012), but also has the additional benefit of allowing direct comparisons among the three lithologies.

The total solute flux from the HF-dominated watershed (Sonadora) is lower than for the other two lithologies (Supplementary Table S7). In addition, stream discharge data indicate that the groundwater component of stream discharge is substantially lower for watersheds containing high proportions of HF (Supplementary Table S7). As determined by the baseflow index (BFI, the proportion of mean annual flow that is from groundwater, tabulated in Supplementary Table S7, see also Supplementary Figure S5), Quebrada Sonadora has the lowest proportion of groundwater of any gauged stream in the Luquillo Mountains. Quebrada Sonadora shares this characteristic with other streams draining significant portions of HF; conversely, watersheds underlain by nearly 100% QD have the highest BFI, followed by Río Mameyes at Puente Roto and the other VC-dominated watersheds.

Solute export rates for the Quebrada Sonadora are 25–50% of those observed for the other study watersheds, and the contribution of groundwater to these fluxes is minimal (Supplementary Table S7, Supplementary Figure S5). In addition, shallow wells and seeps on HF ridges have among the lowest tritium concentrations measured in our study, supporting infiltration through and interaction with a relatively impermeable regolith (Table 5, see also, Supplementary Material). These sites also have relatively low solute concentrations (Figure 2), suggesting that HF regolith reacts slowly with infiltrating water. Despite the relatively high solubility of some of the minerals in the HF, including Ca-rich anorthite and augite, the slow chemical weathering in HF appears to be the result of the characteristically low rate of infiltration associated with the low porosity and low fracture density. All of these characteristics inhibit chemical reactions and formation of deep regolith. These characteristics are summarized schematically in Figure 10.

FIGURE 10

FIGURE 10. Summary schematic of the impact of lithology on stream chemistry. Top row of panels shows the relative landscape surface area as a function of elevation for each lithology. This documents how HF dominates the ridges, QD the mid-elevations, and VC the lowest elevations. The bottom row of colored panels shows conceptual models of regolith structure and hydrology on the three study lithologies. The shaded inset on each regolith profile represents the feldspar weathering front, i.e. the % of feldspar left at each depth from land surface to unweathered parent. Regolith depth decreases from VC to QD to HF and corestones that form are angular (VC, HF) or rounded (QD). The bottom summary shows how primary porosity (porosity of unweathered bedrock), secondary porosity (porosity of weathered bedrock), and drainable porosity (connected porosity) differs for the three lithologies.

In contrast to the Sonadora (HF), the Río Mameyes data imply a strong component of groundwater in the stream at Río Mameyes Puente Roto (MPR) where the river incises bedrock. Apparently, in this section of the stream, groundwater that has flowed in the subsurface beneath the stream in the upper parts of the watershed returns to the stream, increasing the concentrations of bedrock-derived solutes (e.g., Si, Mg, Ca, Sr, SO₄, PO₄). For example, the Mameyes tributaries in the Bisley watersheds (upstream of the MPR site), have been shown to be perched streams flowing on deeply developed regolith (Schellekens et al., 2004; Buss et al., 2013).

Fractures

The picture that emerges from this work is that the differences in longitudinal profiles of stream chemistry are largely caused by differences in mineralogy, fracturing, and the related impacts on groundwater flowpaths within the watersheds. For example, hornfels facies rocks are generally very hard rocks with crystals fitting closely together that do not easily allow the rock to fracture. Thus the river on HF (Sonadora) reflects little influx of groundwater.

On the other hand, sedimentary rocks such as the VC are generally more highly fractured than massive hornfels-facies rock. This is likely because the spacing of vertical fracture joints is generally observed to be smaller in rocks with thinner beds–such as the VC—as compared to rocks with thicker beds or no beds (Gross et al., 1995). This explains why many fractures are observed in the layered sedimentary VC bedrock both in outcrop and in boreholes (Buss et al., 2013). The presence of a pre-existing fracture network in the VC could explain the apparent deep flow of groundwater and return of that water to the Mameyes downstream. Such a network and pattern of flow also can explain why water chemistry moving downstream moves towards that of the deep groundwater on the VC in the Río Mameyes (Figure 2).

In contrast, fractures are only observed around weathered corestones and in vertical fracture zones with stacked corestones in the QD: the rest of the rock beneath or away from the corestones is massive (Figure 4). Groundwater flow through the QD is therefore restricted to the surficial zones of weathering and corestone formation. The flow apparently proceeds through the overlying saprolite and then vertically down into the deep fracture zones (filled with stacks of corestones) or horizontally along the layer comprised of a single corestone between the deep vertical fracture zones (Hynek et al., 2017; Comas et al., 2019). Wherever the water flows through the rindlets around the corestones, it weathers plagioclase and hornblende, picking up rock-derived cations such as Na and Ca. As shown in Figure 8B, the fraction of rock water decreases downstream as the fraction of saprolite water increases in the lower part of the Icacos watershed, explaining the pattern shown in Figure 2. From these observations we infer that eventually, some or almost all of the deeper rock water is pushed up and into the mainstem stream near the Icacos gauge, since the flux values increase at that point (Figure 2, 9). This increase just at the gauge may also be related to the decreasing spacing between the deep vertical fracture zones that has been observed closer to the knickzone (Comas et al., 2019).

Interactions among lithology, chemical denudation and topography

Variation in lithology in the Luquillo Mountains drives changes in porosity/permeability that result in large differences in flow paths. The coloring in Figure 10 schematically indicates the % of the feldspar remaining in the rock compared to the initial abundance in parent. The Figure emphasizes that regolith is thinnest on the HF, thicker on the QD, and thickest on the VC. For example, plagioclase was completely depleted from the entire layer of saprolite on HF but some plagioclase was still present at the land surface because blocky angular bedrock fragments littered the land surface on the ridges (Lebedeva and Brantley, 2017). Thus, the regolith profile on HF is considered to be incompletely developed with respect to primary minerals (Brantley and White, 2009). In an incompletely developed profile, the weathering mineral is not completely weathered away in the regolith even at the land surface. We infer that this lithology defines the highest topography in the LEF because of its low bedrock porosity, its limited tendency to develop secondary porosity, and its resistance to fracturing and physical erosion (Gerrard, 1988; Behrens et al., 2015). These attributes explain why the HF is poorly drained and why bogs are only found in the Luquillo Mountains on ridges underlain by HF (Ogle, 1970). Even though it contains highly weatherable minerals such as anorthite-rich feldspar and augite, the weathering flux in the Sonadora (on HF) is very low compared to both the VC and QD (Supplementary Table S7) as well as other Caribbean catchments (McDowell et al., 1995).

In contrast, on the QD, ample evidence documents that the underlying bedrock is unfractured but that once weathering commences with oxidation of biotite, it induces fracturing (Fletcher et al., 2006; Buss et al., 2008). The Si flux out of the watershed is highest at the top of the watershed where the regolith is thickest on average (Figure 9). Ridgetops on the QD are always saprolite-mantled and thus are completely developed with respect to feldspar: no plagioclase remains in the saprolite at the land surface and no corestones are observed at ridgetop surfaces. Given the relatively unfractured bedrock, meteoric fluid only interacts with feldspar in rindlet zones (White et al., 1998). Dissolved Si concentrations in the streams decrease down gradient because the fraction of rock water is higher near the ridges. At lower elevations, this water is increasingly diluted with saprolite water toward the knickpoint (see Figure 8) because the volume of permeable material in the watershed through which groundwater flows becomes increasingly dominated by saprolite as opposed to the more limited locations of rindletted corestones. Finally, as shown in Figure 9 the elevated fluxes are observed not only high in the watershed, but also toward the outlet at the knickzone: this suggests that groundwater with high solute concentrations from deep long-residence time flowpaths is discharging to the stream channel in the reaches just above the knickzone.

In contrast to the QD and HF, the VC bedrock has a relatively high density of pre-existing fractures due to its sedimentary layering. Weathering-induced fracturing generally does not occur in the VC and the corestones that form are angular instead of spheroidal (Buss et al., 2013). Significant meteoric fluid infiltrates into the fractured rock, and weathering continues beneath the water table, explaining the increasing Si concentration in increasingly deeper water on the VC. Thick regolith on VC and increasing solute concentrations downstream on the Mameyes are all consistent with solute transport during weathering by deeply infiltrating water.

These observations are generally consistent with numerical calculations for weathering of model rocks. In particular, the calculations show that regolith that develops on a rock is thinner when advection fluxes are limited, when other factors are held constant (Bazilevskaya et al., 2013) and also that the regolith thickness increases with increasing fracture density (Lebedeva and Brantley, 2017). This is because solute transport across the weathering reaction front in low-infiltration systems occurs by diffusion rather than advection, and this limits the thickness of regolith that can develop during the residence time that the material weathers. This is schematically indicated in Figure 10: on the HF, fracture density is low, weatherable rock fragments are angular, and the regolith thickness that develops is thin. In contrast, on the other two lithologies, fractures are more prevalent and regolith is correspondingly thicker.

We also argue that above the knickpoint, little to no meteoric fluid passes through the unweathered HF and QD bedrock because of the lack of fractures at depth. For water that does infiltrate HF, the hornfels depletes infiltrating porewaters of oxygen so that when the waters reach the QD at depth, weathering-induced fracturing is not promoted. Only when the HF is removed–for example in the bowl carved by the Río Icacos -- or when most of the iron-containing minerals in the overlying HF are oxidized–such as observed beneath about 7 m depth in borehole EP1 -- do oxygenated waters reach the underlying QD where they can promote oxidation and accelerate spheroidal weathering.

Conclusions

Surface water chemistry was investigated in tandem with groundwater chemistry and regolith architecture in the Luquillo Mountains, Puerto Rico to understand how lithology affects stream chemistry and subsurface water flowpaths. In the Luquillo Mountains, weathering has exhumed a metamorphic aureole, permitting study of how weathering and erosion are affected by lithology.

Three watersheds, each typifying a different lithology, were observed to exhibit different longitudinal patterns of stream chemistry. The hornfels facies contact metamorphic aureole defines the rocky, steep, ridgetops of the mountains. Where the ridges on hornfels have weathered through to the underlying QD, meters of QD- and HF-based saprolite have developed, but the land surface is nonetheless littered with blocky, angular, unweathered hornfels rocks of tens of centimeters in dimension. Thus, the dominant mineral, plagioclase, has not all weathered away, consistent with weathering of a rock that is mostly unfractured (massive) and thus resistant to chemical weathering and physical erosion, and unlikely to act as a permeable porous medium that can support significant through-flow of meteoric water. For these reasons, the Quebrada Sonadora, the watershed in the Luquillo Mountains that has developed on the highest proportion of HF, exhibits the lowest solute concentrations (when compared to waters on QD or VC). The solute concentrations also remain relatively constant downstream as the watershed area increases because little groundwater enters the Sonadora: the hornfels lithology does not host significant groundwater storage or flow.

On the QD in the Río Icacos watershed, solute concentrations decrease downstream as the stream flows away from the high HF ridges towards a large knickzone. This downstream decrease documents the increasing contribution of water flow through saprolite as compared to bedrock. The fraction of stream water that has permeated deeply into the QD is highest in upstream reaches where HF-mantled ridges maintain topography and thick stacks of rindletted corestones of QD have developed. Such corestone development is especially prevalent in deep vertical fracture zones: between these zones, water is presumed to pass through saprolite and along the weathered layer along interfluves which are characterized not by stacks of corestones, but rather by a layer comprised of a single rindletted corestone. The most Si-rich waters apparently emanate from the deeply weathered fracture zones and this water dominates the upper Río Icacos. Water that only flows through shallow, saprolite-dominated regolith increasingly dominates stream chemistry lower in the watershed. In the lowest reaches of Río Icacos just above the prominent knickpoint, water with high solute concentrations re-enters the stream, perhaps because the spacing between the deep vertical fracture zones decreases near the knickpoint (Comas et al., 2019).

Finally, Río Mameyes, the river on VC, is characterized by solute concentrations that increase downstream with watershed area. This pattern reflects the highly fractured nature of the VC. The fracture network allows deep infiltration of groundwater that moves beneath the surface of the upper tributaries, only joining the mainstem of the Mameyes downstream. As groundwater rejoins the Mameyes, solute concentrations increase.

The topography of the LEF is highly affected by its subsurface architecture. If all three lithologies were able to attain a geomorphological steady state such that the total denudation rates for each lithology were equal, one would predict steep slopes on the HF to allow denudation of the rock in the face of its high fracture toughness and lowered weathering extent. Likewise, one might predict lower slopes on the more easily fractured, eroded, and weathered VC. However, because so much of the precipitation that falls on VC infiltrates the rock and enters groundwater, moderately steep slopes must develop so that insoluble minerals that remain after weathering can be eroded away. Finally, one might predict that denudation of the QD is both chemically fast and physically easy because of its soluble minerals and the rock’s high tendency to fracture and erode. These characteristics explain the relatively shallow slopes observed on the QD landscape. The lithology thus helps explain the relative relief of the Luquillo Mountains, namely very steep HF hillslopes at the highest elevations, relatively steep VC hillslopes, and relatively gentle average hillslopes on the QD, as well as the stream chemistries in the forest.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

SH, MB and JO completed field and laboratory work. SH, WM and SB wrote the manuscript. WM and SB supervised the project.

Funding

This material is based on work supported by the U.S. National Science Foundation (NSF) Luquillo Critical Zone Observatory (NSF-LCZO, Grant EAR 1331841 to WM; University of New Hampshire), the Luquillo LTER program supported by NSF (DEB 1239764, DEB 1546686, DEB 1831592), and the USGS Global Change and National Research Programs (NRP).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank Heather Buss (University of Bristol), Miguel Leon and Rafael Jiménez (University of Pennsylvania); Michelle Daly, Jordan Macy, Jody Potter, Jeff Merriam, Geoff Schwaner, and Josh Brown (University of New Hampshire); Carole Johnson, Manuel Rosario Torres, and Martha Scholl (USGS), Carlos R. Colón Rivera, and Christina Henderson (USFS); Wil Mace, Alan Rigby, and Kip Solomon (University of Utah, Dissolved and Noble Gas Lab); and Suelos Inc. (San Juan, PR). Chemical analyses made at Penn State were completed at the Laboratory for Isotopes and Metals in the Environment (LIME). We also thank (and miss) Fred Scatena (1954–2013) whose support and encouragement in the LCZO made this work possible. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2022.779459/full#supplementary-material

References

Aeschbach-Hertig, W., Beyerle, U., Holocher, J., Peeters, F., and Kipfer, R. (2002). Excess Air in Groundwater as a Potential Indicator of Past Environmental Changes, Study of Environmental Change Using Isotope Techniques. Vienna, Austria: IAEA, 174–183.

Google Scholar

Aeschbach-Hertig, W., El-Gamal, H., Wieser, M., and Palcsu, L. (2008). Modeling Excess Air and Degassing by Equilibrium Partitioning with a Gas Phase. Water Resour. Res. 44 (8), W08449. doi:10.1029/2007wr006454

CrossRef Full Text | Google Scholar

Bailey, S. W., McGuire, K. J., Ross, D. S., Green, M. B., and Fraser, O. L. (2019). Mineral Weathering and Podzolization Control Acid Neutralization and Streamwater Chemistry Gradients in Upland Glaciated Catchments, Northeastern United States. Front. Earth Sci. 7 (63). doi:10.3389/feart.2019.00063

CrossRef Full Text | Google Scholar

Banks, E. W., Simmons, C. T., Love, A. J., Cranswick, R., Werner, A. D., Bestland, E. A., et al. (2009). Fractured Bedrock and Saprolite Hydrogeologic Controls on Groundwater/surface-Water Interaction: A Conceptual Model (Australia). Hydrogeol J. 17, 1969–1989. doi:10.1007/s10040-009-0490-7

CrossRef Full Text | Google Scholar

Bazilevskaya, E., Lebedeva, M., Pavich, M., Rother, G., Parkinson, D. Y., Cole, D., et al. (2013). Where Fast Weathering Creates Thin Regolith and Slow Weathering Creates Thick Regolith. Earth Surf. Process. Landforms 38 (8), 847–858. doi:10.1002/esp.3369

CrossRef Full Text | Google Scholar

Behrens, R., Bouchez, J., Schuessler, J. A., Dultz, S., Hewawasam, T., and von Blanckenburg, F. (2015). Mineralogical Transformations Set Slow Weathering Rates in Low-Porosity Metamorphic Bedrock on Mountain Slopes in a Tropical Climate. Chem. Geology. 411, 283–298. doi:10.1016/j.chemgeo.2015.07.008

CrossRef Full Text | Google Scholar

Bhatt, M. P., and McDowell, W. H. (2007). Controls on Major Solutes within the Drainage Network of a Rapidly Weathering Tropical Watershed. Water Resour. Res. 43 (W11402), 1–9. doi:10.1029/2007wr005915

PubMed Abstract | CrossRef Full Text | Google Scholar

Brantley, S. L., Buss, H., Lebedeva, M., Fletcher, R. C., and Ma, L. (2011). Investigating the Complex Interface where Bedrock Transforms to Regolith. Appl. Geochem. 26, S12–S15. doi:10.1016/j.apgeochem.2011.1003.101710.1016/j.apgeochem.2011.03.017

CrossRef Full Text | Google Scholar

Brantley, S. L., Lebedeva, M. I., Balashov, V. N., Singha, K., Sullivan, P. L., and Stinchcomb, G. (2017). Toward a Conceptual Model Relating Chemical Reaction Fronts to Water Flow Paths in hills. Geomorphology 277, 100–117. doi:10.1016/j.geomorph.2016.09.027

CrossRef Full Text | Google Scholar

Brantley, S. L., and Lebedeva, M. I. (2021). Relating Land Surface, Water Table, and Weathering Fronts with a Conceptual Valve Model for Headwater Catchments. Hydrological Process. n/a (n/a), e14010. doi:10.1002/hyp.14010

CrossRef Full Text | Google Scholar

Brantley, S. L., and White, A. F. (2009). “Approaches to Modeling Weathered Regolith,” in Thermodynamics and Kinetics of Water-Rock Interaction. Reviews in Mineralogy and Geochemistry. Editors E. Oelkers, and J. Schott, 435–484. doi:10.1515/9781501508462-012

CrossRef Full Text | Google Scholar

Braun, J.-J., Descloitres, M., Riotte, J., Fleury, S., Barbiéro, L., Boeglin, J.-L., et al. (2009). Regolith Mass Balance Inferred from Combined Mineralogical, Geochemical and Geophysical Studies: Mule Hole Gneissic Watershed, South India. Geochimica et Cosmochimica Acta 73, 935–961. doi:10.1016/j.gca.2008.11.013

CrossRef Full Text | Google Scholar

Brocard, G. Y., Willenbring, J. K., Scatena, F. N., and Johnson, A. H. (2015). Effects of a Tectonically-Triggered Wave of Incision on Riverine Exports and Soil Mineralogy in the Luquillo Mountains of Puerto Rico. Appl. Geochem. 63, 586–598. doi:10.1016/j.apgeochem.2015.04.001

CrossRef Full Text | Google Scholar

Brown, E. T., Stallard, R. F., Larsen, M. C., Raisbeck, G. M., and Yiou, F. (1995). Denudation Rates Determined from the Accumulation of In Situ-produced 10Be in the Luquillo Experimental forest, Puerto Rico. Earth Planet. Sci. Lett. 129, 193–202. doi:10.1016/0012-821x(94)00249-x

CrossRef Full Text | Google Scholar

Buss, H. L., Brantley, S. L., Scatena, F. N., Bazilievskaya, E. A., Blum, A., Schulz, M., et al. (2013). Probing the Deep Critical Zone beneath the Luquillo Experimental Forest, Puerto Rico. Earth Surf. Process. Landforms 38 (10), 1170–1186. doi:10.1002/esp.3409

CrossRef Full Text | Google Scholar

Buss, H. L., Chapela Lara, M., Moore, O. W., Kurtz, A. C., Schulz, M. S., and White, A. F. (2017). Lithological Influences on Contemporary and Long-Term Regolith Weathering at the Luquillo Critical Zone Observatory. Geochimica et Cosmochimica Acta 196, 224–251. doi:10.1016/j.gca.2016.09.038

CrossRef Full Text | Google Scholar

Buss, H. L., Sak, P. B., Webb, S. M., and Brantley, S. L. (2008). Weathering of the Rio Blanco Quartz Diorite, Luquillo Mountains, Puerto Rico: Coupling Oxidation, Dissolution, and Fracturing. Geochimica et Cosmochimica Acta 72, 4488–4507. doi:10.1016/j.gca.2008.06.020

CrossRef Full Text | Google Scholar

Chabaux, F., Blaes, E., Stille, P., di Chiara Roupert, R., Pelt, E., Dosseto, A., et al. (2013). Regolith Formation Rate from U-Series Nuclides: Implications from the Study of a Spheroidal Weathering Profile in the Rio Icacos Watershed (Puerto Rico). Geochimica et Cosmochimica Acta 100, 73–95. doi:10.1016/j.gca.2012.09.037

CrossRef Full Text | Google Scholar

Chapela Lara, M., Buss, H. L., Pogge von Strandmann, P. A. E., Schuessler, J. A., and Moore, O. W. (2017). The Influence of Critical Zone Processes on the Mg Isotope Budget in a Tropical, Highly Weathered Andesitic Catchment. Geochimica et Cosmochimica Acta 202, 77–100. doi:10.1016/j.gca.2016.12.032

CrossRef Full Text | Google Scholar

Chesson, L. A., Tipple, B. J., Mackey, G. N., Hynek, S. A., Fernandez, D. P., and Ehleringer, J. R. (2012). Strontium Isotopes in Tap Water from Coterminous USA. Ecosphere 3 (7). doi:10.1890/ES1812-00122.0012110.1890/es12-00122.1

CrossRef Full Text | Google Scholar

Clarke, W. B., Jenkins, W. J., and Top, Z. (1976). Determination of Tritium by Mass Spectrometric Measurement of 3He. Int. J. Appl. Radiat. Isotopes 27, 515–522. doi:10.1016/0020-708x(76)90082-x

CrossRef Full Text | Google Scholar

Comas, X., Wright, W., Hynek, S. A., Fletcher, R. C., and Brantley, S. L. (2019). Understanding Fracture Distribution and its Relation to Knickpoint Evolution in the Rio Icacos Watershed (Luquillo Critical Zone Observatory, Puerto Rico) Using Landscape‐scale Hydrogeophysics. Earth Surf. Process. Landforms 44 (4), 877–885. doi:10.1002/esp.4540

CrossRef Full Text | Google Scholar

Cox, D. P., Marvin, R. F., M’Gonigle, J. W., McIntyre, D. H., and Rogers, C. L. (1977). Potassium-argon Geochronology of Some Metamorphic, Igneous, and Hydrothermal Events in Puerto Rico and the Virgin Islands. J. Res. U.S. Geol. Surv. 5, 689–703.

Google Scholar

Dixon, J. L., and von Blanckenburg, F. (2012). Soils as Pacemakers and Limiters of Global Silicate Weathering. Comptes Rendus Geosci. 344 (11-12), 597–609. doi:10.1016/j.crte.2012.10.012

CrossRef Full Text | Google Scholar

Dosseto, A., Buss, H. L., and Suresh, P. O. (2012). Rapid Regolith Formation over Volcanic Bedrock and Implications for Landscape Evolution. Earth Planet. Sci. Lett. 337-338, 47–55. doi:10.1016/j.epsl.2012.1005.100810.1016/j.epsl.2012.05.008

CrossRef Full Text | Google Scholar

Fletcher, R. C., and Brantley, S. L. (2010). Reduction of Bedrock Blocks as Corestones in the Weathering Profile: Observations and Model. Am. J. Sci. 310, 131–164. doi:10.2475/03.2010.01

CrossRef Full Text | Google Scholar

Fletcher, R. C., Buss, H. L., and Brantley, S. L. (2006). A Spheroidal Weathering Model Coupling Porewater Chemistry to Soil Thicknesses during Steady-State Denudation. Earth Planet. Sci. Lett. 244 (1-2), 444–457. doi:10.1016/j.epsl.2006.01.055

CrossRef Full Text | Google Scholar

Frost, C. D., Schellekens, J. H., and Smith, A. L. (1998). “Nd, Sr, and Pb Isotopic Characterization of Cretaceous and Paleogene Volcanic and Plutonic Island Arc Rocks from Puerto Rico,” in Tectonics and Geochemistry of the Northeastern Caribbean, G.S.A. Special Paper 322. Editors E. G. Lidiak, and D. K. Larue (Boulder, CO (USA): Geological Society of America). doi:10.1130/0-8137-2322-1.123

CrossRef Full Text | Google Scholar

Gardner, P., and Solomon, D. K. (2009). An Advanced Passive Diffusion Sampler for the Determination of Dissolved Gas Concentrations. Water Resour. Res. 45 (6), n/a–n. doi:10.1029/2008wr007399

CrossRef Full Text | Google Scholar

Gerrard, J. (1988). Rocks and Landforms. London: Unwin Hyman.

Google Scholar

Gioda, A., Mayol-Bracero, O. L., Scatena, F. N., Weathers, K. C., Mateus, V. L., and McDowell, W. H. (2013). Chemical Constituents in Clouds and Rainwater in the Puerto Rican Rainforest: Potential Sources and Seasonal Drivers. Atmos. Environ. 68, 208–220. doi:10.1016/j.atmosenv.2012.11.017

CrossRef Full Text | Google Scholar

Gross, M. R., Fischer, M. P., Engelder, T., and Greenfield, R. J. (1995). Factors Controlling Joint Spacing in Interbedded Sedimentary Rocks: Integrating Numerical Models with Field Observations from the Monterey Formation, USA. Geol. Soc. Lond. Spec. Publications 92, 215–233. doi:10.1144/gsl.sp.1995.092.01.12

CrossRef Full Text | Google Scholar

Gu, X., Kim, H., Hynek, S., Thompson, A., and Brantley, S. L. (2021). Subsurface Particle Transport Shapes the Deep Critical Zone in a Granitoid Watershed. Geochem. Persp. Let. 19, 13–18. doi:10.7185/geochemlet.2127

CrossRef Full Text | Google Scholar

Hayes, N. R., Buss, H. L., Moore, O. W., Kram, P., and Pancost, R. D. (2020). Controls on Granitic Weathering Fronts in Contrasting Climates. Chem. Geology. 535, 19. doi:10.1016/j.chemgeo.2019.119450

CrossRef Full Text | Google Scholar

Hynek, S., Comas, X., and Brantley, S. L. (2017). “The Effect of Fractures on Weathering of Igneous and Volcaniclastic Sedimentary Rocks in the Puerto Rican Tropical Rain forest,” in 15th Water-Rock Interaction International Symposium, Wri-15. Editors J. M. Marques, and A. Chambel (Procedia Earth and Planetary Science), 17, 972–975. doi:10.1016/j.proeps.2017.01.001Proced. Earth Planet. Sci.

CrossRef Full Text | Google Scholar

Jolly, W. T., Lidiak, E. G., Schellekens, J. H., and Santos, H. (1998). “Volcanism, Tectonics, and Stratigraphic Correlations in Puerto Rico,” in Tectonics and Geochemistry of the Northeastern Caribbean, GSA Special Paper 322. Editors E. G. Lidiak, and D. K. Larue (Boulder, CO (USA): Geological Society of America), 1–34. doi:10.1130/0-8137-2322-1.1

CrossRef Full Text | Google Scholar

Jones, L. M., and Kesler, S. E. (1980). Strontium Isotopic Geochemistry of Intrusive Rocks, Puerto Rico, Greater Antilles. Earth Planet. Sci. Lett. 50 (1), 219–224. doi:10.1016/0012-821x(80)90133-8

CrossRef Full Text | Google Scholar

Larsen, M. C., Torres-Sánchez, A. J., and Concepción, I. M. (1999). Slopewash, Surface Runoff and fine-litter Transport in forest and Landslide Scars in Humid-Tropical Steeplands, Luquillo Experimental forest, Puerto Rico. Earth Surf. Process. Landforms 24, 481–502. doi:10.1002/(sici)1096-9837(199906)24:6<481::aid-esp967>3.0.co;2-g

CrossRef Full Text | Google Scholar

Larsen, M. C. (1997). Tropical Geomorphology and Geomorphic Work: A Study of Geomorphic Processes and Sediment and Water Budgets in Montane Humid-Tropical Forested and Developed Watersheds. Ph.D. (Boulder, CO: University of Colorado).

Lawrence, G. B., and Driscoll, C. T. (1990). Longitudinal Patterns of Concentration Discharge Relationships in Stream Water Draining the Hubbard Brook Experimental Forest, New Hampshire. J. Hydrol. 116 (1-4), 147–165. doi:10.1016/0022-1694(90)90120-m

CrossRef Full Text | Google Scholar

Lebedeva, M. I., and Brantley, S. L. (2013). Exploring Geochemical Controls on Weathering and Erosion of Convex Hillslopes: beyond the Empirical Regolith Production Function. Earth Surf. Process. Landforms 38, 1793–1807. doi:10.1002/esp.3424

CrossRef Full Text | Google Scholar

Lebedeva, M. I., and Brantley, S. L. (2017). Weathering and Erosion of Fractured Bedrock Systems. Earth Surf. Process. Landforms 42, 2090–2108. doi:10.1002/esp.4177

CrossRef Full Text | Google Scholar

Li, L., Maher, K., Navarre-Sitchler, A., Druhan, J., Meile, C., Lawrence, C., et al. (2017). Expanding the Role of Reactive Transport Models in Critical Zone Processes. Earth-Science Rev. 165, 280–301. doi:10.1016/j.earscirev.2016.09.001

CrossRef Full Text | Google Scholar

Lucas, L. L., and Unterweger, M. P. (2000). Comprehensive Review and Critical Evaluation of the Half-Life of Tritium. J. Res. Natl. Inst. Stand. Technol. 105 (4), 541–549. doi:10.6028/jres.105.043

PubMed Abstract | CrossRef Full Text | Google Scholar

MacQuarrie, K. T. B., and Mayer, K. U. (2005). Reactive Transport Modeling in Fractured Rock: A State-Of-The-Science Review. Earth-sci. Rev. 72 (3-4), 189–227. doi:10.1016/j.earscirev.2005.07.003

CrossRef Full Text | Google Scholar

Mage, S. M., and Porder, S. (2013). Parent Material and Topography Determine Soil Phosphorus Status in the Luquillo Mountains of Puerto Rico. Ecosystems 16 (2), 284–294. doi:10.1007/s10021-012-9612-5

CrossRef Full Text | Google Scholar

McDowell, W. H., McDowell, W. G., Potter, J. D., and Ramírez, A. (2019). Nutrient export and Elemental Stoichiometry in an Urban Tropical River. Ecol. Appl. 29, e01839. doi:10.1002/eap.1839

PubMed Abstract | CrossRef Full Text | Google Scholar

McDowell, W. H., and Asbury, C. E. (1994). Export of Carbon, Nitrogen, and Major Ions from Three Tropical Montane Watersheds. Limnol. Oceanogr. 39 (1), 111–125. doi:10.4319/lo.1994.39.1.0111

CrossRef Full Text | Google Scholar

McDowell, W. H., Leon, M. C., Shattuck, M. D., Potter, J. D., Heartsill‐Scalley, T., González, G., et al. (2021). Luquillo Experimental Forest: Catchment Science in the Montane Tropics. Hydrological Process. 35. doi:10.1002/hyp.14146

CrossRef Full Text | Google Scholar

McDowell, W. H., Lugo, A. E., and James, A. (1995). Export of Nutrients and Major Ions from Caribbean Catchments. J. North Am. Benthological Soc. 14, 12–20. doi:10.2307/1467721

CrossRef Full Text | Google Scholar

McDowell, W. H., Sánchez, C. G., Asbury, C. E., and Ramos Pérez, C. R. (1990). Influence of Sea Salt Aerosols and Long Range Transport on Precipitation Chemistry at El Verde, Puerto Rico. Atmos. Environ. A. Gen. Top. 24 (11), 2813–2821. doi:10.1016/0960-1686(90)90168-m

CrossRef Full Text | Google Scholar

McDowell, W. H., Scatena, F. N., Waide, R. B., Brokaw, N., Camilo, G. R., Covich, A. P., et al. (2012). “Geographic and Ecological Setting of the Luquillo Mountains,” in A Caribbean Forest Tapestry: The Multidimensional Nature of Disturbance and Response. Editors T. A. C. N. Brokaw, A. E. Lugo, W. H. McDowell, F. N. Scatena, R. B. Waide, and M. W. Willig (New York: Oxford University Press), 72–163. doi:10.1093/acprof:osobl/9780195334692.003.0003

CrossRef Full Text | Google Scholar

Moore, O. W., Buss, H. L., and Dosseto, A. (2019). Incipient Chemical Weathering at Bedrock Fracture Interfaces in a Tropical Critical Zone System, Puerto Rico. Geochimica et Cosmochimica Acta 252, 61–87. doi:10.1016/j.gca.2019.1002.102810.1016/j.gca.2019.02.028

CrossRef Full Text | Google Scholar

Murphy, S. F., Brantley, S. L., Blum, A. E., White, A. F., and Dong, H. (1998). Chemical Weathering in a Tropical Watershed, Luquillo Mountains, Puerto Rico: II. Rate and Mechanism of Biotite Weathering. Geochimica et Cosmochimica Acta 62 (2), 227–243. doi:10.1016/s0016-7037(97)00336-0

CrossRef Full Text | Google Scholar

Murphy, S. F., Stallard, R. F., Scholl, M. A., González, G., and Torres-Sánchez, A. J. (2017). Reassessing Rainfall in the Luquillo Mountains, Puerto Rico: Local and Global Ecohydrological Implications. Plos One. doi:10.1371/journal.pone.0180987

CrossRef Full Text | Google Scholar

Navarre-Sitchler, A. K., Cole, D. R., Rother, G., Jin, L., Buss, H. L., and Brantley, S. L. (2013). Porosity and Surface Area Evolution during Weathering of Two Igneous Rocks. Geochimica et Cosmochimica Acta 109, 400–413. doi:10.1016/j.gca.2013.02.012

CrossRef Full Text | Google Scholar

Ogle, C. J. (1970). “Pollen Analysis of Selected Sphagnum Bog Sites in Puerto Rico,” in A Tropical Rain forest: A Study of Irradiation and Ecology at El Verde. Editors H. T. Odum, and R. F. Pigeon (Puerto Rico, B135–B145.

Google Scholar

Orlando, J., Comas, X., Hynek, S. A., Buss, H. L., and Brantley, S. L. (2016). Architecture of the Deep Critical Zone in the Río Icacos Watershed (Luquillo Critical Zone Observatory, Puerto Rico) Inferred from Drilling and Ground Penetrating Radar (GPR). Earth Surf. Process. Landforms 41, 1826–1840. doi:10.1002/esp.3948

CrossRef Full Text | Google Scholar

Orlando, J. J. (2014). The Anatomy of Weathering Profiles on Different Lithologies in the Tropical forest of Northeastern Puerto Rico: From Bedrock to Clouds. Park, PA (USA): M.S., Pennsylvania State University, University.

Google Scholar

Pett-Ridge, J. C. (2009). Contributions of Dust to Phosphorus Cycling in Tropical Forests of the Luquillo Mountains, Puerto Rico. Biogeochemistry 94 (1), 63–80. doi:10.1007/s10533-009-9308-x

CrossRef Full Text | Google Scholar

Pett-Ridge, J. C., Derry, L. A., and Barrows, J. K. (2009a). Ca/Sr and 87Sr/86Sr Ratios as Tracers of Ca and Sr Cycling in the Rio Icacos Watershed, Luquillo Mountains, Puerto Rico. Chem. Geology. 267, 32–45. doi:10.1016/j.chemgeo.2008.11.022

CrossRef Full Text | Google Scholar

Pett-Ridge, J. C., Derry, L. A., and Kurtz, A. C. (2009b). Sr Isotopes as a Tracer of Weathering Processes and Dust Inputs in a Tropical Granitoid Watershed, Luquillo Mountains, Puerto Rico. Geochimica et Cosmochimica Acta 73 (1), 25–43. doi:10.1016/j.gca.2008.09.032

CrossRef Full Text | Google Scholar

Porder, S., Johnson, A. H., Xing, H. X., Brocard, G., Goldsmith, S., and Pett-Ridge, J. (2015). Linking Geomorphology, Weathering and Cation Availability in the Luquillo Mountains of Puerto Rico. Geoderma 249-250, 100–110. doi:10.1016/j.geoderma.2015.03.002

CrossRef Full Text | Google Scholar

Reid, E. A., Reid, J. S., Meier, M. M., Dunlap, M. R., Cliff, S. S., Broumas, A., et al. (2003). Characterization of African Dust Transported to Puerto Rico by Individual Particle and Size Segregated Bulk Analysis. J. Geophys. Res. Lett. 108 (D19), 1–22. doi:10.1029/2002jd002935

CrossRef Full Text | Google Scholar

Rempe, D. M., and Dietrich, W. E. (2014). A Bottom-Up Control on Fresh-Bedrock Topography under Landscapes. Proc. Natl. Acad. Sci. U.S.A. 111 (18), 6576–6581. doi:10.1073/pnas.1404763111

PubMed Abstract | CrossRef Full Text | Google Scholar

Riebe, C. S., Hahm, W. J., and Brantley, S. L. (2016). Controls on Deep Critical Zone Architecture: a Historical Review and Four Testable Hypotheses. Earth Surf. Process. Landforms 42 (1), 128–156. doi:10.1002/esp.4052

CrossRef Full Text | Google Scholar

Riebe, C. S., Kirchner, J. W., and Finkel, R. C. (2004). Erosional and Climatic Effects on Long-Term Chemical Weathering Rates in Granitic Landscapes Spanning Diverse Climate Regimes. Earth Planet. Sci. Lett. 224, 547–562. doi:10.1016/j.epsl.2004.05.019

CrossRef Full Text | Google Scholar

Schellekens, J., Scatena, F. N., Bruijnzeel, L. A., van Dijk, A. I. J. M., Groen, M. M. A., and Van Hogezand, R. J. P. (2004). Stormflow Generation in a Small Rainforest Catchment in the Luquillo Experimental Forest, Puerto Rico. Hydrol. Process. 18, 505–530. doi:10.1002/hyp.1335

CrossRef Full Text | Google Scholar

Schulz, M. S., and White, A. F. (1999). Chemical Weathering in a Tropical Watershed, Luquillo Mountains, Puerto Rico: III. Quartz Dissolution Rates. Geochimica et Cosmochimica Acta 63 (3-4), 337–350. doi:10.1016/s0016-7037(99)00056-3

CrossRef Full Text | Google Scholar

Seiders, V. M. (1971). Cretaceous and Lower Tertiary Stratigraphy of the Gurabo and El Yunque Quadrangles, Puerto Rico. Virginia, US: United States Geological Survey Bulletin, 1–58.

Google Scholar

S. F. Murphy, and R. F. Stallard (Editors) (2012). Water Quality and Landscape Processes in Four Watersheds in Eastern Puerto Rico, USGS Professional Paper 1789-D (U.S. Geological Survey).

Google Scholar

Shanley, J. B., McDowell, W. H., and Stallard, R. F. (2011). Long-term Patterns and Short-Term Dynamics of Stream Solutes and Suspended Sediment in a Rapidly Weathering Tropical Watershed. Water Resour. Res. 47. doi:10.1029/2010WR009788

CrossRef Full Text | Google Scholar

Smith, A. L., Schellekens, J. H., and Díaz, A.-L. M. (1998). “Batholiths as Markers of Tectonic Change in the Northeastern Caribbean,” in Tectonics and Geochemistry of the Northeastern Caribbea. Editors E. G. Lidiak, and D. K. Larue (Boulder, Colorado (USA): Geological Society of America Special paper), 99–122. doi:10.1130/0-8137-2322-1.99

CrossRef Full Text | Google Scholar

Stallard, R. F. (2012). “Atmospheric Inputs to Watersheds of the Luquillo Mountains in Eastern Puerto Rico, Chapter D,” in Water Quality and Landscape Processes in Four Watersheds in Eastern Puerto Rico, USGS Professional Paper 1789-D. Editors S. F. Murphy, and R. F. Stallard (Reston VA (USA): U.S. Geological Survey), 89–112.

Google Scholar

St. Clair, J., Moon, S., Holbrook, W. S., Perron, J. T., Riebe, C. S., Martel, S. J., et al. (2015). Geophysical Imaging Reveals Topographic Stress Control of Bedrock Weathering. Science 350 (6260), 534–538. doi:10.1126/science.aab2210

PubMed Abstract | CrossRef Full Text | Google Scholar

ten Brink, U. (2005). Vertical Motions of the Puerto Rico Trench and Puerto Rico and Their Cause. J. Geophys. Res. 110, B06404, 06410. doi:10.1029/2004jb003459

CrossRef Full Text | Google Scholar

Turner, B. F., Stallard, R. F., and Brantley, S. L. (2003). Investigation of In Situ Weathering of Quartz Diorite Bedrock in the Rio Icacos basin, Luquillo Experimental Forest, Puerto Rico. Chem. Geology. 202 (3-4), 313–341. doi:10.1016/j.chemgeo.2003.05.001

CrossRef Full Text | Google Scholar

White, A. F., Blum, A. E., Schulz, M. S., Vivit, D. V., Stonestrom, D. A., Larsen, M., et al. (1998). Chemical Weathering in a Tropical Watershed, Luquillo Mountains, Puerto Rico: I. Long-Term versus Short-Term Weathering Fluxes. Geochimica et Cosmochimica Acta 62 (2), 209–226. doi:10.1016/s0016-7037(97)00335-9

CrossRef Full Text | Google Scholar

Yi-Balan, S. A., Amundson, R., and Buss, H. L. (2014). Decoupling of Sulfur and Nitrogen Cycling Due to Biotic Processes in a Tropical Rainforest. Geochimica et Cosmochimica Acta 142, 411–428. doi:10.1016/j.gca.2014.05.049

CrossRef Full Text | Google Scholar