Research ArticleGEOLOGY

Porosity production in weathered rock: Where volumetric strain dominates over chemical mass loss

See allHide authors and affiliations

Science Advances  18 Sep 2019:
Vol. 5, no. 9, eaao0834
DOI: 10.1126/sciadv.aao0834


Weathering in the critical zone causes volumetric strain and mass loss, thereby creating subsurface porosity that is vital to overlying ecosystems. We used geochemical and geophysical measurements to quantify the relative importance of volumetric strain and mass loss---the physical and chemical components of porosity---in weathering of granitic saprolite of the southern Sierra Nevada, California, USA. Porosity and strain decrease with depth and imply that saprolite more than doubles in volume during exhumation to the surface by erosion. Chemical depletion is relatively uniform, indicating that changes in porosity are dominated by processes that cause strain with little mass loss. Strain-induced porosity production at our site may arise from root wedging, biotite weathering, frost cracking, and the opening of fractures under ambient topographic stresses. Our analysis challenges the conventional view that volumetric strain can be assumed to be negligible as a porosity-producing mechanism in saprolite.


Weathering in mountain landscapes creates porosity in rock, saprolite, and soil near Earth’s surface through a combination of volumetric strain and chemical mass loss. Porosity, in turn, provides subsurface storage space and pathways for water (1, 2), enables access to mineral-bound nutrients for organisms (3), and thereby creates habitable substrates for overlying ecosystems (4). The opening of subsurface pore space also helps regulate mineral weathering by moderating fluid flow (5) and the influx of reactants such as dissolved O2 and CO2 (6). By influencing subsurface weathering, porosity production may also affect erosion rates at the surface (7) and thus may play a vital role in landscape evolution (8). Understanding how porosity is produced by volumetric strain and mass loss is therefore important across a broad range of problems in hydrology, biogeochemistry, ecology, and geomorphology (911).

Despite a long history of subsurface weathering studies (4, 7, 9, 10, 1216), the relative importance of physical, chemical, and biological processes has remained controversial. Saprolite is the zone of weathered rock that retains the relative positions of mineral grains (referred to here as “texture”) of the protolith (i.e., parent bedrock) and lies below the layer of bioturbation commonly referred to as “soil” (Fig. 1, A and B). The absence of bioturbation and the retention of protolith texture have commonly been cited as evidence that weathering is isovolumetric in saprolite (17, 18) to the point that a lack of strain has become part of the definition of saprolite in textbooks (19, 20). As a consequence, strain is not commonly measured in studies of subsurface weathering and is instead often assumed to be negligible. However, the original definition of saprolite (21) made no mention of isovolumetric weathering. Rather, it merely specified that the material is “untransported,” meaning that it is neither vertically mixed nor moved down slope, and these conditions are not sufficient to rule out volumetric strain. Dilation and collapse can be induced by physical and chemical processes without vertical mixing or downslope transport that would disrupt protolith texture. For example, in mountain landscapes, erosion at the surface exhumes bedrock from depth, exposing it to gradients in the ambient stress field that can cause existing pores and fractures to open (2225) without changes in the relative positions of mineral grains.

Fig. 1 Conceptual model of variations in porosity, strain, and mass loss in a weathering profile.

From top to bottom: The weathering profile (A) consists of the mobile soil layer (not resolved in our study), saprolite (which is weathered but retains the texture of underlying bedrock), fractured rock, and fresh bedrock at depth. In mountain landscapes, rock is exhumed through the weathering profile by erosion at the surface. The increase in porosity with proximity to the landscape surface (B) can be accounted for by mass loss (C), volumetric strain (D), or a combination of both (E).

Exhumation to the surface may lead to additional porosity production by a number of well-studied chemical, physical, and biological processes that could contribute to strain without disrupting the original protolith texture. For example, mineral expansion during weathering transformations can induce strain in the surrounding rock (9), thereby opening existing pore space and creating new microfractures without causing mixing that would corrupt the protolith texture (4, 9, 26, 27). In addition, frost cracking via segregation ice growth (28) can induce physical strain, or “damage,” both in soil and in the uppermost (near-surface) levels of saprolite, depending on ambient thermal gradients (29, 30). Likewise, biomechanical forcing from root wedging (31, 32) and hyphal growth (33) can expand preexisting pore space without causing vertical mixing. Thus, a variety of chemical, physical, and biological processes can induce volumetric strain in saprolite. This implies that the assumption of isovolumetric weathering in saprolite may be invalid in many landscapes.

However, the assumption is difficult to evaluate because volumetric strain has rarely been measured in studies of subsurface weathering. When it is measured along with chemical mass loss, these two components of weathering can be used together to evaluate their relative importance in creating subsurface porosity (13). For example, if chemical mass loss dominates over physical (volumetric) strain (34), then changes in porosity with depth should be strongly reflected by changes in mass loss (Fig. 1C). Conversely, if volumetric strain dominates over mass loss, then any changes in porosity with depth should be primarily reflected in changes in strain (Fig. 1D). Because minerals exhumed to the surface are typically exposed to increasing throughflow of reactant-rich fluids from Earth’s atmosphere (35), the production of porosity may be more complicated than either of the end-member cases shown in Fig. 1 (C and D), reflecting a combination of physical and chemical processes that cause both volumetric strain and chemical mass loss [e.g., Fig. 1E and (9)].

