Western Interior Plains Confining Unit

The Western Interior Plains confining unit is present in the northwestern and southeastern parts of the study area (fig. 1A) and tends to thin from the study area boundary towards the interior of the study area (Russell and Stivers, 2020). The Western Interior Plains confining unit is mostly composed of Mississippian- to Pennsylvanian-age rocks (fig. 9). The available thickness data for Western Interior Plains confining unit in the study area were compiled and statistically summarized. The median thickness of the Western Interior Plains confining unit is 542 ft within the study area and 551 ft for Oklahoma (table 4). The Western Interior Plains confining unit mostly consists of shales but also of interbedded layers of sandstones and limestones in smaller quantities (Jorgensen and others, 1996; Westerman and others, 2016b). Although the Western Interior Plains confining unit is generally not considered water bearing, fracturing and the weathering of near-surface rocks (less than 300 ft below land surface) of the Western Interior Plains confining unit has increased the permeability of this confining unit in some areas, facilitating shallow groundwater flow in isolated locations (Imes and Emmett, 1994). In Oklahoma, the rocks of the Western Interior Plains confining unit tend to have less shale content and higher sandstone and limestone content compared to the rocks of the Western Interior Plains confining unit in Arkansas, Kansas, and Missouri. Therefore, groundwater yields from the Western Interior Plains confining unit within Oklahoma are often greater than groundwater yields from the Western Interior Plains confining unit outside of Oklahoma (Imes and Emmett, 1994; Hays and others, 2016).

Boone Aquifer

The Boone aquifer (also known as the Springfield Plateau aquifer) is mostly unconfined within Oklahoma. From top to bottom, the rocks that contain the Boone aquifer predominantly consist of the Moorefield Formation, Keokuk Limestone, and Boone Formation (fig. 9). In Oklahoma, the Moorefield Formation consists of calcareous siltstones overlying black shale and oolitic limestone (Huffman, 1958) and is considered as one of the rock units that contain the Boone aquifer. In Arkansas, however, the shale content is greater and thus the Moorefield Formation is considered part of the Western Interior Plains confining unit (Hays and others, 2016). The Keokuk Limestone underlies the Moorefield Formation and is present in Oklahoma, Kansas, and Missouri (fig. 9). The Keokuk Limestone is usually blue-gray and is composed of tripolitic-chert-bearing limestone (Huffman, 1958; Bingham, 1969). The Keokuk Limestone is another important water-bearing unit for the Boone aquifer in northeastern Oklahoma. The Boone Formation (consisting of the Reeds Spring and St. Joe Limestone Members) underlies the Keokuk Limestone and contains various types of crinoidal, oolitic, and chert-rich limestones with occasional beds of chert-free limestone (McKnight and Fischer, 1970). The Boone Formation is the primary water-bearing unit of the Boone aquifer (fig. 10; Imes and Emmett, 1994). The St. Joe Limestone Member typically yields more groundwater than the Reed Springs Member or Keokuk Limestone because it contains much less chert than the Reed Springs Member or Keokuk Limestone (Hays and others, 2016). The Boone aquifer typically yields less than 50 gallons per minute of groundwater (OWRB, 2012a), but some wells have yielded as much as 1,000 gallons per minute of groundwater (Imes and Emmett, 1994).

The thickness of the Boone aquifer is relatively consistent throughout the study area except in areas where the unit has eroded away (Russell and Stivers, 2020). Median thickness of the Boone aquifer is 248 ft in the study area and 238 ft for Oklahoma (table 4). The tops of the rocks that contain the Boone aquifer are exposed at land surface except when the Western Interior Plains confining unit confines the aquifer from about the border between Cherokee County and Sequoyah County towards the Arkansas River and from a few miles northwest of Grand Lake O’ the Cherokees towards the western boundary of the aquifer and the northern extent of the study area (Russell and Stivers, 2020; fig. 1).

Secondary porosity, such as karst development and fracturing, controls the ability of the Boone aquifer to transmit water. The Boone aquifer has developed karst features such as caverns, enlarged fractures, and sinkholes owing to dissolution of soluble limestone along preferential-flow planes such as bedding planes and fracture traces. In some locations, as much as 70 percent of the thickness of the Boone aquifer is associated with chert (Brahana and others, 2009). The presence of chert in the Boone aquifer reduces its susceptibility to karst development and thus reduces yields relative to other aquifers within the Ozark Plateaus aquifer system. However, chert’s susceptibility to fracturing has enhanced permeability in some locations (Imes and Emmett, 1994).

Ozark Confining Unit

The Ozark confining unit limits groundwater flow between the overlying rocks of the Boone aquifer and underlying rocks of the Roubidoux aquifer. The Ozark confining unit mostly consists of Mississippian- and Devonian-age limestones and shales (fig. 9). In Oklahoma, the Ozark confining unit is composed, from top to bottom, of the Northview Shale, the Compton Limestone equivalent, the Woodford Chert, and the Chattanooga Shale. In the study area, the Ozark confining unit is mostly composed of Chattanooga Shale, and in some locations, this shale represents the entirety of the Ozark confining unit (Imes and Emmett, 1994). The Chattanooga Shale is a carbonaceous, pyritic, and fissile black shale that contains minor amounts of white sandstone at the base of the unit (McKnight and Fischer, 1970). Although the lithologic composition of the Ozark confining unit can be highly variable within the study area, leakance through the Ozark confining unit tends to be relatively low even in locations where shale content is low (Imes and Emmett, 1994).

The Ozark confining unit is present within most of the study area except for some small areas around streams where the rocks that compose the Ozark confining unit have eroded away, exposing the rocks that contain the Roubidoux aquifer (fig. 1B). The Ozark confining unit underlies the Boone aquifer and restricts groundwater flow between the Boone and Roubidoux aquifers. The Ozark confining unit is relatively thin compared to the rest of the hydrogeologic units in the study area, and median thickness is 43 ft within the study area and 47 ft within Oklahoma (table 4).

Roubidoux Aquifer

The Roubidoux aquifer is present throughout the entire study area and is mostly confined by the Ozark confining unit except near stream channels where rocks that compose the Ozark confining unit have eroded away. The Roubidoux aquifer is contained in Devonian- to Cambrian-age rocks (primarily limestones, dolomite, and sandstones; fig. 9). The Arbuckle Group is the thickest and highest-yielding unit of the Roubidoux aquifer (Christenson and others, 1990; Hays and others, 2016). The upper formations of the Arbuckle Group (from top to bottom, the Powell Dolomite, Cotter Dolomite, and Jefferson City Dolomite) are composed of dolomite with chert, sandstone, and occasionally shale lenses (Adamski and others, 1995). The Roubidoux Formation underlies the Jefferson City Dolomite (fig. 9) and is composed of sandstones, dolomite, and cherty dolomite (Thompson, 1991). The Gasconade Dolomite underlies the Roubidoux Formation; this unit occasionally contains sandy lenses, which typically yield large amounts of water (Imes and Emmett, 1994). At the base of the Gasconade Dolomite is the Gunter Sandstone Member, which is a fine- to coarse-grained sandstone with some dolomitic content (MacDonald and others, 1977). Underlying the Gunter Sandstone Member are the Eminence and Potosi Dolomites. The Eminence and Potosi Dolomites are similar in composition, consisting of fine- to coarse-grained dolomite with dense chert, drusy quartz, and occasionally glauconitic shale (Caplan, 1957, 1960). The Roubidoux aquifer is underlain by the St. Francois confining unit, which is composed of the Doe Run Dolomite, Derby Dolomite, and Davis Formation, all of Cambrian age. This confining unit restricts flow between the Roubidoux aquifer and the St. Francois aquifer. In the study area, the St. Francois confining unit and St. Francois aquifer are relatively thin (Russell and Stivers, 2020), and groundwater flow between the Roubidoux and St. Francois aquifers is likely negligible (Imes and Emmett, 1994).

The rocks that contain the Roubidoux aquifer dip toward the west from the Ozark uplift in southeastern Missouri towards Arkansas and Oklahoma; the Roubidoux aquifer is deepest in the southern part of the study area (Russell and Stivers, 2020). The aquifer is relatively thick, with a median thickness of 1,300 ft in the study area and 1,281 ft in Oklahoma (table 4). Secondary porosity (fracturing and karst development) is the main control on hydraulic properties (hydraulic conductivity, specific yield, and specific storage). Groundwater yields from the Roubidoux aquifer can be less than 10 gallons per minute and more than 1,000 gallons per minute (OWRB, 2012a).

Faulting

Most of the structural deformation of rocks associated with the Western Interior Plains confining unit, Boone aquifer, Ozark confining unit, and Roubidoux aquifer occurred during the Appalachian-Ouachita orogeny (Mississippian to Permian), when the rocks that contain the Ozark Plateaus aquifer system were compressed, resulting in normal and strike-slip faulting (Orndorff and others, 2001; Hudson and Murray, 2003; Weary and Schindler, 2004; Hudson and others, 2011; Hays and others, 2016). Russell and Stivers (2020) highlighted several regional, northeast-trending faults within the study area. However, offsets along the faults tend to be small relative to the overall thicknesses of the hydrogeologic units (McKnight and Fischer, 1970; Hudson and others, 2011) and do not strongly influence regional groundwater flows within the Ozark Plateaus aquifer system (Imes and Emmett, 1994).

Boone Aquifer

Most of the hydraulic properties for the Boone aquifer have been estimated by using specific-capacity tests and from calibrated values obtained from previously published numerical groundwater-flow models (table 5); hydraulic-property data from previously published aquifer tests conducted by using multiple wells were not available. Horizonal hydraulic conductivity for the Boone aquifer generally ranged from 0.2 to 35.0 ft/d (table 5), but because of heterogeneity, horizontal hydraulic conductivity could be much higher locally in areas of extensive fracturing and karst development. Vertical hydraulic conductivity for the Boone aquifer is difficult to measure, and estimates obtained from groundwater-flow modeling were generally less than those for horizontal hydraulic conductivity. Estimates of vertical hydraulic conductivity ranged from 2.97×10⁻⁴ to 1.0 ft/d (table 5; Czarnecki and others, 2009). Estimates of specific yield (dimensionless) ranged from 0.003 to 0.20. Specific storage estimates ranged from 1×10⁻⁵ to 5×10⁻⁵ ft⁻¹ and were based on values from previous groundwater-flow modeling studies (table 5). Hydraulic properties for the Boone aquifer were generally less variable than those for the Roubidoux aquifer.

Roubidoux Aquifer

Hydraulic properties for the entire Roubidoux aquifer were estimated from several multi-well aquifer tests, specific-capacity tests, and numerical groundwater-flow modeling (table 5). Horizontal hydraulic conductivity for the Roubidoux aquifer generally ranged from 8.64×10⁻⁴ to 86.4 feet per day (ft/d) (table 5; Imes and Emmett, 1994). Hydraulic conductivities in the study area likely represent a smaller range of values (about 0.1 to 31.9 ft/d; table 5; Macfarlane, 2007; Czarnecki and others, 2009). Vertical hydraulic conductivities for the Roubidoux aquifer are difficult to measure, and estimates from groundwater-flow modeling are generally less than estimates for horizontal hydraulic conductivities (table 5). Vertical hydraulic conductivity estimates ranged from 3.28×10⁻⁵ to 3.22 ft/d. Specific yield estimates from groundwater-flow modeling varied from 3×10⁻⁴ to 7.4995×10⁻², and specific storage estimates varied from 1×10⁻⁸ to 5×10⁻⁴ ft⁻¹ (table 5).

