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

Authors: Amjad M. J. Umari; Philip H. Nelson; Chris Fridrich; Gary D. LeCain
Journal: Geofluids
 (2018)

---

## 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.


---

## Body

## 1. Introduction

The existence of underpressures (subhydrostatic hydraulic head values, AppendixA), 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 (Figure1 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.Figure 1
The definition of underpressure.D, depth of measuring point; WT, elevation of water table; GL, elevation of ground level; Up_real, actual underpressure as head drop; Up_approx, head drop from land surface; p, pressure; γ, unit weight of water; Z, elevation of formation.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
′.A hydrodynamic model is deemed appropriate for the confined aquifers of the Great Plains. For example, Belitz and Bredehoeft [3], p. 1356, postulate thatGenerally, 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 h with distance x, elevation y, and time t:(1)
∂
2
h
∂
x
2
+
∂
2
h
∂
y
2
=
S
s
K
∂
h
∂
t
,where K is hydraulic conductivity in m s−1 and S
s 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 (n
=
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 n to model n
+
1, mainly elevation and tilt as calculated in a spreadsheet, are summed with the head values h computed for model n. These augmented values of h then serve as starting values for model n
+
1. This procedure is repeated until n
=
9 and a solution for h 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. AppendixB 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].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]).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
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]).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 (Figure5). 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 lineA
-
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. Uplift during the period of 40 to 27 Ma was assumed to be zero (Appendix B).Figure 6
Geologic map of study area from Nelson et al. ([7], Figure 4) with the 800 km cross-section A
-
A
′.Figure 7
Present-day geologic cross-section along lineA
-
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 8
Estimated cumulative uplift and percent of total uplift as a function of time since 27 Ma, with model numbers and superimposed periods of deposition and erosion.


### 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 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
′.Figure 10
Map of potentiometric surface (hydraulic head contour lines in meters) in rocks of Missourian age, after Nelson et al. ([7], Figure  14C). Hydraulic head ranges from less than 244 m in water wells in Missouri and Oklahoma to greater than 732 m in oil and gas wells in eastern Colorado. Dots represent well control. Control data for the hydrologic model were taken along 800 km line A
-
A
′.Figure 11
Map of potentiometric surface (hydraulic head contour lines in meters) in rocks of Mississippian age, after Nelson et al. ([7], Figure  14E). Hydraulic head decreases from about 366 m in water wells in Missouri and Arkansas to less than 183 m in oil and gas wells in eastern Oklahoma and then increases westward to values greater than 366 m in eastern Colorado. Dots represent well control. Control data for the hydrologic model were taken along 800 km line A
-
A
′.Figure12 presents the hydraulic head values along cross-section 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.Figure 12
Hydraulic head versus distance digitized from head contours. Hydraulic head for six geohydrologic units along lineA
-
A
′, which serve as control data for the hydrologic model of this paper. Head values were digitized from potentiometric maps such as those shown in Figures 9, 10, and 11.


### 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 Figure13.Figure 13
Flow chart of the modeling process.h
E: head specified in east; L
E: Land surface elevation in the east; h
i: interpolated h between h
4.5
M
a and Wolfcampian h contour in E for layers (
1
) to (
7
), and between h
4.5
M
a and Mississippian h contour in E for layers (
8
) to (
10
). Box numbers in the figure after step (
4
) are for labelling and are not strictly sequential.Each flow chart element # is referred to herein as FE # and is related to a section of the paper in Table1.Table 1
Correlation of sections of the paper with FE elements of the flow chart of the modeling process in Figure13 that they refer to.
Section title
Flow chart of the modeling process Figure13, Element #
Abstract
1. Introduction
2. The Geologic Setting
3. Methodology
3.1. Understanding the Geology and Quantifying Structural Changes
3.2. Assembling the present-day geofluid underpressure data
3.3. Creating and Calibrating the Hydrologic Model to Simulate
FE#’s(
3
a
)
,
(
3
b
)
,
…
,
(
3
i
); (4b), (9), (2b), (4a), (5b); (2a),
Structurally Driven Underpressures
(4a), (5a), (7b), (7a)
3.3.1. Uplift, Deposition, And Erosion
FE#(
1
)
3.3.2. Hydraulic Conductivity and Specific Storage
FE# (6), AppendixC, (2a), (8), (5a), (4b), (7a)
Darcy’s equation (2)
FE#’s: (4a), (5a), (4b), (7a)
3.3.3. Boundary Conditions (BCs) for the Sequence of
FE#’s: (2b), (4a), (8), (5b), (4b), (7a), (7b); FE’s:h
E
=
L
E,
Models
h
E
=
h
i arrows of FE# (3)
3.3.4. Initial Heads for the Sequence of Models
FE#’s (3a), (3b), (3i);{
h
i
}’s under FE#’s
(
3
a
)
,
(
3
b
)
,
…
,
(
3
i
)
3.3.5. Transience While Running the Sequence of 9 MODFLOW
FE# (3), FE# (3a)
Models
3.3.6. Sensitivity Analysis
3.3.7. Model Limitations
4. Results and Discussion
5. Summary and Conclusions
A. Theoretical Basis of Underpressure
B. Expanded Geologic Setting
C. Permeability and Hydraulic Conductivity Values from Other Sources
Definitions of Terms
Acknowledgements
ReferencesPart 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 Figure13; 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 z direction (Table 2), 26 “columns” spanning the horizontal x 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 K 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 K 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 h 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.Table 2
Stratigraphic column.
Age
Geohydrologic unit name
Lithologic description
Geohydrologic unit type
Model layer
Range of thickness (feet)
Neogene
Ogallala Formation
Terrestrial gravel, sand, and clay
Aquifer
1
0–300
Neogene-Paleogene
Laramide Basin Fill
Terrestrial gravel, sand, and clay
Aquifer
2
0–2,300
Late Cretaceous
Pierre & Niobrara Formations
Marine shales and carbonates, some sand layers
Aquitard
3
0–6,800
Early Cretaceous
Dakota Formation
Sandstone
Aquifer
4
0–500
Jurassic-Triassic-late Permian
Morrison-Leonardian Formations
Terrestrial sand, silt, and clay with upper Permian shale and evaporites
Aquitard
5
0–1,200
Early Permian
Wolfcampian Formation
Clastic marine sediments with carbonate beds
Aquifer
6
0–1,100
Pennsylvanian
Virgilian-Missourian Formation
Clastic marine sediments with carbonate beds
Aquifer
7
400–1,600
Pennsylvanian
Desmoinesian Formation
Marine carbonates in the east, terrestrial sediments in the west
Aquitard in west, aquifer in east
8
100–700
Early Pennsylvanian
Basal Pennsylvanian Unit
Marine shales with some carbonates
Aquitard
9
100
Early Paleozoic
Mississippian-Silurian-Ordovician-Cambrian Units
Limestones and dolomites with karst weathering, & lesser interbedded shales
Aquifer
10
0–1,100
Precambrian
Precambrian Basement
Granite and gneiss
No-flow boundary
Model lower boundary
NATable 3
The sequence of nine MODFLOW Models (m.y.: million years).
Sequence number of each model
Name of model
Time span (Ma)
Number of years modeled (m.y.)
Figure number showing grid
First
27 Ma
40 to 27 Ma
13
Figure14
Second
18 Ma
27 to 18 Ma
9
Figure15
Third
11 Ma
18 to 11 Ma
7
None
Fourth
7 Ma
11 to 7 Ma
4
None
Fifth
4.5 Ma
7 to 4.5 Ma
2.5
Figure16
Sixth
3.5 Ma
4.5 to 3.5 Ma
1
None
Seventh
3 Ma
3.5 to 3 Ma
0.5
None
Eighth
2 Ma
3 to 2 Ma
1
None
Ninth
0 Ma
2 to 0 Ma
2
Figure7Figure 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.Figure 15
Grid of the second model of the sequence, having the geometry at 18 Ma, used to model the period from 27 to 18 Ma.Figure 16
Grid of the fifth model of the sequence, having the geometry at 4.5 Ma, of maximum deposition, used to model the period from 7 to 4.5 Ma. The vertically expanded part of section from 200 to 600 km in the inset shows details of horizontal hydraulic conductivity structure.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 #’s3
a
,
3
b
,
…
,
3
i 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 (Table3) 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 (
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).Figure 17
Stacked land surfaces of the 9 MODFLOW models showing progressive uplift from 27 Ma to 0 Ma.


