Beyond boundaries: Depositional environment controls on erodibility, process, and form in rivers incising sedimentary bedrock

Bedrock rivers set the relief structure of mountain landscapes, the rate of mass transport out of source regions, and the lower boundary condition for hillslopes. Rates of signal propagation along bedrock channels modulate the response time of drainage basins to environmental perturbations (Whipple et al., 2022). By controlling river and watershed shape as well as sediment transfer across Earth's surface, fluvial bedrock erosion processes influence both modern habitat availability for aquatic organisms (Crowder and Diplas, 2006; Wohl, 2015) and the evolution of species over geologic time (Stokes et al., 2023).

Bedrock erodibility, the conceptual inverse of channel boundary erosion resistance, affects bedrock channel form and process. Erodibility is a somewhat heuristic concept but can be thought of as “some combination of the strength of unfractured rock at the surface and… the degree of fracturing of the rock mass” (Whipple et al., 2022, p. 874), as well as the rock's susceptibility to chemical erosion. The intact rock strength, which governs the erodibility relevant for erosion processes like abrasion that occur at scales smaller than the fracture spacing, depends on mineralogy (Cowie and Walton, 2018), grain size (Fredrich et al., 1990), geologic history such as diagenesis or metamorphism (Cook et al., 2015), and weathering once the rock is exposed at Earth's surface (Hancock et al., 2011; Shobe et al., 2017). Controls on intact strength also affect the degree of fracturing (e.g., Marshall and Roering, 2014; McGinnis et al., 2017), which governs the relative efficiency of abrasion and plucking (Hancock et al., 1998; Whipple et al., 2000). The degree of fracturing further depends on the stress history that the rock experiences during and after its journey to Earth's surface (Leith et al., 2014; Eppes and Keanini, 2017). Chemical erosion susceptibility is set by the solubility of the in-channel bedrock and the saturation index of the stream water with respect to soluble minerals (e.g., Keen-Zebert et al., 2017).

Plentiful field examples show how bedrock erodibility can affect erosion process dominance and channel form. All else being equal, channels tend to be narrower in rock with greater tensile strength and therefore presumably lower erodibility (Bursztyn et al., 2015); the classic explanation for this is that channels in more resistant rock steepen and narrow to enable erosion at the rate of base-level fall (e.g., Yanites and Tucker, 2010). Channels tend to be wider, both in absolute terms and relative to their depth, in rock where closely spaced fractures make bedrock highly erodible by increasing the susceptibility to plucking (Spotila et al., 2015). Plucking is a much more efficient bedrock erosion process than abrasion per unit erosional power in channels with densely spaced fractures (Whipple et al., 2000). Channels may be wider and have lower gradients in these reaches because of their relatively high boundary erodibility (e.g., Wohl and Merritt, 2001; Yanites and Tucker, 2010; Chen and Byun, 2023; Buckley et al., 2024). Observations and theory agree that rock erodibility, which is governed by intact strength, fracture spacing, chemical erosion resistance, and weathering, influences bedrock river process and form, albeit in ways that are incompletely understood (Whipple et al., 2022).

The qualities of rock that govern its fluvial erodibility are inherited from, and contingent on, the outcomes of past geological systems (e.g., Gould, 1989; Phillips, 2021). Examples abound. Knickpoints in Hawaiian channels are collocated with outcrops of a type of lava flow that has relatively low fracture density and is therefore relatively erosion resistant (Raming and Whipple, 2022). Lava type is a function of flow dynamics and therefore the physical structure of the volcanic system (Rowland and Walker, 1990), such that the frequency and size of today's knickpoints are contingent on volcanic factors that may predate the river channel. Marshall and Roering (2014) described a landscape-scale example of erodibility inheritance in sedimentary bedrock. They showed that particular mineral assemblages and diagenetic conditions—the occurrence of fibrous clays and chlorite cements—are associated with relatively resistant, sparsely fractured beds compared to others in the same unit. Application of conceptually similar approaches to the reach-scale morphodynamics of bedrock rivers may enable assessment of how bedrock channels inherit process and form from the depositional environment of the channel-margin bedrock.

Here, we used a field case study to ask: (1) What rock properties seem to exert the most significant control over bedrock channel geometry? (2) How might those rock properties, and therefore how might channel process and form, be inherited from the depositional environment of the bedrock? Working in an Appalachian bedrock channel that crosses a major lithologic contact over a short distance, we combined measurements of rock properties and channel form with site-specific depositional environment interpretations from petrographic analysis of in-channel bedrock. We show that temporal changes in depositional environment, manifested as an up-section progression from an alluvial-fan to shallow-marine environment driven by marine transgression, translate directly into upstream-to-downstream changes in bedrock bed thickness, fracture density, erosion process dominance, and channel shape. We posit a causal link between depositional environment and modern channel form: The (un)steadiness of deposition sets bed thickness, which in turn influences fracture spacing and ultimately fluvial erodibility when the rock appears at Earth's surface as the boundary of a bedrock river.

