当前位置：
X-MOL 学术
›
Lithosphere
›
论文详情
Our official English website, www.x-mol.net, welcomes your
feedback! (Note: you will need to create a separate account there.)
Quantifying Normal Fault Evolution from River Profile Analysis in the Northern Basin and Range Province, Southwest Montana, USA
Lithosphere ( IF 1.8 ) Pub Date : 2021-01-08 , DOI: 10.2113/2021/7866219 Ian P. Armstrong 1 , Brian J. Yanites 1 , Nate Mitchell 1 , Clarke DeLisle 1 , Bruce J. Douglas 1
Lithosphere ( IF 1.8 ) Pub Date : 2021-01-08 , DOI: 10.2113/2021/7866219 Ian P. Armstrong 1 , Brian J. Yanites 1 , Nate Mitchell 1 , Clarke DeLisle 1 , Bruce J. Douglas 1
Affiliation
Over the past few decades, tectonic geomorphology has been widely implemented to constrain spatial and temporal patterns of fault slip, especially where existing geologic or geodetic data are poor. We apply this practice along the eastern margin of Bull Mountain, Southwest Montana, where 15 transient channels are eroding into the flat, upstream relict landscape in response to an ongoing period of increased base level fall along the Western North Boulder fault. We aim to improve constraints on the spatial and temporal slip rates across the Western North Boulder fault zone by applying channel morphometrics, cosmogenic erosion rates, bedrock characteristics, and calibrated reproductions of the modern river profiles using a 1-dimensional stream power incision model that undergoes a change in the rate of base level fall. We perform over 104 base level fall simulations to explore a wide range of fault slip dynamics and stream power parameters. Our best fit simulations suggest that the Western North Boulder fault started as individual fault segments along the middle to southern regions of Bull Mountain that nucleated around 6.2 to 2.5 Ma, respectively. This was followed by the nucleation of fault segments in the northern region around 1.5 to 0.4 Ma. We recreate the evolution of the Western North Boulder fault to show that through time, these individual segments propagate at the fault tips and link together to span over 40 km, with a maximum slip of 462 m in the central portion of the fault. Fault slip rates range from 0.02 to 0.45 mm/yr along strike and are consistent with estimates for other active faults in the region. We find that the timing of fault initiation coincides well with the migration of the Yellowstone hotspot across the nearby Idaho-Montana border and thus attribute the initiation of extension to the crustal bulge from the migrating hotspot. Overall, we provide the first quantitative constraints on fault initiation and evolution of the Western North Boulder fault, perhaps the farthest north basin in the Northern Basin and Range province that such constraints exist. We show that river profiles are powerful tools for documenting the spatial and temporal patterns of normal fault evolution, especially where other geologic/geodetic methods are limited, proving to be a vital tool for accurate tectonic hazard assessments.The evolution of topography within mountainous landscapes is driven by the competition between two major processes: (1) the tectonic uplift of rock [1] and (2) climatically and lithologically controlled erosion [2–4]. When a tectonic perturbation like normal faulting causes a relative drop in the base level, relief is generated, and bedrock rivers respond by incising into the landscape undergoing relative uplift. The transition from the undisturbed upstream relict channel to the actively incising (transient) channel is signified by a discontinuity in the channel slope, known as a knickpoint or knickzone. The spatial distributions of knickpoints and their migration rates are strongly linked to rock uplift rates and therefore act as an archive of tectonic activity within a landscape [5–15].In the past few decades, the widespread study and improvement of bedrock channel morphology analysis, cosmogenic nuclide dating, and river incision modeling techniques have produced insight into transient channel behavior. Combined with the increasing availability of high-resolution digital elevation models (DEMs), these advances have provided the data necessary to extract tectonic information from topography around the globe [4]. For example, Ellis et al. [5] found that multiple series of knickpoints upstream of normal faults in the Pine Forest, Santa Rosa, and Jackson Ranges in Northwestern Nevada were attributed to ~220 m of base level drop since the Pliocene. Boulton and Whittaker [16] found that knickpoints upstream of normal faults in the Hatay Graben in southern Turkey were formed by a 5-fold acceleration in slip rates, and that the half-graben is now bounded by a soft-linked normal fault, proving to have substantial hazard implications for the region. Furthermore, Gallen and Wegmann [17] found that the spatial and temporal changes in bedrock channel steepness (ksn) are consistent with the growth and linkage of the large Ptolemy and South-Central Crete normal fault systems, validating the findings of numerous other studies that show that river profiles will adjust their steepness to keep pace with the rate of rock uplift [13, 18–20]. These studies are just a few of many demonstrating that bedrock channel analyses successfully provide insights into the relative magnitude and timing of tectonic perturbations, even where existing geologic/geodetic constraints are poor [4, 5, 9–11, 14, 15, 17, 20–24].Despite the growing amount of literature using this emerging approach, few studies on the detailed evolution of individual normal fault systems, particularly in the northern end of the Basin and Range province, Southwest Montana, exist. This is due to the fact that the relatively young, slow moving normal fault systems in this region are often inadequately constrained with conventional techniques like thermochronology, which require large amounts of fault throw to exhume thermally reset minerals [25], or with geodetics, which require relatively rapid fault slip rates [26]. Furthermore, thermochronology and geodetics have temporal resolutions of 106-107 and 100-101 years, respectively, hindering the ability to constrain fault slip at the critical time scales of >105-106 years. However, by analyzing bedrock channel morphology and applying stream power incision models [1, 10, 11, 27–29], transient landscapes offer an excellent opportunity to constrain the timing and pattern of base level fall along young fault systems with low slip rates and total throw, even when other geomorphic or geodetic data is lacking [16].This paper aims to use bedrock channel analysis to place constraints on the spatial and temporal slip rates across the Western North Boulder normal fault in the northern end of the Basin and Range province, Southwest Montana (Figures 1(a)–1(c)). Fault slip data in this low strain rate, seismically active region, are sparse and/or poorly constrained, limiting confidence in seismic hazard assessments [30]. Furthermore, the young, actively extending Western North Boulder fault has relatively little fault throw (Johns et al., 1982; Stickney and Bartholomew, 1987), providing an opportunity to observe the beginning/intermediate phases of normal fault evolution and how the surrounding landscape responds in such a scenario. Normal fault mechanics predict elliptical patterns of slip along strike with higher slip rates along the central portion of the fault relative to the distal portions, where the magnitude of slip scales with displacement [5, 9, 17, 30–35]. This serves as a fundamental prediction we may see in our study area, making this region a good place to perform tectonic geomorphology analyses. Moreover, this gradient in tectonic activity provides an opportunity to calibrate and test a stream power model of river incision [36]. To accomplish our analyses, we use (1) channel morphology metrics including channel steepness values, knickpoint distributions, and transient channel incision depths; (2) cosmogenic erosion rates; (3) bedrock strength data; and (4) calibrated stream power incision models to simulate the evolution of transient channels during a sustained period of base level fall. Doing so will shed light on the spatial and temporal dynamics of active Quaternary normal faults within the Northern Basin and Range region, an area where slip rate data are sparse, ambiguous, and/or poorly constrained, therefore limiting seismic hazard assessments. Furthermore, this study will highlight the utility of tectonic geomorphology in constraining fault evolution in young systems with low slip rates and total throw.The Northern Basin and Range province in Southwest Montana is home to a network of semi-isolated, fluvially connected intermontane basins bounded by Quaternary normal faults [37, 38]. Our study area is along the eastern margin of Bull Mountain and the adjacent North Boulder Basin (Figures 1(a)–1(c)). The flanking highlands of these intermontane basins, such as Bull Mountain, are remnants of the Cordilleran fold-and-thrust belt [38] which were primarily formed under three tectonomagmatic events that overlap in space, namely, (1) Sevier fold and thrust belt deformation, (2) calc-alkaline magmatism, and (3) Laramide basement uplifts [37–39]. The basement-cored Laramide uplifts and Sevier fold-and-thrust belt development are the two primary tectonic components responsible for the development of the Cordilleran thrust belt that began in the late Jurassic and continued its development into the late Paleocene to early Eocene [38]. Widespread arc volcanism and emplacement of Cordilleran batholiths also occurred during the late cretaceous to early Eocene time [38, 40, 41], creating the Boulder Batholith (81-73 and 64-61 Ma) and Elkhorn Mountain volcanics (83-80 Ma) that compose Bull Mountain [38, 42–44]. Acting as the roof of the Boulder Batholith, the Elkhorn Mountain volcanics are primarily composed of ignimbrites, dacites, rhyolites, and volcaniclastic sedimentary rocks, while the Boulder Batholith is primarily composed of intrusive bodies of granodiorite and quartz monzonite [41].Following the Sevier fold-and-thrust belt development and Laramide uplifts, regional compression was followed by extension, thought to have initiated by the gravitational collapse of the Cordilleran orogenic wedge that affected Southwest Montana around 45 Ma and created a network of semi-isolated, fluvially connected N-S trending intermontane basins [38, 45–50]. Following the Eocene extension, a period of quiescence during the Oligocene resulted in relatively little tectonic activity and a gap in deposition, which ended during the Late Neogene to Quaternary with the onset of an increase in the extension rate that has continued into the present [5]. Previous work suggests that this increase in Late Neogene to Quaternary normal faulting and seismicity is associated with the crustal bulge from the passage of the Yellowstone hotspot along the nearby Montana-Idaho border, beginning around 6-2 Ma and continuing into the present (Figure 1(c)) [51–55]. This ongoing phase of increased extension further reactivated and segmented the intermontane basins within the Northern Basin and Range province, creating the regional topography observed today (Figure 1(c)) [38, 40, 52].Studies have suggested that slip along the Western North Boulder fault (Figures 1(a)–1(c)) has been ongoing at a rate of <0.2 mm/yr since at least 0.75 Ma and, as a result, induced the incision and deposition of alluvial fans that resulted in as much as 120 m of offset on the Western North Boulder Basin alluvial fan surface (Figure 1(c)) (Johns et al., 1982; U.S. Geological Survey and Montana Bureau of Mines and Geology, Quaternary fault and fold database for the United States). These fans inhibit the persistence of any fault scarps or footwalls on ≥10 m DEMs; therefore, the timing and magnitude of fault slip along the Western Boulder fault are not well constrained. Regardless, combined studies constrain the initiation of increased extension around 6 to 0.75 Ma ([51–53]; Johns et al., 1982; [54, 55]). We aim to analyze the transient channels of eastern Bull Mountain to improve constraints on the spatial and temporal slip rates across the Western North Boulder fault zone and help characterize Northern Basin and Range normal faulting histories. Furthermore, we strive to illustrate the utility of tectonic geomorphology in constraining normal fault evolution in young systems where slip rates and total throw are too low to extract via other geologic methods.Transient landscapes in the Northern Basin and Range province in Southwest Montana provide a record of the region’s extensional tectonic history [56, 57]. Fifteen channels draining along the Western North Boulder fault exhibit a striking contrast in topographic gradient, with a gentle, low relief relict landscape situated upstream at higher elevations and a steep, high relief transient landscape situated at lower elevations (Figures 1(a) and 1(b)). Because the hillslope gradients increase across this transition and show no evidence of glaciation or stream capture that could produce nontectonic pulses of incision [58, 59], we interpret these morphologies to be indicative of a headward-advancing landscape adjustment to a change in the rate of base level fall. Furthermore, the proximal Western North Boulder fault has been actively extending the North Boulder Basin since at least 0.75 Ma (Johns et al., 1982), and this tectonic activity could change the rate of base level fall for these streams. We therefore argue that the transient incision within these channels provides a record of the ongoing extension along the Western North Boulder fault. Channels that do not preserve clear evidence of transience are not shown in Figures 1(a) and 1(b) to maintain visual clarity. For example, our study area does not extend to the southern portion of the Western North Boulder fault and Bull Mountain because no transience is preserved in this area, which could be due to variations in erosional processes, drainage area, bedrock properties, or fault dynamics. Along the 15 channels where transience is ongoing, we aim to improve constraints on the spatial and temporal slip rates across the Western North Boulder fault zone by applying channel morphometrics, cosmogenic erosion rates, bedrock characteristics, and calibrated stream power incision models.We performed a series of channel profile analyses to quantify local channel slopes and their respective drainage area values, spatial distribution of knickpoints, and transient channel incision depths along the 15 channels actively responding to an increase in base level fall (Figures 1(a) and 1(b)). We measured bedrock fracture density and Schmidt hammer rebound values in order to assess variations in rock strength, which can influence the transient behavior of rivers [14, 60–63]. We use a numerical stream power model to constrain the rates of pre- and postbase level fall and river incision parameters for all 15 channels in our study area. Catchment-averaged erosion rates have been calculated to provide further constraints on these values.We extract the elevation, distance, and contributing drainage area from a 10 m DEM (USGS, 2018) for each channel by delineating a Topotoolbox v2 flow direction algorithm in MATLAB that follows a path of 10 m pixels along reaches draining >2 km2 [64, 65]. Elevation data is smoothed over a 5-pixel smoothing window (50-70 m depending on flow path) to remove inherent noise in the river profile [11, 15].To compare the slopes of channels that vary in size, we evaluate channel steepness by dividing channel slopes by the upstream drainage area with an assumed reference concavity (θref). This normalized value is referred to as ksn [15, 62]. We use a fixed reference concavity index that is the average of all 15 relict channel concavities (0.54, Table 1) which is approximate to the concavity values commonly used by other studies in the region [11] and consistent with other empirical studies for river profiles near equilibrium [10, 13, 15, 18, 70].In landscapes responding to an increase in base level fall, low normalized channel steepness (ksn) values at higher elevations are interpreted as the remainder of a relict profile equilibrated to the previous rate of rock uplift, while high ksn values represent the actively incising transient profile that is adjusted to the increased rate of base level fall (Figures 1(a) and 2). The relationship between the relict and transient ksn values have been shown to represent the relative change in the old and new rate of uplift, erosion, or base level fall [9, 14, 17, 18, 71–74].Equation (1) predicts that a temporal change in base level fall will initiate a steepening of the channel gradient at the channel outlet and will progressively migrate upstream until the channel gradient has reached equilibrium with the new rate of base level fall (Figure 2) [1]. The steep, headward propagating reach is known as the transient channel, while the upstream portion is deemed the relict (Figures 2 and 3(a)–3(c)). The knickpoint is the transition from the low eroding relict reach that has not been impacted by the new rate of base level fall to the steep, swiftly eroding transient profile that is adjusting to the new rate of base level fall (Figures 2 and 3(a)–3(c)). Knickpoints are identified as breaks in a river’s slope area relationship (Figures 3(a)–3(c)) [15]. These breaks can manifest as either a sharp discontinuity in channel slope or as the top of a convex-upward reach called a knickzone [11]. We define knickpoints as the lowest elevation of the relict reach where it diverges from a straight line in Χ-space (Figure 3(b)) using a MATLAB script developed by [11] containing functions from Topotoolbox v2 [64, 65].The upstream propagation of knickpoints following a new rate of base level fall can be measured in both vertical and horizontal travel distances relative to the channel outlet (Figures 3(a) and 3(b)) [9, 11, 14]. The vertical and horizontal knickpoint migration rates have different relationships with U, n, and K, where a knickpoint’s horizontal celerity is dependent on the new uplift rate (U) only when n>1 and vertical celerity are dependent on K only when n >1 [11].Transient channel incision depths reflect the total base level fall since fault initiation minus relict exhumation (Figure 2) [9–11, 76, 77] which is independent of catchment size [78]. We quantify the transient channel incision depths to provide constraints on the magnitude of base level fall along the outlet of each stream profile. We measure incision depths by projecting the relict reach downstream of the knickpoint (paleoriver) by using the relict reference concavity (0.54) and Ksn values in equation (4) (e.g., [77]) and taking the elevation difference from the paleo- to modern stream outlet, as shown in Figure 3(b).Through these equations, steepness indices, knickpoint distributions, and incision depths can be used to infer the spatial and temporal changes in the rate of base level fall within a given region, given that variations in climate and bedrock strength are minimal in relation to the new rates of base level fall [14].The concentration of 10Be in quartz rich fluvial sediment provides long-term (103−106 yrs) spatially averaged basin-wide erosion rates [18, 79–81] and is commonly used to provide constraints on the magnitude and timing of base level fall [10, 18]. 10Be in quartz is produced in proportion to the residence time of quartz grains in the uppermost few meters of the Earth’s surface; therefore, larger concentrations of 10Be suggest longer in situ residence time and a slower erosion rate and vice-versa [18]. We calculated relict and transient channel catchment-averaged erosion rates from cosmogenic 10Be in stream sediment from basins 2, 4, 5, 10, and 14, which span the length of the Western North Boulder fault. Samples were collected in the main-stem relict and transient channels that show no sign of upstream landslides or glaciation [18, 82–84], therefore ensuring that the sediment is representative of the entire watershed. Our samples were also collected in locations with variations in upstream lithology and therefore quartz concentration (e.g., the granite-dominated relict landscape contains a higher quartz concentration than the quartz-depleted downstream transient landscape of gabbro, diorite, and rhyolite), which will be discussed as a source of uncertainty in our erosion rate calculations. Quartz separation and purification were carried out according to the methods of [85]. The isolation and measurement of 10Be present within the quartz were carried out at the Purdue University PRIME AMS laboratory. Basin average cosmogenic nuclide production and erosion rates were calculated through the CAIRN method [86] and were corrected for topographic shielding.Bedrock strength and fracture density have an important role in the pace of the transient channel response. These properties cause variations in the bedrock erodibility parameter K in equation (1) and therefore modulate knickpoint migration rates [14, 60–62, 87]. To account for these effects, we collected nearly 5000 Schmidt hammer rebound values and over 150 bedrock fracture density values for each of the main lithological units along the 15 transient channels (Selby 1980). Schmidt hammer measurements were taken in groups of 60 consecutive readings distributed across bedrock surfaces that show minimal signs of weathering and fracturing [60]. Fracture density measurements were taken in the same proximal locations of Schmidt hammer measurements. Fracture density values were measured by counting the number of fractures or joints that intersect a meter-long transect of bedrock. Fracture density is therefore reported here as the number of fractures per meter of exposed bedrock. At each location, we used a horizontal and vertical transect to measure fracture density.We use a first-order upwind 1-D stream power-based river incision model to constrain the spatial and temporal patterns of fault slip along the Western North Boulder Fault. The goal of the model is to determine the rate and initiation time of base level fall required for the knickpoints to travel to their current locations. We test a wide range of uplift rates (or rates of base level fall), slope exponent (n) values, and bedrock erodibility (K) values in our simulations (Table 2). We use a weighted sum of square deviation (WSSD2) misfit analysis to calibrate acceptable ranges of these values.Assuming that the relict landscape is in steady state (i.e., relict uplift rates are equal to relict erosion rates), we use the CRN-derived relict and basin-wide erosion rates from basins 4, 5, 10, and 13 to constrain a range of 3 relict uplift (Uinitial) rates (2×10−5, 5×10−5, 8×10−5 m/yr; Table 2). In each model run, the relict reference concavity (0.54), Uinitial, K, and n values are applied to reach a steady-state profile (i.e., paleoriver profile) within all 15 channels. Starting at the equilibrium state, the downstream end node of the paleoriver profile is lowered to its modern position through time (Figures 4(a)–4(c)). We test a wide range of postacceleration uplift rates (Ut) that range from Uinitial to 1×10−3 m/yr, which encapsulates the estimated rates of base level fall for the Western North Boulder fault (>0.2 mm/yr, Johns et al., 1982) and other studies in the Northern Basin and Range province [33, 82]. For consistency, we focus on the results from the mean Uinitial (5×10−5 m/yr) simulations. We report the results with the minimum and maximum Uinitial simulations in the supplementary materials (available here).The exponent for channel slope (n) in equations (8) and (9) strongly influences transient channel behavior [1, 11] and is postulated to reflect dominant incision processes, with plucking consistent with n values of ~0.67 to 1 and abrasion consistent with n values of ~1.5 [67]. Because these incision processes are controlled by bedrock properties like tensile strength [3], compressive strength [60], and fracture density [61, 88], it is important that we recognize potential variations in n in our lithologically diverse study area. Therefore, we test all simulations with n values of 0.67, 1, and 1.5.We test a wide range of bedrock erodibility (K) values (1×10−7 to 1×10−4 m0.28/yr for n=0.67; 1×10−8 to 1×10−5 m-0.08/yr for n=1; 1×10−9 to 1×10−6 m-0.62/yr for n=1.5) that capture the range of calibrated stream power erodibility values found in other studies (i.e., [10, 19, 62, 89–91]). We assume that erodibility is constant along an individual river profile, which is reasonable for uniform lithology [10, 13]; yet, the lithological differences among different drainage basins in the study area (Figure 1(b)) will be expected to produce a larger amount of variability within the model results. A summary of the parameters used within each model run is displayed in Table 2. Each model run contains a specific n and Uinitial value and 100 intervals of base level fall (U) and bedrock erodibility (K) values, resulting in 10,000 profile simulations.We use the CRN-derived erosion rates for E and the mean relict and transient channel steepness (ksn) values in equation (10) to solve for K. We compare these K values to those constrained in the incision simulations to further assess the confidence in our results.Several relationships between the measured channel metrics and simulation results should exist (i.e., observed knickpoint horizontal travel distances in Χ-space scale with the modeled timing of base level fall) if the 15 transient channels in the study area behave in a manner consistent with the stream power model of a step function in base level fall (Table 1) [11, 17, 94]. For these relationships, we assume spatially and temporally constant K, n, and m. It is important to note that horizontal and vertical knickpoint migration rates have different relationships with U, 𝑛, and 𝐾 (Table 2) [11]. We report the relationships in Table 2 to assess if the behavior of the 15 channels in our study area is consistent with the stream power equation.The pattern of normalized channel steepness (ksn) mimics the distribution of hillslope gradients proximal to channels, where low channel steepness values are found in the flat, upstream relict landscape, while high channel slopes are found in the steep, narrow transient reaches of each channel (Figures 1(a) and 1(b)). Steepness increases by an average factor of 2.9 from relict (109±44 m1.08 average) to transient (286±96 m1.08 average) reaches (Table 3), which should roughly reflect the difference between the relict (Uinitial) and new (U) rates of base level fall depending on the slope exponent (n) value [15]. Most relict and transient channel segments appear to be linear on Χ plots (Figure 5(d)), confirming that the reference concavity (θ) of 0.54 is a reasonable estimate of the concavity of the relict and transient channels along Bull Mountain [17, 21, 75].All 15 transient channels contain a single knickpoint with a slope-break morphology (Figures 5(a)–5(d)). We observe that the knickpoints within most of the profiles contain a more stretched, elongated knickzone that is consistent with n values<1, rather than a sharp slope break knickpoint zone that is found with n values ≥1 (Figure 5(d)) [11, 21, 36].Knickpoint elevations and transient channel incision depths have a general pattern consistent with normal fault displacement—higher values towards the center of the Western North Boulder fault (Figures 5(a)–5(c)). Incision depth in basin 1 at the southern end is ~124 m, which increases to 497 m incision in basin 8 and then decreases to 99 m incision by basin 15 (Figures 5(a) and 5(b), Table 3). There is some variability from the general trend, including relatively high (basin 8, 497 m) and low (basin 6, 207 m; basin 7, 236 m; and basin 12, 59 m) incision depths (Figure 5(b)). These variations in knickpoint elevation could be attributed to localized variations in the temporal and spatial patterns of base level fall (i.e., due to normal fault linkage) [17], noise within the DEM [15], or deviations in relict steepness values that lead to inaccurate paleoriver reconstructions. Nonetheless, this elliptical pattern of transient channel incision depths suggests that the duration (i.e., initiation timing) and/or magnitude of base level fall is the highest near the center of the fault around basin 8 and minimizes towards the north end near basin 14, which is consistent with the general model of the slip distribution along a normal fault [17, 33].The knickpoint travel distances in Χ-space are similar (~2 km) from basins 1 to 9 and then decrease northward along strike from basins 9 to 15 to ~1 km (Figure 5(c)). There is some variability from this general trend, including relatively low (basin 2, 1237 m; basin 10, 893 m) and high (basin 15, 1175 m) horizontal travel distances in Χ-space (Figure 5(c)). These variations could be explained through local variations in the timing of base level fall initiation, channel narrowing effects [9, 14, 23, 95–98] or variations in the erosional resistance of bedrock [60, 62]. Because the horizontal travel distance of knickpoints formed at the same time should scale with the drainage area [28, 71, 99, 100] and therefore plot at the same Χ values [11, 17, 21, 94], these data suggests that the initiation of base level fall is similar in the south to middle portion of the fault and then relatively later in the northern end of the fault.Beryllium concentrations are consistent with the expectation of the general pattern of hillslope and river steepness. Concentrations in the relict landscape of basins 4, 5, and 10 range from ~28 to 78×104 atoms/g SiO2, resulting in erosion rates of 0.019±0.001 to 0.049±0.002 mm/yr (Table 4, Figure 5(a)). Concentrations taken at the mouth of basins 2, 4, 5, 10, and 14, which include both relict and transient topography (basin-wide), range from ~34 to 61×104 atoms/g SiO2, resulting in erosion rates of 0.022±0.002 to 0.042±0.001 mm/yr (Table 4, Figure 5(a)). Using a drainage area weighted mixing model [101] along with rates of 0.019±0.001 to 0.049±0.002 mm/yr for the relict, we back out an estimated range of 0.021±0.002 to 0.04±0.001 mm/yr for the transient/adjusted portion of the landscape in basins 4, 5, and 10 (Table 4, Figure 5(a)).Overall, no observable erosion rate trends exist from the northern to southern relict or adjusted basins. Relict, adjusted, and basin-wide erosion rates do not show a significant pattern of increasing with channel steepness (Figure 5(a)), and we observe no relationship between the erosion rate and lithology. For example, in basin 4, the relict erosion rate (0.049±0.002 mm/yr) is higher than the transient/adjusted (0.031±0.001 mm/yr), deviating from the expectation of the increasing erosion rate within increasing channel steepness. We discuss these discrepancies below.The Schmidt hammer values show little variation northward along strike, with median values having a narrow range from 57 to 65 (Figures 6(a) and 7). This suggests that all major lithological units have similar compressive strength ([60, 62]; Selby, 1980). Although incision into bedrock streams occur through breakage in tension, not compression, studies have shown that the compressive strength of rock is strongly correlated to the tensile strength of bedrock [60]. Therefore, we assume that the basin-wide tensile strength values also show no considerable variation northward along strike.The basin-wide fracture density values, however, display considerable variation, with median values ranging from ~0.8 to 22 fractures every meter (Figures 6(b) and 7. We observe far fewer fractures in the Boulder Batholith granite and similar intrusive units than in the extrusive ignimbrites and dacite porphyries (Figure 7).The stream power model does a reasonable job reproducing all transient channel profiles in our study area (e.g., basin 10 in Figure 4) given the constraints from the erosion rate and channel morphology data. The weighted sum of square deviation (WSSD2) misfit analysis results show that the percentage of acceptable fits for the simulations using a step function in the base level fall rate is significantly higher than the simulations with a linearly increasing rate of base level fall (Figure 8). This result might be expected given that the steepness of adjusted reaches in all 15 basins is near constant (Figure 5(d)) and therefore reflects more of a step function rate of base level fall [17], rather than temporally increasing. Therefore, we focus on the step function in base level fall simulations for the rest of the analysis but report the linear increase results in the supplemental.In the step function in base level fall simulations, the fault slip initiation times vary depending on the relict base level fall (Uinitial) rate. The simulations with the mean Uinitial rate (5×10−5 m/yr) resulted in a maximum initiation time of ~6 Ma for n=1 (Figure 9). Using the minimum Uinitial rate (2×10−5 m/yr), the model simulations predict the oldest initiation time (max of ~8.5 Ma for n=1), while simulations with the max Uinitial rate (8×10−5 m/yr) produced the youngest initiation time (max of ~4.5 Ma for n=1) (Supplemental Figures S5 and S6). This suggests that the lower the Uinitial rate, the earlier the estimated initiation of new base level fall. For consistency, we focus on the step function in base level fall simulations with the mean Uinitial rate (5×10−5 m/yr) for the rest of the analysis.Figure 9 illustrates the range of acceptable model fits for each basin and each value of n. For the simulations using a step function in base level fall, WSSD2 results show that the simulations with n values of 0.67, 1, and 1.5 generated similar relative patterns for the timing of the increase in along-strike base level fall initiations (Figure 9). The magnitude of the base level fall for each of the 3 n value simulations, however, differs significantly, where the smaller the n value results in earlier initiation of base level fall (Figure 9). This makes sense as the erosion rate, which controls the rate of river response to the new base level fall, is less sensitive to changes in slope with smaller values of n [11].All three simulations have a general along-strike trend of older initiations of base level fall in the southern to middle basins (1 to 8) along the Western North Boulder fault (median initiations from 8.8 Ma to 4.5 Ma for n=0.67; 6.2 to 2.5 Ma for n=1; 3.9 to 1.4 Ma for n=1.5). The results in the northern basins (9, 10, 13 to 15) suggest more recent fault initiation (median initiations from 2.6 Ma to 0.5 Ma for n=0.67; 1.5 Ma to 0.4 Ma for n=1; 0.7 Ma to 0.2 Ma for n=1.5). Basins 11 and 12 deviate from this trend with older fault initiation times similar to the southern/middle basins (median initiations from 4.6 to 4.2 Ma for n=0.67; ~2.5 Ma for n=1; 1.7 Ma to 1.5 Ma for n=1.5). The simulations with n=1 generated the highest percent of acceptable fits for all basins aside from 14 and 15, with the highest percent of fits being near 30% and the lowest percent of fits around 0.3% (Figure 9). However, even though the specified n value changes the geometry of the knickzone and transient profile (e.g., [11]), these geometric changes were not large enough to produce significant variations in the percent of acceptable fits, given the resolution of the DEM (Figures 4(a)–4(c) and 9).We assess several morphologic stream power relationships between the measured channel metrics and simulation results (Table 2) to see if the 15 transient channels behave in a manner consistent with the stream power model of a step function in base level fall (Figures 10(a)–10(c)). Knickpoint horizontal travel distances in Χ-space have a weak pattern (R2=0.34) of increasing with the median simulated initiation of base level fall (Figure 10(c)). Basins 7 (2 km travel distance), 8 (2.4 km travel distance), and 9 (2.1 km travel distance) have the largest deviation from the general trend and contain some of the highest travel distances in the study area.The average knickpoint horizontal migration rates in Χ-space have no significant correlation (R2=0.04) with their respective basin-wide median bedrock fracture density values (Figure 10(b)). The fracture density values also do not have any significant relationship with the median modeled bedrock erodibility values (Figure 10(b)). The mean transient channel normalized steepness (ksn) values have a weak pattern (R2=0.32) of increasing with the median modeled uplift rates (or rates of base level fall) (Figure 10(c)). Again, basins 7, 8, and 9 have the largest deviation from the general trend and contain some of the highest estimated median uplift rates (3.8×10−4 m/yr, 2.5×10−4 m/yr, and 2×10−4 m/yr, respectively) and fracture densities in the study area (Figures 6(c), 7, and 10(b). The fact that these three basins have the highest fracture densities in the region suggests that complex incisional processes could be occurring along their channels (i.e. ,spatially variant K and n values) causing them to deviate from stream power relationships assuming spatially uniform K, n, and m. Other factors such as drainage divide instability (potentially suggested by the differences in shape and size of these three basins as well as surrounding basins) could also explain the source of these outliers [102].The constrained bedrock erodibility (K) values for all 15 basins range from ~6.8×10−7 to 5.6×10−6 m0.28/yr for the incision simulations with an n value of 0.67, ~2.2×10−7 to 2.4×10−6 m-0.08/yr with an n value of 1, and ~1.1×10−8 to 3.2×10−7 m-0.62/yr with an n value of 1.5 (Figure 11). These K values, particularly those from the n=1 simulations, do not show any significant spatial trends (Figure 10(a)–10(c)). These ranges of constrained K values are consistent with the K values calculated from the CRN-derived erosion rates in equation (10) and calibrated stream power K values reported in other studies (Figure 11) [10, 19, 62, 63, 67, 89–91, 103, 104]. The observed overlap between the two methods of constraining K (i.e., from incision simulations and CRN-derived erosion rates) as well as other published studies increases our confidence in these results.Rivers are powerful tools for understanding the spatial and temporal patterns of normal fault evolution, especially where geologic and geodetic data is limited. Our results are consistent with the widely observed patterns of normal fault slip distribution/evolution (i.e., higher slip displacement along the central portion of the fault relative to the distal portions). Furthermore, our cosmogenic erosion rates and base level fall simulations provide the first quantitative spatiotemporal constraints of fault slip initiation and evolution along the Western North Boulder fault, possibly the farthest North Basin in the northern Basin and Range province that such constraints exist.Overall, our analysis of a wide range of models using a step function in base level fall suggests that fault slip accelerated along the central portion of the Western North Boulder fault sometime between ~8.8 and 3.9 Ma. The slip initiated near basin 7 and propagated to the north and south, followed by the initiation of slip towards the south near basins 1-8 between ~6.5 and 1.4 Ma and further slip initiation near basins 9-15 in the northern end of Bull Mountain between ~2.6 and 0.2 Ma (Figure 9). The timing of base level fall strongly depends on each channel’s slope exponent (n) (Figure 9), as variations in n cause the wide range of slip initiation values we report here.Here, we discuss which range of fault initiation values are we most confident with. Combined studies have suggested that the ongoing period of increased extension within the Northern Basin and Range province and along the Western North Boulder fault initiated sometime around 6 to 0.75 Ma and is continuing into the present ([51–53]; Johns et al., 1982; [54, 55]). The modeled timing of slip initiation for the n=1 and n=1.5 simulations are from ~6.2 to 2.5 and ~3.9 to 1.4 Ma, respectively, which all fit within our existing constraints. However, the n=1 simulations produced a significantly higher number of acceptable fits than the n=1.5 simulations. This makes sense given that the n=1.5 simulations produced sharp slope break knickpoints (rather than more broad, stretched knickzones seen with the n=0.67 and n=1 simulations) that are not represented in any of the channels. Therefore, we place the most confidence in the n=1 simulations with fault initiation times of ~6.2 to 2.5 Ma. This time range coincides well with the migrations of Yellowstone hotspot over the Idaho-Montana border from ~6-2 Ma [51–55]. Therefore, we attribute the increased extension along the Western North Boulder fault to the crustal bulge from the migrating hotspot.We built an interpretive reconstruction of the Western North Boulder fault evolution using (1) incision depth values, (2) normal fault length-displacement scaling relationships, and (3) median values of base level fall rate and initiation from the step function in base level fall simulations with n=1. Previous work suggests that large normal fault zones like the Western North Boulder fault can form by the nucleation and propagation of isolated fault segments, followed by an intermediate stage of overlap and linkage of these segments (i.e., soft-linked faults) that eventually evolve into larger hard-linked faults that kinematically behave like isolated faults [5, 31, 35]. We adhere to these 3 stages of normal fault evolution in our interpretive reconstruction, starting with individual fault segments nucleating at the outlet of basins followed by the propagation and linkage of these segments into a large hard linked fault (i.e., the Western North Boulder fault).Studies suggest that a rapid increase in fault slip rate occurs following the linkage of two isolated faults such that it reestablishes the displacement-length ratio ([9, 16, 17, 30, 32, 33, 105, 106], 2016; [35, 78]). This action of relatively rapid slip in zones of fault linkage is applied in our interpretive reconstruction (Figure 12).While the two end members of fault evolution (isolated and hard-linked faults) maintain similar displacement-length ratios for the entire fault system [32, 35, 106], studies show that the dynamics of fault linkage and displacement readjustment can develop a wide range of displacement-length ratios (Dmax~0.001-0.1 L) [30, 107–110]. We fit displacement-length contours by hand, assuming the ratio of Dmax~0.03 L throughout the beginning (i.e., isolated faults) stages of fault evolution; yet, during the intermediate (i.e., fault linkage) and final (i.e., hard-linked faults) stages, we find it difficult to maintain this scaling and therefore relax the ratio to a range of Dmax~0.01-0.03 L in order to maintain fault scaling relationships. Throughout all phases of fault evolution, and particularly in the zones of fault linkage, we fit the displacement-length contours to best match the constrained present distribution of total fault displacement values (i.e., black circles in Figure 12) while obeying the predicted range of scaling ratios. An illustrative guide on how we created our interpretive reconstruction is shown in the supplemental (Figure S8).The largest incision depth is 497 m in basin 8 (Figure 5(b)), which is at least 140 m larger than its neighboring basins 7 (236 m) and 9 (355 m). The relict channel reconstruction for basin 8 has significant uncertainty due to its relatively small relict channel. Based on fault scaling relationships, we expect that the incision depth of basin 8 is overestimated. In order to obey the displacement-length ratio (equation (12)), we thus lower the incision depth of basin 8 to 400 m (462 m fault displacement), which still serves as the zone of maximum displacement and provides a more consistent displacement relationship pattern with neighboring basins (Figure 12).Our interpretive fault evolution reconstruction begins with the nucleation of isolated fault segments that increase in displacement, propagate at the fault tips, and link together through time, developing relay ramps in the zones of linkage that accumulate extra displacement in order to maintain the range of displacement-length ratios in the resulting hard-linked faults (Dmax~0.01-0.03 L) (Figure 12). Studies suggest that two sets of knickpoints in each channel would be expected if fault initiation was being recorded, in addition to a later linkage event [9, 78]. We do not observe multiple knickpoints in any of the channels; therefore, either (1) no linkage took place along the outlet of the channels, (2) the knickpoints formed by linkage consumed the older knickpoints, which can occur if n>1 [11], or (3) the fault initiation took place in a zone of linkage (i.e., basin 14, Figure 12).The fault at basin 1 propagates only a few kilometers southward, which serves as a plausible reason why no transience is observed in rivers south of basin 1. Other plausible reasons include the fact that different rock types such as carbonates and sandstones also exist in rivers south of basin 1, potentially modulating transient channel response. However, given the consistency in fault length and our observed displacement distribution, we prefer the model in which the fault is pinched out just south of basin 1. Overall, given the listed assumptions above, the spatial and temporal slip patterns of the Western North Boulder fault were recreated assuming a step function in base level fall and an n value of 1.With an assumed normal fault geometry of 60°, we calculate a range of slip rates for each basin with initiation times from the step function in base level fall and n=1 simulations (Table 5). Slip rates across all 15 basins range from 0.02 to 0.45 mm/yr (Table 5). The highest slip rates take place along basin 14 (0.27-0.45 mm/yr), the outlet of which is assumed to be in a zone of fault linkage (Figure 12), and therefore, this is a reasonable result. Relatively high slip rates also take place along basins 8 (0.17-0.22 mm/yr) and 9 (0.32-0.41 mm/yr), which are in the central portions of the fault, a region that is anticipated to contain relatively high slip rates [9, 30]. Furthermore, all of the calculated slip rates are comparable to the USGS estimated slip rates of the Western North Boulder Fault (<0.2 mm/yr, Johns et al., 1982) and other normal faults in the region, including the Beaver Creek Fault (<0.2 mm/yr), Tobacco Root Fault (<0.2 mm/yr), Bridger Fault (<0.2 mm/yr) (U.S. Geological Survey and Montana Bureau of Mines and Geology, Quaternary fault and fold database for the United States), and the Lemhi, Lost River, and Beaverhead faults (≤0.3 mm/yr) [33], thus increasing our confidence in the results.Our range of transient/adjusted landscape erosion rates (0.021±0.002 to 0.04±0.001 mm/yr) falls within the low end of the calculated slip rates (0.02 to 0.45 mm/yr). We conclude that these erosion rates may be underestimates due to the limited extent of the transient/adjusted landscape and variable lithology within the basins (e.g., presence of nonquartz yielding units such as gabbro within the transient/adjusted landscape).Our results show that the percent of acceptable fits is significantly higher for the step function in base level fall rather than the linearly increasing rate of base level fall (Figure 8). We expect the large percent of acceptable fits for this step function in base level fall for the following reasons: (1) all channels contain a single knickpoint with a slope-break morphology (Figures 5(a) and 5(d)) that displays a consistent change in channel steepness above and below the knickpoint (Figure 1(a)), which typically develops in response to a step function in tectonic forcing [70]. (2) The consistency of adjusted channel ksn values and their linearity on Χ-plots (Figures 5(a) and 5(d)) suggests that either the rate of base level fall is near constant (i.e., step function) from its initiation to the present [11] or it is hidden by our reference concavity and assumption of spatially uniform uplift. (3) All 15 knickpoints do not correlate with a specific rock type (Figures 5(d) and 7), suggesting that knickpoint formation is not attributed to lithology. Therefore, we attribute the most probable formation of the slope-break knickpoints and pattern of ksn values in the transient and relict reaches to a two-phase base level fall history along Western North Boulder Fault, where the relict reaches are representative of the relatively low, old rate of base level fall, while the transient reaches are adjusting to the step function of a relatively faster rate of base level fall (Figures 1(a) and 1(b)). We acknowledge that there could be temporal variations in base level fall; however, any such variability is likely less than the difference between the modern and relict base level fall and is not resolvable with our river profile reconstruction approach.The relationships among knickpoint horizontal travel distance in Χ-space, initiation of base level fall, and bedrock properties (Figures 10(a)–10(c); Table 1) warrant some discussion. Specifically, knickpoint horizontal travel distance in Χ-space and initiation of base level fall have a weak correlation (Figure 10(a)); however, there is some variation from this trend (i.e., basins 7, 8, and 9). Deviation from this correlation is reasonable, as studies have shown that even after knickpoints are normalized for the drainage area, there is a trend for higher knickpoint retreat rates with higher slip rates, which may be driven by dynamic channel narrowing effects [9, 14, 23, 95–97, 111].In the simulations using a step function in base level fall and an n value of 1, our constrained bedrock erodibility (K) values range between ~2.2×10−7 to 2.4×10−6 m-0.08/yr (Figure 11). These values have a strong pattern of increasing with the normalized knickpoint horizontal travel distances (Figure 10(b)). This intuitively makes sense, as the higher the bedrock erodibility (K), the faster the knickpoint can migrate upstream [11]. If our assumptions are correct (i.e., step function change in base level fall and no dynamic channel width adjustments), then we have confidence that the modeled bedrock erodibility (K) values reflect the actual values in nature. This claim is reinforced by the fact that our constrained K values fall within range of the K values calculated from our CRN-derived erosion rates and the calibrated stream power K values reported in other studies (Figure 11).Because the slope exponent n values are postulated to reflect dominant incision processes, with plucking consistent with n values of ~0.67 to 1 and abrasion consistent with n values of ~1.5 [67], the poor correlation between the knickpoint travel distances (normalized by respective fault initiation times for n =1 simulations) and respective fracture density values (Figure 10(b)) suggests that either (1) abrasion (correlated with Schmidt rebound values) is a more dominant process than plucking (correlated with fracture density), (2) the fracture density values are not truly representative of the basins, or (3) the Schmidt rebound and fracture density data collected are not sufficient enough to deduce the prominent erosional processes taking place within the channels (Figures 10(a) and 10(b)). The latter of these hypotheses is also supported by the poor correlation between mean transient channel ksn values and median uplift rates (Figure 10(c)). Even with the confidence in our modeled K values, we clearly do not have a strong metric for erodibility (e.g., fracture spacing). We acknowledge that bedrock strength, fracture density, and other factors such as grain size, sediment load, bed cover, drainage area, n value, drainage capture, and drainage divide instabilities vary in space and time during transient incision, further complicating dominant erosional processes and the rate of knickpoint retreat in our study area [23, 36, 60, 62, 95, 111–114].There are several assumptions that are necessary to obey the fault scaling relationships in our Western North Boulder fault reconstruction. First, basin 7 has a significantly older initiation of base level fall (>2 Ma) than its neighboring basins 6 and 8 (Figure 9). To obey fault scaling relationships, we recognize that the initiation of base level fall is unreasonably older than surrounding basins due to its low modeled bedrock erodibility value, causing a slower rate of knickpoint migration and therefore an older initiation of base level fall. This low bedrock erodibility value is a result of its relatively high normalized steepness value, which may be due to drainage capture from surrounding basins. This is suggested by the differences in shape and size of neighboring basins (Figures 1(a) and 1(b) and 5(a)) [58]. If we assume that its fault initiation is similar to neighboring basins 6 and 8, we reconstruct more realistic fault scaling relationships. Secondly, basins 2, 3, and 4 all share the same outlet, yet vary in the initiation of base level fall due to differences in knickpoint horizontal travel distances, elevations, and transient channel incision depths (Figures 5(a)–5(d)). These differences in knickpoint metrics could be due to drainage capture or lithological variations (Figure 7). While we acknowledge variability exists from basin to basin, the general trend of fault slip and initiation follows expectations from the fault slip distribution and timing in the Northern Basin and Range province and Southwest Montana.Furthermore, we acknowledge that the cosmogenic erosion rate data only integrates over ~10 k.y. but the landscape has been evolving over millions of years. If the cosmogenic erosion rate data underestimates the true long-term erosion rates, bedrock erodibility values would be larger and therefore, the modeled initiations of base level fall would be younger. Conversely, an overestimate of long-term erosion rates would cause the modeled initiation of the base level to be older. Nonetheless, the general pattern of initiation times would be the same.The erosion rate values also contain considerable error due to the limited extent of the transient/adjusted landscape and varying lithologies and therefore quartz concentrations within the basins. This is most clearly seen in basin 4, where the granite-dominated relict landscape contains a higher quartz concentration and erosion rate than the quartz-depleted downstream transient landscape of gabbro and diorite. Here, it is likely that most of the quartz within the transient landscape were shed from the relict landscape, and therefore, the corresponding erosion rates are not truly reflective of their surrounding topography. This source of error applies to all basins with variable lithology and limited extent of the transient/adjusted landscape.While absolute values of the fault reconstruction in Figure 12 have uncertainty, we are confident that the general pattern is robust. Moreover, we provide the first quantitative constraints on fault initiation and evolution of the Western North Boulder fault; perhaps, the farthest north basin in the Northern Basin and Range province that such constraints exist. We illustrate the utility of river profile analysis in documenting normal fault evolution, especially in young systems where slip rates and total throw are too low to extract via other geologic methods.The eastern margin of Bull Mountain in the Northern Basin and Range province contains a steep, dissected transient landscape that is most likely attributed to the recent base level fall of the Western North Boulder fault. We show that transient channel morphometrics, rock strength metrics, cosmogenic erosion rates, and calibrated stream power incision models can constrain the temporal and spatial patterns of base level fall along the Western North Boulder fault. We show that with a step function in base level fall and an n value of 1, the Western North Boulder fault started as individual fault segments along the middle to southern regions of Bull Mountain that nucleated around 6.2 to 2.5 Ma, respectively. This was followed by the nucleation of other fault segments in the northern region around 1.5 to 0.4 Ma. Through time, these faults linked together to span over 40 km, with a maximum fault slip of 462 m (basin 8) in the central portion of the fault. Fault slip rates range from 0.02 to 0.45 mm/yr along strike and are similar to estimates for other faults in the region.Our results augment existing constraints of Quaternary fault slip histories in the Northern Basin and Range province. We find that the timing of fault initiation coincides well with the migration of the Yellowstone hotspot across the Idaho-Montana border. Therefore, we confidently attribute the increased extension along the Western North Boulder fault to the crustal bulge from the migrating hotspot. However, the base level fall simulations show that the slope exponent (n) value can significantly change the estimated timing of fault initiation required to match observed knickpoint locations, highlighting the importance of this parameter when using the stream power equation. Furthermore, our range of calibrated bedrock erodibility (K) values are comparable to those calculated from the cosmogenic erosion rates and calibrated stream power K values reported in other studies, increasing confidence in our results. Overall, we show that rivers are powerful tools for documenting the spatial and temporal patterns of normal fault evolution, especially where other geologic/geodetic methods are limited, proving to be a vital tool for accurate tectonic hazard assessments.This study’s data are available in the IUScholarWorks digital repository at Indiana University (http://hdl.handle.net/2022/25807).The authors declare no conflict of interest regarding this publication.This work has been conducted with financial support from the National Science Foundation grants 4827341 and 1727139, the Army Research Laboratory grant W911NF-17-1-0248, the Geological Society of America Robert K. Fahnestock grant, the Indiana University Field Station Vitaliano grant, Indiana University Graduate and Professional Student Government, and the Indiana University Department of Earth and Atmospheric sciences. Yanites and Mitchell was supported by the National Science Foundation grant EAR-1727139. A large portion of data in this study was collected on the properties of the Carey, Dawson, and McCauley families, Ms. Norby, Ms. Warthen, and Dr. Pontiac. Cosmogenic nuclide concentrations were measured at the Purdue University PRIME laboratory. This work greatly benefited from thoughtful reviews by Sean Gallen, Sarah Boulton, and an anonymous reviewer; however, any errors or inaccuracies presented here are solely the responsibility of the authors.The attached supplementary material file contains 14 figures displaying calibrated stream power incision model results using different uplift rates through time, a guide to our interpretive reconstruction of the Western North Boulder fault evolution, and field photos showing some examples of the relict and transient reaches in the studied basins. Figure descriptions are located below each respective figure.
更新日期：2021-01-08
京公网安备 11010802027423号