#### 3.3.1. Uplift, Deposition, and Erosion

Uplift, deposition, and erosion were incorporated into each of the sequences of MODFLOW models. Table4 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.Table 4
Percent and cumulative percent of total uplift, relative to 27 Ma, for the 9 MODFLOW model times.
Geologic time
Cumulative percent of total uplift at mountain front (assumed constant for allx) (CpTU)
Percent of total uplift (PTU)
27 Ma
0
13.125
18 Ma
13.125
13.125
11 Ma
26.250
11.250
7 Ma
37.500
15.625
4.5 Ma
53.125
15.625
3.5 Ma
68.750
9.375
3 Ma
78.125
12.500
2 Ma
90.625
9.375
0
100Figure7 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 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 Figure14, 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 (Figure8). It was further assumed that the maximum deposition, at 4.5 Ma, was 183 m atx
=
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 atx
=
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 Figure16. 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 ofx, 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 (S
s), horizontal hydraulic conductivity (
K
h
), and vertical hydraulic conductivity (
K
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.


#### 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 Figure12 are highly sensitive to hydraulic conductivity K, 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 K.For further adjustment ofK values, a useful tool is provided by Darcy’s equation relating flux q to the hydraulic conductivity K and the hydraulic gradient, d
h
/
d
x, which aided us during calibration:(2)
q
h
=
-
K
h
·
d
h
d
x
,where K
h is the horizontal hydraulic conductivity and h is hydraulic head. Equation (2) indicates that, for a specific q
h, K
h and d
h
/
d
x are inversely proportional. That is, as the medium becomes more conductive and K increases, the hydraulic gradient decreases and the plot of h versus x becomes shallower or less steep. Conversely, a steep hydraulic gradient indicates a tight or low hydraulic conductivity medium. Based on this, the K
h and K
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 (Figure12) 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 K
v
=
K
h
/
10 in most of the geohydrologic units (Table 5), a selection that is the default setting in MODFLOW. However, K
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, K
v, is large, the hydraulic head plots of h 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 K
v, is needed.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). In units (1)–(7), the vertical hydraulic conductivity K
v is set to one-tenth the value of the horizontal hydraulic conductivity K
h. In units (6) and (7), a low-K choke separates more permeable units in the west and east. In unit (8), K is higher in the east than in the west.
Geohydrologic unit
K
h (m d−1)
K
v (m d−1)
(
1
) Ogallala
1 × 10−6
K
h/10
(
2
) Tertiary sediments
1 × 10−6
K
h/10
(3) Pierre
1 × 10−6
K
h/10
(4) Dakota
1 × 10−6
K
h/10
(5) Morrison-Leonardian
1 × 10−6
K
h/10
West
Choke
East
(6) Wolfcampian—choke from 407 to 482 km
5 × 10−4
1 × 10−4
5 × 10−4
K
h/10
(7) Virgilian-Missourian—choke from 375 to 393 km
5 × 10−4
1 × 10−4
5 × 10−4
K
h/10
(8) Desmoinesian
LowK west of 315 km
1 × 10−9 west of 315 km
2 × 10−10
HighK east of 315 km
5 × 10−3 east of 315 km
2 × 10−10
(9) Basal Pennsylvanian
1 × 10−9
2 × 10−10
(10) Mississippian
5 × 10−3
2 × 10−10TheK
h 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, K
h 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 K values from the model calibration.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 m2. 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 K (m/d) = 1.28E-11 ×k (ft-squared), at standard temperature and pressure.Figure12 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 K
h of the Desmoinesian geohydrologic unit. The model assigns a K
h value of 1 × 10−9 m d−1 for the western section of the Desmoinesian and a K
h 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,K
h is one-fifth that of K
h 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).Figure 19
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 k
v
/
k
h ranges from 0.01 to 1.0.K
h 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,S
s, value of 1 × 10−6 m−1 was used throughout the model. The storage coefficient, S, was calculated as S
s× 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,Q
w, or flux, q
w, 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 (q) 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 (Figures2 and 6), where the confined units crop out and all layers are exposed to the atmosphere. From the definition of the hydraulic head, h, in (A.1), we see that when the pressure is atmospheric, p becomes 0 and h
=
z. 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 h
E
=
L
E arrows in Figure 13). The same constant head BC (h
=
z) 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 (Figures9 and 11; [7]); see h
E
=
h
i 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 toC
W
E 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, h
i, 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, ash
E
=
h
i, but with h
i using calculated and actual contoured h values of the Mississippian instead of the Wolfcampian, with a C
M
E 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,Q
w, S
s, and K (K
h and K
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 Q
w, S
s, and K values could have produced calculated present-day hydraulic head values as close, or closer, to the actual present-day hydraulic head values. The “optimal” Q
w, S
s, and K (K
h and K
v) values were used in obtaining the final modeling results presented in this paper.The sensitivity of the calculatedh values to assumed flux (q
w), S
s, and K (K
h and K
v) 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 (h), 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, h at the top is 0 because p
=
0 (no water above it) and z
=
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 D, p
/
γ
=
D and z
=
-
D; therefore, h
=
D
-
D
=
0. So, the initial hydraulic head h is 0 everywhere in the aquifer at 70 Ma and is taken to represent mean sea level. Referring to the h values at all the cells in the model symbolically as the “vector” {
h
}, we can say that h
i
=
{
0
}, where i means “initial.”The land surface for the first model, the 27 Ma model shown in Figure14 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), z will be increased by 250 m in the west, decreasing linearly to 50 m in the east. Therefore, the initial {
h
} 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 h 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{
h
} values at 40 Ma, as described above, and the boundary conditions and K and S parameters presented earlier, MODFLOW calculates {
h
} 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 h
27
M
a.For the 18 Ma model shown in Figure15, a sudden (instantaneous) additional rise, z, equal to the uplift from 27 to 18 Ma, U
p
l
i
f
t
27–18 Ma, takes place along the length of the model, at 27 Ma, creating an equal rise in h, according to (A.1). The incremented {
h
} values, h
27
M
a
+
U
p
l
i
f
t
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,h
27
M
a + U
p
l
i
f
t
27–18 Ma, at 27 Ma, and the boundary condition presented earlier, MODFLOW calculates {
h
} for each of the 9 million years from 27 to 18 Ma that the 18 Ma model represents, ending at 18 Ma with h
18
M
a. The above process is repeated until MODFLOW calculates h
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.


