Simulating the Evolution of Fluid Underpressures in the Great Plains, by Incorporation of Tectonic Uplift and Tilting, with a Groundwater Flow Model

1Earth Systems Modeling Branch, Integrated Modeling and Prediction Division, Water Mission Area, US Geological Survey (USGS), Denver Federal Center, P.O. Box 25046, MS 406, Lakewood, CO 80225, USA 2Central Energy Resources Science Center, US Geological Survey (USGS), Denver Federal Center, P.O. Box 25046, MS 964, Lakewood, CO 80225, USA 3Central Energy Resources Science Center, US Geological Survey (USGS), Denver Federal Center, P.O. Box 25046, MS 980, Lakewood, CO 80225, USA 4Earth Systems Processes Division, Water Mission Area, US Geological Survey (USGS), Denver Federal Center, P.O. Box 25046, MS 406, Lakewood, CO 80225, USA


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: (1) the history of uplift and tilting, which is recorded in the tectonostratigraphy of the Great Plains, and (2) 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 ( = 1, . . ., 9), 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  + 1, 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  + 1.This procedure is repeated until  = 9 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.

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 0 500 (km) Figure 3: Index map showing key elements of the Eaton Alvarado ridge in relation to political boundaries (dash-dot lines) and physiographic province boundaries (thin solid lines); some of the latter provinces are grouped together and labelled.The region around the crest of the Alvarado ridge is outlined by a heavy solid line (Figure 2 of [11]).
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.

Understanding the Geology and Quantifying Structural
Changes.A geologic map of the study area and a crosssection along line A-A  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.Uplif t during the period of 40 to 27 Ma was assumed to be zero (Appendix B).

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 crosssection A-A  (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.

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.  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 threedimensional 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 (3a), (3b), . . ., (3i) 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-millionyear, 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 (1) the increased tilting of the land surface, (2) the increased tilting of all strata, and (3) 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).

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, x, 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 twodimensional cross-section of Figure 14), whereas the presentday land surface of Figure 7 is concave upwards.The upward 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 x = 340 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 x = 340 km to 0 m at both east (x = 810 km) and west (x = 0 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 x, 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 ( V ) values, until the calculated present-day hydraulic head values approximated the present-day hydraulic head values.This is described in the next two sections.

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, dℎ/dx, which aided us during calibration: where  ℎ is the horizontal hydraulic conductivity and ℎ is hydraulic head.Equation (2) indicates that, for a specific  ℎ ,  ℎ and dℎ/dx are inversely proportional.That is, as the medium becomes more conductive and  increases, the hydraulic gradient decreases and the plot of ℎ versus x becomes shallower or less steep.Conversely, a steep hydraulic gradient indicates a tight or low hydraulic conductivity medium.Based on this, the  ℎ and  V 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 etc.: Figure 13: Flow chart of the modeling process.ℎ  : head specified in east;   : Land surface elevation in the east; ℎ  : interpolated ℎ between {ℎ} 4.5 Ma and Wolfcampian ℎ contour in  for layers (1) to (7), and between {ℎ} 4.5 Ma and Mississippian ℎ contour in  for layers (8) to (10).
Box numbers in the figure after step (4) are for labelling and are not strictly sequential.
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  V =  ℎ /10 in most of the geohydrologic units (Table 5), a selection that is the default setting in MODFLOW.However,  V 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,  V , is large, the hydraulic head plots of ℎ versus x 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  V , 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, lowpermeability 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 largescale 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 onetwentieth 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.

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 m 3 d −1 of water (109,500 m 3 per 1 m.y.) enters each of the Wolfcampian (cross-sectional area 15,239 m 2 ) and Desmoinesian (cross-sectional area 18,287 m 2 ) and 0.0007 m 3 d −1 (255,500 m 3 per 1 m.y.) of water enters the Virgilian-Missourian (cross-sectional area 33,526 m 2 ).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 x 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 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 (6), 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  ME 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  V ), 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 presentday hydraulic head values.The "optimal"   ,   , and  ( ℎ and  V ) values were used in obtaining the final modeling results presented in this paper.
The sensitivity of the calculated ℎ values to assumed flux (  ),   , and  ( ℎ and  V ) values are presented in Sensitivity Analysis.

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).).In units ( 1)-( 7), the vertical hydraulic conductivity  V is set to one-tenth the value of the horizontal hydraulic conductivity  ℎ .In units ( 6) and ( 7), a low- choke separates more permeable units in the west and east.In unit (8),  is higher in the east than in the west.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  = 0 (no water above it) and  = 0 (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, ℎ = − = 0. 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 {ℎ}  = {0}, where  means "initial."

Geohydrologic unit 𝐾
The land surface for the first model, the 27 Ma model shown in Figure 14     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 {ℎ} 27 Ma .
For the 18 Ma model shown in Figure 15, a sudden (instantaneous) additional rise, , equal to the uplift from 27 to 18 Ma, Uplift 27-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, {ℎ} 27 Ma +Uplift 27-18 Ma , represent the initial head values for the 18 Ma model at 27 Ma.
With the above initial head values for the 18 Ma model, {ℎ} 27 Ma + Uplift 27-18 Ma , 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 {ℎ} 18 Ma .The above process is repeated until MODFLOW calculates {ℎ} 0 , 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.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.

Transience While Running the
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.

27-Ma land surface
Present-day land surface Present-day measured heads (3) Mississippian  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.

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  V ) 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. (3)

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).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  = 160 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 crosssection, showing the early development of underpressure.The eastern boundary condition (ℎ = surface elevation) prevents the development of any underpressure at  = 800 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; (2) increasing head; (3) increasing underpressure in the west; and (4) 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 presentday 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.

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 MOD-FLOW 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 presentday heads, were obtained.Hydraulic conductivity values and western specified-flux values were varied until the presentday 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 Geofluids 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.

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 by where  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, GL−WT, 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 GL − ℎ.

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 uplif t 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.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 ∼70to ∼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 A-A  is at the mountain front).Numerous workers have observed the predominantly finegrained 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 concaveupward 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 differentage 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 twodimensional 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  V / ℎ was set to 0.01 for the shelf carbonate facies.Other values of  V / ℎ 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.A cell west of cell 2, for illustrative purposes cell 2:

