Purpose and Scope

The purpose of this report is to describe a hydrologic investigation of the Washita River aquifer in southern Oklahoma. This description includes (1) an updated summary of the hydrogeology with a definition of the hydrogeologic framework (including an updated geographic extent) of the aquifer, (2) development of conceptual and calibrated numerical groundwater-flow models for the aquifer representing the 1980–2017 study period, and (3) results of simulations of groundwater-availability scenarios. The groundwater-availability scenarios use the calibrated numerical groundwater-flow model to (1) estimate the EPS pumping rate that ensures a minimum 20-, 40-, and 50-year life of the aquifer, (2) quantify the potential effects of projected well withdrawals on groundwater storage over a 50-year period, and (3) simulate the potential effects of a hypothetical (10-year) drought on groundwater storage. The calibrated numerical groundwater-flow model and groundwater-availability scenarios were archived and released in a USGS data release (Rogers and others, 2023).

The entire Washita River aquifer in Oklahoma consists of four administrative sections, or reaches, which are managed independently by the OWRB. Reach 1 extends from the Texas border to downstream from Foss Reservoir in western Oklahoma (fig. 1; Ellis and others, 2020), and reach 2 is hydraulically connected to and managed as part of the Rush Springs aquifer in western Oklahoma (fig. 1; Ellis, 2018a); both reach 1 and reach 2 are outside the scope of this report. The modeled area of this report is a 279.6-square-mile (mi²) (178,938-acre) subarea of the 411.5-mi² extent of alluvium and terrace deposits (Heran and others, 2003; Miller and Stanley, 2004; Chang and Stanley, 2010; Stanley and Chang, 2012) of reaches 3 and 4 of the Washita River aquifer in southern Oklahoma (fig. 1). Reach 3 extends from near Anadarko, Okla., to Alex, Okla., and includes 113.6 mi² (72,697 acres, or 40.6 percent) of the modeled area. Reach 4 extends from near Alex to south of Davis, Okla., and includes 166.0 mi² (106,241 acres, or 59.4 percent) of the modeled area. The aquifer becomes narrow and thin in smaller tributary valleys, so this investigation focused on the main body of the aquifer along the Washita River main stem and selected large tributaries (modeled aquifer area, fig. 1). Though sometimes referred to as the “Washita River alluvial aquifer” (Schipper, 1983; OWRB, 2019b), the aquifer is referred to as the “Washita River aquifer” in this report because it consists not only of alluvium deposits but also of terrace and dune deposits. Selected sections in this report are modified from Ellis (2018a, b), Ellis and others (2020), and Smith and others (2021). Though the study areas are different, the organization and wording of this report are largely based on Smith and others (2021).

Description of Study Area

The Washita River aquifer is a long, narrow, unconfined aquifer that includes alluvium and terrace deposits of the Washita River and tributaries in Caddo, Grady, Garvin, McClain, and Murray Counties, Oklahoma (fig. 1). The Washita River is perennial in the study area except during the most extreme droughts as documented by historical USGS streamgaging records for the stream obtained from the USGS National Water Information System (NWIS) database (USGS, 2020; table 1). The Washita River generally flows from northwest to southeast for about 207 miles (mi) (Horizon Systems Corporation, 2015) in the study area. The largest tributaries in the study area are Little Washita River and Sugar, Rush, and Wildhorse Creeks (fig. 1). Many of the smaller tributaries to the Washita River in the study area are dammed by hundreds of Natural Resources Conservation Service (NRCS) floodwater-retarding structures (fig. 2), which were mostly built in the 1950s through 1970s (U.S. Army Corps of Engineers, 2021). Floodwater-retarding structures were designed to impound runoff from headwater catchments and slowly release it downstream with the purpose of reducing the frequency and magnitude of damaging flood flows on the Washita River main stem (Smith and Esralew, 2010). The Washita River main stem, however, is not dammed in the study area.

Land-cover data for the Washita River aquifer study area were obtained from the CropScape database (National Agricultural Statistics Service, 2019) (fig. 3), which includes land-cover characteristics at 30-meter (m) resolution for land overlying the Washita River aquifer for the 2010–17 period. During this period, land-cover type was primarily cropland (44.1 percent), grass or pasture (36.3 percent), forest and shrubland (9.3 percent), or developed (7.7 percent). Winter wheat was the major crop-cover type in the study area, accounting for 57.6 percent of cropland by area. Alfalfa or hay (22.0 percent) was the next largest crop-cover type by area. Corn and soybeans accounted for 7.5 and 6.6 percent of cropland by area, respectively. Cotton (1.3 percent) and fallow or idle cropland (1.2 percent) were the only other crop-cover types that accounted for more than 1 percent of cropland by area (fig. 3). Crop-cover types typically change in response to economic conditions and hydrologic factors, but the percentages of total cropland cover and individual crop types did not change substantially during the 2010–17 period.

Climate Characteristics

The climate of the study area is classified as subhumid with generally southerly prevailing winds (Patterson, 1984). Historical data from selected climate stations in southern Oklahoma (Climate Division 8) have been quality assured and summarized monthly as part of the Global Historical Climatology Network (National Centers for Environmental Information, 2021). These data were used to calculate and graph annual and monthly temperature and precipitation statistics for the study area (fig. 4, table 2). A locally weighted scatterplot smoothing line (Cleveland, 1979) was used to delineate periods of below- and above-mean annual temperature and precipitation (fig. 4). The mean annual temperature during 1895–2020 was 62.2 degrees Fahrenheit (°F) (fig. 4, table 2), and the mean annual temperature differed by about 2 °F from northwest to southeast across the study area (Oklahoma Climatological Survey, 2021a, b). The mean annual precipitation during 1895–2020 was 38.0 inches per year (in/yr) (fig. 4, table 2), and the mean annual precipitation increased by about 7.5 in/yr from northwest to southeast across the study area (Oklahoma Climatological Survey, 2021a, b). The mean annual precipitation for the 1980–2017 study period was 40.1 in/yr (fig. 4, table 2). The 1981–2002 period was an unprecedented wet period in the available record; above-mean precipitation was recorded in 18 of these 22 years (82 percent of the years), contributing to a higher mean precipitation during the study period compared to the mean for the period of record (fig. 4A). May is typically the wettest month, and January is typically the driest month (fig. 5A).

Multiyear to decadal droughts are common for the study area (fig. 4). The 1929–41 (Dust Bowl; Egan, 2006) and 1952–56 drought periods were among the most severe in Oklahoma in the 20th century; two shorter, less severe drought periods also occurred in the late 20th century during 1961–72 and 1976–81 (Tortorelli, 2008; Shivers and Andrews, 2013). The 21st century began with the 2002–6 and 2010–14 drought periods (Tortorelli, 2008; Shivers and Andrews, 2013). The most severe droughts on record developed from extended periods of below-mean precipitation paired with above-mean temperature.

Streamflow Characteristics and Trends

Daily streamflow data, with various and sometimes interrupted periods of record, were recorded at selected USGS streamgages in the study area (figs. 1 and 2) and summarized for the 1980–2017 study period (table 3). Streamflow measured at streamgages is the sum of runoff and base flow originating upstream; base flow is the component of streamflow that is supplied by the discharge of groundwater to streams (Barlow and Leake, 2012). For this report, streamflow-hydrograph data obtained from the NWIS database (USGS, 2020) were separated into runoff and base-flow components by using the standard Base-Flow Index code (Wahl and Wahl, 1995) in the USGS Groundwater Toolbox (Barlow and others, 2015). The Base-Flow Index code uses the minimum streamflow in a moving n-day window as a basis for hydrograph separation; for consistency, a 6-day window was used for all streamgages in this report. The 6-day window was selected by testing multiple n-day windows for the study period at USGS streamgage 07328100 Washita River at Alex, Okla. (hereinafter referred to as the “Alex gage”) and plotting the resulting mean base-flow percentage against n; a slope change was evident at 6 days. Base flow, as computed using this method, accounts for about 35–48 percent of annual streamflow for the period of record at streamgages on the Washita River main stem (table 3). This percentage, known as the base-flow index (BFI) (Wahl and Wahl, 1995; Barlow and others, 2015), generally increased over the period of record at streamgages on the Washita River (fig. 6A–D). When summarized annually using the Kendall (1938) tau test with the Theil-Sen (Sen, 1968) slope estimator and an alpha value of 0.0500 set as the significance level, the upward trends in the BFI were found to be statistically significant (probability value [p-value] = 0.0054) at USGS streamgage 07326500 Washita River at Anadarko, Okla. (hereinafter referred to as the “Anadarko gage”), statistically significant (p-value = 0.0248) at the Alex gage, statistically significant (p-value = 0.0000) at USGS streamgage 07328500 Washita River near Pauls Valley, Okla. (hereinafter referred to as the “Pauls Valley gage”), and statistically significant (p-value = 0.0083) at USGS streamgage 07331000 Washita River near Dickson, Okla. (hereinafter referred to as the “Dickson gage”). Upward trends in base flow (fig. 6A–D) were found to be not statistically significant (p-value = 0.0772) at the Anadarko gage, statistically significant (p-value = 0.0111) at the Alex gage, statistically significant (p-value = 0.0163) at the Pauls Valley gage, and statistically significant (p-value = 0.0177) at the Dickson gage. The overall pattern of upward trends in the BFI and base flow may have developed following the construction of more than 570 NRCS floodwater-retarding structures (dams; fig. 2) on tributaries of the Washita River within the study area. NRCS dam construction in this area began in 1948, peaked at 65 dams per year in 1963, greatly diminished to about 1 dam per year in the 1980s, and concluded in the early 1990s (U.S. Army Corps of Engineers, 2021). NRCS dams temporarily impound runoff and release it over several days with the purpose of reducing peak flows on the river main stem. The flows that are slowly and steadily released by these dams are indistinguishable from groundwater-derived base flows when they reach the main stem and are mostly classified as base flows by the Base-Flow Index code. Though most upward trends in the BFI and base flow were statistically significant for the streamgage periods of record, no statistically significant trends in base flow or BFI were detected for the 1980–2017 study period at streamgages on the Washita River main stem. One interpretation of the upward trends in the BFI and base flow for the period of record is that the construction of the floodwater-retarding structures caused a step change in apparent base flows in the Washita River.

Surface-water use in the Washita River aquifer study area consists of permitted diversions along the main stem of the Washita River and its associated tributaries (OWRB, 2019d). At the end of the study period (2017), 60 active, nontemporary permitted diversions on the Washita River main stem allowed for a maximum total withdrawal of 13,957 acre-ft/yr. All of those permitted diversions were for irrigation with the exception of three permitted diversions that were used for public supply, and all but 15 of those permitted diversions were in reach 3 of the Washita River aquifer. At the end of the study period (2017), 29 active, nontemporary permitted diversions on Washita River tributaries allowed for a maximum total withdrawal of 9,749 acre-ft/yr.

Groundwater Use

The OWRB permits and regulates groundwater use of more than 5 acre-ft/yr for domestic and agricultural purposes and groundwater use for irrigating more than 3 acres of land for growing gardens, orchards, or lawns in Oklahoma (Oklahoma Statutes §82–1020.1[2] [Oklahoma State Legislature, 2021a; OWRB, 2014]; Oklahoma Statutes §82–1020.3 [Oklahoma State Legislature, 2021c]). Groundwater-use data since 1980 are self-reported annually to the OWRB by permitted users (Christopher Neel, OWRB, written commun., 2019). In 2017, about 61 long-term temporary (regular) groundwater-use permits and about 66 prior-right groundwater-use permits were active for the Washita River aquifer (OWRB, 2019d). Each permit may be served by multiple wells that share the allocated groundwater use. Most groundwater-use permits for the Washita River aquifer were allocated for irrigation (38.0 percent) and public supply (61.4 percent), with some allocated for other uses (0.6 percent), including agricultural (nonirrigation); commercial; industrial; mining; power; recreation, fish, and wildlife; or uncategorized uses (figs. 7 and 8, table 4) (OWRB, 2019d; Rogers and others, 2023). Groundwater use for domestic supply (self-supplied directly to a home by a private well) was assumed to be a negligible part of the total groundwater use because the population of the land overlying the aquifer was predominantly rural or concentrated in cities supplied water by municipal or rural water districts. For the purposes of this report, all groundwater use was assumed to be consumptive use (that is, none of the groundwater that is withdrawn returns to the aquifer as recharge).

