The causal mechanisms for the onset and patterns of post-Miocene erosion of the western Great Plains remain the subject of an enthusiastic debate concerning the roles of climatically modulated geomorphic parameters and tectonic rock uplift as drivers of long-term erosion. This study distinguishes between these drivers on the plains of New Mexico and Colorado, where post-Miocene erosion and late Cenozoic volcanism of the Jemez lineament have produced distinctive modern landscapes characterized by deep bedrock canyons and inverted, basalt-capped mesas. The 40Ar/39Ar ages of basalt-capped paleosurfaces define an episodic eruption history in the Raton-Clayton and Ocate volcanic fields and help to quantify patterns, amounts, and rates of differential erosion. Several data sets indicate patterns of NE-trending geologic features that require explanation, including: (1) crude volcanic “belts” of similar age, (2) parallel NE-oriented erosional escarpments retreating toward the NW, (3) differential denudation rates increasing systematically NW from a NE-trending hinge line of “low to no” erosion on the Great Plains, (4) a NE-trending zone of broad (50–100 km) convexities in stream profiles identified by an analysis of strath terraces and channel steepness (ksn), (5) reorganization of stream networks that took advantage of an apparent relative base-level fall in the SE, and (6) an ∼150-km-long, 40Ar/39Ar-dated composite paleosurface that has been tilted 64 millidegrees/Ma since 3.4 Ma. Our synthesis of new 40Ar/39Ar geochronology with new calculations of regional surface denudation, channel steepness, and tilt rates shows that post-Miocene patterns of landscape evolution are best interpreted as related to dynamic uplift along the NE-trending Jemez lineament, with second-order impacts from climate-driven geomorphic variables. We interpret this dynamic uplift to be ultimately due to changing mantle structure and buoyancy with associated crustal melt flux.


Drivers of Post-Miocene Landscape Change

The western Great Plains have seen a dramatic change since the late Miocene from an environment of widespread net aggradation associated with Ogallala Group deposition to one of widespread net incision and landscape lowering that continues to the present. Evidence for deep, post–6 Ma incision has also been observed throughout the Rocky Mountains and Colorado Plateau, and an enthusiastic debate has developed in the literature regarding the causal mechanisms driving this phenomenon. The debate largely hinges on the relative roles of tectonic rock uplift (Trimble, 1980; McMillan et al., 2002, 2006; Leonard, 2002; Eaton, 1986, 2008; Karlstrom et al., 2012; Cather et al., 2012) and climatically modulated geomorphic parameters such as changes in water and sediment fluxes (Molnar and England, 1990; Gregory and Chase, 1992, 1994; Chapin, 2008; Wobus et al., 2010) and drainage reorganization as drivers of incision. Because the physical effects of tectonic versus climate end-member models are difficult to distinguish in the geologic record, workers face the challenge of identifying criteria that differentiate climate-driven patterns of landscape evolution from patterns driven by tectonic rock uplift. Furthermore, studies advocating either end-member model must identify their mechanism and demonstrate that the timing of their mechanism coincides with the onset of rapid landscape change.

Recent approaches to this question for the Great Plains region have included analyses of warped paleosurfaces (McMillan et al., 2002; Leonard, 2002; Duller et al., 2012), regional patterns of incision (McMillan et al., 2006), sediment transport models (Wobus et al., 2010), and the timing of regional-scale climatic and sedimentation events as proxies of upstream erosion (Chapin, 2008). All of these studies have concluded that both climate and tectonics may have played important roles in the turnaround from net aggradation to net lowering on the Great Plains since ca. 6 Ma, but their relative roles remain controversial.

New Approaches to an Old Hypothesis

In settings across the western United States, especially in the Rocky Mountains and Colorado Plateau, the hypothesis for late Cenozoic tectonic rock uplift has been revitalized by increasingly well-resolved EarthScope tomographic images of mantle-velocity structure (Dueker et al., 2001; Frederiksen et al., 2001; Goes and Van der Lee, 2002; Schmandt and Humphreys, 2010) that reveal strong multiscale heterogeneity. Large seismic velocity contrasts seen in these images are variously interpreted in terms of different scales and mechanisms of mantle flow (Moucha et al., 2009; van Wijk et al., 2010; Crow et al., 2011; Levander et al., 2011), but all suggest that the sharp mantle-velocity variations observed in tomographic images represent large changes in temperature and rheology sufficient to drive young and ongoing mantle flow. Mantle flow, convective heat transfer, and melt flux, in turn, are modeled to drive perhaps 0.5–1 km of differential surface uplift (Van Wijk et al., 2010). These seismic velocity structures have been linked with tectonic and magmatic activity, and recent studies have also begun to document a relationship between the position of mantle-velocity anomalies and spatially variable rates of long-term river incision and tilting of paleosurfaces (McMillan et al., 2006; Eaton, 2008; Crow et al., 2011; Karlstrom et al., 2012). Thermochronology and incision studies of the Colorado River system have postulated that incision of the western slope of the Rocky Mountains in the last 6–10 Ma may be due to mantle-driven uplift of the Rockies relative to the Colorado Plateau (Karlstrom et al., 2008, 2012).

The goal of this study is to better understand the relative roles of various causal mechanisms in driving landscape evolution on the western Great Plains. This is accomplished using integrated data sets to decipher whether one, as a primary driver, may dominate over others. The causal mechanisms we consider include: (1) tectonic rock uplift in response to magmatic inflation and/or mantle dynamics; (2) climate change events (e.g., the shift toward large global climate oscillations in the late Cenozoic; Molnar, 2004); and (3) geomorphic parameters (e.g., drainage reorganization, regional evaporite collapse events, and isostatic responses to differential erosion).

Important new data sets in this study include a synthesis of published (Olmsted and McIntosh, 2004) and previously unpublished (Stroud, 1997) 40Ar/39Ar basalt ages from the Raton-Clayton and Ocate volcanic fields to create a comprehensive eruption chronology for the northeastern portion of the Jemez volcano-tectonic lineament. We also present new data on landscape denudation rates based on topographically inverted, basalt-capped mesas, a quantitative synthesis of normalized channel steepness of river profiles, and a new analysis of drainage reorganization in this region. Spatial patterns of erosion are compared to mantle-velocity structure and volcanic age domains to address the hypothesis of late Cenozoic, mantle-driven uplift across the Jemez lineament as presented by Wisniewski and Pazzaglia (2002).

To clarify terminology of this study, we consider “tectonic mechanisms” to include the broad range of deformational and magmatic processes that cause rock deformation and change lithospheric buoyancy structure. For landscape changes, we use “incision” for bedrock incision of streams and “exhumation” for the related displacement of rocks relative to the surface, and we use “denudation” for erosional landscape lowering at a variety of scales. In this terminology: rivers incise, rocks are exhumed, and landscapes are denuded. We use the term “dynamic topography” to indicate topography that is modified by active and/or young tectonics (younger than 10 Ma in our case), and we use the term “mantle-driven dynamic topography” for topography that has been modified by changes in mantle flow and buoyancy within and below the plate (Karlstrom et al., 2012).


Jemez Volcanic and Geologic Lineament

The Jemez lineament (Fig. 1A) was originally defined as a NE-trending alignment of late Cenozoic basaltic-volcanic centers (Mayo, 1958). This zone is also described as a long-lived intracontinental tectonic and magmatic boundary in the lithosphere (Magnani et al., 2004) that crosses four major physiographic and tectonic provinces: the Colorado Plateau, Rio Grande rift, Southern Rocky Mountains, and Great Plains. The southern edge of the lineament is coincident with the southern edge of a lithospheric transition zone between the Proterozoic Yavapai (1.8–1.7 Ga) and Mazatzal (1.7–1.6 Ga) Provinces (Fig. 1B). It has been suggested that this transition zone may contain fertile rocks and may have been a favorable conduit for magma ascent in the Cenozoic, thus providing a control on the distribution of magmatism at the surface (Magnani et al., 2004). This zone has a long history of reactivation (Aldrich, 1986) and is proposed to be characterized by a >100-km-wide zone of active uplift of Earth’s surface (Wisniewski and Pazzaglia, 2002; Karlstrom et al., 2005).

Jemez Mantle Anomaly

Figure 1B shows the relative P-wave velocity (%Vp) below New Mexico and adjacent areas as derived from teleseismic P-wave residual measurements using the ∼70-km-spacing EarthScope Transportable Array (TA) and all available PASSCAL experiment data rendered at 80 km depth (Schmandt and Humphreys, 2010). The late Cenozoic volcanic fields of the Jemez lineament are observed to be underlain by a broad, NE-trending, low-velocity anomaly that is distinct from, and crosscuts, the N-trending low-velocity domains associated with the Rio Grande rift and Southern Rocky Mountains. The sharp velocity contrasts that border the Jemez anomaly at a depth of 80 km (both below the Great Plains and SE margin of the Colorado Plateau) are among the sharpest in the western United States: up to 4% (Vp) over lateral distances of ≤100 km. These contrasts are robust features of these regularized tomographic inversions, which, being preferentially smoothed, may underestimate true velocity variation amplitudes. Using standard velocity-temperature scaling relations (e.g., Cammarano et al., 2003; Jackson and Faul, 2010), these differences in velocity may correspond to as much as 500–700 °C of thermal contrast over lateral distances less than 100 km. If a melt phase is present in the low-velocity domains shown in Figure 1B (e.g., Sine et al., 2008; Schmandt and Humphreys, 2010), the thermal contrasts could be less than 500 °C, but rheologic contrasts would be enhanced.