Multi-well aquifer tests provide a robust methodology to estimate the hydraulic properties of an aquifer (Freeze and Cherry, 1979; Heath, 1983). For the Roubidoux aquifer, multi-well aquifer tests were conducted in Pittsburg, Kans. (Stramel, 1957; Macfarlane, 2007) and in Ottawa County, Okla. (Reed and others, 1955; Christenson and others, 1990).

Two multi-well aquifer tests were conducted in Pittsburg, Kans., for the Roubidoux aquifer (fig. 1; Macfarlane, 2007; Stramel, 1957). Estimated horizontal hydraulic conductivity from the first multi-well aquifer test ranged from 25.5 to 31.9 ft/d (table 5; Macfarlane, 2007). This estimate of hydraulic conductivity is typically higher than estimates of hydraulic conductivity from specific-capacity tests for the northern part of the study area (Macfarlane, 2007) and is likely higher than hydraulic conductivity within the study area (Imes and Emmett, 1994; Clark and others, 2018, 2019) or the Ozark Plateaus aquifer system (Pugh, 2008). Hydraulic conductivity was not calculated for the second study (Stramel, 1957), and estimated transmissivities were about two times greater than those of the first study (Macfarlane, 2007). A Theis solution (Heath, 1983) was used to estimate transmissivities for the multi-well aquifer test in the second study, and a Cooper-Jacob solution (Cooper and Jacob, 1946) was used to estimate transmissivities in the first study. A Cooper-Jacob solution is likely more appropriate for karst aquifers (Kresic, 2013).

Two studies analyzed multi-well aquifer tests for the Roubidoux aquifer in Ottawa County, Okla. The first study found that a Theis solution did not have strong agreement with data from the aquifer tests (Reed and others, 1955). Based on estimates of saturated thickness and transmissivities from this study, hydraulic conductivity ranged from about 9 to 13 ft/d (Reed and others, 1955). The second study re-analyzed these same multi-well aquifer tests assuming a leaky confining unit (Christenson and others, 1990). Estimates of transmissivities for the second study were about an order of magnitude less than those of the first study (Reed and others, 1955; Christenson and others, 1990; table 5). Estimates from the first study (Reed and others, 1955) were thought to be more representative of hydraulic conductivity for the Roubidoux aquifer near Ottawa County because these values are more similar to hydraulic conductivity estimated by using specific-capacity tests within the area (Macfarlane, 2007) and from previous groundwater modeling studies (Imes and Emmett, 1994; Czarnecki and others, 2009; Clark and others, 2018; table 5). Early-time groundwater levels were used to estimate hydraulic properties in the second study (Christenson and others, 1990). If karst features (such as conduits, vugs, or caves) were present in the aquifer, groundwater responses indicated by early-time groundwater levels would mostly represent water drained from the karst features and would be less representative of water drained from the aquifer matrix (Kresic, 2013). Using early-time groundwater levels with a Theis solution would result in lower estimates of hydraulic properties compared to using late-time groundwater levels for aquifers with karst development (p. 112 of Kresic, 2013). Thus, estimates for hydraulic properties from the first study (Reed and others, 1955) are likely more representative of the Roubidoux aquifer than estimates from the second study (Christenson and others, 1990).

Dissolved-Solids Concentration Analysis

An analysis of dissolved-solids concentrations was done to assess potential salinity migration from evaporite-rich rock units in western Oklahoma and Kansas into the aquifers in the study area, and a geochemical analysis of dissolved ions was done to better understand the hydrologic connection between the Boone and Roubidoux aquifers. A saline transitional zone for the deep aquifers in the study area (which are typically contained in rocks of Paleozoic age) exists along the western boundary of the Ozark Plateaus aquifer system (fig. 12; Macfarlane and Hathaway, 1987; Jorgensen and others, 1993; Imes and Emmett, 1994; Macfarlane, 2010; Nelson and others, 2015). The identified saline transitional zone extended further eastward in Kansas and Missouri compared to Oklahoma (Macfarlane, 2010; Nelson and others, 2015). Groundwater within the saline transitional zone tends to be less saline than deeper groundwater residing in Paleozoic rocks west of the study area but slightly more saline than the groundwater residing in Paleozoic rocks in the study area (figs. 1A; Nelson and others, 2015). Dissolved-solids concentrations from groundwater in the Boone and Roubidoux aquifers within Oklahoma were mostly less than 5,000 mg/L (fig. 12).

Dissolved-solids concentrations in groundwater-quality samples collected from wells completed in the Roubidoux aquifer were usually higher than those in samples collected from wells completed in the Boone aquifer (fig. 13). Dissolved-solids concentrations listed in the GMAP database for groundwater-quality samples from wells completed in the Boone aquifer in Oklahoma ranged from 122 to 472 mg/L, with a mean concentration of 236 mg/L (OWRB, 2017). Dissolved-solids concentrations in groundwater-quality samples from wells completed in the Roubidoux aquifer ranged from 136 to 1,281 mg/L, and the mean and median concentrations were 474 and 330 mg/L, respectively. None of the dissolved-solids concentrations listed in the GMAP database were greater than the OWRB (2017) regulatory threshold of 5,000 mg/L for freshwater. Some concentrations for dissolved solids, however, did exceed the EPA’s secondary drinking water standard of 500 mg/L for dissolved-solids concentrations (EPA, 2020a).

Within Oklahoma, hardness for groundwater-quality samples from wells completed in the Boone aquifer ranged from 100 to 326 mg/L of calcium carbonate, with a mean of 195 mg/L of calcium carbonate (OWRB, 2017). Hardness for groundwater-quality samples from wells completed in the Roubidoux aquifer ranged from 38 to 435 mg/L of calcium carbonate, with a mean of 141 mg/L of calcium carbonate.

The EPA sets primary drinking water standards for various contaminants (EPA, 2020a). The maximum contaminant level (MCL) is the permitted contaminant concentration limit for treated drinking water (EPA, 2020a). Although groundwater-quality samples collected from wells completed in the Boone or Roubidoux aquifer are not the same as treated drinking water, the EPA drinking water standards are still useful benchmarks for assessing groundwater quality. One groundwater-quality sample from the Boone aquifer exceeded the MCL for nitrate (10 mg/L; EPA, 2020a). Three groundwater-quality samples from the Roubidoux aquifer contained dissolved radium-226 and radium-228 concentrations exceeding the EPA’s MCL (5 picocuries per liter; OWRB, 2017). No other analyzed constituents exceeded the MCL standards (OWRB, 2017).

Dissolved Major-Ion Concentration Analysis

Dissolved major-ion concentration data from the GMAP program were plotted on a trilinear Piper diagram (Piper, 1944) to evaluate the hydrochemical facies, or dominant water types, for groundwater sampled from each aquifer by determining the relative abundance of major ions (fig. 13). Calcium, magnesium, sodium, potassium, sulfate, fluoride, chloride, bicarbonate, and carbonate were the major ions considered in the analysis. Groundwater-quality samples were excluded if the error for the milliequivalent ionic balance of the constituents was 10 percent or more. The hydrochemical facies were used to better understand the hydrogeologic connection between the Boone and Roubidoux aquifers and identify any groundwater-quality or groundwater-flow patterns.

The groundwater-quality samples were characterized as mostly calcium-bicarbonate dominant and sodium sulfate-chloride dominant (fig. 13). Groundwater-quality samples from wells completed in the Roubidoux aquifer generally had higher dissolved-solids concentrations compared to those from the Boone aquifer. When dissolved-solids concentrations were higher, groundwater-quality samples typically contained relatively large sodium and chloride ion concentrations. Groundwater-quality samples from wells completed in the Boone aquifer were mostly calcium-bicarbonate dominant; magnesium ions were present in lesser proportions than in the groundwater samples from the Roubidoux aquifer. Magnesium, calcium, and bicarbonate ions were likely from dissolution of the aquifer matrix by groundwater. The dissolved major-ion concentrations in groundwater-quality samples collected from wells completed in the Roubidoux aquifer were generally more dominated by magnesium and sodium ions compared to samples collected from wells completed in the Boone aquifer. Major-ion concentration differences indicate minimal groundwater flow between the Boone and Roubidoux aquifers in many areas. However, the hydrogeologic connection between the Boone and Roubidoux aquifers is difficult to conclusively determine from dissolved major-ion concentrations because both aquifers are contained in rocks composed of similar geology and geochemical characteristics.

Most of the samples for analysis of groundwater quality in the Roubidoux aquifer were collected in the northwestern part of the study area (fig. 12). Major-ion concentrations indicate that the water types of groundwater-quality samples in the northwest were generally sodium-sulfate dominant. These groundwater-quality samples were collected from wells that were close to the transitional saline zones identified by Macfarlane (2010) and Nelson and others (2015). The source of the sodium chloride ions was attributed to eastward flow of saline groundwater into the Roubidoux aquifer. A few samples collected from the Boone and Roubidoux aquifers in the southern part of the study area were sodium-chloride dominant. Salinity in these areas could potentially result from encroachment of saline groundwater from the west or the south (Macfarlane, 2010; Nelson and others, 2015). Groundwater salinity is high (more than 30,000 mg/L dissolved-solids concentration) for wells completed in Paleozoic-age rocks south and west of the Arkansas River (fig. 12; Nelson and others, 2015). Areas of groundwater containing elevated dissolved-solids concentrations (more than 5,000 mg/L) are isolated, and saline groundwater is not likely transported by groundwater flow, because regional groundwater-flow directions are generally towards the more saline groundwater in the west and the south (figs. 10, 11).

Dominant water types of samples collected from wells completed in different aquifers in the study area are depicted in figure 14 in relation to the spatial extent of the different hydrogeologic units in the study area. One groundwater-quality sample collected from a well completed in the Roubidoux aquifer was classified as calcium-magnesium sulfate-chloride water type. Eight groundwater-quality samples (seven from the wells completed in the Roubidoux aquifer and one completed in the Boone and Roubidoux aquifers) were classified as sodium sulfate-chloride water type. The remaining 44 groundwater-quality samples (from 31 wells completed in the Boone, 7 wells completed in the Roubidoux aquifer, and 6 wells completed in both the Boone and Roubidoux aquifers) were classified as calcium-magnesium bicarbonate water type. Typically, sulfate ions in groundwater in the Ozark Plateaus aquifer system in northern Oklahoma, eastern Kansas, and southwestern Missouri are derived from sulfate-bearing minerals within shales of the confining units of the Boone aquifer (Imes and Emmett, 1994; Adamski and others, 1995). Sulfate-bearing minerals are abundant in the Tri-state mining district (fig. 1; Spruill, 1987; Christenson, 1995; Petersen and others, 1998; DeHay, 2003; DeHay and others, 2004); sulfate ions in groundwater in the Tristate mining district derive mostly from sphalerite, galena, and chalcopyrites formed from hydrothermal fluids transported along fractures and fault planes (Wenz and others, 2012). High concentrations of sulfate ions in groundwater-quality samples from the Roubidoux aquifer likely indicated that the Boone and Roubidoux aquifers had a strong hydrogeologic connection where these samples were collected. Alternatively, many groundwater wells in the Picher mining district were completed in both the Boone and Roubidoux aquifers but have since been capped (IT Corporation, 1985). Abandoned wells open to both aquifers could also create a pathway for sulfate-rich groundwater to flow from the Boone aquifer into the Roubidoux aquifer.