Mean annual groundwater use for the 1980–2017 study period was 1,214.1 acre-ft/yr, which includes 745.8 acre-ft/yr (61.4 percent) for public supply and 461.8 acre-ft/yr (38.0 percent) for irrigation (table 4) (Rogers and others, 2023). Mean annual groundwater use for the period of record (1967–2017) was 1,590.7 acre-ft/yr, which includes 861.6 acre-ft/yr (54.2 percent) for irrigation and 709.7 acre-ft/yr (44.6 percent) for public supply (table 4) (Rogers and others, 2023). The decrease in groundwater use, particularly in irrigation groundwater use after 1980, is likely a result of the OWRB implementing changes in methods of groundwater-use reporting. Since 1980, the OWRB has required irrigation permit holders to estimate and annually report their groundwater use in terms of number of applications and number of inches of water applied (Christopher Neel, OWRB, written commun., 2019). Prior to 1980, however, inches of water applied during irrigation were not required to be reported. As a result, the OWRB adopted rules to estimate the number of inches of water applied for pre-1980 data on the basis of the number of water applications. This estimate results in what appears to be a step-change decrease in irrigation groundwater use after 1980 (fig. 8). In 2017, regular permits accounted for about 58 percent and prior-right permits accounted for about 42 percent of reported groundwater use from the Washita River aquifer (OWRB, 2019d).

Mean annual groundwater use for the 1980–2017 study period was 309.9 acre-ft/yr in reach 3 and 904.2 acre-ft/yr in reach 4 of the Washita River aquifer (table 4) (Rogers and others, 2023). Mean annual groundwater use for the period of record 1967–2017 was 406.7 acre-ft/yr in reach 3 and 1,184.0 acre-ft/yr in reach 4 of the Washita River aquifer (table 4) (Rogers and others, 2023). Groundwater-use permits for public supply were predominantly for wells near Lindsay, Okla., or Pauls Valley, Okla., in reach 4 of the Washita River aquifer, whereas groundwater-use permits for irrigation were for wells scattered throughout the Washita River aquifer (fig. 7). Reach 4 accounted for 74.5 percent of the groundwater use for the 1980–2017 study period and 74.4 percent of the groundwater use for the period of record 1967–2017.

Historical groundwater yields from wells completed in the Washita River aquifer generally decreased from northwest to southeast as the aquifer materials thin; groundwater yields of about 350 gallons per minute (gal/min) were reported near Lindsay and groundwater yields of about 200 gal/min were reported near Davis in the 1950s (Leonard and others, 1958). Driller-estimated groundwater-yield values, which have been mostly recorded since the 1980s, averaged about 70 gal/min; about 50 percent of the values recorded were greater than 25 gal/min (the median value), and about 9 percent of the values recorded were greater than 200 gal/min (OWRB, 2019a). Historical groundwater-yield values in the Washita River aquifer obtained from the NWIS database (USGS, 2020) averaged about 240 gal/min. The median groundwater-yield value was 175 gal/min, and a groundwater-yield value of 1,200 gal/min was reported at one well northwest of Pauls Valley (USGS station 344625097181001). These groundwater-yield statistics likely were biased high because high-yield irrigation and public-supply wells were the most common type of well to have measured discharges (USGS, 2020).

Alluvium and Terrace Deposits

Alluvium and terrace deposits of Quaternary age are the principal water-bearing units of the Washita River aquifer (fig. 9). These deposits are made up of alluvial gravel, sand, silt, and clay (figs. 9 and 10; Leonard and others, 1958). Most terrace deposits are hydraulically connected to the valley-bottom alluvium, but some are hydraulically disconnected from the other terrace deposits and the valley-bottom alluvium and thus are unable to store water (Leonard and others, 1958). The disconnected terrace deposits are most evident north of Anadarko and northwest of Pauls Valley (fig. 9). The alluvium and terrace deposits are typically about 2–3 mi wide, but where the Washita River Valley constricts, those deposits are less than 1 mi wide (fig. 11). The Washita River Valley is narrowest where it approaches the Arbuckle Uplift (figs. 9 and 11E). In that area, the alluvium and terrace deposits are only about 0.25 mi wide, which the OWRB regards as the southern boundary of the Washita River aquifer (fig. 9; OWRB, 2019a).

Bedrock Units

Bedrock units of Permian age underlie all of reach 3 (fig. 9A) and most of reach 4 (fig. 9B) of the Washita River aquifer. In the northwesternmost part of the aquifer, west of Anadarko, the Whitehorse Group directly underlies the geologic units that contain the Washita River aquifer (fig. 9). The Whitehorse Group consists of the Rush Springs Formation and the Marlow Formation (fig. 10) but is shown in figure 9 as an undivided unit. The Rush Springs Formation is described as orange-brown, cross-bedded fine-grained sandstone, with some dolomite and gypsum beds (Ellis, 2018a; Neel and others, 2018). The thickness of this geologic unit ranges between 186 and 560 ft in the study area (Ellis, 2018a). Stratigraphically below the Rush Springs Formation is the Marlow Formation, which consists of orange-brown, fine-grained sandstone and siltstone interbedded with mudstone, gypsum, and dolomite (Kent and others, 1984; Ellis, 2018a). The Marlow Formation directly underlies the Washita River aquifer near Anadarko and along Sugar Creek (Kent and others, 1984).

Beneath the White Horse Group is the Permian-age El Reno Group, which consists of the Dog Creek Shale, the Chickasha Formation, and the Duncan Sandstone (fig. 10) but is shown in figure 9 as an undivided unit (Heran and others, 2003). The Dog Creek Shale underlies the Washita River aquifer in the western part of Grady County and consists of red, brown, and green gypsiferous shales interbedded with siltstone, sandstone, and dolomite (Kent and others, 1984; Heran and others, 2003; Ellis, 2018a). The Dog Creek Shale has an estimated thickness of 200 ft and is described as having very low permeability (Mogg and others, 1960; Ellis, 2018a), which serves to limit downward flow in the Washita River aquifer. The Chickasha Formation and Duncan Sandstone make up the stratigraphically lower part of the El Reno Group. The Chickasha Formation is characterized as a red-brown mudstone conglomerate with some shale and siltstone and is estimated to be about 200 ft thick (Kent and others, 1984; Ellis and others, 2017). The Duncan Sandstone has a similar thickness to the Chickasha Formation and the Dog Creek Shale and is characterized as a light gray to reddish cross-bedded fine-grained sandstone interbedded with shale (Heran and others, 2003; Ellis and others, 2017). The Duncan Sandstone directly underlies the Washita River aquifer as reach 3 transitions to reach 4 (Heran and others, 2003). The units of the El Reno Group contain the El Reno minor aquifer in the study area.

The Permian-age Hennessey Group directly underlies the Washita River aquifer in the northern part of Garvin County (Kent and others, 1984; fig. 9). The Hennessey Group is made up of the Bison Formation, Purcell Sandstone, and Fairmont Shale (fig. 10) but is shown in figure 9 as an undivided unit. The Bison Formation, estimated to be 50–90 ft thick in the study area, makes up the upper part of the Hennessey Group and is characterized as gray to red-brown shale containing calcium carbonate (Heran and others, 2003). The Purcell Sandstone is a red-brown, fine- to coarse-grained sandstone with some shale and mudstone conglomerate (Heran and others, 2003). The Purcell Sandstone is about 90–150 ft thick in the study area and lies directly above the Fairmont Shale. The Fairmont Shale is similar in color to the Purcell Sandstone, characterized as red-brown and blocky. The Fairmont Shale thickness is 40–80 ft in the study area.

The Permian-age Sumner Group, which includes the Garber Sandstone and Wellington Formation, underlies the Washita River aquifer northwest of Pauls Valley (figs. 9 and 10). The Garber Sandstone is characterized as red-brown fine-grained sandstone with a thickness of about 100–150 ft in the study area (Heran and others, 2003). The Wellington Formation is composed of red-brown shale with several 20–30-ft bituminous sandstones at the base and a thickness of about 100–200 ft decreasing southeastward in the study area (Heran and others, 2003). The Garber Sandstone and Wellington Formation contain the Garber-Wellington aquifer (also known as the Central Oklahoma aquifer, Mashburn and others, 2013) in the northeastern part of the study area, where the Garber Sandstone and Wellington Formation are thicker and more permeable. However, the Garber Sandstone and Wellington Formation are not considered to be a major aquifer near the Washita River aquifer (Hart, 1974; Patterson, 1984).

Underlying the Sumner Group is the Pontotoc Group, a mostly Permian-age red-brown to grayish shale with arkosic sandstone and limestone conglomerates (fig. 10). The Pontotoc Group directly underlies the Washita River aquifer in reach 4 where it approaches the Arbuckle Mountains (fig. 9) and has a relative thickness of 300–500 ft in the study area (Heran and others, 2003). The Pontotoc Group is not sufficiently permeable to yield large quantities of water in most locations and is not considered to be a major water-bearing unit in the study area; however, well yields of 60 gal/min have been reported in a small area in southwestern Garvin County (Hart, 1974).

Pre-Permian Paleozoic-age units make up the Arbuckle Mountains and underlie the southeasternmost part of the Washita River aquifer (fig. 9). These units are composed of gray-brown limestone with shale and some bituminous sandstone and limestone conglomerate (fig. 10) and were greatly faulted and folded by the Arbuckle Uplift (Reeds, 1927; Ye and others, 1996). Some of the pre-Permian Paleozoic-age units form the Arbuckle-Simpson aquifer in southern Oklahoma. It is uncertain, however, whether that aquifer is in direct hydraulic connection with the Washita River aquifer as defined in this report. Igneous units (Wichita Granite and Raggedy Mountain Gabbro Groups [fig. 9], consisting of diallage and olivine gabbro, anorthosite, perthite leucogranites, and diorite composing a layered intrusion [fig. 10]) are present within the study area southwest of Davis; however, these igneous units are not in direct hydraulic connection with the Washita River aquifer.

Groundwater Quality

Groundwater-quality data were collected in the study area during July 28–August 25, 2014, by the OWRB (OWRB, 2017). Groundwater was sampled from 14 wells in the alluvium and terrace deposits of the Washita River aquifer study area (fig. 9; site information is available in Rogers and others [2023]) and analyzed for selected constituents including physical properties (pH, temperature, and specific conductance), major ions, nutrients, and trace metals. Total dissolved solids (TDS) concentrations in the groundwater samples ranged from 215 to 1,070 milligrams per liter (mg/L) with a median concentration of 523 mg/L and a mean concentration of 594 mg/L. The U.S. Environmental Protection Agency (2017) has established a secondary drinking water standard of 500 mg/L for TDS. The State of Oklahoma, however, acknowledges a beneficial domestic use for general use (class II) groundwater with TDS concentrations of less than 3,000 mg/L and limited use (class III) groundwater with TDS concentrations of 3,000–5,000 mg/L (OWRB, 2015). Specific conductance ranged from 316 to 1,879 microsiemens per centimeter at 25 degrees Celsius (µS/cm at 25 °C) with a median value of 955 µS/cm at 25 °C and a mean of 1,051 µS/cm at 25 °C. The hardness ranged from 142 to 1,470 mg/L; the median hardness value was 493 mg/L, and the mean hardness value was 553 mg/L. All groundwater samples were classified as hard water (hardness greater than 180 mg/L) except for two samples that were collected from terrace deposits of northern Garvin County, which had hardness values of 142 and 170 mg/L.