#### 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. Figure20 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 K 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.Table 6
Analytical solution for decay time to steady state in the west as a function of hydraulic conductivity, based on equations in Carslaw and Jaeger [17].
K (m d−1)
K (m s−1)
Decay time (million years)
1.00 × 10−2
1.16 × 10−7
0.149
1.00 × 10−3
1.16 × 10−8
1.49
1.00 × 10−4
1.16 × 10−9
14.9
1.00 × 10−5
1.16 × 10−10
149
1.00 × 10−6
1.16 × 10−11
1,490Table 7
Sensitivity analysis (OM: order of magnitude, Og: Ogallala, Lbf: Laramide Basin Fill, Dak: Dakota, Wolfc: Wolfcampian, Virg: Virgilian).
K
h in m/d (K
v = 0.1 K
h)
Q in m3/1,000 yr
S
s
h in Wolfcampian at 200 km, in m
h in Wolfcampian at 200 km, % change
See figure
Base case
474.5
1 × 10−6
119.90
0
Og, Lbf, Dak
1 × 10−6
Wolfc, Virg,
5 × 10−4
Desmoinesian,
5 × 10−3
Mississippian
5 × 10−3
1
/
2 OM increase in K of aquifers:
89.21
−25.6
Figure23
Wolfc, Virg,
1 × 10−3
Desmoinesian,
1 × 10−2
Mississippian
1 × 10−2
1
/
2 OM decrease in K of aquifers:
200.03
+66.83
Figure24
Wolfc, Virg,
1 × 10−4
Desmoinesian,
1 × 10−3
Mississippian
1 × 10−3
2 OM higherK in Og, Lbf, Dak
1 × 10−4
109.167
−8.95
Figure27
DoubleQ in west
949.0
152.05
+26.81
Figure21
HalveQ in west
237.3
96.73
−19.33
Figure22
1 OM higherS
s
1 × 10−5
200.67
+67.36
Figure25
1 OM lowerS
s
1 × 10−7
74.96
−37.48
Figure26Figure 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 Figures20–33.Figure 21
Sensitivity to doubling the flux in the west.Figure 22
Sensitivity to decreasing the flux in the west by one-half.Figure 23
Sensitivity to increasingK by half an order of magnitude in aquifers.Figure 24
Sensitivity to decreasingK by half an order of magnitude in aquifers.Figure 25
Sensitivity to increasingS
s by an order of magnitude.Figure 26
Sensitivity to decreasingS
s by an order of magnitude.Figure 27
Sensitivity to increasingK by 2 orders of magnitude in layers (
1
), (
2
), and (
4
).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 29
Results of running the second model, the 18 Ma model, representing the 9-million-year period from 27 to 18 Ma, compared to present-day heads. Model curves(
1
), (
2
), and (
3
) (Wolfcampian, Desmoinesian, and Mississippian, resp.) are approaching the present-day data curves. Red triangles at centroids of the grid columns mark the initial zero heads at 70 Ma.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
Results from the 0 Ma model, representing the 2-million-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 32
Progression of Wolfcampian head solutions from 27 to 0 Ma.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 cross-section but has not reached steady state in the central and western parts of the model.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 Table3), 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 Table5 for values of K), we computed perturbed h profiles corresponding to perturbed values of flux on the western boundary (Figures 21 and 22). Doubling the flux produces perturbed h values greater than the base case, and halving the flux produces perturbed h values less than the base case. At x
=
200 km, the changes in h range from +5 m in the Mississippian unit to +32 m in the Wolfcampian unit. The relatively low changes in h 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 x-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 parametersK and S
s (Figures 23–26). The key parameter in (1) is the ratio K
/
S
s, suggesting that perturbations in K and S
s should produce opposite effects, and this is indeed the case. An increase in K drops h below the base case (Figure 23), and so does a decrease in S
s (Figure 26). Conversely, a decrease in K and an increase in S
s both produce increases in h (Figures 24 and 25).A third sensitivity test examines the effect of a much higherK in the uppermost layers of the model (Figure 27). K of Pierre Shale is maintained at K = 10−6 and the overlying and underlying units are increased to 10−4. Figure 27 shows that if K 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, h 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 Table7.