Overall, the water type of groundwater-quality samples collected from wells completed in the Boone aquifer tended to be calcium-bicarbonate dominant, and the water type of samples collected from wells completed in the Roubidoux aquifer tended to be calcium-bicarbonate or sodium sulfate-chloride dominant (fig. 14). In areas where the sampling results indicate that the Roubidoux aquifer contained groundwater that was sodium sulfate-chloride dominant, groundwater-quality samples collected from wells completed in the overlying Boone aquifer were not sodium-chloride dominant, indicating that mixing between the aquifers was likely minimal in these areas. Mixing is also likely minimal where the apparent groundwater age is on the order of only a few years. Apparent groundwater ages vary considerably in the study area, ranging from 3 to thousands of years (Jurgens and others, 2022). Groundwater-quality samples collected from wells completed in the Roubidoux aquifer contained higher ionic proportions of magnesium, likely because of the thicker dolomite-bearing formations in the rock units that contain the Roubidoux aquifer compared to those that contain the Boone aquifer (Imes and Emmett, 1994). Magnesium represented more than 20 percent of the ionic composition in some of the groundwater-quality samples collected from wells completed in the Boone aquifer (fig. 13). Dolomite is present in the Boone aquifer, and higher magnesium content could indicate a higher presence of dolomite in these areas. The different hydrochemical patterns for the waters of the Boone and Roubidoux aquifers indicate mixing of groundwater between the two aquifers is likely minimal at a regional scale but may be appreciable at a local scale.

Rorabaugh Method

The Rorabaugh method (Rorabaugh, 1964; Rutledge, 1998) is a hydrograph method for estimating basin-scale recharge. Rorabaugh-method-estimated recharge is proportional to the change in discharge between the pre-storm and post-storm recession curves at the critical time (the time after a peak where a linear recession is observed on a semi-log, daily discharge hydrograph) and is calculated by using the following equation:

For the Rorabaugh method, USGS streamgages were selected if the watershed drained mostly areas overlying the Boone and Roubidoux aquifers, and, for streamgages on the same stream, the farthest downstream streamgage that did not violate the assumptions or best practices of the Rorabaugh method (Rutledge, 2000). The Rorabaugh method is most appropriate for longer (annual) timeframes, basins with steeper relief, streams with watershed areas of 1–500 square miles, and streams that are gaining flow from the groundwater system (Rutledge, 2000). Some streamgages on smaller order streams were within the same watershed as larger order streams; however, long-term recharge was estimated by using a weighted mean based on the watershed area upstream from each streamgage. Compared to larger order streams, smaller order streams contribute a smaller component to the mean recharge because they drain smaller watershed areas, and any bias for the watershed-area-weighted mean annual recharge estimation is likely minimal where a greater number of streamgages were used within the selected watersheds.

RORA is a computer program within the USGS Groundwater Toolbox (Barlow and others, 2015) named after the developer of the Rorabaugh method (Rorabaugh, 1964) that allows the user to estimate recharge from streamflow by using a partially automated approach of the Rorabaugh method. For this program, streamflows and watershed areas from the USGS Groundwater Toolbox are downloaded from NWIS (USGS, 2021a, b). The RORA program was applied to selected streamgages (table 7) for continuous periods of record within the study period (1980–2017). To estimate long-term recharge, mean annual RORA-program estimated recharge values for each streamgage were averaged by computing a weighted mean; the weights consisted of the watershed area for a given streamgage divided by the sum of the watershed areas for all of the streamgages used for the RORA program. Watershed-area-weighted mean RORA-program estimated recharge (hereinafter referred to as “RORA-program estimated recharge”) was estimated to be 10.62 in/yr (tables 6, 7).

RORA-program estimated recharge was similar to recharge estimated by Hays and others (2016) by using the SWB method for the entire Ozark Plateaus aquifer system (mean annual recharge of 11 in/yr). Hays and others (2016) estimated recharge ranging from 0.4 to 15 in/yr (about 30 percent of precipitation) for the entire Ozark Plateaus aquifer system. Hydraulic properties within the study area tend to be slightly lower than other areas of the Ozark Plateaus aquifer system (Imes and Emmett, 1994), and thus recharge is likely less for the study area compared to the entire Ozark Plateaus aquifer system. For the study period (1980–2017), the RORA-program estimated recharge of 10.62 in/yr constituted about 23 percent of mean annual precipitation (45.7 in.; figs. 3A, 8E). RORA-program estimated recharge of 10.62 in/yr is similar to the 10.5 in/yr of recharge estimated for the Boone aquifer by the OWRB (OWRB, 2012a) and is in the range of 5 to 15 in/yr from Dugan and Peckenpaugh (1985) within the study area (fig. 1A).

RORA-program estimated recharge was used for conceptual model recharge estimates by applying the estimate to the land-surface areas where the rocks that contain the Boone and Roubidoux aquifers are exposed at land surface. Additionally, RORA-program estimated recharge was compared to recharge estimates from other methods used in this study.

Water-Table Fluctuation Method

The WTF method is a recharge estimation technique that assumes recharge is the main factor causing short-term, daily to weekly rises in water-table altitude following precipitation events (Healy and Cook, 2002). This method is best suited for groundwater wells instrumented with continuous recorders (hereinafter referred to as “groundwater gages”) with hydrographs that display rapid rises in response to precipitation (Healy and Cook, 2002). The WTF method estimates recharge at a groundwater well by using the following equation:

Several groundwater gages were installed in the Boone and Roubidoux aquifers (fig. 1B). Two groundwater gages—well 5g (USGS well 355741094461001; table 1) in Adair County and well 9g (USGS well 361317094594101; table 1) in Delaware county (fig. 1B)—were completed in the Boone aquifer and were selected for the WTF method based on their suitability with respect to the assumptions inherent to the WTF method (fig. 15A, B). Other groundwater gages were excluded from the WTF method for various reasons. For example, if the groundwater gage was near a stream, it could be difficult to determine if rises were caused by precipitation or by groundwater inflows. Groundwater gages were also excluded if they were near karst development or chert, which could potentially affect vertical recharge (figs. 1B; for example, groundwater gages 6g, 12g, 15g; 15C) (p. 125 in Kresic, 2013; Murdoch and others, 2016). Other reasons to exclude gages were if the hydrographs displayed a subdued, delayed, or nonexistent response to precipitation (for example, groundwater gages 3g, 10g, and 14g; fig. 15D), or the well at the gage was completed under confined conditions (for example, groundwater gages 2g, 8g, 11g, and 13g; fig. 15E). In these instances, the WTF method would not be suitable for the groundwater gages because recharge was not rapidly reaching the water table or rises were not likely solely caused by recharge (if, after each precipitation event, groundwater levels rose to approximately the same level regardless of the quantity of precipitation). A dimensionless specific yield value of 0.01 was used for the WTF method based on specific yield estimates from previous studies of the Boone aquifer (table 5).

Mean WTF-estimated recharge for 2018–19 was 6.9 in/yr (10.6 percent of mean annual precipitation; table 8). This value was about two-thirds as much as the mean RORA-program estimated recharge (10.62 in/yr; table 7). Many factors, such as groundwater levels not responding to precipitation or an underestimation of specific yield, could contribute to low estimates of recharge. Also, the periods analyzed were not the same; data collected during 2018–19 were used for the WTF method, and data collected during 1980 to 2017 were used for the RORA method (periods of record for some streamgages did not span this entire period). Additionally, precipitation for the WTF analyzed period was typically higher than normal (fig. 3A).

WTF-estimated recharge was not used to develop estimates of long-term recharge for the conceptual model for several reasons. The WTF-estimated recharge was likely less reliable than the estimate obtained by using the RORA method because specific yield for the Boone aquifer is likely highly variable and the period of record of groundwater gages used in the WTF method was short (2 years).

Soil-Water-Balance Method

For the 1980–2017 study period, spatially distributed recharge within the study area was estimated by using the SWB method (Westenbroek and others, 2010). The computer code in the SWB method requires climatological (precipitation and temperature) and landscape (soil-water storage capacity, hydrologic soil group, land-surface gradient, and land-cover type) characteristics as input files in order to estimate recharge. Various sources were used to compile data for the input files to the SWB computer code; these input files and the output files are included in the associated data release (Trevisan and others, 2024). The SWB method uses a gridded data structure to determine the amount of infiltration resulting from precipitation that exceeds the storage capacity of the plant root zone (the depth below ground surface to which plant roots extend) and evapotranspiration demand.

The SWB method is used to simulate recharge by applying a modified version of the Thornthwaite-Mather method (Thornthwaite and Mather, 1957). The grid used for the SWB method consisted of 354 rows and 261 columns discretized into cells; each cell covered an area of 2,000 by 2,000 ft. The Hargreaves and Samani (1985) method for a reference latitude range of 34.5–35.5 degrees was used to estimate evapotranspiration for the SWB recharge simulation. Land-cover types were used to assign curve numbers (related to antecedent runoff potential) and plant root-zone depths for each cell. The default SWB plant root-zone depths ranged from 0 to 4.5 ft (Westenbroek and others, 2010). Because the root-zone depths are not well known and are highly sensitive to SWB recharge estimation, the root-zone depths were adjusted to calibrate to the conceptual model recharge (RORA-program estimated recharge). Root-zone depths from Westenbroek and others (2010) were for plants in Wisconsin; however, root-zone depths for Oklahoma are largely unknown, and root zones in the shales and limestones of the confining units and aquifers in the study area likely preclude deep root penetration (Schenk, 2008; Estrada-Medina and others, 2013). Root-zone depths were therefore reduced to better match conceptual model recharge (RORA-program estimated recharge). To match RORA-program estimated recharge, only the part of the study area where the Western Interior Plains confining unit was not present was used to calculate mean recharge for each array in the SWB simulation (1980–2017). The final root-zone depths were scaled to 27 percent of the original root-zone depths from Westenbroek and others (2010).

Daily precipitation, minimum temperature, and maximum temperature grids from 1980 to 2017 were obtained from the Daymet database (version 4; Thornton and others, 2020) and were used as inputs to the SWB code. Annual precipitation in the Daymet database (fig. 16C) was slightly higher than that from the PRISM Climate Group (2022) (fig. 3A); however, daily gridded climate data were not available from PRISM Climate Group (2022) for the entire 1980–2017 study period. Tabular summaries of the Daymet (Thornton and others, 2020) and PRISM Climate Group (2022) annual and monthly precipitation and annual and monthly SWB recharge can be found in Trevisan and others (2024). Soil-water storage capacity and hydrologic soil group data were obtained from the Gridded Soil Survey Geographic Database (U.S. Department of Agriculture, 2021) and sampled to the SWB model grid resolution. The land-surface gradient was obtained by using the D8 method (O’Callaghan and Mark, 1984) and was calculated from a one-third arc-second DEM (USGS, 2019) after resampling to the SWB model grid resolution using a mean land-surface altitude for each cell and filling sinks. Land-surface gradient is used by the SWB code to route runoff between cells. Land-cover types were obtained from the National Land Cover Database (Multi-Resolution Land Characteristics Consortium, 2021) and resampled to the SWB model grid resolution by using the most common land-cover type within each cell.