Although it has not previously been shown in the literature, the combined effects of volumetric strain (ε) and mass loss (τb) on subsurface porosity (ϕw) can be expressed in Eq. 1, which is derived from mass balance principles summarized in Supplementary Text.ϕw=1τb+1ε+1(1)

In Eq. 1, the term τb, referred to hereafter as “bulk tau,” is the abundance-weighted sum of the mass transfer coefficients (τi,j) for all chemical elements present in the bedrock. Bulk tau can be expressed in terms of immobile element concentrations as shown in Eq. 2 (see Supplementary Text for derivation).τb=j=1nCj,pτi,j=Ci,pCi,w1(2)

Here, C is the concentration, subscript j refers to a specific chemical element in either the protolith (subscript p) or the weathered material (subscript w), and n is the total number of chemical elements present in the protolith. The subscript i refers to an element that is conserved (i.e., immobile) during chemical losses of other more soluble elements. The rightmost term in Eq. 2 shows that the overall mass loss of a soil or saprolite can be calculated from the enrichment of an immobile element relative to its concentration in the protolith. When enrichment is negligible, bulk tau is roughly 0, reflecting negligible mass loss. Conversely, when enrichment is very high (i.e., Ci,wCi,p), bulk tau approaches −1, reflecting large mass losses approaching 100% of the material present in the protolith.

Immobile element enrichment ratios can also help quantify volumetric strain as shown in Eq. 3 (13).ε=VwVp1=ρb,pCi,pρb,wCi,w1(3)

Here, V is the volume and ρb is the bulk density. Positive strain indicates dilation, negative strain indicates contraction, and ε = 1 (i.e., 100% strain) implies a doubling in protolith volume during weathering (13).

Collectively, Eq. 1 to 3 and the three scenarios in Fig. 1 provide a framework for understanding variations observed in weathering profiles. Scenario 1 represents the conventional view that saprolite weathers without strain (Fig. 1C). Although strain has rarely been measured in studies of saprolite weathering, several examples of isovolumetric saprolite weathering have been reported in the literature. For example, measurements of bulk density and geochemistry from the Luquillo Critical Zone Observatory, Puerto Rico, suggest that strain is negligible across nearly the entire weathering profile (34). Similarly, volumetric strain in saprolite at Panola Mountain, Georgia, USA, was judged to be small compared to chemical mass losses (36).

Here, we demonstrate the critical importance of strain in porosity production within granitic saprolite using data from a new study of subsurface weathering at the Southern Sierra Critical Zone Observatory (SSCZO). Extensive existing measurements of bulk geochemistry (7, 37, 38) and geophysics (39) make the SSCZO ideal for our analyses. We focused on a hillslope in the Providence Creek drainage (Fig. 2, A to D), which is underlain by Dinkey Creek Granodiorite (40), was not glaciated during the Pleistocene (41), is canopied by a mixed conifer forest dominated by white fir, and receives an average of ~110 cm/year of precipitation mostly as snow (42). Forests at the ~2-km elevation of our site support year-round growth despite a long summer dry season (43). This implies that forest productivity is sustained by water stored in the subsurface (44), underscoring the importance of understanding porosity production in saprolite at this site. Regolith (the combination of soil, saprolite, and fractured bedrock that mantle unweathered bedrock at depth) is typically thicker than 10 m (39), and erosion rates are generally less than 50 mm ka−1 (thousand years) (7, 38) despite variability implied by knickpoints in stream channels and the surrounding topography (8). This indicates that residence times in the weathering profile are commonly greater than 200 ka (8), allowing for long exposure to weathering and thereby potentially explaining the high porosities (>40%) observed in near-surface saprolite at our site (39). Building on existing geophysical observations of subsurface weathering (Fig. 2E), we present new measurements of bulk geochemistry and show that volumetric strain dominates over mass loss in saprolite weathering along a forested slope at our site. We also present a predictive model that uses the bulk geochemical data to explain 92% of the variance in geophysical surveys of seismic velocities. In addition, we show how the data can be coupled with seismic refraction surveys conducted at the surface to predict both vertical and lateral variations in mass loss and volumetric strain over hillslope scales.

Fig. 2 Geophysical survey location and results.