#### 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 ofS
s and K (K
h and K
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, S
s, and K 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 thex direction covering a horizontal distance of 800 km and 10 layers in the z direction covering a vertical distance of 4 km, and aquifer properties (K and S
s) 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 horizontalx (length) and vertical z (depth) directions and that no flow takes place in the horizontal y direction.


## 3.1. Understanding the Geology and Quantifying Structural Changes

A geologic map of the study area and a cross-section along lineA
-
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. Uplift during the period of 40 to 27 Ma was assumed to be zero (Appendix B).Figure 6
Geologic map of study area from Nelson et al. ([7], Figure 4) with the 800 km cross-section A
-
A
′.Figure 7
Present-day geologic cross-section along lineA
-
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 8
Estimated cumulative uplift and percent of total uplift as a function of time since 27 Ma, with model numbers and superimposed periods of deposition and erosion.


## 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 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
′.Figure 10
Map of potentiometric surface (hydraulic head contour lines in meters) in rocks of Missourian age, after Nelson et al. ([7], Figure  14C). Hydraulic head ranges from less than 244 m in water wells in Missouri and Oklahoma to greater than 732 m in oil and gas wells in eastern Colorado. Dots represent well control. Control data for the hydrologic model were taken along 800 km line A
-
A
′.Figure 11
Map of potentiometric surface (hydraulic head contour lines in meters) in rocks of Mississippian age, after Nelson et al. ([7], Figure  14E). Hydraulic head decreases from about 366 m in water wells in Missouri and Arkansas to less than 183 m in oil and gas wells in eastern Oklahoma and then increases westward to values greater than 366 m in eastern Colorado. Dots represent well control. Control data for the hydrologic model were taken along 800 km line A
-
A
′.Figure12 presents the hydraulic head values along cross-section 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.Figure 12
Hydraulic head versus distance digitized from head contours. Hydraulic head for six geohydrologic units along lineA
-
A
′, which serve as control data for the hydrologic model of this paper. Head values were digitized from potentiometric maps such as those shown in Figures 9, 10, and 11.


## 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 Figure13.Figure 13
Flow chart of the modeling process.h
E: head specified in east; L
E: Land surface elevation in the east; h
i: interpolated h between h
4.5
M
a and Wolfcampian h contour in E for layers (
1
) to (
7
), and between h
4.5
M
a and Mississippian h contour in E for layers (
8
) to (
10
). Box numbers in the figure after step (
4
) are for labelling and are not strictly sequential.Each flow chart element # is referred to herein as FE # and is related to a section of the paper in Table1.Table 1
Correlation of sections of the paper with FE elements of the flow chart of the modeling process in Figure13 that they refer to.
Section title
Flow chart of the modeling process Figure13, Element #
Abstract
1. Introduction
2. The Geologic Setting
3. Methodology
3.1. Understanding the Geology and Quantifying Structural Changes
3.2. Assembling the present-day geofluid underpressure data
3.3. Creating and Calibrating the Hydrologic Model to Simulate
FE#’s(
3
a
)
,
(
3
b
)
,
…
,
(
3
i
); (4b), (9), (2b), (4a), (5b); (2a),
Structurally Driven Underpressures
(4a), (5a), (7b), (7a)
3.3.1. Uplift, Deposition, And Erosion
FE#(
1
)
3.3.2. Hydraulic Conductivity and Specific Storage
FE# (6), AppendixC, (2a), (8), (5a), (4b), (7a)
Darcy’s equation (2)
FE#’s: (4a), (5a), (4b), (7a)
3.3.3. Boundary Conditions (BCs) for the Sequence of
FE#’s: (2b), (4a), (8), (5b), (4b), (7a), (7b); FE’s:h
E
=
L
E,
Models
h
E
=
h
i arrows of FE# (3)
3.3.4. Initial Heads for the Sequence of Models
FE#’s (3a), (3b), (3i);{
h
i
}’s under FE#’s
(
3
a
)
,
(
3
b
)
,
…
,
(
3
i
)
3.3.5. Transience While Running the Sequence of 9 MODFLOW
FE# (3), FE# (3a)
Models
3.3.6. Sensitivity Analysis
3.3.7. Model Limitations
4. Results and Discussion
5. Summary and Conclusions
A. Theoretical Basis of Underpressure
B. Expanded Geologic Setting
C. Permeability and Hydraulic Conductivity Values from Other Sources
Definitions of Terms
Acknowledgements
ReferencesPart 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 Figure13; 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 z direction (Table 2), 26 “columns” spanning the horizontal x 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 K 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 K 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 h 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.Table 2
Stratigraphic column.
Age
Geohydrologic unit name
Lithologic description
Geohydrologic unit type
Model layer
Range of thickness (feet)
Neogene
Ogallala Formation
Terrestrial gravel, sand, and clay
Aquifer
1
0–300
Neogene-Paleogene
Laramide Basin Fill
Terrestrial gravel, sand, and clay
Aquifer
2
0–2,300
Late Cretaceous
Pierre & Niobrara Formations
Marine shales and carbonates, some sand layers
Aquitard
3
0–6,800
Early Cretaceous
Dakota Formation
Sandstone
Aquifer
4
0–500
Jurassic-Triassic-late Permian
Morrison-Leonardian Formations
Terrestrial sand, silt, and clay with upper Permian shale and evaporites
Aquitard
5
0–1,200
Early Permian
Wolfcampian Formation
Clastic marine sediments with carbonate beds
Aquifer
6
0–1,100
Pennsylvanian
Virgilian-Missourian Formation
Clastic marine sediments with carbonate beds
Aquifer
7
400–1,600
Pennsylvanian
Desmoinesian Formation
Marine carbonates in the east, terrestrial sediments in the west
Aquitard in west, aquifer in east
8
100–700
Early Pennsylvanian
Basal Pennsylvanian Unit
Marine shales with some carbonates
Aquitard
9
100
Early Paleozoic
Mississippian-Silurian-Ordovician-Cambrian Units
Limestones and dolomites with karst weathering, & lesser interbedded shales
Aquifer
10
0–1,100
Precambrian
Precambrian Basement
Granite and gneiss
No-flow boundary
Model lower boundary
NATable 3
The sequence of nine MODFLOW Models (m.y.: million years).
Sequence number of each model
Name of model
Time span (Ma)
Number of years modeled (m.y.)
Figure number showing grid
First
27 Ma
40 to 27 Ma
13
Figure14
Second
18 Ma
27 to 18 Ma
9
Figure15
Third
11 Ma
18 to 11 Ma
7
None
Fourth
7 Ma
11 to 7 Ma
4
None
Fifth
4.5 Ma
7 to 4.5 Ma
2.5
Figure16
Sixth
3.5 Ma
4.5 to 3.5 Ma
1
None
Seventh
3 Ma
3.5 to 3 Ma
0.5
None
Eighth
2 Ma
3 to 2 Ma
1
None
Ninth
0 Ma
2 to 0 Ma
2
Figure7Figure 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.Figure 15
Grid of the second model of the sequence, having the geometry at 18 Ma, used to model the period from 27 to 18 Ma.Figure 16
Grid of the fifth model of the sequence, having the geometry at 4.5 Ma, of maximum deposition, used to model the period from 7 to 4.5 Ma. The vertically expanded part of section from 200 to 600 km in the inset shows details of horizontal hydraulic conductivity structure.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 #’s3
a
,
3
b
,
…
,
3
i 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 (Table3) 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 (
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).Figure 17
Stacked land surfaces of the 9 MODFLOW models showing progressive uplift from 27 Ma to 0 Ma.


