The uplift history of the Colorado Plateau has been debated for over a century with still no unified hypotheses for the cause, timing, and rate of uplift. 40Ar/39Ar and K/Ar dating of recurrent basaltic volcanism over the past ∼6 Ma within the Virgin River drainage system, southwest Utah, northwest Arizona, and southern Nevada, provides a way to reconstruct paleoprofiles and quantify differential river incision across the boundary faults of the Colorado Plateau–Basin and Range boundary. We compare differential incision data with patterns of channel steepness, bedrock erodibility, basaltic migration, and mantle velocity structure to understand the birth and evolution of the Virgin River system.

New detrital sanidine ages constrain the arrival of the Virgin River across the Virgin Mountains to less than 5.9 Ma. Virgin River incision rates and amounts show an eastward stair-step increase in bedrock incision across multiple N-S–trending normal faults. Using block incision values away from fault-related flexures, average bedrock incision rates are near zero since 4.6 Ma in the Lower Colorado River corridor, 23 m/Ma from 6.8 to 3.6 Ma in the Lake Mead block, 85 m/Ma from 3 to 0.4 Ma in the combined St. George and Hurricane blocks, and 338 m/Ma from 1 to 0.1 Ma in the Zion block. Steady incision within each block is documented by incision constraints that span these age ranges. We test two end-member hypotheses to explain the observed differential incision magnitudes and rates along the Virgin River system over the past ∼5 Ma: (1) as a measure of mantle-driven differential uplift of the Colorado Plateau relative to sea level; or (2) due to river integration across previously uplifted topography and differential rock types with down-dropping of Transition Zone blocks but no post–5 Ma uplift.

We favor headwater uplift of the Colorado Plateau because basalt-preserved paleoprofiles indicate that eastern fault blocks have been the “active” blocks that moved upwards relative to western blocks with little base-level change of the lower Colorado River corridor in the past 4.6 Ma. Block-to-block differential incision adds cumulatively such that the Zion block (Colorado Plateau edge) has been deeply incised 880–1200 m (∼338 m/Ma) over the 2.6–3.6 Ma period of Hurricane fault neotectonic movement, which has a slip magnitude of 1100 m. Mantle-driven uplift is implicated by a strong correlation throughout the Virgin River drainage between high normalized channel steepness (ksn) and low underlying mantle velocity, whereas there is a weaker correlation between high ksn and resistant lithologies. Basaltic volcanism has migrated northeastward at a rate of ∼18 km/Ma parallel to the Virgin River between ca. 13 and 0.5 Ma, also suggesting a mantle-driven mechanism for the combined epeirogenic uplift of the western Colorado Plateau, recurrent slip on its bounding faults, and headward propagation and differential incision of the Virgin River. Thus, we interpret the Virgin River to be a <5 Ma disequilibrium river system responding to ongoing upper-mantle modification and related basalt extraction that has driven ∼1 km of young (and ongoing) surface uplift of the western Colorado Plateau.

The timing and processes of Colorado Plateau (CP) uplift from sea level at 70 Ma to its present 2 km average elevation have been debated for over a century. CP uplift likely occurred in three stages: Laramide, mid-Tertiary, and post–10 Ma (Karlstrom et al., 2012), but the relative magnitude of each uplift episode is debated. Neogene uplift (post–10 Ma) of the western CP has been suggested by differential incision studies of the Colorado River (CR) through Grand Canyon (Karlstrom et al., 2007; Karlstrom et al., 2008; Crow et al., 2014). However, others explain the observed differential incision by means of geomorphic controls such as variable rock strengths (Pederson and Tressler, 2012; Bursztyn et al., 2015; Darling and Whipple, 2015) and/or cyclic climate changes (Small and Anderson, 1998; Molnar, 2004; Chapin, 2008), rejecting the young plateau uplift hypothesis.

We use the Virgin River (VR) fluvial system (Fig. 1) to further evaluate the hypothesis that the CP has been uplifted relative to sea level in the past 5 Ma. The East and North forks of the VR form the headwaters of the river system at an elevation of ∼2.5 km. The VR and its tributaries carve deep bedrock canyons (Zion, Parunuweap, and Virgin River Gorge) and cross several major normal fault systems (Piedmont, Grand Wash, Washington, and Hurricane; Fig. 1) before ultimately joining the CR in what is now the Overton Arm of Lake Mead, at a pre-dam elevation of 225 m (Birdseye, 1924). The common base level at the confluence of the CR and the VR since ca. 5 Ma that we document here allows comparisons between Grand Canyon and Zion Plateau incision along two different major rivers that have carved impressive young canyons across the CP–Basin and Range (BR) boundary zone.

We test the uplift versus no-uplift hypotheses by collecting and evaluating several data sets. (1) Nested river profiles of the VR and its major tributaries and knickpoints provide information about relationships between topography, faults, rock type, and other features that may help explain profile geometry (Fig. 2; Table 1). (2) The near continuum of 40Ar/39Ar and K/Ar dated basaltic volcanism from ca. 15 Ma to present (Fig. 3) provides timing constraints to evaluate differential river incision. Following Hamblin (1984), we use basalts that overlie river gravels to reconstruct paleoriver profiles of different ages that can be compared to the modern profile. A spatial trend of basaltic magmatism is also used to indicate relationships between mantle melting and surface responses (Supplemental Table S11). (3) A summary of documented displacement along major normal faults provides insight into the evolution of the VR landscape (Fig. 1; Table 2). (4) Approximately 100 new calculated incision rates are presented in the context of a synthesis of published incision rates (Hamblin et al., 1981; Willis and Biek, 2001; Hayden and Sable, 2008; Crow, 2012) to provide a complete differential incisional history of the VR and its tributaries (Table S2 [footnote 1]). (5) Correlations among normalized river steepness (ksn; Wobus et al., 2006; Kirby and Whipple, 2012), incision rate and magnitude, upper-mantle tomography (Schmandt and Humphreys, 2010), rock type, and precipitation are used to attempt to discriminate possible deep (mantle-driven uplift) and shallow (rock type and climate) controls on river evolution. (6) A sediment provenance analysis using detrital zircons and detrital sanidines of ancestral VR gravels is used to help date the birth of the VR system and examine its provenance (Table S3 [footnote 1]).

This study has three primary goals that contribute to the VR story. The first goal of this paper is to constrain the timing and processes of the VR birth. The birth of the through-going VR (integration of a major river draining from the CP to near its present base level) has been interpreted to be ca. 6–4 Ma based on studies of the VR depression filled with the fluvial upper Muddy Creek Formation (Williams, 1996; Dickinson et al., 2014). These interpreted ages are tested in this paper using detrital sanidine dating to search for young grains representing a maximum depositional age of the first-arriving, far-traveled, VR-derived gravels of the upper Muddy Creek Formation.

A second goal of this paper is to reconstruct the incision history of the VR drainage system. Differential incision studies are aided by the nearly continuous record of late Cenozoic basaltic volcanism (ranging from 14.19 to 0.12 Ma; Fig. 3) that allows for accurate dating of paleoriver deposits and ancient landscapes. This study builds on and extends similar studies by Hamblin et al. (1981), Hamblin (1984), Willis and Biek (2001), and Crow (2012). We also build on recent work supporting the ca. 6 Ma Grand Canyon hypothesis as proposed by Karlstrom et al. (2014), Darling and Whipple (2015), and Winn et al. (2017).

A third goal is to understand temporal and spatial relationships between incision and the locus of basaltic volcanism to help evaluate mechanisms of potential uplift. Some papers have hypothesized young and ongoing mantle-driven CP uplift, but proposed mechanisms have differed from: edge-driven, upper-mantle convection driving uplift (van Wijk et al., 2010), asthenospheric upwelling and buoyancy change driving uplift (Karlstrom et al., 2008; Crow et al., 2014), delamination of lithosphere and asthenospheric return flow around the Escalante anomaly (Levander et al., 2011), and whole mantle flow driving dynamic topography (Moucha et al., 2009). A long-term goal of our research is to evaluate when the CP’s ∼2 km modern elevation and its high-relief river gorges formed.