We identified the Dry Fork, a river in the Allegheny Mountains of West Virginia, eastern United States (Fig. 1), as a promising natural laboratory in which to assess the ways in which the depositional environment of in-channel bedrock governs modern channel process and form. The Allegheny Mountains are tectonically quiescent and south of the Last Glacial Maximum extent. Knickpoints in the region are often lithologically controlled (e.g., Chilton and Spotila, 2022) as rivers cross different sedimentary units deposited in the Appalachian Basin (Bjerstedt and Kammer, 1988), indicating that local rock properties might exert a strong control on channel shape in this region. While hillslope-derived boulders can impose a nonlocal lithologic influence on channel form in steep landscapes like this one (Chilton and Spotila, 2020; Shobe et al., 2021), we did not observe deposits of hillslope-derived boulders in our study area.

The Dry Fork watershed lies immediately west of the Eastern Continental Divide, which separates basins draining to the Atlantic Ocean from those draining to the Gulf of Mexico, and flows to the Gulf of Mexico via the Cheat, Monongahela, Ohio, and Mississippi Rivers. It drains the Allegheny Mountain section, a transitional physiographic province connecting the folded rocks of the Appalachian Valley and Ridge to the flatter-lying units of the Appalachian Plateau. Mean annual precipitation derived from interpolation of precipitation gauge data from 1895 to 2023 is ~1330 mm (PRISM Climate Group, 2023). The Dry Fork is a fifth-order stream at our study sites, which were located ~5 km upstream from U.S. Geological Survey (USGS) discharge gauge #03065000 at Hendricks, West Virginia (drainage area 873 km²). Gauge records from 1940 through January 2024 record annual peak flows ranging from 202 to 2831 m³/s, with an average annual peak flow of 461 m³/s.

The Dry Fork watershed today is heavily forested with Appalachian hardwoods. Like most Appalachian rivers, the Dry Fork watershed experienced extensive deforestation during the nineteenth and early twentieth centuries due to timber harvesting and may also have hosted tie drives—downstream floating of railroad ties—that influenced flow, sediment, and erosion dynamics during that time (see, for example Wohl, 2001). A railroad operated in the Dry Fork valley until the mid-twentieth century (Hicks, 1963); the railroad grade appears to impinge on the channel in narrow reaches of the valley.

We studied two reaches, or short sections, of the Dry Fork between which the stream crosses the contact between two geologic units (Figs. 1 and 2). The upstream study reach, which we refer to as the sandstone reach, is ~160 m long and drains a 753 km² watershed (Fig. 2A). The downstream study reach, which we refer to as the carbonate reach, is ~100 m long and drains 783 km² (Fig. 2B). The two reaches are ~2.3 km apart, and a small tributary, Red Run, enters the Dry Fork from the east between the two reaches. Reach-averaged slopes calculated from 1-m-resolution, light detection and ranging (LiDAR)–derived digital elevation data are ~0.003 m/m in the sandstone reach and ~0.005 m/m in the carbonate reach.

The two study reaches lie on either side of the contact between the Upper Devonian–Lower Mississippian Price Formation (upstream, down section) and Mississippian Greenbrier Group (downstream, up section; Fig. 1B). The Price Formation (ca. 355–347 Ma) represents a deltaic complex with a mix of marine, deltaic, and alluvial plain deposits (Kammer and Bjerstedt, 1986; Bjerstedt and Kammer, 1988; Read and Eriksson, 2016). It exhibits great spatial variability in composition, but in our study area, it is an upward-coarsening siltstone and sandstone unit culminating in a quartz-pebble conglomerate that marks the top of the formation. The Greenbrier Group (ca. 340–325 Ma), a set of shallow-marine carbonates, unconformably overlies the Price Formation in our study area. It contains a wide diversity of carbonate facies (Rittenhouse, 1949; Read and Eriksson, 2012; see also Results and Discussion sections). Although to our knowledge no previous workers have made local paleoenvironmental interpretations for our study area specifically, the broad interpretation for the transition from Price to Greenbrier deposition in the region is one of marine transgression due to increased tectonic subsidence (Bjerstedt and Kammer, 1988; Read and Eriksson, 2012). Both units were folded as part of the North Potomac syncline, the hinge of which strikes NE-SW through the study area (Fig. 1B); dips reach up to ~10° in our study reaches.

The study area represents a space-for-time substitution in which temporal changes in ancient depositional environment over millions of years (down section to up section) are reflected in spatial changes (upstream to downstream) in channel-margin rock type along the modern Dry Fork. Although we primarily focused this study on the differences between the sandstone and carbonate bedrock units, each of these reaches has lithologic variability of its own. The upward-coarsening Price Formation contains a fine sandstone that dominates the upstream half of the sandstone reach and a coarse sandstone that underlies the downstream half (Fig. 2A). In the carbonate reach, we recognized four distinct subunits exposed along the channel margin (Fig. 3). We did not divide the carbonate reach into multiple sites because the carbonate channel margin is everywhere composed of multiple subunits. This sequence of depositional environments, now expressed as significant variations in channel-margin bedrock over short spatial scales, enables testing of relationships among depositional environment, channel-margin bedrock properties, and bedrock river process and form.

