Geophysical imaging reveals topographic stress control of bedrock weathering

See allHide authors and affiliations

Science  30 Oct 2015:
Vol. 350, Issue 6260, pp. 534-538
DOI: 10.1126/science.aab2210

Bedrock weathering runs to the hills

Fractures in bedrock drive the breakdown of rock into soil. Soil makes observations of bedrock processes challenging. St. Clair et al. combined a three-dimensional stress model with geophysical measurements to show that bedrock erosion rates mirror changes in topography (see the Perspective by Anderson). Seismic reflection and electromagnetic profiles allowed mapping of the bedrock fracture density. The profiles mirror changes in surface elevation and thus provide a way to study the critical zone between rock and soil.

Science, this issue p. 534; see also p. 506


Bedrock fracture systems facilitate weathering, allowing fresh mineral surfaces to interact with corrosive waters and biota from Earth’s surface, while simultaneously promoting drainage of chemically equilibrated fluids. We show that topographic perturbations to regional stress fields explain bedrock fracture distributions, as revealed by seismic velocity and electrical resistivity surveys from three landscapes. The base of the fracture-rich zone mirrors surface topography where the ratio of horizontal compressive tectonic stresses to near-surface gravitational stresses is relatively large, and it parallels the surface topography where the ratio is relatively small. Three-dimensional stress calculations predict these results, suggesting that tectonic stresses interact with topography to influence bedrock disaggregation, groundwater flow, chemical weathering, and the depth of the “critical zone” in which many biogeochemical processes occur.

Weathered bedrock and soil are part of the life-sustaining layer at Earth’s surface commonly referred to as the “critical zone.” Its thickness depends on the competition between erosion, which removes weathered material from the land surface, and weathering, which breaks down rock mechanically and chemically and thus deepens the interface between weathered and fresh bedrock (13). Erosion at the surface can be studied both directly by observation (4) and indirectly using isotopic tracers (5). In contrast, weathering at depth is generally obscured by overlying rock and soil, making it more difficult to study. In this Report, we combine results from landscape-scale geophysical imaging and stress-field modeling to evaluate hypotheses about the relationship of subsurface weathering to surface topography.

One hypothesis is that weathering at depth is regulated by the hydraulic conductivity and porosity of bedrock and the incision rate of river channels into bedrock, which together determine how rapidly groundwater flowing toward the river channels can drain bedrock pores of chemically equilibrated water (1). Another hypothesis, based on reactive transport modeling, predicts that the propagation of the weathered zone depends on the balance between mineral reaction kinetics and groundwater residence times (2, 6). These mechanisms are not mutually exclusive, and both emphasize the importance of fluid flow through bedrock. This implies that the development of fracture systems, which provide preferred fluid-flow paths in bedrock (7, 8), should help regulate the downward propagation of weathering.

Bedrock fractures may be inherited from past tectonic or thermal events, but the common observation that fracture abundance declines with depth (9) suggests that near-surface processes are capable of generating new fractures or reactivating existing ones. Examples of potential near-surface fracturing mechanisms include wedging due to growth of ice crystals (10), salt crystals (11), or roots (12); volume expansion due to weathering reactions of certain minerals (13); and stress perturbations associated with surface topography (14). Some of these mechanisms require specific conditions that only occur in certain geographic regions, rock types, or parts of the subsurface. For example, ice weathering requires repeated freeze-thaw cycles and is therefore limited to depths of a few meters, and volume expansion due to mineral weathering is most effective in rocks containing abundant biotite or hornblende. In contrast, surface topography always perturbs the bedrock stress field. Here we focus on this potentially widespread control on bedrock fractures.

Bedrock stress fields can be evaluated by considering the effects of local topography on gravitational stresses and regional tectonic stresses. Theoretical calculations have shown that the presence of topographic features such as ridges and valleys can cause vertical and lateral variations in the shallow subsurface stress field (15) and that these stress perturbations can be large enough to alter bedrock fracture patterns (14, 16). Some observational evidence suggests that topographic stresses influence near-surface fracture distributions (9, 17, 18) and groundwater flow (1921), but it is unknown whether topographic stresses systematically affect the spatial distribution of bedrock weathering. In this study, we used geophysical surveys of seismic velocity and electrical resistivity to image the thickness of the weathered zone in three landscapes with similar topography but different tectonic stress conditions, and we tested whether the geometry of the weathered zone in each landscape matches the modeled topographic stress field.