In addition to the groundwater-quality data collected in 2014 from 14 wells by the OWRB (OWRB, 2017), groundwater-quality data from 3 wells sampled during 1947–53 were obtained from the NWIS database (USGS, 2020) (fig. 9). Groundwater from those wells was analyzed for selected constituents including physical properties (pH, temperature, and specific conductance), major ions, and nutrients. TDS concentrations in the groundwater samples ranged from 432 to 7,960 mg/L with a median concentration of 716 mg/L and a mean concentration of 1,471 mg/L. Specific conductance ranged from 195 to 6,990 µS/cm at 25 °C with a median value of 1,120 µS/cm at 25 °C and a mean value of 2,115 µS/cm at 25 °C. The hardness ranged from 160 to 900 mg/L; the median value was 622 mg/L, and the mean value was 219 mg/L. All groundwater samples were classified as hard water except for two samples, which had hardness values of 160 and 180 mg/L, from western Grady County.

Major-ion groundwater-quality data were converted to units of milliequivalents per liter and excluded from further analysis if the difference in the cation-anion balance for the sample was greater than 10 percent; one sample was excluded from the OWRB (2017) data, and two samples were excluded from the USGS (2020) data. The data from the remaining 13 OWRB (2017) and 12 USGS (2020) samples were plotted on a Piper (1944) diagram for visualization of groundwater types and trends (fig. 12). The northern parts of the Washita River aquifer are adjacent to Permian-age bedrock transitioning in the southern parts to the Pontotoc Group, Permian and pre-Permian Paleozoic-age red-brown to grayish shale with arkosic sandstone and limestone conglomerates near the Arbuckle Mountains (figs. 9 and 10). The groundwater samples show a downgradient transition in anions from a sulfate-rich water to a bicarbonate-rich water (fig. 12). Bicarbonate dominates the anions, particularly in the downgradient samples, resulting in the calcium-magnesium bicarbonate water type for a majority of the samples. A calcium-magnesium sulfate-chloride water type was identified in seven of the samples. These samples occurred mostly in the northwestern part of the aquifer (fig. 9A), where interactions with the Permian-age bedrock contribute the sulfate anion through dissolution of gypsum beds (fig. 10).

Aquifer Extent and Thickness

The geographic extent of the Washita River aquifer (fig. 1) was updated from the OWRB (2019b) by using information obtained from relatively new 1:100,000-scale geologic maps (Miller and Stanley, 2004; Chang and Stanley, 2010; Stanley and Chang, 2012). An older 1:250,000-scale geologic map (Heran and others, 2003) was used for a part of the aquifer near Chickasha, Okla., for which finer scale geologic maps were unavailable. Compared to the relatively coarse scale of the older 1:250,000-scale geologic map of the study area, the finer 1:100,000-scale geologic maps showed a narrower extent of the alluvium and terrace deposits that contain the aquifer, so the updated aquifer area presented in this report is smaller than previously documented.

Where present, the top of the Washita River aquifer was defined as the land-surface altitude obtained from a 10-m (horizontal resolution) digital elevation model (DEM) (USGS, 2015) with filled depressions. The altitude of the base of the Washita River aquifer, which was equivalent to the top of the bedrock, was contoured at a 50-ft interval from bedrock depths obtained from drillers’ lithologic logs, well-completion reports, and test-hole data (OWRB, 2019a; USGS, 2020) and direct-push test holes (fig. 13; USGS, 2020). For each of these data sources, the altitude of the base of the aquifer was calculated by subtracting the measured bedrock depth from the land-surface altitude. For consistency, the land-surface altitude was obtained from the 10-m DEM, even when the data source provided a land-surface altitude.

Selected wells and direct-push test holes in the Washita River aquifer were assumed to fully penetrate the aquifer (USGS, 2020; fig. 13), and their depths were assumed to be equal to the bedrock depth. Bedrock depths from these wells and test holes were given the highest priority because they had each been visited and checked by the authors of this report. Other wells and test holes included reports with drillers’ lithologic logs (OWRB, 2019a; fig. 13) that were searched for terms representing consolidated Permian-age bedrock units (such as “redbed” and “bedrock”). The altitude associated with the first occurrence of these terms in the logs was used as the altitude of the aquifer base. For logs that indicated that the well completely penetrated the Washita River aquifer (for example, the well log did not include a bedrock-related term), the altitude associated with the hole depth listed on the lithologic log was assumed to be the maximum (highest) possible altitude of the aquifer base at that location.

The mean thickness of alluvium and terrace deposits in the Washita River aquifer was about 54 ft, estimated by using the land-surface altitude and the aquifer-base altitude defined in this report. The mean thickness of alluvium and terrace deposits in reach 3 was about 56 ft and in reach 4 was about 53 ft.

Potentiometric Surface and Saturated Thickness

Potentiometric-surface maps show the altitude at which the water level would have stood in tightly cased wells at a specified time; the potentiometric surface is usually contoured or spatially interpolated from synoptic water-table-altitude measurements in many wells across an aquifer extent. Potentiometric-surface maps are used to indicate the general directions of groundwater flow in an aquifer. Groundwater generally flows perpendicular to potentiometric contours in the direction of decreasing contour altitudes (Freeze and Cherry, 1979).

In this report, the Washita River aquifer water table was assumed to be equivalent to the potentiometric surface because the aquifer is unconfined (Kent and others, 1984). The potentiometric surface for the Washita River aquifer was mapped by using groundwater-level measurements from selected wells and subtracting the measured depth to water from the land-surface altitude derived from a light detection and ranging (lidar) DEM (USGS, 2019). The selected wells were domestic and irrigation wells that were unused at the time of measurement. Methods described in Cunningham and Schalk (2011) were used to obtain groundwater-level measurements. A January 2020 potentiometric surface (fig. 14) for the Washita River aquifer was contoured primarily from synoptic groundwater-level altitudes measured during January 21–23, 2020, in 66 wells and supplemented with selected water-level altitudes measured mostly during February 2018 and December 1986–January 1987 in 24 wells (USGS, 2020; sites not listed in table 1; instead, site information is available in Rogers and others [2023]). Water-level altitudes measured in 2018 and 1986–87 were used to add control in areas where there was not a measurement in 2020. The mean saturated thickness of the Washita River aquifer was about 47 ft. The mean saturated thickness of reach 3 was about 50 ft and of reach 4 was about 45 ft (Rogers and others, 2023). Local flow in the Washita River aquifer was generally from topographically high areas toward the Washita River and major tributaries with regional flow in the aquifer generally from northwest to southeast (fig. 14).

Hydraulic and Storage Properties of Aquifer Materials

The distribution and variability of textural and hydraulic properties of aquifer materials, especially the horizontal hydraulic conductivity, were assumed to be the primary controls on groundwater flow in the Washita River aquifer. Multiple methods were used to estimate the range and central tendency of horizontal hydraulic conductivity values in the aquifer. These methods included in-place estimation in test holes, a multiwell aquifer test, and summary of data in drillers’ lithologic logs.

Horizontal Hydraulic Conductivity Estimated in Test Holes

A direct-push Geoprobe hydraulic profiling tool (HPT; Geoprobe Systems, 2015) was used at six Geoprobe test-hole sites (table 5, fig. 1) during June 24–July 15, 2021, to obtain in-place estimates of horizontal hydraulic conductivity by using the McCall (2010) method. This method can be used in the saturated zone to estimate horizontal hydraulic conductivity values between 0.1 and 150 feet per day (ft/d). During direct-push drilling, the HPT injects water into the aquifer material while logging electrical conductivity, injection pressure (corrected for hydrostatic gradient in the saturated zone), and injection flow rate at 0.05-ft depth intervals (Geoprobe Systems, 2015). The Direct Image Viewer (version 3.0) software package (Geoprobe Systems, 2020) was used to calculate discrete horizontal hydraulic conductivity values with depth by using the ratio of the injection flow rate and the hydrostatic-corrected injection pressure. The HPT test holes were drilled to a depth of refusal, which was assumed to be the bedrock contact. The 5,595 discrete (every 0.05 ft) HPT horizontal hydraulic conductivity measurements from all six test holes ranged from 0.1 to 89.0 ft/d (fig. 15) with a mean of 12.2 ft/d (table 5; USGS, 2021). These horizontal hydraulic conductivity data may be underestimated, though; small algal mats from the injection-water tank became entrained in the injection line during logging, which may have partially obstructed flow through the injection line, possibly increasing the measured injection pressure and making the aquifer appear less permeable. Progressive obstruction of the injection line may explain why the first two drilled test holes (TH11 and TH12, table 5, fig. 15A) had greater mean horizontal hydraulic conductivity measurements (22.1 and 26.7 ft/d, respectively) than did the others. Most of the HPT test holes showed generally decreasing electrical conductivity and generally increasing horizontal hydraulic conductivity with depth, indicating that the aquifer materials tend to fine upward from fine to medium sand to silt and clay (fig. 9.2 in Boggs [1995]).

Table 5.    

Hydraulic properties calculated from Geoprobe hydraulic profiling tool logs at direct-push Geoprobe test holes in the Washita River aquifer, southern Oklahoma, 2020–21.

[Data from U.S. Geological Survey (USGS; 2020, 2021). M/D/Y, month/day/year; NAVD 88, North American Vertical Datum of 1988; HPT, hydraulic profiling tool; Kh, horizontal hydraulic conductivity; --, not available or not applicable]

Table 5.    Hydraulic properties calculated from Geoprobe hydraulic profiling tool logs at direct-push Geoprobe test holes in the Washita River aquifer, southern Oklahoma, 2020–21.

Map name (fig. 1)	Field identifier	USGS station identifier	USGS site name	Latitude, in decimal degrees	Longitude, in decimal degrees	Reach	Drilling date (M/D/Y)	Land-surface altitude, in feet above NAVD 88	Well or hole depth, in feet	HPT logging date (M/D/Y)	HPT depth, in feet below land surface	HPT depth to water, in feet below land surface	HPT mean Kh, in feet per day
TH01	WYNNEWOOD3	343537097085801	02N-01E-36 DCC 1	34.5936	−97.1494	4	10/5/2020	826	40.6	7/13/2021	33.2	3.4	5.7
TH02	PAULSVALLEY1	344011097114301	02N-01E-03 CCB 1	34.6697	−97.1953	4	10/5/2020	864	48.6	--	--	--	--
TH03	PAULSVALLEY2	344053097101001	03N-01E-35 CDD 1	34.6814	−97.1696	4	10/2/2020	851	46.6	--	--	--	--
TH04	LINDSAY3	344842097283401	04N-03W-23 AAA 1	34.8118	−97.4762	4	9/11/2020	958	46.9	--	--	--	--
TH05	LINDSAY1	344916097311401	04N-03W-16 BDA 1	34.8212	−97.5205	4	10/2/2020	961	63.2	7/14/2021	63.0	15.6	6.8
TH06	LINDSAY4	345122097294001	05N-03W-34 DDD 1	34.8562	−97.4946	4	10/2/2020	956	74.1	--	--	--	--
TH07	BRADLEY1	345402097434901	05N-05W-16 DCB 1	34.9005	−97.7302	4	10/6/2020	1,019	49.2	7/14/2021	48.5	4.3	7.1
TH08	BRADLEY3	345422097421801	05N-05W-14 CBB 1	34.9061	−97.7051	4	10/9/2020	1,013	90.2	--	--	--	--
TH09	ALEX3	345749097484101	06N-06W-26 CBB 1	34.9637	−97.8114	3	10/1/2020	1,038	59.2	--	--	--	--
TH10	ALEX2	345817097470602	06N-06W-24 DCC 2	34.9715	−97.7849	3	10/1/2020	1,038	43.4	--	--	--	--
TH11	ALEX1	345829097484101	06N-06W-23 CCB 1	34.9748	−97.8115	3	10/9/2020	1,043	57.5	6/28/2021	58.9	19.2	22.1
TH12	CHICKASHA2	350513098021301	07N-08W-10 CCD 1	35.0870	−98.0369	3	10/1/2020	1,115	63.6	6/24/2021	63.6	2.8	26.7
TH13	ANADARKO2	350513098115501	07N-10W-12 DDD 1	35.0869	−98.1987	3	10/8/2020	1,164	69.2	7/15/2021	68.5	7.4	4.8
TH14	CHICKASHA1	350606098022501	07N-08W-03 CCC 1	35.1017	−98.0403	3	10/9/2020	1,123	49.3	--	--	--	--
TH15	ANADARKO1	350610098114701	07N-09W-06 CCC 1	35.1028	−98.1965	3	10/8/2020	1,165	98.0	--	--	--	--
TH16	CHICKASHA3	350700098022501	07N-08W-04 AAA 1	35.1167	−98.0403	3	10/1/2020	1,114	35.2	--	--	--	--
TH17	ANADARKO3	350701098121901	07N-10W-01 ABB 1	35.1170	−98.2054	3	10/8/2020	1,172	92.7	--	--	--	--
TH18	CHICKASHA4	350753098022001	08N-08W-28 DDD 1	35.1314	−98.0389	3	10/9/2020	1,117	47.3	--	--	--	--
Mean											55.9	8.8	12.2