We investigated potential relationships among rock composition, strength, and discontinuity spacing, as well as channel geometry and erosion process dominance, to understand how changes in depositional environment during Mississippian time might govern the shape and erosion dynamics of today's channel. We used field measurements of channel cross-section form, bed thickness, and fracture spacing, laboratory measurements of rock strength, and hand-sample and thin-section petrography. The Price Formation and Greenbrier Group crop out across the region, vary in composition across space and depth, and have not been described at our study locations (e.g., Rittenhouse, 1949; Bjerstedt and Kammer, 1988; Read and Eriksson, 2012). Identification of minerals and fossils, along with estimation of their relative abundance, through thin-section petrography is essential to inferring the depositional environment of the rocks exposed in our study reaches and assessing connections between depositional history and modern bedrock erosion.

We used hand-sample and thin-section petrography to classify the different rock types exposed along the Dry Fork channel margin and to infer a depositional environment for each type. Hand samples were collected from the channel margin (Fig. 2) using a sledgehammer and chisel; Wagner Petrographic (Lindon, Utah, USA) prepared thin sections from the hand samples. We examined a total of 12 large-format (51 × 76 mm) thin sections: two from the fine sandstone reach, two from the coarse sandstone reach, and two from each of the four subunits we identified in the carbonate reach. Thin sections were partially stained with Alizarin Red S to facilitate distinction between calcite and dolomite. The thin-section epoxy is dyed blue, which can distinguish white or clear grains and cement from pore space. Thin sections were examined using an Olympus BX53 binocular microscope with a magnification range of 6.3–63× and transmitted, reflected, and polarized light. We identified and estimated relative abundances of visible minerals, cements, and fossils to establish the depositional environment of the rock in our study reaches specifically, rather than rely on the formation-scale interpretations for the relevant units available in past literature (Rittenhouse, 1949; Bjerstedt and Kammer, 1988).

We quantified hand-sample-scale rock strength using the point-load testing technique, which measures the load required to break centimeter-scale rock samples (e.g., Marshall and Roering, 2014; Chilton, 2021; Chilton and Spotila, 2022). The point-load index correlates with, and can be converted to, both compressive and tensile strength (Zhang, 2017), but it is best thought of as a relative strength index rather than an absolute strength value (Marshall and Roering, 2014; Chilton and Spotila, 2022). We followed the point-load testing procedures reported in Chilton and Spotila (2022; see their section 3.3.2 and their supporting information S1 for detailed testing procedures) and the International Society of Rock Mechanics guidelines (Franklin, 1985), as summarized in the following paragraph.

Samples were collected from the channel margin (Fig. 2) using a sledgehammer and chisel. Because it was not possible to closely control the size of samples chiseled from the bedrock, we used a rock saw to cut any samples too large to fit into the testing device to approximate cubes ~5 cm on a side. We submerged all samples in water for at least 4 d before testing to avoid strength differences due to varying water content. We then tested the samples by crushing them between the platens of a RocTest PIL-10 point-load testing device. Once failure occurred in the sample, we corrected the point-load index for sample size to yield a sample size–independent point-load index in units of megapascals (Chilton and Spotila, 2022). If breaking a sample yielded a fragment large enough to be tested as its own sample, we tested the strength of that fragment. We obtained 57 strength measurements from the sandstone reach and 121 from the carbonate reach. The greater number of measurements in the carbonate was purposeful because the carbonate unit had greater intra-unit heterogeneity (four identifiable subunits) than the sandstone unit (two subunits).

We measured bed thickness in the field using a tape measure. The limited subaerial exposure of bedding plane edges due to low banks and an absence of knickpoints precluded making large numbers of bed thickness measurements, but we were able to measure 17 beds in the fine sandstone unit, 18 beds in the coarse sandstone unit, and 25 beds in the carbonate unit. We collected images from uncrewed aerial vehicle (UAV) flights at each site. We processed the UAV images and global positioning system (GPS) data using the open-source Web Open Drone Map software to build orthomosaics with 2 cm resolution using photogrammetry. We then used the orthoimages to measure the fracture spacing of each rock unit. We quantified fracture spacing by drawing sampling lines totaling 100 m in length in the sandstone and 100 m in the carbonate and then measuring the distance along those lines between fractures that the lines crossed (Fig. 2). Fracture spacing measurements were intended to capture the major, throughgoing fracture sets visible in Figures 2 and 3. Fracture measurements made in rock exposed above the water surface were validated in the field with a tape measure; this was essential to capturing fractures spaced less than 4 cm (2 pixels in the orthomosaic) apart from one another, which could not be resolved from the imagery we collected. This procedure resulted in 52 fracture spacing measurements in the fine sandstone unit, 26 in the coarse sandstone unit, and 25 in the carbonate unit. We also collected the orientation of each intersecting fracture using a Brunton compass. There was not sufficient exposure of each carbonate subunit to make bed thickness and fracture spacing measurements for each, so we collected these measurements for the carbonate unit as a whole.