A calculation of the stresses beneath an idealized topographic profile in the presence of regional tectonic stress illustrates how topographic stress might influence the geometry of the weathered zone (Fig. 1). We calculated the total stress field as the sum of the ambient stress due to gravity and tectonics and the stress perturbation due to topography (22). From the stress field, we calculated two scalar quantities that act as proxies for two mechanisms that could influence the abundance of fractures. As a proxy for shear fracturing or shear sliding on existing fractures, we calculated the failure potential (Φ) (23), defined as (σmc – σlc)/(σmc + σlc), where σ indicates a stress, and the subscripts denote the most compressive (mc) and least compressive (lc) principal stresses, with compression being positive. As a proxy for opening-mode displacement on fractures, we used the magnitude of σlc. Because we do not have prior knowledge of the orientations or abundances of existing fractures in a given landscape, we used both quantities to represent the propensity for generating open fractures. A larger Φ indicates that new shear fractures are more likely to form and that existing fractures oriented obliquely to the σmc and σlc directions are more likely to dilate as a result of sliding on the rough fracture surfaces (24). A smaller σlc indicates that fractures oriented nearly perpendicular to the σlc direction will be more likely to open (20, 25). Given that open fractures permit water to flow more quickly through bedrock and enhance the rate of chemical weathering, we expect that zones of larger Φ or smaller σlc correspond to zones of more weathered bedrock.

Fig. 1 Hypothesized influence of topography and tectonic stress on weathering zone geometry.

(A to C) Cross-sectional plots of failure potential (Φ) computed with a numerical stress model (22). Arrows indicate the magnitude of ambient horizontal compression. (A) shows a horizontal surface under weak compression, (B) shows ridges and valleys under weak compression, and (C) shows ridges and valleys under strong compression. The magnitude of the horizontal compression at the mean elevation of the land surface is 0.25ρga in (A) and (B) and 2ρga in (C), where ρ is rock density, g is gravitational acceleration, and a is the horizontal distance from the valley bottom to the ridge top in (B) and (C). This corresponds to σ* values of 0.25 in (B) and 2 in (C) (Eq. 1). Topography is horizontally periodic. (D to F) The same scenarios as in (A) to (C), but showing the magnitude of the least compressive principal stress, σlc. (G to I) Conceptual diagrams of bedrock fracture abundance in the three scenarios.

In a scenario with a horizontal land surface and a stress field determined by both gravity and a low, uniform ambient horizontal compression, Φ declines and σlc increases with depth beneath the land surface (Fig. 1, A and D). If topography is added and the ambient horizontal compression remains low, contours of Φ and σlc generally parallel the surface (Fig. 1, B and E). This surface-parallel pattern occurs because σlc is determined primarily by the overburden in this scenario, and therefore it increases with depth similarly beneath ridges and valleys (fig. S1A). However, if the ambient horizontal stress becomes strongly compressive—so that σmc at shallow depths parallels the topographic surface, and σlc becomes small relative to σmc (fig. S1, E and F)—then contours of Φ and σlc resemble mirror images of the land surface, plunging below ridges and rising beneath valleys (Fig. 1, C and F). This mirror-image pattern occurs because σlc becomes much smaller beneath ridges than beneath valleys, whereas the fractional decrease in σmc is relatively minor. In all three scenarios, spatial variations in σlc are the dominant cause of the spatial variations in Φ (Fig. 1, A to F, and fig. S1).