### 3.3.1. Uplift, Deposition, and Erosion

Uplift, deposition, and erosion were incorporated into each of the sequences of MODFLOW models. Table4 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.Table 4
Percent and cumulative percent of total uplift, relative to 27 Ma, for the 9 MODFLOW model times.
Geologic time
Cumulative percent of total uplift at mountain front (assumed constant for allx) (CpTU)
Percent of total uplift (PTU)
27 Ma
0
13.125
18 Ma
13.125
13.125
11 Ma
26.250
11.250
7 Ma
37.500
15.625
4.5 Ma
53.125
15.625
3.5 Ma
68.750
9.375
3 Ma
78.125
12.500
2 Ma
90.625
9.375
0
100Figure7 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 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 Figure14, 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 (Figure8). It was further assumed that the maximum deposition, at 4.5 Ma, was 183 m atx
=
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 atx
=
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 Figure16. 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 ofx, 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 (S
s), horizontal hydraulic conductivity (
K
h
), and vertical hydraulic conductivity (
K
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.


### 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 Figure12 are highly sensitive to hydraulic conductivity K, 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 K.For further adjustment ofK values, a useful tool is provided by Darcy’s equation relating flux q to the hydraulic conductivity K and the hydraulic gradient, d
h
/
d
x, which aided us during calibration:(2)
q
h
=
-
K
h
·
d
h
d
x
,where K
h is the horizontal hydraulic conductivity and h is hydraulic head. Equation (2) indicates that, for a specific q
h, K
h and d
h
/
d
x are inversely proportional. That is, as the medium becomes more conductive and K increases, the hydraulic gradient decreases and the plot of h versus x becomes shallower or less steep. Conversely, a steep hydraulic gradient indicates a tight or low hydraulic conductivity medium. Based on this, the K
h and K
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 (Figure12) 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 K
v
=
K
h
/
10 in most of the geohydrologic units (Table 5), a selection that is the default setting in MODFLOW. However, K
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, K
v, is large, the hydraulic head plots of h 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 K
v, is needed.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). In units (1)–(7), the vertical hydraulic conductivity K
v is set to one-tenth the value of the horizontal hydraulic conductivity K
h. In units (6) and (7), a low-K choke separates more permeable units in the west and east. In unit (8), K is higher in the east than in the west.
Geohydrologic unit
K
h (m d−1)
K
v (m d−1)
(
1
) Ogallala
1 × 10−6
K
h/10
(
2
) Tertiary sediments
1 × 10−6
K
h/10
(3) Pierre
1 × 10−6
K
h/10
(4) Dakota
1 × 10−6
K
h/10
(5) Morrison-Leonardian
1 × 10−6
K
h/10
West
Choke
East
(6) Wolfcampian—choke from 407 to 482 km
5 × 10−4
1 × 10−4
5 × 10−4
K
h/10
(7) Virgilian-Missourian—choke from 375 to 393 km
5 × 10−4
1 × 10−4
5 × 10−4
K
h/10
(8) Desmoinesian
LowK west of 315 km
1 × 10−9 west of 315 km
2 × 10−10
HighK east of 315 km
5 × 10−3 east of 315 km
2 × 10−10
(9) Basal Pennsylvanian
1 × 10−9
2 × 10−10
(10) Mississippian
5 × 10−3
2 × 10−10TheK
h 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, K
h 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 K values from the model calibration.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 m2. 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 K (m/d) = 1.28E-11 ×k (ft-squared), at standard temperature and pressure.Figure12 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 K
h of the Desmoinesian geohydrologic unit. The model assigns a K
h value of 1 × 10−9 m d−1 for the western section of the Desmoinesian and a K
h 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,K
h is one-fifth that of K
h 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).Figure 19
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 k
v
/
k
h ranges from 0.01 to 1.0.K
h 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,S
s, value of 1 × 10−6 m−1 was used throughout the model. The storage coefficient, S, was calculated as S
s× 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,Q
w, or flux, q
w, 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 (q) 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 (Figures2 and 6), where the confined units crop out and all layers are exposed to the atmosphere. From the definition of the hydraulic head, h, in (A.1), we see that when the pressure is atmospheric, p becomes 0 and h
=
z. 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 h
E
=
L
E arrows in Figure 13). The same constant head BC (h
=
z) 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 (Figures9 and 11; [7]); see h
E
=
h
i 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 toC
W
E 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, h
i, 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, ash
E
=
h
i, but with h
i using calculated and actual contoured h values of the Mississippian instead of the Wolfcampian, with a C
M
E 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,Q
w, S
s, and K (K
h and K
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 Q
w, S
s, and K values could have produced calculated present-day hydraulic head values as close, or closer, to the actual present-day hydraulic head values. The “optimal” Q
w, S
s, and K (K
h and K
v) values were used in obtaining the final modeling results presented in this paper.The sensitivity of the calculatedh values to assumed flux (q
w), S
s, and K (K
h and K
v) 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 (h), 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, h at the top is 0 because p
=
0 (no water above it) and z
=
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 D, p
/
γ
=
D and z
=
-
D; therefore, h
=
D
-
D
=
0. So, the initial hydraulic head h is 0 everywhere in the aquifer at 70 Ma and is taken to represent mean sea level. Referring to the h values at all the cells in the model symbolically as the “vector” {
h
}, we can say that h
i
=
{
0
}, where i means “initial.”The land surface for the first model, the 27 Ma model shown in Figure14 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), z will be increased by 250 m in the west, decreasing linearly to 50 m in the east. Therefore, the initial {
h
} 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 h 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{
h
} values at 40 Ma, as described above, and the boundary conditions and K and S parameters presented earlier, MODFLOW calculates {
h
} 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 h
27
M
a.For the 18 Ma model shown in Figure15, a sudden (instantaneous) additional rise, z, equal to the uplift from 27 to 18 Ma, U
p
l
i
f
t
27–18 Ma, takes place along the length of the model, at 27 Ma, creating an equal rise in h, according to (A.1). The incremented {
h
} values, h
27
M
a
+
U
p
l
i
f
t
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,h
27
M
a + U
p
l
i
f
t
27–18 Ma, at 27 Ma, and the boundary condition presented earlier, MODFLOW calculates {
h
} for each of the 9 million years from 27 to 18 Ma that the 18 Ma model represents, ending at 18 Ma with h
18
M
a. The above process is repeated until MODFLOW calculates h
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.