Mean annual SWB recharge was 10.52 in. (fig. 16C) for parts of the study area where the Western Interior Plains confining unit was not present and was about 22 percent of Daymet mean annual precipitation (Thornton and others, 2020) over the 1980–2017 study period (48.65 in.; fig. 16C). Mean monthly and mean annual SWB recharge values were slightly lower in Oklahoma (0.78 and 9.31 in., respectively; Trevisan and others, 2024). Mean annual SWB recharge was typically less than 9 in. and ranged from 0 to 32.76 in. (fig. 17).

The locations with the highest mean annual SWB recharge were generally in the central and northeastern parts of the study area where rocks that contain the Boone and Roubidoux aquifers are exposed at land surface (fig. 17). SWB recharge tends to be lowest in the northwestern and southern parts of the study area. Where the Western Interior Plains confining unit is the uppermost hydrogeologic unit, mean annual SWB recharge was typically less than 9 in., and frequently less than 6 in. (fig. 17). All open-water features (lakes and ponds) identified in the study area did not receive any recharge because during the SWB simulation, recharge is set to zero over open water by the SWB code. Precipitation that falls on open water is intercepted and thus does not directly contribute to recharge (Westenbroek and others, 2010).

Conceptual Model Seepage Estimation

Conceptual model seepage was determined by using several streamgages at the most downstream location within each watershed. Streamgages along dammed reaches (reaches downstream from dams where the stream is not free flowing) were not used in the analysis. Base flow estimated at streamgages was assumed to be equal to net seepage occurring upstream within the watershed of each streamgage.

Watersheds were delineated for each streamgage by using the USGS StreamStats tool (Ries and others, 2017). Small parts of some of the watersheds included areas where the Western Interior Plains confining unit was present, but seepage was assumed to be entirely from the Boone or Roubidoux aquifers for this analysis. Base flows were scaled by the percentage of watershed area within the study area to determine conceptual model seepage. Areas without streamgages were assumed to contribute base flows equal to the mean annual long-term (1980–2017) base flow per area of watershed for all the analyzed streamgages (table 9). Total estimated seepage was scaled to the areal extent of the surficial exposure of the rocks that contain each aquifer to estimate conceptual model seepage for the Boone and Roubidoux aquifers. These outcrop areas were then scaled to the part of the study area in Oklahoma in order to estimate the total seepage for Oklahoma.

Seepage estimates for the Boone and Roubidoux aquifers represented 58.32 percent of recharge, which was near the middle of previously estimated percentages of stream seepage to recharge (20–90 percent; table 6; Imes and Emmett, 1994; Hays and others, 2016).

Synoptic Streamflow-Measurement Analysis

Synoptic streamflow measurements for streams that drain the Boone and Roubidoux aquifers were compiled for seepage analysis (fig. 18; USGS, 2021a, b). The data were from January 31 to February 1, 2018, a period when runoff was likely small and streamflows were likely supported by base flows. Additionally, continuously gaged streamflow data obtained from USGS streamgages within the study area were used to estimate seepage (water gained or lost by the stream from groundwater inflows or outflows, respectively) for several stream segments. Seepage was normalized to the length of the stream segment for which discharge was measured to compare seepage between stream segments of different lengths (Hays and others, 2016; Ellis, 2018; Ellis and others, 2020). Synoptic streamflow measurements and continuously gaged streamflow records are available in NWIS (USGS, 2021a, b).

When a synoptic streamflow measurement was made directly upstream or downstream from a streamgage, the nearest continuously gaged streamflow record from that streamgage was used in the stream seepage computation. If no streamflow measurement was made directly upstream or downstream of a streamgage, the nearest continuously gaged streamflow record corresponding to the mean synoptic streamflow measurement time was used in the stream seepage computation. By using this approach, streamflow was assumed constant between the time of each measurement such that travel times of water down a stream were assumed to be insignificant for determining seepage. For determining upstream streamflow, if multiple synoptic measurements or streamgages on multiple tributaries were upstream of a synoptic measurement, then the sum of streamflow from those streamgages or synoptic measurements was subtracted from the downstream streamflow value. Tributaries with no measurements or streamgages were assumed to contribute no seepage. Stream lengths were calculated by using the National Hydrography Dataset (NHD; EPA, 2020b). Seepage values were estimated by taking the change in streamflow divided by the stream lengths.

Seepage was normalized and classified into three categories: (1) gaining, (2) neutral, or (3) losing (fig. 18). Gaining stream segments were those for which normalized seepage was greater than 0.5 cubic foot per second per mile ([ft³/s]/mi), neutral stream segments were those for which normalized seepage was at or between 0.5 and −0.5 (ft³/s)/mi, and losing stream segments were those for which stream normalized seepage was less than −0.5 (ft³/s)/mi. Streams in the study area tend to be gaining; however, three segments in the study area were losing (fig. 18). Five stream segments were neutral. Losing stream segments tended to be shorter compared to the gaining stream segments. Seepage rates from the seepage analysis were not used for the conceptual model seepage estimates because they represent a short period of time and streamflows were measured following a year of above normal precipitation (2017). Therefore, the seepage values estimated by using the streamflows were not likely representative of long-term seepage for the study area.

Recharge

Recharge to the groundwater-flow model was simulated by using the Recharge (RCH) package (Harbaugh, 2005). The RCH package requires a user-specified recharge rate for each simulated cell and applies that rate over the specified time interval of each stress period. Monthly estimated aquifer recharge rates from the SWB code were used as inputs to the groundwater-flow model. Recharge was specified to occur in the top active cell, which is in the first layer of the model for every cell in layer 1. Recharge was initially set to a multiplier of 1.0 and varied within the range of 0.5 to 1.2 during calibration for each stress period. The multiplier is the value specified in the RCH package that is multiplied with each value in the recharge array to produce a new recharge array used by MODFLOW-NWT.

Lateral Groundwater Flow

Lateral groundwater flow across the study area boundary was simulated by using the General Head Boundary (GHB) package (Harbaugh, 2005). The GHB package is a groundwater-altitude-dependent boundary that requires the input of GHB altitude and conductance values for each GHB package cell. Conductance is computed by multiplying the hydraulic conductivity by the saturated cross-sectional area of a GHB package cell orthogonal to groundwater flow, and then dividing this sum by the distance between cell centroids. Groundwater flow between GHB cells and active cells in the model is calculated as the difference in groundwater-altitude values between the GHB cell and the individual cell in the model multiplied by the conductance of the GHB cell (Harbaugh, 2005).

GHBs were placed along the northern edge of the study area and along the eastern edges of the study area in each layer (figs. 20, 21A–C). Groundwater altitudes for each layer of the groundwater-flow model were estimated by using groundwater-altitude observations near the study area boundaries (USGS, 2021a, b) and by using potentiometric-surface maps (Imes and Emmett, 1994; Nottmeier, 2015). Groundwater altitudes for the GHB boundary were then calibrated by using bounds of 50 ft above and below the initial GHB altitude; groundwater altitudes were then adjusted to ensure the bounds were not more than model land-surface altitude or less than the altitude of the base of the layer. Initial GHB conductance values were estimated by using initial hydraulic conductivity estimates and the estimated saturated thickness for each GHB cell. GHB groundwater altitude and conductance were adjusted during calibration.

Lakes

The Time-Variant Specified-Head (CHD) package (Harbaugh, 2005) was used to simulate lake seepage to and from the major lakes and the hydrogeologic units in the study area (OWRB, 2021): Grand Lake O’ the Cherokees, Lake Eucha, Spavinaw Lake, Lake Hudson, Fort Gibson Lake, Webbers Falls Reservoir, Tenkiller Ferry Lake, and Robert S. Kerr Reservoir (fig. 1). Lake seepage for smaller lakes in the model area was assumed negligible.

The CHD package requires the input of beginning and ending lake-surface altitudes (hereinafter referred to as “lake stage”) for each stress period. For each model stress period, the beginning lake stage was set to the monthly lake stage from the stress period corresponding to the previous month, and the ending lake stage was set to the monthly lake stage from the month corresponding to a given stress period. For the steady-state stress period, beginning and ending lake stages were set to conservation pool altitude (USACE, 2021).

For most lakes, monthly lake stage was obtained from the USACE Water Control Data System (USACE, 2021) unless otherwise specified. Monthly lake stage was only available from November 1994 to the present. Lake stage before November 1994 was set to conservation pool altitude (USACE, 2021).

For Grand Lake O’ the Cherokees, lake-stage values inferred from the rating at USGS station 07190000 were used to estimate daily lake-storage values (USGS, 2021a, b). The rating depicts the relation between lake stage and lake storage and was computed in the 1940s; the full rating was published in Hunter and others (2020). The rating consists of 1-ft increments of lake stage and corresponding lake storage. The method used for estimating lake stage included calculating a linear equation for each 1-ft increment in lake stage consisting of a slope (change in lake stage, in feet, divided by the change in lake storage, in acre-feet) and an intercept (offset, in feet). Lake stage was estimated by finding the smallest storage value in the rating that was greater than the published lake storage. The slope and intercept at this section of the rating were then used to estimate a lake stage. There was some overlap between lake stage published in NWIS and the implementation of the rating used to publish lake storage. The overlapping data were used to determine a root mean square error of approximately 0.10 ft of lake stage. The statistics for the lake stage are included in the data release associated with this report (Trevisan and others, 2024). Daily lake-stage estimates were converted to monthly values by taking the mean of the daily lake-stage estimates for each month.

Available lake-stage values for Lake Eucha and Spavinaw Lake were also compiled for use in the CHD package. Historical lake-stage values prior to October 2007 were not available from the USACE for these lakes. Daily lake stages for these lakes were obtained from NWIS (USGS, 2021a, b) from October 2007 through December 2017. Daily lake stage was averaged for each month to obtain a monthly lake-stage value. Prior to October 2007, lake-stage values were set to conservation pool altitudes for Lake Eucha and Spavinaw Lake (USACE, 2021). When summarized in the simulated groundwater budget, the net groundwater flow between the lakes (CHD) and the aquifer was added to the seepage water-budget category.

Streams

The Streamflow Routing (SFR) package (version 2; Niswonger and Prudic, 2005) was used to simulate stream seepage between the stream and aquifer in accordance with Darcy’s law (Heath, 1983) and downstream routing of base flow in accordance with Manning’s equation (Chow, 1959). Stream components were separated into stream segments (portions of the stream that only contain inflows at one end and outflows at the other) and reaches (individual cells of the model along each segment).

Stream data representing the stream network in the study area were obtained from the NHD (EPA, 2020b). Streams with stream orders greater than one were used for the simulation to reduce the stream density. Minor adjustments were made to the representation of stream channels when a stream in an active cell crossed inactive cells to ensure that the stream network was continuous through the active model domain. Additionally, short stream segments (typically covering one to two cells of the model) and divergences (loops) in the NHD stream network were removed.