Where contours of Φ and σlc are generally parallel to the land surface, we expect to find a zone of abundant open fractures and weathered bedrock with a roughly uniform thickness, underlain by less weathered bedrock with fewer open fractures (Fig. 1, G and H). Where contours of Φ and σlc mirror the land surface, we expect to find thick weathered zones beneath ridges and thin zones beneath valleys (Fig. 1I). This conceptual model assumes that the geometry of the weathered zone reflects the present-day stress field. In an eroding landscape, a parcel of rock experiences a time-varying stress field as the eroding surface draws nearer. However, given that Φ increases and σlc decreases toward the surface (provided that the ambient horizontal tectonic stresses are compressive) (Fig. 1, A to F), the bottom boundary of the weathered zone is likely to be demarcated by the contour of the lowest Φ or highest σlc at which fractures open enough to accelerate chemical weathering. After crossing this boundary, rocks will only experience higher Φ and lower σlc and should become even more weathered.

A dimensionless ratio of tectonic stresses to gravitational stresses that accounts for topographic perturbations captures the difference between the surface-parallel scenario (Fig. 1, B, E, and H) and the surface-mirroring scenario (Fig. 1, C, F, and I). This ratio is defined asEmbedded Image (1)where σt is the magnitude of the horizontal tectonic stress perpendicular to the ridges and valleys, which we refer to as the tectonic compression; b/a is the characteristic slope of a topographic feature with height b and horizontal length a; ρ is rock density; and g is gravitational acceleration. Our model predicts that a landscape with σ* < ~1, corresponding to low tectonic compression or widely spaced ridges and valleys, should have a zone of large Φ, small σlc, and weathered rock that parallels the surface topography (Fig. 1, B, E, and H). Conversely, a landscape with σ* > ~1, corresponding to high tectonic compression or closely spaced ridges and valleys, should have a zone of large Φ, small σlc, and weathered rock that thickens under ridges and thins under valleys (Fig. 1, C, F, and I). The surface-mirroring pattern is most pronounced for landscapes with b/a > 0.05.

To test our conceptual model, we investigated three field sites in the United States with similar topography but different ambient tectonic stress regimes. Gordon Gulch (Fig. 2B), in the Boulder Creek Critical Zone Observatory, Colorado (3), is underlain by Precambrian gneiss and has weak tectonic compression (σt = 1 MPa; σ* = 0.2). Calhoun Critical Zone Observatory, South Carolina (Fig. 2C) (26), is underlain by Neoproterozoic or Cambrian granite gneiss and has stronger tectonic compression (σt = 4.7 MPa; σ* = 1.2). Pond Branch, Maryland (Fig. 2D) (27), is underlain by Precambrian schist and has the strongest tectonic compression of the three sites (σt = 8.4 MPa; σ* = 3.2). At each site, we selected a valley with a cross-sectional relief of 20 to 40 m and used a three-dimensional stress model to calculate the topographic perturbation to the ambient stress field in the vicinity of the valley (22). Our model calculations include ~1 to 15 km2 of the topography surrounding our field sites (Fig. 2). Ambient stress magnitudes, directions, and depth gradients were constrained from compilations of regional measurements (fig. S2) (22).

Fig. 2 Maps of study sites.

(A) Locations of the three U.S. sites. Colors indicate elevation, with greens being lowest and browns being highest. (B) Gordon Gulch, Colorado. (C) Calhoun, South Carolina. (D) Pond Branch, Maryland. Geophysical survey lines are numbered. White lines mark the transects shown in Fig. 3. Elevation data are from the U.S. National Elevation Dataset and the National Center for Airborne Laser Mapping.

Transects of Φ and σlc across the valleys (Fig. 3, A to F) show patterns consistent with the idealized calculations shown in Fig. 1, A to F. A zone of high Φ and low σlc, with roughly uniform depth, parallels the land surface at Gordon Gulch (σ* = 0.2) (Fig. 3, A and D), whereas the zones of high Φ and low σlc at Calhoun (σ* = 1.2) and Pond Branch (σ* = 3.2) are deeper beneath ridges and shallower beneath valleys (Fig. 3, B, C, E, and F). Sensitivity tests accounting for uncertainties in the ambient stress estimates and variations in topography across the sites confirm that these patterns are robust (fig. S3 to S5) (22).