Hydrologic Boundaries

Hydrologic boundaries in the conceptual model represent actual sources (inflows) and sinks (outflows) of water to and from the aquifer. Boundaries that act as both inflows and outflows may be referred to as “net inflows” or “net outflows” depending on which flow component dominates.

Recharge

Recharge is the predominant inflow to the Washita River aquifer. Recharge is defined in this report as the amount of precipitation that infiltrates from the land surface through the unsaturated zone and reaches the water table over a given time. This definition of recharge includes irrigation return flows to groundwater. Other processes such as stream seepage or lateral groundwater flow from adjacent hydrogeologic units are not considered recharge and are accounted for separately in the conceptual-model water budget. Recharge rates are controlled by many factors including precipitation rate, land-surface gradient, soil and sediment permeability, evapotranspiration rates, and vegetation cover type. Although recharge rates are difficult to measure because of high spatial and temporal variability, methods that include environmental tracers, physical measurements, streamflow-hydrograph techniques, and computer codes can be used to estimate recharge rates. For this report, a groundwater-hydrograph-based water-table fluctuation (WTF) method (Healy and Cook, 2002) was used to estimate localized recharge rates for 2017–19, and a code-based water-balance-estimation technique (Westenbroek and others, 2010) was used to estimate spatially distributed recharge rates for the 1980–2017 study period.

Water-Table Fluctuation Method

The WTF method was the primary method used to estimate recharge to the Washita River aquifer. The WTF method assumes that rises in groundwater levels in unconfined aquifers that occurred during a relatively short period (hours to a few days) can be attributed to recharge arriving at the saturation zone following a period of precipitation. This method is most appropriately applied to groundwater wells with shallow water tables and hydrographs that display sharp rises in groundwater levels after precipitation (Healy and Cook, 2002). The WTF method cannot account for a steady rate of recharge or recharge from sources other than precipitation. Annual recharge (R), in inches per year, was estimated by using the following equation:

R

=

Sy

× Σ(

Δh

/

Δt

)

(2)

where

Sy

is the specific yield (dimensionless);

Δh

is the rise in groundwater-level altitude, in inches; and

Δt

is the change in time, in years.

Water-level hydrographs from three USGS continuous water-level recorder wells (W01, W02, and W05, table 1; USGS, 2020) were used to estimate annual recharge during 2017–19. W01 was located near Lindsay, W02 was located near Chickasha, and W05 was located near Pauls Valley (fig. 1). The water-level hydrographs from other continuous water-level recorder wells in the study area were not analyzed because they were affected by nearby pumping (in the case of W03, fig. 17A) or because they did not display sharp water-level rises (in the case of W04 and W06, fig. 17B). Annual precipitation data for the 2016–19 period of analysis were obtained from the Chickasha (CHIC) and Pauls Valley (PAUL) climate stations (fig. 1; Oklahoma Mesonet, 2019). Using the daily precipitation record from the nearest climate station and a specific yield of 0.10 from Hemann (1985), the annual recharge estimates for wells W02 and W05 in 2017 were 5.3 inches (in.; about 13.5 percent of the station annual precipitation) and 7.4 in. (about 18.5 percent of the station annual precipitation), respectively (table 7). The annual recharge estimates for wells W02 and W05 in 2018 were 9.0 in. (about 23.7 percent of the station annual precipitation) and 11.5 in. (about 22.8 percent of the station annual precipitation), respectively. In a year-long period from April 1, 2018, through March 31, 2019, the annual recharge estimate at well W01 was 10.3 in. (about 22.6 percent of the station annual precipitation). When those recharge estimates are normalized by the annual mean precipitation of 40.1 in/yr for the 1980–2017 study period and averaged, the resulting station-averaged mean annual recharge for 2017–19 in reach 3 and 4 of the Washita River aquifer was 8.1 in/yr (about 20.2 percent of mean annual precipitation for the 1980–2017 period). Multiplied by the 178,938-acre modeled area and unit converted, the WTF-calculated mean annual recharge for 2017–19 was 120,783 acre-ft/yr; this value was used as the conceptual-model recharge during 1980–2017 (table 6).

Table 7.    

Summary of recharge amounts estimated using the water-table fluctuation method for the Washita River aquifer, southern Oklahoma, 2017–19.

[Continuous water-level recorder wells W03 and W04 were not suitable for analysis with the water-table fluctuation method (Healy and Cook, 2002). Dates for the periods are in month-day-year (MM-DD-YYYY) format; --, data not available or not applicable]

Table 7.    Summary of recharge amounts estimated using the water-table fluctuation method for the Washita River aquifer, southern Oklahoma, 2017–19.

	U.S. Geological Survey continuous water-level recorder well (fig. 1, table 1)
	W02	W01	W05
Annual mean precipitation 1980–2017, in inches per year, southern Oklahoma, Climate Division 8 (National Centers for Environmental Information, 2021)	40.1	40.1	40.1
Climate station (fig. 2, table 1)	CHIC	PAUL	PAUL
Estimated specific yield	0.10	0.10	0.10
Period 1: 01-01-2017 to 12-31-2017
Station annual precipitation, in inches	39.4	--	40.1
Sum of water-level rises, in feet	4.4	--	6.2
Recharge, in inches per year	5.3	--	7.4
Recharge, percent of annual precipitation	13.5	--	18.5
Recharge, in inches per year, normalized to mean annual precipitation during 1980–2017	5.4	--	7.4
Period 2: 01-01-2018 to 12-31-2018
Station annual precipitation, in inches	37.9	--	50.5
Sum of water-level rises, in feet	7.5	--	9.6
Recharge, in inches per year	9.0	--	11.5
Recharge, percent of annual precipitation	23.7	--	22.8
Recharge, in inches per year, normalized to mean annual precipitation during 1980–2017	9.5	--	9.2
Period 3: 04-01-2018 to 03-31-2019
Station annual precipitation, in inches	--	45.5	--
Sum of water-level rises, in feet	--	8.6	--
Recharge, in inches per year	--	10.3	--
Recharge, percent of annual precipitation	--	22.6	--
Recharge, in inches per year, normalized to mean annual precipitation during 1980–2017	--	9.0	--
Mean annual recharge during 2017–19, in inches per year, normalized to mean annual precipitation 1980–2017			8.1

Soil-Water-Balance Code

The Soil-Water-Balance code (SWB; Westenbroek and others, 2010) was used to estimate the amount and spatial distribution of daily groundwater recharge to the Washita River aquifer for each month of the 1980–2017 study period. SWB uses a modified Thornthwaite and Mather (1957) soil-water-balance method on a gridded data structure to compute the daily amount of infiltration (minus losses) that exceeds the storage capacity of the plant root zone. The soil-water-balance equation (modified from Westenbroek and others [2010]) has the form of the following equation:

R

= (

P

+

S

+

R_i

) – (

Int

+

R_o

+

P_et

) –

ΔSm

(3)

where

R

is recharge, in inches per day;

P

is precipitation, in inches per day;

S

is snowmelt, in inches per day;

R_i

is surface runoff inflow, in inches per day;

Int

is plant interception, in inches per day;

R_o

is surface runoff outflow, in inches per day;

P_et

is potential evapotranspiration, in inches per day; and

ΔSm

is the change in soil moisture, in inches per day.

SWB requires climate and landscape characteristic data inputs including precipitation, temperature, soil-water storage capacity, hydrologic soil group, surface-water flow direction, and land-cover type (Westenbroek and others, 2010). Each of these inputs was assigned to a user-specified grid that was 600 by 500 cells of 800 by 800 ft each. The landscape inputs were assumed to remain constant during the study period, but climate data inputs varied daily. Daily climate data including precipitation and minimum and maximum air temperature were obtained for 14 selected climate stations (fig. 2, table 1; National Centers for Environmental Information, 2022) in and near the study area and gridded using inverse-distance-weighted (Esri, 2021b) interpolation. Accumulated snowmelt was added on the basis of the daily minimum and maximum air temperatures. Potential evapotranspiration was calculated by using the Hargreaves and Samani (1985) method for a reference latitude range of 34.5–35.4 degrees. Surface runoff was routed downslope by using a flow-direction grid derived from a 10-m DEM (USGS, 2015). Depressions in the DEM were filled by using the ArcGIS Fill tool (Esri, 2021a) to ensure correct routing of surface runoff and to eliminate areas of internal drainage that can result in unrealistically high amounts of recharge. Soil properties (soil-water storage capacity and hydrologic soil group) were derived from the Soil Survey Geographic database (NRCS, 2022), which is an inventory of generalized soil characteristics. Land-cover types (Fry and others, 2011; Multi-Resolution Land Characteristics Consortium, 2011; fig. 3) were used in conjunction with hydrologic soil groups to partition daily precipitation into plant interception (Int) and surface runoff (R_i and R_o) components and assign plant root-zone depths. The root-zone depths for grass/pasture and cropland (the dominant land-cover types for land overlying the aquifer; fig. 3) varied with soil texture but ranged from about 0.8 to 1.5 ft (40 percent of the values used by Westenbroek and others [2010] for permeable glacial deposits in Wisconsin). The soil-water storage capacity, analogous to specific yield in the saturated zone, multiplied by the root-zone depth determined the maximum volume of water available in the root zone. Changes in soil moisture (ΔSm) exceeding the soil-water storage capacity were assumed to be recharge (R) to the saturated zone. Larger root-zone depths resulted in increased evapotranspiration of water from the root zone and decreased recharge; smaller root-zone depths resulted in decreased evapotranspiration of water from the root zone and increased recharge. Recharge from irrigation was not simulated by SWB but was assumed to be negligible given the relatively small amount of irrigation groundwater use in the study area. The SWB model for the Washita River aquifer study area is included in Rogers and others (2023).