### 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. Figure20 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 K 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.Table 6
Analytical solution for decay time to steady state in the west as a function of hydraulic conductivity, based on equations in Carslaw and Jaeger [17].
K (m d−1)
K (m s−1)
Decay time (million years)
1.00 × 10−2
1.16 × 10−7
0.149
1.00 × 10−3
1.16 × 10−8
1.49
1.00 × 10−4
1.16 × 10−9
14.9
1.00 × 10−5
1.16 × 10−10
149
1.00 × 10−6
1.16 × 10−11
1,490Table 7
Sensitivity analysis (OM: order of magnitude, Og: Ogallala, Lbf: Laramide Basin Fill, Dak: Dakota, Wolfc: Wolfcampian, Virg: Virgilian).
K
h in m/d (K
v = 0.1 K
h)
Q in m3/1,000 yr
S
s
h in Wolfcampian at 200 km, in m
h in Wolfcampian at 200 km, % change
See figure
Base case
474.5
1 × 10−6
119.90
0
Og, Lbf, Dak
1 × 10−6
Wolfc, Virg,
5 × 10−4
Desmoinesian,
5 × 10−3
Mississippian
5 × 10−3
1
/
2 OM increase in K of aquifers:
89.21
−25.6
Figure23
Wolfc, Virg,
1 × 10−3
Desmoinesian,
1 × 10−2
Mississippian
1 × 10−2
1
/
2 OM decrease in K of aquifers:
200.03
+66.83
Figure24
Wolfc, Virg,
1 × 10−4
Desmoinesian,
1 × 10−3
Mississippian
1 × 10−3
2 OM higherK in Og, Lbf, Dak
1 × 10−4
109.167
−8.95
Figure27
DoubleQ in west
949.0
152.05
+26.81
Figure21
HalveQ in west
237.3
96.73
−19.33
Figure22
1 OM higherS
s
1 × 10−5
200.67
+67.36
Figure25
1 OM lowerS
s
1 × 10−7
74.96
−37.48
Figure26Figure 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 Figures20–33.Figure 21
Sensitivity to doubling the flux in the west.Figure 22
Sensitivity to decreasing the flux in the west by one-half.Figure 23
Sensitivity to increasingK by half an order of magnitude in aquifers.Figure 24
Sensitivity to decreasingK by half an order of magnitude in aquifers.Figure 25
Sensitivity to increasingS
s by an order of magnitude.Figure 26
Sensitivity to decreasingS
s by an order of magnitude.Figure 27
Sensitivity to increasingK by 2 orders of magnitude in layers (
1
), (
2
), and (
4
).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 29
Results of running the second model, the 18 Ma model, representing the 9-million-year period from 27 to 18 Ma, compared to present-day heads. Model curves(
1
), (
2
), and (
3
) (Wolfcampian, Desmoinesian, and Mississippian, resp.) are approaching the present-day data curves. Red triangles at centroids of the grid columns mark the initial zero heads at 70 Ma.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
Results from the 0 Ma model, representing the 2-million-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 32
Progression of Wolfcampian head solutions from 27 to 0 Ma.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 cross-section but has not reached steady state in the central and western parts of the model.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 Table3), 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 Table5 for values of K), we computed perturbed h profiles corresponding to perturbed values of flux on the western boundary (Figures 21 and 22). Doubling the flux produces perturbed h values greater than the base case, and halving the flux produces perturbed h values less than the base case. At x
=
200 km, the changes in h range from +5 m in the Mississippian unit to +32 m in the Wolfcampian unit. The relatively low changes in h 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 x-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 parametersK and S
s (Figures 23–26). The key parameter in (1) is the ratio K
/
S
s, suggesting that perturbations in K and S
s should produce opposite effects, and this is indeed the case. An increase in K drops h below the base case (Figure 23), and so does a decrease in S
s (Figure 26). Conversely, a decrease in K and an increase in S
s both produce increases in h (Figures 24 and 25).A third sensitivity test examines the effect of a much higherK in the uppermost layers of the model (Figure 27). K of Pierre Shale is maintained at K = 10−6 and the overlying and underlying units are increased to 10−4. Figure 27 shows that if K 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, h 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 Table7.