Definitions of Terms
A cell east of cell 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

Figure 2 :
Figure 2: Physiographic map of study area, in the Great Plains, from Nelson et al. ([7], Figure 1).The figure shows land surface elevation contours (m), state lines, location of structural features, sedimentary basins (shaded), and the 800 km cross-section line A-A  .

Figure 4 :
Figure 4: Cross-section of the Eaton Alvarado ridge along the 36 ∘ N parallel.The cross-section shows the Rio Grande rift incised into the Alvarado ridge (Figure 3 of [11]).

Figure 5 :
Figure 5: Distribution of Miocene and Pliocene sedimentary rocks (shaded) on the Eaton Alvarado ridge.These sediments lie in a chain of deep medial grabens (like in the Rio Grande rift in NM and southern CO) and in thin blankets on the flanking rises of the ridge.The region around the crest of the Alvarado ridge is outlined by a heavy solid line.The large continuous shaded area in the east, running from Nebraska down to the Texas Panhandle, is the Ogallala Formation (Figure 4 of [11]).

Figure 6 :
Figure 6: Geologic map of study area from Nelson et al. ([7], Figure 4) with the 800 km cross-section A-A  .

Figure 7 :Figure 8 :
Figure 7: Present-day geologic cross-section along line A-A  (Figure 6) with names of geohydrologic units.Modified from Nelson et al. [7].Vertical lines are both well locations and internal model boundaries.

Figure 9 :Figure 10 :Figure 11 :Figure 12 :
Figure 9: Map of potentiometric surface (hydraulic head contour lines in meters) in rocks of Wolfcampian age, after Nelson et al. ([7], Figure 14A).Hydraulic head ranges from less than 366 m in water wells in northeastern Kansas to greater than 610 m in oil and gas wells in eastern Colorado and southwestern Nebraska.Dots represent well control.Control data for the hydrologic model were taken along 800 km line A-A  .