Streambed slopes and arbolate sums (the sum of the length of all upstream stream segments) from the NHD dataset were used for stream slopes and to estimate stream widths by using the arbolate-sum regression equation from appendix 2 on page 266 in Feinstein and others (2010). Stream widths for the Arkansas, Neosho, and Spring Rivers were reduced by 30 percent to more closely match the stream widths measured from aerial photography (Microsoft Corporation, 2021). Streambed widths were further adjusted during calibration within bounds of plus or minus 20 percent of initial estimates. Streambed altitudes were set to mean model-cell-resolution altitude based on the one-third arc-second DEM (USGS, 2019). The streambed altitudes derived from the one-third arc-second DEM altitudes (instead of the NHD streambed altitudes) were used to be consistent with the altitudes used for the rest of the groundwater-flow model.

The SFR package requires streambed altitudes to decrease in the downstream direction; streambed altitudes were adjusted where necessary by starting at stream segments with no inflows and moving reach-by-reach (cell-by-cell) downstream along the stream network to the outflow. When streambed altitudes increased downstream, they were removed and then linearly interpolated between the streambed altitudes of adjacent reaches. Additionally, some SFR streambed altitudes were manually adjusted by segment to better match calibration targets while maintaining an altitude of at least the minimum sampled altitude from the one-third arc-second DEM for each SFR cell.

Inflows to the SFR network were estimated by using estimated base flow at several USGS streamgages north of the study area (fig. 6C). Not all streamgage records spanned the entire study period (1980–2017). Streamgage records from streamgages north of the study area that did not span the entire study period were estimated by using the standardization by means method (Farmer and Vogel, 2013), where monthly flows are estimated by multiplying flows at a reference streamgage by the ratio of mean flow for the streamgage to be estimated divided by mean flow of the reference streamgage. Monthly inflows (as base flows) were estimated for USGS streamgage 07184500 Labette Creek near Oswego, Kans. (streamgage 3s) (table 1; fig. 1A) (using the sum of mean monthly base flows between USGS streamgage 07184000 Lightning Creek near McCune, Kans. [streamgage 2s] and USGS streamgage 07183500 Neosho River near Parsons, Kans. [streamgage 1s] as the reference base flow) and USGS streamgage 07186055 Cow Creek near Scammon, Kans. (streamgage 8s) (using mean monthly base flows for USGS streamgage 07186000 Spring River near Waco, Mo.[streamgage 7s] as the reference base flow). Estimated monthly base flows are shown in figure 22.

Permitted surface-water diversion locations were obtained from the OWRB (2019b). The Grand River Dam Authority (GRDA) has the jurisdiction to regulate surface-water diversions in northeastern Oklahoma, which is an area that is about one-third of the part of the study area that is in Oklahoma (fig. 7A). The permits issued by the GRDA were not used for the groundwater-flow model because they were not publicly available. The OWRB-issued permits used in the SFR package were less than 1 percent of total base flows, and the permits issued by GRDA were likely only a small percentage of base flows. Annual surface-water diversions were scaled by OWRB (2012b–d) water-use monthly multipliers for each monthly stress period based on the use type associated with the permit and the basin where each surface-water diversion was located (figs. 7, 23). When summarized in the simulated groundwater budget, the net groundwater flow between the streams (SFR) and the aquifer was added to the seepage water-budget category.

Seeps and Springs

For the groundwater-flow model, seeps and springs in the study area were the USGS-identified springs (USGS, 2019), which were represented by using the Drain (DRN) package (Harbaugh, 2005). The DRN package functions by removing groundwater when simulated groundwater altitudes exceed user-specified DRN altitudes; groundwater is removed at a rate proportional to the product of a user-specified DRN conductance and the difference between the simulated groundwater altitude and DRN groundwater altitude (eq. 6–10 in Harbaugh, 2005). Land-surface altitudes were point-sampled from the one-third arc-second DEM (USGS, 2019) at spring locations to use as altitudes for the DRN package. When multiple springs were located in the same cell, the mean of all sampled altitudes of all the springs in the cell was used for the DRN package altitudes for that cell. Initial DRN conductance values were set to 8,000 foot squared per day (ft²/d) and were further adjusted during calibration. When summarized in the simulated groundwater budget, the net groundwater flow between the springs (DRN) and the aquifer was added to the seepage water-budget category.

Groundwater Withdrawals

Groundwater withdrawals were simulated by using the Well (WEL) package (Niswonger and others, 2011). Groundwater withdrawals were distributed evenly among wells that shared the same groundwater permit. Annual groundwater withdrawals were discretized to monthly stress-period withdrawals based on monthly water use demand distributions for the six Oklahoma Comprehensive Water Plan water-management planning basins (basin numbers 46, 47, and 79–82) in the study area (OWRB, 2012a–d; figs. 7, 23). For wells outside of Oklahoma, but in the study area, the nearest surface-water basin in Oklahoma was used. Groundwater withdrawals were not simulated if they were located in cells used for the CHD or GHB packages; these packages could potentially simulate limitless amounts of water, and thus groundwater withdrawals would have minimal influence on simulated groundwater flow in the locations where GHB cells were used. Excluded groundwater withdrawals were a small component of total groundwater use and were mostly along GHBs outside of Oklahoma. For the first steady-state stress period, groundwater withdrawals were set to mean annual values for the 1980–2017 model period.

Automated Model Calibration

Automated model calibration was done by using parameter estimation software version 4 (PEST++ version 4; White and others, 2020a), which was modified from PEST++ version 3 (Welter and others, 2015). PEST++ version 4 contains a suite of tools including the PEST++ (PESTPP) iterative ensemble smoother (PESTPP-IES; White, 2018) software tool that was used to minimize the objective function. For minimizing the objective function, the PESTPP-IES software tool employs a Levenberg-Marquardt algorithm with an iterative ensemble smoother (Chen and Oliver, 2013). The objective function (Φ) is defined as:

PESTPP-IES minimizes the objective function by running multiple sets of ensembles where model parameters are iteratively adjusted within a user-defined range of values. In this manner, automated calibration can assess a broader combination of parameters faster than manually adjusting each parameter value. Additionally, PESTPP-IES uses singular value decomposition for regularizing solution subspace (Doherty, 2015; White and others, 2020b). Singular value decomposition improves the solution stability by identifying and adjusting the most sensitive parameters affecting observations and reduces the computation time required to complete automated calibration. Automated calibration was done by running ensembles in parallel on a high-performance computing (HPC) cluster using the parallel run manager built into the PEST++ suite (White and others, 2020b).

Calibration Parameters

During automated calibration, 10,552 parameters categorized into 29 parameter groups (table 10) were adjusted; the data associated with these parameters are published in the companion data release (Trevisan and others, 2024). These parameter groups include GHB conductance, GHB groundwater altitude, SFR parameters (streambed hydraulic conductivity, streambed thickness, and stream width), RCH monthly array multipliers, DRN conductance, and hydraulic parameters (hydraulic conductivity, horizontal-to-vertical hydraulic conductivity ratio, specific storage, and specific yield). Calibration parameters were adjusted to help reduce the objective function and more closely match calibration targets.

Pilot points are gridded points used to distribute parameters more easily throughout the model extent and estimate parameter heterogeneity without prior knowledge about the spatial parameter distribution during automated calibration (Doherty, 2003) (fig. 20). During automated calibration, pilot-point parameters are adjusted. Then, an exponentially weighted kriging procedure is performed to interpolate parameters to the groundwater-flow model grid resolution. Pilot-point parameters were used for the four hydraulic parameters in each model layer: (1) horizontal hydraulic conductivity, (2) the vertical anisotropy, (3) specific storage, and (4) specific yield. For simplicity in the calibration process, each hydraulic parameter group and pilot points within that parameter group pertained to a specific hydrogeologic unit (for example, pilot points in group hk1, hk2, hk3, and hk4 pertained to hydraulic conductivities for the Western Interior Plains confining unit, Boone aquifer, Ozark confining unit, and Roubidoux aquifer, respectively; table 10). Parameter bounds for most hydraulic properties (table 5) were estimated by Imes and Emmett (1994) and Clark and others (2018, 2019); however, hydraulic conductivities for the Boone and Roubidoux aquifers were estimated by using information from other publications. For the Boone aquifer, horizontal hydraulic conductivity was estimated by Reed and Czarnecki (2006) and Czarnecki and others (2009) as indicated in table 5. For the Roubidoux aquifer, the upper bound of horizontal hydraulic conductivity was estimated from a multi-well aquifer test (Macfarlane, 2007); because multi-well aquifer tests are not typically performed in low-water-yielding areas of an aquifer, the lower bound was set to 0.001 ft/d based on iterative testing of different hydraulic parameters during manual calibration. Groundwater flow to streams was restricted to areas around streams where the Ozark confining unit was the uppermost hydrogeologic unit, and simulated groundwater altitudes were high within those areas (more than 500 ft above the top of the model [land surface]). Thus, in locations where the Ozark confining unit was the uppermost hydrogeologic unit in layer 1, the Roubidoux aquifer hydraulic parameters were used as an estimate for realistic near-surface hydraulic parameters. Near-surface hydraulic parameters for shallower rocks in the Ozark confining unit were likely greater than those for rocks at depth in the Ozark confining unit because near-surface rocks can be prone to fracturing (particularly for the more fissile rocks contained in the Ozark confining unit, which readily fracture along bedding planes) and are often more permeable compared to deeply buried rocks of this unit (Imes and Emmett, 1994).

Calibration Targets

Calibration targets consisted of three observation groups: groundwater-altitude observations, base-flow observations, and conceptual model recharge values for the Boone and Roubidoux aquifers. For automated calibration, calibration targets were weighted when determining the objective function (the sum of squared weighted residuals). That is, residuals are multiplied by weights, squared (multiplied by themselves), then summed. Weighting is necessary to control how much each observation contributes to the objective function. Observations often have different units of measure, and using unweighted values would cause the automated calibration process to favor observations with greater residuals.

For automated calibration, each observation group was weighted to contribute approximately equally to the objective function (sum of squared residuals), except for the conceptual model recharge components which were given a weight to scale the objective function two orders of magnitude less than the objective function contributions of the other observation groups. This decision was made because the conceptual model budget components were estimated by making simplifying assumptions and using mean estimates, and because the data were sparse for large parts of the study area. Also, initial starting values for recharge parameters closely matched the conceptual model targets, and, by weighting two orders of magnitude less than other observation groups, recharge could be adjusted without automated calibration favoring conceptual model recharge over the other observed calibration data.

After the initial weights were calculated for each group, the individual observation weights were further adjusted on the basis of several criteria. For groundwater-altitude observations, those that were within Oklahoma were weighted about double of those outside of Oklahoma. Additionally, groundwater-altitude observations were weighted proportionally to the inverse of the number of observations at each location to provide a more robust spatial calibration because a substantial number of observations (10 or more) were available for 47 wells, and few observations (less than 10) were available for 508 wells. Weighting each measurement equally would bias the objective function towards wells with more measurements. Groundwater-altitude observations were weighted equally for the Western Interior Plains confining unit and the Boone and Roubidoux aquifers. Manual calibration revealed that Western Interior Plains confining unit groundwater-altitude observations were sensitive to net leakage to or from the Boone aquifer. Western Interior Plains confining unit groundwater-altitude observations were important for constraining net leakage from the Boone aquifer and were often in locations with few Boone aquifer groundwater-altitude observations.

