Abstract
Underpressures (subhydrostatic heads) in the Paleozoic units underlying the Great Plains of North America are a consequence of Cenozoic uplift of the area. Based on tectonostratigraphic data, we have developed a cumulative uplift history with superimposed periods of deposition and erosion for the Great Plains for the period from 40 Ma to the present. Uplift, deposition, and erosion on an 800 km geologic cross-section extending from northeast Colorado to eastern Kansas is represented in nine time-stepped geohydrologic models. Sequential solution of the two-dimensional diffusion equation reveals the evolution of hydraulic head and underpressure in a changing structural environment after 40 Ma, culminating in an approximate match with the measured present-day values. The modeled and measured hydraulic head values indicate that underpressures increase to the west. The 2 to 0 Ma model indicates that the present-day hydraulic head values of the Paleozoic units have not reached steady state. This result is significant because it indicates that present-day hydraulic heads are not at equilibrium, and underpressures will increase in the future. The pattern uncovered by the series of nine MODFLOW models is of increased underpressures with time. Overall, the models indicate that tectonic uplift explains the development of underpressures in the Great Plains.
1. Introduction
The existence of underpressures (subhydrostatic hydraulic head values, Appendix A), in the geologic units that underlie the Great Plains [1, 2] of North America has been well documented. Appendix A presents the theoretical basis for defining and quantifying underpressure.
These underpressures are subsurface zones where the hydraulic head of confined aquifers is less than the elevation of the water table (Figure 1 and Appendix A). Belitz and Bredehoeft [3] provide a summary of underpressures identified in the Texas Panhandle [4], the San Juan Basin in New Mexico [5], the Denver Basin in Nebraska [6], and numerous other examples. Nelson et al. [7] used water levels measured in wells and oil and gas drillstem tests to quantify the underpressures in the Paleozoic units along a cross-section from eastern Colorado to eastern Kansas (Figure 2). The data and analysis identified large underpressures in the west, equivalent to several hundred meters of water. As noted by Nelson et al. [7] and Belitz and Bredehoeft [3], underpressure values are greatest in the west and are least where the geologic units crop out in the east.
A hydrodynamic model is deemed appropriate for the confined aquifers of the Great Plains. For example, Belitz and Bredehoeft [3], p. 1356, postulate that
Generally, subnormal fluid pressures might be found in any subaerial, topographically tilted, structural basin capped by a thick sequence of low-permeability rocks (i.e., shale or evaporites). The tilt can provide the topographic driving force for the fluid flow. The low-permeability cap can provide insulation from the elevation head of the water table, and the structure can provide the mechanism for reducing permeability in the basin deep that allows for better hydrologic connection to low-elevation outcrops than to high-elevation outcrops.
They present a hydrodynamic model of the Denver Basin and conclude that the underpressures are due to groundwater flow combined with limited recharge. Sorenson [8] concluded that the underpressures in the Panhandle-Hugoton gas field of Kansas, Oklahoma, and Texas were due to fluid flow and controlled by the reservoir outcrop elevations in eastern Kansas. Analysis of groundwater flow in the Dakota Sandstone aquifer of Kansas and southeastern Colorado concluded that the underpressure values in the Dakota Aquifer were due to groundwater discharge at the unit outcrop in eastern Kansas [9]. Nelson et al. [7] also concluded that the control of the hydraulic head values in Paleozoic strata in Colorado and Kansas was also exerted by the outcrops in eastern Kansas and Oklahoma.
A review of the literature finds that the work to date has provided some insight on our conceptual and numerical understanding of the existence and development of underpressure in the Great Plains. However, no one has attempted to incorporate the structural changes that may have provided the driving force. The key to a numerical groundwater flow model that would be able to represent the development of underpressure over geologic time in the Great Plains is, in our view, the inclusion of the following: the history of uplift and tilting, which is recorded in the tectonostratigraphy of the Great Plains, and the hydrologic framework of the Great Plains. Over the last 70 million years, the area that is now the Great Plains has undergone the natural geologic processes of uplift, tilting, subsidence, erosion, and deposition. Major structural changes include the Laramide orogeny and formation of the Rocky Mountains, subsidence of the Denver Basin, uplift of the Colorado Rockies, and erosion and exposure of the Paleozoic units in eastern Kansas. The area of the Great Plains has transitioned from the bottom of a large Cretaceous inland sea, the Cretaceous Western Interior Seaway, to a sea of grass elevated to over a mile above sea level at the base of the Rocky Mountains.
We hypothesize that underpressure is a condition driven by structural changes (primarily differential uplift) that occurred during the Cenozoic Era, resulting in exposure of outcrops in the eastern Great Plains and tilt of the land surface. To investigate this hypothesis and incorporate the structural changes in the Great Plains over the last 40 million years, we use a numerical groundwater flow model called MODFLOW [10], to solve the two-dimensional diffusion equation (the “groundwater flow equation” of hydrologists), and find the variation of hydraulic head with distance , elevation , and time :where is hydraulic conductivity in m s−1 and is specific storage, with dimensions of m−1. These terms and others used in this paper are presented under the Definitions of Terms.
Because MODFLOW utilizes a fixed geometry and does not incorporate time-varying structural changes, we elected to run a sequence of nine numerical models (), each representing a limited time interval, such as 27 to 18 Ma, with the geometry appropriate for that time interval. The structural changes required to transition from model to model , mainly elevation and tilt as calculated in a spreadsheet, are summed with the head values computed for model . These augmented values of then serve as starting values for model . This procedure is repeated until and a solution for is obtained at 0 Ma. Equation (1) is solved with uniform time steps of one year to ensure proper numerical solution of the transient groundwater flow, diffusion equation (see (1)), whereas the structural changes incorporated in the nine models are implemented over durations ranging from 1 to 9 million years depending on whether the uplift as a function of time was steep or shallow, to keep the total uplifts within all time intervals comparable.
2. The Geologic Setting
The sequence of nine structural models rests upon the history of uplift of the Great Plains, which is recorded in its tectonostratigraphy. Appendix B gives further details.
Eaton [11, 12] interpreted the Great Plains as the eastern flank of a huge antiformal uplift that is centered on the Rio Grande rift and includes the Colorado Plateau as its western flank (Figures 3, 4, and 5). The Rio Grande rift, in turn, is partly coincident with the southern Rocky Mountains, and the current high elevation of the Rocky Mountains dates only to uplift of this antiformal bulge, which entirely postdates the original formation of these mountains in the Laramide (~70 to ~40 Ma) orogeny. Eaton [11] called this giant domal uplift the Alvarado ridge. From shortly after the end of the Laramide orogeny (at ~40 Ma) to just before the beginning of formation of the Rio Grande rift (at ~27 Ma), subduction continued along the west coast of North America; however, contractional uplift of the southern Rocky Mountains had ceased. The onset of Rio Grande rifting to the southwest of the study area (near the central Colorado/New Mexico border), at ~27 Ma, was accompanied by an early stage of uplift [13].
The interval from ~18 to ~4.5 Ma on the Great Plains was characterized by the deposition of the Ogallala Formation, which is a vast apron of debris shed from the crest of the Alvarado ridge onto the Great Plains, and which is an important water-table aquifer in a major part of the Great Plains (Figure 5). Deposition of the Ogallala therefore reflects the uplift of the Alvarado ridge [14]. The Ogallala contains cobbles derived from the Rocky Mountains as far east as the Colorado-Kansas border, an indicator of its deposition on an already tilted and partially uplifted surface. Presumably, continued uplift and tilting during deposition of the Ogallala eventually created grades steep enough to end deposition and initiate erosion. Strong tilting of the Ogallala in the western part of the Great Plains indicates that at least about half of the uplift and tilting of the Great Plains occurred after 4.5 Ma. Current slopes across the Great Plains average ~10−3. This, combined with the sustained pattern of downcutting throughout the Great Plains since ~4.5 Ma, reflects a much stronger pattern of uplift since ~4.5 Ma than in the preceding ~27 to ~4.5 Ma. Whereas it is difficult if not impossible to quantify this, it appears likely that at least about half of the total uplift of the Great Plains occurred in the interval from ~4.5 Ma to the present.
3. Methodology
3.1. Understanding the Geology and Quantifying Structural Changes
A geologic map of the study area and a cross-section along line are presented in Figures 6 and 7, respectively. From the geologic evidence presented in Appendix B, we have estimated an uplift history for the Great Plains for the period of 27 Ma to the present plotted as cumulative uplift with time, along with superimposed periods of deposition and erosion (Figure 8). Uplift during the Laramide orogeny (70–40 Ma) was incorporated as 250 m of uplift at the Rocky Mountain Front and 50 m of uplift in the east, relative to the Cretaceous Western Interior Seaway. Uplift during the period of 40 to 27 Ma was assumed to be zero (Appendix B).
3.2. Assembling the Present-Day Geofluid Underpressure Data
Nelson et al. [7] present six maps of present-day hydraulic heads in the Paleozoic rocks that underlie the Great Plains. These values are derived from pressures measured in drillstem tests in oil and gas wells, combined with water levels measured in water wells. The maps extend northward from the Anadarko Basin in Oklahoma to southern Nebraska and eastward from eastern Colorado to eastern Kansas and Oklahoma. On all six maps, the potentiometric surfaces approach land elevation in the vicinity of the geologic unit outcrops in the eastern part of the mapped area. However, as land elevation rises to the west (Figure 2), the potentiometric surfaces drop below land surface, and underpressures develop. Figures 9, 10, and 11 present the potentiometric surfaces from three of the maps presented in Nelson et al. [7]. The figures indicate that the hydraulic gradient is smaller in central and eastern Kansas (wide contour spacing) and larger in western Kansas and eastern Colorado (narrow contour spacing). The smaller hydraulic gradient in central and eastern Kansas may be attributed to hydraulic conductivity values greater than those of western Kansas and eastern Colorado.
Figure 12 presents the hydraulic head values along cross-section (Figure 6) for the six potentiometric surfaces in Nelson et al. [7]. Figure 12 indicates that the hydraulic head data fall into two groups: the Wolfcampian, Virgilian, and Missourian profiles form one group with higher head values and higher gradients, and the Mississippian and Cambrian-Ordovician-Silurian profiles form a second group with lower head values and lower gradients. The Desmoinesian profile matches the low-head, low-gradient character in the east but displays high head values in the west. Replication of these features is a primary goal of the hydrologic model described in this paper.
3.3. Creating and Calibrating the Hydrologic Model to Simulate Structurally Driven Underpressures
The creation and calibration of the hydrologic model is best summarized in a flow chart as shown in Figure 13.
Each flow chart element # is referred to herein as FE # and is related to a section of the paper in Table 1.
Part one of the model creation and calibration process includes incorporating tectonic uplift, deposition, and erosion.
The sequence of MODFLOW models created for this work (represented conceptually by FE # (3) in the flow chart of Figure 13; and listed in Table 3) is based on the geologic cross-section running from northeastern Colorado to eastern Kansas by Nelson et al. [7]. The line of the cross-section is shown in Figures 2 and 6 and the cross-section in Figure 7. The cross-section cuts through, from top to bottom, the Cenozoic strata [Tertiary Laramide Basin Fill], the Mesozoic strata [Cretaceous Pierre and Niobrara Shales and Dakota Sandstone; Jurassic sediments], and the Paleozoic strata [Permian (Upper Permian, Lower Permian: Wolfcampian), Pennsylvanian (Virgilian, Missourian, Desmoinesian, Basal Pennsylvanian), Mississippian, Silurian, Ordovician, and Cambrian formations]. The geologic time line starts 40 Ma, the end of the Laramide orogeny and the beginning of a quiescent period lasting until the onset of uplift of the Great Plains at 27 Ma. The results from model 1 represent an equilibrium potentiometric regime for the geometry existing at 27 Ma. Wellbore information from 27 oil and gas wells (N. J. Gianoutsos, USGS, written commun., 2014), which is background to data published in Nelson et al. [7], was used to define the tops and bottoms of geologic formations (Figure 7). These data were also used to calculate formation thicknesses in each borehole using a Microsoft Excel spreadsheet. The MODFLOW preprocessor program ModelMuse [15] was used to create the model geometry. Figure 14 presents the first in the sequence of 9 cross-sectional models, representing the period from 40 to 27 Ma. The MODFLOW model has 10 “layers” spanning the vertical direction (Table 2), 26 “columns” spanning the horizontal direction, and 1 “row” perpendicular to the page (because we are representing three-dimensional reality by a 2-dimensional section in the plane of a report page). Each layer-column-row element of the model is referred to as a “cell.” The thicknesses of the formations in each borehole, as calculated above, define the elevations of the tops and bottoms of these formations/layers in Figure 14. A stratigraphic column (Table 2) presents the ages, lithologies, hydrogeologic unit, and thicknesses of these 10 model layers (see Appendix B for background). Even though the Ogallala, Laramide Basin Fill, and Dakota are shown in Table 2 to be aquifers, they were treated as aquitards, with a of 1 × 10−6 m/d, essentially to remove them from the analysis and emphasize the flow regimes through the deeper units where head data are available. The Ogallala, Laramide Basin Fill, and Dakota are only included for the sake of geologic completeness and for the option to expand the modeling to the upper horizons if desired in the future. The sensitivity analysis Table 7 (to be introduced later under Sensitivity Analysis) shows that if in the Ogallala, Laramide Basin Fill, and Dakota were to be increased by 2 orders of magnitude, the conceptual model of isolation of the deeper units would be compromised and in the Wolfcampian would change by −8.95 percent. Figures 14, 15, and 16 present the model grids for the periods ending at 27, 18, and 4.5 Ma, respectively.
To investigate the hypothesis that underpressure is a condition driven by structural changes that occurred during the Cenozoic Era, such as exposure of outcrops in the eastern Great Plains and tilt of the land surface, the sequence of MODFLOW models was developed (FE #’s of the flow chart in Figure 13; Table 3). This sequence of models represents our hydrologic modeling of the specific geologic events of uplift, deposition, and erosion, shown in Figure 8. The goal of the modeling was that the calculated heads for the ninth and last model of the sequence approximate the actual present-day hydraulic head data values observed in the field.
This model sequence (Table 3) consists of 9 multi-million-year, coarse-grid cross-sectional models with 10 layers that simulate the period of uplift, deposition, and erosion from 27 Ma to the present, represented in Figure 8. The structural deformation ends with a net uplift of 1.6 km in the west along the Rocky Mountain Front relative to eastern Kansas. Models are named for the end of a time period; so the 27 Ma model simulates the 13-million-year period from 40 to 27 Ma. The progression of uplift, deposition, and erosion is portrayed in four cross-sections (Figures 14, 15, 16, and 7), starting with the 27 Ma model and concluding with the 0 Ma model, which also serves as the present-day geologic cross-section (Figure 7). The primary effects through time are the increased tilting of the land surface, the increased tilting of all strata, and the higher elevation of the entire section. Elevation gain, tilt, and the development of a concave surface with time are also apparent on the stacked land surfaces of the nine models (Figure 17).
3.3.1. Uplift, Deposition, and Erosion
Uplift, deposition, and erosion were incorporated into each of the sequences of MODFLOW models. Table 4 presents the percent of total uplift (PTU) and cumulative percent of total uplift (CpTU) calculated for the 9 model intervals. As shown in Figure 8, post-Laramide uplift occurs from 27 Ma to the present, superimposed first by a period of deposition of the Ogallala, from 18 to 4.5 Ma, and then a period of erosion from 4.5 Ma to the present.
Figure 7 represents the present-day surface elevations. Figure 14 represents the land surface elevation at 27 Ma. The difference in elevation between these two surfaces represents the total uplift (TU), which varies with distance, , from the Golden Fault at the Rocky Mountain Front. The 27 Ma land surface is a tilted plane (seen as the straight line top of the two-dimensional cross-section of Figure 14), whereas the present-day land surface of Figure 7 is concave upwards. The upward concavity represents increased uplift as the Rocky Mountains are approached from the east.
The uplift process assumes the land surface at 27 Ma, shown in Figure 14, is a plane (straight line in two dimensions) rising from 50 m ASL in the east to 250 m ASL in the west. Because the uplift from 27 to 18 Ma is 13.125 percent of total uplift (PTU in Table 4), 13.125 percent of total uplift calculated above is added to the 27 Ma land surface of Figure 14 to obtain the 18 Ma land surface of Figure 15. This process is repeated for 18 to 11 Ma, 11 to 7 Ma, 7 to 4.5 Ma, and so on, to calculate the uplift at each stage of the modeling sequence, ending with the 0 Ma land surface shown in the geologic cross-section of Figure 7.
The model also includes periods of deposition. Deposition of the Ogallala takes place from 18 to 4.5 Ma (Figure 8). It was further assumed that the maximum deposition, at 4.5 Ma, was 183 m at km, near the midpoint and thickest location of the Ogallala Formation at 4.5 Ma, shown in Figure 16. We then assumed a linear spatial interpolation of Ogallala thickness from 183 m at km to 0 m at both east ( km) and west ( km) ends of the model, at 4.5 Ma. We then linearly interpolated with time from a thickness of 0 m at 18 Ma to the Ogallala thickness at 4.5 Ma to obtain the Ogallala deposition thickness at the intermediate times between 18 and 4.5 Ma, namely, 11 and 7 Ma. The Ogallala deposition was added to the uplifted surfaces of the 11, 7, and 4.5 Ma models to obtain the final surfaces of these models.
Uplift and deposition, as explained above, define the 4.5 Ma land surface shown in Figure 16. The model also includes erosion that occurred from 4.5 Ma to the present. The land surface at 4.5 Ma was eroded to the present-day (0 Ma) land surface shown in Figure 7. To calculate erosion, the 0 Ma land surface is subtracted from the 4.5 Ma land surface to obtain the erodible thickness (ET), which is a function of, that is, ET(x). This ET(x) has to be removed from the top of the model from 4.5 Ma to the present. So, we interpolate the ET(x) with time between the values calculated at 4.5 Ma and 0 m at 0 Ma. These interpolated values of ET(x) are subtracted from the uplifted surfaces of the 3.5, 3, 2, and 0 Ma models to obtain the final surfaces.
The second part of the calibration process includes running the sequence of models using variable rates of flux at the western boundary and variable specific storage (), horizontal hydraulic conductivity , and vertical hydraulic conductivity values, until the calculated present-day hydraulic head values approximated the present-day hydraulic head values. This is described in the next two sections.
3.3.2. Hydraulic Conductivity and Specific Storage
The simulations of the heads, horizontal gradients, and the relation among the hydraulic head profiles displayed in Figure 12 are highly sensitive to hydraulic conductivity , which can vary over many orders of magnitude. We have relied on previous modeling studies (Appendix C) to establish a range of reasonable values of .
For further adjustment of values, a useful tool is provided by Darcy’s equation relating flux to the hydraulic conductivity and the hydraulic gradient, , which aided us during calibration:where is the horizontal hydraulic conductivity and is hydraulic head. Equation (2) indicates that, for a specific , and are inversely proportional. That is, as the medium becomes more conductive and increases, the hydraulic gradient decreases and the plot of versus becomes shallower or less steep. Conversely, a steep hydraulic gradient indicates a tight or low hydraulic conductivity medium. Based on this, the and values were varied to make the calculated heads as close to actual values as possible.
The Mississippian geohydrologic unit represents the combined strata of Cambrian-Ordovician-Silurian as well as rocks of Mississippian age. The data (Figure 12) show that these lowermost units must be considerably more permeable (shallow hydraulic gradient) than the Virgilian-Missourian and Wolfcampian units (steeper hydraulic gradient) and, based on the separation of the hydraulic head values, must be hydraulically isolated. Hydraulic isolation is achieved by assigning low vertical permeability in the Desmoinesian and basal Pennsylvanian units. Permeability normal to bedding is generally considerably less than parallel to bedding. Consequently, we set in most of the geohydrologic units (Table 5), a selection that is the default setting in MODFLOW. However, was set to 2 × 10−10 m d−1 in the lowermost three layers to provide vertical isolation and to prevent leakage into underlying basement rocks. If the vertical hydraulic conductivity, , is large, the hydraulic head plots of versus for the 10 layers of the model tend to overlie each other. In order to articulate different behaviors in the vertically stacked 10 layers, vertical isolation, or low , is needed.
The values for the Mississippian geohydrologic unit range from 10−4 to 10−1 m d−1 within the upper part of the distribution shown in Figure 18, as described in Signor et al. [16]. After some experimentation, for the Mississippian geohydrologic unit was set to 5 × 10−3 m d−1, which falls within the range cited in Signor et al. [16]. Table 5 presents the final values from the model calibration.
Figure 12 shows an abrupt change in the hydraulic head of the Desmoinesian geohydrologic unit where the hydraulic head values transition from the shallower Wolfcampian-Virgilian-Missourian to the lower Mississippian and Cambrian-Ordovician-Silurian geohydrologic units. This abrupt change is due to a change in of the Desmoinesian geohydrologic unit. The model assigns a value of 1 × 10−9 m d−1 for the western section of the Desmoinesian and a value of 5 × 10−3 m d−1 for the eastern section (Table 5). This, combined with vertical isolation, enabled us to simulate how the Desmoinesian head plot resembles the Wolfcampian-Virgilian-Missourian head plots in the west and the Mississippian plot in the east in the present-day head data.
Two geohydrologic units—Virgilian-Missourian and Wolfcampian—exhibit small changes in slope of the hydraulic head at a distance of 400 km. To replicate these transitions in the Virgilian-Missourian and Wolfcampian, low-permeability segments referred to as chokes are inserted in the model. Within the chokes, is one-fifth that of on either side of the choke (Table 5). The low-permeability chokes may represent losses or discontinuities in large-scale hydraulic continuity resulting from the transition from mainly marine deposition in the east to exclusively terrestrial deposition in the west in upper Paleozoic rocks. The hydraulic conductivity of 5 × 10−4 m d−1 representing the permeable sections of the Virgilian-Missourian and Wolfcampian units is about one-tenth that of values found throughout Kansas in strata of Pennsylvanian-Jurassic age (Figure 18) and is one-twentieth of the value of 10−2 m d−1 used for shelf carbonates in the Palo Duro Basin (Figure 19).
of the upper five geohydrologic units was set to a low value of 1 × 10−6 m d−1. This choice is reasonable for the Upper Cretaceous Pierre Shale and Upper Permian-Jurassic units and does assure isolation of the lower units from land surface. However, it is unrealistic for the permeable Upper Cretaceous Dakota Formation, which was not modeled in this study. The broad horseshoe shape of the Dakota outcrop is difficult to model with a two-dimensional grid, and potentiometric maps indicate that the Dakota is not hydraulically connected to the Wolfcampian because of the presence of the low-permeability Upper Permian-Jurassic unit [7].
For the storage characteristics (compressibility of the medium and fluid), a specific storage, , value of 1 × 10−6 m−1 was used throughout the model. The storage coefficient, , was calculated as × layer thickness.
3.3.3. Boundary Conditions (BCs) for the Sequence of Models
The models’ western boundary is defined by the Rocky Mountain Front. It is modeled as a specified-flux boundary to represent flow from the west across the Golden Fault. By altering the flux until the head gradient matched the data (Darcy’s equation (2) shows that there is a direct relationship between flux and hydraulic gradient), it was found that 0.0003 m3 d−1 of water (109,500 m3 per 1 m.y.) enters each of the Wolfcampian (cross-sectional area 15,239 m2) and Desmoinesian (cross-sectional area 18,287 m2) and 0.0007 m3 d−1 (255,500 m3 per 1 m.y.) of water enters the Virgilian-Missourian (cross-sectional area 33,526 m2). These values were applied throughout the modeling sequence.
It should be noted that when the amount of specified flow rate into the western boundary, , or flux, , is increased, the hydraulic head in the geohydrologic unit receiving the flux increases and the plot of head versus becomes steeper in the west (based on Darcy’s equation (2)). The optimal values of flux () obtained by the calibration process were used as the western BCs for the sequence of models.
The models’ eastern boundary lies in east-central Kansas (Figures 2 and 6), where the confined units crop out and all layers are exposed to the atmosphere. From the definition of the hydraulic head, , in (A.1), we see that when the pressure is atmospheric, becomes 0 and . So, for units exposed to the atmosphere in the east, the hydraulic head is equal to the land elevation. Consequently, the eastern BC for the 27, 18, 11, 7, and 4.5 Ma models of the sequence is a constant head BC set equal to the land elevation for these models in the east (see arrows in Figure 13). The same constant head BC () in the east was set for all layers in the model, which implies hydrostatic conditions at the eastern boundary.
However, for the 3.5, 3, 2, and 0 Ma models, the eastern BC values were calculated by interpolation between the 4.5 Ma model hydraulic head values and present-day hydraulic head values obtained from potentiometric surface maps of the Wolfcampian and Mississippian (Figures and ; [7]); see arrows in Figure 13. The Wolfcampian at the east end of the model crops out at a land surface elevation of about 305 m ([7], Figure 14A).
The Wolfcampian head is, therefore, fixed to of 305 m in the last column of the model in the east at 0 Ma. Wolfcampian head values were then interpolated between 304.785 m and the Wolfcampian head of 310.906 m in the easternmost column calculated at 4.5 Ma by the 4.5 Ma model. These intermediate temporally interpolated values, , of the Wolfcampian heads at 3.5 Ma, 3 Ma, and 2 Ma were then used to set the constant head BC in the east for the Wolfcampian unit for the 3.5, 3, 2, and 0 Ma models. These same interpolated values of the Wolfcampian heads were used to set the constant head BC in the east for the top seven layers (which include the Wolfcampian as layer , Table 2) of the 3.5, 3, 2, and 0 Ma models, assuming equal heads in these layers.
The eastern BC for the bottom 3 layers in the 3.5, 3, 2, and 0 Ma models was established in a similar fashion, as , but with using calculated and actual contoured values of the Mississippian instead of the Wolfcampian, with a value of 183 m at 0 Ma, based on contoured values of Figure 11.
It should be noted that when the head is fixed to a specific value at the eastern boundary, MODFLOW varies the flux at the boundary, as needed, to maintain the specified head value, according to Darcy’s equation (see (2)). This affects the head gradient terminating at the boundary.
The values of specified flow rate in the west, , , and ( and ), resulting in calculated present-day hydraulic head values closest to the actual present-day hydraulic head values, are referred to as “optimal” values. Even though we refer to these values as “optimal,” it should be noted that this solution is not unique, and other combinations of , , and values could have produced calculated present-day hydraulic head values as close, or closer, to the actual present-day hydraulic head values. The “optimal” , , and ( and ) values were used in obtaining the final modeling results presented in this paper.
The sensitivity of the calculated values to assumed flux (), , and ( and ) values are presented in Sensitivity Analysis.
3.3.4. Initial Heads for the Sequence of Models
The initial heads at 70 Ma, prior to any uplift, represent the hydrostatic head of the Cretaceous Western Interior Seaway that was assumed to overlie most of the study area. That hydrostatic head is 0 m if we take mean sea level to be the vertical datum.
As was presented earlier, hydraulic head (), in units of L, is given as shown in (A.1).
It can be seen from (A.1) and Figure 1 that, in the aquifer at 70 Ma, which is saturated by the sea up to its top, at the top is 0 because (no water above it) and (the aquifer top is defined as sea level, the datum specifically for the purpose of the discussion in this section). For a point in the aquifer at a depth , and ; therefore, . So, the initial hydraulic head is 0 everywhere in the aquifer at 70 Ma and is taken to represent mean sea level. Referring to the values at all the cells in the model symbolically as the “vector” , we can say that , where means “initial.”
The land surface for the first model, the 27 Ma model shown in Figure 14 used to model the 13-million-year period from 40 to 27 Ma, is raised 250 m ASL in the west (to accommodate the Florissant fossil data) and 50 m ASL in the east, relative to the Cretaceous Western Interior Seaway. We assume that, at modeling time 0 years (40 Ma), there is a sudden (instantaneous) uplift in the west of 250 m and in the east of 50 m. According to (A.1), will be increased by 250 m in the west, decreasing linearly to 50 m in the east. Therefore, the initial at modeling time 0 years (40 Ma) is from the initial heads in the aquifer at 70 Minearly decreasing values from 250 m in the west to 50 m in the eas. The values are the same in all cells of the same column but decrease from 250 m for column 1 in the west to 0 m for column 26 in the east.
With these initial head values at 40 Ma, as described above, and the boundary conditions and and parameters presented earlier, MODFLOW calculates for each of the 13 million years (40 to 27 Ma) that the first model, the 27 Ma model, represents, ending at 27 Ma with .
For the 18 Ma model shown in Figure 15, a sudden (instantaneous) additional rise, , equal to the uplift from 27 to 18 Ma, , takes place along the length of the model, at 27 Ma, creating an equal rise in , according to (A.1). The incremented values, , represent the initial head values for the 18 Ma model at 27 Ma.
With the above initial head values for the 18 Ma model, + , at 27 Ma, and the boundary condition presented earlier, MODFLOW calculates for each of the 9 million years from 27 to 18 Ma that the 18 Ma model represents, ending at 18 Ma with . The above process is repeated until MODFLOW calculates , present-day values at 0 Ma. In each time step, MODFLOW calculates the response to a step-function increase in head imposed at the beginning of that time step.
3.3.5. Transience While Running the Sequence of 9 MODFLOW Models
When running the first of the 9 MODFLOW models, the 27 Ma model representing the 13-million-year period from 40 to 27 Ma, the transient solution at 30 Ma was close enough to the steady-state solution (which was calculated separately) that the 30 Ma transient solution was considered the final solution for the 27 Ma model. Figure 20 shows that the 27 Ma model solution (in the Wolfcampian) is largely at steady state at approximately 34.5 Ma (5.5 million years of modeled time from 40 Ma). This modeled duration to steady state is bracketed by the analytic solution given in Carslaw and Jaeger ([17], Eq. 3.3.14) and presented in Table 6. Because of the Wolfcampian, from Table 5, varies between 1 × 10−4 and 5 × 10−4 m/d, the analytic solution indicates from Table 6 that the decay time until steady state is between 14.9 and 1.49 million years.
The remaining 8 MODFLOW models, the 18, 11, 7, 4.5, 3.5, 3, 2, and 0 Ma models, were run for their full durations of 9, 7, 4, 2.5, 1, 0.5, 1, and 2 million years, respectively (see Table 3), whether they reached steady state or not. The models with longer modeling periods, specifically the 18, 11, 7, and 4.5 Ma models, reached steady state within their modeling periods. The models with shorter modeling periods, specifically the 3.5, 3, 2, and 0 Ma models, did not reach steady state. The hydraulic heads at the end of a modeling period were used to create the initial heads for the following period, as described earlier. It did not matter whether the heads at the end of a modeling period were at steady state or not.
3.3.6. Sensitivity Analysis
Using model 1 (40 to 27 Ma) with optimal parameters as a base case (see Table 5 for values of ), we computed perturbed profiles corresponding to perturbed values of flux on the western boundary (Figures 21 and 22). Doubling the flux produces perturbed values greater than the base case, and halving the flux produces perturbed values less than the base case. At 0 km, the changes in range from +5 m in the Mississippian unit to +32 m in the Wolfcampian unit. The relatively low changes in in the Mississippian unit are due to its westward pinch-out and lack of direct hydraulic connection to flux at the western boundary (Figure 7). Perturbations tend to zero at the eastern boundary where the solution is constrained by the constant head boundary condition. The base case flux values determined for the western boundary are quite low: 109,500, 109,500, and 255,500 m3 per million years for the Wolfcampian, Desmoinesian, and Virgilian-Missourian units, respectively. If these fluxes were to penetrate uniformly into a horizontal prism and fill only one percent of rock volume, the -length of penetration would be 719, 599, and 762 m, or less than one kilometer per million years. In other words, the total penetration, assuming that one percent of rock volume is utilized, would be about 20 to 30 km for the 40 m.y. duration of flux. Our model suggests that fluxes from the west are of little significance volumetrically even though the rates are highly constrained.
A second concern is the sensitivity of the model to its fixed parameters and (Figures 23–26). The key parameter in (1) is the ratio , suggesting that perturbations in and should produce opposite effects, and this is indeed the case. An increase in drops below the base case (Figure 23), and so does a decrease in (Figure 26). Conversely, a decrease in and an increase in both produce increases in (Figures 24 and 25).
A third sensitivity test examines the effect of a much higher in the uppermost layers of the model (Figure 27). of Pierre Shale is maintained at = 10−6 and the overlying and underlying units are increased to 10−4. Figure 27 shows that if in the Ogallala, Laramide Basin Fill, and Dakota were to be increased by 2 orders of magnitude, the conceptual model of isolation of the deeper units would be compromised, in the Wolfcampian would change by −8.95 percent (Table 7), and the heads in the deeper aquifers would not be matched as well.
Results of the sensitivity analysis given above are summarized in Table 7.
3.3.7. Model Limitations
As stated under Creating and Calibrating the Hydrologic Model to Simulate Structurally Driven Underpressures, even though we refer to values of and ( and ) obtained by the sequence of models as “optimal,” it should be noted that the sequence of models and their solution presented here is not unique, and other combinations of flux, , and values could have produced calculated present-day hydraulic head values as close, or closer, to the actual present-day hydraulic head values.
The model geometry is rather simple, or coarse, with only 26 columns in the direction covering a horizontal distance of 800 km and 10 layers in the direction covering a vertical distance of 4 km, and aquifer properties ( and ) are assumed to be homogenous in each model cell. Likewise, the one-million year time steps of the diffusive hydrologic model could be shortened as the uplift history becomes more refined. Shorter time steps and keeping the diffusive and uplift (model) time steps the same would result in better estimates of equilibration times.
The 2-dimensional vertical nature of the models assumes only flow in the horizontal (length) and vertical (depth) directions and that no flow takes place in the horizontal direction.
4. Results and Discussion
We have selected four of the nine models to show the buildup of land surface and hydraulic head with time: the 27, 18, 4.5, and 0 Ma models, which sample the periods of nondeposition, deposition of the Ogallala, and erosion (Figure 8). Figures 28–31 display the land surface and hydraulic head at the end of each time period, as well as the present-day land surface and hydraulic head profiles. The sequence of models utilizes the optimal values of flux from the west and hydraulic conductivity values, as previously discussed. The four figures show the progression from zero head values at 40 Ma to the near-match of underpressured head values at 0 Ma for the Wolfcampian, Desmoinesian, and Mississippian hydrologic units.
At 27 Ma, head values for the three hydrologic units are approximately 100 m at km, less than the land surface elevation of 200 m (Figure 28). Thus, heads at 27 Ma are roughly half of land surface along most of the cross-section, showing the early development of underpressure. The eastern boundary condition ( = surface elevation) prevents the development of any underpressure at km. At 18 Ma (Figure 29), the amount of underpressure, equivalent to the separation between land surface and head, increases in the west and separation among the three modeled heads—Wolfcampian, Desmoinesian, and Mississippian—is more visible. These four trends—(1) higher land surface elevation; increasing head; increasing underpressure in the west; and increasing differentiation among the three hydrologic units—are progressively more manifest at the end of 4.5 Ma (Figure 30) and 0.0 Ma (Figure 31). In addition, we observe a concavity developing in the modeled Wolfcampian and Desmoinesian hydraulic head curves at 24 Ma, similar to what is seen in the present-day data. At 0 Ma (Figure 31), the modeled hydraulic head values of the Wolfcampian, Desmoinesian, and Mississippian units approximate their corresponding present-day hydraulic head values.
Figure 32 presents the progression of seven of the nine modeled Wolfcampian hydraulic head values for the 27 Ma through 0 Ma models (for visible clarity, the 3.5 and 3 Ma models are not shown). Figure 32 shows that, through the series of nine models, the Wolfcampian modeled head values gradually approach, and eventually approximate, the present-day Wolfcampian measured hydraulic head values, also shown in Figure 31. Note that the hydraulic head values for the 2 Ma model are larger than the 0 Ma values. The decrease in the hydraulic head values from 2 to 0 Ma is the result of fixing the hydraulic head values on the eastern boundary to the geohydrologic unit elevations, as described under Boundary Conditions (BCs) for the Sequence of Models.
Figure 33 presents the modeled hydraulic head values with time for 3 model columns (6, 11, and 23; located at 140, 410, and 740 km) of the Wolfcampian layer for the period from 2 to 0 Ma. Figure 33 shows that the hydraulic head values at 0 Ma have not reached steady state. This is especially clear when compared to Figure 20, where steady state does occur during the modeled period at all locations. This result is significant because it indicates that present-day hydraulic heads are not at equilibrium and that, consequently, underpressures are going to increase in the future (beyond 0 Ma). The pattern uncovered by the series of 9 MODFLOW models is of increased underpressures with each subsequent model, as discussed earlier when comparing the underpressures in Figures 28, 29, 30, and 31. Overall, the models indicate that the tectonic uplift plays a primary role in the development of underpressures in the Great Plains.
5. Summary and Conclusions
From the geologic evidence, we have constructed an uplift history for the Great Plains plotted as cumulative uplift as a function of time from 27 Ma to the present, with superimposed periods of deposition and erosion.
The uplift history is embedded in 9 two-dimensional cross-sectional models, based on a geologic cross-section running from northeast Colorado to southeast Kansas. Tectonic uplift, deposition, and erosion were incorporated as changes in the surface elevation, geohydrologic unit elevation, and hydraulic heads at the beginning of each MODFLOW run. Fluxes were specified on the western boundary and constant heads were specified on the eastern boundary. The hydraulic heads at the beginning of each MODFLOW run (the initial heads) were assumed to be the values from the previous model run plus the increase in elevation due to tectonic uplift. This process was repeated until the calculated heads of the 0 Ma model, that is, the calculated present-day heads, were obtained. Hydraulic conductivity values and western specified-flux values were varied until the present-day calculated heads were as close as possible to known (measured) present-day heads.
The results from the sequence of MODFLOW models indicate a gradual increase in hydraulic heads and development of underpressures. Early models, such as the first model (40 to 27 Ma), indicate that the modeled Wolfcampian, Desmoinesian, and Mississippian heads are substantially different from the corresponding present-day Wolfcampian, Desmoinesian, and Mississippian hydraulic head data. Later models, such as the 7 to 4.5 Ma model, indicate that the modeled hydraulic head values of the Wolfcampian, Desmoinesian, and Mississippian units continue to rise, though not as rapidly as the land surface. Consequently, comparison of the modeled hydraulic head values to the 4.5 Ma land surface shows a continued increase in the underpressures and also indicates development of a downward concavity of the Wolfcampian and Desmoinesian hydraulic heads in the west, similar to what is seen in the present-day data. The final model (2 to 0 Ma) indicates that the modeled hydraulic head values of the Wolfcampian, Desmoinesian, and Mississippian units approximate their corresponding present-day hydraulic head values. The model also provides a good representation of the present-day underpressures measured in the Wolfcampian, Desmoinesian, and Mississippian units. The modeled and measured hydraulic head values also indicate that the underpressures increase to the west.
The 2 to 0 Ma model indicates that the present-day hydraulic head value of the Wolfcampian geohydrologic unit has not reached steady state. This result is significant because it indicates that present-day hydraulic heads are not at equilibrium and, therefore, underpressures are going to increase in the future. The pattern uncovered by the series of 9 MODFLOW models is of increased underpressures with time. Overall, the models indicate that tectonic uplift explains the development of underpressures in the Great Plains.
Appendix
A. Theoretical Basis of Underpressure
To understand what underpressure means, the following definition of hydraulic head is needed. Hydraulic head , in units of L, is given bywhere is the hydraulic pressure at a location [F L−2], is the unit weight of water [F L−3], and z is the elevation (“elevation head”) above a datum [L] (L represents units of length and F units of force). In (A.1), has units of L and is the “pressure head,” representing the height of a column of water that produces the pressure at its bottom for normal hydrostatic pressure conditions. Equation (A.1), therefore, says that the hydraulic head equals pressure head plus elevation head.
For a hydrostatic system, where the pressure head at a point is equal to the depth of submergence of the point, the hydraulic head is constant throughout, including at the water table. Alternately, at any point in an aquifer, if the fluid pressure is equal to normal hydrostatic pressure (i.e., pressure head is equal to the depth of submergence), the hydraulic head at that point is equal to the hydraulic head at the water table.
Underpressure, on the other hand, is the condition where pore fluid pressure (hydraulic pressure) is below normal hydrostatic pressure for a point in the aquifer (the pressure : depth ratio is less than that of a hydrostatic column; [7]). Hence, it follows that, for a depth where fluid pressure is below normal hydrostatic pressure, the hydraulic head is less than the hydraulic head at the water table, as depicted in Figure 1. Note that, in Figure 1, is the pressure head at the observation point relative to the potentiometric surface, not the water table.
Thus, the following definition can be made: underpressure is the condition where the hydraulic head is less than the hydraulic head at the water table. The amount that the hydraulic head at a point is less than the hydraulic head at the water table is the amount of underpressure, represented by the term Up_real in Figure 1.
Additionally, if the vertical distance between the water table elevation (WT) and land surface elevation (GL), namely, , as shown in Figure 1, is negligible relative to GL like in the case of deep aquifers with shallow water tables, then the hydraulic head at the water table is approximately equal to the land surface elevation. In these cases, which are assumed in this report, the above definition of underpressure can be restated as follows: underpressure is the condition where the hydraulic head is less than land surface elevation. The amount by which is less than GL is the approximate amount of underpressure, depicted as Up_approx in Figure 1. The data and model results in this paper are presented with reference to hydraulic head (h) and land surface elevation (GL), from which Up_approx can be calculated as .
B. Expanded Geologic Setting
The central Great Plains are a vast, almost monotonically eastward-sloping terrain that stretches from the front of the southern Rocky Mountains eastward to the Mississippi River at the meandering eastern and northeastern boundaries of Nebraska and Kansas, respectively (Figure 2). Elevations range from about 200 meters above sea level in the east to almost 2000 meters locally at the western limit of the Great Plains. Slopes within this tilted terrain are dominantly low (gentle) but become progressively steeper to the west. The resulting form is very gently concave upward relative to the steep but nonvertical axis of this very long wavelength and very low amplitude tilted synclinal fold (Figure 7). No closed low is present in the center of this noncontractional fold because the tilting greatly exceeds the folding; thus, it is not a true syncline.
Eaton [11, 12] interpreted the Great Plains as the eastern flank of a huge antiformal uplift (Figures 3 and 4) that is centered on the Rio Grande rift and includes the Colorado Plateau as its western flank. The Rio Grande rift (Figures 4 and 5) in turn is partly coincident with the southern Rocky Mountains (Figure 4), and the current high elevation of this mountain range dates only to uplift of this antiformal bulge, which entirely postdates the original formation of these mountains in the Laramide (~70 to ~40 Ma) orogeny. Eaton [11] called this giant domal uplift the Alvarado ridge, a name that has not caught on but which nevertheless represents an important observation. Timing constraints indicate that uplift of the Alvarado ridge followed a similar history throughout its geographic range [11, 12], indicating it is a unitary tectonic (epeirogenic) landform—despite the overwhelming focus in the literature on the uplift of the Colorado Plateau and the resulting incision of the Grand Canyon.
Rocks exposed across the central Great Plains include the upper Cenozoic (~18 to ~4.5 Ma) Ogallala alluvial conglomerate, Laramide fills of foreland basins (such as the Denver Basin, Figure 2), Upper Cretaceous marine shales and shoreline sandstones, and, in eastern Kansas, Paleozoic marine limestones and related rocks (Figure 6). Studies of the Ogallala Formation indicate that it has been strongly tilted in the western Great Plains and only weakly tilted in the eastern plains [18, 19], indicating that it is folded in a fashion that mimics the overall geomorphology of the Great Plains.
Moreover, several major subsurface structures (basins and anticlines) are present at depth under the Great Plains (Figure 2), and none of these have anything but erosional geomorphic expression, including the Laramide basins, indicating that the uplift of the plains is entirely late Cenozoic (after 40 Ma). The lack of geomorphic expression of the buried structures under the Great Plains is particularly clear in cross-section (such as in Figure 7). Further, the generally steep front of the southern Rocky Mountains is entirely an erosional landform in its current expression. The history of uplift of the Great Plains is recorded in its tectonostratigraphy.
B.1. Tectonostratigraphy and Hydrostratigraphy
B.1.1. Early Paleozoic
Cambrian through Mississippian rocks under the central Great Plains are marine limestones and dolomites and lesser interbedded shales deposited when shallow seas covered nearly all of the midcontinental United States. An interval of uplift and erosion at the end of the Mississippian resulted in removal of the lower Paleozoic rocks from east-central Colorado westward to the front of the southern Rocky Mountains and beyond. Over the remainder of the study area, where these rocks are still present, the early Paleozoic rocks have nearly ubiquitous development of a karstic weathering profile at the top.
This karstic horizon is the most transmissive aquifer/reservoir rock under much of the central Great Plains. Moreover, the whole early Paleozoic section appears to act as a single continuous geohydrologic unit as indicated by potentiometric continuity [7]. The uninterrupted nature of the karstic profile results in excellent continuity of this unit over much of the study area. None of the interbedded shales appears to act as a significant confining unit, at least on a regional basis. Because of the westward pinch-out under the Plains, the lower Paleozoic rocks receive no recharge from the west. Further, these strata are capped by a basal Pennsylvanian section which is a confining unit, strongly restricting recharge from overlying strata.
B.1.2. Late Paleozoic
Pennsylvanian and Permian strata under the central Great Plains are mainly marine rocks in the east changing westward to almost exclusively terrestrial rocks in the far west at the front of the southern Rocky Mountains and for at least tens of kilometers to the east of that front. Under at least half of the central Great Plains, the Pennsylvanian-Permian rocks exhibit a complex intertonguing of marine and terrestrial rocks. The purely terrestrial rocks are clastic sediments derived from Pennsylvanian-Permian uplifts that formed at the time, mainly in central and western Colorado, which are known as the Ancestral Rocky Mountains. Because of the abundant clastic-sediment supply that existed, a majority of the upper Paleozoic marine rocks are also clastic sediments, although they contain usually thin, but very widespread, carbonate beds. Extensive and locally thick evaporite deposits are present at the marine/terrestrial contact in many parts of the study area, reflecting episodic marine incursions.
The Pennsylvanian-Permian interval was one characterized by numerous marine transgressions and regressions, reflecting repetitive episodic sea-level fluctuations related mainly to glaciations but also to tectonic events in some cases. Because of the very low relief of much of the midcontinent in this interval, the marine shoreline moved hundreds of kilometers westward and then back eastward in these transgressions and regressions, respectively. Consequently, the marine carbonate horizons are aquifers/reservoir rocks with regional-scale physical and hydrologic continuity in many cases, much like the karstic lower Paleozoic hydrologic unit.
Unlike the lower Paleozoic rocks, however, the upper Paleozoic sections contain a number of confining units, mainly shales. On a regional basis, two of these confining units are important. One of these separates the Wolfcampian, Virgilian, and Missourian aquifers from the lower part of the Pennsylvanian, as indicated by differing pressure regimes [7]. The second, part of the Upper Permian section, consists almost entirely of shale and evaporites and is a confining unit restricting recharge to the upper Paleozoic strata from above. Further, recharge from the west is restricted by the westward facies change to exclusively terrestrial rocks, both because the latter are low-permeability rocks in general and because the permeable beds within the section lack continuity across the marine/terrestrial transition.
Many of the arches and basins that are buried under the central Great Plains, such as the Las Animas arch (Figure 2), were active during the late Paleozoic. Consequently, the upper Paleozoic section, and especially the carbonate beds in this section, thin over this arch and other related features. This thinning over buried uplifts is accompanied by an increase in the clastic content of carbonate beds in the section, resulting in lower permeability in these areas.
B.1.3. Early and Middle Mesozoic
The Triassic, Jurassic, and Lower Cretaceous section in the midcontinent is very thin in general and is discontinuous in many areas. This section consists of terrestrial sediments which are generally low in permeability. Hydraulically, this section is continuous with the Upper Permian shale and evaporite section and acts as part of this confining unit.
B.1.4. Late Cretaceous
The vast majority of Upper Cretaceous strata under the Great Plains are black shales, which were deposited in the Cretaceous Western Interior Seaway. These shallow-water pelagic deposits (mainly the Pierre and Niobrara Shales, Figure 7) are underlain and capped by horizons which record a major transgression and a major regression near the beginning and end of the Late Cretaceous, respectively. These events resulted in deposition of widespread marine shoreline sandstones over the area of the Great Plains—the Dakota and Trinidad Sandstones and their stratigraphic equivalents [20]. These marine shoreline deposits as well as the shallow-water nature of the thick intervening black shales prove that, for tens of millions of years prior to the onset of the ~70- to ~40-Ma Laramide orogeny, the area of the Great Plains was close to sea level.
The Upper Cretaceous black shale is very thick over all but the eastern quarter of the study area (Figure 7). This shale is the most effective confining unit in the study area, strongly restricting recharge to all of the aquifers/reservoir rocks that underlie it. The Dakota Sandstone is an aquifer that is much thinner and less permeable than the underlying Paleozoic aquifers/reservoir rocks. The extreme vertical exaggeration of the cross-section (Figure 7) precludes showing that this unit is bent up to the surface immediately east of the front of the Rocky Mountains, where it crops out in a narrow hogback. The Dakota Sandstone receives minor recharge from this small-footprint exposure but receives no recharge from the mountains to the west because it does not quite extend to the mountains.
B.1.5. The Laramide Orogeny
Laramide strata of the Great Plains are the fills of foreland basins, such as the Denver Basin, which formed in the ~70 to ~40 Ma Laramide orogeny—during the contractional uplift of the Rocky Mountains and the local formation of the Golden Fault at the mountain front (not shown in Figure 7, but western end of geologic cross-section is at the mountain front). Numerous workers have observed the predominantly fine-grained nature and thus low-energy depositional environment of these terrestrial basin fills. This is consistent with other evidence (e.g., [21]) indicating that these foreland basins, and thus the area of the future Great Plains itself, remained near sea level during the initial contractional uplift of the adjacent Rocky Mountains. The Laramide strata are present only in the westernmost part of the study area (Figure 7) and contain no major aquifers. Where present, this section adds extra confinement to the deep aquifers/reservoir rocks under the Plains, beyond that provided by the Pierre and Niobrara Shales.
B.1.6. Post-Laramide Time
From shortly after the end of the Laramide orogeny (at ~40 Ma) to just before the beginning of formation of the Rio Grande rift (at ~27 Ma), subduction continued along the west coast of North America; however, contractional uplift of the southern Rocky Mountains had ceased. This was an interval of widespread subduction-related volcanism in and around the southern Rocky Mountains [22]. Eaton [11] has inferred that uplift of the Alvarado ridge began during the interval of ≤40 to ~27 Ma; however, he provides no evidence to support that interpretation. In the few cases where data exist that can be directed at testing this minor element of Eaton’s [11] interpretation, no support was found for onset of uplift of the Alvarado ridge, including the variable uplift/tilting of the Great Plains, before initiation of Rio Grande rifting at ~27 Ma.
Late Eocene/early Oligocene sediments (~34 Ma) in the area of Florissant, Colorado, are currently at ~2500 m elevation but were only ~300 m above sea level at the time of deposition, based on an analysis of plant fossils [23]. This fossil evidence is somewhat controversial but provides supportive evidence that major uplift of the Alvarado ridge did not occur until after the interval of subduction-related volcanism during ≤40 to ~27 Ma. Lava vesicle paleoaltimetry from flows of different ages on different parts of the Alvarado ridge [24] provides more definitive evidence and is consistent with the above conclusion.
The onset of (Rio Grande) rifting to the southwest of the study area (near the central Colorado/New Mexico border) at ~27 Ma was accompanied by an early stage of uplift and of resulting erosion [13]. Erosion in the highest part of the Great Plains in the southern part of the study area—across most of northeastern New Mexico and in a small part of adjacent southeastern Colorado—extended to ~18 Ma, when deposition of the Ogallala Formation began. In contrast, in at least the early part of the same interval, the more northerly Great Plains, in the vicinity of the Colorado/Wyoming border, remained a depocenter, mainly for fine-grained deposits of the Eocene to Oligocene White River and late Oligocene to early Miocene Arikaree Formations, which were evidently deposited in a swampy environment [13]. The formation of the Rio Grande rift itself was diachronous, with onset starting in southern New Mexico at ~30 Ma and progressing northward to the area of the Colorado/Wyoming border at ~25 Ma (R. A. Thompson, USGS, oral comm., 2015).
Erosion/deposition patterns on the Great Plains indicate that uplift/tilting of the eastern flank of the Alvarado ridge followed a northward-migrating diachronous pattern that was similar to that of the onset of rifting. Whereas onset of uplift of the Great Plains was slightly staggered behind that of the Rio Grande rift, the time of maximal uplift of all parts of the Alvarado ridge were strongly staggered behind the predominant time of formation of the rift. For example, the tectonostratigraphy of the Great Plains indicate that approximately half of the uplift and tilting occurred from latest Miocene time (~6 Ma) to the present, as discussed below. Whereas a small number of data points, mainly from the Grand Canyon, have been interpreted as indicating a much earlier uplift of the western part of the Alvarado ridge (at ~70 Ma; [25]), most workers find this conjectural interpretation to be irreconcilable with the majority of available data, most of which are more reliable and robust than the questionable uplift data from the Grand Canyon (e.g., see [26]).
The roughly east-west cross-section that we have constructed as part of our framework model (Figure 7) lies approximately halfway between the Colorado/New Mexico and Colorado/Wyoming border areas in which the ~27 to ~18 Ma interval was one of erosion and one of deposition, respectively. On the line of section itself, which was chosen to be along the section line published in Nelson et al. [7], there is no evidence as to whether erosion or deposition was occurring in this interval, and thus we have assumed nondeposition in this interval along the line of section as the simplest assumption from available evidence.
The interval from ~18 to ~4.5 Ma on the Great Plains was characterized by deposition of a vast alluvial apron covering most or all of the Plains. These alluvial deposits are known as the Ogallala Formation, which is an important water-table aquifer in a major part of the Plains. The early stages of uplift of the Alvarado ridge precipitated formation of this alluvial apron [14]. Continued uplift and tilting eventually created grades steep enough to end deposition of this apron and to initiate erosion. Strong tilting of the Ogallala in the western part of the Great Plains indicates that at least about half of the uplift and tilting followed the transition from deposition to erosion.
The major area of preservation of the Ogallala alluvial apron is in the central part of the Great Plains. The base of the formation in this central area is a concave-upward surface that is almost certainly a tectonic fold (see Figure 7). As discussed above, the Great Plains are not simply tilted eastward, but they have also been bent into a concave-upward surface during epeirogeny as a result of the nonlinear nature of the uplift. Additionally, the major mass of preserved Ogallala Formation probably has escaped erosion at least partly because it is relatively protected from erosion within the central part of the tectonic concavity. Near the Colorado/Wyoming border, the Ogallala Formation is locally preserved in the far western part of the Great Plains in a feature known as the “gangplank.” Tiny remnants of the Ogallala Formation are also found capping small knolls in eastern Kansas [27]. Thus, it appears that the Ogallala Formation was deposited over virtually the entirety of the Great Plains during an early stage of uplift of the Alvarado ridge, based on evidence from McMillan et al. [28].
The vast majority of the erosion after 27 Ma on the Great Plains evidently occurred from ~4.5 Ma to the present—after deposition of the Ogallala Formation ceased. A number of different, mainly Quaternary, formations were deposited locally on the Great Plains during this interval, mainly as the caps of terraces formed during downcutting. Erosion rates can be estimated from the relative heights of these different-age deposits to yield estimates of the change in erosion rate during the interval of ~4.5 Ma to the present. However, the rates so obtained are all within the error bars of estimation (E. Taylor, USGS, oral comm., Denver), and thus the simplest assumption that satisfies existing data is a constant erosion rate since ~4.5 Ma.
Current slopes across the Great Plains average ~10−3 and, in at least the steeper western half of these plains, are well in excess of the 10−4 to 10−3 slopes [28] estimated during the time of deposition of the Ogallala Formation. This and the sustained pattern of downcutting throughout the Great Plains since ~4.5 Ma reflect a much stronger pattern of uplift since ~4.5 Ma than in the preceding ~27 to ~4.5 Ma. Whereas it is difficult if not impossible to quantify this, it appears likely that at least about half of the total uplift of the Great Plains occurred in the interval of ~4.5 Ma to the present, which is the same approximate timing that has been estimated for the western half of the Alvarado ridge (the Colorado Plateau; [26]) and within the area of the Rocky Mountains [29, 30].
Paleoelevation estimates based on basalt vesicle studies [24] indicate that about 50 percent of total uplift of the Alvarado ridge occurred from ~27 to ~5 Ma and the remaining half from ~5 Ma to the present. Further, studies in the vicinity of the “gangplank” show that late Pliocene sediments are not discernibly tilted, although only strong tilting is discernible. Thus, we have concluded that uplift peaked in the Pliocene and has been lesser in the Quaternary.
C. Permeability and Hydraulic Conductivity Values from Other Sources
In this appendix, values of permeability (and equivalent hydraulic conductivity values) chosen for two previous hydrologic modeling studies of the Great Plains are summarized. These values can be compared with values chosen for the model.
The underpressured deep brine aquifers of the Palo Duro Basin in the Texas Panhandle were studied with a two-dimensional model to understand the causes of underpressuring and the nature of deep groundwater circulation [31]. Figure 19 shows the permeability values used in their model for Triassic, Permian, and Pennsylvanian strata. The most important rock unit of the deep basin brine aquifer, the Permian-Pennsylvanian shelf carbonate facies, was assigned a value of 0.013 m d−1 based on values for Wolfcampian carbonates, Pennsylvanian carbonates, and pre-Pennsylvanian strata. The ratio of vertical to horizontal permeability was set to 0.01 for the shelf carbonate facies. Other values of range from 0.01 in a salt dissolution zone to 0.1 in a fluvial-lacustrine system to 1.0 in low-permeability salt and basinal systems.
In considering the selection of permeability values for a regional study of the Great Plains, Signor et al. [16] produced a series of maps of for each of their aquifer systems. Values of were computed from porosity values obtained from geophysical logs. Values from their maps, taken along the cross-section shown in Figure 7, are summarized in Figure 18. The lower units of the Western Interior Plains Aquifer System of Signor et al. corresponds to rock units of Ordovician and Silurian age, while the upper units correspond to rock units of Mississippian age. In the present paper, both the lower and upper units of the Western Interior Plains Aquifer System are combined into the Mississippian and Cambrian-Ordovician-Silurian hydrologic unit. The third grouping of Figure 18, the Western Interior Plains Confining System, comprises rock units ranging from Pennsylvanian to Jurassic in age, which are treated in the present paper as five distinct hydrologic units: Basal Pennsylvanian, Desmoinesian, Missourian-Virgilian, Wolfcampian, and Upper Permian-Jurassic. The Lower Cretaceous is named the Dakota hydrologic unit in this paper, but it was not modeled.
Definitions of Terms
| : | Vector (list) of values |
| : | Vector of zeros |
| : | Gamma; unit weight of water [F L−3] |
| Actual present-day hydraulic head values: | Actual present-day hydraulic head data values observed in the field or derived from field observations |
| ASL: | Above sea level |
| BC: | Boundary condition |
| c1: | Subscript indicating cell 1 |
| c2: | Subscript indicating cell 2 |
| Calculated present-day hydraulic head values: | Hydraulic head values calculated for the ninth MODFLOW model in the sequence, at 0 Ma (the present) |
| cell: | Basic element of each MODFLOW model at the intersection of a layer and column |
| cell 1: | A cell west of cell 2, for illustrative purposes |
| cell 2: | A cell east of cell 1, for illustrative purposes |
| : | Head contour value of the Mississippian in the east |
| Column: | Horizontally discrete element of each MODFLOW model. Each MODFLOW model is divided into 26 columns spanning the horizontal direction |
| Conductivity: | Qualitative reference to hydraulic conductivity |
| CpTU: | Cumulative percent of total uplift |
| Cretaceous Western Interior Seaway: | The seaway, interior to the North American continent, connected to the present-day North Sea and Pacific and Atlantic oceans, which prevailed during the Cretaceous and submerged present-day Colorado and half of Kansas |
| Cumulative percent of total uplift (CpTU): | The total uplift (TU) at a particular geologic time, since 27 Ma, as a percent of total uplift |
| : | Head contour value of Wolfcampian in the east |
| d: | Day |
| : | Used for depth under land surface in discussion of Figure 1 |
| : | Darcy, unit of |
| Dak: | Dakota Formation |
| Data: | Present-day measured values of Wolfcampian, Desmoinesian, and Mississippian hydraulic head values |
| Data curve: | Wolfcampian, Desmoinesian, or Mississippian curve of present-day measured hydraulic head values as a function of |
| : | Hydraulic gradient or rate of change of with respect to . See (2) |
| : | East |
| East: | The right side of all figures in the report is aligned with the geographic east |
| Eq.: | Equation |
| Erodible thickness, ET: | The thickness of the model to be eroded from 4.5 Ma (max value) to the present, in m |
| F: | Dimensional unit of force, as in [F] |
| FE, FE#: | Flowchart Element, or element #, in Figure 13 |
| Flux: | , flow rate [L3 T−1] per unit area [L2], so equal to [L T−1] |
| Formation: | The same as “geohydrologic formation” |
| Formula: | This terminology was used exclusively in this paper in relation to formulas for use with ModelMuse, a preprocessor to MODFLOW. A formula, for example, can be an algebraic equation describing how is a function of . Further description of these formulas can be found in Winston [15] |
| ft: | Foot or feet |
| : | Acceleration due to gravity |
| Gamma: | See |
| Geohydrologic formation: | A grouping of rocks considered to have the same geohydrologic parameters and given the geologic name, or age, of one or more of the rocks in the grouping, the same as “geohydrologic unit.” Each layer of a MODFLOW model represents a geohydrologic formation (see Table 2) |
| Geohydrologic parameters: | The values of and used for the model sequence |
| Geohydrologic unit: | The same as “geohydrologic formation” |
| Geometry: | The geometry of a model is the layout of columns and layers of a model; see “MODFLOW model” |
| Great Plains: | Physiographic region of North America, bordered by the Southern, Middle, and Northern Rocky Mountain physiographic regions to the west and the Central Lowland physiographic region to the east |
| Grid: | The grid of a model, the same as “geometry” |
| : | Hydraulic head, [L], m |
| : | The subscript , as in , or , indicates the horizontal direction |
| : | Hydraulic head, , at cell 1 |
| : | Hydraulic head, , at cell 2 |
| : | Head specified in (east) |
| [distinct from ]: | Temporally interpolated head between at 4.5 Ma, calculated by the 4.5 Ma model, and at 0 Ma for layers through (7), and between at 4.5 Ma and at 0 Ma for layers (8) through (10). This is used as for the 3.5, 3, 2, and 0 Ma models |
| : | Initial for a model run, distinct from |
| : | Vector, or list, of hydraulic heads for all the cells in a model |
| : | at 18 Ma |
| : | at 27 Ma |
| Head: | The same as “hydraulic head” |
| Head plot: | Plot of head, , versus , the same as “hydraulic head plot” |
| Hydraulic gradient: | See |
| Hydraulic head: | ; summation of “elevation head” and “pressure head,” and in (A.1) |
| Hydraulic head plot: | A plot of hydraulic head, , versus |
| : | Subscript indicating “initial” in and “temporally interpolated” in |
| : | Hydraulic conductivity (depends on medium, , and fluid properties): includes and . [L T−1], m d−1 |
| : | permeability [L2] depends only on the medium (not fluid) properties; it includes horizontal, , and vertical, , components |
| : | Horizontal hydraulic conductivity |
| km: | Kilometer |
| : | Vertical hydraulic conductivity |
| L: | (1) Dimensional unit of length, as in [L], [L−1], etc. |
| Laramide Basin Fill: | Layer , second from top, of the model. Age: the combined Paleogene and Neogene Periods (formerly the Tertiary Period). Terrestrial gravel, sand, and clay (also called Tertiary Sediments) |
| Layer: | Vertically discrete element of each MODFLOW model. Each MODFLOW model is divided into 10 layers spanning the vertical direction and represents one or more geohydrologic unit(s) |
| Lbf: | Laramide Basin Fill |
| : | Land surface elevation in the east. |
| Leonardian: | Upper part of the Permian, combined with the lower part of the Jurassic (Morrison formation) in layer (5) |
| m: | Meters |
| m3: | Cubic meters |
| Ma: | Mega annum, or millions of years ago |
| model: | Same as “MODFLOW model” |
| ModelMuse: | Preprocessor program for MODFLOW: the user inputs values of MODFLOW model parameters, such as , and ModelMuse creates files in the format required for MODFLOW input files |
| Model sequence: | The same as “sequence of MODFLOW models” |
| Modeling time: | Time counted from the beginning of modeling, which starts at 40 Ma. So, modeling time of 0 is 40 Ma; a modeling time of 13 million years is 27 Ma because there are 13 million years from 40 to 27 Ma, and so on |
| MODFLOW: | Groundwater flow modeling program of the US Geological Survey |
| MODFLOW model: | One MODFLOW cross-sectional model in the sequence of models. For this report, it is a two-dimensional cross-sectional (vertical) representation of the geologic cross-section of the study area, divided into “layers” spanning the vertical direction and “columns” spanning the horizontal direction |
| m.y.: | Million years |
| Neogene Period: | Late part of the former Tertiary Period, consisting of the Miocene and Pliocene Epochs |
| no.: | Number |
| Og: | Ogallala Formation |
| Ogallala Formation: | Model layer , top layer. Age: Neogene period (old Miocene and Pliocene epochs). Terrestrial gravel, sand, and clay |
| OM: | Order of magnitude |
| Optimal values of underflow and : | Values of underflow and ( and ) that make the calculated present-day hydraulic head values as close as possible to the actual present-day hydraulic head values |
| : | Pressure, [F L−2] |
| Paleogene Period: | Early part of the former Tertiary Period, consisting of the Paleocene, Eocene, and Oligocene Epochs |
| Percent of total uplift: | PTU. Uplift between MODFLOW model times, such as between the 27 Ma and 18 Ma models, represented as a percent of total uplift |
| Potential: | Driving force in physics where a physical quantity like heat or mass moves from high to low potential values. In groundwater flow, the potential is the hydraulic head |
| Potentiometric surface map: | Map of the three-dimensional potential surface, represented as contours lines of equal potential, or equal hydraulic head |
| Present-day: | Describes current conditions, at 0 Ma, as opposed to in the geologic past |
| PTU: | Percent of total uplift |
| : | Flow rate [L3 T−1] |
| : | Flow rate specified at the west end of the model |
| : | Flux; flow rate [L3 T−1] per unit area [L2], so equal to [L T−1] |
| : | Flux specified at the west end of the model |
| : | Flux in the horizontal, , direction |
| resp.: | Respectively |
| : | Storage coefficient (storativity). The volume of water, L3, released from a unit area of the aquifer, L2, per unit drop in head, L. So, units of (L3 L−2) L−1 = 1, that is, dimensionless |
| Sequence of [MODFLOW] models: | A sequence of 9 MODFLOW models representing the structural, depositional, and erosional changes from 70 Ma to the present, the same as “model sequence” |
| : | Specific storage = per unit height of aquifer, L. So, units of 1 L−1, m−1 |
| STP: | Standard temperature and pressure |
| : | Time in millions of years, can be 27 Ma, 18 Ma, and so forth |
| T: | Dimensional unit of time, as in [T], [L T−1], and so forth in days, d, and transmissivity |
| The present: | Current, at 0 Ma |
| Total Uplift | (TU): Total uplift from 27 Ma to the present |
| TU: | Total uplift |
| Underflow: | Groundwater flow entering an aquifer from another aquifer bordering it |
| Underpressure: | The amount that the hydraulic head at a location is lower than the water table. For this report, this is approximated by the amount that the hydraulic head is lower than land surface |
| USGS: | US Geological Survey |
| : | the subscript , as in , indicates the vertical direction |
| Vector : | List of values |
| Vertical datum: | Reference for elevation in the vertical direction, taken in this report to be mean sea level |
| Virg: | Virgilian geohydrologic unit |
| : | West |
| West: | The left side of all figures in the report is aligned with the geographic west |
| Wolfc: | The Wolfcampian, Permian, geohydrologic unit |
| -direction: | The horizontal direction from west to east |
| : | Distance to the east of the Golden Fault (in m). Graphically shown in km |
| : | at cell 1 |
| : | at cell 2 |
| -direction: | The horizontal direction from south to north. For this 2-dimensional vertical model, there is no flow in the -direction and the model width in the -direction is 100 m |
| yr: | Year |
| : | Elevation head, or elevation ASL, in (A.1) |
| -direction: | The vertical direction, from low elevation to high elevation. |
Disclosure
Philip H. Nelson, Chris Fridrich, and Gary D. LeCain are Emeritus.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The authors are grateful to the US Geological Survey Carbon Sequestration Program, which funded this work for several years, and for its director, Peter Warwick who also reviewed the manuscript and provided valuable input. They thank Geoffrey Delin of the USGS for his thorough review and request to clearly present our model sensitivity analysis.