The SWB-estimated mean annual recharge for the 1980–2017 study period was about 4.9 in/yr, or about 12.2 percent of the mean annual precipitation of 40.1 in/yr during 1980–2017 (fig. 18A). The minimum and maximum SWB-estimated annual recharge values were about 1.1 in/yr in 2003 and 11.6 in/yr in 2015, respectively. Recharge efficiency (mean monthly recharge as a percentage of mean monthly precipitation) was greatest from November to February, when evapotranspiration was at a minimum; recharge in these months was greater than 20 percent of precipitation (fig. 18B). Recharge efficiency was lowest in July and August, when evapotranspiration was at a maximum; recharge in these months was less than 6 percent of precipitation. Spatially, mean annual recharge for the study period was greatest in areas of alluvium near the Washita River and tributaries (fig. 19). The SWB-estimated mean annual recharge generally increased from northwest to southeast in the Washita River aquifer; the mean annual recharge for reach 4 was about 14 percent greater than mean annual recharge for reach 3 of the Washita River aquifer. Increased recharge from reach 3 to reach 4 is correlated with increased precipitation from northwest to southeast in the study area (Oklahoma Climatological Survey, 2021a, b).

Mean annual recharge rates estimated from SWB were compared to published estimates of mean annual recharge for the Washita River aquifer. Kent and others (1984) estimated a mean annual recharge of 3.72 in/yr, or about 12 percent of the mean annual precipitation of 32.2 in/yr during a 20-year study period (July 1, 1973–July 1, 1993). That estimate was applied across the entire Washita River aquifer in Oklahoma and included reaches that are outside the scope of this report. Schipper (1983) estimated a mean annual recharge of 3.17 in/yr, or about 13 percent of the mean annual precipitation (24.90 in/yr), during the same 20-year study period for a groundwater model of reach 1 of the Washita River aquifer (in western Oklahoma). Ellis and others (2020) estimated a mean annual recharge for reach 1 of the Washita River aquifer (in western Oklahoma) of 3.15 in/yr, or about 14 percent of the mean annual precipitation of 23 in/yr, during a 36-year study period (1980–2015). The mean annual recharge rate of 4.9 in/yr that was estimated by using SWB for the study area in this report is thus within the range of published estimates of recharge for the Washita River aquifer, particularly when calculated as a percentage of the applicable mean annual precipitation.

The SWB-estimated mean annual recharge was, however, about 40 percent less than the WTF-estimated mean annual recharge used for the conceptual model. SWB is a numerical model, and SWB-estimated recharge must, therefore, be checked against and calibrated to field-based measurements or estimates of recharge like those estimated using WTF. In other reports (Smith and others, 2017, 2021; Ellis, 2018a; Ellis and others, 2020), SWB-estimated recharge was calibrated by decreasing root-zone depths until the SWB-estimated mean annual recharge approximated the conceptual-model recharge. That method was also applied in this report to a point but was halted when the mean root-zone depths began to decrease to less than 1.2 ft. Additional adjustment of SWB-estimated recharge was accomplished during calibration of the numerical model as detailed in the “Numerical Groundwater-Flow Model” section of this report.

Streambed Seepage

Base flow can be measured directly in streams when the runoff component of streamflow is at or near zero cubic feet per second (Garner and Bills, 2012). When base-flow measurements are collected at multiple locations over a short span of time, they are called synoptic base-flow (seepage-run) measurements. These measurements can be used to calculate net streambed seepage and classify stream reaches as gaining (exhibiting an increase in base flow between the upstream and downstream endpoints of the reach) or losing (exhibiting a decrease in base flow between the upstream and downstream endpoints of the reach). Streambed seepage rates were calculated as the difference in measured base flows (adjusted for tributary inflows) divided by the stream reach length between measurement locations.

Synoptic base-flow measurements were collected by using the methods of Rantz and others (1982) during a period of minimal runoff on January 23, 2018. These measurements, along with data from active USGS streamgages, served to delineate tributary inflow and base-flow conditions across the aquifer in January 2018. Tributaries along the Washita River alluvium contributed base flows ranging from 0.21 to 32.3 cubic feet per second (ft³/s; fig. 20). Tributaries with unmeasured base flows, many of which were visually observed as having minimal or no flows, were assumed to contribute no base flows. The base-flow measurements on the Washita River main stem indicated that there were gaining and losing reaches of the stream. The greatest base-flow gain (12.6 cubic feet per second per mile) occurred over a 5-mi reach south of Davis, and the greatest base-flow loss (–4.8 cubic feet per second per mile) occurred over a 9-mi reach near Lindsay.

Net streambed seepage terms for the conceptual model were assumed to be net outflows from the aquifer and were estimated from mean annual base flows computed by using the Base-Flow Index code (Wahl and Wahl, 1995) in the USGS Groundwater Toolbox (Barlow and others, 2015) on selected streamgages in the study area. Net streambed seepage in reach 3 for the 1980–2017 study period (49,041 acre-ft/yr, table 6) was roughly estimated as the mean annual base flow at the Alex gage minus the mean annual base flows at the Anadarko gage and selected tributary streamgages: USGS streamgage 07327000 Sugar Creek near Gracemont, Okla.; USGS streamgage 07327550 Little Washita River east of Ninnekah, Okla.; and USGS streamgage 07328070 Winter Creek near Alex, Okla. (fig. 2, tables 1, 3; tributary base flows are summarized in Rogers and others [2023]). Net streambed seepage in reach 4 during 1980–2017 (75,130 acre-ft/yr, table 6) was estimated as the mean annual base flow at the downgradient end of the aquifer (estimated using data from the Dickson gage) minus the mean annual base flows at the Alex gage and tributary streamgages: USGS streamgage 07328180 North Criner Creek near Criner, Okla.; USGS streamgage 07329000 Rush Creek at Purdy, Okla.; USGS streamgage 07329700 Wildhorse Creek near Hoover, Okla.; and USGS streamgage 07329780 Honey Creek below Turner Falls near Davis, Okla. (fig. 2, tables 1, 3; tributary base flows are summarized in Rogers and others [2023]). When available, streamflow data from the 1980–2017 study period were used to estimate mean annual base flows for the conceptual model. Streamflow data from outside the study period were used to estimate mean annual base flows when other data were unavailable. The total net streambed seepage was 124,171 acre-ft/yr and accounted for 87.8 percent of total outflows in the conceptual-model water budget (table 6).

Conceptual-Model Water Budget

The conceptual-model water budget (table 6) summarized mean water flows exchanged between each hydrologic boundary and the Washita River aquifer for the 1980–2017 study period. Components of the water budget were estimated from analyses of available data; however, where data were not available, assumptions were made by using published analogs as described in the “Hydrologic Boundaries” section. Recharge accounts for 85.5 percent of the conceptual-model inflows to the Washita River aquifer, and net streambed seepage accounts for 87.8 percent of the outflows from the Washita River aquifer. Net lateral groundwater flow (14.5 percent of inflows) and saturated-zone evapotranspiration (11.3 percent of outflows) were the only other components estimated to be greater than 5 percent of inflows or outflows. Well withdrawals accounted for less than 1 percent of conceptual-model water-budget outflows.

Spatial and Temporal Discretization

The model domain (fig. 21) of the Washita River aquifer was spatially discretized into 500 rows; 600 columns; 12,179 active cells measuring 800 by 800 ft each; and a single convertible layer (Harbaugh, 2005) based on the hydrogeologic framework described in this report. The cell size was chosen to minimize model-processing time while still representing the variability of properties being simulated. The chosen cell size also ensured that the narrowest part of the aquifer was represented by no fewer than three cells; model instability and water-level volatility often occur in narrow parts of the modeled area where groundwater flow is focused into fewer cells. The single convertible layer represented the undifferentiated Quaternary-age alluvium and terrace deposits with variable thickness determined from the hydrogeologic framework; the underlying bedrock was not represented as a layer. The altitude of the top of the aquifer was multiplied by 1.01 in the numerical model to prevent confined aquifer conditions that occur as a side effect of model discretization when the simulated water-table altitude exceeds the altitude of the top of the aquifer. This multiplier, which adds about 10 ft to the altitude of the top of the aquifer, has no effect on most parts of the model where the simulated water-table altitude is below the altitude of the top of the aquifer.

The active modeled area (fig. 21) was initially derived from the Washita River aquifer extent (modified from OWRB [2019b]) as defined in the hydrogeologic framework of the Washita River aquifer. The active modeled area was expanded or contracted in some areas to ensure that each active cell was in connection with at least one other active cell. The active modeled area was further modified by inactivating some model cells in narrow areas along selected Washita River tributaries to improve model stability.

The numerical model was temporally discretized into 456 monthly transient stress periods (each with two time steps to improve model stability) representing the 1980–2017 study period (the numerical-modeling period [hereinafter referred to as the “modeling period”] coincides with the 1980–2017 study period). An initial steady-state stress period represented mean annual inflows to, and outflows from, the aquifer. The steady-state solution was used as the initial condition for subsequent transient stress periods, as well as some groundwater-availability scenarios. The numerical model was constructed by using length and time units of feet and days, respectively.

Simulation of Hydrologic Boundaries and Hydraulic Parameters

Hydrologic boundaries in the numerical model (fig. 21) define where and how water may enter or leave the model and include specified-flux and head-dependent boundaries (Harbaugh and others, 2000). Specified-flux boundaries were used to simulate recharge and well withdrawals. Head-dependent boundaries were used to simulate streambed seepage, saturated-zone evapotranspiration, and lateral groundwater flow exchanged with adjacent (upstream and downstream) alluvium and terrace deposits of the Washita River aquifer and with the surrounding bedrock units. When available, hydrologic data (along with data-based assumptions and analogs) were used to estimate or constrain precalibration model inputs for each hydrologic boundary.

Well Withdrawals

Well withdrawals for the 1980–2017 study period were simulated by using the Well package (Harbaugh and others, 2000). Annual reported groundwater use for each permit was evenly distributed among all well locations recorded for that permit (OWRB, 2019d). In the case of permits for which there were no recorded wells, a single well was placed in the center of the land parcel recorded for the permit. Annual reported groundwater use was temporally split into model stress periods by using the mean monthly water demand distribution (fig. 22; OWRB, 2012a) from Oklahoma Comprehensive Water Plan water-management planning basins 14–16 (fig. 2; OWRB, 2019c). The monthly water demand for irrigation was greatest in the summer months; July–September accounted for about 76 percent of irrigation groundwater use (fig. 22). The monthly water demand for public supply also was greatest in the summer months; however, July–September only accounted for about 34 percent of public-supply groundwater use (fig. 22).

Storage and Hydraulic Properties

The Upstream Weighting package of MODFLOW-2005 (Niswonger and others, 2011) was used to represent storage and hydraulic properties of the Washita River aquifer. Initial storage and hydraulic property values were similar to those used in numerical models for other alluvial aquifers in Oklahoma (Ellis and others, 2017, 2020; Smith and others, 2017) and identical to those used in the numerical model for the Salt Fork Red River aquifer (Smith and others, 2021). The storage and hydraulic property values were adjusted during model calibration but were held spatially uniform and temporally constant through all stress periods in all simulations. The specific yield (sy3 and sy4, table 8) and specific storage (ss, table 8) of the Washita River aquifer were initially set as 0.10 (dimensionless) and 1×10^–5 inverse foot, respectively. Horizontal hydraulic conductivity (hkpp001–hkpp146, table 8) was represented by using 146 pilot points (1–146, fig. 23; Doherty, 2010) placed at 12,800-ft intervals (1 at every 16th column and row) within a 10,560-ft buffer of the aquifer boundary. (Pilot points are regularly spaced control points used to support an interpolation of cell-based horizontal hydraulic conductivity values across an active modeled area.) The pilot points were assigned an initial horizontal hydraulic conductivity value of 20 ft/d and interpolated before each model run to create the horizontal hydraulic conductivity array read by the numerical model (Doherty, 2010). During calibration, pilot-point horizontal hydraulic conductivity values were allowed to vary independently between 10 and 100 ft/d (Rogers and others, 2023).

Solver Settings and Budget Percentage Discrepancies