Fig. 3 Comparisons of topographic stress and seismic velocity.

Modeled failure potential (A to C), modeled magnitude of the least compressive stress (D to F), and measured P-wave velocity (G to I) for Gordon Gulch [(A), (D), and (G)], Calhoun [(B), (E), and (H)], and Pond Branch [(C), (F), and (I)]. Transect locations are marked in white in Fig. 2. The topographic surface was smoothed for stress calculations (17). Vertical exaggerations are 3.6× at Gordon Gulch, 3.7× at Calhoun, and 2.3× at Pond Branch.

We performed seismic refraction and electrical resistivity surveys along the same transects to investigate whether topographic stresses influence the geometry of the weathered zone (22). The seismic velocity structure reflects variations in rock damage such as chemical weathering and fracturing (28), both of which contribute to an increase in porosity and a decrease in seismic velocity (29). Velocity in fractured rocks can also increase with higher water saturation, which increases the effective bulk modulus of the composite medium (28, 30), or higher confining stress, which closes microfractures formed during the rock’s thermal and tectonic history (31). We observed steep velocity gradients with depth, with P-wave velocities (vP) increasing downward from ~0.3 km/s near the surface to ~4.0 to 4.5 km/s at 10 to 40 m depth (Fig. 3, G to I), consistent with a transition from soil regolith to unweathered bedrock (28). The rapid increase in velocity with depth is likely also due in part to closure of fractures with increasing confining stress (31). At depths greater than the 4-km/s velocity contour, fractures are likely to be effectively closed, and the velocity gradients are greatly reduced. Borehole observations, including velocity measurements at fine depth increments (fig. S7B) and fractures mapped in optical image logs (fig. S10) (22), support the interpretation that faster velocities correspond to less abundant open fractures.

The variations in seismic velocity are highly similar in shape to the variations in modeled Φ and σlc, particularly in the zones of low velocity (vP < 4 km/s; Fig. 3, G to I), which generally correspond to zones of high Φ (Φ > 0.75; Fig. 3, A to C) and low σlclc < 0.5 MPa; Fig. 3, D to F). The 4-km/s velocity contour rises to within 10 m of the surface beneath valley bottoms and plunges to depths of 40 m beneath ridgetops at the sites with strong tectonic compression and σ* > 1 (Fig. 3, H and I). This is very similar to the mirror-image pattern of Φ and σlc contours that is predicted by the stress model (Fig. 3, B, C, E, and F). In contrast to these sites, both the 4-km/s velocity contour and the zones of high Φ and low σlc roughly parallel the land surface at Gordon Gulch (Fig. 3, A, D, and G), which has weak tectonic compression and σ* < 1. The velocity fields contain some fine-scale structures that do not appear in the modeled stress fields, but the stress models do predict some secondary features of the velocity structure. For example, the zones of high Φ and low seismic velocity at Pond Branch are both thicker beneath the right-hand side of the valley bottom (Fig. 3, C and F), and the zone of low σlc shows this asymmetry for slightly smaller values of σ* (fig. S4D). The bases of these three zones slope down to the right beneath both ridges at Calhoun (Fig. 3, B, E, and H).

At all three sites, Φ and σlc correlate well with seismic velocity (fig. S8). The trends of σlc versus vP for the three sites are very similar (fig. S8, D to F), whereas the trends of Φ versus vP vary (fig. S8, A to C), particularly at Gordon Gulch. This may indicate that opening of fractures oriented perpendicular to the σlc direction is a stronger control on seismic velocity than shear fracturing or sliding at these sites, or it may reflect differences in rock properties among sites that influence the relationship between stress and seismic velocity. In general, the strong correspondence between Φ, σlc, and seismic velocity at these sites suggests a robust influence of topographic stress on the physical state of weathered bedrock in the subsurface.