We used a Leica TS-7 total station to survey 10 channel cross sections in each of the three rock units (fine sandstone, coarse sandstone, and carbonate; Fig. 2). We surveyed across the active channel, which we defined as ending where exposed bedrock and alluvium gave way to soil-mantled forest floor. Our goal was to quantify the shape and surface roughness of each cross section to assess whether there were diagnostic differences in channel form that might have arisen from changes in erosion mechanism (e.g., plucking vs. abrasion) and/or efficiency due to lithologic characteristics. We discuss qualitative patterns in channel width and depth among cross sections but note that quantitative analysis of width and depth would ideally be based on cross-section surveys that span many channel widths in the downstream direction to average across morphologic patterns like pool-riffle sequences and meander bends. Limited bedrock exposure prevented collection of morphologic data at these length scales, so we focused on scales relevant to determining erosion mechanism and efficiency (~1–10 m). We used two techniques to combine the 10 surveyed cross sections from each lithology and quantify channel form at the sub-cross-section scale: analysis of cross-section hypsometry and cross-section roughness. To avoid potential bias in both analyses due to variable point spacing in the surveys, we first linearly interpolated our cross-section surveys to 0.1 m point spacing. We also tested monotonic cubic interpolation and found no significant differences in the results. This interpolation is reasonable given that we attempted to capture all significant observed topographic features when surveying.

Hypsometry at the scale of river cross sections is the frequency distribution of elevation along the surveyed channel boundary (e.g., Strahler, 1952; Finnegan and Balco, 2013). Hypsometric analysis can identify differences in cross-channel shape, for example, between channel cross sections that are V-shaped and those that might have a thalweg incised below flat portions of the channel bed (Finnegan and Balco, 2013, their fig. 2). We binned the elevation values of all 10 interpolated cross sections from each rock unit into 0.1 m increments to yield a single distribution of elevation frequency within channel cross sections in that rock unit, relying on the combined data to average across the variability from individual cross-section surveys. For cross sections with surveyed elevations that differed between the river-left and river-right banks (typically due to difficulty in identifying the edge of the active channel or other field constraints), we truncated the cross sections such that only elevations represented on both banks were considered. This analysis resulted in three elevation frequency distributions—one for each unit—that could be compared statistically to identify whether there were quantitative differences in the distribution of elevation above the thalweg among lithologies.

The roughness of channel cross sections might indicate the erosion processes acting on the channel margin. All else being equal, we might expect channels eroding primarily through the plucking of fracture-bounded blocks of rock to exhibit greater cross-channel roughness than those eroding through other mechanisms such as abrasion (e.g., Spotila et al., 2015, their fig. 7B) or dissolution. To quantify differences in cross-section roughness, we calculated the frequency of inflection points (transitions from positive to negative cross-stream slope or vice versa) in our cross sections. We then resampled the interpolated cross sections at point spacings ranging from 0.1 m to 10 m, each time counting the frequency of inflection points (units of m^–1) in the resampled cross section. This method captures the effects of both small (decimeter-scale) and large (decameter-scale) roughness elements. When the cross sections are resampled at dense point spacing, inflection points that occur due to small-scale roughness dominate, while only inflection points due to large-scale roughness remain at sparser point spacing. Roughness values were averaged across the 10 cross sections surveyed in each rock unit.

Fine sandstone reach. The fine sandstone samples of the Price Formation are dominated by quartz grains, which range from angular to rounded and exhibit close contacts between grains in some locations (Figs. 4A, 5A, and 5B). Isopachous intergranular hematite cement surrounds the grains. Organic material occurs in flecks throughout samples and is concentrated in some intergranular spaces. There is little evident porosity. The fine sandstone unit coarsens upward at the centimeter scale, as seen in the two large-format thin sections we analyzed, from a very fine to a fine sandstone, and the proportion of quartz grains to lithic mudstone grains also increases up section. The stratigraphically lower sample is well sorted, while the upper sample is moderately well sorted. The upper sample also contains more intergranular hematite and clay cements. First-generation hematite cements and the presence of organics may suggest subaerial deposition. The fine sandstone is a lithic arenite sandstone that may have formed in an alluvial fan or fluvial system.

Coarse sandstone reach. The coarse sandstone unit of the Price Formation consists of moderate- to high-sphericity, subangular to subrounded, medium sand to pebble quartz grains, medium sand-sized lithic mudstone grains, and very well-rounded coarse sand to pebble quartzite grains (Figs. 4B, 5C, and 5D). Minor feldspar grains make up ~5% of the grains. Quartz overgrowth cement surrounds ~30% of quartz grains. Hematite cement is common and forms isopachous coatings on the remaining grains as well as their overgrowths. Green clay cement is minor and most visible where mudstone lithic grains have been plucked from the sample during thin-section preparation; clay cement appears to intermingle with hematite coatings. Hematite cement completely fills the few intergranular void spaces and forms well-developed coatings on grain surfaces. Close grain contacts are common, especially between quartz matrix grains. Dark-brown organic material is also common and concentrated at grain boundaries. Rock porosity is low: <5% of the thin section area. The reworked subangular-subrounded, moderately sorted grains with hematite cement are consistent with subaerial deposition by fluvial and/or alluvial processes. Pits in thin sections indicate that mudstone lithic grains were plucked during thin-section preparation, suggesting weak cementation. The similar size and rounding of the quartz and mudstone sand grains indicate a separate source from the well-rounded quartzite grains. The coarse sandstone unit is a pebbly lithic arenite.

The continuous changes in grain size and composition from the lower portion of the fine sandstone unit through the upper portion of that unit and culminating in the overlying coarse sandstone unit probably indicate that these uppermost sandstones of the Price Formation represent a single depositional environment. We interpret that the two sandstone units represent a subaerially deposited alluvial fan.