Most of the settings for the Newton solver were unchanged from suggested input values (Winston, 2018), which were converted to model units of feet and days. The head tolerance was increased to 1.5×10^–3 ft, and the flux tolerance was increased to 5.0×10⁴ cubic feet per day. These settings improved model stability while keeping the budget percentage discrepancies in the solution under 1.0 percent for all 457 stress periods; the largest budget percentage discrepancy (by absolute value) was about 0.16 percent in stress period 411 (Rogers and others, 2023).

Calibration

Model calibration is the process of systematically changing selected model input values (parameters) within limits to improve the fit between model-simulated data and observed data (calibration targets). The preferred calibration results (1) minimize the differences (residuals) between simulated and observed data and (2) conform to the predetermined conceptual model. The calibration process for the numerical model included both manual and automated adjustment of parameters. The manual calibration approach was primarily focused on aligning the numerical-model water budget, especially the recharge and streambed seepage components, to the conceptual-model water budget. The manual calibration also involved selected structural elements of the model, such as boundary elevations, that are not easily adjusted by automated calibration methods. The automated calibration approach was focused solely on minimizing residuals and used the PEST++ iterative ensemble smoother (White, 2018) to reduce run times associated with the calibration of highly parameterized models.

Calibration Targets

The suite of calibrated parameter values was evaluated on the basis of the minimization of an objective function that quantifies the difference between simulated and observed values included in the model (White, 2018). The objective function was calculated as the sum of squared weighted residuals for calibration targets in five observation groups: water-table-altitude observations, base-flow observations at the Alex gage, base-flow observations at the Pauls Valley gage, base-flow observations (estimated) at the downgradient end of the model, and the conceptual-model recharge rate (table 9). Observation weights were required to prevent calibration targets of any one observation group from dominating the objective function and, therefore, becoming the sole focus of parameter adjustments during automated calibration (Hill and Tiedeman, 2007; Doherty, 2010). The observation weights (1.0, 5.0 × 10^–7, 1.7 × 10^–7, 1.8 × 10^–7, 2.5 × 10^–2 for water-table-altitude observations, base-flow observations at the Alex gage, base-flow observations at the Pauls Valley gage, base-flow observations [estimated] at the downgradient end of the model, and conceptual-model recharge rate, respectively) were uniform for all observations from the transient simulation within each observation group; the base-flow observations from the steady-state simulation were given increased weights (200 times greater) relative to the base-flow observations from the transient simulation. The observation weights were chosen so that, at the beginning of automated calibration (that is, precalibration), the observation groups accounted for 24.0, 10.5, 55.3, 10.2, and 0.0 percent of the objective function, respectively (table 9). This weighting scheme reflects an assumption that base-flow observations (totaling 76.0 percent of the objective function) are more important than water-table-altitude observations for the purpose of determining the parameter values that are most influential in estimating groundwater availability (the primary subject of this report) in an unconfined alluvial aquifer. Base flows are generally more sensitive to changes in specific yield than they are to changes in water-table altitudes. Water-table altitudes are generally more sensitive to changes in horizontal hydraulic conductivity than to changes in specific yield.

Table 9.    

Components of the objective function for the automated calibration of the numerical groundwater-flow model of the Washita River aquifer, southern Oklahoma, 1980–2017.

[Objective function is calculated as the sum of squared weighted residuals. SS, steady-state simulation; TR, transient simulation; Min, minimum; Max, maximum; RMSE, root-mean-square error; OWRB, Oklahoma Water Resources Board; USGS, U.S. Geological Survey; E, base-10 exponent (for example, 1.7E-07 means 1.7 × 10⁻⁷)]

Table 9.    Components of the objective function for the automated calibration of the numerical groundwater-flow model of the Washita River aquifer, southern Oklahoma, 1980–2017.

Observation group	Source1	Number of observations	Weight multiplier	Precalibration residuals						Calibrated residuals
		SS TR	SS TR	Min	Mean	Max	Interquartile range	RMSE	Objective function component (percent in parentheses)	Min	Mean	Max	Interquartile range	RMSE	Objective function component (percent in parentheses)
Water-table altitude, in feet (headobs)	Water-table-altitude observations (OWRB, 2019a; USGS, 2020)	63 130	2.0E+00 2.0E+00	−22.5	−6.5	17.2	9.1	18.6	66,641 (24.0)	−19.5	−4.1	10.5	6.8	13.8	36,522 (24.9)
Base flow, in cubic feet per second (gageobs1)	Alex gage (USGS, 2020)	1 432	1.0E-04 5.0E-07	−1,238.8	14.0	970.2	95.6	8.2	28,993 (10.5)	−1,228.1	−7.7	981.8	93.9	7.9	27,109 (18.5)
Base flow, in cubic feet per second (gageobs2)	Pauls Valley gage (USGS, 2020)	1 456	3.4E-05 1.7E-07	−1,207.7	105.4	2,549.3	213.2	18.3	153,413 (55.3)	−1,224.2	45.6	2,489.8	213.8	10.8	53,763 (36.6)
Base flow, in cubic feet per second (gageobs3)	Downgradient end of model; adjusted from Dickson gage (USGS, 2020)	1 456	3.6E-05 1.8E-07	−1,457.2	8.8	4,952.7	234.1	7.9	28,231 (10.2)	−1,433.3	−8.8	4,932.0	234.1	8.0	29,173 (19.9)
Conceptual-model recharge, in acre-feet per year (budget)	Conceptual model (table 6)	1 1	2.0E-02 2.0E-02	353.7	353.8	353.8	0.0	7.1	100.1 (0.0)	−829.8	−416.5	−3.2	413.3	11.7	275.4 (0.2)
Total		1,542		Precalibration objective function:					277,378	Calibrated objective function:					146,842

¹

Alex gage refers to USGS streamgage 07328100 Washita River at Alex, Okla.; Pauls Valley gage refers to USGS streamgage 07328500 Washita River near Pauls Valley, Okla.; and Dickson gage refers to USGS streamgage 07331000 Washita River near Dickson, Okla. (fig. 1, table 1).

Water-Table-Altitude Observations

The Head Observation (HOB) package (Hill and others, 2000) was used to compare simulated water-table altitudes to observed water-table altitudes in the Washita River aquifer for the 1980–2017 modeling period. Observed water-table altitudes were calculated from depth-to-water measurements and the land-surface altitude obtained from a lidar DEM (USGS, 2019). For consistency, the land-surface altitude was obtained from the lidar DEM, even when the data source provided a land-surface altitude. Sites with water-table altitude observations were screened for location errors, and observations were filtered such that only one observation per cell per stress period was included in the HOB package. Only 130 water-table-altitude observations (OWRB, 2019a; USGS, 2020) from 63 wells were included for the transient simulation in the HOB package (fig. 24A). Water-table-altitude observations from the period 2018–20 were assigned to 2017 in the HOB package because few historical observations were available from the modeling period. Some water-table-altitude observations from the period 2018–20 were assigned to the same month and day in 2017, but a set of 51 synoptic water-table-altitude observations from January 2020 were assigned to the same day in December 2017. Water-table-altitude observations from the transient simulation were averaged by well to derive 63 calibration targets for the steady-state simulation because no water-table-altitude observations were available in 1980, at the beginning of the modeling period.

Base-Flow Observations

Base-flow observations were available for the full modeling period (456 months) at the Pauls Valley gage on the Washita River and for most of the modeling period (432 months) at the Alex gage (fig. 24B) (USGS, 2020) on the Washita River within the modeled area (fig. 1). Base-flow observations also were available for the full modeling period (456 months) at the Dickson gage (fig. 24B) (USGS, 2020), which is downstream from the modeled area on the Washita River (fig. 2). Base-flow observations from the Dickson gage were decreased by 30 percent (to account for drainage that originated from the tributaries outside of the modeled area) and used as estimated base-flow observations for the Washita River at the downgradient end of the model. The monthly mean base flows for these three streamgages, accounting for 1,344 observations (table 9), were used as calibration targets for the transient simulation (Rogers and others, 2023). The mean annual base flows for the three streamgages (table 3), accounting for three observations (table 9), were used as calibration targets for the steady-state simulation. Base-flow observations from other selected streamgages were not used as calibration targets because they (1) were previously used to define inflows for SFR2 or (2) had periods of record that were too short to be representative of streamflow conditions during the modeling period.

Calibration Results

Calibration results were evaluated on the basis of the reduction of the objective function and the general fit of the calibrated numerical-model water budget to the predetermined conceptual-model water budget. Automated calibration reduced the calibrated objective function total by about 47 percent from the precalibration objective function (table 9). Most of that reduction resulted from reducing base-flow residuals at the Pauls Valley gage. Base-flow residuals for the Alex gage and the downgradient end of the model were mostly reduced during a round of manual calibration prior to calculation of the precalibration objective function.

Observation-Sensitivity Analysis

As part of the calibration process, an observation-sensitivity analysis was performed with the iterative parameter estimation software package PEST (Doherty, 2010) to identify which parameters were most effective (and most ineffective) at reducing the objective function. PEST is used to measure the changes in calibration-target residuals resulting from 1-percent changes in each parameter and records those changes in a Jacobian matrix with dimensions equal to the number of observations times the number of parameters (Doherty, 2010). Parameters with the greatest effects on calibration-target residuals have the greatest observation sensitivities (Rogers and others, 2023).

Parameters that were closely linked, either spatially or temporally, to observations typically had the greatest observation sensitivities in the numerical model. Parameters with the greatest observation sensitivities were the saturated-zone evapotranspiration reference altitude multiplier (evttop), the steady-state recharge-rate multiplier (rch001), the altitude of the upstream GHB (ghdup), the steady-state saturated-zone evapotranspiration-rate multiplier (evtavg), and the streambed depth for SFR2 zone 3 (sbd03) (Rogers and others, 2023; table 8). Horizontal hydraulic conductivity pilot points 15, 26, 71, 77, 78, 79, and 89 (hkpp015, hkpp026, hkpp071, hkpp077, hkpp078, hkpp079, and hkpp089) were the only horizontal hydraulic conductivity pilot points with sensitivities ranking in the top 100 (of 753 calibration parameters). The horizontal hydraulic conductivity pilot points with the highest ranking sensitivities were in parts of the aquifer with large spatial and temporal concentrations of water-table-altitude observations. The specific yield for reach 3 (sy3) and the specific yield for reach 4 (sy4) ranked 72 and 28 (of 753 calibration parameters), respectively. The steady-state recharge-rate multiplier parameter (rch001) determined how much recharge was applied in the steady-state simulation, which included about one-third of the water-table-altitude observations (Rogers and others, 2023). The specific yield parameters (sy3 and sy4) controlled the volume of water released (from storage) as streambed seepage to streams across the entire aquifer, which in turn partially controlled how much base flow was simulated in the Washita River at the Alex gage, Pauls Valley gage, and downgradient end of the model.

To graphically summarize the observation-sensitivity analysis, observation-group sensitivities were calculated as the weighted sum of the absolute Jacobian matrix output for each parameter group (fig. 25). Because recharge parameters and horizontal hydraulic conductivity parameters influence every cell in the model, those parameter groups had the greatest observation-group sensitivities. The storage parameter group, which included specific yield and specific storage parameters, had the least observation-group sensitivities.

Calibrated Parameter Values