We measured electrical resistivity to determine whether the zones of high Φ, low σlc, and low seismic velocity also have high water content and more weathered rock. Rocks become less electrically resistive with increasing water saturation, connectedness of pores, and clay content (32). The resistivity images show structures similar to those in the seismic velocity images and reveal additional details about the hydrology of the weathered zone. For example, resistivity at Calhoun increases with depth, and we identified three regions that roughly correlate with the seismic velocity structure (Fig. 4). Where seismic velocities exceed 4 km/s, resistivity is generally high (>4000 ohm·m), consistent with low water content and poorly connected pores in intact bedrock. The upper 10 m of the weathered zone has low resistivity (100 to 700 ohm·m) and velocities < 1.5 km/s, consistent with a highly porous, water-saturated regolith. Between these zones is a zone of intermediate resistivity (700 to 4000 ohm·m) and velocity (1.5 to 4.0 km/s), consistent with fractured, slightly weathered bedrock with intermediate water content. We additionally observed a resistive body roughly 10 m beneath the ridgetop that has higher seismic velocity than the surrounding rock and is in close proximity to a region of locally reduced Φ and elevated σlc (Fig. 3, B and E).

Fig. 4 Electrical resistivity and seismic velocity at the Calhoun site.

The transect is the same as in Fig. 3, B, E, and H. Resistivity (A) and P-wave velocity (B) increase with depth, probably reflecting the closing of fractures. Both images show the surface-mirroring pattern predicted by the stress model (Fig. 3, B and E). Vertical exaggeration is 3×.

The geophysical observations support the hypothesis that topographic stress influences bedrock weathering by modifying the abundance of open fractures. This interpretation is consistent with the borehole sonic velocity log at Calhoun (fig. S7B) and the borehole image logs at Gordon Gulch (fig. S10) (22), both of which show lower velocities where fractures are more abundant. However, it is important to consider how factors other than fracture abundance might influence seismic velocity and resistivity. Higher water saturation can elevate seismic velocity, but our resistivity maps show the opposite: a near-surface layer with high water content but lower velocity (Fig. 4). Stress-induced pore pressure variations are also an unlikely explanation for the observed seismic velocity variations, because the velocity field resembles the least compressive stress rather than the mean stress. A more likely scenario is that areas with lower velocities and lower resistivity have experienced a higher degree of chemical weathering. This possibility does not conflict with the interpretation that fractures influence velocity and resistivity variations; more abundant, open, and connected fractures should promote faster groundwater flow and chemical weathering, which would enhance the connection between topographic stress and bedrock damage.

Our results suggest that topographic perturbations to gravitational and tectonic stress fields control the distribution of open fractures, permeability, and water storage potential in a weathered subsurface zone with a thickness comparable to the topographic relief. We propose that the base of the weathered zone is the depth above which stresses allow abundant fractures to grow or open, providing pathways for water and biota to invade the bedrock and initiate weathering. The stress field thus sets the physical stage for the interacting geochemical (2), hydrological (1), biological (33), and climatic (3) processes that drive weathering and nutrient cycling in Earth’s critical zone, which in turn help regulate the global carbon cycle (6). Lateral variations in stress-mediated weathering may also lead to variations in bedrock erodibility and the production rate of mobile soil. Topographic stress may ultimately govern the evolution of landforms through its influence on subsurface weathering and, thus, on erosion at the surface.

Supplementary Materials

Materials and Methods

Figs. S1 to S10

Tables S1 to S3

References (3455)


  1. Materials and methods are available as supplementary materials on Science Online.
  2. Acknowledgments: This study was supported by the NSF Experimental Program to Stimulate Competitive Research (award EPS-1208909), the U.S. Army Research Office (grant W911NF-14-1-0037), NSF Earth Sciences (awards EAR-0724960 to the University of Colorado–Boulder and EAR-1331846 to Duke University), and the NSF Critical Zone Observatory Network. Poly3D software was provided by Schlumberger. 3DMove software was provided by Midland Valley Software. Geophysical data presented in this paper are available through the Wyoming Center for Environmental Hydrology and Geophysics data repository at The authors thank the editor and three anonymous reviewers for their suggestions.
View Abstract

Stay Connected to Science

Navigate This Article