Locations of SSCZO in California (A) and P301 catchment (thick blue outline) in the headwaters of Providence Creek (B) showing orientation of geophysical survey (red line). Southernmost point on survey located at 37.067657°N, 119.194315°W, and 2016-m elevation above mean sea level (projection: World Geodetic System of 1984). Thin blue lines are 10-m elevation contours. N, north; CA, California; ID, Idaho; NV, Nevada; UT, Utah, OR, Oregon; AZ, Arizona. (C) Bottom of a pit excavated in P301, near the geophysical survey, showing roots in soil profile. Measuring tape divisions are in decimeters. Bottom of the pit shows transition from dark, organic-rich, bioturbated soil to lighter, organic-poor, decomposed but nontransported saprolite. (D) Subsamples from 10 of the 82 outcrop samples used to calculate average bulk geochemistry of the protolith in (38). These billets show the typical scale, texture, and mineralogy of interlocking grains in the Dinkey Creek Granodiorite. Seismic velocity tomogram (E) across forested hillslope and meadow with locations of Geoprobe cores (CZG-1, CZG-6, and CZG-7) showing P-wave velocity contours (thin black lines) at 1, 2, 3, 4, and 5 km/s. (F) Porosity calculated from rock physics model using travel time tomography and fractional saturation measured from cores (see Materials and Methods) for saprolite, defined here as material with a P-wave velocity of <1.2 km/s (see text). Stippling shows depths with P-wave velocities >1.2 km/s. Seismic velocities are lowest and porosities are correspondingly highest under the forested ridge. (Photo credit: C. Riebe).


Variations in bulk density, porosity, mass loss, and strain

Our geochemical analyses focused on saprolite obtained from three push cores collected along a previously studied geophysical transect. We report measurements of bulk density, porosity, and bulk geochemistry in data file S1, which also includes values of strain and mass transfer coefficients that were calculated using Eqs. 2 and 3. Across 26 saprolite samples, bulk density increases with depth (r2 = 0.62, P < 0.0001; table S1), ranging from 0.96 g cm−3 at the surface to 1.71 g cm−3 at the bottom of the deepest core (Fig. 3A). This corresponds to a decrease in porosity from 0.65 to 0.35 (Fig. 3B), which is calculated according to Eq. 4 (Materials and Methods). Meanwhile, bulk tau, calculated from Eq. 2, ranges from −0.006 to −0.35 but does not vary systematically with depth (r2 = 0.02, P = 0.53; table S1), fluctuating around a mean ± SE of −0.22 ± 0.01 (Fig. 3C); this reflects the relatively narrow range in saprolite zirconium (Zr) concentrations [130 to 199 parts per million (ppm)] and the absence of systematic variability in Zr with depth. Thus, the variations in porosity with depth (Fig. 3B) are not well explained by variations in mass loss (Fig. 3C), i.e., porosity and bulk tau are not well correlated (Fig. 4A) (r2 = 0.07, P = 0.19; table S1). In contrast, volumetric strain increases by nearly a factor of 3, ranging from 31% at the bottom of the core to 112% at the surface (Fig. 3D), indicative of more than a doubling of protolith volume during exhumation. Therefore, strain is strongly correlated with porosity (Fig. 4B) (r2 = 0.62, P < 0.0001; table S1), indicating that changes in porosity are well explained by changes in volumetric strain with depth (cf. Fig. 3, D and B). The positive correlation between strain and porosity is consistent with what we would expect from Eq. 1 when mass loss does not vary systematically with depth. This result is nonetheless unexpected, given the conventional view of saprolite as bedrock that has weathered in place without volumetric strain.

Fig. 3 Bulk physical and chemical properties from Geoprobe cores.

Bulk density (A), porosity (B), mass loss (C), and volumetric strain (D) as a function of depth (±10 cm) from three Geoprobe cores at the SSCZO. Triangles, squares, and circles represent data from cores CZG-1, CZG-6, and CZG-7, respectively. Error bars estimated using Gaussian error propagation from ±1 SE of variability in physical properties and element concentrations of both the protolith and the weathered material. Bulk tau (C) is roughly uniform with depth, while porosity (B) and volumetric strain (D) both decrease with depth in the cores.

Fig. 4 Controls on porosity.

Cross-plots of porosity, bulk tau (A), and volumetric strain (B). Triangles, squares, and circles represent data from cores CZG-1, CZG-6, and CZG-7, respectively. Error bars estimated using Gaussian error propagation from ±1 SE of variability in physical properties and element concentrations of both the protolith and the weathered material. The strong trend in porosity with volumetric strain (and lack thereof with bulk tau) indicates that variations in porosity can be largely explained by variations in volumetric strain.

The elemental tau values (τi,j) of Na and Ca are more strongly coupled with depth than bulk tau (Fig. 5 and table S1), consistent with some progressive weathering of feldspar with increasing proximity to the surface. However, because Na and Ca make up a small fraction of the protolith (together averaging just 9.3% by weight as oxides), the trend in bulk tau—which is the chemical component of porosity in Eq. 1—is dominated by the lack of correlation between τi,j and depth for elements like Si, Al, and Fe (Fig. 5 and table S1), which, as oxides, collectively make up ~74%, on average, of the protolith by weight. This highlights the importance of using bulk tau, rather than τi,j, in analyses of porosity production by chemical mass loss, particularly when the element does not constitute a large fraction of the bedrock.

Fig. 5 Mass transfer coefficients of individual elements in saprolite.

Triangles, squares, and circles represent data from cores CZG-1, CZG-6, and CZG-7, respectively. All estimates use Zr as the immobile reference element. Labels in panels show weight percentages of elemental oxides in protolith. Error bars estimated using Gaussian error propagation from ±1 SE of variability in physical properties and element concentrations of both the protolith and the weathered material. Although CaO and Na2O show statistically significant correlations with depth (see table S1), together they make up only 9.3% of the underlying bedrock. Meanwhile, more abundant oxides such as SiO2, Al2O3, and Fe2O3 show no clear trend with depth, thus explaining the lack of a clear trend in bulk tau with depth (Fig. 3C).