Several geodynamic models that scale the observed velocity variations to density structure and use reasonable mantle flow laws suggest that these large velocity contrasts may reflect Neogene and ongoing upper-mantle convective flow (Schmandt and Humphreys, 2010) of perhaps 1–5 cm/yr. This flow pressure could, in turn, drive surface uplift of several hundreds of meters (Moucha et al., 2009; van Wijk et al., 2010). Alternatively, or in addition, melt flux through dike and sill networks in these low-velocity regions could cause changes in mass distribution and induce differential isostatic adjustments that could influence topography.

Tectonic History of the Great Plains

The Great Plains consists of a broad topographic ramp that connects the uplifted regions of the Rockies with the low-elevation midcontinent. The cross-sectional geometry of the ramp follows a power-law function decay curve, consistent with thermal cooling over a half-space and suggests broad epeirogenic uplift associated with the generally N-S–oriented Rockies and Rio Grande rift (Roy et al., 2004). This geometry is similar in cross section to mid-ocean ridges, leading to models linking topography to thermal decay symmetrically away from a N-S “ridge” axis called the Alvarado Ridge (Eaton, 1986, 2008). However, both the timing and mechanisms of the proposed uplift are poorly understood for this region and probably interact with dynamic processes such as convective mantle flow and magma transfer. Timing of development of the Great Plains ramp is important for unraveling the relative importance of Laramide, Oligocene, and Neogene components of Rocky Mountain uplift (Karlstrom et al., 2012; Cather et al., 2012).

The interior plains region of North America spent nearly half a billion years, between ca. 570 Ma and 70 Ma, covered by shallow seas that deposited a thick blanket of flat-lying sedimentary rocks (Trimble, 1980). Beginning at ca. 75 Ma, the Laramide orogeny uplifted much of the western United States, including the Great Plains surface, which sloped eastward from the nascent Rocky Mountains. For much of the 70 Ma that followed, the Great Plains acted as zone of net aggradation of sediments being shed eastward off the uplifting Rockies. The Ogallala Group (deposited between ca. 13 and 5 Ma in this region; Cather et al., 2012) is a large coalescent lithostratigraphic body of sandstones and gravels deposited during numerous cut-and-fill cycles. Since the end of the Miocene, net fluvial incision has dominated and continues to the present.

Raton Section

Our field laboratory lies in the Raton section of the western Great Plains (Fig. 2A), a physiographically distinct region of the interior plains located adjacent to the Southern Rocky Mountains that spans >44,000 km2 along the New Mexico–Colorado border at 37°N (Fig. 2A). Here, late Cenozoic volcanism of the Jemez lineament and post-Miocene incision have combined to create modern landscapes characterized by deep bedrock canyons, topographically inverted basalt-capped mesas, and appreciably higher topographic roughness than adjacent areas of the Great Plains (Coblentz and Karlstrom, 2011). From east to west, elevations increase from ∼1000 to over 2500 m, and bedrock geology consists primarily of mildly deformed Mesozoic sedimentary rocks, which are unconformably overlain by remnants of Neogene fluvioclastic units (e.g., the Miocene Ogallala Group). The landscapes in this region are dominated by the flows and vents of the Raton-Clayton and Ocate volcanic fields, which have been episodically active since ca. 10 Ma (Stroud, 1997; Olmsted and McIntosh, 2004; this study). These fields represent the northeasternmost expression of Jemez lineament volcanism.

The diverse volcanic landscapes of the Raton section have been carved into locally high relief by post-Miocene erosion of the Great Plains surface, resulting in the formation of basalt-capped mesas and volcanic edifices that often stand high above the surrounding plains (Figs. 2B and 2C), especially in the northwest. For instance, mesa-capping flows near Raton, New Mexico (“Ra” in Fig. 2A) rest nearly 550 m above the Canadian River. The gradual southeasterly slope of the Raton section is strikingly punctuated in two places by the NE-trending Canadian and Cimarron erosional escarpments (Figs. 2A and 2D). These erosional features rise 300–400 m from base to top and extend for straight-line distances of ∼160 km and 80 km, respectively, imparting to the region a staircase-like topography of retreating cliffs held up by relatively resistant strata (Fig. 3).

Elsewhere, rivers have carved spectacular bedrock canyons into the Mesozoic and late Cenozoic sedimentary rocks that underlie much of the Great Plains. The Canadian River, for example, has carved a 100-km-long, 400-m-deep, and 1-km-wide canyon from Springer to Sabinoso, New Mexico (Figs. 2A and 2E). For comparison, these dimensions make this relatively unknown feature longer, deeper, and nearly as wide as the often-visited Grand Canyon of the Yellowstone at Yellowstone National Park in Wyoming. Other impressive canyons can be found along the Mora and Dry Cimarron Rivers (Fig. 2A). This region is especially well suited to studies of landscape evolution given the abundant age constraints on basalt-capped mesas and evidence for widespread erosion in the late Cenozoic.


Field Methods

The Raton-Clayton volcanic field spans ∼20,000 km2, with eruptive centers in New Mexico and southern Colorado. A general geographic overview of the field is shown in Figure 4, along with place names that will be used frequently hereafter. Sixty-four basalt samples were collected from localities throughout the Raton-Clayton volcanic field to provide a representative sampling of most of the existing volcanic features. At each sample locality, one to two hand samples were collected from a volcanic feature, such as a flow or cone, for 40Ar-39Ar dating. Hand samples consisted of 1–2 kg of fresh, nonvesicular lava or the central portions of volcanic bombs.

40Ar/39Ar Analytical Methods

Hand samples were initially crushed and sieved as a first step in obtaining groundmass concentrates. Groundmass concentrates were prepared from the 1.18–0.850 mm fraction by handpicking fragments devoid of visible phenocrysts under a binocular microscope. For seven samples, concentrates were also prepared from the 300–250 µm, 250–180 µm, and 180–150 µm sieve fractions. The nonmagnetic and highly magnetic portions of all samples were removed with a Frantz magnetic separator prior to irradiation in an attempt to eliminate xenocrysts from the groundmass concentrates.

The sample grains were treated in a solution of dilute HCl and then ultrasonically cleaned in distilled water. Samples were placed in aluminum irradiation trays with alternating flux monitors of the Fish Canyon Tuff sanidine (27.84 Ma relative to Mmhb-1 hornblende at 520.4 Ma; Samson and Alexander, 1987) and were irradiated between 1 and 10 h in the reactor at Texas A&M or the Ford reactor at the University of Michigan. The samples were analyzed 2 to 5 mo following irradiation in order to allow decay of 37Ar, which can degrade mass spectrometer performance. Moreover, three additional experiments were run. In the first of these, some samples were prepared in much larger quantities (400–500 mg) in an effort to increase precision by increasing the 36Ar signal. In the second experiment, variable grain sizes (850–150 µm) were used in an attempt to eliminate excess argon derived from xenocrysts present in the sample. A third experiment was conducted to measure the presence, if any, of excess 40Ar in olivine present in sample JS-94–4. Results and inferences from these additional experiments are reported in Stroud (1997).

Analyses were performed at the New Mexico Geochronology Research Laboratory at the New Mexico Institute of Mining and Technology. When samples were run, this facility included a MAP 215–50 mass spectrometer attached to a computer-automated all-metal argon extraction line. FC-2 sanidine flux monitors were fused by a CO2 laser for 15 s, and reactive gases were removed by SAES GP-50 getters prior to gas expansion into the mass spectrometer. A double-vacuum Mo resistance furnace was used to heat the groundmass concentrate samples. Samples were incrementally heated in typically 11 steps from 500 °C to 1650 °C. Heating times of 8 to 10 min were simultaneous with reaction with an NP-10 getter. Additional gas cleanup was conducted by reaction with a GP-50 getter for 5 to 10 min. Isotopic data are corrected for extraction line blank, which for furnace analyses ranged from 3 × 10–16 to 3 × 10–15 moles 40Ar and from 2 × 10–18 to 9 × 10–18 moles 36Ar.

40Ar/39Ar Results and Volcanic Episodicity

The 40Ar/39Ar results from the Raton-Clayton volcanic field add to the large 40Ar/39Ar data set published for the Ocate volcanic field by Olmsted and McIntosh (2004). The combination of these data sets and the high precision of 40Ar/39Ar ages allow for a more detailed interpretation regarding the timing and areal extent of volcanism in northeastern New Mexico. Results will be discussed in the context of a five-phase division of volcanism for the region that synthesizes a previously recognized three-phase division for the Raton-Clayton volcanic field (Calvin, 1987, and references therein) and a 14-phase division for the Ocate volcanic field (Olmsted and McIntosh, 2004).

Table 1 summarizes the eruptive history of the volcanic field using 40Ar/39Ar ages for the Raton-Clayton field that were initially reported by Stroud (1997). Stroud reported only plateau ages, and, in some cases, we have reevaluated Stroud’s (1997) original preferred ages and chosen to report isochron ages instead. This was done for about one third of the samples, where isochron data indicated the presence of excess argon, such that the isochron age is a more accurate estimate of the correct eruption age of the sample. Use of the isochron ages rather than the plateau ages makes little difference to our erosion rate calculations and interpretations. All age spectra and isochron plots (Figs. DR2–DR72) are located in the Supplemental File.1 Ocate volcanic field ages were reported in Olmsted and McIntosh (2004) and are also discussed herein and enumerated in Table DR1 in the Supplemental File (see footnote 1). Figure 5 shows age probability plots and histograms of age-frequency, and Figures 6 and 7 show the locations and ages for samples from both volcanic fields.