We observed four distinct strata exposed within the carbonate reach (Fig. 3) and labeled these subunits C1–C4 from bottom to top of the exposed section. The C2–C4 subunits are of roughly equal thickness (0.5–1 m); the thickness of C1 could not be measured because its lower boundary was not exposed. Conformable contacts between units are gradational over several centimeters.

The bottommost strata, C1, is a mottled gray and brown, poorly sorted, siliciclastic-carbonate graywacke (Figs. 4C and 6A). It is composed of sand- to pebble-sized carbonate lithic fragments in a quartz silt matrix. The quartz silt grains are well sorted, angular, and of moderate sphericity. Sand- and pebble-sized carbonate intraclasts consist of carbonate mud and skeletal fragments of crinoids, brachiopods, sponges, and trilobites (Figs. 4C and 6A). Diagenetic features include calcite pendant cements and recrystallized crinoid columnals. The poorly defined grain boundaries on the larger carbonate lithic grains suggest only partial cementation prior to subsequent transport, and the fine grain size of the siliciclastic material compared to the underlying Price Formation suggests transport and reworking of eroded material prior to deposition. We interpret C1 to be the basal carbonate unit of the Greenbrier Group, formed as the underlying Price Formation clastics were eroded. Due to its poor sorting, presence of organic material, lack of hematite cement in clastic grains, and location within the stratigraphic framework of the area, we interpret C1 as a likely shallow fluvial or tidal channel deposit.

Unit C2 overlies C1 and has a higher proportion of carbonate material. C2 contains ~60%–70% skeletal grains, 20% peloids, and 10%–20% quartz silt (Figs. 4D and 6B). Like C1, C2 hosts a diverse fossil assemblage including gastropods, trilobites, echinoderms, brachiopods, algae, and sponges. Some skeletal grains have been partially dissolved and recrystallized. There are two generations of blocky calcite intergranular cement. Geopetal fill in gastropods and intragranular pores of other skeletal grains consists of fossils, mud, and quartz silt. Geopetal fill indicates some degree of reworking; this is supported by the presence of skeletal hash. The presence of fine quartz silt and absence of other clastic material point to an environment with selective addition of clastic material, possibly transported by wind from an already-sorted source such as a beach. We interpret C2 as a skeletal packstone/grainstone deposited in an open lagoon or shallow shelf.

Unit C3 is a fossiliferous wackestone (Figs. 4E and 6C) composed of brachiopod and dasycladacean algae fragments, as well as peloids and carbonate mud. Quartz silt grains comprise ~5%–10% of thin-section area. The skeletal fragments likely indicate high-energy transport, but the carbonate mud and peloids may point to deposition in quiescent waters. We interpret the C3 unit as a shallow-marine shelf deposit.

Unit C4, the uppermost strata of the carbonate unit in the field area, is composed primarily of ooids and coated grains with subrounded quartz silt and recrystallized skeletal fragments, including gastropods and sponge spicules, at their cores (Figs. 4F and 6D). Small numbers of uncoated grains, including peloids, crinoid fragments, and uncoated quartz silt grains, were also observed. The grains are cemented by carbonate intergranular cement. We call C4 an ooidal packstone and infer that it was deposited on a shoal at the edge of a shallow-marine carbonate platform.

From upstream to downstream and down section to up section in our study area, the bedrock of the modern Dry Fork channel margin transitions from sandstones deposited in an alluvial fan or fluvial system to carbonates deposited in marginal marine and marine settings. We interpret this depositional sequence to represent a marine transgression, consistent with interpretations from past work (e.g., Bjerstedt and Kammer, 1988).

We made 174 point-load strength measurements, 57 distributed between the coarse and fine sandstones and 117 distributed among the four carbonate subunits (Fig. 7). Using Kruskal-Wallis H tests with Dunn's post hoc tests and Bonferroni corrections for multiple comparisons implemented in the scikit-posthocs Python package (Terpilowski, 2019), there were few significant differences in strength among the individual subunits we sampled. Significant differences were mainly caused by the high strength of the fine sandstone unit, which was significantly stronger than the C1 and C3 carbonate subunits (p = 0.001 and p < 0.001, respectively). The coarse sandstone may be stronger than carbonate subunit C3 (p = 0.022), but not than any other carbonate subunits. The carbonate subunits were not different from one another except for possibly C3 and C4 (p = 0.038); the two sandstones did not have significantly different strengths (p > 0.05). When the carbonate units were combined into a single data set—a reasonable approach given that multiple subunits underlie each of our carbonate cross sections—the fine sandstone was stronger than the carbonate (p < 0.001), but strength differences between the fine and coarse sandstone and between the coarse sandstone and carbonate were not significant (p > 0.05).

Combining all carbonate strength data and comparing them against the combined coarse and fine sandstone data revealed significant differences between the two units overall (Mann-Whitney U test; p < 0.001). Even for this aggregated comparison, however, differences between median point-load strength indices for the two populations were only ~1 MPa, which is fairly minor relative to (1) the variability in strength within each unit we sampled (Fig. 7) and (2) intra-unit strength variability found by a prior study in this region (Chilton and Spotila, 2022).