### 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 ofS
s and K (K
h and K
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, S
s, and K 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 thex direction covering a horizontal distance of 800 km and 10 layers in the z direction covering a vertical distance of 4 km, and aquifer properties (K and S
s) 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 horizontalx (length) and vertical z (depth) directions and that no flow takes place in the horizontal y direction.


## 3.3.1. Uplift, Deposition, and Erosion

Uplift, deposition, and erosion were incorporated into each of the sequences of MODFLOW models. Table4 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.Table 4
Percent and cumulative percent of total uplift, relative to 27 Ma, for the 9 MODFLOW model times.
Geologic time
Cumulative percent of total uplift at mountain front (assumed constant for allx) (CpTU)
Percent of total uplift (PTU)
27 Ma
0
13.125
18 Ma
13.125
13.125
11 Ma
26.250
11.250
7 Ma
37.500
15.625
4.5 Ma
53.125
15.625
3.5 Ma
68.750
9.375
3 Ma
78.125
12.500
2 Ma
90.625
9.375
0
100Figure7 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 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 Figure14, 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 (Figure8). It was further assumed that the maximum deposition, at 4.5 Ma, was 183 m atx
=
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 atx
=
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 Figure16. 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 ofx, 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 (S
s), horizontal hydraulic conductivity (
K
h
), and vertical hydraulic conductivity (
K
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.


## 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 Figure12 are highly sensitive to hydraulic conductivity K, 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 K.For further adjustment ofK values, a useful tool is provided by Darcy’s equation relating flux q to the hydraulic conductivity K and the hydraulic gradient, d
h
/
d
x, which aided us during calibration:(2)
q
h
=
-
K
h
·
d
h
d
x
,where K
h is the horizontal hydraulic conductivity and h is hydraulic head. Equation (2) indicates that, for a specific q
h, K
h and d
h
/
d
x are inversely proportional. That is, as the medium becomes more conductive and K increases, the hydraulic gradient decreases and the plot of h versus x becomes shallower or less steep. Conversely, a steep hydraulic gradient indicates a tight or low hydraulic conductivity medium. Based on this, the K
h and K
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 (Figure12) 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 K
v
=
K
h
/
10 in most of the geohydrologic units (Table 5), a selection that is the default setting in MODFLOW. However, K
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, K
v, is large, the hydraulic head plots of h 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 K
v, is needed.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). In units (1)–(7), the vertical hydraulic conductivity K
v is set to one-tenth the value of the horizontal hydraulic conductivity K
h. In units (6) and (7), a low-K choke separates more permeable units in the west and east. In unit (8), K is higher in the east than in the west.
Geohydrologic unit
K
h (m d−1)
K
v (m d−1)
(
1
) Ogallala
1 × 10−6
K
h/10
(
2
) Tertiary sediments
1 × 10−6
K
h/10
(3) Pierre
1 × 10−6
K
h/10
(4) Dakota
1 × 10−6
K
h/10
(5) Morrison-Leonardian
1 × 10−6
K
h/10
West
Choke
East
(6) Wolfcampian—choke from 407 to 482 km
5 × 10−4
1 × 10−4
5 × 10−4
K
h/10
(7) Virgilian-Missourian—choke from 375 to 393 km
5 × 10−4
1 × 10−4
5 × 10−4
K
h/10
(8) Desmoinesian
LowK west of 315 km
1 × 10−9 west of 315 km
2 × 10−10
HighK east of 315 km
5 × 10−3 east of 315 km
2 × 10−10
(9) Basal Pennsylvanian
1 × 10−9
2 × 10−10
(10) Mississippian
5 × 10−3
2 × 10−10TheK
h 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, K
h 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 K values from the model calibration.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 m2. 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 K (m/d) = 1.28E-11 ×k (ft-squared), at standard temperature and pressure.Figure12 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 K
h of the Desmoinesian geohydrologic unit. The model assigns a K
h value of 1 × 10−9 m d−1 for the western section of the Desmoinesian and a K
h 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,K
h is one-fifth that of K
h 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).Figure 19
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 k
v
/
k
h ranges from 0.01 to 1.0.K
h 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,S
s, value of 1 × 10−6 m−1 was used throughout the model. The storage coefficient, S, was calculated as S
s× 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,Q
w, or flux, q
w, 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 (q) 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 (Figures2 and 6), where the confined units crop out and all layers are exposed to the atmosphere. From the definition of the hydraulic head, h, in (A.1), we see that when the pressure is atmospheric, p becomes 0 and h
=
z. 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 h
E
=
L
E arrows in Figure 13). The same constant head BC (h
=
z) 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 (Figures9 and 11; [7]); see h
E
=
h
i 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 toC
W
E 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, h
i, 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, ash
E
=
h
i, but with h
i using calculated and actual contoured h values of the Mississippian instead of the Wolfcampian, with a C
M
E 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,Q
w, S
s, and K (K
h and K
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 Q
w, S
s, and K values could have produced calculated present-day hydraulic head values as close, or closer, to the actual present-day hydraulic head values. The “optimal” Q
w, S
s, and K (K
h and K
v) values were used in obtaining the final modeling results presented in this paper.The sensitivity of the calculatedh values to assumed flux (q
w), S
s, and K (K
h and K
v) 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 (h), 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, h at the top is 0 because p
=
0 (no water above it) and z
=
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 D, p
/
γ
=
D and z
=
-
D; therefore, h
=
D
-
D
=
0. So, the initial hydraulic head h is 0 everywhere in the aquifer at 70 Ma and is taken to represent mean sea level. Referring to the h values at all the cells in the model symbolically as the “vector” {
h
}, we can say that h
i
=
{
0
}, where i means “initial.”The land surface for the first model, the 27 Ma model shown in Figure14 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), z will be increased by 250 m in the west, decreasing linearly to 50 m in the east. Therefore, the initial {
h
} 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 h 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{
h
} values at 40 Ma, as described above, and the boundary conditions and K and S parameters presented earlier, MODFLOW calculates {
h
} 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 h
27
M
a.For the 18 Ma model shown in Figure15, a sudden (instantaneous) additional rise, z, equal to the uplift from 27 to 18 Ma, U
p
l
i
f
t
27–18 Ma, takes place along the length of the model, at 27 Ma, creating an equal rise in h, according to (A.1). The incremented {
h
} values, h
27
M
a
+
U
p
l
i
f
t
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,h
27
M
a + U
p
l
i
f
t
27–18 Ma, at 27 Ma, and the boundary condition presented earlier, MODFLOW calculates {
h
} for each of the 9 million years from 27 to 18 Ma that the 18 Ma model represents, ending at 18 Ma with h
18
M
a. The above process is repeated until MODFLOW calculates h
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.