As a proxy for eruption volume, we used geographic information system (GIS) and New Mexico State geologic maps to estimate the land-surface area covered by volcanic units belonging to each phase. Table 2 enumerates absolute and fractional surface areas, and Figure 8 compares these results graphically. General errors of ±20% were applied to all values to account for erosion of flows since emplacement, buried/unmapped flows, and uncertainties in identifying phase boundaries between some adjacent flows. We relied on the detailed flow maps of Olmsted and McIntosh (2004) and Stroud (1997) to delineate phase boundaries between adjacent flows and used dashed polygons where the ages of flows are unconstrained by 40Ar/39Ar ages and/or detailed mapping. Flow areas in dashed polygons are split evenly between the volcanic phases to which they likely belong.

Phase One: Early Volcanism (10–5.25 Ma)

The earliest episode of eruption involved volcanism over relatively small areas between ca. 10 Ma and 5.25 Ma. Vents for these eruptions were located dominantly in the NW portions of both volcanic fields, with some volcanism in the far south. Near Raton, New Mexico, eruptions lasted from 9.04 ± 0.04 Ma to 7.73 ± 0.03 Ma. These alkali olivine basalts filled local topographic lows, probably valley floors and plains, at the time of eruption. The Red Mountain rhyodacite domes on Johnson Mesa were extruded ca. 7 Ma (Scott et al., 1990), and then small feldspathoidal basalt eruptions occurred on Johnson Mesa around 6 Ma (Scott et al., 1990). Flows from this phase account for 14.8% (694 km2) of the total surface area of volcanic rocks in the study area (Fig. 8).

Phase Two: 5.25–4 Ma Volcanism

Eruptions from 5.25 to 4 Ma occurred primarily in the northern parts of both volcanic fields, but they have long lava runouts that extend to the south and east. In the Raton-Clayton volcanic field, vents on Mesa de Maya erupted at ca. 5.10 ± 0.03 Ma. Eruption of the Yates olivine basalt flows continued from phase one until at least 4.67 ± 0.04 Ma. Flows from this phase account for 19.9% (932 km2) of the total surface area of volcanic rocks in the study area. The cumulative surface area of flows in the Raton-Clayton field roughly doubled during this phase, while the cumulative surface area of flows in the Ocate field nearly quadrupled (Fig. 8). Nevertheless, total surface area of flows remained higher in the Raton-Clayton volcanic field.

Phase Three: 3.8–2.58 Ma Volcanism

In the Raton-Clayton volcanic field, eruptions from 3.8 Ma to 2.58 Ma occurred at multiple nonadjacent locations. The only stratovolcano in the field, Sierra Grande (Fig. 4), began erupting at 3.8 ± 0.2 Ma and continued until 2.77 ± 0.04 Ma. Ages obtained from the andesite flows of Sierra Grande indicate an active eruption cycle shortly before eruptive activity near Clayton, New Mexico. The Bartlett Mesa basalts (near Raton, New Mexico; Fig. 4) erupted around 3.37 ± 0.10 Ma. The 3.58 ± 0.02 Ma Pine Butte dome is the only other andesite unit in the Raton-Clayton volcanic field. Volcanism at the Don Carlos hills occurred around ca. 3.5 Ma. In the Ocate volcanic field, vents from this phase are dominantly located in the south; flows of this age also cover the largest areal extent (576 km2) of any phase in the Ocate field. Flows from both fields emplaced during this phase account for 27.7% (1299 km2) of the total surface area of volcanic rocks in the region.

Phase Four: 2.58–1.45 Ma Volcanism

Eruptions from 2.58 Ma to 1.45 Ma occurred primarily in the south and east of both volcanic fields. In the Raton-Clayton field, this phase includes the Clayton flows, an areally extensive amalgamation of flows with long (up to 40 km) runouts. The Clayton volcanics consist of olivine basalts interpreted as having erupted from several fissures, which were subsequently punctuated by younger cinder cones of later eruptions (Dungan et al., 1989). The final eruption of the Clayton phase was at Mount Dora, a shield volcano lying between Clayton and Sierra Grande, at 2.03 ± 0.13 Ma. In the Ocate field, the largest set of eruptions occurred in the far south at Maxson crater, where ca. 1.57 Ma flows poured into the paleo–Mora River canyon and traveled down the Mora River channel over 60 km until coming to rest near the Canadian escarpment (Fig. 7). Flows from this phase account for 29.4% (1379 km2) of the total surface area of volcanic rocks, making it the largest single phase in the region.

Phase Five: Waning Volcanism (Post–1.45 Ma)

The most recent eruptive episode lasted from 1.45 Ma to ca. 50 ka. In the Raton-Clayton field, flows from this phase include the Capulin basalts and scattered feldspathoidal flows (Fig. 6). Capulin basalts are found in the western and central part of the volcanic field and occupy modern valley floors and plains. The exception to this is in the western part of the field, where flows form topographically inverted mesas 30–165 m above local base levels. These western flows had previously been considered as part of an older phase because of their high relief (Scott et al., 1990), but their petrologic properties (Dungan et al., 1989) and 40Ar/39Ar age range place them within this youngest phase. In the Ocate field, flows from this phase originated at vents located in a relatively small geographic area located in the central part of the field. Flows from this phase account for 8.1% (381 km2) of the total surface area of volcanic rocks in the study area, making it the smallest phase overall. It is noteworthy that the surface area of flows of this age in the Raton-Clayton field (335 km2) is more than seven times greater than the surface area of flows in the Ocate field (45 km2) for the same period.

In summary, the Raton-Clayton and Ocate volcanic fields both show evidence for episodic activity since ca. 10 Ma. In every volcanic phase, the calculated eruption area for the Raton-Clayton field is greater than that for the Ocate field, and, as a result, the Raton-Clayton field accounts for 68.8% of the total regional eruption area (Table 2; Fig. 8). Using eruption area as a proxy for magma flux, magmatism in the Ocate field apparently peaked during phase three (3.8–2.58 Ma), whereas the Raton-Clayton field saw increased activity as late as phase four (2.58–1.45 Ma), before magmatism in both fields decreased in phase five (post–1.45 Ma). Despite the long history of volcanic activity in both fields, 57.1% of all volcanic surface area appears to have been formed in the proportionally small time frame of ∼2.4 Ma covered by phases three and four (3.8–1.45 Ma). The intensity of magmatism since 1.45 Ma has been markedly lower than any other volcanic phase examined.


The 40Ar/39Ar geochronology data set for the Raton-Clayton and Ocate volcanic fields shows both NW- and NE-trending patterns of alignment of vents and flows and changing distribution of volcanic eruptions through time. Figure 9 shows both volcanic fields where flows have been grouped by age according to the five-phase division that we synthesized from Stroud (1997) and Olmsted and McIntosh (2004). Figure 9B shows a compilation of all mapped vents from these volcanic fields. Vents with 40Ar/39Ar ages are shown as dots with a black bull’s-eye, while the eruptive phase to which all other vents belong has been inferred from the flow mapping of Stroud (1997) and Olmsted and McIntosh (2004). Figure 9B (and also Fig. 4) shows that vents align in NW-trending clusters within the overall NE-trending Jemez lineament. These are interpreted as an expression of feeder dikes that followed NW-trending faults and other crustal weaknesses related to the underlying Ancestral Rocky Mountains (of Pennsylvanian age) and Sierra Grande arch. We suggest that approximately synchronous volcanism was concentrated in crude NE-trending alignments, or belts, during a given eruptive phase (Fig. 9B). Phase one may exhibit two parallel belts, whereas volcanism of phases two through five was concentrated in between, giving the region a NE-trending zonation, subparallel to the overall Jemez lineament. This pattern suggests a weak volcanic migration trend toward the center of the Jemez lineament in both fields. Although both fields overlapped in time and shared a large eruption episode in phase three (3.80–2.58 Ma), Ocate volcanism then began to wane during early phase four (2.58–1.45 Ma), while Raton-Clayton volcanism remained very active, and many spatially extensive flows were emplaced (Fig. 8). This suggests a NE migration in magmatism and volcanism along the Jemez lineament.


Basalt-capped mesas form “inverted topography” in the sense that flows that once occupied the lowest parts of the local landscape are preserved today as the highest topographic mesas because of basalt’s relative resistance to erosion. Topographically inverted flows in the Raton-Clayton and Ocate volcanic fields are widely distributed in terms of both age and geography, thus providing an excellent opportunity to examine rates and patterns of denudation on a regional scale. Previous studies also recognized this potential and utilized erosion data in their reports on the physiographic evolution of one field or the other (O’Neill, 1988; Stroud, 1997; Olmsted and McIntosh, 2004). However, this study is the first to integrate all data using a consistent measurement technique on elevated basalt flows (n = 28) within both volcanic fields to evaluate patterns of differential denudation rates across the region.

Earlier workers defined broad erosion surfaces based on several (low, intermediate, and high) composite surfaces that bevel across gently folded Paleozoic and Mesozoic rocks. These surfaces are marked by gravel lags containing abundant rounded cobbles of Precambrian granite, gneiss, and pegmatite originating from the Sangre de Cristo Mountains to the west and were interpreted to be upland pediment surfaces related to Ogallala deposition (O’Neill, 1988). Basalt flows that cap the gravels, in general, are older on the higher-elevation surface, and successively younger flows cap progressively lower surfaces. In places, however, young flows also erupted at high elevations to cover older basalts. The older surface is capped by flows ranging from 8.3 to 5.5 Ma and is interpreted to record landscape stability until ca. 5.5 Ma (O’Neill, 1988).