Average bulk tau in the core samples indicates that ~22% of the protolith mass was lost below the deepest limits of our cores. This contributed substantially to porosity observed in the cores (see Eq. 1), but the lack of variation in bulk tau across depth in our cores shows that changes in porosity in saprolite are dominated by volumetric strain in the top ~11 m beneath the landscape surface. Thus, the patterns observed here are most consistent with scenario 3 for porosity production in the critical zone (Fig. 1E), with a mass loss evident at the base of saprolite and increasing strain throughout the saprolite. However, in this case, we do not know where the 22% mass loss occurs because our cores do not reach into fractured bedrock, much less unweathered protolith.

Immobility of Zr, Ti, and Nb

In our analysis of strain and mass loss (Figs. 3 to 5), we assume that Zr is immobile, recognizing that Zr-bearing mineral phases such as zircon and baddeleyite are highly insoluble. While this assumption is common in solid-phase mass balance studies (14, 34, 45, 46), even highly refractory Zr-bearing minerals can suffer weathering losses (47), leading to underestimates of both bulk and elemental tau values. In addition, tau can be either underestimated or overestimated in profiles where Zr-poor or Zr-rich mineral phases, derived, e.g., from exogenous dust (48) or as by-products of weathering (49), have been preferentially transferred via eluviation across depths within the saprolite.

To evaluate whether or not Zr weathering or eluviation are confounding factors in our analysis, we compared patterns in bulk tau from Zr with bulk tau calculated using titanium (Ti) and niobium (Nb), which have sometimes been used as immobile tracers of weathering (47, 50). Zr enrichments produce the largest bulk tau values (Fig. 6A), consistent with Zr being the most immobile of the three prospective tracer elements. Significant Zr weathering is therefore unlikely. The absence of confounding effects on Zr concentrations is further supported by the lack of trends in element enrichment for Zr, Ti, and Nb, which also show considerable overlap (Fig. 6B and table S1), without mismatches in trend that would be evident if eluviation of dust or weathering by-products preferentially influenced one element more than the others. In addition, the chemical enrichments of Zr, Ti, and Nb do not decrease with depth and thus do not follow the trend that we would see if there had been significant eluviation of Zr-, Ti-, or Nb-poor secondary minerals from shallow to deeper levels in the profile. On the other hand, dust inputs have been shown to have a strong effect on ecosystem nutrient budgets in the region (51, 52), and downward eluviation of exogenous Zr-, Ti-, and Nb-poor dust from the surface, if present, would tend to suppress weathering-related enrichments in these elements near the surface, thus creating a more uniform pattern of bulk tau versus depth, similar to what we observe in Figs. 3C and 6A. Nevertheless, we can rule out dust as a significant confounding factor in our analyses because the uniform patterns in Zr, Ti, and Nb (Fig. 3C, Fig. 6, and table S1) are not parsimonious with contamination by exogenous dust fluxes: To discredit the finding that strain dominates porosity production in saprolite, the Zr, Ti, and Nb concentrations in the dust and the amount of added dust would need to have coincided to eliminate downward decreases in weathering-related enrichment of Zr, Ti, and Nb across each of the measured cores. In addition, the downward flux of added dust would need to have outpaced the conversion of saprolite (and added dust) to soil, which would at least partly counteract eluviation if present. So many coincidences across multiple saprolite-forming processes seem highly unlikely. A more parsimonious explanation is that the bulk tau values from the Zr concentrations robustly reflect a pattern of uniform chemical mass losses with depth. This, in turn, supports the conclusion that the observed increase in total porosity toward the surface dominantly reflects an increase in strain—the physical component of porosity.

Fig. 6 Evaluation of other potentially immobile elements.

(A) Bulk tau was calculated in three ways, assuming alternatively that Ti, Nb, and Zr are immobile elements. Bulk tau from Nb and Ti are generally lower than bulk tau from Zr, suggesting that Zr is the most immobile element in saprolite at our study site. Triangles, squares, and circles represent data from cores CZG-1, CZG-6, and CZG-7, respectively. (B) Immobile element enrichment in saprolite relative to protolith normalized to Zr enrichment for Nb (red) and Ti (blue). The lack of any clear trends with depth suggests that our estimates of bulk tau are not strongly influenced either by translocation of clays within the profiles or by infiltration of exogenous dust (see text).

Physical weathering predicted from seismic velocities