The calibrated parameter values (Rogers and others, 2023) selected for the numerical model were the combined result of manual and automated calibration approaches. Many of the calibrated parameter values were at the minimum or maximum bounds specified in the PEST control file (Rogers and others, 2023). Calibrated parameter values at bounds indicate that the objective function could likely be further reduced by expanding those bounds; doing so, however, would divert the numerical model further from the conceptual model or increase model instability. For the recharge parameter group (rch, table 8), 333 of 457 calibrated recharge-rate multiplier values were assigned to either the minimum or maximum bounds; 156 of those values were at the minimum bound (1.40), and 177 were at the maximum bound (2.00) (Rogers and others, 2023). For the horizontal hydraulic conductivity parameter group (hyd, table 8), 53 of 146 calibrated parameter (pilot point) values were assigned to either the minimum or maximum bounds; 15 of those values were at the minimum bound (10.0 ft/d), and 38 were at the maximum bound (100.0 ft/d) (Rogers and others, 2023). However, the minimum value interpolated in the calibrated horizontal hydraulic conductivity array was increased to 20.0 ft/d (fig. 23) to alleviate problems with flooding (cells with water levels above land surface) in low-lying areas; 32 of the calibrated parameter (pilot point) values were less than or equal to 20.0 ft/d. Most of the calibrated parameter values in the evapotranspiration group (evt, table 8) were at bounds, including the saturated-zone evapotranspiration extinction (root-zone) depth, the steady-state saturated-zone evapotranspiration-rate multiplier, and 7 of the 12 transient (monthly) saturated-zone evapotranspiration-rate multipliers. For the streamflow routing parameter group (sfr, table 8), all but 2 of the streambed hydraulic conductivity zone calibrated parameter values (sbk01–sbk28), all but 4 of the streambed thickness zone calibrated parameter values (sbb01–sbb28), all but 8 of the streambed depth zone calibrated parameter values (sbd01–sbd28), and all but 12 of the streambed width zone calibrated parameter values (sbw01–sbw28) were at bounds.

Comparison of Simulated and Observed Calibration Targets

Conforming to MODFLOW convention (Harbaugh, 2005), calibration-target residuals in this report were calculated as observed minus simulated values; positive residuals indicate lower simulated than observed values, and negative residuals indicate higher simulated than observed values. The mean calibrated residual for water-table-altitude observations was –4.1 ft (table 9), indicating that, on average, simulated water levels were slightly higher than observed water levels. The combined (steady-state and transient simulation) water-table-altitude root-mean-square error (RMSE) was 13.8 ft, and the interquartile range was 6.8 ft (table 9, fig. 26B).

The water-table-altitude residuals with the greatest magnitudes were associated with wells near Alex, in the downgradient part of reach 3 and the upgradient part of reach 4 of the Washita River aquifer (figs. 26A and 27). The residuals near Alex were mostly negative, indicating that simulated water-table altitudes were higher than observed water-table altitudes (fig. 27). In that part of the model, given the weights applied to observations in the objective function (table 9), there was a tradeoff between minimizing the water-table-altitude residuals and minimizing the base-flow residuals at the Alex streamgage.

Annual water-table-altitude observations representing more than 5 years were measured in three wells (9428, 25848, and 25851, table 1, figs. 26C and 27B; OWRB, 2019a; USGS, 2020). Well 9428 was in the upgradient part of reach 4 of the Washita River aquifer near Alex. Simulated water-table altitudes generally matched the magnitude and pattern of observed water-table altitudes in that well. Wells 25848 and 25851 were in reach 4 of the Washita River aquifer near Maysville, Okla. The simulated water-table altitudes for well 25848 were comparable in pattern but were less variable in magnitude than the observed water-table altitudes. The simulated water-table altitudes for well 25851 generally matched the pattern of but were about 6 ft higher than the observed water-table altitudes.

The mean calibrated residual for base-flow observations at the Alex gage was −7.7 ft³/s (table 9), indicating that, on average, simulated base flows were slightly higher than observed base flows. The combined (steady-state and transient simulations) Alex gage base-flow RMSE was 7.9 ft³/s, and the interquartile range of residuals was 93.9 ft³/s. The mean calibrated residual for base-flow observations at the Pauls Valley gage was 45.6 ft³/s (table 9), indicating that, on average, simulated base flows were lower than observed base flows. The combined Pauls Valley gage base-flow RMSE was 10.8 ft³/s, and the interquartile range of residuals was 213.8 ft³/s. The mean calibrated residual for base-flow observations at the downgradient end of the model was −8.8 ft³/s (table 9), indicating that, on average, simulated base flows were slightly higher than observed base flows. The combined base-flow RMSE at the downgradient end of the model was 8.0 ft³/s, and the interquartile range of residuals was 234.1 ft³/s. Although simulated base flows generally matched the pattern of observed base flows, they often failed to reach the high and low extremes (fig. 28). For example, the simulated base flows failed to reach the highest observed base flows in 2007 at the Pauls Valley gage and at the downgradient end of the model (fig. 28B, C). Likewise, the simulated base flows failed to reach the low (sometimes 0 ft³/s) observed base flows during the 2011–14 drought (fig. 4) at all locations on the Washita River, especially the downgradient end of the model (fig. 28). The failure to simulate extreme values in the observed data is typical of groundwater-flow models (Ellis and others, 2017; Smith and others, 2017, 2021) and probably results from numerical-model discretization, or the necessary simplification of spatially and temporally variable hydrologic processes that actually occur (Mandelbrot, 1983).

Calibrated Numerical-Model Water Budget

The calibrated numerical-model water budget (table 10, fig. 16B) shows simulated mean annual inflows and outflows for the 1980–2017 modeling period; a subaccounting for reach 3 and reach 4 of the Washita River aquifer was computed by using the ZONEBUDGET utility (Harbaugh, 1990). Simulated recharge (121,613 acre-ft/yr) was the primary inflow for the calibrated numerical model and accounted for 88.9 percent of inflows in the water budget. Net lateral groundwater flow accounted for the other 11.1 percent of inflows. Net streambed seepage (54.3 percent of outflows in the water budget) was the largest outflow for the calibrated numerical model; saturated-zone evapotranspiration accounted for 43.5 percent, and well withdrawals accounted for 0.7 percent of outflows (table 10, fig. 16B). The inflow components of the calibrated numerical-model water budget compare well with those of the conceptual-model water budget (table 6, fig. 16). The saturated-zone evapotranspiration and net streambed seepage outflow components, however, differed substantially between the calibrated numerical-model and conceptual-model water budgets (table 6, fig 16). These components of the conceptual model, especially the saturated-zone evapotranspiration component, were poorly constrained by available field data. No field data were available to quantify saturated-zone evapotranspiration. The types of field activities required to effectively model saturated-zone evapotranspiration are expensive to conduct, and the resulting data are difficult to integrate at the scale of the study area. Historical base-flow observations at the downgradient end of the model were estimated from the Dickson gage, which included base-flow contributions from a relatively large tributary (Caddo Creek) and periodic releases from Lake of the Arbuckles (fig. 2); thus, the budgeted outflows for reach 4 were more poorly constrained than those for reach 3. The conceptual-model net streambed seepage also may be biased high because of the inclusion of some outflows from floodwater-retarding structures that could not be distinguished and separated from groundwater-derived base flows.

Table 10.    

Calibrated numerical-model water budget of simulated mean annual inflows and outflows for the Washita River aquifer, southern Oklahoma, 1980–2017.

[All values are in units of acre-feet per year; --, not applicable. Totals may not equal sum of components because of independent rounding. Net water-budget totals are calculated as the difference between outflow and inflow; negative values indicate outflow from the aquifer, and positive values indicate inflow to the aquifer. Net streambed seepage, net lateral groundwater flow, and net change in groundwater storage represent the net effect of aquifer inflows and outflows. Negative net change in groundwater storage indicates gain of groundwater storage to the aquifer, which is reported as an aquifer outflow in the numerical groundwater-flow model mass balance. Reach 3 and Reach 4 are upgradient and downgradient, respectively, from the vicinity of U.S. Geological Survey streamgage 07328100 Washita River at Alex, Okla. (fig. 1)]

Table 10.    Calibrated numerical-model water budget of simulated mean annual inflows and outflows for the Washita River aquifer, southern Oklahoma, 1980–2017.

		Reach 3	Reach 4	Total	Percentage of water budget
Modeled area	in cells	4,948	7,231	12,179	--
	in acres	72,697	106,241	178,938	--
	in percent	40.6	59.4	100.0	--
Transient simulation
Water-budget component	Inflow
	Recharge	46,073	75,540	121,613	88.9
	Net lateral groundwater flow	11,331	3,810	15,141	11.1
	Total inflow	57,404	79,350	136,754	100.0
	Outflow
	Net streambed seepage	−33,220	−41,068	−74,288	54.3
	Saturated-zone evapotranspiration	−23,045	−36,401	−59,446	43.5
	Well withdrawals	−206	−795	−1,001	0.7
	Net change in groundwater storage	−1,082	−933	−2,014	−1.5
	Total outflow	−57,553	−79,197	−136,749	100.0

Simulated annual net changes in groundwater storage reflect climatic trends during the 1980–2017 modeling period (fig. 29). Groundwater storage, for the purposes of this report, is the volume of water that would drain from the aquifer under gravity; this definition excludes the water that adheres to the aquifer matrix and water that remains locked in isolated or unconnected pore spaces. Simulated groundwater storage was calculated by multiplying the calibrated specific yield (Rogers and others, 2023) by the simulated saturated thickness of each active model cell (see section “Simulated Saturated Thickness and Groundwater Storage”). The cumulative net change in simulated groundwater storage (surplus or deficit), referenced to the simulated groundwater storage from the steady-state simulation (811,037 acre-feet [acre-ft]), was between −5 and 5 percent in most years. The greatest simulated groundwater storage surplus (10.8 percent) occurred in 2015 after a year of record recharge for the period 1980–2017. The greatest simulated groundwater storage deficit (−8.4 percent) occurred in 2016. The 2016 recharge was below but not exceptionally below the mean for the period 1980–2017; the recharge during 6 other years from 1980 to 2017 was less than the 2016 recharge. The large change in simulated groundwater storage from 2015 to 2016 could be explained by increased inflow from daily recharge received late in 2015 that was not fully discharged as streambed seepage and saturated-zone evapotranspiration until early in 2016. Whatever the cause, the occurrence of the greatest simulated groundwater storage surplus and deficit in consecutive years illustrates the potential for large changes in groundwater storage over short time frames in the study area; between 2015 and 2016, the simulated groundwater storage changed by about 155,000 acre-ft. The 2016 groundwater storage deficit was fully regained in 2017, and the transient simulation modeling period ended with a slight (0.4 percent) groundwater storage surplus as compared to the steady-state simulation (fig. 29).

Simulated Saturated Thickness and Groundwater Storage

The simulated saturated thickness of the Washita River aquifer was determined by subtracting the altitude of the aquifer base from the simulated water-table altitude at the end of the modeling period (2017). The simulated saturated thickness was more than 100 ft along the assumed paleochannel of the Washita River near Verden, Okla. (fig. 30). The simulated mean aquifer thickness (sum of saturated and unsaturated thicknesses) was 54.28 ft (Rogers and others, 2023), and the simulated mean saturated thickness was 47.12 ft (table 11). The simulated mean saturated thicknesses for reaches 3 and 4 were 50.03 and 45.13 ft, respectively (table 11). A simulated mean transmissivity of 2,276 ft²/d was calculated as the mean of the calibrated horizontal hydraulic conductivity multiplied by the saturated thickness of each cell; the simulated mean transmissivities for reaches 3 and 4 were 2,061 and 2,423 ft²/d, respectively. The simulated groundwater storage of the Washita River aquifer at the end of the modeling period (2017) was 869,053 acre-ft; 436,807 acre-ft (50.3 percent) of that total was in reach 3, and 432,247 acre-ft (49.7 percent) of that total was in reach 4 (table 11).

Equal-Proportionate-Share Pumping Rates