The intermediate surface was named the Urraca surface in the Ocate field (after Urraca Mesa; Fig. 1; Ray and Smith, 1941) and was correlated with ridge tops in the Park Plateau in the Raton area (O’Neill, 1988). This surface is several tens of meters lower and is capped by flows that range from 5.0 to 4.0 Ma (O’Neill, 1988). In places, the Urraca surface is locally convex upward (Fig. 3B) and was interpreted by O’Neill (1988) to record local warping of an originally concave-up pediment after ca. 4–5 Ma. The lowest surface is overlain by 3.3–3.0 Ma basalts (e.g., Charette Mesa; Fig. 7) and is graded to the modern Great Plains surface at the SE part of the Ocate field. By the time of the 1.57 ± 0.02 Ma Maxson crater eruption (Olmsted and McIntosh, 2004), modern drainage patterns were in place, and this flow followed the Mora River to the confluence of the Canadian at a time when the canyons of both these rivers had incised to about half their present depth. This flow documents post–1.57 Ma average bedrock incision rates of 60–70 m/Ma (Wisniewski and Pazzaglia, 2002) and provides the best timing control for post–1.57 Ma warping of Canadian River terraces.

Calculating Long-Term Denudation Rates

Here, denudation refers to local erosion in the vertical dimension, as measured directly adjacent to basalt-capped mesas. Denudation rates were calculated by measuring the elevation difference between the estimated bottom of a mesa-capping basalt flow and the base of the mesa, and then dividing by the age of the flow. We estimated the elevation of the local base of the mesas from digital topographic profiles, defined qualitatively as the first prominent decreasing break in slope along the immediate mesa flanks. The elevation difference between a basalt-capped paleosurface and the local base of a mesa was then measured in the field using a laser rangefinder. This method provides minimum erosion values and hence minimum erosion rates compared to using elevation of more distant stream channels, for example. The 40Ar/39Ar ages on flows come from Olmsted and McIntosh (2004) and this study. Errors on rates of denudation were calculated by propagating the uncertainties associated with 40Ar/39Ar eruption ages and field measurements of elevation via addition in quadrature. Error calculation is further discussed in the Supplemental File (see footnote 1).

In cases where multiple ages and elevations exist for a single mesa (e.g., Charette Mesa), we used the age and elevation of the flow nearest to our area of interest. We avoided the method of measuring denudation to the nearest major drainage because this procedure typically results in significant overestimates of denudation due to the deep, local incision of many stream networks into the modern Great Plains surface. It is therefore important to distinguish canyon incision rates from longer-wavelength incision, or denudation, to draw meaningful comparisons between the elevations of modern and paleo–Great Plains surfaces. Canyon incision rates for the Canadian River are discussed in Wisniewski and Pazzaglia (2002). Note that “average” denudation rates are averaged over time, not catchment area.

Spatial Patterns in Long-Term Denudation

The results of our denudation calculations can be seen in Figure 10A and are enumerated in Table 3, where labeled/listed flows include only those that were accessible for elevation measurements in the 2011 field season. Figure 10A shows that long-term average erosion rates vary widely depending on the location of the flow, with maximum rates observed in the northwest, for example, at Yankee and Bartlett Mesas in the Raton-Clayton field, where rates are >100 m/Ma since 3.4 Ma. Minimum long-term rates are seen in the southern and eastern parts of both volcanic fields, where basalt flows that range in age from ca. 1.6 Ma to 4.7 Ma show little to no erosion since they were emplaced. Our long-term denudation rates usually only differ by a few meters per million years from estimates made by Stroud (1997) for the Raton-Clayton field, and they are typically lower than those presented by Olmsted and McIntosh (2004) for the Ocate field (in one case, our rate is lower by >20 m/Ma). Although these differences are sometimes large, they were also expected given our use of more local base levels, and they nevertheless verify the same general spatial patterns that were evident in previous studies.

The spatial pattern of long-term denudation is striking and shows systematic variability, with the highest denudation rates to the NW and lowest to the SE. Rates decrease systematically in the SE direction toward a NE-trending zone, or hinge line, of low to no erosion. Based on these data, this hinge line is drawn in Figure 10A. It is important to note that, because basalts occupy most of the area between drainages, as well as paleodrainages, these denudation numbers are good approximations of the landscape denudation of different parts of the Great Plains surface itself.

Figure 10B shows graphically that denudation rates decrease toward the SE in both volcanic fields (r2 = 0.60), but that the decrease is more pronounced in the Ocate field (open circles), where there is an almost-linear (r2 = 0.82, not regressed in Fig. 10B) increase in long-term denudation rate as a function of a flow’s projected position along the transect B–B′, which is roughly parallel to the Canadian River. The relationship between the total magnitude of denudation (rather than the long-term rate) and flow position along B–B′ is also relatively strong, where r2 = 0.49 for both fields (Fig. 10C). These plots clearly demonstrate that the high variability of long-term surface denudation in our study area has a systematic distribution, with higher rates to the NW, decreasing to low values along a NE-trending hinge line located in the SE part of both fields. This spatial trend is not well developed in transects oriented in other directions; for example, see Figure DR1 (Supplemental File [see footnote 1]) for an unsystematic plot of denudation values projected onto a W-E–oriented transect line.

The strength of this spatial relationship is further emphasized when we compare denudation rates and magnitudes to the 40Ar/39Ar age of a flow. In steadily eroding landscapes, for example, age-elevation plots would be expected to show average denudation rates through time, or linear trends, for periods of semisteady erosion, regardless of position along the river profile. Instead, Figure 10D shows that for both volcanic fields, there is very little correlation (r2 = 0.02) between any long-term denudation rate and the age of a mesa-capping flow. Figure 10E shows that the relationship between the local magnitude of denudation and flow age is weaker (r2 = 0.38, for data from both fields) than the relationship between denudation magnitude and flow position along the B–B′ transect (r2 = 0.49; Fig. 10C). These results differ from the situation in uniformly eroding landscapes where older mesas are expected to be always located at higher points on the land surface (O’Neill, 1988; Stroud, 1997; Olmsted and McIntosh, 2004). The plots in Figure 10 imply that differential denudation rates are strongly dependent on the NW-SE position of flows on the landscape and that they are relatively independent of basalt flow age.

Temporal Patterns in Long-Term Denudation

We examined temporal variations in long-term denudation rates at four locations where dated basalt flows of variable ages reside in close proximity to one another. The subset of data used in the calculation of rates through time is tabulated in Table DR2 in the Supplemental File (see footnote 1) and shown graphically in Figure 11.

Near Raton, New Mexico (red box, Fig. 10A), four closely spaced flows show that average denudation rates were relatively low (∼10 m/Ma) between ca. 9 Ma and 3.4 Ma, but rates rose sharply (90–114 m/Ma) between ca. 3.4 Ma and the present. In the Mesa Larga region of the Raton-Clayton volcanic field (green box, Fig. 10A), three flows demonstrate that denudation rates between ca. 6.9 Ma and 0.75 Ma were slow (∼30 m/Ma) but accelerated to ∼80 m/Ma between 0.73 Ma and the present day. Although there is a long time gap in the Mesa Larga data set from ca. 6.1 Ma to 0.73 Ma, the available data permit an onset of rapid denudation sometime after 4 Ma, as observed near Raton, New Mexico. Both sites in the Raton-Clayton volcanic field document denudation that has increased unsteadily since the late Miocene, from rates measured in the tens of meters per million years to modern rates measured near 100 m/Ma. Slow rates of denudation before ca. 3.4 Ma in the Raton-Clayton volcanic field are consistent with the existence of a stable, intermediate-level Urraca–Park Plateau surface in this area until at least 4.0 Ma (O’Neill, 1988) and deposition of Ogallala Group sediment that may have continued to some extent in the southern High Plains until ca. 4.5 Ma (Chapin, 2008, and references therein).

Temporal variations in denudation rate in the Ocate volcanic field are constrained at two locations. In the Charette Lakes area (yellow box, Fig. 10A), five flows indicate that the average denudation rate was ∼116 m/Ma between ca. 4.6 Ma and 3.2 Ma and decreased to ∼37 m/Ma after ca. 3.2 Ma. This immediate area does not contain older flows to document the onset of this rapid denudation, but our data are compatible with O’Neill’s (1988) argument for a period of relative stability from 8.3 to 5.5 Ma, during which his high-level Ocate surface formed. This is shown in Figure 11, where we have extrapolated the average denudation rate for the Charette Lakes region (116 m/Ma) up to the elevation of Apache Mesa West at 2350 m (yellow dashed line), which O’Neill identified as part of his high-level Ocate surface.

Near Wagon Mound, at the southern end of the Ocate field, five closely spaced mesas indicate that the average denudation rates were around ∼44 m/Ma between ca. 6.1 Ma and 3.0 Ma and then decreased to ∼16 m/Ma thereafter. Denudation rates and magnitudes near Wagon Mound are low overall because the location is near to the low to no erosion hinge line shown in Figure 10A. At this location, the 6 Ma basalts on Los Mesas del Conjelon rest on the high-level surface of O’Neill (1988), which we have shown with schematic blue dashed lines in Figure 11.

Olmsted and McIntosh (2004) speculated about the difference in landscape denudation history between the two volcanic fields and noted that the period of relatively rapid denudation in the Ocate field seemed to have been earlier (6–4 Ma) than the onset in the Raton-Clayton field (ca. 3.4 Ma near Raton, New Mexico). This time window of 6–3.4 Ma is indicated with the shaded gray bar in Figure 11. Furthermore, they noted that rates in the Ocate field decreased in the last ∼4 Ma in this region, while rates increased in the Raton-Clayton field. They also suggested that the high denudation rates observed from ca. 3.4 Ma to the present near Raton, New Mexico, may reflect headward erosion of streams, because these flows are located near the headwaters of the modern Canadian River. However, the flows near Raton are not substantially closer to the headwaters of modern major drainages in the Sangre de Cristo Mountains than many of the flows in the Ocate volcanic field, and the Mesa Larga flows are substantially farther away. Furthermore, a simple headward erosion model does not explain how a period of moderate to high rates of denudation in upstream portions of the Canadian catchment (e.g., at Charette Lakes between 4.6 and 3.2 Ma) preceded a period moderate to high rates of canyon incision further downstream at the confluence of the Canadian and Mora Rivers (60–70 m/Ma since 1.57 Ma; Wisniewski and Pazzaglia, 2002). This model may also have trouble explaining the higher denudation magnitudes to the NW in both volcanic fields and would not predict the observed bowed river terraces in the Canadian River (Wisniewski and Pazzaglia, 2002) and convex gravel-capped pediment surfaces in the Ocate volcanic field (O’Neill, 1988).