We investigated two types of discontinuities in the rock units in our study area: bedding planes and fractures (Fig. 8). Bed thickness averaged ~45 cm in the carbonate units (range of 3–200 cm; n = 25), which to the extent we could see did not vary consistently among the C1–C4 subunits. Beds were thinnest in the coarse sandstone unit with an average thickness of ~3 cm (range of 1–10 cm; n = 18). The fine sandstone unit had intermediate bed thickness, with an average thickness of ~10 cm (range of 2–21 cm; n = 17). There were significant differences in bed thickness among the three groups: Post hoc tests showed that the carbonate beds are significantly thicker than the coarse sandstone beds (p < 0.001) and may be significantly thicker than the fine sandstone beds (p = 0.036). Differences in bed thickness between the two sandstone units were not significant (p > 0.05). When measurements were combined for the two sandstone reaches, carbonate bed thickness was significantly greater than sandstone bed thickness (p < 0.001).

Measurements of the spacing of fractures that were perpendicular to the bedding revealed that the carbonate unit hosts two distinct fracture sets (visible in Fig. 2): one in subunit C1, typically exposed at the thalweg in our study reach, and one in subunits C3 and C4, exposed along the banks. Average fracture spacing for these two sets was 178 cm (range of 53–278 cm; n = 10) and 313 cm (range of 108–707 cm; n = 15), respectively (Fig. 8). Average fracture spacing was 44 cm (range of 3–200 cm; n = 26) in the coarse sandstone and 96 cm (range of 24–268 cm; n = 52) in the fine sandstone (Fig. 8). There were significant differences among the four groups. The bank assemblage of carbonate fractures had greater spacing than both the fine (p < 0.001) and coarse (p < 0.001) sandstone fractures. The thalweg assemblage had significantly greater fracture spacing than the coarse sandstone (p < 0.001) but not than the less densely fractured fine sandstone (p > 0.05). The fine sandstone had greater fracture spacing than the coarse sandstone (p = 0.002), while the two carbonate assemblages did not have significantly different spacing (p > 0.05). When combining the two carbonate assemblages to get a mean fracture spacing for the carbonate unit as a whole, each of the three rock units had different fracture spacing from the others (fine sandstone greater than coarse sandstone: p = 0.001; carbonate greater than fine sandstone: p < 0.001; carbonate greater than coarse sandstone: p < 0.001). Combining the two sandstone units and the two carbonate fracture assemblages each into a single data set revealed that the carbonate had significantly greater fracture spacing than the sandstone (p < 0.001). In addition to differences in fracture spacing, fracture orientation was much more consistent in the carbonate units than the sandstone units (Fig. 9).

Cross sections in the carbonate reach have a U-shape in which a relatively flat thalweg is bounded by steep banks (Fig. 10). Sandstone reaches display in some cases an inner-gorge morphology in which a narrow, V-shaped notch is inset into the broader cross section and in some cases a flat but rough channel bottom (Fig. 10). Sandstone channels tend to be wider and shallower, while carbonate cross sections are narrower and deeper. However, analysis of width and depth must remain qualitative given the restricted bedrock exposure in our study area, which precluded us from measuring cross sections along a distance of many channel widths in each lithology. We attempted to quantify differences between cross-section morphology in the different rock units at the erosion process scale by using two morphologic metrics: cross-section hypsometry (distribution of elevations above the thalweg) and cross-section roughness (frequency of inflection points).

The carbonate cross sections showed a peak in elevation frequency within 0.5 m of the thalweg and a decline in frequency with additional height above the thalweg (Fig. 11A). A minor peak in frequency occurred at ~3 m above the thalweg, where surveys captured the transition to the top of the channel banks (Figs. 10 and 11A). Coarse sandstone cross sections showed an initial peak between 0 and 0.1 m and a second, broad peak centered on 0.5–0.6 m, followed by a decline in frequency with height (Fig. 11B). Fine sandstone cross sections showed a bimodal hypsometric distribution, having a broad peak in elevation frequency centered on 0.3–0.4 m and a minor peak centered on 1.2–1.3 m (Fig. 11C). The carbonate elevation frequency distribution was significantly different from the coarse sandstone distribution (p < 0.001) and possibly the fine sandstone distribution (p = 0.028). The two sandstones showed significantly different elevation frequency distributions (p < 0.001). The hypsometry of the carbonate cross sections differed from that of the sandstone cross sections when the fine and coarse sandstone data were lumped together (p < 0.001).

Cross-section surveys (Fig. 10) suggest that the sandstone unit might generate rougher channel boundaries than the carbonate unit. Inflection-point analysis showed that the sandstone units have more frequently occurring inflection points across a wide range of sampling intervals (Fig. 12). The plateau in the roughness data at decimeter-scale sampling intervals represents the inflection frequency in the raw data because the sampling interval was less than or equal to the average point spacing in our surveys. The sandstone inflection-point data sets inevitably collapse toward the carbonate data set as the sampling interval grows toward a value approximating the channel width. The increased roughness of the sandstone cross sections is most pronounced at sampling intervals below ~3 m, consistent with the qualitative impression given by the cross sections (Fig. 10).