Base-flow observations were initially weighted such that each streamgage contributed approximately equally to the objective function. During the calibration process, the objective function weighting for base-flow observations were subsequently doubled for streamgages on streams that predominantly drained areas overlying the Boone and Roubidoux aquifers within Oklahoma to improve calibration results for those areas.

Groundwater-Altitude Observations

For the study area, 2,274 water-level observations for 1979 through 2017 are available in NWIS (USGS, 2021a, b). These water-level observations and their corresponding well depths were converted to groundwater-altitude observations and well-bottom altitudes, respectively, each measured from the land-surface altitude by using a one-third arc-second DEM (USGS, 2019). Observations were assigned to aquifers on the basis of well-bottom altitude and base-of-aquifer altitude (Russell, 2020). When conflicts arose between the aquifer assigned by using well-bottom altitude and the geologic unit coded for the measurement (based on the geologic unit associated with the aquifer from Russell and Stivers [2020]), publications using groundwater-altitude observations at those wells were used to identify the correct aquifer in which those wells were completed (Marcher and Kenny, 1983; Goemaat and others, 1984; Christenson, 1995; DeHay, 2003; DeHay and others, 2004). Discrepancies that resulted from using well-bottom altitudes likely arose from the relatively coarse (2,000-ft by 2,000-ft) resolution of the base-of-aquifer altitude and uncertainty in the DEM and base-of-aquifer altitudes.

Simulated groundwater altitudes were extracted from the groundwater-flow model by using the Head Observations (HOB) package for MODFLOW-NWT (Hill and others, 2000). The HOB package is used to assign a simulated groundwater-altitude observation to a cell and requires a user-specified row offset and column offset to place the observation in the cell. In addition to spatial interpolation, the HOB package also temporally interpolates between stress periods by using a simulated groundwater altitude that is based on a user-specified time offset.

Groundwater-altitude observations for 1979 were used for the steady-state simulation because precipitation during 1979 approximated the long-term (1885–2020) mean annual precipitation (fig. 3A) and did not overlap with any of the transient simulated periods. Because the steady-state simulation lacks a time component, steady-state groundwater-altitude observations were averaged when multiple observations were made in the same well in 1979. This reduced the total number of groundwater-altitude observations from 2,274 to 2,217, leaving 38 total observations in 1979 that were used for the steady-state simulation (fig. 24). Groundwater-altitude observations during 1980–2017 were used for the transient simulation. The number of groundwater-altitude observations for the transient stress periods totaled 2,179 (fig. 24). No groundwater altitudes were identified for the Ozark confining unit. For the period 1885–2020, groundwater-altitude measurements were most frequently made during January through April (fig. 25). March was the month with the most measurements (398), and November had the fewest (88) (fig. 25).

Conceptual Model Budget Components

Components of the conceptual model were used to match simulated mean annual budget values (1980–2017). Each budget component was used for manual calibration, but for automated calibration, only long-term recharge for the Boone and Roubidoux aquifers was used. Recharge was the only parameter used for automated calibration because most other budget components were already calibrated to other calibration targets or were minor components of the conceptual model (table 6). Because MODFLOW-NWT reports groundwater flows for the entire simulation, subareas for each hydrogeologic unit were extracted by using the ZONEBUDGET utility (fig. 26A–D; Harbaugh, 1990).

Calibration Results

Calibration results for this report were evaluated on the goodness-of-fit between the simulated and the observed calibration targets (table 11). The primary statistic for assessing the calibration results was the reduction in the total objective function. Automated calibration reduced the total objective function by about 96 percent from calibration inputs (from 1,340,001 to 65,145; table 11). The conceptual model recharge objective function component value was approximately double the precalibration objective function component value (from 6,667 to 12,851). Base-flow and groundwater-altitude objective function components were reduced by about 95 and 99 percent, respectively (from 666,667 to 24,953 and 666,667 to 27,340, respectively).

Sensitivity Analysis

The Jacobian matrix represents a first-order relation between parameters and observations (White and others, 2020b). A groundwater-flow model sensitivity analysis was done by using the PEST++ Gauss-Levenburg-Marquardt (PESTPP-GLM) program (White and others, 2020b) to calculate a Jacobian matrix with dimensions of 1,163 parameters by 13,166 observations. The sensitivity analysis consisted of adjusting parameters by 1 percent and calculating a derivative (the change in the observation residual divided by the change in the parameter). Thus, higher derivatives within the Jacobian matrix indicate that a change in a given parameter affects an observation more than other parameters with smaller values in the Jacobian matrix (less sensitivity). The derivatives in the Jacobian matrix were with respect to the log of that parameter because adjustable parameters were log transformed using PESTPP-GLM (White and others, 2020b).

To complete the Jacobian matrix, the groundwater-flow model was simulated n times, where n is the total number of parameters, and another simulation was done without adjusting any parameters. Owing to the large number of adjusted parameters, the sensitivity analysis was done by parameter group, and the analysis was conducted in parallel on an HPC cluster and summarized based on parameter groups (fig. 27). Composite parameter sensitivities (eq. 5.3.1 on p. 73 of Doherty, 2015) were calculated to normalize observation sensitivities and parameter groups. For this report, sensitivity with respect to the groundwater-flow model refers to composite parameter sensitivity. A composite parameter sensitivity is the weighted, normalized column for each parameter from the Jacobian matrix with respect to the number of observations with nonzero weights. The composite parameter sensitivity is calculated by using the equation (all terms are considered dimensionless):

cps

= [

J^T Q J

]

^1/2

/

n

(4)

where

cps

is the composite parameter sensitivity;

J^T

is the transpose of the Jacobian matrix;

Q

is the matrix of weights assigned to each observation;

J

is the Jacobian matrix; and

n

is the total number of observations with nonzero weights.

Groundwater-altitude observations were consistently the least sensitive to parameter adjustment, the recharge multipliers were most sensitive, and base-flow observations were moderately sensitive during parameter adjustment (fig. 27A–C). Because recharge was only affected when adjusting the recharge multipliers (rch), all other parameter groups were insensitive (composite sensitivity of zero) for conceptual-model recharge observations (fig. 27A). Base-flow observations (fig. 27B) were generally more sensitive compared to groundwater-altitude observations (fig. 27C) during parameter adjustment. Generally, for parameters that corresponded to specific aquifers, observations were often more sensitive to parameter adjustment for the Boone and Roubidoux aquifers (parameter groups generally associated with layers 2 and 4, respectively) compared to adjusting the same parameters for the Western Interior Plains confining unit and Ozark confining unit (parameter groups generally associated with layers 1 and 3, respectively).

Observations were most sensitive to adjusting the recharge multipliers and the GHB altitudes for the Boone and Roubidoux aquifers (ghbhd2 and ghbhd4, respectively; fig. 27B, C), but hydraulic conductivity for the Roubidoux aquifer (hk4) was slightly more sensitive than the GHB altitudes for the Roubidoux aquifer with respect to groundwater altitudes (fig. 27C). Hydraulic conductivity pilot points for the Boone aquifer (hk2) were the next most sensitive parameters for the base-flow and groundwater-altitude observations. Base-flow and groundwater-altitude observations were generally less sensitive to parameter adjustment for the Ozark confining unit compared to other parameters. Thus, leakage from the Boone aquifer to the Roubidoux aquifer through the intervening Ozark confining unit may have high uncertainty. For the Ozark confining unit, observations were about as sensitive to the vertical anisotropy parameters as they were to hydraulic conductivity parameters.

Study Area Calibrated Parameters

The parameters for the groundwater-flow model (Trevisan and others, 2024) were calibrated manually and automatically by using various packages designed to work with MODFLOW-NWT as described in the “Simulation of Groundwater Flow” section of this report.

Calibrated SFR streambed hydraulic conductivity (sfrhk) and streambed thickness (sfrthick) ranged from 0.357 to 10 ft/d and 0.333 to 6 ft, respectively (table 10). The bounds for these values were not determined by field data collected as part of this study but were based on values from previous studies for other groundwater-flow models in Oklahoma (table 10; Smith and others, 2017, 2021; Ellis, 2018; Ellis and others, 2020). Because the Boone and Roubidoux aquifers that underlie the streams in the study area are bedrock aquifers, streambed thickness was assumed to be slightly less than the streambed thickness typical for streams overlying alluvial aquifers. GHB conductance values for the Boone aquifer (ghbhcd2 in table 10) ranged from 0.16 to 67,400 ft²/d, and conductance values for the Roubidoux aquifer (ghbcd4 in table 10) ranged from 3.22 to 130,000 ft²/d. DRN (representing springs) conductance values ranged from 984 to 10,000 ft²/d (drn in table 10). Conductance values for springs are difficult to determine, and springs for the groundwater-flow model lacked calibration data. Moderate conductance values were assumed for springs because they were considered minor outflows from the aquifers. RCH package array multipliers for SWB recharge grids ranged from 0.598 to 1.2 (rch in table 10).

For the Boone aquifer, horizontal hydraulic conductivity values ranged from 0.647 to 19.1 ft/d (hk2 in table 10), and for the Roubidoux aquifer, horizontal hydraulic conductivity values ranged from 0.054 to 31.9 ft/d (hk4 in table 10). The horizontal hydraulic conductivity values for both aquifers were within the range of horizontal hydraulic conductivity of 8.64×10⁻⁴ to 86.4 ft/d estimated in previous studies for the Ozark Plateaus aquifer system (Imes and Emmett, 1994; Hays and others, 2016; Clark and others, 2018, 2019).

Calibrated vertical hydraulic conductivity values for the Boone and Roubidoux aquifers were generally higher than those of previous studies (tables 5, 10); however, manual calibration revealed that using the vertical hydraulic conductivity values in table 5 would restrict drainage such that simulated base flows would be much greater than base-flow observations, groundwater altitudes would be hundreds of feet higher than the maximum measured groundwater altitude, and leakage between the Boone and Roubidoux aquifers would also be much less than the conceptual model leakage. Many groundwater-flow models were simulated during relatively drier periods when recharge was estimated to only be about 4 to 5 in/yr and as small as an order of magnitude less than the conceptual model recharge values (Reed and Czarnecki, 2006; Czarnecki and others, 2009; Clark and others, 2019). Differences among the parameters that were adjusted, the time periods analyzed, and the locations of the simulations between the groundwater-flow model of this report and the groundwater-flow models of previous reports (Reed and Czarnecki, 2006; Czarnecki and others, 2009; Clark and others, 2018, 2019) could also affect the differences in calibration of vertical conductivity. The spatial variability of vertical hydraulic conductivity (fig. 28C, D) was similar to that of horizontal hydraulic conductivity (fig. 28A, B).

Comparison of Simulated and Observed Calibration Targets