An alternative explanation for this pattern is that the differences in the timing of relatively rapid denudation are determined by magma flux and uplift. By this model, the decrease in denudation rates after 3 Ma in the Ocate volcanic field may reflect the waning of magmatism after phase three (3.8–2.58 Ma; Fig. 8). In contrast, the onset of relatively rapid denudation rates in the Raton-Clayton field by 3.4 Ma coincides with the increase in magmatism there during phase four (2.58–1.45 Ma; Fig. 8). However, better constraints on temporal variation of denudation rates at different localities are needed, and other data sets that compare spatial patterns of denudation with patterns of basaltic volcanism and velocity variations in the potential mantle source should also be considered.

Another possible interpretation of the data in Figure 11 is that there are no actual temporal variations in the onset of rapid denudation between the Raton-Clayton and Ocate fields. Given the large time gaps between dated flows in some areas, the data are permissive of a period of high denudation rates (≤100 m/Ma) that began at all four locations at about the same time (4–4.6 Ma; blue bar in Fig. 11) and could have lasted for as long as 1 Ma. We do not strongly advocate this interpretation, however, as it represents one special case that is not well constrained with the rest of the available data.


In addition to measuring landscape erosion, we build on the work of Wisniewski and Pazzaglia (2002), who observed the presence of a broad (100 km) convexity in the profile of the Canadian River and its fluvial terraces. They inferred that the cause for this profile morphology was epeirogenic arching related to dynamic support of the lithosphere by warm, buoyant asthenosphere. Since the effects of epeirogeny should occur at a regional scale over long wavelengths, their interpretation raises the question of whether similar convexities can be seen in the profiles of other regional streams. Following the methods of Whipple et al. (2007), we used channel steepness indices to quantitatively compare the longitudinal stream profiles of the Canadian River with other drainages in our study area. The goal of this analysis was to test the hypothesis that profiles of rivers crossing the Jemez lineament may exhibit oversteepened, disequilibrium profiles that could reflect some combination of differential surface uplift, variable bedrock strength, and/or transient knickpoints in the river profile.

Channel steepness analysis is based on a power-law scaling between local channel slope (S) and upstream catchment area (A):
where ks and θ are indices of steepness and concavity, respectively. This relationship between channel slope and catchment area has been observed in a variety of tectonic and climatic settings (Flint, 1974; Hack, 1973; Howard and Kerby, 1983). In order to make direct comparisons of channel slopes among different rivers, a reference concavity θref is widely used to derive a dimensionless normalized channel steepness value, ksn (Wobus et al., 2006). We used the common reference concavity of θref = 0.45 in our calculations (Burbank and Anderson, 2011) and 90-m-resolution digital elevation models (DEMs) from the National Elevation Data Set. Additional information regarding user-defined parameters in the digital analysis of DEMs is recorded in Table DR3 in the Supplemental File (see footnote 1).

The results of this analysis (Fig. 12A) indicate that the normalized steepness (ksn) of the Canadian River is low (<60) in the Southern Rocky Mountains (Sangre de Cristo Range), but it steps up to moderate values (60–80) as it approaches the Cimarron escarpment at Raton. After the river passes Raton, ksn returns to low values until Springer, where it steps up to moderate-to-high (>120) values that persist throughout the Canadian River canyon. This steep canyon reach approximately coincides with the bowed reach inferred from the terrace study of Wisniewski and Pazzaglia (2002) (Fig. 12C). Finally, ksn returns to low values downstream of the Canadian escarpment, punctuated only by artificial steps associated with Conchas and Ute Lakes. Other streams in the study area exhibit low normalized steepness (ksn) everywhere, except in the Southern Rocky Mountains and along a broad, NE-trending region associated with the Canadian escarpment, highlighted with a white ellipse in Figure 12A. The pattern of convexities is generally consistent with a base-level fall or headwater uplift signal being recorded primarily in the larger drainages, but steep reaches in the headwaters of the Raton-Clayton and Ocate volcanic fields, in areas well away from trunk streams, may be more consistent with headwater uplift.

Normalized channel steepness is expected to vary with uplift rate and substrate erodibility (Duvall et al., 2004; Kirby and Whipple, 2003; Lague, 2003; Snyder et al., 2000; Whipple and Tucker, 1999), and recent studies have documented systematic increases in ksn as erosion rates increase in landscapes exhibiting a relatively uniform precipitation and lithology (DiBiase et al., 2010; Ouimet et al., 2009). In these studies, channels are argued to be locally adjusted to the differential trunk river incision (Ouimet et al., 2009) and differential rock uplift (DiBiase et al., 2010). Although lithology could play a major role in the steepness of channels in some drainages in our study area (e.g., the Maxson Crater basalts occupy the Mora River canyon for tens of kilometers), there does not seem to be any lithological marker that is particularly resistant to erosion that could be consistently responsible for the high-ksn reaches in other drainages. The fact that these knickpoints are evident in the modern river profiles, as well as in terrace geometries dating back as far as ca. 1.6 Ma (Wisniewski and Pazzaglia, 2002), suggests that some larger-scale control may be causing the regional convexity and its persistence through time in the form of bowed terraces. Figure 3A shows that the Dakota Sandstone holds up the Canadian River escarpment and is generally in the same map position as the high-ksn reaches (Fig. 12B). However, in detail (Fig. 12C), the steep reaches do not closely correspond to the Dakota lithology and, instead, steep reaches are often present in soft Jurassic rocks. Over the 1.57 Ma time scale of the bowed terraces, the river would have re-established itself in equilibrium with the softer rock substrate. Thus, similar to Wisniewski and Pazzaglia (2002), our preference is to consider other possibilities to explain the combined data sets. The geometry of the NE-trending oval containing steep river reaches (Fig. 12A) is located above the NE-trending Jemez mantle-velocity gradient (Fig. 1B) and suggests the possibility of dynamic, differential uplift to help explain the position of these knickpoints (e.g., Karlstrom et al., 2012).


Other NE-trending landscape elements in this part of the Great Plains include two regionally prominent erosional escarpments, the Canadian and Cimarron escarpments (Fig. 2A). The origins and ages of these can be inferred, in part, using the positions of dated volcanic flows and other stratigraphic relationships. The Cimarron escarpment extends for 80 km between Cimarron and Raton, New Mexico, forming the southeastern edge of the Park Plateau (Fig. 2A). The Park Plateau is a surface (capped by gravels in places) beveled into early Paleogene rocks of the Raton Basin, a Laramide sedimentary depocenter that began filling during the retreat of the Late Cretaceous epeiric seaway and continued to fill until the middle Paleocene (Cather, 2004). Since the Cimarron escarpment roughly mimics the preserved limits of the Raton Basin (Cather, 2004), its present position and orientation could be related to the architecture of the Raton Basin itself. However, there is evidence instead that the Cimarron escarpment is a relatively young feature and could not have existed before the middle Pliocene, based on the geometry of ca. 3.4 Ma Bartlett Mesa flows located near its northeastern terminus. These flows erupted onto a low-relief surface with no indication of cascades into steep canyons or down escarpment faces. Moreover, ridge-top elevations in the Park Plateau correlate with the elevations of 4–5 Ma basalt-capped surfaces in the Ocate volcanic field (Fig. 2A). This suggests that prior to ca. 4.5 Ma, a continuous low-relief surface may have extended between the ridge tops of the modern Park Plateau and the northeastern parts of the Ocate volcanic field, the Urraca surface of O’Neill (1988), which now has convexities in it due to post–4–5 Ma deformation (O’Neill, 1988).

The more southerly Canadian escarpment was described by Wisniewski and Pazzaglia (2002) as a relatively youthful landform based on its high topography and lack of large consequent drainages. More specific constraints on the age of the Canadian escarpment are as follows: (1) The escarpment could not have existed before ca. 4.5 Ma, the estimated end of Ogallala Group deposition on the southern High Plains (Chapin, 2008), because preserved remnants of Ogallala Group strata rest at similar elevations on the Canadian escarpment, Mesa Rica (to the south), and on the Caprock escarpment (even further south). Even though the Ogallala Group is likely a series of inset cut-and-fill deposits, lacking a single common strath, the regional base of the Ogallala section is subplanar on the scale of the study area, and there was likely a semicontinuous, gently sloping, low-relief Miocene paleo–Great Plains surface that extended across what is now the Canadian River valley (Fig. 13; Cather et al., 2012). (2) Circa 4.7 Ma Yates flows outcrop at the very edge of the modern Canadian escarpment, yet there are no cascades to suggest that they poured down a paleo-escarpment face. (3) The 1.57 Ma Maxson Crater basalts (Olmsted and McIntosh, 2004) flowed down the paleo–Mora River and Canadian River canyons, which were approximately half as deep as their modern counterparts. Since depth of the canyon and the height of the escarpment are equal at the canyon mouth near Sabinoso (Fig. 2A), this could indicate that the escarpment was also only half as tall at the time of emplacement of the Maxson Crater flows, as it is now. By extrapolating a steady rate of Canadian River canyon incision between 1.57 Ma and the present, we can estimate that canyon carving had begun by at least ca. 3 Ma. This time might also be a good estimate for the establishment of erosional escarpments. Finally, (4) the channelized Maxson Crater flows extend ∼10 km SE of the modern escarpment/canyon mouth (basalts near Sabinoso, New Mexico; Fig. 2A), but they appear to have been channelized at the time of emplacement. This suggests ∼10 km of lateral escarpment retreat since 1.57 Ma, at a rate of ∼6 km/Ma.