The observation that bulk tau, the chemical component of porosity, is roughly uniform across the weathering profile implies that we can predict strain, the physical component of porosity, everywhere that total porosity is known. In other words, we can use the subsurface geophysical image of total porosity values (Fig. 2F) to generate a two-dimensional cross section of subsurface strain values (Fig. 7A). To accomplish this, we applied a uniform bulk tau, equal to the measured average of −0.22, at each depth and distance along the transect to predict values of subsurface strain from the porosity cross section using Eq. 7 (Materials and Methods). Strain predicted in this way varies from nearly 0% under the meadow at the base of the slope to ~100% (representing a doubling of volume relative to protolith) at some depths under the forested ridge (Fig. 7A). However, we restrict our analysis to material that can be realistically modeled using Hertz-Mindlin contact theory and only report strain at depths with P-wave velocities less than 1.2 km/s (Fig. 2), a threshold that marks the transition between saprolite and fractured granitic bedrock elsewhere in the U.S. West (53). To determine whether the strain predictions are realistic, we compared them to observations from the cores at CZG-1, CZG-6, and CZG-7 (Fig. 7B). The predicted strain profile (dashed line) closely follows the running average of observed strain values (solid line). This implies that our combined geophysical and geochemical approach yields realistic predictions of strain in saprolite (Fig. 7C).

Fig. 7 Spatial distribution of strain under a forested slope.

(A) Volumetric strain estimated from seismic velocities under the assumption that chemical mass loss is uniform (22% of the original protolith across the profile). We restrict our analysis to saprolite, with P-wave velocities <1.2 km/s (see text); stippling shows depths with P-wave velocities >1.2 km/s. Predicted strain is highest under the forested ridge and lowest under the meadow. (B) Depth-averaged volumetric strain from cores (dashed line calculated using locally weighted scatterplot smoothing (LOWESS) with a 0.5-m window) closely follows trend predicted from seismic velocities (solid line), matching well in a predicted versus observed plot (C). The line in (C) is a reference 1:1 line of perfect agreement.

As a separate check on our ability to partition geophysics-based porosity values into their physical and chemical components, we predicted seismic velocities for each core using the observed values of volumetric strain and mass loss from the geochemical measurements (Fig. 8). Considered individually, the contributions to porosity from strain and mass loss are unable to generate velocities low enough to match the observed velocity profile (cf. cyan, purple, and red lines in Fig. 8). However, when contributions to porosity from both strain and mass loss are considered together, we observe a much closer match between predicted and observed velocities (cf. black and red lines in Fig. 8), confirming the importance of both factors and quantifying their relative influence on inferred seismic velocity structure of the subsurface.

Fig. 8 Variations in predicted and observed seismic velocities at core locations.

Observed average velocities (black line) at core locations compared with velocities predicted from three rock physics models. Model predicted from measured variations in bulk tau (light blue) is offset from observations by roughly 0.5 km/s from observations across all depths. Model predicted from measured variations in volumetric strain (dark blue) shows less offset, but a combination of both bulk tau and volumetric strain in the rock physics model provides the closest match (red).


Several mechanisms could help explain the observed variations in strain and thus porosity as a function of depth in saprolite. For example, the biosphere extends many meters deep into saprolite (54) and may often be responsible for the inception of mineral weathering at depth (32). Roots of coniferous trees present at our site can extend as deep as 7 m beneath the landscape surface (55). Given that root density generally decreases with increasing depth (56), the decrease in strain with depth could partly reflect diminishing effects of root wedging. This might help explain the higher strain and deeper extent of physical weathering and saprolite under the forested ridge relative to the adjacent meadow (Fig. 7A). Measurements of the lateral and vertical distribution of roots and microorganisms would help test this hypothesis.

Another possible explanation for the observed strain patterns is frost cracking due to segregation ice growth when temperatures fall within the “frost-cracking window” (28). Models of temperate mountain sites like the SSCZO show that frost cracking is most intense in the top 2 m and is negligible at depths greater than 4 m (29). This suggests that the influence of frost cracking is too shallow to produce the patterns in Fig. 7. However, the thick saprolite (39) and slow erosion rate (7, 8, 38) of our site imply 105-year regolith residence times (8, 52) that likely integrate effects of past periglacial processes (39, 57), when damage due to frost cracking may have extended deeper (29, 30). Thus, frost cracking could help explain some of the observed variations in strain.

Another commonly recognized mechanism of volumetric strain in granite weathering is stress induced by expansion of sheet silicates such as biotite and smectite clays (9, 26, 27). For example, expanding clay within the saprolite could help explain the observed near-surface increase in strain if clay is more abundant at shallower depths within the cores. However, we found no evidence of expanding clays in x-ray diffraction patterns of saprolite sampled from our site (e.g., fig. S2). On the other hand, biotite accounts for up to 28% of the granitic protolith in the southern Sierra Nevada (58), and biotite oxidation has been shown to cause expansion that can disaggregate crystalline bedrock along adjacent mineral grain boundaries (9), creating porosity that promotes throughflow of reactive fluids and thereby enhances chemical weathering (9, 26, 27). Hence, it seems likely that biotite oxidation plays a role in the patterns of subsurface strain observed at our site.