For this study, residuals were calculated as observed values minus simulated values. Thus, negative values indicate simulation overprediction, and positive values indicate simulation underprediction. The mean calibrated groundwater-altitude residual was −86.2 ft, indicating that, on average, simulated groundwater altitudes were greater than observed groundwater altitudes (table 11). Groundwater-altitude residuals for the Boone aquifer were not generally overpredicted or underpredicted spatially (figs. 29B, 30A), but simulated groundwater-altitude residuals for the Roubidoux aquifer were generally overpredicted spatially (figs. 29B, 30B). Structural features (such as fracturing, presence of chert, and karst development) were not represented in the model because of a lack of available field data, but these features may be contributing to overprediction in many of the simulations.

Long-term (1980–2017) simulated and observed groundwater altitudes generally followed similar patterns. The calibrated groundwater-flow model generally overpredicted groundwater altitudes (fig. 31). For the Boone aquifer, simulated groundwater altitude was often more variable than observed groundwater altitudes, however the opposite was true for the Roubidoux aquifer (fig. 31). For the Roubidoux aquifer, simulated groundwater altitudes were overpredicted for lower groundwater-altitude observations (less than 600 ft; wells 16g, 31g, and 35g, fig. 31). Well 35g is one example where groundwater altitudes were influenced by nearby groundwater withdrawals in this location (figs. 3D, 31). Observed groundwater altitudes were frequently influenced by groundwater withdrawals in the study area.

The calibrated numerical groundwater-flow model slightly overpredicted base flows (fig. 29C, D). The mean base-flow residual was −1.95 ft³/s, indicating that, on average, base flows were overpredicted (table 11). Simulated base flows in the central and northern parts of the study area generally agreed closely with observed base flows (USGS streamgages 07187600 Spring River near Baxter Springs, Kans. [streamgage 10s], 07188000 Spring River near Quapaw, Okla. [streamgage 11s], 07195500 Illinois River near Watts, Okla. [streamgage 30s], and 07196000 Flint Creek near Kansas, Okla. [streamgage 34s]) (table 1; figs. 1A, 32). Peak base flows and base flows towards the southern part of the study area were likely difficult to match because of model discretization effects such as the absence of simulated alluvial aquifers at the coarse (2,000-ft by 2,000-ft) cell size used in the model.

Calibrated Groundwater-Flow Model Water Budget

The calibrated groundwater-flow model water budget (table 13; fig. 33) includes mean annual inflows and outflows for the model period 1980–2017. For the calibrated groundwater-flow model, simulated recharge (2,520,887 and 1,420,898 acre-ft/yr for the study area and for the part of the study area in Oklahoma, respectively) was the primary inflow for the Boone aquifer, and simulated leakage was the primary inflow for the Roubidoux aquifer (912,606 and 103,529 acre-ft/yr for the study area and for the part of the study area in Oklahoma, respectively). For both aquifers, seepage (to streams, springs, and lakes) was the largest outflow. Lateral flow across the study area boundary was a larger component of outflows than seepage for the Roubidoux aquifer. Otherwise, lateral flow across the study area and Oklahoma boundaries were relatively minor components of the water budget. Groundwater use accounted for less than 1 percent of total outflows for the Boone aquifer (in the study area and in the part of the study area in Oklahoma) and 2.27 percent (study area) and 5.02 percent (Oklahoma) of total outflows for the Roubidoux aquifer. Simulated groundwater use for the Boone aquifer in Oklahoma (−2,613 acre-ft/yr; table 13) was slightly higher than for the conceptual model (−2,608 acre-ft/yr; table 6). The borders that Oklahoma shares with Kansas, Missouri, and Arkansas cross some cells in the model; therefore, some small groundwater withdrawals in surrounding States were included when calculating the groundwater withdrawals for Oklahoma with the ZONEBUDGET utility. The difference accounted for about 0.19 percent of total groundwater use. Also, simulated groundwater use was less than the conceptual-model groundwater use for the study area because wells that were present along the GHB boundary were not included in the simulation.

Table 13.    

Long-term simulated water budget for the Boone and Roubidoux aquifers for the entire study area in northeastern Oklahoma, southeastern Kansas, southwestern Missouri, and northwestern Arkansas and for the part of the study area in Oklahoma, 1980–2017.

[Positive and negative values indicate flows into and out of the aquifer, respectively. Values not in parentheses are in acre-feet per year. Values in parentheses are percentages of total inflows or outflows. If water-budget category was an inflow (positive) then the percentage is referenced to total inflows, otherwise the percentage was referenced to outflows. Percentages may not sum to 100 percent due to rounding]

Table 13.    Long-term simulated water budget for the Boone and Roubidoux aquifers for the entire study area in northeastern Oklahoma, southeastern Kansas, southwestern Missouri, and northwestern Arkansas and for the part of the study area in Oklahoma, 1980–2017.

Water-budget category	Study area						Study area in Oklahoma
	Boone aquifer		Roubidoux aquifer		Net total		Boone aquifer		Roubidoux aquifer		Net total

Table 13.    Long-term simulated water budget for the Boone and Roubidoux aquifers for the entire study area in northeastern Oklahoma, southeastern Kansas, southwestern Missouri, and northwestern Arkansas and for the part of the study area in Oklahoma, 1980–2017.

Water-budget category	Study area						Study area in Oklahoma
	Boone aquifer		Roubidoux aquifer		Net total		Boone aquifer		Roubidoux aquifer		Net total
Recharge	2,520,887	(95.68)	53,389	(5.53)	2,574,276	(93.53)	1,420,898	(100)	6,881	(6.23)	1,427,779	(100)
Seepage	−1,896,569	(75.13)	−420,977	(43.45)	−2,317,546	(89.8)	−1,219,919	(85.75)	−95,821	(85.67)	−1,315,739	(91.95)
Leakage1	−734,385	(29.09)	912,606	(94.47)	178,221	(6.47)	−199,075	(13.99)	103,529	(−92.56)	−95,546	(6.68)
Lateral groundwater flow across boundary	113,875	(4.32)	−525,899	(54.28)	−412,025	(15.96)	−1,051	(−0.07)	−10,416	(−9.43)	−11,467	(−0.8)
Groundwater use	−7,470	(0.3)	−22,003	(2.27)	−29,472	(1.14)	−2,613	(0.18)	−5,616	(5.02)	−8,229	(0.58)
Net change in groundwater storage	−3,125	(0.12)	−2,884	(0.3)	−6,008	(0.23)	−1,234	(0.09)	−1,443	(1.29)	−2,677	(0.19)
Total inflow	2,634,762		965,995		2,752,497		1,420,898		110,409		1,427,779
Total outflow	−2,524,549		−968,879		−2,580,822		−1,422,657		−111,852		−1,430,981

¹

Leakage consists of the sum of leakage to the Western Interior Plains confining unit and Ozark confining unit.

For the Boone aquifer, most of the component flows of the conceptual model were similar to those of the simulated water budget (tables 6 and 13). Simulated recharge was about 3.5 percent (study area) and 1.3 percent (Oklahoma) less than conceptual model recharge (tables 6 and 13). Simulated seepage was about 20 percent (study area) and 56 percent (Oklahoma) greater than conceptual model seepage (tables 6 and 13). Simulated Boone aquifer leakage was about 251 percent greater for the simulation compared to the conceptual model for the study area (tables 6 and 13). Leakage was largely unknown and was specified as a mean leakage estimated for the entire Ozark Plateaus aquifer system. Simulated leakage from the Ozark confining unit, when normalized to each subarea, was still within the leakance bounds (multiplied by the volume of the Ozark confining unit) estimated by Christenson and others (1990) in Ottawa County, Okla. (table 5; Trevisan and others, 2024).

The components of the conceptual model did not match the simulated water budget as closely for the Roubidoux aquifer compared to the Boone aquifer (tables 6 and 13). Most of the components of the conceptual model for the Roubidoux aquifer were poorly represented by available field data. Simulated recharge, however, did match closely to that of the conceptual model. Seepage was estimated by using base-flow estimates from streamgages. Most streamgages drained either only the Boone aquifer or mostly the Boone aquifer; therefore, conceptual-model seepage estimates might be more representative of the Boone aquifer than the Roubidoux aquifer. Seeps also tend to occur in areas of the aquifer with more conduit development (Imes and Emmett, 1994), which would result in greater groundwater discharge.

Simulated leakage was about 251 percent and 336 percent greater than conceptual model leakage within the study area for the Boone and Roubidoux aquifers, respectively, but only about 73 percent greater and 10 percent less than conceptual model leakage within the Oklahoma part of the study area for the Boone and Roubidoux aquifers, respectively (tables 6 and 13). Leakage was largely unknown and specified as a mean leakage estimated for the entire Ozark Plateaus aquifer system. Simulated leakage from the Ozark confining unit, when normalized to each subarea, was still within the leakance bounds (multiplied by the volume of the Ozark confining unit) estimated by Christenson and others (1990) in Ottawa County, Okla. (table 5; Trevisan and others, 2024).

Simulated lateral groundwater flow was greater than conceptual model lateral groundwater flow by about 251 percent and 336 percent within the study area for the Boone and Roubidoux aquifers, respectively, and by about 99 and 90 percent within Oklahoma for the Boone and Roubidoux aquifers, respectively (tables 6 and 13). Lateral groundwater flow was not estimated by field data but was estimated by balancing groundwater flow in the conceptual model. For the study area, increased simulated leakage (compared to the conceptual model) accounts for the increased inflows, and thus, increased simulated seepage and lateral groundwater flow across the study area and Oklahoma boundaries. Changes in groundwater storage represent a larger component of the simulated groundwater budget for the Roubidoux aquifer than for the Boone aquifer (fig. 33A–D).

Simulated annual net groundwater storage changes reflect climatic variability during the transient-simulation model period 1980–2017 (fig. 33). Periods with loss in groundwater storage generally occur during years with below normal precipitation, and periods with gain in groundwater storage generally occur during years with above normal precipitation for both aquifers (fig. 33; fig. 8E). For both aquifers, the largest declines in storage occurred during 1980–83, 2005–07, and 2011–14, which were drier periods (fig. 3A). Cumulative net change in storage during the 10-year period 1990–99 was slightly positive for the Boone and Roubidoux aquifers (fig. 33).

Simulated Saturated Thickness and Groundwater Storage Availability

Simulated saturated thicknesses were calculated at the end of December 2017 for the Boone and Roubidoux aquifers. Patterns in simulated saturated thickness were generally similar to patterns in aquifer thickness. For example, the Boone and Roubidoux aquifers thicken towards the southern part of the study area where the saturated thicknesses were often greatest (fig. 34; Russell, 2020). The thinnest part of each aquifer with the smallest saturated thicknesses were generally in the northwestern part of the study area (fig. 34; Russell, 2020).

Boone aquifer saturated thickness ranged from 0 to 1,373 ft for the study area (fig. 34A). Areas of the Boone aquifer with the highest saturated thickness were mostly where the aquifer was confined (mostly 200 to 300 ft of saturated thickness; fig. 34A). In the unconfined part of the Boone aquifer, saturated thickness was often less than 200 ft (fig. 34A). Saturated thickness was generally less around streams and lakes except for streams and lakes in the northwestern part of the study area (fig. 34A). Saturated thickness greater than 650 ft was present where the Western Interior Plains confining unit has eroded and were also present in the thickest part of the Boone aquifer in the southeastern part of the study area. Saturated thickness was also less in an area within Craig County (where the aquifer is relatively thin in the northwestern part of the study area).

