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).

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).

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.

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).

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.

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.

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).

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).

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.

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.

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.

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).

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).

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.

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.

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.

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).

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.

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”.

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.

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).

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.

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).

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).

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.