Changes in lithostatic pressure during exhumation as well as interactions between topographic and regional stress fields could also explain some of the observed patterns in strain. Emplacement of the Sierra Nevada Batholith at >4-km depth (59) and subsequent exhumation due to rock uplift and surface erosion have caused substantial changes in lithostatic pressure through time. The reduction in pressure during exhumation likely caused fracturing of the protolith (24), with implications for physical and chemical weathering of saprolite (10). In addition, modern topographic and regional stresses can influence the distribution of subsurface fractures: The integrated modern stress field produces vertical and lateral variations in subsurface stresses that can mimic seismic velocity patterns (25) and thus may predict subsurface porosity distributions. The observed lateral variations in seismic velocities at our site (Fig. 2) and elsewhere in the region (39) broadly match predictions for sites with weak regional compression or extension (25), which is, in turn, broadly consistent with the weak north-south compressional and stronger east-west tensional stresses that have been predicted for the region (60). Regardless of far-field stresses, the inferred increase in strain toward the surface is consistent with the general prediction from topographic stress models, i.e., that fractures increase in frequency and aperture toward the surface (24, 25).

Together, our observations and analyses suggest that a combination of processes may be responsible for the observed variations in subsurface strain along our profile. The biotite expansion and pressure reduction mechanisms can account for our deepest measurements of strain, which are difficult to explain by the other, more surficial mechanisms of frost cracking and root wedging. Closer to the surface, all the proposed mechanisms may conspire to promote the upward increase in volumetric strain. Hence, many factors likely contribute to the observed and predicted patterns in Figs. 2, 3, 4, and 7. Although we are unable to quantify the relative importance of different strain-producing mechanisms without additional process-based research, our analysis shows that physical weathering contributes substantially to subsurface water holding capacity, which is vital to sustaining mountain ecosystems during the region’s summer dry seasons and through its extended droughts (44, 61).

While we can identify several plausible strain-producing mechanisms, the lack of any clear trends in bulk and elemental tau with depth (Figs. 3C and 5) is more enigmatic. In saprolite containing reactive mineral phases, geochemical modeling predicts vertical gradients in element concentrations due to progressive dissolution following the downward propagation of weathering fronts (5, 6, 35). Bulk and elemental tau values should only be uniform in profiles where the element of interest has been completely lost via mineral dissolution (i.e., with tau = −1, corresponding to 100% loss). Yet, across our cores, bulk tau shows no trend with depth and indicates that just 22% of the protolith has been lost because of dissolution (Fig. 3C). Moreover, none of the major rock-forming elements (Si, Al, K, Ca, Na, Fe, or Mg) are completely depleted from the profiles (Fig. 5), and the most abundant of them (Si, Al, and Fe) have tau values that are roughly uniform with depth (Fig. 5 and see table S1 for regression statistics). This suggests that chemical erosion has been shut down in saprolite, even though it retains many phases that are commonly assumed to be reactive in geochemical modeling studies. According to reactive transport modeling, the presence of supposedly reactive phases implies that dissolution has not reached a “local equilibrium” limit (11, 35) and that it is instead “kinetically limited.” However, kinetically limited dissolution should produce clear vertical gradients in bulk and elemental tau (6), which we do not observe. One explanation for this apparent contradiction is that the host minerals of many of the elements in the saprolite have been chemically depleted to the point that no further chemical erosion is possible. The roughly constant degree of chemical depletion in the saprolite indicates that the rate of mass loss from the profile can increase only if there is an increase in the rate of fresh, unreacted mineral supply to the weathering zone via the conversion of protolith to saprolite at depth. This condition is analogous to “supply limited” chemical erosion (11, 6264), which has been recognized at catchment scales across many granitic sites around the world (62, 64).

Our results show that volumetric strain dominates over mass loss in porosity production within saprolite along a mountain hillslope at the SSCZO. They therefore challenge the conventional view of saprolite as rock that has weathered in place without volumetric strain. Our analysis also challenges the prediction that the downward propagation of weathering fronts should produce vertical gradients in inferred chemical loss for elements that have not been completely depleted from the profile. We hypothesize that variations in volumetric strain at our site are driven, in part, by gradients in subsurface stress, which are present everywhere across Earth’s terrestrial surface because of interactions between topography and tectonics. Subsurface gradients in the intensity of other strain-producing processes such as biotite weathering, frost cracking, and root wedging are also likely important at our site and at other mountain granitic sites around the world. We therefore suggest that the dominance of physical over chemical weathering measured at our site may be widespread in other landscapes. More measurements are needed to test this hypothesis. Studies of chemical losses are common, but they do not always quantify volumetric strain, thus ignoring a crucial indicator of physical weathering. Our analyses outline a coupled geochemical and geophysical framework for quantifying both mass loss and volumetric strain---the chemical and physical components of porosity---over hillslope scales by combining total subsurface porosity from seismic refraction surveys with immobile element enrichment data from bulk geochemical analyses. Our work highlights the value of coupling direct observations from drilling with indirect geophysical measurements in investigations of subsurface weathering and porosity.


Experimental design

Our analysis builds on existing information about weathering and surface processes at the SSCZO, a site of more than a decade of intensive research on hydrology, forest ecology, and biogeochemistry across the rain-snow transition of the Sierra Nevada, California (3739, 42, 43, 52). To quantify spatial variations in chemical mass losses (τb) in saprolite, we analyzed material from three cores collected using a Geoprobe 6610DT direct-push, dual-speed auger (39). Density, porosity, and water content were already measured in previous work (39). This allowed us to leverage existing data on physical properties in our analysis of connections between bulk geochemistry, density, and porosity (Eqs. 1 to 3). Because the cores were obtained from points along a previously published seismic refraction profile (Fig. 2), we were also able to connect the one-dimensional view of chemical and physical properties from the cores with a two-dimensional perspective on subsurface weathering from geophysics using rock physics modeling.