For the Roubidoux aquifer, saturated thickness ranged from approximately 77 to 2,701 ft for the study area and the Oklahoma part of the study area (fig. 34B). Saturated thickness was relatively low in northeastern Cherokee County and northwestern Adair County. From about central Craig County to northwestern Ottawa County, saturated thickness was the lowest (500 ft or less).

Mean saturated thickness and transmissivity were calculated for the Boone and Roubidoux aquifers based on each active cell within Oklahoma (table 12). Mean transmissivity for the Roubidoux aquifer (7,591 ft²/d) was approximately an order of magnitude greater than transmissivity for the Boone aquifer (902 ft²/d), largely because the Roubidoux aquifer simulated saturated thickness (mean of 1,487 ft) was much thicker than that of the Boone aquifer (mean of 224 ft; table 12). The differences in saturated thickness between the Boone and Roubidoux aquifers resulted in greater simulated groundwater storage for the Roubidoux aquifer (50,226,846 acre-ft) than for the Boone aquifer (9,656,595 acre-ft) despite higher calibrated specific yield and specific storage for the Boone aquifer (table 12).

Boone Aquifer

The EPS simulations for the Boone aquifer were conducted by using the PESTPP-GLM program (White and others, 2020b). The PESTPP-GLM program iteratively adjusted the groundwater withdrawal rate in the WEL package, conducted the simulation, and recorded the percentage of aquifer with saturated thickness less than 15 ft after solving for the EPS target after 20, 40, and 50 years. The 20-, 40-, and 50-year EPS rates for the Boone aquifer were 1.10, 0.98, and 0.96 acre-feet per acre per year ([acre-ft/acre]/yr), respectively, when simulating normal recharge conditions (that is, recharge equal to the long-term [1980–2017] mean annual SWB recharge for the study area; table 14). Given the 2,799,816 simulated acres for the Boone aquifer in Oklahoma, these rates correspond to annual yields of approximately 3,071,199; 2,749,343; and 2,674,636 acre-ft/yr for the 20-, 40-, and 50-year normal recharge scenarios, respectively. These rates exceed the reported 2018 groundwater use for the Boone aquifer (fig. 8A) by about three orders of magnitude. Decreasing and increasing recharge by 10 percent resulted in a change in the EPS rate of about 7–8 percent. At the end of the 20-year EPS simulation for the Boone aquifer (2037), simulated saturated thickness was lowest in the central to northern part of the study area (fig. 35). Streamflows were mostly depleted at the end of the 20-year EPS simulation for the Boone aquifer using normal recharge (fig. 35). Base flows at the end of the scenario were highest in streams with simulated inflows (Neosho and Spring Rivers) and the Elk River (in the northern part of the study area; fig. 1). Base flows were less than 50 ft³/s for other streams in the study area. Because the mean SWB recharge for 1980–2017 was used for the EPS simulation, the Boone aquifer was near equilibrium with EPS groundwater withdrawals and mean recharge because the change in groundwater storage was minimal between each of the final simulated years (fig. 36).

Roubidoux Aquifer

The EPS simulations for the Roubidoux aquifer required large groundwater withdrawal rates that rapidly dewatered many parts of the aquifer and caused simulation instability, resulting in the EPS target not being achieved. Instead, the EPS simulations for the Roubidoux aquifer were conducted by gradually increasing the EPS rates and graphing the percentage of the aquifer with less than 15 ft of saturated thickness for many simulations (fig. 37). These simulations were conducted by using the PEST++ parameter sweep program (PESTPP-SWP; White and others, 2020b) on an HPC cluster. The PESTPP-SWP program runs multiple simulations by adjusting parameters based on user-specified parameters. After each simulation, the simulated EPS rate and the percentage of the aquifer with less than 15 ft of saturated thickness (hereinafter referred to as the “PCT15”) were recorded as observations.

Another observation was an EPS rate (hereinafter referred to as the “leakage-adjusted EPS rate”) that was estimated by subtracting leakage from the Boone aquifer from the amount of simulated groundwater withdrawn from the Roubidoux aquifer. The leakage-adjusted EPS rate represents an EPS rate that assumes no leakage from the Boone aquifer and was determined by finding the maximum simulated (actual) EPS rate after correcting for leakage. Attempts were made to project the EPS rates to the EPS target by fitting exponential regressions. The actual EPS rates, which were calculated by dividing the total volume of groundwater withdrawn by the number of years and the area of the aquifer, plateaued (and declined slightly) even as the applied EPS rate (the rates applied in the WEL package) increased because of a reduction in groundwater availability and the smoothing conducted by the MODFLOW-NWT WEL package (fig. 37).

The Boone aquifer was mostly drained (PCT15 of about 99 percent for the maximum simulated EPS rate; Trevisan and others, 2024), resulting in small inflows to the Roubidoux aquifer. Depleting the Boone aquifer also caused most streams to transition to no flow for the 20-year EPS groundwater-withdrawal scenario for the Roubidoux aquifer (fig. 38); most streams also transitioned to no flow for the 20-year EPS groundwater-withdrawal scenario for the Boone aquifer (fig. 35). The EPS rate projected with exponential regressions was within 5 percent of this maximum simulated rate; however, the regressions were consistently overpredicting the EPS rate based on the observations for the simulations using higher simulated EPS rates. Because exponential regression equations poorly fit the data and the amount of groundwater withdrawn did not increase when the applied EPS rate was increased for larger applied EPS rates, the EPS rate was estimated from the maximum actual EPS rate (fig. 37A–I).

The 20-, 40-, and 50-year EPS groundwater withdrawal rates for the Roubidoux aquifer were 1.76, 1.34, and 1.25 (acre-ft/acre)/yr, respectively, when simulating normal recharge (table 15). Given the 2,937,557 simulated acres for the Roubidoux aquifer in Oklahoma, these rates correspond to annual yields of about 5,155,653; 3,930,404; and 3,667,315 acre-ft/yr, for the 20-, 40-, and 50-year normal recharge scenarios, respectively. The 20-, 40-, and 50-year leakage-adjusted EPS rates were 0.60, 0.32, and 0.26 (acre-ft/acre)/yr, respectively, when simulating normal recharge (table 15). For the Roubidoux aquifer in Oklahoma, the leakage-adjusted EPS rates corresponded to annual yields of about 1,759,718; 932,514; and 750,852 acre-ft/yr, respectively. These rates were greater than 2018 groundwater use by about two to three orders of magnitude (fig. 8B). Decreasing and increasing recharge by 10 percent resulted in a change in the EPS rate of about 6–7 percent and almost no change in the leakage-adjusted EPS rate. The EPS rates for the Roubidoux aquifer are likely less sensitive to recharge than the EPS rates for the Boone aquifer because the Roubidoux aquifer is mostly confined and recharge is a smaller component of the water budget for the Roubidoux aquifer despite the generally greater Roubidoux aquifer thickness compared to that of the Boone aquifer (table 4).

By using normal recharge and the highest applied 20-year EPS groundwater withdrawal rate for the Roubidoux aquifer, saturated thickness and base flow were mapped to visualize the conditions of the aquifer after the simulation (fig. 38). The highest applied 20-year EPS groundwater withdrawal rate simulation was used because the greatest reduction in groundwater altitudes would result from the simulation using this rate. Although the EPS groundwater withdrawal target rate was small, the aquifer was drained rapidly, a larger percentage (likely more than 50 percent) of the aquifer contained 50 ft of saturated thickness or less, and base flow was depleted from most of the streams (fig. 38). The EPS groundwater withdrawal rates which could produce the EPS groundwater withdrawal target would likely be close to the highest applied 20-year EPS groundwater withdrawal rate. Streams that were not dry were in the northern part of the model that contained SFR inflows; SFR inflows were likely sustaining base flows and saturated thickness in these locations. The highest actual EPS rates were selected as the EPS rate for each 20-, 40-, and 50-year EPS and simulated using normal recharge to illustrate changes in groundwater storage for the Roubidoux aquifer (fig. 39). Similar to the Boone aquifer (fig. 36), the Roubidoux aquifer was potentially in equilibrium with mean 1980–2017 recharge because decreases in groundwater storage were near zero at the ends of the EPS simulations (fig. 39). Cumulative percent changes in groundwater storage were slightly less than those for the Boone aquifer EPS simulations and could potentially be greater if the EPS simulations could be simulated to a PCT15 of 50 percent.

Comparison of Equal-Proportionate-Share and 2018 Groundwater-Withdrawal Rates

EPS simulations are a theoretical construct assuming the entire aquifer is fully developed with regularly spaced wells (one well per cell or one well per approximately 91.8 acres of simulated area), each withdrawing groundwater at rates often more than 100 times greater than 2018 reported groundwater use. Some parts of the aquifers may be more developed than others, but reported groundwater use for each aquifer was far less than groundwater use simulated in the EPS simulations. In Oklahoma, reported groundwater use for 2018 was approximately 2,279.87 acre-ft/yr and 6,136.16 acre-ft/yr for the Boone and Roubidoux aquifers, respectively (fig. 8A, B; Trevisan and others, 2024); 2018 reported groundwater use was equivalent to groundwater-withdrawal rates of approximately 0.0008 and 0.002 (acre-ft/acre)/yr for the Boone and Roubidoux aquifers, respectively, when equally distributed over the 2,799,816- and 2,937,557-acre spatial extents of the Boone and Roubidoux aquifers, respectively. Both rates are approximately 1 percent or less of the EPS rates estimated in this report.

Boone Aquifer

For the Boone aquifer, groundwater storage after 50 years with no groundwater withdrawals was 9,812,658 acre-ft, or 5,993 acre-ft (0.06 percent) greater than the groundwater storage at the end of the 2017-groundwater-withdrawal simulation (table 16). Groundwater storage at the end of the 50-year period with the increasing demand groundwater withdrawals was 9,803,276 acre-ft, or 3,389 acre-ft (0.03 percent) less than the storage at the end of the 2017-groundwater-withdrawal scenario. Groundwater storage after 50 years of groundwater withdrawals at the mean groundwater withdrawal for the 1980–2017 study period was 9,807,904 acre-ft, or 1,239 acre-ft (0.01 percent) greater than the groundwater storage at the end of the 2017-groundwater-withdrawal scenario.

Roubidoux Aquifer

For the Roubidoux aquifer, groundwater storage after 50 years with no groundwater withdrawals was 50,309,351 acre-ft, or 6,820 acre-ft (less than 0.01 percent) greater than the groundwater storage at the end of the 2017-groundwater-withdrawal scenario (table 16). Groundwater storage at the end of the 50-year period with the increasing demand groundwater withdrawals was 50,298,792 acre-ft, or 3,739 acre-ft (less than 0.01 percent) less than the storage at the end of the 2017-groundwater-withdrawal scenario. Groundwater storage after 50 years of groundwater withdrawals at the mean rate for the 1980–2017 study period was 50,303,211 acre-ft, or 680 acre-ft (less than 0.01 percent) greater than the groundwater storage at the end of the 2017-groundwater-withdrawal scenario.