River geometries provide a sensitive record of the potential interplay between tectonics (e.g., differential rock uplift) and geomorphic adjustments like headward erosion and stream capture. Stream piracy, for example, can be a cause for differential incision and/or a symptom of tectonic movements, and hence the evolution of drainages in this region is a key data set to test alternate models. We combined new 40Ar/39Ar ages from the Raton-Clayton and Ocate volcanic fields with preexisting drainage evolution scenarios from other workers (Seni, 1980; O’Neill, 1988; Wisniewski and Pazzaglia, 2002; Galloway et al., 2011; Cather et al., 2012) to provide new constraints on the evolution of the stream network since the late Miocene. Figure 14 shows our proposed drainage reconstruction at four time intervals: the late Miocene, early and middle Pliocene, and Pleistocene. During the late Miocene at ca. 6 Ma (Fig. 14A), aggrading streams still flowed E-SW off the Southern Rocky Mountains (O’Neill, 1988) toward the distal fans of the Ogallala distributive fluvial systems near the Texas–New Mexico border (Seni, 1980). The nascent Raton-Clayton and Ocate volcanic fields consisted of only a few eruptive centers located primarily in the headwater regions of major streams. Ogallala Group river gravels below the 5.1 Ma Mesa de Maya and Black Mesa outflow (locations shown in Fig. 4) provide strong evidence that the headwaters of the ancestral Dry Cimarron River were located in the Southern Rocky Mountains and flowed eastward during this time (Stroud, 1997; this study). Likewise, what is here termed the ancestral Vermejo River flowed E-SE and supplied the Dalhart-Amarillo sedimentary fan lobe in the west Texas Panhandle (Seni, 1980). Several basalt-capped gravel surfaces near Angel Fire slope toward the southeast (O’Neill, 1988), suggesting that the ancestral Mora and Pecos Rivers took a more southeasterly route off of the mountains before continuing on toward the alluvial fans to the east. In the Miocene, there would have been no N-S integration of either the Canadian or Pecos River systems, as evidenced by the presence of elevated Ogallala Group gravels far to the east of these N-S drainages that contain Precambrian clasts originating from the Rocky Mountains. The dominantly E-flowing arrangement of drainages during the Miocene mimics the geometry of the N-S–oriented uplift of the Rocky Mountains and reflects Miocene and earlier regional uplift components.

By the early Pliocene (Fig. 14B), widespread Ogallala Group deposition had largely ceased (not later than 4.5 Ma; Chapin, 2008), and streams had begun to incise the Great Plains surface (the “turnaround” of McMillan et al., 2006). In some cases, new canyons/drainages were cut south of, and parallel to, basalt-filled paleovalleys/depressions, potentially indicating a relative base-level fall to the southeast or headwater uplift in the northwest. In the north, for example, eruptions at Mesa de Maya filled the channel of the ancestral Dry Cimarron River at ca. 5.1 Ma with a long (46 km) basalt runout, instigating the stream to shift to the south edge of the basalt flow to continue incising. Elsewhere, eruptions continued to occur near the headwaters of the ancestral Mora River and also along the ancestral Vermejo River, as evidenced by the ca. 5.3–4.7 Ma Yates flows. We observed Ogallala gravels located stratigraphically above the 4.7 Ma Yates flows, indicating that the ancestral Vermejo River was still fed by its headwaters in the Rocky Mountains until after 4.7 Ma. Farther to the east, the base of the Ogallala Group had already experienced substantial subsidence (up to 100–200 m) in the Texas Panhandle due to the large-scale dissolution and collapse of Permian evaporite deposits, which may have begun in the middle Miocene (Gustavson, 1986). This extended collapse history is recorded in thickening sequences of Ogallala deposits, tilted lacustrine deposits, dated sink holes, and faulted fluvial terraces, and may continue to the present day based on the high solute load of the modern Canadian River (Gustavson, 1986).

Around the middle Pliocene (3.6 Ma; Fig. 14C), the ancestral Canadian and Dry Cimarron Rivers had probably begun to carve their deep canyons, and the Canadian River system headwaters had begun to migrate northward, likely by headward erosion from its confluence with the ancestral Mora, and capture the ancestral Vermejo River. This capture event is inferred by the onset of rapid landscape denudation after ca. 4.5 Ma in the Charette Lakes area, which could have been driven by increased incision as tributary streams were integrated with the master Canadian River system. New volcanic eruptions continued to occur throughout the region around this time.

By the Pleistocene (2.2 Ma; Fig. 14D), canyons and escarpments were probably well developed (O’Neill, 1988; Wisniewski and Pazzaglia, 2002) but had not likely attained more than half of their modern depth/height, assuming steady-state incision along the Canadian River (Wisniewski and Pazzaglia, 2002; described in the following section). The migration of the Canadian River headwaters had likely continued northward, capturing the headwaters of the ancestral Dry Cimarron River and cutting it off from a large portion of its catchment area, just as many smaller streams draining the southeastern Raton-Clayton field had previously been cut off from their mountain headwaters as the Canadian River system became integrated. The Dry Cimarron River capture event is inferred from the order-of-magnitude increase in landscape denudation rates near Raton, New Mexico, at ca. 3.4 Ma (Fig. 11). By the Pleistocene, stream networks largely resembled their modern configurations, with the S-flowing Canadian River being well established as the master stream draining the eastern flank of the southernmost Rocky Mountains. Similarly, the ancestral, E-flowing Pecos River had been diverted toward the south sometime after 3.6 Ma (Wisniewski and Pazzaglia, 2002).

The primary question arising from this set of observations is: What combination of tectonic and geomorphic forcing factors may have driven drainage reorganization from a series of ESE-flowing streams draining the Southern Rocky Mountains during the late Miocene to a system dominated by the SSE-flowing Canadian River? One option is base-level fall associated with regional evaporite collapse downstream of our study area. Large-scale dissolution of Permian bedded salts along the modern Canadian River in west Texas has been documented by Gustavson (1986), who argued that dissolution has led to significant subsidence (100–200 m) of the base of the Ogallala Group since the middle Miocene and may have exerted some control over the location and morphology of the modern Canadian River valley. It is clear that a relative base-level drop of this magnitude has the strong potential to influence the rates and patterns of erosion observed upstream in the Raton section and must be considered carefully.

There are two lines of reasoning, however, which likely preclude evaporite collapse in west Texas from exerting a dominant control of the patterns of erosion observed in this study. First, the magnitude of subsidence due to evaporite dissolution is relatively small. It seems implausible that a 200 m drop in base level could result in the carving of canyons and escarpments that are twice as deep. Moreover, this subsidence does not explain the spatial patterns of long-term average denudation that increase toward the headwaters from the NE-trending, low- to no-erosion hinge line (Fig. 10A).

Second, the timing of evaporite collapse does not closely correspond with the timing of observed erosion upstream. Gustavson (1986) noted that a significant portion of dissolution-induced subsidence may have occurred since the middle Miocene—during the deposition of Ogallala Group and before the turnaround from net aggradation to net erosion—as evidenced by thickening Ogallala sedimentary packages into the zone of proposed subsidence. In other words, much of the documented dissolution and subsidence may have occurred during a period when E-flowing rivers were aggrading, long before the development of the Canadian River canyon and escarpment, and it therefore seems unlikely that downstream evaporite collapse can be implicated as a dominant driver of spatial and temporal erosion patterns in this part of the Great Plains.

A second option involves headwater uplift due to broad differential rock uplift along the NE-trending Jemez lineament. A NE-trending uplift that overprinted the earlier N-S Rocky Mountain uplift could have influenced the change from E- to SE-flowing rivers. Evidence potentially supporting this interpretation is shown in Figure 13 (modified from Cather et al., 2012), which shows the reconstructed 4.5 Ma base of the Ogallala Group in N-S cross section drawn at 103.5 W, parallel to the Rocky Mountain front. This basal contact appears to be arched above the Jemez low-velocity mantle as it bevels across warped Mesozoic substrate. Given that the Ogallala aggradational surface developed during the time of E-flowing rivers, this warping of that surface, plus the drainage reorganization, could be due to tectonic uplift along the NE-trending Jemez lineament (Fig. 14D). Likewise, the NE-trending zone of highksn knickpoints (described earlier) in this model may be explained as transient knickpoints due to adjustment of rivers to steeper headwaters as the escarpments also migrated up the tectonically modified topographic slope. In this model, the dynamic support and uplift above the Jemez lineament set into motion a series of base-level falls, drainage rearrangement, and stream captures that have been incising the western High Plains in an unsteady and nonuniform way since the end of the Miocene.


In the northern part of the Raton-Clayton volcanic field, near Raton, New Mexico, the apparent onset of relatively rapid denudation at ca. 3.4 Ma (Fig. 11) was preceded by a period of low (∼10 m/Ma) surface denudation rates beginning no later than 9 Ma. This period coincides, at least partly, with the deposition of Miocene Ogallala Group sediments (until at least 4.5 Ma; Chapin, 2008) and suggests a long-lived period of relative surface stability in the northern portion of the study area between 9 and 3.4 Ma. Thus, by examining only the basalt flows in the northern part of the study area, which were emplaced during this period of relative surface stability, these flows can be considered as capping a mappable, composite paleosurface that persisted from ca. 9 Ma to 3–4 Ma. From west to east, these flows include: Bartlett Mesa (3.37 ± 0.10 Ma), Barela Mesa (9.04 ± 0.04 Ma), Johnson Mesa (7.73 ± 0.03 Ma), Mesa de Maya (5.05 ± 0.02 Ma), and Black Mesa (5.10 ± 0.03 Ma), which is a long, sinuous outflow from the Mesa de Maya basalts (for mesa locations, see Fig. 4). Moreover, because most of these flows (with the possible exception of Bartlett Mesa) are underlain by thin Ogallala river gravels containing Precambrian clasts originating from the Rocky Mountains, this composite surface may approximately mimic the Ogallala transport slope in the northern part of our study area. This composite 9–3.4 Ma surface in the Raton-Clayton field is best visualized if each flow is projected into a W-E profile along the transect C–C′ (Fig. 10A).