Sample collection and preparation

The three cores—CZG-1, CZG-6 and CZG-7—reached depths of 10.2, 5.7, and 11.4 m, respectively, corresponding to the depth of refusal of the Geoprobe during coring (see Fig. 2 for locations). In each case, the depth of refusal was above or coincident with the base of saprolite (Fig. 2), which, in other granitic landscapes, corresponds to a P-wave velocity threshold of ~1.2 km/s (53). All our samples were taken from below the mobile soil layer, which is typically thinner than 1 m at our site (7).

To prevent moisture loss from saprolite before analysis, plastic core tubes were sealed with vinyl endcaps and parafilm. In the laboratory, 1-m core sections were subsampled in 10-cm increments. To avoid contamination from infill of the borehole sidewall that occurs between coring intervals, we used only the bottom 10 cm from each 1-m core section. In this way, we obtained a total of 29 subsamples for bulk geochemical analyses.

The core material resembles saprolite observed in nearby roadcuts: It differs from overlying soil in that it is devoid of visible organic matter and displays crystalline rock texture (i.e., interlocking grain boundaries in cross section), suggesting that protolith texture (Fig. 2D) and structure have not been disrupted by vertical or lateral mixing. It would not have been possible to collect our samples using push coring techniques if they were not at least somewhat friable. Thus, we are highly confident that samples from our cores were properly classified as saprolite.

Geophysical data collection

Seismic refraction data acquisition was described in detail in previous work (39). We reiterate highlights here because the data were vital to our new analyses. The survey consisted of two 24-channel Geometrics Geode systems with 40-Hz vertical component geophones spaced at 5-m intervals. Seismic source energy was generated every ~15 m along the profile by a 12-pound sledgehammer striking a stainless steel plate. In some instances, source energy was generated from a 12-gauge shotgun blank fired into a 1- to 2-m-deep auger hole.

Quantifying physical properties

To quantify bulk density, porosity, and volumetric water content, the mass of each subsample was measured before and after drying it in an oven at 105°C for 24 hours. Bulk density was calculated using the dry mass and the known volume of the coring cylinder. Volume estimates were corrected for compaction on the basis of the amount of core recovered relative to the distance cored by the Geoprobe (39), a standard approach in push-core acquisition. Assuming an average mineral grain density (ρg) of 2.65 g cm−3, porosity was estimated from bulk density (ρb) using Eq. 4.ϕ=1ρbρg(4)

The mineral grain density assumed here should be reasonable for saprolite derived from a granodiorite and containing very little organic matter, as is the case in our samples. Volumetric water content was estimated from the difference in mass between the original and dried sample weights.

Geochemical analyses

Equations 2 and 3 show that ratios of immobile element concentrations in protolith and weathered material can be used to quantify volumetric strain (ε) and bulk mass loss (τb). We measured concentrations of major, minor, and trace elements in representative aliquots of our subsamples using x-ray fluorescence (XRF). The aliquots were first powdered to <50 μm in a tungsten carbide grinding pot using a SPEX Shatterbox. To drive off water and any organic matter, which constituted a tiny fraction of our saprolite samples, we incinerated powders at 550°C for 12 hours. Each powdered aliquot (1 g) was carefully massed and mixed with lithium tetraborate at a ratio of 1:9. The sample and flux mixture was fused into a glass disc for XRF analysis of major and minor element-oxide concentrations using standard techniques. We measured trace element concentrations using XRF analysis of pressed powders infused with SPEX UltraBind. All XRF analyses were performed at the University of Wyoming using a 4-kW wavelength-dispersive XRF spectrometer. The geochemical data are recorded in data file S1.

We quantified average protolith bulk geochemistry using data from 82 representative samples of fresh outcrops of Dinkey Creek Granodiorite, the bedrock underlying the study site. These samples, which were collected from surface outcrops and analyzed in previous work (38), were vital to estimating the chemical mass loss that has occurred in the saprolite sampled from the cores (Eq. 2). Their average concentrations of Zr, Ti, and Nb in protolith were 129.2 ± 3.5 ppm, 5.94 ± 0.15 mg/g, and 8.8 ± 0.2 ppm (mean ± SE), respectively (data file S1). The uncertainties in mean protolith element concentrations were propagated through the calculations of strain and mass transfer coefficient using Gaussian error propagation.

Geophysical analyses

We created a seismic velocity tomogram for the geophysical survey data using conventional techniques. The timing of first arrivals was picked manually from waveform data with high signal-to-noise ratios; noisy wave traces, although rare, were excluded from the analysis. We used a least-squares MATLAB-based tomographic inversion code (25) to convert the travel time data from the picking analysis into a subsurface map of seismic velocities (Fig. 2E) that closely resembles the model for the same seismic refraction data published in (39)