Compressional forces of the Sevier and Laramide orogenies followed by BR extension established the present physiographic framework of the VR region by ca. 17 Ma (Biek et al., 2010). This was the time of formation of the Grand Wash fault and establishment of significant topographic relief across the CP-BR boundary (Faulds et al., 2001). Best et al. (1980), Wenrich et al. (1995), and Nelson and Tingey (1997) observed an eastward-propagating sweep of basaltic volcanism in which basalts get younger and also become more asthenospheric in Nd composition to the east (Crow et al., 2011). The combined data are interpreted to suggest that the CP lithosphere is being thinned and replaced by asthenosphere as North America moves SW (absolute velocity) over warm mantle. This has been envisioned as the East Pacific Rise mantle domain by Moucha et al. (2009), edge-driven convection around the CP margin (Karlstrom et al., 2008; van Wijk et al., 2010; Crow et al., 2011), or lithospheric delamination tectonism and its surface response (Levander et al., 2011).

Similar to magmatism, normal faulting propagated west to east into the CP (Table 2; Pearthree, 1998; Karlstrom et al., 2007). Similarly, the Wasatch fault system that forms the western CP boundary to the north transitions southward into several faults that represent a southward continuation of the intermountain seismic zone (Smith and Arabasz, 1991), termed the Utah Transition Zone (Wannamaker et al., 2001). The Sevier-Toroweap and Hurricane faults act as the western neotectonic boundary of the CP with older (ca. 17 Ma; Faulds et al., 2001) normal faults of the Grand Wash fault zone forming the physiographic western boundary of the CP (Brumbaugh, 1987). The ∼90 km distance between the Grand Wash and Sevier-Toroweap faults represents a transition zone between the BR and CP provinces (Fig. 1).

Hamblin et al. (1981) first proposed calculating differential incision across the Hurricane fault as a measure of CP tectonic uplift; he calculated a differential uplift of 364 m/Ma across the Grand Wash and Hurricane faults using four incision points. Hamblin (1984) used these data and basalt-preserved paleoprofiles across faults to infer that the eastern (footwall) blocks were the “active” blocks. Willis and Biek (2001) added new 40Ar/39Ar ages of basalts and 12 new incision points along the VR and found 307 m/Ma of differential uplift across the Hurricane and Washington faults. Similar analyses in the Grand Canyon (Pederson et al., 2002; Karlstrom et al., 2007, 2008; Crow et al., 2014) also attempt to quantify differential neotectonic uplift. Howard et al. (2015) made the case that the CR was graded to near sea level by 4.5 Ma, with major aggradation of Bullhead gravels from 4.5 to 3.5 Ma (Howard et al., 2015; Crow et al., 2016). This reinforces the concept that differential incision along the CR system and major tributaries such as the VR can be used to infer uplift of the CP (Lucchitta, 1979; Karlstrom et al., 2007; Karlstrom et al., 2008) instead of downstream subsidence (Pederson et al., 2002; Ott et al., 2018).

Muddy Creek Formation

The Muddy Creek Formation has been separated into upper and lower units (Swenberg, 2012). The lower Muddy Creek Formation, located west of the Grand Wash fault, generally consists of locally derived, fine-grained, basin-fill deposits ranging in age from 11 to 4 Ma (Pederson, 2008; Dickinson et al., 2014). The upper Muddy Creek Formation (Muddy Creek conglomerate of Williams, 1996) contains coarse-grained, well-rounded, gravel clasts that were deposited ca. 4 Ma. Williams (1996) mapped a sharp contact of well-rounded gravels and cobbles above fine-grained sands of the Muddy Creek Formation ∼1 mile northwest of Mesquite, Nevada, near the mouth of the Virgin River Gorge. He interpreted this to represent the rapid arrival of the ancestral VR across the Virgin Mountains, which he considered to have been a previous drainage divide. Swenberg (2012) found similar outcrops of well-rounded, far-traveled gravels above angular, locally derived gravels near the northernmost extent of Lake Mead.

The birth of the VR drainage system is defined here as the age of arrival of the first coarse-grained, far-traveled gravels to appear in the Mesquite Basin (eastern basin adjacent to the Piedmont fault within the Virgin River Depression; Fig. 1). Pederson (2008), in agreement with Longwell (1928), performed petrographic analyses of fine-grained sand of lower Muddy Creek Formation and concluded that none of the sediments represented material derived from the CP and hence were not paleo-CR deposits but, instead, were locally derived. Provenance analyses on dated zircons of the Muddy Creek Formation provide evidence for the arrival of an ancestral river draining the CP across the Virgin Mountains at 6–4 Ma (Muntean, 2012; Dickinson et al., 2014).

Profile Analysis

Longitudinal profiles of the VR and ten of its major tributaries were extracted from 10 m digital elevation models (DEMs) using ArcGIS. Distances of river profiles were calculated starting at the historic (pre-dam) confluence of the VR and CR as mapped by Birdseye (1924). Observation of knickpoints and varying stream gradients were quantified using Topotoolbox in MATLAB to calculate normalized channel steepness values (ksn) from a U.S. Geological Survey (USGS) 30 m DEM of the entire drainage basin (Schwanghart and Kuhn, 2010; Schwanghart and Scherler, 2014).

Values for ksn give the gradient of the streambed normalized by the upstream area at each location along the channel:
where S is the gradient, A is the upstream drainage basin area, and Θref is the reference concavity index. This method uses upstream contributing area as a proxy for discharge so that shallower gradients of downstream channels with high discharges can be compared to steeper gradients of headwater channels with lower discharges.