Figure 15B shows a topographic profile along C–C′ starting in the Rocky Mountains and decreasing in elevation onto the High Plains to the east. The projections of Miocene–Pliocene basalts (red) are also shown, which were erupted during the period of relative landscape stability and define a gently concave-up profile. This profile morphology is characteristic of fluvially carved valleys and/or paleo–Great Plains pediment surfaces and is possibly continuous with the one of the Ocate surfaces defined by O’Neill (1988). The average slope of this paleosurface is ∼8.5 m/km, i.e., nearly twice as steep as the average slope of the modern Dry Cimarron River (∼4.6 m/km, blue), which runs parallel to the line of transect. There are two options that may explain the separation of the blue profile and the red profile: (1) The older (red) profile was originally carved at the gradient of the younger (blue) profile, and it has since been rotated to a steeper gradient by uplift concentrated to the west; or (2) the older (red) profile was carved at the gradient shown, followed by an increase in concavity driven by climate change and changes hydrologic parameters (e.g., Wobus et al., 2010; Zaprowski, 2005) that produced lower modern gradients, as shown by the blue profile.

Option 1 was supported by sedimentological studies in the Ogallala Group in the northern Great Plains (McMillan et al., 2002; Duller et al., 2012). If we follow this model and assume that the ca. 3.4 Ma Ogallala depositional surface formed at the same average slope as that of the modern Dry Cimarron, we calculate that this surface may have tilted 64 ± 7 millidegrees/Ma (millidegree = 10–3 degrees) since ca. 3.4 Ma. Starting at our interpreted tilt axis in Figure 15B, this seemingly insignificant angular tilt rate becomes nontrivial toward the western extent of our age control, where it translates into an apparent differential uplift of 650 m since 3.4 Ma. Data used in the calculation of tilt rates are provided in Table DR4 in the Supplemental File (see footnote 1).

The assumption regarding the original gradient of the Ogallala depositional surface is important because it directly influences our estimate for the magnitude of apparent tilting and differential uplift since ca. 3.4 Ma in the Raton-Clayton field. While the original depositional gradient of the Ogallala Group certainly could have been steeper than we suggest, there are two primary reasons for concluding that the original transport slope may have actually been shallower than that of the modern Dry Cimarron River, thus rendering our calculations as conservative estimates of apparent tilting. First, our assumed original Ogallala transport slope is steeper than estimates made by previous workers, which were determined using paleohydraulic methods to be in the range of 0.8–3 m/km at similar distances from the Rocky Mountain front on the High Plains of southern Wyoming and Nebraska (McMillan et al., 2002; Duller et al., 2012) and 1.7 m/km on the Colorado Piedmont (Leonard, 2002). Second, the modern Dry Cimarron River almost certainly has a much smaller catchment today than it did at ca. 3.4 Ma, when its headwaters were supplying Precambrian clasts from ∼100 km to the west in the Rocky Mountains. Since river channel slope is inversely related to upstream catchment area via Equation 1, it is likely that the average slope of the modern river is steeper now than it was in the Pliocene. Moreover, other obstacles persist when explaining channel slope differences via a climate-driven increase in concavity (option 2): The convex bowed river terraces in the Canadian River (Wisniewski and Pazzaglia, 2002), for example, require the major trunk stream to preserve a long-term disequilibrium profile; the convexities in sub-basalt, gravel-capped paleosurfaces in the Ocate volcanic field (O’Neill, 1988) indicate surface deformation since the end of the Miocene; and large-scale geometry of the sub–Ogallala Group unconformity shows broad warping across the Jemez lineament from north to south (Fig. 13).

Thus, option 1 (the surface tilting hypothesis) seems to be compatible with existing data, while option 2 (increased channel concavity hypothesis) does not. However, even if the apparent tilt described here represents a physical rotation of the Great Plains surface, there are still at least two plausible explanations for this rotation since ca. 3.4 Ma. The first is some form of tectonic tilting, which has traditionally been inferred to explain these types of geometries (e.g., Epis et al., 1980; Trimble, 1980; Eaton, 1986, 2008). An alternative model involves passive isostatic uplift in response to post-Miocene erosion that is driven by changes in climate (Molnar and England, 1990; Gregory and Chase, 1992, 1994; Wobus et al., 2010). The modeling results of Lazear et al. (2011), however, suggest that erosional flexural isostasy only accounts for ∼300 m of differential isostatic rebound since 10 Ma (15 millidegrees/Ma) along the C–C′ transect (Fig. 15A). This magnitude amounts to less than half of our estimated maximum differential uplift, leaving ∼350 m of uplift since ca. 3.4 Ma (34 millidegrees/Ma) that must be explained by some tectonic forcing. Given the long wavelength of apparent tilting (at least 400 km) and the position of this paleosurface above the Jemez low-velocity mantle anomaly (Fig. 15C), dynamic, mantle-driven uplift seems to be a likely cause for the observed differential uplift. Moreover, the interpreted tilt axis is coincident with the NE-trending hinge line of low to no erosion shown in Figure 10A.

Apparent late Cenozoic tilting has also been observed elsewhere on the western Great Plains, and Figure 16 synthesizes the results of studies that have previously proposed tilts of similar magnitudes. McMillan et al. (2002) measured tilting of Ogallala Group sediments far to the north of our study area in Wyoming and western Nebraska and inferred a maximum of 540 m of uplift/tilting of the Rockies since 18 Ma. Although angular rates were not explicitly calculated by McMillan et al. (2002), their estimates are equivalent to ∼9 millidegrees/Ma since 18 Ma. Leonard (2002) also made estimates of Ogallala Group tilting along transects spanning the Colorado Piedmont that yielded a maximum of 600 m/Ma of tilt since 17.5 Ma along the Canadian-Arkansas interfluve. Leonard’s estimated tilt is also equivalent to ∼9 millidegrees/Ma. Still other workers have noted a tilted middle Cenozoic partial annealing zone (PAZ), defined by apatite fission-track thermochronology of surface and borehole samples, that indicates a rate of tilting equivalent to ∼22 millidegrees/Ma since ca. 28 Ma (House et al., 2003; Leonard et al., 2002).

Notably, the angular tilt rate that we have calculated in our study area is ∼3–7 times greater than the rate-equivalent values calculated in previous studies elsewhere on the Great Plains. Furthermore, our proposed tilts occur across a NE-trending hinge line parallel to the Jemez lineament, rather than across the N-S Rocky Mountain trend, and they are considerably younger (averaged over 3.4 Ma compared to 17.5–28 Ma). If the amount of tilting estimated by McMillan et al. (2002) and Leonard (2002) took place since the end of the Miocene (5 Ma, similar to our study area), however, then their angular tilt rates would increase to 35 millidegrees/Ma, in closer agreement with our results. This may indicate that either more tilting has occurred within our study area, or that the rate of tilting of the entire western Great Plains has accelerated in late Cenozoic time, with much of the observed tilting occurring since the end of the Miocene.

Both McMillan et al. (2002) and Leonard (2002) have also argued against the model that climate-driven erosional isostasy can account for the full magnitudes of tilting based on results of flexural models, which replicate only about one half of observed tilts. Furthermore, these authors suggested various tectonic drivers to account for the remaining tilt, including the northward propagation of the Rio Grande rift or mantle-driven uplift. As described earlier herein, we suggest that the observed patterns of surface lowering and apparent tilting of paleosurfaces are consistent with post-Miocene mantle-driven uplift. However, because our study area is located in a special position above the Jemez low-velocity mantle anomaly, it is difficult to speculate about the influence of dynamic uplift on the study areas of McMillan et al. (2002) and Leonard (2002) farther to the north. Recent work by Duller et al. (2012) in Wyoming and Nebraska, however, does suggest that that short-lived dynamic uplift may be responsible for tilting in these study areas as well.


Another test of Wisniewski and Pazzaglia’s (2002) hypothesis that mantle-driven dynamic uplift may be affecting incision of the Canadian River canyon is to compare spatial patterns of long-term denudation with the geometry of mantle-velocity domains. Figures 17A and 17B were created by using GIS to query the digital tomographic model of Schmandt and Humphreys (2010) for the P-wave velocity at 80 km depth directly underneath the denudation point locations listed in Table 3. These plots of denudation rate versus underlying P-wave velocity show that the highest denudation rates calculated for the Raton section are concentrated above mantle with the lowest seismic velocities. These locations are in the NW part of the area (reference Fig. 1B), where possible dynamic components of uplift would be expected to be greatest. Southeast decreases in denudation rates and magnitudes toward the hinge line of low to no erosion (Fig. 10A) are accompanied by a SE increase in mantle velocity. Moreover, the erosion hinge line is subparallel to a sharp mantle-velocity boundary, and this observation is compatible with the notion that it may also be a tectonic hinge line above the edge of an uplifting region directly above a low-mantle-velocity domain.