Channel form changes from upstream to downstream as the Dry Fork crosses from sandstone to carbonate rock. In the fine sandstone, and in the coarse sandstone immediately downstream, we observe wider, shallower cross sections with broad, flat beds (Fig. 10) that have peaks in elevation density well above the thalweg (Fig. 11) and rough bedrock boundaries (Fig. 12). The carbonate units only a couple of kilometers downstream show narrower, deeper, U-shaped cross sections (Fig. 10) with monotonic declines in elevation density with height above the thalweg (i.e., no incision of the thalweg below the remainder of the bed; Fig. 11) and smoother bedrock boundaries across decimeter to meter scales (Fig. 12). We also observe steeper slopes in the carbonate reach relative to the sandstone reach. These two distinct morphologies are broadly consistent with those observed by others working in Appalachian bedrock rivers—Spotila et al. (2015) named them “incision plain” and “channel neck,” respectively (their fig. 11)—with the caveat that our study sites and theirs differ substantially in scale (their sites along the New River drain approximately an order of magnitude more area than ours along the Dry Fork).

The differences in morphology as the channel crosses between rock types could be due to any of the three facets of bedrock erodibility: the strength of the unfractured rock, the degree of fracturing (Whipple et al., 2022), or the rock's susceptibility to chemical erosion. Rock strength measurements show that, at the scale of samples testable with the point-load method (~5 cm), the sandstone units are generally stronger than the carbonates. However, differences in point-load index between the two units (~1 MPa) are modest given the range of variability in Appalachian sedimentary rocks (>2.5 MPa; Chilton and Spotila, 2022). Further, we would not expect wider, shallower, flat-bottomed channels in stronger rock and narrower, deeper, U-shaped channels in weaker rock (e.g., Wohl and Merritt, 2001; Yanites and Tucker, 2010; Chen and Byun, 2023; Buckley et al., 2024), especially in relatively low-sediment-supply conditions like those prevailing in Appalachian rivers. We therefore suggest that, as in other Appalachian rivers (Chilton and Spotila, 2022), the strength of the unfractured bedrock is not the primary control on channel morphology in this case.

Like Spotila et al. (2015), we interpret the change from incision plain to channel neck morphology to occur due to a change in erosion process dominance. Our sandstone reaches appear to be plucking-dominated; jagged edges of bedding planes and block-shaped voids are common, and there are few flutes or other abrasion forms (Figs. 2C, 2D, and 13). The carbonate reach appears to be dominated by, or at least to reflect greater contributions of, abrasion and chemical erosion. We observed rounding of large, in situ bedrock blocks and pitting that probably results from dissolution (Fig. 3B). Measurements of bed thickness and fracture spacing support this interpretation; sizes of potentially pluckable blocks are centimeters on a side in the sandstones but several decimeters on a side in the carbonates. Because plucking tends to be the most efficient bedrock erosion mechanism when blocks are sufficiently small (Whipple et al., 2000), the Dry Fork may be able to maintain lower slopes and a wider, shallower, flat-bottomed cross section in the sandstone reaches, where flows frequently exceed the threshold to pluck centimeter-scale blocks, while steepening and excavating a deeper, narrower cross section in the carbonate reaches, where blocks are an order of magnitude larger and only pluckable during rare floods. This interpretation is consistent with many studies that have investigated the influence of bedrock properties on channel form (e.g., Scott and Wohl, 2019, and references therein) but does not explain the genesis of those properties. We used our petrology-derived depositional environment interpretations to ask: What factors, in the sequence of geologic events leading up to the appearance of these rock units in the bed of the Dry Fork, control the spacing of discontinuities that we now observe to be governing channel process and form?