To perform correlations among different data sets (ksn, gradient, mantle velocity anomalies, rock type, and annual precipitation), data were collected at 1 km intervals along the entire drainage system. This study performs an analysis of the correlation between upper-mantle velocity anomalies and ksn values along the VR and its tributaries using 30 m DEM instead of previously used 90 m DEM (Crow, 2012). We use mantle tomography from Schmandt and Humphreys (2010) with velocity values below different geographic locations extracted from the 2010 3D tomographic model ( In this data set, mantle velocity is measured as a percent change from a regional average mantle velocity. We refer to relative mantle velocity anomaly data (%) simply as “mantle velocity.”

Rock type was also taken into consideration by classifying the rock type underlying each data point into slope- and cliff-forming lithologies. The slope and cliff formers were identified based upon the unit descriptions found in geologic maps throughout the region (Billingsley and Workman, 2000; Billingsley and Wellmeyer, 2003; Sable and Hereford, 2004; Beard et al., 2007; Ludington et al., 2007; Biek et al., 2010). This study also compares tensile strength measurements taken from the CR drainage (Bursztyn et al., 2015) and channel steepness within the VR drainage (this study). Bursztyn et al. (2015) provided average tensile strengths for seven formations (Navajo, Kayenta, Shinarump, Moenkopi, Kaibab, Esplanade, and Redwall) found in common between the two drainages. Formation average tensile strengths measured along the CR were then compared to calculated formation average ksn data along the VR. This analysis assumes tensile rock strength of a particular formation is consistent across the region in adjacent drainages. The underlying mantle velocities at the same data point locations that overlie the seven formations were also used to search for a correlation between mantle velocity and channel steepness.

We use modern mean annual precipitation around the VR region (National Atlas of the United States, 2010) as a proxy for spatial variations in climate to compare to channel steepness data sets. Climate has changed through time and space; however, this is the best available data set that can represent spatial climatic variations at the scale of the VR drainage.

Magmatic Sweep

Previously published 40Ar/39Ar and K/Ar basalt ages throughout the eastern VR drainage system were compiled into Table S1 (footnote 1). K/Ar and 40Ar/39Ar ages were recalculated using a modern Fish Canyon Tuff standard age of 28.201 Ma (Kuiper et al., 2008) and a modern decay constant of 5.463 × 10−10/a (Min et al., 2000). Some published K/Ar ages contained insufficient amounts of information to perform the calculations and therefore were not recalculated using modern decay constants (Table S1 [footnote 1]); however, small (∼1%) age corrections would be insignificant given the low precision of the measurements. We use the oldest vent located within a 0.5° × 0.5° grid for analysis of migration direction and rates (similar to Crow et al., 2011). Magmatic migration rates were calculated by measuring distances between volcanic sources (e.g., cinder cones and vents) and dividing these distances by the difference in age between flows. The Lovell Wash member basalt flow (13.3 Ma; Harlan et al., 1998) is the oldest flow in the study area with a known vent location and is located in the southwestern corner of the study area (Fig. 1). This basalt flow is used as the location and age from which all distances to other vents are measured and elapsed time to other flows is calculated.

Differential Incision

Our incision measurements are designed to understand bedrock incision in the erosional landscape of the CP. We calculate heights of bedrock straths above the modern river and use best age constraints on overlying river deposits to give average bedrock incision rates in m/Ma This method makes the assumption that instantaneous incision initiates once a lava flow crystallizes. Therefore, our calculated incision rates treat basalts like fill terraces in the sense that the time to re-incise back to prior strath positions gets integrated into long-term bedrock incision averages. By calculating bedrock incision at the million-year time scale, we are able to average out 100 ka glacial-interglacial cycles of incision and aggradation over the past several million years to arrive at better estimates of river and canyon incision. To analyze rock type controls, we assume all unit descriptions in published geologic maps that include “cliff-forming” rock types represent resistant lithologies, while “slope-forming” rock types represent less resistant lithologies.

We use recalibrated K/Ar and 40Ar/39Ar ages of basalts that cap perched gravels of the ancestral VR to quantify the most reliable or “preferred” incision rates (Table S2 [footnote 1]). Other basalt flows have elongated and sinuous outcrop geometries indicating they flowed down paleodrainages, even if river gravels have not been mapped directly beneath the flow. The majority of heights used in incision rate calculations for the VR and its tributaries were estimated from 1:24,000-scale quadrangle maps and Google Earth, while others were measured in the field using a laser range finder. Following the methods of Crow et al. (2014), heights measured from quadrangle maps and a laser range finder are given uncertainties of ±10 m and ±2 m, respectively. Heights at “preferred” locations were measured from the top of the bedrock strath beneath the ancestral river gravels to the current river or tributary elevation to quantify bedrock incision. “Approximate” rates from flows were measured from the contact between the basalt flow and the underlying bedrock to the current river or tributary elevation as used by Willis and Biek (2001; Grant Willis, 2016, personal commun.).

Structures such as hangingwall anticlines and footwall rebound can form above both listric and planar faults (Seixas et al., 2015) that significantly dampen measured incision rates immediately next to faults. Crow et al. (2014) estimated a distance of 10 km from fault traces to be sufficient to eliminate differential incision observed along the Colorado River due to these localized flexural influences and more closely record regional block uplift amounts. Along the VR, we observe a stabilization of incision rates ∼5 km away from major faults. Thus, average incision rates within structural blocks were calculated using the slope of a best-fit linear trend line among the incision points that lie ≥5 km away from major faults. In the Zion block, decreased incision rates near the headwaters of VR tributaries were excluded from regional uplift calculations because these points represent geomorphic controls, because headwater reaches do not have enough stream power to fully incise and represent block uplift. A decrease in incision rates near headwaters is expected and observed in all tributaries. One outlying data point along Fort Pearce Wash with an abnormally low incision rate (incision point 1, Table S2 [footnote 1]) in comparison to other rates on the Zion Block was also excluded from the average block calculations. The St. George and Hurricane blocks were combined when calculating averages due to the lack of basalt flows that exist in the Hurricane block >5 km from any major fault. In the Hurricane block, only one existing flow complies with these criteria, and therefore a trend line could not be drawn to calculate the slope.

Provenance Analysis

We collected three samples for detrital mineral dating from first-arriving ancestral VR gravels of the upper Muddy Creek Formation. Two of these samples are from locations that had previously been analyzed for detrital zircons (Forrester, 2009; Muntean, 2012; Dickinson et al., 2014) and by clast counts (Forrester, 2009; Muntean, 2012). 40Ar/39Ar dating of sanidine grains instead of zircons allows for precise ages and can help understand the exact source of ancestral river deposits (Hereford et al., 2016). From these samples, sanidine grains were concentrated using heavy-liquid mineral separation techniques, and ∼150 sanidine grains were handpicked from each density separate while immersed in wintergreen oil. Sanidine grains were distinguished from plutonic K-feldspar by their optical clarity and lack of microtextures when viewed in the wintergreen oil. 40Ar/39Ar dating of individual sanidine grains was performed by single-crystal laser fusion with a CO2 laser and the gases were measured on an ARGUS VI noble gas mass spectrometer at the New Mexico Geochronology Research Laboratory. Additional details can be found in Table S3 (footnote 1).

Profile Analysis

The nested longitudinal profile of the VR and its major tributaries and knickpoints are shown in Figure 2 and are divided into three groups of rivers: (1) main-stem VR (including Deep Creek and North and East Forks), (2) Colorado Plateau –Transition Zone (CP-TZ) tributaries, and (3) BR tributaries. These groups are described here and will be interpreted in terms of age and history later in this paper.

Crow (2012) noted a difference in average ksn values between Deep Creek (∼130) and East Fork (∼90) tributaries of the upper VR and explored possible causes for differences in steepness. We analyzed these and two additional tributaries, La Verkin Creek and Fort Pearce Wash, for possible spatial associations with climatic, rock type, and tectonic parameters. Figure 4 shows color-coded profiles for mean annual precipitation (Fig. 4A). La Verkin Creek, Deep Creek, and East Fork Virgin River are all south-flowing tributaries sourced in higher-elevation, wetter regions near Zion National Park and show similar precipitation amounts; Fort Pearce Wash drains the lower and drier Hualapai Plateau on the south side of the main-stem VR. There is no observed correlation of higher normalized channel steepness with annual precipitation. Figure 4B shows changing lithology along the profiles in terms of rock erodibility as manifested by cliff-forming versus slope-forming rock types. Resistant cliff-forming lithologies seem to be associated with knickpoints and steeper tributaries to some extent. Deep Creek has a steep gradient and primarily consists of cliff formers, while Fort Pearce Wash has a shallow gradient and primarily consists of slope formers. However, both East Fork (shallow gradient) and La Verkin Creek (steep gradient) have a mixture of lithologies. Figure 4C shows mantle p-wave velocity anomalies. Mantle velocity differs greatly below the steeper gradients of La Verkin Creek and Deep Creek, which are underlain by low-velocity mantle and the shallower gradients of Fort Pearce Wash and East Fork Virgin River, which are underlain by high-velocity mantle.

A regional plot of ksn versus Vp along the main-stem VR and the CP-TZ tributaries is shown in Figure 5A. Three separate trends can be regressed when looking at maximum ksn values in three subsets of data. All three regressions show an increase in maximum ksn as mantle velocity decreases, similar to that shown by Crow (2012). Highest ksn points reflect steepest portions of knickpoints as the VR and its tributaries cross major faults or lithologic contacts. The intermediate regression reflects steep reaches between knickpoints, and the bulk subset of data shows modest increase in channel steepness in areas above low-velocity mantle. Figure 5B considers lithologic erodibility and shows that, rather than high ksn values being associated with cliff-forming bedrock, there is a wide range of ksn values in both slope and cliff-forming rock types. Figure 5C also shows a lack of strong correlation between normalized channel steepness and annual precipitation amounts. The precipitation data set gives annual inches of precipitation to the nearest five inches. For example, 0–10 in/a includes all locations labeled with 0, 5, or 10 in/a. A summary table of ksn values on the main-stem VR and the CP-TZ tributaries (Fig. 5D) shows that the average ksn in tributaries underlain by lower-velocity mantle (δVp < −1.9) is ∼70% greater (100.85) than the average ksn in areas underlain by higher-velocity mantle (δVp > −1.89), where ksn averages 59.07. Average ksn is also higher in cliff-forming (90.72) versus slope-forming (70.15) rock types; however, this is only a ∼30% increase. Differences among areas of high precipitation (15–30 in/a) and low precipitation (0–10 in/a) show only a 5% increase from average ksn values of 82.18 and 77.37, respectively. Two sample t-tests indicate the difference of means found within the mantle velocity and lithology data is statistically significant; however, the precipitation means are not. As a non-parametric alternative, rank sum tests were also performed on the data sets and resulted in similar p-values to those calculated from the two-sample t-tests (mantle velocity = 2.2 × 10−16; lithology = 2.2 × 10−16; precipitation = 0.858). Therefore, ksn is inversely correlated with mantle velocity and positively correlated with lithology; the stronger correlation is with mantle velocity. Figure 6A plots average rock tensile strength from Bursztyn et al. (2015) against average ksn values for seven Paleozoic and Mesozoic units within the VR drainage and shows no correlation (R2 = 0.004, p-value = 0.893). However, Bursztyn et al. (2015) acknowledge the unexpectedly low tensile strength measured in the massively-bedded, resistant, cliff-forming, Navajo Sandstone. Perhaps more data are needed to properly represent rock strengths. Within their data set, approximately ten measurements were made at each sample location with one calculated average tensile strength. Figure 6A shows the minimum and maximum tensile strength values for formations in which samples were collected at more than one location. In contrast to the poor correlation seen in Figure 6A, Figure 6B shows a strong correlation (R2 = 0.98, p-value = 1.63 × 10−5) between mantle velocity and ksn regionally when Vp is plotted against binned average ksn for different mantle velocities in 0.25% bins.

Magmatic Sweep

An analysis of ages and locations of basalt flows in northwestern Arizona, southern Nevada, and southwestern Utah shows a northeast-trending migration of Cenozoic volcanism (Fig. 7). Many flows have long runouts along paleodrainages such that we plotted only vent locations. Contouring the oldest vent ages (yellow stars) within each 0.5° × 0.5° grid shows a strong northeastward-younging trend that is also seen when all known vent ages (red dots) are contoured. Best et al. (1980) and Wenrich et al. (1995) calculated a NE-migration rate of 12 km/Ma in western Grand Canyon. This study includes basalt flows throughout the Lake Mead area and the eastern lobe of the VR drainage (SW Utah) to produce new migration rates. Our migration rates calculated from all vent locations and oldest vent locations (first-arriving flows) in each subregion are higher and range from 15.8 to 17.8 km/Ma, with R2 values of 0.738 and 0.965 and p-values of 2.2 × 10−16 and 4.14 × 10−7, respectively.

Differential Incision

Incision rates vary spatially along the VR and its major tributaries at both local and regional scales. Small transient knickpoints may have swept through the drainage system throughout the past 5 Ma; however, our calculated incision rates show steady incision through time within each structural block. Figure 8 shows incision rates along the profile for two time periods, 0.9–0.2 Ma and 4–1 Ma. Even the shorter-term rates are averaged over >200 ka and hence span more than one glacial-interglacial oscillation; thus, we interpret them to reflect realistic estimates of bedrock incision (cf. Karlstrom et al., 2013; Pederson et al., 2013). Variations in bedrock incision (differential incision) occur dramatically at a local scale, immediately across faults. For example, Volcano Mountain lava flow (0.353 Ma; Sanchez, 1995) was emplaced across the Hurricane fault and provides an excellent example of fault-dampened incision (Pederson et al., 2002). Incision rates in the hangingwall near the Hurricane fault decrease from 110 m/Ma at 5 km west of the fault (incision point 64; Fig. 8), to 0 m/Ma directly adjacent to the fault (still in the hangingwall; incision point 105; Fig. 8) over the same time period. We use the term “apparent” incision rates when flexure near the fault appears to account for variation in incision rates. Our observations are compatible with Crow et al. (2014), who noted such variations due to fault-related flexure are greatest within ∼10 km from the trace of the fault. Calculated incision rates of the same flow (Volcano Mountain) increase to 304 m/Ma on the immediate upthrown footwall (incision point 106; Fig. 8) directly east of the fault. But on the east side, the calculated incision rate of 304 m/Ma (incision point 106; Fig. 8) is similar to the average regional block incision rates (338 m/Ma) of data points farther east, >5 km away from the fault trace, suggesting limited footwall flexure where the VR crosses the Hurricane fault. In contrast, in La Verkin Creek tributary ∼20 km north of the VR, apparent incision rates along the footwall of the Hurricane fault increase to 782 m/Ma (incision point 66; Table S2 [footnote 1]), suggesting substantial footwall flexure at this location.

Given the large number of high-quality incision points made possible by numerous dated basalts in the VR area (Table S2 [footnote 1]), it is possible to separate out differences in average incision rate within blocks using only rates >5 km away from major faults and regress these through time in each block (Fig. 9B). The regression lines in Figure 9B show quasi-steady incision rates within each individual block. Average incision rates, excluding local fault-related deformation and headwater erosion effects, show a regional increase eastward from block to block: from 0 m/Ma in the Lower Colorado River Corridor (Crow et al., 2016), to 23 m/Ma in the Lake Mead block, to 85 m/Ma in the combined St. George and Hurricane block, to 338 m/Ma in the Zion block (Fig. 9A).

Summing the differential incision between blocks multiplied by the duration of fault slip gives a total magnitude of differential uplift between the CP and BR. Duration of fault slip on the youngest fault, the Hurricane fault, is constrained by the 3.6 Ma Bundyville basalt, which has subequal throw as the underlying Moenkopi strata (Billingsley and Workman, 2000). Biek et al. (2010) agree that “most normal faults on the Shivwits and Uinkaret Plateaus became active after 3.6 and possibly after 2.6 Ma.” Using a fault duration of 3.6 Ma gives 1217 m of differential incision (23 m/Ma × 3.6 Ma + 62 m/Ma × 3.6 Ma + 253 m/Ma × 3.6 Ma = 1217 m; Fig. 9), whereas using 2.6 Ma gives 879 m (23 m/Ma × 2.6 Ma + 62 m/Ma × 2.6 Ma + 253 m/Ma × 2.6 Ma = 879 m).

Basalt Paleoprofiles

Basalts flowed into and preserved segments of paleodrainages that inform the evolution of the tributary system for the VR. Grand Wash drains the eastern Virgin Mountains and was filled with basalt at 4.71 Ma (Beard et al., 2007). Basalts overlie gravels that are well rounded, with clasts similar to those found in the upper Muddy Creek conglomerate in the Virgin Depression, including yellow and purple quartzites, yellow volcanics, and chert litharenites that resemble clasts from the Canaan Peak Formation, which outcrops on the east flank of the Pine Valley Mountains (Goldstrand, 1992, 1994). The paleoprofile of the Grand Wash basalts (Figs. 10A and 10B) is subparallel to the present-day profile of Grand Wash. Howard and Bohannon (2001) show that this basalt flow has been downfolded due to the formation of a hangingwall anticline along the Wheeler fault. This tributary incision rate, located >5 km away from the fault trace, varies from 11 m/Ma to 23 m/Ma and is compatible with the Lake Mead block average rate of 23 m/Ma. Crow et al. (2016) reported a higher incision rate of ∼50 m/Ma at the southernmost extent of this flow (Fig. 10B).

The 3.7 Ma Black Rock Mountain basalt erupted from a vent near the divide between VR and paleo–Grand Wash. It flowed northward into Virgin Gorge and southward into Grand Wash Trough and crossed the Grand Wash fault, thus suggesting no post–3.7 Ma fault slip at this location (Figs. 10C and 10D). Figure 10D compares the basalt flow profile to the nearest modern drainage profile. The profiles generally form concave-up longitudinal profiles, and the paleodivide, represented by multiple vents, has remained at a similar location and elevation for the past 3.7 Ma. The slope of the southern drainage has remained similar, and there have been minor amounts of relief production (∼100 m) over the past 3.7 Ma, giving tributary incision rates of 27 m/Ma compared to the Lake Mead block average rate of 23 m/Ma seen in Figure 9. The slope of the northern drainage has steepened with high-relief production in the Virgin River Gorge (∼850 m) over the past 3.7 Ma.

Provenance Analysis

The upper Muddy Creek conglomerate (first-arriving VR gravels) consists of yellow and purple quartzite, volcanic, carbonate, black chert, and chert litharenite clasts. These gravels overlie fine-grained sediments of the Muddy Creek Formation and record the birth of a major high-energy VR system entering the Mesquite basin (Williams, 1996; Forrester, 2009; Muntean, 2012). In one location, a small outcrop (∼15 m × 15 m) of 4.1 Ma basalt near Mesquite, Nevada, appears to underlie these gravels and may imply a 4.1 Ma date for the birth of the VR (Williams, 1996). However, the geologic context of this basalt seemed equivocal to us due to its small size, whether it is a flow or small vent, causing uncertainty about its age relationship to upper Muddy Creek gravels.

Detrital sanidine (DS) samples were collected throughout the Virgin Depression from sand lenses within both the upper and lower Muddy Creek deposits. These data were combined with detrital zircon (DZ) samples collected from previous studies (Forrester, 2009; Muntean, 2012) to form a complete provenance analysis. DS and DZ sample locations are mapped in Figure 11. Potential source locations and ages of the detrital grains are mapped and plotted in Figures 11 and 12A, respectively. Age distribution plots (Fig. 12B) of zircons show distinctly different curves between lower Muddy Creek (sandstone and siltstone) and upper Muddy Creek paleo-VR (conglomerate and sandstone) samples (Forrester, 2009; Muntean, 2012). The lower Muddy Creek curve has a broad peak at 19 Ma with smaller peaks at 14 and 17 Ma. The upper Muddy Creek zircon samples show a sharp peak ca. 20 Ma with a few younger grains at 13 Ma. A Kolmogorov-Smirnov (K-S) test (Dickinson and Gehrels, 2008) between the two data sets gives a calculated p-value of 0.00036, indicating a >95% confidence level that the parent sources of the upper and lower Muddy Creek deposits are statistically distinguishable. Therefore, the upper and lower Muddy Creek deposits came from different sources.

We analyzed DS from three samples of first-arriving VR gravels, of which DZ analyses at two of these locations have previously been published (Fig. 11). An age distribution curve of the three combined DS samples shows four distinct peaks at 13.72 Ma, 18.68 Ma, 20.56 Ma, and 23.8 Ma (Fig. 12C). The precise ages of the two highest peaks define two discrete sources as opposed to the single peak at 20 Ma shown by the DZ data (Fig. 12B). Two grains dated in this study at ca. 5.9 Ma give a maximum depositional age (Dickinson and Gehrels, 2009) of these first-arriving VR gravels.

Profile Analysis

Of the ten identified knickpoints (Figs. 1 and 3; Table 1): A, B, G, and H are interpreted to represent steady-state knickpoints caused by repeated young slip across major faults. Knickpoints C, D, I, and J are interpreted to represent lithologically controlled steps in the longitudinal profile; C is believed to represent the contact of Kayenta and Navajo Formations in Parunuweap Canyon; and D, I, and J are thought to represent basalts in the channel (Fig. 1; Sable and Hereford, 2004; Biek et al., 2010). Knickpoint E formed as a result of the Sentinel landslide, ca. 4.8 ka (Grater, 1945; Castleton et al., 2016). The Sentinel landslide and the resultant knickpoint may hide a knickpoint similar to C, which should be found on the North Fork Virgin River at the Navajo and Kayenta contact.

Possible tectonic influences on channel steepness include both faults and epeirogenic uplift above mantle low-velocity zones. Fault-related knickpoints across the Hurricane and Toroweap faults are not seen in Grand Canyon (Karlstrom et al., 2012; Crow et al., 2014) possibly due to higher stream power of the CR and/or lower slip rates. Landslides are common in rapidly uplifting environments and steep canyons. The largest known landslide in the Zion Plateau area is the Sentinel landslide (Grater, 1945; Castleton et al., 2016). Castleton et al. (2016) shows the knickpoint produced by the Sentinel landslide has been greatly reduced in the past 4800 years as the VR has incised through ∼130 m of the original ∼180 m of debris. Overall, we interpret landslide effects on tributary profiles and incision rates at the millions-of-years timescale to be minor in this region.

VR tributary profiles show evidence of a disequilibrium landscape with deeply incised slot canyons and convex-up knickpoints on the Zion Plateau. Deep, narrow canyons form where the VR North and East Forks incise through the resistant, cliff-forming, Navajo Sandstone. Rock erodibility clearly influences channel steepness (Pederson and Tressler, 2012; Bursztyn et al., 2015). However, our ksn results show a much stronger correlation with low-velocity mantle (Fig. 6). Low-velocity mantle is suggestive of the presence of partial melts in buoyant and rheologically weaker and hotter mantle, which acts as a potential cause for uplift (Karlstrom et al., 2008; Sine et al., 2008; Schmandt and Humphreys, 2010). Therefore, we propose that both rock type and ongoing uplift above upwelling mantle influence the VR profiles with the latter being the dominant driver for the observed difference in steepness between North Fork and/or Deep Creek and East Fork tributaries (Fig. 4C; Crow, 2012).

Magmatic Sweep

The migration path of basaltic volcanism followed a general northeast direction at a rate of ∼18 km/Ma (18 mm/a) from the Lake Mead basalts to the East Fork Virgin River basalts. This parallels the general southwest vector of North American absolute plate motion at a similar rate of ∼20 mm/a relative to the asthenosphere (Minster and Jordan, 1978; Gordon and Jurdy, 1986). The strong negative correlation of steepest normalized stream gradients with underlying low-velocity mantle and the sweep of young basaltic volcanism that involves mixed lithosphere and asthenosphere sources (Crow et al., 2011) leads us to interpret that the sweep of magmatism is related to the southwest movement of the North American plate over upwelling asthenosphere, causing lithospheric removal and/or modification, basalt extraction, and buoyancy-driven uplift above zones of low-velocity mantle.

Differential Incision

Spatially variable and temporally steady incision rates are among key observations that are helpful in determining potential drivers of differential incision. Following the reasoning of Hamblin (1984), we interpret the observed steady incision rates through time measured from the central portions of each block to represent persistent east-side-up regional uplift (Fig. 9). Temporally unstable incision rates, if observed, would suggest transient knickpoints that would be better explained by geomorphic controls. We do observe transient knickpoints at headwater locations and interpret these as locations where the VR is extending its reach into the uplifting Zion Plateau. Incision rates along the uppermost extents of the North and East forks of the VR are 0 m/Ma, but the rates rapidly increase and then stabilize to an average rate for the Zion block only ∼5–10 km downstream.

We report differential incision magnitudes over 3.6 and 2.6 Ma time spans as maximum (1217 m) and minimum (879 m) values, respectively, but we interpret 1217 m of incision as the more accurate magnitude. Both the Lake Mead and St. George and Hurricane blocks show steady rates over 4 Ma timeframes, whereas steady incision is documented only back ∼1 Ma for the Zion block. But there are geologic arguments that steady uplift and/or incision at 338 m/Ma (Fig. 9) can be extended back 3.6 Ma as the VR has headward eroded since it first extended into the Zion block at the Hurricane fault. The VR is presently at an elevation of ∼940 m where it crosses the Hurricane fault. Fault throw measured from central portions of blocks (away from the effects of fault-related flexure) is ∼1100 m (Biek et al., 2010). Assuming slip for the past 3.6 Ma, net fault throw measured from the top of the Navajo Sandstone (1100 m) and calculated footwall incision magnitudes (253 m × 3.6 Ma = 911 m) are in close enough agreement to suggest that the VR first incised through rock that is now at an uplifted elevation of ∼2000 m. The Jurassic Carmel Formation lies at this elevation, ∼20 km east of the Hurricane fault at the “West Temple” (Fig. 1), which is the uppermost stratigraphic layer of the southernmost Zion Canyon. The Carmel and overlying Iron Springs Formations are also the uppermost stratigraphic layers found immediately west of the Hurricane fault. Therefore, an ancestral “Zion Canyon” may have initiated at the Hurricane fault ca. 3.6 Ma when the VR eroded into the Zion block at the stratigraphic level of the Carmel Formation such that the high relief of modern Zion Canyon reflects upstream canyon development in the stronger Navajo Sandstone accompanied by ∼20 km of retreat of Navajo Sandstone (White) cliffs.

Thus, our interpretation of the combined data sets is that a young (4–5 Ma) VR is headward eroding into an uplifting Zion Plateau. Differential incision magnitudes are a proxy for differential block uplift across faults and across the region. Because incision rates can be shown to be quasi-steady over the past ∼6 Ma (Fig. 9B), differential incision rates as well as differential incision magnitude are both proxies for differential uplift. The upper-crustal faulting and the sweep of magmatism implicate tectonic influences on VR evolution, and both are interpreted here to be manifestations of mantle modification processes that are driving melt extraction and buoyancy modification near a step in the lithosphere-asthenosphere boundary that is located beneath the southwestern CP.

Basalt Paleoprofiles

The subparallel profiles in Grand Wash Trough (Figs. 10A and 10B) support the hypothesis of regional active uplift of the Lake Mead Block relative to base level, since base level (CR) has been graded to sea level since ca. 4.5 Ma (Howard et al., 2015). Howard et al. (2015) argue for ∼200 m of regional offset across the Fortification fault system in the Black Mountains located ∼15 km west of the CR-VR confluence (Fig. 1) over the past 4–5 Ma. In contrast, Figure 9 shows an incision rate difference (across this fault) of 23 m/Ma, which would accumulate to 82 m over 3.6 Ma or 104 m over 4.5 Ma (about half the reported regional offset). This discrepancy is relatively small (<100 m over 4 Ma) and may be due to fault-related flexure and/or an added epeirogenic footwall uplift component.

The 3.7 Ma Black Rock Mountain paleoprofile analysis (Figs. 10C and 10D) shows similar results. At the time of this flow, the northern and southern drainages had very similar slopes, which may imply similar base levels, uplift rates, and/or drainage areas. The stark contrast in relief generation between the northern and southern profiles since 3.7 Ma is explained by the northern profile’s connection to the VR (local base level). Since 3.7 Ma, the VR greatly increased its drainage area and scale, while the southern drainage has had no significant change. Since the initiation of the Hurricane fault at 3.6 Ma, the VR has also experienced a large increase in active headwater uplift, relative to smaller amounts of regional uplift observed in the Lake Mead block. The uplifting headwaters steepened the river profiles and increased stream power, which then generated the observed relief.

Provenance Analysis

Provenance analysis is important for interpreting when and how the VR became integrated. The well-rounded, far-traveled clasts found in Grand Wash Trough beneath a 4.7 Ma basalt (Beard et al., 2007) are interpreted to represent reworked Canaan Peak gravels from the eastern flank of Pine Valley Mountains or perhaps from now buried or eroded outcrops that once existed on the western flank. The combination of yellow volcanics, chert litharenites, black argillites with white quartz veins, and maroon and yellow quartzite clasts is diagnostic of the Canaan Peak Formation (Goldstrand, 1992, 1994). The presence of these clasts in the ancestral Grand Wash leads us to believe a paleoriver (ancestral Virgin?) entered Grand Wash Trough from the north ca. 4.7 Ma en route to the newly established Colorado River. Similar clasts found near Mesquite, Nevada, imply an ancestral Virgin exited the Virgin Gorge <4.1 Ma. This age constraint is based on a basalt outcrop dated at 4.1 Ma by Williams (1996); the outcrop is located near the base of a cliff of first-arriving VR gravels.

An age distribution plot of detrital zircons with sharp peaks ca. 19 Ma led Dickinson et al. (2014) to indicate that the ca. 4 Ma arrival of the VR had headwaters in the Pine Valley laccolith (20.5 Ma; Hacker et al., 2007). Dickinson et al. (2014) proposed that the distribution plots would have a much broader peak if the source was from BR igneous centers such as the Indian Peak and Caliente caldera complexes, since they have a much wider range of ages (ca. 33–11 Ma; Best et al., 1993, 2013) (Figs. 11 and 12A). We interpret our analysis of upper versus lower Muddy Creek detrital zircons to indicate a change in headwater location. The lower, internally drained, Muddy Creek had a source in the northern BR igneous centers (primarily the Caliente caldera complex), which is supported by a broader age distribution (Fig. 12B). The upper Muddy Creek detrital zircons, first-arriving VR gravels, show a sharp peak at ca. 20 Ma, implying the capture of an ancestral VR through the VR Gorge with headwaters at Pine Valley Mountains, which contain the 20.5 Ma Pine Valley intrusion.

The more precise data set of detrital sanidine ages allows us to pinpoint source regions and resolve the question discussed in Dickinson et al. (2014) regarding a 19 Ma peak in the zircon age distribution plots. An age distribution curve of the upper Muddy Creek detrital sanidine grains shows two distinct peaks at 18.68 Ma and 20.56 Ma. We interpret the 20.56 Ma age to represent the influx of sediments derived from the Pine Valley laccolith and latite (20.5 Ma; Hacker et al., 2007) located east of Grand Wash fault. Other peaks at 18.68 Ma and 23.8 Ma represent the Caliente caldera complex (Best et al., 1993, 2013), while the 13.72 Ma peak is sourced from Mineral Mountain (14–10.2 Ma; Hacker et al., 2007), part of the Iron Axis intrusives and extrusives (Fig. 12). Therefore, the rapidly arriving, far-traveled gravels that represent the upper Muddy Creek had a mixed population of grains sourced from (1) Beaver Dam wash draining the Caliente caldera complex and Mineral Mountain and (2) an ancestral VR draining the Pine Valley laccolith. The large peak of Pine Valley age grains found among the lowest (first-arriving) gravels within the upper Muddy Creek is interpreted to indicate that an incipient VR headward eroded from Mesquite Basin and tapped into a well-established drainage system with headwaters in the high-elevation Pine Valley laccolith. A small headward-eroding stream into the Pine Valley laccolith would show a gradual increase of 20.5 Ma age grains from lower (KCW17-4 and KCW17-6) to stratigraphically higher Muddy Creek samples (KCW17-19); such an increase is not seen. The two youngest grains (ca. 5.9 Ma) provide a much younger maximum depositional constraint than the previously published youngest detrital zircon grains of ca. 11 Ma (Muntean, 2012). This indicates that the Virgin River Gorge is a relatively young canyon, similar in age to the Grand Canyon.

Our model for VR evolution (Fig. 13) accommodates the following constraints. (1) Detrital zircon data suggest CP (or CP-equivalent formations of the BR)–derived sediments entered the Virgin Depression and Overton Arm basins by at least 6 Ma (Muntean, 2012; Dickinson et al., 2014). (2) A paleoriver (likely the paleo-VR) entered the GWT from the north 4.7 Ma, depositing reworked Canaan Peak gravels from an inferred source east of Pine Valley Mountains. (3) A sharp contact exists in the Overton Arm, Mormon, and Mesquite basins between the informally named lower and upper Muddy Creek Formations. Lower Muddy Creek Formation is fine-grained siltstone interpreted as internally drained basin deposits. Upper Muddy Creek Formation contains coarse gravels that represent the arrival of the VR (Williams, 1996; Swenberg, 2012) at 6–4 Ma (Dickinson et al., 2014; this study). (4) The northward and southward flows of the Black Rock Mountain basalt (3.7 Ma) suggest an ancestral river gorge and drainage divide existed 3.7 Ma in a similar location to the present-day gorge and divide (Figs. 10C and 10D; this study). (5) A 3 Ma tuff intercalated with Muddy Creek Formation in Glendale Basin, west of Mormon Mesa (Fig. 1), was interpreted to indicate that the BR tributaries (White River and Meadow Valley Wash) did not become integrated with the VR until after 3 Ma (Dickinson et al., 2014). (6) Differential incision data are interpreted to suggest a headward progression with slow incision rates since 6–4 Ma near the CR-VR confluence, suggesting stable base level (with the CR graded to sea level) since then in the lower Colorado River corridor; incision rates increase upstream in steps across faults (Fig. 9). (7) The shape of the modern VR watershed (Fig. 1) is unique with two distinct lobes, west and east, which appear to have been two separate drainage basins divided by the Beaver Dam and Virgin Mountains until integration through the Virgin Gorge. (8) The pattern of south-flowing (fault-controlled) tributaries entering an east-propagating main stem is suggested by modern geometries (Fig. 1).

Pre–6 Ma Landscape

Prior to the integration of the Colorado River to the Gulf of California and the incision of Grand Canyon, the Lake Mead region consisted of large internally drained basins throughout the late Miocene (Fig. 13A). Basins relevant to this research include Grand Wash Trough, Virgin Depression (Mesquite and Mormon basins), Overton Arm, Temple Bar, Greggs, and Glendale basins. These basins formed at the major onset of extension ca. 17 Ma (Faulds et al., 2016). The Hualapai Limestone provides a record of a spring-fed lake (Lake Hualapai) in Grand Wash Trough that persisted from 13 Ma until the 6–4.5 Ma arrival of the CR (Spencer et al., 2001; Crossey et al., 2015). The Grand Wash cliffs were the major topographic feature providing relief in the area; the Pine Valley Mountains acted as the major topographic feature in the north.

Integration and Birth of the Virgin River

Reworked Canaan Peak gravels in Grand Wash Trough indicate that streams drained the Pine Valley Mountains by ca. 4.7 Ma (Fig. 13B). Thick deposits of gypsum in the northern Grand Wash Trough provide evidence for internal drainage until ca. 5 Ma (Faulds et al., 2016). A stream entering Grand Wash from the north likely eroded headward until reaching the Canaan Peak Formation exposed on the southeastern flanks of the Pine Valley Mountains. Also, note a stream flowing from the north into the Mesquite Basin at ca. 6 Ma, bringing in ca. 24–11 Ma detrital zircon and/or sanidine grains from the Caliente caldera complex. This explains the 23.8 Ma sanidine grain peak age (Fig. 12C). Caliente caldera complex grains could also be coming from tephra deposits throughout the region (Best et al., 2013). Figure 13C shows that headward erosion from the Virgin Depression across the Piedmont and Grand Wash faults produced the Virgin Gorge and the capture of the Grand Wash drainage by ca. 4 Ma. Stream capture led to the deposition of reworked Canaan Peak gravels at the mouth of the Virgin River Gorge and within the Overton Arm basin. From 4.0 to 3.6 Ma, the large influx of sediment from incision of the Virgin Gorge causes major aggradation within the Mesquite, Mormon, and Overton Arm basins similar to the Bullhead Alluvium located in the lower Colorado River (Howard et al., 2015). At ca. 3.6 Ma (Fig. 13D), slip along the Hurricane fault initiates as an upper-crustal response to asthenospheric upwelling beneath the edge of the CP, as tracked by migrating basaltic volcanism. Combined footwall uplift and epeirogenic doming uplifts the headwaters, increasing the slope, which increases stream power and triggers a transition from aggradation to incision within the lower VR drainage. By 2 Ma (Fig. 13E), the main stem of the upper VR has incised headward across the Hurricane fault and continues to propagate eastward into the uplifting Zion Plateau.

Therefore, a few conclusions can be drawn from the proposed paleogeographic evolution. First, recent and ongoing differential uplift of the Zion Plateau has been the major driver for drainage evolution in this region in the past 5 Ma. Second, headward erosion triggered by uplift along the Piedmont fault was the primary mechanism for the birth of the present-day VR and the transition from internally drained basins to a major through-flowing drainage system.

Steady versus Non-Steady Incision

Steady incision over long time periods at a given location on a river often takes place in tectonic settings where there is a persistent forcing such that bedrock incision rates are subequal to uplift rates (Whipple, 2004). In contrast, variable incision at a given location and transient incision pulses are often driven by geomorphic forcings such as integration or piracy events or rock strength variation (Cook et al., 2009; Darling and Whipple, 2015). The “Sadler effect” (Sadler, 1981; Gardner et al., 1987) expresses a potential timescale bias in interpreting sediment accumulation rates over long periods—for example, if unconformities are not adequately characterized. This concept has also been applied to river incision because bedrock incision processes can be unsteady due to river dynamics and potential long hiatuses in river incision such that incision records can be unsteady even over long time intervals (Finnegan et al., 2014). Our data show an increase in average incision rates with a decrease in average timescale from west to east, block to block, but this is true at overlapping timescales and is interpreted to reflect differential incision of tectonic blocks rather than an unrecognized time-dependent “Sadler effect.”

A potential problem in calculating incision rates is the use of a constantly changing modern river profile (e.g., incision, planation, and aggradation) as the datum for calculating incision rates (Gallen et al., 2015). However, the calculation of rates through time via strath-to-strath measurements in a given reach (e.g., Crow et al., 2014), use of paleoprofiles frozen beneath far-traveled basalt flows, and construction of paleoprofiles from dated strath elevations along the entire river system (e.g., Karlstrom et al., 2008) obviate arguments about the relatively minor variations in incision interpretations that may occur, for example, as a result of using present river level versus present depth to bedrock as a datum (Karlstrom et al., 2007). In this paper, similar to Karlstrom et al. (2013), by emphasizing long time periods that average out the known complexities introduced by glacial-interglacial aggradation and/or incision cycles, we think the resulting differential bedrock incision data set is robust.

Fault-Dampened Incision

The fault-dampened hypothesis (Pederson et al., 2002) is that incision rates and/or magnitudes on upthrown blocks are subequal to the sum of incision rates and/or magnitudes on downthrown blocks plus fault slip rates and/or magnitudes. As discussed above, this hypothesis is modified in this paper (and in Crow et al., 2014) by using data points at sufficient distance from faults (>5 km) to minimize effects of hangingwall anticlines and footwall uplift. Using incision data points from central positions in fault blocks, accompanied by fault slip measured in the same way (Anderson and Christenson, 1989), is considered here to accurately reflect block uplift across faults. In our case, since the lower CR side of the faults preserves a sea level datum, the differential incision is used to infer headwater surface uplift of the Colorado Plateau.

Examining individual faults, uncertainties remain about the total slip and slip history across different segments of the Hurricane and other CP-bounding faults. An estimate of total east-up throw across the Hurricane fault in the Virgin Canyon area, using offset of the top of the Navajo Sandstone at sufficient distance from the fault to be representative of block interiors (i.e., to eliminate local fault-related flexure), is ∼1100 m (Anderson and Christenson, 1989; Biek et al., 2010). Much of this displacement took place after 3.6 Ma as shown by subequal throw of the Bundyville basalt and the underlying Moenkopi strata. Washington, Grand Wash, Piedmont, and Wheeler faults also have Pliocene to Quaternary slip such that the overall CP-BR boundary is a complex zone of anastomosing neotectonic fault strands with displacement migrating inboard (Moore, 1972; Bohannon, 1993; Biek et al., 2010).

Figure 9 shows that the observed differential incision between centers of St. George and Hurricane versus Zion blocks (253 m/Ma) multiplied by an estimated 3.6 Ma of proposed fault slip would equal 911 m of block displacement at a rate of 253 m/Ma. This is similar to the measured amount of total block displacement across the Hurricane fault in this region of 1100 m and is in the range of young fault slip rates of 210–570 m/Ma calculated from offset basalts ranging in age from 350 to 850 ka (Biek et al., 2010). We conclude that the fault-dampened incision model (Karlstrom et al., 2007, modified by Crow et al., 2014) explains the combined differential incision and fault slip history data to a first order.

More refined data and geodynamic models are needed to test the mechanics of cumulative fault slip plus isostatic unloading and uplift of footwalls. Thompson and Parsons (2017) addressed this problem for the BR Province where multi-seismic uplifts of the footwall seem too small to explain the breadth of mountain blocks. They propose a model where elastic rebound causes unloading of the footwall and inflow of lower crust and mantle that, in turn, cause postseismic broad uplift of the footwall. But to explain the observation that both footwall and hangingwall blocks have been incised relative to basalt-preserved paleoprofiles, western CP models (both Grand Canyon and VR) need to also consider the effects of additional buoyancy from mantle-driven subregional epeirogenic block uplift interacting with upper-crustal faulting, where each amplifies the other.

No-Uplift Hypotheses for the Colorado Plateau

In the debate about the scale of post–10 Ma surface uplift of the CP, an alternative hypothesis (e.g., Bursztyn et al., 2015) is that CP uplift is not required to explain the data. Instead, a “no post–10 Ma uplift” interpretation might posit that the elevation difference between the CP and the BR was established during or before slip on the Grand Wash fault in the Miocene, ca. 17 Ma, and is being now incised by the ca. 5 Ma integrated CR system. In this hypothesis, similar to Pelletier (2010), a headward wave of incision with evolving profiles is proposed which shows a migrating knickpoint with fastest incision rates as the knickpoint passes (Cook et al., 2009; Abbott et al., 2015). However, this hypothesis predicts non-steady rather than the steady incision we observe over millions of years all along the profile. The no-uplift interpretation also provides no explanation for the combination of documented fault slips, the sweep of basaltic volcanism, and association of steep stream profiles underlying low-velocity mantle; and it seems less complete. It might be argued that progressively higher young incision rates and magnitudes for inboard blocks reflect Pliocene–Pleistocene downdropping of each block to the base level set by the lower Colorado River west of the Fortification fault (Callville fault of Howard et al., 2015) envisioning fixed elevations of CP headwaters and BR base level since 5 Ma but still accommodating km-scale post–5 Ma faulting across the transition zone. In this view, high relief of the Virgin Gorge and Zion Canyon might be envisioned as being generated by fault slip, with incision amounts subequal to cumulative fault slip, and by upstream migrating waves of incision that reflect base-level fall with no plateau uplift relative to sea level.

However, criteria for determining whether western (hangingwall) blocks moved down or eastern (footwall) blocks moved up, relative to sea level, were summarized by Hamblin (1984). He showed that when both blocks record net incision relative to the basalt-preserved paleoprofiles of different ages, this indicates that both blocks persistently went up relative to the paleoprofiles. Our data also indicate long-term steady incision at different rates across faults (Fig. 9B) with no observed long-term aggradation on downdropped western blocks. This is further supported regionally by evidence farther downstream, west of the Fortification fault (Fig. 9), that the lower Colorado River has been a low-gradient regional river graded to sea level since at least ca 4.5 Ma when deposition of the Bullhead Alluvium began (Howard et al., 2015; Crow et al., 2016). There has been very little aggradation (or incision) of the modern Colorado River profile relative to 5 Ma paleoprofiles (model A of Karlstrom et al., 2008). Thus, like Hamblin (1984), we prefer the interpretation that the eastern blocks were active uplifting blocks across each of the faults and that cumulative slip across the faults reflects mainly headwater uplift of the CP rather than base-level fall of western blocks.

Mechanisms for Mantle-Driven Uplift

The sweep of basalts found within the VR drainage is a compelling data set that links observed differential incision at the surface to subsurface mantle activity through time. Tomographic data of upper-mantle velocity anomalies at 80 km depth show a low-velocity anomaly ring around the western CP boundary and a high-velocity anomaly known as the Escalante anomaly (Schmandt and Humphreys, 2010). Low mantle velocities underlie the majority of the northeastern portion of the VR drainage system including the entirety of Zion National Park. Figure 14 shows a correlation between the location of low-velocity mantle and increased incision rates. We interpret the correlation between high incision rates and low relative mantle velocity to suggest that mantle buoyancy is a contributor to differential incision and differential surface uplift.

Details of a mantle mechanism that could be driving uplift are incompletely understood, but lithospheric delamination and edge-driven convection have both been proposed. The preferred mechanism must be able to explain both the migration of basaltic magmatism toward the center of the CP (Best et al., 1980; Wenrich et al., 1995; Roy et al., 2009; Crow et al., 2011) and an increase in asthenospheric melt through time (Crow et al., 2011). A model indicating delamination below the Escalante anomaly, supported by Levander et al. (2011), suggests thermochemical convection as upwelling mantle generates an intrusion of basaltic partial melts into the base of the CP, increasing negative buoyancy and thermally weakening the mantle lithosphere. Decreased viscosity and increased density of the mantle lithosphere allow a “drip” to form and delaminate the mantle lithosphere and perhaps some of the lower crust. The failure within the lowermost crust along a localized surface, as observed in the receiver function images (figure 3 of Levander et al., 2011), allows for the replacement of lithospheric mantle and lower crust with hot buoyant asthenosphere. The introduction of asthenosphere as the mantle lithosphere continues to tear away (delaminate) from the lower crust is suggested by Levander et al. (2011) to explain the observed migration of magmatism. However, the dipping structure in figure 3a of Levander et al. (2011) opens to the northeast and therefore implies that the hot asthenosphere first encounters the thinned lithosphere in the northeast and migrates to the southwest. Thus, in detail, this model (figure 4 of Levander et al., 2011) does not support the NE direction of basalt migration observed at the surface (Fig. 7).

A second model explaining uplift of the CP proposes small-scale, upper-mantle, edge-driven convection. The CP lithosphere is modeled to have a thickness of 120–140 km, whereas the BR lithosphere has a thickness of 60–80 km (West et al., 2004). Van Wijk et al. (2010) use numerical modeling of an assumed step in lithospheric thicknesses to better understand the convective processes that take place. In this model, a thermally induced Rayleigh-Taylor instability forms at the large step with hot asthenosphere juxtaposed with colder mantle lithosphere. The hydrated and weakened mantle lithosphere begins to drip off the base of the CP, and thinned lithosphere is replaced by hot buoyant asthenosphere, which drives surface uplift. The edge-driven convection model applied to a thicker CP lithosphere as it moves SW over warm asthenosphere may best explain the northeastward migration of magmatism and increased uplift of the Zion Plateau (Fig. 15).

The differential incision history of the 5–4 Ma Virgin River drainage system is interpreted to provide evidence for ∼880–1200 m of surface uplift of the CP. Documentation of steady incision for each reach, and removal of local effects of fault-dampened incision directly adjacent to faults, allow us to equate differential incision magnitude with uplift magnitude. Averaged incision rates show an eastward-propagating, stair-stepping increase across structural blocks: 0 m/Ma in the lower Colorado River corridor, 23 m/Ma in the Lake Mead block, 85 m/Ma in the St. George and Hurricane blocks, and 338 m/Ma in the Zion block (headwaters). Differences in incision magnitude accumulate upstream to ∼880–1200 m of uplift in the western CP, relative to the lower Colorado River corridor, depending on whether Zion block has been incising for 2.6 or 3.6 Ma. Thus, ∼30%–50% of the total surface uplift of the CP since 70 Ma occurred in the past ∼5 Ma. Increased incision rates and oversteepened channel segments (e.g., North Fork Virgin River) correlate better with areas of underlying low-velocity mantle than with variances in substrate lithology or precipitation. An observed NE-propagation of basaltic magmatism at rates of ∼18 km/Ma is similar to North American absolute plate motion (∼20 km/Ma). We conclude that NE-propagating, upper-mantle convection drove differential uplift of the western CP in the past 5 Ma and that the birth and evolution of the VR provide evidence for large-scale landscape response to mantle-driven uplift.

Analytical support was in part from National Science Foundation (NSF) Division of Earth Sciences (EAR) grants EAR-1348007 from the Tectonics Program and EAR- 1545986 from the Sedimentary Geology Division. We acknowledge the long-term efforts by the Utah Geological Survey to perform 40Ar/39Ar dating of all basalt flows in southwest Utah. Their robust and precise data set greatly influenced this research. We thank Andres Aslan, Bob Biek, Keith Howard, and an anonymous reviewer for their constructive criticism, which greatly improved this manuscript. Informal review and discussion with Kelin Whipple also helped improve the paper. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

1Supplemental Tables. Table S1 contains compilation of all dated basalt flows throughout Virgin River region. Table S2 contains acquired data used in the calculation of incision rates. Table S3 contains analytical data, methods, and instrumentation regarding the dating of detrital sanidine grains. Please visit or access the full-text article on to view the Supplemental Tables.
Science Editor: Raymond M. Russo
Guest Associate Editor: Andres Aslan
Gold Open Access: This paper is published under the terms of the CC-BY-NC license.