EPS scenarios for reaches 3 and 4 of the Washita River aquifer were run for periods of 20, 40, and 50 years. The 2017 simulated water table from the calibrated numerical model was used as the starting water table in each EPS scenario. Annual stress periods were used in these scenarios instead of monthly stress periods to simplify the analysis, improve model stability, and lessen run times. Model inputs for recharge, saturated-zone evapotranspiration, and most stream inflows to the active modeled area were configured as the means of each annual period used in the calibrated numerical model; stream inflows from the Alex gage were added for the Washita River at the boundary between reach 3 and reach 4 (which has no specified Washita River inflow in the calibrated numerical model) so that results from reach 3 and reach 4 could be compared. The General-Head Boundary package (Harbaugh and others, 2000) was disabled for the EPS scenarios to prevent the input of unrealistically high simulated flows from adjacent bedrock units. All simulated wells from the calibrated numerical model, including wells with prior-right permits, were discarded, and a hypothetical well was placed in each active cell (covering about 14.7 acres). The aquifer top multiplier was reduced to 1.00 for the EPS scenarios; these extreme pumping scenarios are sensitive to the aquifer top altitude because well pumping is automatically reduced when the ratio of saturated thickness to aquifer thickness (aquifer top minus aquifer base altitude) approaches zero. The adjustment of specified well pumping is necessary because the full pumping rate is not possible in cells with minimal saturated thickness. The MODFLOW-2005 Well package PHIRAMP variable (Niswonger and others, 2011), which specifies the ratio of saturated thickness to aquifer thickness at which to begin reducing pumping rates, was set to 0.05 (5 percent of aquifer thickness) for the EPS scenarios. Lower PHIRAMP settings were tested but caused model instability.

PEST++ (White and others, 2020), a software package for iterative parameter estimation, was used to determine the EPS pumping rate for the selected period (20, 40, or 50 years). In each PEST++ iteration, the hypothetical wells were assigned a uniform pumping rate for the selected period. If a saturated thickness of at least 5 ft was determined for more than 50 percent of the active cells at the end of an iteration, successive iterations with an increased uniform pumping rate were performed until a saturated thickness of less than 5 ft was determined for at least 50 percent of the cells.

The EPS pumping rates for reach 3 and reach 4 were computed independently, although they were simulated with the same model. The EPS scenarios were sensitive to the rate of Washita River main-stem inflows, which can be highly variable and would likely cease at full EPS development. Therefore, an alternate version of the EPS scenarios without the Washita River main-stem inflows was provided in this report for comparison. To account for potential climate variability, additional EPS scenarios were run with recharge increased and decreased by 10 percent.

PHIRAMP complicates the reporting of EPS pumping rates because the actual (reduced) EPS pumping rate is not uniform across cells and may be substantially less, on average, than the nominal EPS pumping rate specified in the Well package file (Harbaugh, 2005). The nominal and actual EPS pumping rates are provided in table 12, but only the actual EPS pumping rates are discussed in this section of the report. With Washita River main-stem inflows, the 20-, 40-, and 50-year EPS pumping rates under normal recharge conditions were about 3.08 (acre-ft/acre)/yr for reach 3 and about 3.80 (acre-ft/acre)/yr for reach 4 (table 12). Given the respective 72,697- and 106,241-acre modeled areas, these EPS pumping rates correspond to annual yields of about 224,000 acre-ft/yr for reach 3 and about 404,000 acre-ft/yr for reach 4. For the 20-, 40- and 50-year EPS scenarios with Washita River main-stem inflows, decreasing and increasing recharge by 10 percent resulted in a 1.5–2.0-percent change in the EPS pumping rate (table 12). Without Washita River main-stem inflows, the 20-year EPS pumping rate under normal recharge conditions was about 1.14 (acre-ft/acre)/yr for reach 3 and about 1.23 (acre-ft/acre)/yr for reach 4 (table 12). Without Washita River main-stem inflows, the 40- and 50-year EPS pumping rates under normal recharge conditions were about 1.13 (acre-ft/acre)/yr for reach 3 and about 1.23 (acre-ft/acre)/yr for reach 4 (table 12). Given the respective 72,697- and 106,241-acre modeled areas, these EPS pumping rates correspond to annual yields of about 82,000 acre-ft/yr for reach 3 and about 131,000 acre-ft/yr for reach 4. For the EPS scenarios without Washita River main-stem inflows, decreasing and increasing recharge by 10 percent resulted in about a 5.3–5.9-percent change in the EPS pumping rate (table 12).

The EPS scenarios with Washita River main-stem inflows reached equilibrium (less than 0.1 percent annual change in storage) after 6 years, so the 20-, 40-, and 50-year scenarios were nearly identical (fig. 31A–C). At the end of all EPS scenarios, the simulated groundwater storage had decreased by about 59 percent, and simulated base flows on the Washita River at the Alex gage, the Pauls Valley gage, and the downgradient end of the model had all substantially decreased (fig. 31). The Washita River likely would be dry (0 ft³/s) at all locations if the mean main-stem and tributary base flows had not been supplied at the SFR2 inflow locations (figs. 31D–F and 32). After 20 years of EPS pumping, only areas along streams with specified SFR2 inflows retained more than 5 ft of saturated thickness (fig. 32); areas along those streams remained saturated because they received streambed seepage supplied by tributaries of the Washita River (SFR2 inflows).

The EPS pumping rates computed for the Washita River aquifer with Washita River main-stem inflows were higher than those computed for similar alluvial aquifers in Oklahoma (Ellis and others, 2017, 2020; Smith and others, 2017, 2021). Many hydrologic factors influence the EPS pumping rate, but one of the most important factors is the presence or absence of large areas of elevated terrace deposits. These areas of terrace deposits, which often have minimal saturated thickness as compared to alluvium deposits along the center of the river valley, quickly go dry at the beginning of the EPS scenarios. Because the Washita River aquifer as defined in this report includes few elevated terrace deposits, the saturated alluvium deposits along the center of the river valley must be drained to a greater degree with a greater EPS pumping rate to make half of the aquifer reach a saturated thickness of less than 5 ft.

With Washita River main-stem inflows, the actual EPS pumping rate for the Washita River aquifer (both reaches) stabilized at about 59 percent of the nominal EPS pumping rate (fig. 33A). Without Washita River main-stem inflows, the actual EPS pumping rate for the Washita River aquifer (both reaches) stabilized at about 71 percent of the nominal EPS pumping rate (fig. 33B).

EPS scenarios represent an extreme theoretical construct in which the aquifer is fully developed with regularly spaced wells (one well per model cell, which is one well for every 14.7 acres of modeled area), with each well pumping at a high rate. This level of development is unlikely to occur, and these EPS scenarios bear no resemblance to the current (2022) level of development in the Washita River aquifer in terms of well distribution and pumping rates. Some parts of the aquifer may be more developed than others, but the current level of development in the aquifer as a whole is far smaller than the level of development simulated in the EPS scenarios. The total groundwater use for 2017 (1,402.1 acre-ft/yr, table 4) (Rogers and others, 2023), if divided by the modeled area (178,938 acres), is equivalent to an aquifer-averaged pumping rate of about 0.0078 (acre-ft/acre)/yr—less than 1 percent of the 20-year EPS pumping rate for the Washita River aquifer with or without Washita River main-stem inflows.

Projected 50-Year Pumping

Projected 50-year pumping scenarios were used to simulate the effects of modified well withdrawal rates on groundwater storage in the Washita River aquifer and base flows in the Washita River. Monthly stress periods were used in these scenarios. Well withdrawals in these scenarios, unlike in the EPS scenarios, used historical well withdrawal rates (or multipliers on historical well withdrawal rates) and historical well locations used in the transient simulation of the calibrated numerical model. The effects of modified well withdrawals were evaluated by quantifying differences in groundwater storage and base flow in four 50-year scenarios, which applied (1) no groundwater pumping, (2) 2017 pumping rates, (3) increasing demand pumping rates at simulated wells, and (4) mean pumping rates for the study period (1980–2017). The increasing demand pumping rates assumed a cumulative 45.7-percent increase (0.756-percent compounding annual increase) in pumping over 50 years based on 2010–60 demand projections for southern Oklahoma (OWRB, 2012b). Other model stresses were configured as the means of each monthly stress period from the calibrated numerical model, and the scenarios assumed that future climate conditions were comparable to the 1980–2017 study period. The simulated water table from the end of the calibrated numerical modeling period (2017) was used as the starting water table in each projected pumping scenario covering the period 2018–67.

Because well withdrawals were a small part (less than 1 percent) of the calibrated numerical-model water budget (table 10), changes to the well pumping rates had little effect on simulated Washita River base flows (fig. 34) and groundwater storage (table 13) in the Washita River aquifer. For convenience, all groundwater storage changes were referenced to the groundwater storage at the end of the 2017-pumping-rate scenario (876,165 acre-ft, table 13). Groundwater storage after 50 years with no pumping was 876,593 acre-ft, or 428 acre-ft (0.049 percent) greater than the groundwater storage at the end of the 2017-pumping-rate scenario; this groundwater storage increase is equivalent to a mean water-table-altitude increase of 0.023 ft. Groundwater storage at the end of the 50-year period with the increasing demand pumping rates was 876,030 acre-ft, or 135 acre-ft (0.015 percent) less than the storage at the end of the 2017-pumping-rate scenario; this groundwater storage decrease is equivalent to a mean water-table-altitude decline of 0.007 ft. Groundwater storage after 50 years of pumping at the mean rate for the study period (1980–2017) was 874,807 acre-ft, or 1,358 acre-ft (0.155 percent) less than the groundwater storage at the end of the 2017-pumping-rate scenario; this groundwater storage decrease is equivalent to a mean water-table-altitude decline of 0.074 ft. Though changes in groundwater storage were minimal, some cells near well boundaries showed increased development of simulated cones of depression. For example, local simulated drawdown in some cells was about 3 ft at the end of the 50-year period with the increasing demand pumping rates (Rogers and others, 2023).

Hypothetical 10-Year Drought

A hypothetical 10-year drought scenario was used to simulate the potential effects of a prolonged period of reduced recharge on groundwater storage. The period January 1993–December 2002 was chosen as the hypothetical drought period because it had a mean annual recharge rate (8.3 in/yr) similar to that of the 1980–2017 study period (8.1 in/yr, table 7). Drought effects were quantified by comparing the results of the drought scenario to those of the calibrated numerical model (no drought) at the end of the hypothetical drought period (2002). To simulate the hypothetical drought, recharge in the calibrated numerical model was reduced by 50 percent during the hypothetical drought period (1993–2002), and upstream inflows from the Washita River and simulated tributaries were reduced by 75 percent, which was comparable to the reduction in Washita River base flows at the Anadarko gage during the 2010–14 drought period (fig. 6A).

Simulated groundwater storage with drought at the end of the hypothetical 10-year drought period (stress period 277) was 826,141 acre-ft, which is 39,487 acre-ft (4.6 percent) less than the simulated groundwater storage of the calibrated numerical model (with no drought) at the end of the drought period (865,627 acre-ft, table 14, fig. 35). This decrease in simulated groundwater storage is equivalent to a mean water-table-altitude decline of 2.2 ft (table 14). The largest simulated water-table-altitude declines occurred in the terrace areas most upgradient from the Washita River (Rogers and others, 2023). The simulated saturated thickness of areas near the Washita River and major tributaries changed little during the hypothetical drought, but the simulated base flow in streams in those areas decreased rapidly. After 12 months of the hypothetical drought, simulated base flows in the Washita River at the Alex gage, Pauls Valley gage, and downgradient end of the model had all decreased by greater than 60 percent as compared to the calibrated numerical model (fig. 36). At the end of the hypothetical 10-year drought period (120 months), simulated base flows at the Alex gage, Pauls Valley gage, and downgradient end of the model had decreased by about 72, 68, and 71 percent, respectively (fig. 36). After the drought period, simulated base flows mostly recovered (returned to less than a 1-percent decrease in simulated base flow) after 2 years (fig. 36), but a groundwater storage deficit of about 1 percent remained until the end of the simulation (fig. 35).