## 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. Figure20 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 K 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.Table 6
Analytical solution for decay time to steady state in the west as a function of hydraulic conductivity, based on equations in Carslaw and Jaeger [17].
K (m d−1)
K (m s−1)
Decay time (million years)
1.00 × 10−2
1.16 × 10−7
0.149
1.00 × 10−3
1.16 × 10−8
1.49
1.00 × 10−4
1.16 × 10−9
14.9
1.00 × 10−5
1.16 × 10−10
149
1.00 × 10−6
1.16 × 10−11
1,490Table 7
Sensitivity analysis (OM: order of magnitude, Og: Ogallala, Lbf: Laramide Basin Fill, Dak: Dakota, Wolfc: Wolfcampian, Virg: Virgilian).
K
h in m/d (K
v = 0.1 K
h)
Q in m3/1,000 yr
S
s
h in Wolfcampian at 200 km, in m
h in Wolfcampian at 200 km, % change
See figure
Base case
474.5
1 × 10−6
119.90
0
Og, Lbf, Dak
1 × 10−6
Wolfc, Virg,
5 × 10−4
Desmoinesian,
5 × 10−3
Mississippian
5 × 10−3
1
/
2 OM increase in K of aquifers:
89.21
−25.6
Figure23
Wolfc, Virg,
1 × 10−3
Desmoinesian,
1 × 10−2
Mississippian
1 × 10−2
1
/
2 OM decrease in K of aquifers:
200.03
+66.83
Figure24
Wolfc, Virg,
1 × 10−4
Desmoinesian,
1 × 10−3
Mississippian
1 × 10−3
2 OM higherK in Og, Lbf, Dak
1 × 10−4
109.167
−8.95
Figure27
DoubleQ in west
949.0
152.05
+26.81
Figure21
HalveQ in west
237.3
96.73
−19.33
Figure22
1 OM higherS
s
1 × 10−5
200.67
+67.36
Figure25
1 OM lowerS
s
1 × 10−7
74.96
−37.48
Figure26Figure 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 Figures20–33.Figure 21
Sensitivity to doubling the flux in the west.Figure 22
Sensitivity to decreasing the flux in the west by one-half.Figure 23
Sensitivity to increasingK by half an order of magnitude in aquifers.Figure 24
Sensitivity to decreasingK by half an order of magnitude in aquifers.Figure 25
Sensitivity to increasingS
s by an order of magnitude.Figure 26
Sensitivity to decreasingS
s by an order of magnitude.Figure 27
Sensitivity to increasingK by 2 orders of magnitude in layers (
1
), (
2
), and (
4
).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 29
Results of running the second model, the 18 Ma model, representing the 9-million-year period from 27 to 18 Ma, compared to present-day heads. Model curves(
1
), (
2
), and (
3
) (Wolfcampian, Desmoinesian, and Mississippian, resp.) are approaching the present-day data curves. Red triangles at centroids of the grid columns mark the initial zero heads at 70 Ma.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
Results from the 0 Ma model, representing the 2-million-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 32
Progression of Wolfcampian head solutions from 27 to 0 Ma.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 cross-section but has not reached steady state in the central and western parts of the model.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 Table3), 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 Table5 for values of K), we computed perturbed h profiles corresponding to perturbed values of flux on the western boundary (Figures 21 and 22). Doubling the flux produces perturbed h values greater than the base case, and halving the flux produces perturbed h values less than the base case. At x
=
200 km, the changes in h range from +5 m in the Mississippian unit to +32 m in the Wolfcampian unit. The relatively low changes in h 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 x-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 parametersK and S
s (Figures 23–26). The key parameter in (1) is the ratio K
/
S
s, suggesting that perturbations in K and S
s should produce opposite effects, and this is indeed the case. An increase in K drops h below the base case (Figure 23), and so does a decrease in S
s (Figure 26). Conversely, a decrease in K and an increase in S
s both produce increases in h (Figures 24 and 25).A third sensitivity test examines the effect of a much higherK in the uppermost layers of the model (Figure 27). K of Pierre Shale is maintained at K = 10−6 and the overlying and underlying units are increased to 10−4. Figure 27 shows that if K 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, h 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 Table7.


## 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 ofS
s and K (K
h and K
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, S
s, and K 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 thex direction covering a horizontal distance of 800 km and 10 layers in the z direction covering a vertical distance of 4 km, and aquifer properties (K and S
s) 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 horizontalx (length) and vertical z (depth) directions and that no flow takes place in the horizontal y 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 (Figure8). 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 atx
=
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 cross-section, showing the early development of underpressure. The eastern boundary condition (h = surface elevation) prevents the development of any underpressure at x
=
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.Figure32 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.Figure33 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.



---

*Source: 3765743-2018-02-26.xml*