Because discontinuity spacing appears to be the key rock property influencing channel morphology in our study area, we hypothesize that the bed thickness inherited from the depositional environment of the sedimentary bedrock controls modern process and form. Much prior work has established that bed boundaries are discontinuities (e.g., Shang et al., 2016) along which blocks may be plucked (Fig. 2C; Scott and Wohl, 2019; Chilton and Spotila, 2022), and that bed thickness controls the spacing of fractures that develop perpendicular to the bedding during deformation (e.g., Ladeira and Price, 1981; Gross et al., 1995; Bai and Pollard, 2000; Chemenda, 2022). All else remaining equal, thinner beds tend to develop more densely spaced fractures than thicker beds develop. This relationship can be nonlinear; fracture spacing is most sensitive to bed thickness when beds are thin (decimeter scale) but becomes less sensitive to bed thickness when beds are a meter or more thick (e.g., Chemenda, 2022, their fig. 1). Although the quantitative form of this relationship, to our knowledge, has not been established for the units we investigated, it is likely that the thin (centimeter-scale) beds of the sandstone units allowed the development of more densely spaced fractures (decimeter-scale fracture spacing) than did the thicker (decimeter-scale) beds of the carbonate units (meter-scale fracture spacing) because fracture density is very sensitive to bed thickness for thin beds. The sandstones may therefore have become more densely fractured than the carbonates as both units experienced tectonic stresses from the folding of the Allegheny Mountain section. It is also possible that denser fractures in the sandstone relative to the carbonate arise in part because fracturing can occur at lower strain in rocks with greater stiffness (Gross et al., 1995), a property that tends to correlate with rock strength (D'Andrea et al., 1965).

Differences in channel form between rock units in our study area arise from changes in the relative efficacy of plucking, which depends on discontinuity density and therefore bed thickness. Tracing the origins of modern bedrock channel form to the bedrock's depositional environment requires understanding why bed thickness differs between units. Our hypothesis is that bed thickness in the units we studied is largely a function of the dynamism, or lack thereof, of the depositional environment. Petrologic evidence suggests that the sandstone units were deposited in an alluvial-fan or other fluvial environment that might have experienced rapid channel migration and/or frequent avulsions (e.g., Leenman and Eaton, 2021). Each bed might therefore only represent a relatively short time between changes in the location of deposition. Conversely, the petrology of the carbonates suggests a marginal marine to marine depositional environment with water depth progressively increasing from C1 to C4. In these environments, deposition at any point may have been more continuous, resulting in thicker beds and therefore lower fracture density.

Investigation of how micron- to meter-scale rock properties relate to reach-scale bedrock channel form and erosion processes reveals one potential pathway by which the depositional environment of in-channel bedrock might influence modern river dynamics (Fig. 13). We suggest that the depositional environment leaves its imprint on the modern landscape in a two-step chain of causality (Fig. 13). The relative dynamism of the depositional environment sets bed thickness, which controls the development of the fracture network as rock units undergo tectonic stresses. Fracture density and bed thickness then govern the relative ease of plucking, which controls bedrock erosion efficiency and channel shape.

Our results expand upon past work suggesting that bed thickness in sedimentary rock represents an important link between the rock's depositional environment and its resistance to fluvial erosion (Spotila et al., 2015; Chilton and Spotila, 2022). We have hypothesized that bed thickness might exert an additional control on fluvial processes and channel geometry by influencing the extent to which tectonic stresses precondition bedrock for erosion through fracturing (e.g., Molnar et al., 2007; Roy et al., 2016). These findings complement previously documented diagenetic controls on rock strength, fracture spacing, and rock erodibility (Marshall and Roering, 2014) and suggest that it may be possible to establish a priori quantitative relationships that allow prediction of geomorphic process variability from widely available sources of information about bedrock properties (e.g., geologic map databases).

To connect bedrock depositional history with erosional process and form in bedrock channels, we studied a river reach crossing a lithologic contact that represents a transition from alluvial-fan to shallow-marine deposition. Our case study shows potential causal links across >300 m.y., relating modern channel form to the depositional environment of the in-channel bedrock. Beginning at the modern channel margin and working backward in time, we find that:

1. Discontinuity spacing, rather than hand-sample-scale rock strength, may set modern bedrock erosion process dominance and channel form. This is consistent with past work in suggesting a scaling break—at the discontinuity spacing below which channels become plucking-dominated—distinguishing a regime in which channel form reflects the strength and/or solubility of the intact rock from one in which channel form reflects discontinuity spacing.

2. Denser discontinuities preferentially occur in the sandstone units we studied over the carbonates, probably due to the thinner beds in the sandstones. Given the known positive relationship between bed thickness and spacing of fractures in rock, we suggest that the sandstones host denser discontinuities both because the bedding planes themselves are discontinuities and because tectonic stresses develop denser fracture networks in thinly bedded rocks.

3. Depositional environment interpretations from thin-section petrography suggest that the sandstones were deposited in a more dynamic, unsteady depositional environment (alluvial fan) than the carbonates (shallow marine). The change in depositional environment may have caused the differences in bed thickness, and therefore discontinuity spacing, erosion process dominance, and channel form, that we observe between the two units.

Given the wide variability in the geomechanical and chemical properties of the rocks currently eroding at Earth's surface, the specific mechanisms by which the history of channel-margin bedrock is imprinted on modern river form are likely highly variable. Our results may apply most to landscapes where, like the Appalachian Plateau, rivers incise into relatively lightly deformed sedimentary bedrock in a tectonically quiescent setting. This study and others like it suggest the possibility that, in such landscapes, knowledge of bedrock history gained from geologic maps or similar data sets could be used to predict variability in geomorphic processes. Rather than viewing the properties of Earth-surface materials as boundary conditions, we can think of them as a rich inheritance from past geologic processes that sets the template for today's surface processes, ecosystems, and natural resources.

The data and code underlying this study are permanently archived in a public repository available at https://doi.org/10.5281/zenodo.13840741 (Colaianne and Shobe, 2024). The findings and conclusions in this publication are those of the authors and should not be construed to represent any official U.S. Department of Agriculture or U.S. government determination or policy. Any use of trade names is for descriptive purposes only and does not imply endorsement by the U.S. government. This study was funded by a National Aeronautics and Space Administration West Virginia Space Grant Consortium graduate fellowship and a West Virginia University (WVU) Geology graduate student grant (both to N.J. Colaianne), a WVU Research and Scholarship Advancement award (to C.M. Shobe), and a National Science Foundation EAR postdoctoral fellowship (to K.D. Chilton). Thanks go to Chris Belland, Sam Bower, Maya Bradford, and Dietrich Kuhlmann for field and laboratory assistance, and to Jaime Toro for helpful input. Julia Carr, one anonymous reviewer, the associate editor, and Science Editor Andrea Hampel provided perceptive comments on the manuscript.