The subsurface velocity map is prone to errors in picking the timing of first arrivals (estimated to be ±1 ms) and variations arising from assumptions in the starting model. However, the resulting uncertainties in P-wave structure are likely limited to 5 to 10% for this dataset based on sensitivity analyses presented in (39). Anisotropy is another potential source of error in the inferred P-wave velocity structure. Although rock and therefore saprolite textures are likely isotropic at our site (40), it is difficult, on the basis of available data, to quantify anisotropy that might result from subsurface fracturing. A previous study of granitic saprolite elsewhere in the western United States showed that anisotropy in P-wave velocity structure decreased with depth from 20 to 25% at 2 m to 10 to 15% at 4 m (65).

Rock physics model of porosity

To estimate porosity from the P-wave velocity model, we used rock physics relationships following previous work (39) but, in this case, explicitly incorporated average fractional saturation measured from cores using Brie’s fluid-mixing model (66). First, we calculated the elastic moduli of individual minerals with a Voigt-Reuss-Hill average (67) and then calculated the bulk and shear moduli of the saprolite using Hertz-Mindlin contact theory with a critical porosity of 0.38 and a Hashin-Strickman lower boundary (39, 68). Next, we used the bulk modulus of the fluid (66) in Gassman’s (69) equation to find the bulk and shear moduli of partially saturated saprolite. Last, we used the resulting estimates of physical properties in Eq. 5 to calculate P-wave velocity over a range of saturation and porosity values.U=K+43Gρavg(5)

Here, U is the P-wave velocity, K is the bulk modulus, G is the shear modulus, and ρavg is the average of the mineral and fluid densities. Velocity is also influenced by overburden pressure and thus depth beneath the surface. However, for a given depth, mineral assemblage, and range of fractional saturation values, P-wave velocity should have a specific relationship with porosity (e.g., fig. S1). Thus, we can estimate porosity from the model at each depth using velocities from travel time tomography and fractional saturation measured from Geoprobe cores. To estimate porosity at depths below the lower limits of the cores where no fractional saturation measurements are available, we used linear extrapolation from values measured in the top 11 m. However, to restrict our analysis to material that can be realistically modeled using Hertz-Mindlin contact theory (which applies to disaggregated granular material such as granitic saprolite), we only report porosity at depths with U < 1.2 km/s (Fig. 2), which should limit our analysis to saprolite, based on the measured transition between saprolite and fractured granitic bedrock in the Laramie Range, Wyoming, USA (53). The velocity threshold, U = 1.2 km/s, also corresponds to the deepest limits of CZG-7 (Fig. 2F), implying that our geochemical measurements span the entire thickness of saprolite at that core’s location.

The agreement between predicted and observed patterns of subsurface strain (Fig. 7) suggests that the porosity model yields realistic results despite several assumptions that introduce potential for uncertainty. For example, modal mineralogy, which influences both the bulk and shear moduli in Hertz-Mindlin contact theory, was assumed here to be 50% feldspar, 25% quartz, and 25% clay [after (39)]. The assumed clay content is broadly consistent with the observed bulk tau values (i.e., −0.22 ± 0.01) and is in the middle of the range of values explored in a previously reported sensitivity analysis of the dataset used here. The calculated sensitivity of porosity to assumptions about clay content decreases with depth from ±7.5% at 5 m to ±2% at 15 m for clay contents ranging from 0 to 65% (39). Modeled porosity is also sensitive to degree of saturation, but this likely contributes minimally to errors in predicted porosity given that our model incorporates direct observations of saturation from the cores.

Linking geophysics and geochemistry

Equation 1 shows that porosity can be quantified by combining estimates of chemical mass loss (quantified with Eq. 2) and volumetric strain (quantified with Eq. 3). By rearranging Eq. 1, we can also express τb and ε explicitly.τb=εϕwεϕw(6)ε=τb+1ϕw11(7)

Equations 6 and 7 show that either bulk mass loss or strain can be inferred by combining subsurface porosity distributions (inferred from geophysics and rock physics modeling) with estimates of strain or bulk mass loss (inferred from bulk geochemistry).


Supplementary material for this article is available at

Supplementary Text

Fig. S1. Rock physics model of porosity.

Fig. S2. Representative x-ray diffraction pattern from separated clay fraction.

Table S1. Regression statistics.

Data file S1. Physical properties, bulk geochemistry, mass transfer coefficients, and volumetric strain.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank N. Swoboda-Colberg, M. Schell, and C. Fobes for help with sample preparation and analysis and H. Buss and J. Kaszuba for insightful discussions. Funding: J.L.H., C.S.R., W.S.H., and B.A.F. were supported by NSF EPS-1208909. C.S.R. and P.C.H. were supported by NSF EAR-1331939. J.L.H. acknowledges support from Dickinson College Research and Development. Author contributions: J.L.H. and C.S.R. contributed equally to manuscript preparation. J.L.H., C.S.R., and W.S.H. contributed equally to study design. J.L.H. and C.S.R. conducted or oversaw all geochemical analyses. W.S.H. and C.S.R. designed and oversaw the geophysical data acquisition. B.A.F. conducted all rock physics modeling of porosity. P.C.H. conducted or oversaw Geoprobe sample collection and all physical analyses of cores. All authors discussed the results and commented on the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article