Although erosional patterns are generally consistent with epeirogenic uplift above the Jemez lineament, climate changes likely interacted with tectonic forcing over the last 5–6 Ma, and it is important to evaluate whether climate forcing may have been a dominant control on landscape evolution. There are several postulated changes in post-Miocene climate that broadly correspond to changes in denudation rates that some authors have argued could provide a control on the turnaround from aggradation to incision of the western Great Plains. Chapin (2008) argued that opening of the Gulf of California by ca. 6.4 Ma resulted in an intensification of the North American Monsoon and increased erosion in the semiarid Southern Rocky Mountains–western Great Plains area (Mann and Meltzer, 2007; Tucker et al., 2006) since the latest Miocene. This proposed climate shift, however, predates the onset of both rapid denudation and drainage reorganization, thus making it difficult to tie landscape change in our study area to this climate change episode.

Wobus et al. (2010) listed several possible climatic shifts that may have contributed to post-Miocene landscape evolution of the Great Plains by driving increases in erosion: (1) The global δ18O record shows the beginnings of an ∼2‰ negative shift started near the Miocene-Pliocene boundary (5.3 Ma) (Zachos et al., 2004); this may have caused a shift toward more snowmelt-dominated hydrologic systems, which modified the type and extent of vegetative cover in the Rockies and/or changed atmospheric circulation patterns and the frequency distribution of storms. (2) The carbon record from the High Plains begins a regional ∼6‰ positive shift at the end of the Miocene (Fox and Koch, 2003), consistent with an increase in aridity to the east of the Rockies. (3) The terrestrial record shows the first major advance of ice sheets into the central United States around 2.5 Ma (e.g., Balco et al., 2005), which was most likely accompanied by changes in alpine glaciation in the headwater systems as well and could have modified water and sediment fluxes in streams draining the Southern Rocky Mountains.

While the advance of ice sheets at 2.5 Ma postdates the observed onset of relatively rapid denudation in both the Raton-Clayton and Ocate volcanic fields, climate changes recorded by oxygen and carbon isotopic records at the end of the Miocene stand as possible mechanisms for increased erosion. However, the timing of these two mechanisms only loosely correlates with the onset of rapid denudation in the Ocate volcanic field and predates rapid denudation in the Raton-Clayton field (at 3.4 Ma) by as much as 1.9 Ma. Most importantly, such climate forcing would be expected to be regional and therefore seems insufficient to explain the spatial and temporal patterns of differential erosion observed across our relatively small study area. One of the largest obstacles to climate-driven landscape models, as discussed by Wobus et al. (2010), is that there is no direct evidence of the ways in which documented late Cenozoic climatic changes would have influenced the hydrology of the Great Plains river systems. That is, even if climate events were temporally well correlated with the turnaround from aggradation to erosion, there is still uncertainty about how those changes would impact the pattern of erosion on the plains. Our view is that the best way to distinguish between patterns of river incision that are largely the result of climatically driven geomorphic parameters versus those that would result from differential tectonic uplift is to examine temporal and spatial patterns of change, which need to correlate with the timing of postulated climate changes, have temporal scales of differential erosion that correspond with climate cycles (e.g., river terraces), and a have a spatial variation that can be explained by regional climate patterns.

Wobus et al. (2010) attempted to do this using a mathematical sediment transport model to identify distinctive patterns of landscape incision. Their model suggested that sustained rock uplift creates river channels that are everywhere steeper than their initial gradients, while incision driven by an increase in discharge or a decrease in sediment supply creates channels of lower steepness. Wobus et al. (2010) dismissed the paleohydraulic analysis by McMillan et al. (2002), which suggested channels have been steepened, citing uncertainty and high margins of error in their analysis, and argued that the gradients of the Arkansas and South Platte Rivers have decreased since the late Miocene, which would be consistent with their predicted effects for climate-driven incision. However, paleogradients of these rivers are unknown to the extent that these conclusions are equivocal. For the Canadian River system, our data support the arguments and approaches of McMillan et al. (2002), Leonard (2002), and Duller et al. (2012), which suggest that the paleo–Great Plains surfaces and river channels have experienced a period of steepening since the late Miocene. Additionally, our channel steepness analysis further supports this concept by showing that rivers and erosional escarpments have been actively responding to disequilibrium conditions via knickpoint formation rather than becoming better graded and decreasing in slope over this time period.

Spatial patterns of drainage reorganization are also difficult to explain if one implicates climate-driven geomorphic parameters as the dominant driver of post-Miocene erosion on this part of the Great Plains, with little or no tectonic component. The reorganization of E-flowing rivers during the deposition of Ogallala Group sediments to the SE-flowing Canadian River system can be explained by base-level fall or headwater uplift that generated NW-SE topographic gradients. This geometry is readily explained by post-Miocene epeirogeny above the NE-trending Jemez lineament, but it seems less well explained by evaporite collapse, other base-level fall models, or climate-influence landscape modification.


Our preferred model for the onset of post-Miocene erosion involves magmatic inflation of the Raton-Clayton and Ocate volcanic fields in response to longer-wavelength mantle driven dynamic uplift. Our evaluation of surface evidence demonstrates a strong connection between mantle-velocity structure and a wide range of landscape features. These NE-trending landscape features include: (1) crude volcanic “belts” of similar age, (2) retreating erosional escarpments, (3) differential denudation measured from a NE-trending hinge line of low-to-no erosion, (4) a NE-trending zone of broad (50–100 km) convexities in stream profiles identified by normalized channel steepness index (ksn) analyses, (5) a reorganization of stream networks from E-flowing streams during the Miocene to a modern, SE-flowing Canadian River system that takes advantage of a relative base-level fall in the southeast, and (6) an 40Ar/39Ar-dated composite paleosurface that indicates total tilting of 64 millidegrees/Ma (with 34 millidegrees/Ma of tectonic tilting) since ca. 3.4 Ma. These local trends are overprinted on and stand in contrast to the overall N-S-trending orientation of features associated with the Southern Rocky Mountains and the broader “Alvarado Ridge” as defined by Eaton (1986) and redefined by Eaton (2008).

We argue that the crucial data connecting landscape response with a dynamic tectonic forcing are the 40Ar/39Ar ages on volcanic flows, which point to large-volume basaltic magmatism (by proxy of eruption surface area) above the Jemez lineament from 3.8 to 1.45 Ma. This magmatism suggests substantial melt flux into and through the crust over a several-million-year time frame that closely coincides with the onset of rapid landscape denudation in these volcanic fields. Our suggestion that the onset of voluminous magmatism (and perhaps onset of denudation) was a bit younger in the Raton-Clayton field than the Ocate field raises the possibility of NE propagation of the Jemez lineament magmatism and associated uplift. This needs further investigation, but at the scale of this study, the spatial and temporal association of mantle, magmatic, and landscape features strongly suggests that mantle-driven surface uplift has shaped the Great Plains landscape over the last 5 Ma.

Neogene climate change events likely also played important roles in late Cenozoic landscape evolution, but apparently at different temporal and spatial scales relative to the landscape features we observe. For example, postulated climate changes summarized by Chapin (2008) and Wobus et al. (2010) do not aptly explain patterns, magnitudes, or timing of long-term differential denudation and tilting. Moreover, climate-driven erosional isostasy does not fully explain the timing or magnitudes tilting, as the flexural models of McMillan et al. (2002), Leonard (2002), and Lazear et al. (2011) show that this effect accounts for only ∼50% of apparent tilts, thus necessitating a tectonic component of uplift to explain the remainder. Evaporite collapse mechanisms for base-level fall likely operated in the past (e.g., Gustavson, 1986), but their magnitude and timing in west Texas do not account for the observed differential denudation relative to a NE-trending hinge line or the onset of more rapid denudation rates since ca. 5 Ma. Instead, we envision salt collapse as a process that is entwined with tectonically-instigated landscape evolution. Similarly, the change from E-flowing to SE-flowing drainages is spatially consistent with a hypothesis for uplift along the NE-trending Jemez lineament, but it is unexplained by climate change scenarios.

Mechanisms for increasing mantle buoyancy remain poorly understood, but they seem to be related to the geologically-recent flux of basaltic magma along the Jemez lineament. Although asthenospheric flow pressures at the base of the lithosphere are modeled to cause hundreds of meters of surface uplift over many hundreds of kilometers (e.g., Moucha et al., 2009; van Wijk et al., 2010), they are not strongly favored here because of the relatively short wavelength of the Jemez zone (a few hundred kilometers). Uplift caused by magmatic inflation in dike and sill networks in the crust and mantle lithosphere below the Raton-Clayton and Ocate volcanic fields is perhaps a more likely explanation for the observed surface patterns, considering the scale of the Jemez zone. Thus, our preferred interpretation for the tilted paleosurface and patterns of landscape evolution in our western Great Plains study area involves erosion in response to mantle-driven dynamic topography, generated by an interrelated combination of flow (minor) and buoyancy change (major) in the crust and mantle lithosphere since the Miocene due to magma flux from the deeper mantle.

This research was generously funded by the National Science Foundation (NSF) via the New Mexico Alliance for Minority Participation’s “Bridge to the Doctorate” program (EHR-1026412, BDVIII, for A. Nereson via L. Crossey). Funding for K. Karlstrom and M. Heizler came from the NSF Continental Dynamics Program CREST experiment (EAR-0607808). Brandon Schmandt provided tomographic images used in this study. We extend our thanks to the commercial ranches and families for field access and their generous hospitality. We thank Eric Leonard, Frank Pazzaglia, and Mike Williams for thorough, insightful, and challenging reviews that greatly improved this manuscript.

1Supplemental File. PDF file of Tables DR1–DR4 and Figures DR1–DR72. If you are viewing the PDF of this paper or reading it offline, please visit http://dx.doi.org/10.1130/GES00837.S1 or the full-text article on www.gsapubs.org to view the Supplemental File.