( 3 ) 8 )
Run the model sequence once (4b) If h difference is as small as possible (4a) If h difference can be reduced (0) Creating and calibrating a hydrologic model (1) Uplift, deposition, erosion (4) Compare with measured {h} Appendix C: Permeability and Hydraulic Conductivity Values from Other Sources (Use Darcy's equation in (5a) and (5b) to see effect of trial K and Q values on dh/dx (2a) Choose 1st guess of K ℎ , K  , S s values (2b) Choose 1st guess of flow rates in W, Q w

Figure 14 :Figure 15 : 9 ElevationFigure 16 :Figure 17 :
Figure 14: Grid of the first model of the sequence, having the geometry at 27 Ma, used to model the period from 40 to 27 Ma.

/ 10 ( 7 ) 10
Wolfcampian-choke from 407 to 482 km 5 × 10 −4 1 × 10 −4 5 × 10 −4  ℎ Virgilian-Missourian-choke from 375 to 393 km 5 × 10 −4 1 × 10 −4 5 × 10 −4  ℎ /10 (8) Desmoinesian Low K west of 315 km 1 × 10 −9 west of 315 km 2 × 10 −10 High K east of 315 km 5 × 10 −3 east of 315 km 2 × 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 {0 from the initial heads in the aquifer at 70 Ma} + {Linearly decreasing values from 250 m in the west to 50 m in the east}.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.GeofluidsHydraulic conductivity (m/day) Sequence of 9 MODFLOW Models.When running the first of the 9 MODFLOW models,

Figure 19 :
Figure19: Permeability values selected by Senger and Fogg[31] for hydrologic units in a study of the Palo Duro Basin, Texas.Values range over seven orders of magnitude, from one darcy (D) to 0.1 microdarcies (uD).The left axis label shows equivalents for hydraulic conductivity in m d −1 using a conversion factor of 1 D = 0.835 m d −1 .STP, standard temperature and pressure; Perm-Penn, Permian-Pennsylvanian.The ratio of vertical to horizontal permeability  V / ℎ ranges from 0.01 to 1.0.

Figure 20 :Figure 21 :Head
Figure 20: Head versus time for Wolfcampian cells in columns 6, 11, and 23 of the 27 Ma model, showing that steady state is reached around 34.5 Ma, after 5.5 million years of simulated time.Head in meters above sea level for Figures 20-33.

Figure 22 :
Figure 22: Sensitivity to decreasing the flux in the west by one-half.

Figure 23 :Figure 24 :
Figure 23: Sensitivity to increasing  by half an order of magnitude in aquifers.
Sensitivity to increasing S s by 1 order of magnitude

Figure 25 :
Figure 25: Sensitivity to increasing   by an order of magnitude.

Figure 26 :
Figure 26: Sensitivity to decreasing   by an order of magnitude.
). Doubling the flux produces perturbed ℎ values greater than the base case, and halving the flux produces perturbed ℎ values less than the base case.At  = 200 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 pinchout and lack of direct hydraulic connection to flux at the

Figure 28 :Figure 29 :
Figure 28: Modeled head results of running the first model, the 27 Ma model representing the 13-million-year period from 40 to 27 Ma, compared to present-day heads.Model curves (1), (2), and (3) (Wolfcampian, Desmoinesian, and Mississippian, resp.)eventually approach the present-day data curves.Red triangles at centroids of the grid columns mark the initial zero heads at 70 Ma.

Figure 30 :
Figure 30: Results of running the fifth model, the 4.5 Ma model, representing the 2.5-million-year period from 7 to 4.5 Ma, compared to present-day heads.Model curves (1), (2), and (3) (Wolfcampian, Desmoinesian, and Mississippian, resp.) are further approaching the present-day data curves.Red triangles at centroids of the grid columns mark the initial zero heads at 70 Ma.

Figure 31 :Figure 32 :
Figure 31: Results from the 0 Ma model, representing the 2million-year period from 2 to 0 Ma and present-day head values.Model curves (1), (2), and (3) (Wolfcampian, Desmoinesian, and Mississippian, resp.)approximate the present-day data curves.Red triangles at centroids of the grid columns mark the initial zero heads at 70 Ma.

Figure 33 :
Figure 33: Head versus time for Wolfcampian cells in columns 6, 11, and 23 of the 0 Ma model, representing the period from 2 to 0 Ma.At 0 Ma, the solution has equilibrated at the eastern part of the crosssection but has not reached steady state in the central and western parts of the model.

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.
t o r o f z e r o s : 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 s c r i p t i n d i c a t i n g c e l l 1 c2: S u b s c r i p t i n d i c a t i n g c e l l 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:

Table 1 :
Correlation of sections of the paper with FE elements of the flow chart of the modeling process in Figure13that they refer to.

Table 3 :
The sequence of nine MODFLOW Models (m.y.: million years).

Table 4 :
Percent and cumulative percent of total uplift, relative to 27 Ma, for the 9 MODFLOW model times.

Table 5 :
Hydraulic conductivities of the (10) geohydrologic units (layers) of each model in the sequence.Hydraulic conductivities in meters per day (m d −1 [16]or et al.[16], Midcontinent, 1996 Figure 18: Hydraulic conductivity values used by Signor et al. [16].Permeability values in square feet from four maps by Signor et al. [16] for hydrologic units in a study of the Great Plains of the United States.Mapped values increase from west to east, indicated by the open bars.Values between the Kansas/Colorado border and the Nemaha uplift lie within the solid (shaded) bars.Names of hydrologic units used by Signor et al. appear vertically in italics; equivalent names used in this paper appear vertically in normal font.
Scale for hydraulic conductivity units (m d −1 ) used in this paper appears at left, using a conversion of 1 m d −1 = 1.19 × 10 −12 m 2 .Values given by Signor and others are permeability in ft-squared, taken from maps (Figures 21, 22, 29, and 36 in Signor et al. [16]).Porosity values from geophysical logs were used to compute permeability estimates using the Kozeny-Carman equation.Conversion used is  (m/d) = 1.28E-11 × k (ft-squared), at standard temperature and pressure.

Table 6 :
[17]ytical solution for decay time to steady state in the west as a function of hydraulic conductivity, based on equations in Carslaw and Jaeger[17].