Surface Water

Primary surface-water features in the Cache region are the southerly flowing White River, Cache River, Bayou De View, and L’Anguille River (fig. 2A). The White River, with an average daily discharge of about 21,334,139 cubic meters per day (m³/d) (8,720 cubic feet per second [ft³/s]), is the largest stream in the region and flows along the western side of the regions, draining into the Mississippi River in southeast Arkansas (Schrader and others, 2005). The Cache River, located east of the White River, which has an average daily discharge of about 2,446,576 m³/d (1,000 ft³/s) and flows south from near Egypt, Ark., to its confluence with the White River in the Grand Prairie region (fig. 2A and 2B). Bayou De View flows through the central part of the region, discharging to the Cache River in the Grand Prairie region, and the L’Anguille River flows south through the eastern part of the region, just west of Crowley’s Ridge, which is the headwaters region for tributaries to the L’Anguille River (fig. 2A).

Primary surface-water features in the Grand Prairie region include the White River in the western region, Wattensaw Bayou, Bayou Two Prairie, and Bayou Meto in the central Grand Prairie region, and the Arkansas River that flows southeast from the Ozark Plateau to the Mississippi River in the southeastern corner of the region (fig. 2B). Peak flows in this area occur from March through April with low flows occurring October through November. Pumping of groundwater for irrigation has reduced base flows and increased the frequency of low flows for many rivers and streams in this area (Yasarer and others, 2020). The flow alteration is greatest in October owing to the cumulative withdrawal of groundwater (Yasarer and others, 2020). The drainage system in the MAP has changed in the last 80 years owing to changes in land use as well as drainage improvement projects, flood controls, the construction of levees and ditches, and the deepening and straightening of streams.

Groundwater and Hydrogeologic Units

Groundwater is present in several major hydrogeologic units in the northeastern Arkansas region of the MERAS (Hart and others, 2008). These major units, from shallowest to deepest, include the MRVA, Vicksburg-Jackson confining unit, the upper Claiborne aquifer, the middle Claiborne confining unit, the middle Claiborne aquifer, the lower Claiborne confining unit, the lower Claiborne aquifer, the middle Wilcox aquifer, the lower Wilcox aquifer, and the Midway confining unit (Hart and others, 2008). Detailed descriptions of each hydrogeologic unit for Arkansas is available in Kresse and others (2014). Hydraulic properties of primary aquifers and confining units in the Cache River and Grand Prairie regions from previous studies are summarized in table 1.

The MRVA is the shallowest aquifer and primary source of fresh groundwater in the MAP (Kresse and others, 2014). Groundwater flow in the MRVA is generally from north to south but has changed in some regions owing to intensive pumping. In the Cache region, owing to decades of intensive groundwater pumping for irrigation and aquaculture, a cone of depression developed between Crowley’s Ridge and the Cache River (fig. 2A; Czarnecki, 2010; McGuire and others, 2019). Similarly, in the Grand Prairie region, groundwater declines have been recorded since the onset of irrigation pumping in the early 1900s with significant declines resulting in the depletion of the MRVA in the Grand Prairie CGWA by the late 20th century (fig. 2B; Clark and others, 2011a). Depletion of the MRVA primarily owing to irrigation pumping has caused increased water use in the deeper Sparta aquifer. Although not well quantified, low yields from wells screened in the MRVA in the Grand Prairie region have led farmers to install wells screened in the deeper Sparta aquifer so they can continue irrigating their crops (Kresse and others, 2014).

The thickness of the MRVA varies from 0 to 60 m throughout the Cache and Grand Prairie regions (Clark and Hart, 2009). The transmissivity of the MRVA in northeastern Arkansas ranges from 2,591 to 19,782 meters squared per day (m²/day) (Broom and Lyford, 1981). The MRVA horizontal hydraulic conductivity (Kh) ranges from 30.5 to 223 meters per day (m/d), specific yield (SY) values for unconfined conditions were about 0.19 to 0.3 (unitless), and specific storage (SS) values for confined conditions were 1.2x10^-6 to 6x10^-3 1/m (Broom and Lyford, 1981; Czarnecki and others, 2003; Clark and Hart, 2009). Vertical anisotropy of hydraulic conductivity (Kvani), the ratio between Kh and vertical hydraulic conductivity (Kv) (Kh:Kv, unitless), ranged from 10 to 100 (table 1; Reed, 2003; Clark and Hart, 2009).

The Vicksburg-Jackson confining unit, where it exists in northeastern Arkansas in the Grand Prairie region, reduces the hydraulic connection between the MRVA and underlying middle Claiborne aquifer. The Kh, Kvani, and SS of the Vicksburg-Jackson confining unit is about 0.3 m/d; 1,475 m/d; and 1.1 x10^-6 1/m, respectively (Clark and Hart, 2009). The range of values from several sources is summarized in table 1.

The Tertiary-age aquifers below the MRVA include the Claiborne and Wilcox aquifers (Kresse and others, 2014). In northeastern Arkansas, the primary water use in the middle Claiborne aquifer occur in the undifferentiated Memphis Sand and Sparta Sand (Kresse and others, 2014). The Cache River Formation and Carrizo Sand are present in northeastern Arkansas (Clark and Hart, 2009). Locally and in various hydrologic studies, the middle Claiborne aquifer is referred to as the Memphis aquifer, the Sparta aquifer, or the Sparta-Memphis aquifer; these various names can add to the naming confusion, as described in Kresse and others (2014). Hereafter, in this report, the Sparta, Memphis, and Sparta-Memphis aquifers are referred to as the “middle Claiborne aquifer.” The middle Claiborne confining unit, above the middle Claiborne aquifer, had a Kh value of about 0.05 m/d and a Kvani of about 565 (same as table 1) (Clark and Hart, 2009). The middle Claiborne aquifer, in northeastern Arkansas, had Kh values from 1.1 to 44.5 m/d (Clark and Hart, 2009); middle Claiborne aquifer Kvani was about 207.2 and SS was about 1x10^-6 to 2.6x10^-5 1/m. The lower Claiborne aquifer had a Kh value of about 7.6 m/d and Kvani was about 23 (Clark and Hart, 2009). The middle Wilcox aquifer Kh was about 0.7 m/d, the Kvani was about 617, and SS was 1.1x10^-6 (Clark and Hart, 2009). The lower Wilcox aquifer Kh was about 7.5 m/d, the Kvani was about 28, and SS was 1.1x10^-6 (Clark and Hart, 2009). The Midway confining unit, consisting of thick marine clays, underlies both regions, is assumed to contribute insignificant amounts of water to the overlying aquifer units, and is considered to be a no-flow boundary where it contacts the Wilcox aquifer system across the study area in northeastern Arkansas (Clark and Hart, 2009). Property values and their ranges from previous studies are also summarized in table 1.

Cache and Grand Prairie Model Domains

A combination of the 2016 MRVA potentiometric surface map from McGuire and others (2019) and major streams were used to delineate the Cache and Grand Prairie domains for the numerical models (hereafter referred to as the “Cache model domain” and the “Grand Prairie model domain”; fig. 2A, B). Each model domain was delineated to ensure focus areas of each numerical model were minimally affected by boundary conditions and to maximize natural no-flow boundaries where groundwater flow was parallel to the domain boundaries (Anderson and others, 2015). The focus areas for each model domain were the primary cones of depression (areas of largest water-level declines compared to the surrounding water table) underlying the Cache River and Grand Prairie CGWAs (fig. 2A, B).

The northern boundary of the Cache model domain extends from the Ozark Plateau and is just west of the MERAS fall line to just east of Crowley’s Ridge (fig. 2A). The groundwater-flow direction (perpendicular to equipotential lines in fig. 2A) along the northern Cache model domain boundary is more parallel to the boundary, flowing eastward until near Crowley’s Ridge where groundwater flow turns perpendicular and flows south along Crowley’s Ridge (fig. 2A). Flow parallel to the domain boundary indicates a no-flow boundary. Selection of the southern domain boundary applies the same principle as the northern boundary where groundwater flow is parallel to the domain boundary, approximating a no-flow regime. The eastern edge of the domain was chosen as the easternmost extent of Crowley’s Ridge, which acts as a barrier to groundwater flow as seen by the head gradient difference in the MRVA across the Crowley’s Ridge (fig. 2A; McGuire and others, 2019). The western edge of the domain was chosen where there was some parallel groundwater flow in the southwestern region (fig. 2A).

The northern boundary of the Grand Prairie model domain extends from the fall line between the MERAS and Ozark Plateau (fig. 2B) to St. Francis County (fig. 2B). The northern domain boundary is north of a cone of depression in northern Prairie County and about 20 km north of the primary cones of depression in the CGWA. The southern domain boundary, which extends from the forested lands in Cleveland County east such that the Arkansas River was inside the domain until near its confluence with the Mississippi River (fig. 2B). The eastern edge of the domain was chosen to include most of the White River through Monroe and Arkansas Counties (fig. 2B; McGuire and others, 2019). The western edge of the domain was chosen to encompass the Arkansas River where it flowed out of the Ozark Plateau and into the MERAS (fig. 2B).

Hydrostratigraphy from Geophysical Data

The resistivity data from AEM and waterborne geophysical surveys were used to further the understanding of the hydrostratigraphy of the Cache and Grand Prairie model domains to improve the conceptual hydrogeologic framework developed in Hart and others (2008). A brief description follows, but for more detail please refer to Minsley and others (2021). In the Cache model domain, the shallow alluvial system exhibited areas of low and high resistivity values horizontally and vertically throughout the MRVA, which indicated the MRVA was more heterogeneous than was represented in previous models that had specified single or uniform Kh, Kv, SS, and SY values. The middle Claiborne aquifer is a higher resistivity unit directly underlying the MRVA in the eastern region of the Cache extending under Crowley’s Ridge (fig. 4A). At a depth of about 40 to 70 m below the land surface, underlying Crowley’s Ridge in the central and southern area of the Cache model domain, the Cook Mountain Formation comprises the middle Claiborne confining unit and underlies the unconsolidated sediment cap of Crowley’s Ridge. The transition area between the prominent Midway confining unit and the middle Claiborne aquifer is the Wilcox aquifer system, which consists of the middle Wilcox and lower Wilcox aquifers; however, these are difficult to differentiate in the Cache model domain and can be considered a single aquifer system. The Midway confining unit was prominent in the resistivity data where it subcropped beneath the MRVA in the west and central MAP with an eastward dip (table 1.1).

The Grand Prairie model domain, like the Cache model domain, exhibited a heterogeneous resistivity field in the MRVA that provided more detail on the potential Kh and Kv values throughout the aquifer. The middle Claiborne aquifer directly underlies the MRVA in the northern half of the region (Minsley and others, 2021). The Vicksburg-Jackson confining unit separates the MRVA and middle Claiborne aquifer in the southern half of the region; however, based on the resistivity classes, there are some small areas of higher resistivity within the Vicksburg Jackson confining unit that may allow for some hydraulic connection between the MRVA and middle Claiborne aquifer (fig. 4B). For the deepest areas resolved by the AEM data, parts of the middle Claiborne aquifer and lower Claiborne confining unit exhibit low resistivities, which represent less hydraulically conductive material.

Groundwater Inflows, Outflows, and Storage Changes

Inflows to both model domains included recharge from precipitation, recharge as seepage from streams, recharge return flows from irrigation, cross-boundary groundwater flow from upgradient adjacent areas, and releases from groundwater storage. Recharge from precipitation was estimated to be 0–5.73 in/yr with some small areas exceeding 20 in/yr (Reed, 2003; Clark and Hart, 2009; Reitz and others, 2017). Recharge as seepage from streams is affected by streambed properties such as Kv, can vary along a stream reach, and was estimated to be about 15–25 percent of recharge from precipitation, or as much as about 1.5 in. (Reed, 2003). Recharge as return flow from traditional gravity flood or furrow irrigation water has not been well quantified by previous studies but generally occurs where shallow clay layers in the MRVA are absent (Yaeger and others, 2018).

Outflows to both regions include irrigation and aquaculture pumping from wells, pumping from municipal wells, base flow discharge to streams, cross-boundary groundwater flow to downgradient adjacent areas, and replenishment of groundwater storage. The primary outflow is groundwater pumping for irrigation and to a lesser extent for aquaculture, and municipal from wells screened in the MRVA (Clark and others, 2011a). Evapotranspiration from groundwater by plant roots that reach the water table was negligible owing to widespread water-level declines in the 20th century, and it was not considered in past groundwater-flow models developed in the region (Reed, 2003; Clark and Hart 2009; Rashid and others, 2015).

Groundwater pumping is widespread in both model domains; however, to reduce the use of groundwater for crop irrigation and limit groundwater declines in areas with less saturated thickness in the MRVA, some surface-water irrigation exists in the Cache and Grand Prairie model domains (Yaeger and others, 2018). Surface-water irrigation occurs as diversions of streamflow to on-farm irrigation reservoirs that range in size from less than 10 acres to more than 300 acres (Yaeger and others, 2018). Diversions can occur from small bayous to large rivers. Construction and utilization of on-farm irrigation reservoirs has been ongoing since the mid-1970s with an increasing abundance of reservoirs through recent times (Yaeger and others, 2018). These on-farm storage reservoirs are shallow, typically have a depth of 2.4–3 m, and are filled to capacity in the wetter winter and early spring months (January through March). Farmers typically use the water stored in these reservoirs to furrow or flood irrigate their crops in the drier summer months (June through August) when evapotranspiration exceeds precipitation (Yaeger and others, 2018). Farms that have irrigation reservoirs typically also have irrigation wells that are assumed to only be used to irrigate their crops if the reservoirs do not have enough water stored because of a dry winter and spring when diverting streamflow was not possible.

Hydrologic Boundaries

Hydrologic boundaries, with respect to the groundwater-flow system, exist where groundwater is exchanged with the landscape or surface-water system, or where groundwater flows in or out of the model domain. Hydrologic boundaries can be used to determine the extent of groundwater flow simulated in a numerical model. The hydrologic boundaries for the Cache and Grand Prairie models included aquifer boundaries and streams. The western extent of the MERAS was the main hydrologic boundary for northwestern and west-central part of the Cache model domain and for the northwestern and western part of the Grand Prairie model domains. The interface between the MRVA and the Ozark Plateau aquifer system was a boundary in the northwestern section of both model domains (fig. 5A, B). The Ozark Plateau aquifer system has a moderate hydraulic connection with the MRVA where groundwater generally flows from the higher elevation Ozark Plateau aquifer to the lower elevation MRVA; however, the contribution of the Ozark Plateau aquifer to the overall MERAS water budget is negligible (Clark and Hart, 2009; Clark and others, 2018).

Horizontal boundaries of groundwater flow exist along all four sides of each model domain except where the Ozark Plateau aquifer is absent. Specified fluxes, derived from MERAS 3 cell-by-cell fluxes, were assigned to the outermost model cells for all model layers using the Well (WEL) package (fig. 5; Harbaugh, 2005). Within model cells where the Ozark Plateau aquifer system intersected the MRVA, specified heads were assigned using the General Head Boundary (GHB) package (fig. 5A, B; Harbaugh, 2005). Hydraulic conductances assigned to the GHB cells were based on values in the Ozark model interface with the MRVA from Clark and others (2018).

The upper vertical boundary was the land surface, which was defined using a digital elevation model (U.S. Geological Survey, 2016). The lower vertical boundary was the top of the Midway confining unit represented as a vertical no-flow groundwater boundary in both models because the Midway confining unit is not hydrologically connected with the MERAS (Clark and Hart, 2009; Traylor and Weisser, 2024).

Spatial and Temporal Discretization

The Cache and Grand Prairie models were spatially discretized horizontally into 500-meter x 500-meter orthogonal cells on a grid with the following vertical layers, rows, and columns:

Vertical layering schemes used the same approach as described in Leaf and others (2023) whereby layers 1 through 8 in both models used the AEM data at 5-m constant thicknesses to represent the MRVA, layers 9 through 15 used increasing uniform thickness layers to represent the Vicksburg-Jackson confining unit and, if present, the Claiborne aquifers and confining units. Below layer 15, the MERAS framework layers from Hart and others (2008) were used to represent the bottom surfaces of Wilcox aquifers and confining units (fig. 4A, B).

The bottom of layer 8, a depth of 40 m (eight layers at 5-m thickness) was represented using the updated bottom of the MRVA, which was produced by James and Minsley (2021). Where the MRVA was less than a depth of 40 m, the remaining layers to the 40-m depth were assigned as pass-through cells of zero thickness (refer to “Regional and Inset Model Layering” section in appendix 1 of Leaf and others [2023] for a more detailed description). Layers 9 through 15 increased in thickness for both models at a multiplier of 1.5 where layer 15 was the lowest extent of the AEM data resistivity facies classes. The framework layers from Hart and others (2008) delineated the bottom elevations for layers 16 through 18 of the Cache model and layers 16 through 19 of the Grand Prairie model (fig. 4A, B).

Based on the resistivity facies classes from the AEM data, the MRVA thickness in both models varied from about 20 to 40 m (fig. 4A, B). Most cells were active in layers 1 through 4, because the MRVA was present almost everywhere at those depths. Pass-through cells, a new feature in MODFLOW 6, were assigned to cells where units did not exist; however, the overlying and underlying cells were active. The MRVA progressively thinned and became absent in layers 5–8 (fig. 6A, B). In layer 8, for both models, less than 50 percent of the cells were active (no pass-through) cells (fig. 6A, B).

Temporal discretization for both models included 148 stress periods of varying lengths to simulate steady state and transient hydrologic conditions from January 1, 1900, through December 31, 2018. The first stress period simulated predevelopment steady-state conditions (without pumping), followed by six multiyear stress periods that simulated transient conditions from January 1, 1900, to April 1, 2007. The multiyear stress periods start and end dates aligned with the previous MERAS models, which had stress periods defined to account for major changes in pumping history (fig. 10 in Clark and Hart, 2009). Monthly stress periods simulated the transient hydrologic conditions from April 1, 2007, through December 31, 2018 (table 2.1; table 2.1 shows end date of January 1, 2019, because that is the stop time of the model, but that date is not simulated).

MODFLOW Inputs and Configuration

Recharge inputs to the Cache and Grand Prairie models were simulated using the Recharge (RCH; Harbaugh, 2005) package. Recharge was specified as transient array-based datasets using the values from a soil-water-balance (SWB; Westenbroek and others, 2018) model of the MERAS developed for the MAP project (Westenbroek and others, 2021). The recharge outputs from the SWB model, which had a cell resolution of 1,000 m to match the MERAS 3 resolution, were clipped and resampled to the size of the Cache and Grand Prairie model domains (Traylor and Weisser, 2024). The unsaturated zone was not simulated in the Cache model or Grand Prairie model because that functionality was not available in MODFLOW 6 at the time of model development, and it was not necessary to address the objective of the study.

Streams were simulated in both models with the Streamflow Routing (SFR; Niswonger and Prudic, 2005) package except for the Black River, White River, and Little Red River (not shown) in the Cache model and the Arkansas River in the Grand Prairie model, which were simulated with the River (RIV; Harbaugh, 2005) package (fig. 5A). Stream reaches simulated in each model were delineated from the NHDPlus version 2 dataset (McKay and others, 2012) using python scripts with SFRmaker (Leaf and others, 2021) and MODFLOW setup (Leaf and Fienen, 2022) as described in the workflow from Leaf and others (2023). The SFR package simulated total stress period averaged streamflow in both models through the simulation of groundwater/surface-water exchange, accumulation of flow from upstream reaches and stream inflows from the edge of the model, and inflows of runoff generated by the SWB model from Westenbroek and others (2021) (fig. 5A, B). Daily runoff values from the SWB model were aggregated to stress period averaged values and assigned using the methods described in Leaf and others (2023).

Monthly stream inflow values were assigned using monthly average streamflow measured at streamgages where available. Inflow points where streamgages were absent or the measured streamflow datasets did not span the entire simulation period were filled with estimated streamflows from a random-forest regression surface-water model by Dietsch and others (2023). The Cache model included five inflow points where only one streamgage, the Cache River at Egypt, Ark. (USGS streamgage 07077380) was close enough to the model boundary and had available data to specify inflows to the SFR package network (fig. 5A). The other four inflow points used estimated monthly average streamflows from the random-forest regression surface-water model (Dietsch and others, 2023). The Grand Prairie model included 23 inflow points and, like the Cache model, only one inflow point used measured streamflows from a streamgage, the Cache River at Cotton Plant, Ark. (USGS streamgage 07077555); this streamgage did not have records prior to 2008, so estimated streamflows from the random-forest regression surface-water model were used to assign inflow values prior to 2008.

Diversions of streamflow from the SFR package were estimated in both models to simulate the removal of stream water to fill irrigation reservoirs (fig. 5A, B; Yaeger and others, 2017). Diversion data records on volumes or locations by each user were not available, so the diversion volumes, rates, and locations were estimated based on general reservoir filling practices and a remote sensing survey, which provided an inventory of irrigation reservoir locations (Yaeger and others, 2017; Michele Reba, U.S. Department of Agriculture-Agricultural Research Service, written commun., 2022). The SFR package diversion reaches were assigned based on proximity to the closest SFR reach to each reservoir. Diversion volumes for each irrigation reservoir were estimated based on the area of each irrigation reservoir from the inventory dataset in Yaeger and others (2017) with the assumption of an 8-ft filling depth. Diversion rates were estimated by equally dividing the diversion volume across typical filling months of February through May (Yaeger and others, 2017).

The SFR package physical properties, which were required in the package data file and were not generated from the NHDPlus version 2 dataset, included streambed thickness, Manning’s roughness coefficient, and streambed hydraulic conductivity (rhk) (Langevin and others, 2022; Traylor and Weisser, 2024). The rhk values were specified as leakance because streambed thicknesses were assigned a uniform value of 1 m for all reaches; streambed thickness is difficult to estimate on a large scale and a value of 1 is commonly used in models, with all the variability in leakance represented by rhk (Leaf and others, 2023). Initial values of streambed leakance were assigned using the values from a waterborne electrical resistivity survey of streams and lakes in the MAP by Adams and others (2019) or the geometric mean of the entire waterborne electrical resistivity survey dataset (Adams and others, 2019) where waterborne surveys were not completed. Streambed top elevations for each reach were assigned using a mean value from a digital elevation model dataset that was resampled to the model grid discretization (U.S. Geological Survey, 2019).

The RIV package was used for those major streams outside the focus CGWAs: the Black River, White River, and Little Red River (not shown) in the Cache model, and the Arkansas River in the Grand Prairie model (fig. 5A, B). The RIV package was used for these streams because they were far enough away from the cones of depression that the simulation of total routed streamflow by the SFR package was not necessary. The RIV package, a head-dependent boundary condition, simulated the groundwater/surface-water exchange based on the hydraulic head gradient among the RIV package stage and the head in underlying groundwater cell and a specified riverbed conductance term (Traylor and Weisser, 2024). The Cache model RIV package stages were interpolated from U.S. Army Corps of Engineers (USACE) stage records at one streamgage along the Black River in the model domain (fig. 5A; USACE streamgage Black River at Elgin Ferry, Ark.) and two USACE streamgages on the White River (fig. 5A; USACE streamgages White River at Augusta, Ark., and White River at Newport, Ark.) (U.S. Army Corps of Engineers, 2020). The Grand Prairie model RIV package stages were interpolated from USACE stage records at three streamgages along the Arkansas River in the model domain (fig. 5B; USACE streamgage Arkansas River at Pendleton, Ark.; Arkansas River at Pine Bluff, Ark.; Arkansas River at David D Terry L&D below Little Rock, Ark.) (U.S. Army Corps of Engineers, 2020).

Groundwater pumping was simulated with the WEL package (Harbaugh, 2005) in both models for irrigation, aquaculture, and municipal wells, and, only for the Grand Prairie model, thermoelectric wells using MODFLOW-setup methods described in Leaf and others (2023). For stress periods 2 through 6 (January 1, 1900, through April 1, 2007), the MERAS 2.2 model (Haugh and others, 2020) was used to specify the pumping rates and locations for all types of pumping. Irrigation and aquaculture, municipal, and thermoelectric pumping was simulated for the monthly stress periods 7 through 148 (April 1, 2007, through December 31, 2018) using MODFLOW-setup (Leaf and Fienen, 2022) and included in the WEL package input datasets. For the monthly stress periods, outputs from the AIWUM version 1.1 (Bristow and Wilson, 2023), which is summarized in the “Water Use” section of this report, were used to specify the locations and rates of irrigation and aquaculture pumping by crop type. When available, site-specific water use data (SWUDS) from 2000 through 2018 (U.S. Geological Survey, 2020) were used to simulate municipal pumping, and national estimates of thermoelectric pumping were provided by Diehl and Harris (2014) and Harris and Diehl (2019a,b). For monthly stress periods (April 1, 2007, through December 31, 2018), when SWUDS and thermoelectric pumping data were not available, the pumping rates were assigned in the WEL package based on pumping rates that were closest in time to the missing data gap.

Vertical layer assignments for irrigation and aquaculture pumping were determined using a geostatistical estimation dataset by Torak (2021), which estimated the irrigation production zone interval of the MRVA. The AEM survey associated with the MAP project was used to delineate a new bottom of the MRVA surface, included in the Cache and Grand Prairie models as the bottom of layer 8, so there is a difference between the two surfaces that may affect pumping placement vertically in the models. Most of the irrigation and aquaculture pumping were initially assigned to model layers 1–10 in both models. However, in the Grand Prairie region, owing to decades of groundwater declines in the MRVA, irrigation wells have gone dry and caused farmers to drill deeper past the MRVA and into the middle Claiborne aquifer (locally the Sparta aquifer) in search of a reliable source of groundwater (Kresse and others, 2014). This practice of installing wells screened in the middle Claiborne aquifer is not well quantified or documented, which precluded the accurate assignment of irrigation and aquaculture pumping in the middle Claiborne aquifer. To capture this practice in the Grand Prairie model, a fraction of irrigation and aquaculture pumping from wells in layers 5–8 was moved to the middle Claiborne aquifer model layer 15 if present and if resistivity facies classes were 3–7, indicating aquifer material. This was done to improve the simulation of the Sparta aquifer in the Grand Prairie model and capture a conceptual element of the system that has not been simulated in previous models. Vertical layer assignments for the nonagricultural pumping were specified using open intervals or estimated production zones from Knierim and others (2019).

Aquifer properties (Kh, Kv, SS, and SY) in both models were specified in the Node Properties File and the Storage (STO) package (Langevin and others, 2022) by way of zones that corresponded to the resistivity facies classes derived from the AEM survey data in Minsley and others (2021) and the original MERAS framework from Hart and others (2008) below the extent of the AEM data. Initial values for Kh, Kvani, SS, and SY were based on previous studies that included calibrated models discussed in the “Previous Studies” and “Groundwater and Hydrogeologic Units” sections of this report. Using the prior knowledge on aquifer properties and the resistivity classes that relate to fine-to coarse-grained material, the Kh values assigned in ascending order related to the resistivity facies class such that the lowest values of Kh corresponded to the lowest facies resistivity class (class 1) and the highest Kh values corresponded to the highest facies resistivity class (class 7). For non-AEM layers (layers 10 through 19), values from previous studies were used initially, then modified during calibration.

The Newton-Raphson formulation within MODFLOW 6 solved the groundwater-flow system of equations as specified in the integrated model solution (Langevin and others, 2022). The “moderate” option for the Cache and Grand Prairie models was employed owing to the simulation of at least one unconfined layer and minimal convergence issues with this setting (Langevin and others, 2022). Simulation run times on a modern laptop computer for the calibrated Cache and Grand Prairie models (refer to “Groundwater-Flow Model Calibration” section) were about 4 and 10 minutes, respectively. The Observations (OBS) package in MODLFOW 6 was also used to retrieve simulated groundwater altitudes for all layers at cells that intersected observation well locations (Langevin and others, 2022).

Calibration Approach

The Cache and Grand Prairie models were both calibrated using the PESTPP–IES version 5 code to automate the adjustment of parameters and improve the fit between measured observations or calibration targets and simulated equivalent model outputs. Manual trial-and-error adjustment of some parameters was done prior to automated calibration to improve the history matching and the understanding of the groundwater systems for each model. These initial values for the automated parameter estimation step are referred to as the “prior” parameter distribution. PESTPP–IES produces a “posterior” parameter distribution for a user-specified number of realizations for each iteration of the calibration (White, 2018). The “posterior” distribution of parameter values are different for each realization and each iteration of the calibration. At the conclusion of an automated PESTPP–IES calibration run, an individual realization from a specific iteration can be chosen to represent the best single combination of “calibrated” parameter values. Further, the ensemble from that chosen “calibrated” iteration can be used to represent the uncertainty or range of parameter values if a single realization is not adequate to address the needs of the study. In this study, the “base” realization was chosen as the single combination of “calibrated” parameter values because it is most similar to the traditional least squares parameter estimation algorithm used in the traditional Parameter Estimation software (Doherty, 2005).

The general calibration approach and workflow for the Cache and Grand Prairie models used the Delta inset model calibration approach from Leaf and others (2023); however, some changes were required in the Cache and Grand Prairie models to achieve acceptable calibration results. Similarities with the Delta inset model calibration approach include utilization of methods in the pyemu Python module to create the PEST++ framework around each model (White and others, 2016). This framework included set-up of template and instruction files and creation of a PEST control file. Some differences included parameter schemes, types of observation, and weighting of observation groups in automated calibration.

Automated calibration was performed using PESTPP–IES (White, 2018; White and others, 2020). The IES in PESTPP–IES is a tool that combines the Gauss-Levenberg-Marquardt (GLM) least-squares parameter estimation algorithm (used in PEST++–GLM; White and others, 2020) and Monte Carlo using the PEST++ interface to reduce an objective function (phi [Φ]), which is the sum of squared and weighted residual (measured minus simulated) values (Welter and others, 2015). PESTPP–IES is advantageous because it approximates the Jacobian (sensitivity) matrix by calculating correlations between an ensemble of parameters and observations, which allows for very high-dimensional parameterization without the computational burden associated with PEST++–GLM (White and others, 2020; Hunt and others, 2021) during an automated calibration run. The IES has been used to successfully calibrate groundwater-flow models in several recent studies (Corson-Dosch and others, 2022; Fienen and others, 2022; Ellis and others, 2023; Leaf and others, 2023) using thousands more parameters than traditional PEST from Doherty (2005). Therefore, the Cache and Grand Prairie models employed high dimensional parameterization schemes of 13,740 and 30,436 parameters, respectively.

Observations

Observations, used as calibration targets in the Cache and Grand Prairie model calibration process, inform PESTPP–IES during the parameter upgrade process. The fit between observations and their simulated equivalent values is an indication of a model’s ability to accurately simulate historic conditions and helps to illuminate bias in the results. In PEST++, similar observations can be assigned to groups where each group’s contribution to Φ can be assessed. The Cache and Grand Prairie models included 13 and 18 observation groups with nonzero weighted observations, respectively, and were grouped based on location and type of observation (tables 2 and 3). The observation groups (tables 2 and 3) for both models included the following:

Most of the primary and derivative observation groups were based on those used in Leaf and others (2023).

Weights were assigned to each observation in both models based on each model’s goal so that PESTPP–IES prioritized the fit between certain observations and simulated equivalent values during the calibration process. Weights were initially set using the error-based weighting scheme from Hill and Tiedeman (2007) based on approximated measurement error and variance in observations. Then, weights were adjusted by observation group to represent their relative contribution to Φ; similar balancing approaches have been applied in other models (Leaf and others, 2023; Traylor and others, 2023). Each observation group’s contribution to Φ was based on the importance of the observation group to better simulate groundwater levels in and around the cones of depression within the CGWAs, particularly for the final 9 years of the simulation period (January 1, 2010, through December 31, 2018). The Grand Prairie observation group weights were balanced to improve the fit in the MRVA and middle Claiborne aquifers (Traylor and Weisser, 2024).

The Grand Prairie model observations were also assigned a standard deviation that was used to create the noise observations during the automated calibration with PESTPP–IES, which creates noise observation data to properly assess the effect of measurement error on parameter uncertainty (White and others, 2020; Traylor and Weisser, 2024). For preliminary automated calibration runs, PESTPP–IES generated noise observations that included some unrealistic values such as negative streamflows or very low groundwater-level altitudes. These unrealistic values interfered with the approximation of parameter-to-observation gradients and limited the ability of PESTPP–IES to improve the fit to the observations and provide a robust uncertainty assessment. Therefore, standard deviations of observations were specified, which enabled PESTPP–IES to create realistic noise observations and produce an acceptable calibration result.

Groundwater Levels

The Cache model was calibrated to 2,534 stress period length head observations from 195 wells located throughout the model domain (fig. 7A; table 2). Observation groups for primary groundwater levels included head and priority_wells. The secondary (or derivative) observation groups were head_disp, head_sdiff, head_spinup, head_tdiff, and priority_wells_disp (table 2). Detailed descriptions of each observation group and the number of observations for each group are included in table 2. Although the middle Claiborne aquifer exists or subcrops in the Cache model domain, separate observation groups were not necessary for the MRVA and middle Claiborne aquifer as the head observation group was reflective of conditions in the MRVA and middle Claiborne aquifer because they are hydrologically connected in this area.

The Grand Prairie model was calibrated to 5,778 stress period length head observations from 340 wells located throughout the model domain (fig. 7B; table 3). The Grand Prairie model included more observation groups for primary groundwater levels than the Cache model. Additional groups of groundwater-level observations were needed to track the residuals in the MRVA and middle Claiborne aquifers separately because the middle Claiborne aquifer subcropped under much of the Grand Prairie model domain and had varying levels of hydrologic connection with the overlying MRVA. Also, a substantial number of observation wells were screened in the middle Claiborne aquifer owing to its importance as a municipal and agricultural water source (Kresse and others, 2014). Therefore, monthly measured groundwater-level observations were grouped by aquifer code (aqfr_cd; 110ALVM for MRVA and 124SPRT for middle Claiborne wells), provided in the National Water Information System database (U.S. Geological Survey, 2020), and by location in the priority Grand Prairie CGWA (fig. 1; table 3). These groups by aquifer and priority area also included the derivative observation groups with suffix “_disp” and “_tdiff.” Observation wells with other MRVA related aquifer codes such as 112TRRC and 112MRVAA were not included in the MRVA and middle Claiborne groups; however, they were included in the priority_wells group (table 3).

Calibration Parameters

Both models were parameterized using the same workflow as Leaf and others (2023), which included a combination of custom python scripts that utilized the pyEMU and Flopy python packages, following the methods in White and others (2020) (White and others, 2016; Bakker and others, 20239; Traylor and Weisser, 2024). A mix of fine- and coarse-scale spatial and temporal constants; pilot points; and zonal, gridded, and multiplier parameters were specified for inputs to provide PESTPP–IES adequate flexibility to adjust parameters and improve the spatial and transient simulation of the groundwater flow in each model during the calibration process. Many of the parameters in both models were the same; however, the Grand Prairie model included some additional parameters to aid in better simulating the groundwater system and stresses in the model domain.

Aquifer properties for both models were parameterized as “direct” values of Kh, Kvani, SY, and SS for zones corresponding to the resistivity facies classes from Minsley and others (2021), where the MRVA and non-MRVA facies were combined, and the hydrostratigraphic framework unit from Hart and others (2008). Further, Kh and Kvani properties were parameterized with pilot point multipliers by model layer (1–18 for Cache and 1–19 for Grand Prairie). SY pilot points were specified for model layers 1–10. Pilot points were placed every 10 cells for layers 1–15 and every 20 cells for layers 16–18 in both models. Initial values of “direct” parameters were assigned based on previous studies with groundwater-flow models in the regions (Broom and Lyford, 1981; Czarnecki and others, 2003; Clark and Hart, 2009) (fig. 8A, B). Pilot point multipliers of Kh, Kvani, and SY were assigned an initial value of 1.0 and adjusted during calibration.

Recharge from the SWB model was parameterized into spatial zone multipliers based on surface connectivity of the upper 15 m of AEM resistivity zones from Minsley and others (2021): pilot point multipliers with a spacing of 23 and 19 cells for the Cache and Grand Prairie models, respectively; and constant grid-scale multipliers for each stress period (fig. 8A, B). The Cache and Grand Prairie models included 18 and 48 surface connectivity zone multipliers, respectively; 41 and 116 pilot point multipliers, respectively; and 148 stress period multipliers each. The surface connectivity zones were used to parameterize recharge because it was assumed that the upper 15 m of the alluvial aquifer material affected the downward movement of water from the root zone to the water table. The surface connectivity dataset showed zones of lower and higher connection, where lower connection was indicated by the presence of more clay and higher connection was indicated by more sand (Minsley and others, 2021). The surface connectivity zones within the Cache model were further partitioned into zones that overlay MRVA and Crowley’s Ridge separately to avoid parameter compensation in two different physiographic regions (fig. 8A). Pilot points were specified to account for spatial heterogeneity in recharge that was not captured by the SWB model or the surface connectivity zones, which provided a finer-scale parameter to PESTPP–IES. Model grid-scale parameters for each stress period accounted for temporal trends in recharge not captured by the SWB model. Initial values for all recharge parameters were set to 1.0 for each multiplier; however, stress period multipliers were manually adjusted prior to automated calibration and discussed in the “Cache Model Calibration Process” and “Grand Prairie Model Calibration Process” sections. The lower and upper bounds for each group of recharge parameters were set to allow realistic flexibility during the automated calibration process so that recharge values would stay near or within the range of values from previous studies. Recharge parameter upper and lower bounds are listed in tables 4 and 5.

Cache and Grand Prairie model stream parameters included 601 and 1,496 SFR package rhk multipliers; 740 and 3,403 multipliers for SFR package inflows; 1,332 and 3,848 multipliers for runoff to the SFR package; and two and one conductance multipliers for RIV package cells, respectively. The rhk multipliers were specified for each NHDPlus version 2 stream segment unique identifier (Moore and Dewald, 2016) and surface connectivity zones from Minsley and others (2021). There is one stream segment unique identifier per hydrologic unit code level 12, which represented the fine-scale subwatershed parameterization. The surface connectivity zones (fig. 8A, B) were assumed to affect stream/aquifer interactions and acted as the coarse-scale parameters for the rhk (fig. 9A, B). Runoff generated by the SWB model and added to the SFR package network was parameterized as multipliers by surface connectivity zone and by stress period to account for coarse-scale uncertainty in runoff and temporal dynamics that may not have been entirely captured by the SWB model. The lower and upper bounds on the runoff multipliers were assigned values of 0.2 and 1.5 to allow for adjustment of parameters during automated calibration within reasonable limits. Inflows to the SFR package were parameterized with multipliers by stress period to account for the uncertainty in the estimated flows from the random-forest surface-water model and monthly average streamflows processed from daily measured values at streamgages near the model boundary. The Grand Prairie model stream parameters were specified like those from the Cache model, except that rhk values were also parameterized by the 11 coarse-scale hydrologic unit code level 8 watersheds to allow for more coarse-scale adjustment during automated calibration (fig. 9B). The conductance term for the Cache model RIV package was parameterized into two grid-scale multipliers for the cells representing the Black River and White River (fig. 9A). Similarly, the Grand Prairie model RIV package conductance was parameterized with a single grid-scale multiplier for the cells representing the Arkansas River (fig. 9B).

Conductance for the GHB package cells, where boundaries of the models intersected the Ozark Plateau aquifer system, were parameterized with grid-scale multipliers to account for uncertainty in the values obtained from Clark and others (2018). The Cache model included two GBH multipliers and a north and a south multiplier, and the Grand Prairie model included one GBH multiplier (fig. 9A, B). Lower and upper bounds for the conductance multipliers in both models were set to 0.01 and 100, respectively, to provide several orders of magnitude flexibility for PESTPP–IES during the automated calibration process.

Pumping rates in the WEL package were parameterized for both models by crop type, data source, possible reductions in pumping owing to surface-water access, and stress period. The Grand Prairie model included additional WEL package pumping rate parameters by combined model layer and AEM framework zone, and a separate parameter to adjust the fraction of MRVA pumping moved to the middle Claiborne aquifer to account for the uncertainty in agricultural pumping from the middle Claiborne aquifer in the Grand Prairie region, noted by McKee and Clark (2003). A python script moved a fraction of pumping from MRVA layers 5 through 8 to active cells designated with AEM framework zone two through seven in layer 15, which represented aquifer material in the middle Claiborne aquifer (Traylor and Weisser, 2024). Crop type pumping multipliers for both models included aquaculture, corn, cotton, nonagricultural, other, rice, and soybeans. The Cache model combined the unique stress periods with each crop type to create 889 crop types by stress period parameters, but the Grand Prairie model did not combine stress period multipliers with crop types. Grid scale pumping rate multipliers for both models were specified by SWUDS, Mississippi Embayment Regional Aquifer System version 2 (MERAS 2), and AIWUM data sources to provide flexibility and independence in parameters during the automated calibration process for pumping rates with different data sources. The Grand Prairie model combined data source with stress periods to create 412 “data source by stress period” multipliers; preliminary calibration runs indicated that the crop type by stress period multipliers scheme (used for the Cache model) was not adequate to capture the temporal trends in pumping rates.

Farms where irrigation reservoirs supplied stored surface water to irrigate crops during the summer were assumed to use less groundwater. Groundwater pumping cells in the WEL package that intersected potential irrigated cells supplied by surface water from irrigation reservoirs were assigned pumping rate multipliers to reduce pumping. The surface-water irrigated areas adjacent to irrigation reservoirs were unknown at the scale of both models; therefore, a dataset was constructed to approximate the crop area for which irrigation reservoirs supplied water. Using the irrigation reservoir dataset from Yaeger and others (2017), model cells that were within 200 m of a 10–30-acre reservoir, within 500 m of a 30–100-acre reservoir, and within 1,000 m of a reservoir greater than 100 acres were assigned a new variable boundname designation of “sw_irr” in the WEL package files, indicating they may potentially receive stored surface water; cells that did not intersect the irrigation reservoir area dataset were assigned the boundname “no_sw_irr.” These boundnames in each WEL package file were used to parameterize pumping in both models where “no_sw_irr” was assigned a multiplier of 1.0 and fixed during the calibration, and where “sw_irr” received a multiplier of 0.7 for the Cache model and 0.001 for the Grand Prairie model. The multipliers differed for both models owing to uncertainty in irrigation practices in both areas based on anecdotal evidence and the state of the MRVA in both areas. Reductions in groundwater pumping were assumed to preserve the aquifer and use surface water, if available, which was more common in the Grand Prairie area because the MRVA was more depleted in this region compared to the Cache area. Further, the climate was assumed to be a factor in using groundwater where surface water was available; groundwater was assumed to be used in drier than normal summer months.

Parameter-elevation Regressions on Independent Slopes Model (PRISM) data were used to define dry, average, and wet stress periods in the Grand Prairie model, where average was defined as the 40th to 60th percentile in monthly precipitation from 1999 through 2017 (PRISM Climate Group, 2022). The initial parameters were assigned using the monthly climate type from the PRISM data. For “sw_irr” pumping cells, if the June, July, or August stress periods were dry (less than 40th percentile of mean precipitation), they were assumed to use groundwater and their multipliers were set to 1.0; average and wet stress periods were assumed to use only surface water from the irrigation reservoirs and therefore were assigned a multiplier of 0.001, which effectively reduced pumping in the cell to near zero. Note that log normal parameters in PEST++ cannot be zero, but the difference between zero pumping and 0.001 of original pumping is effectively the same in the model. These multipliers were left as adjustable given the uncertainty in practices in the region but allowed for a more reasonable starting point for the automated calibration with PESTPP–IES (Traylor and Weisser, 2024).

Cache Model Calibration Process

Preliminary calibration runs for both models indicated that PESTPP–IES, which is a gradient-based algorithm like the traditional GLM, was substantially affected by the prior parameter ensemble of values that it generated (White, 2018; White and others, 2020). Therefore, some parameters were manually adjusted to improve the prior ensemble of parameter values for the automated PESTPP–IES calibration process. The Cache model calibration process included some manual adjustment of coarse parameters prior to automated calibration with PESTPP–IES. Manual adjustments were made to areal recharge and well pumping; specifically, Cache model recharge stress period multipliers were reduced from 1.0 to 0.2 to reduce recharge to within previously reported values (Reed, 2003; Clark and Hart, 2009; Reitz and others, 2017). Well pumping by crop-type stress period multipliers and well pumping by data source multipliers were increased from 1.0 to 1.2.

Grand Prairie Model Calibration Process

The Grand Prairie model calibration process underwent a series of preliminary calibration attempts focused on achieving an acceptable fit between groundwater-level observations and simulated equivalent values in the MRVA. Although these preliminary “MRVA-focused” automated calibration runs yield very good results within the MRVA, they yielded poor and unacceptable results in the middle Claiborne aquifer. The “MRVA-focused” calibration attempts weighted MRVA groundwater-level observations an order of magnitude higher than middle Claiborne aquifer observations and did not include the script that moved some fraction of MRVA groundwater pumping to the middle Claiborne aquifer. The implementation of the python script that moved some fraction of MRVA pumping to layer 15 of the model (middle Claiborne aquifer), coupled with balancing of the objective function such that middle Claiborne aquifer groundwater-level observation groups were weighted similar to MRVA observations, aimed to better align the numerical model with the conceptual model and reduce bias in simulated middle Claiborne aquifer groundwater levels while maintaining an acceptable fit to observations in the MRVA, which resulted in the “balanced” version of the Grand Prairie model. Additional preliminary calibration attempts with the “balanced” model did not yield acceptable fit to the observed groundwater levels in the MRVA when using the prior parameter values that were generated by PESTPP–IES for the “MRVA-focused” model. Finally, the best parameters of the base realization from iteration two for the “MRVA-focused” calibration that were coincident with “balanced” model parameters were used as the prior values for the final automated calibration with PESTPP–IES for the “balanced” model.

Groundwater-Flow Model Calibration Results

This section of the report presents calibration results for the Cache and Grand Prairie models that include ensemble (all realizations) and “base” realization results for each model’s fit to the measured or estimated observations and their final “posterior” parameter distributions and values. A model’s fit provides a metric for assessing a model’s ability to simulate past conditions; the objective function Φ is an indication of model fit. Model fit to observations, sometimes called history matching, can identify spatial and or temporal bias in the calibrated model based on the residuals (Leaf and others, 2023). Residuals are the difference between observations and their simulated equivalent values where a positive residual indicates a model’s undersimulation bias and a negative value indicates a model’s oversimulation bias (Traylor and others, 2023). Each model’s simulated groundwater levels were extracted using the OBS package and streamflows were extracted using the SFR observations option in the SFR package (Langevin and others, 2022).

The assessment of parameter fields is also an important aspect of calibration results, wherein the final parameter values obtained from the calibration process are compared to reasonable ranges of previously published values based on the current conceptual understanding of the hydrologic system. New insights gained from the calibration process about the conceptual model are discussed in subsequent sections of this report.

Ensemble Objective Function

The automated calibration step for the Cache model, using PESTPP–IES, was run through three iterations with a prior ensemble of 500 realizations. The automated calibration step for the Grand Prairie model, using PESTPP–IES, was run through four iterations with a prior ensemble of 800 realizations. Additional prior realizations were specified for the Grand Prairie model to facilitate exploration of parameter space based on a more uncertain conceptual model of prior parameters. The conceptual model focuses on uncertainty in agricultural pumping and anisotropy of Kv below the MRVA, particularly the Vicksburg-Jackson confining unit and middle Claiborne aquifer. The reduction in Φ across all iterations for the Cache and Grand Prairie models is shown in figure 10A an 10B. The Cache model reduction in Φ was about a 73-percent change between iteration zero and iteration two (table 6). The reduction in Φ was accompanied by a reduction in the standard deviation of Φ from values larger than the mean Φ to values more than 50 percent less than mean Φ, which was an indication of a reduction in model uncertainty (table 6). The Grand Prairie model reduction in Φ was 59.2 percent between iteration zero and iteration one (fig. 10B). The standard deviation in Φ between iteration zero and one decreased by from 79 percent of total Φ for iteration zero to about 24 percent of total Φ for iteration one (table 7). The ensemble of Φ for the Grand Prairie model excluded Φ values greater than 1x10⁵ for iteration one because values greater than that were determined to produce unacceptable fits to observations (fig. 10B).

Table 6.    

Summary statistics for the ensemble of objective functions from the PEST++ Iterative Ensemble Smoother calibration of the Cache model.

[--, not applicable]

Table 6.    Summary statistics for the ensemble of objective functions from the PEST++ Iterative Ensemble Smoother calibration of the Cache model.

Objective Function Statistic	Iteration
	0a	1	2b	3
Mean of ensemble objective function (Φ) value	48,820	15,178	13,190	13,058
Standard deviation ensemble objective function (Φ) value	81,640	6,531	5,317	5,303
Change in mean from previous iteration, percent	--	−68.9	−13.1	−1.0
Change in mean from iteration 0a to iteration twob, in percent	--	--	−73.0	--

^a

Prior ensemble.

^b

Posterior ensemble.

Table 7.    

Summary statistics for the ensemble of objective functions from the PEST++ Iterative Ensemble Smoother calibration of the Grand Prairie model.

[--, not applicable]

Table 7.    Summary statistics for the ensemble of objective functions from the PEST++ Iterative Ensemble Smoother calibration of the Grand Prairie model.

Objective Function Statistic	Iteration
	0a	1b	2	3	4
Mean of ensemble objective function (Φ) value	124,267	50,668	40,254	36,459	36,674
Standard deviation ensemble objective function (Φ) value	97,689	12,111	8,506	8,809	11,101
Change in mean from previous iteration, percent	--	−59.2	−20.6	−9.4	0.6
Change in mean from iteration 0a to iteration oneb, in percent	--	−59.2	--	--	--

^a

Prior model.

^b

Posterior model.

Calibrated Model Fit

The Cache model fit to observations was assessed using the base realization from iteration two of the PESTPP–IES calibration run. In the priority CGWA, absolute groundwater levels and spatial and temporal trends in groundwater levels were simulated with minimal bias and exhibited a good fit to measured observations (fig. 11A–F). The measured versus simulated one to one plots and histograms of residuals for each observation group are shown in figures 3.1–3.26, and the measured versus simulated changes in groundwater-level plots from the calibrated base realization of iteration two for the Cache model are shown in figures 4.1 and 4.2. The mean absolute residual for the priority_wells observation group was 1.58 m, which indicated a good fit to observed groundwater levels in the Cache CGWA (table 8). The fit to observed groundwater levels outside the Cache CGWA (head observation group) was also good, with a mean absolute residual of 2.41 m (table 8). Mean groundwater-level residuals for priority_wells and head observation groups were both positive, 0.33 and 1.12 m, respectively, which indicated a slight undersimulation bias (table 8). Spatial and temporal trends in simulated groundwater levels were minimal in the priority CGWA with some small over- and undersimulation, likely caused by uncertainty in agricultural water use from year to year. Most of the larger residuals in the head observation group were located outside of the CGWA, mainly in the north-central region of the model domain or east of Crowley’s Ridge (fig. 12A). Temporal trends in simulated groundwater levels during the primary calibration period from January 1, 2010, through December 1, 2018 (priority_wells_disp and head_disp observation groups) exhibited minimal bias with a mean residual of −0.01 m and 0.17 m, respectively (table 8; fig. 11A–F).

Table 8.    

Calibration statistics summary of residuals for the base realization of iteration two for the Cache model. [Left]

[Residual, observed value minus simulated value; Absolute residual, absolute value of observed value minus simulated value; m3/d, cubic meter per day; m, meter]

Table 8.    Calibration statistics summary of residuals for the base realization of iteration two for the Cache model. [Left]

Parameter groupa (units)	Count	Residuals
		Standard deviation	Minimum	25th percentile	Mean	75th percentile	Maximum
flux_estimated (m3/d)	352	309,461.21	−1,855,420.00	62,637.38	98,889.08	262,212.25	636,250.00
flux_estimated_sdiff (m3/d)	441	731,460.46	−4,218,250.00	−789,971.60	−547,164.94	−39,828.00	401,650.00
flux_estimated_tdiff (m3/d)	352	472,672.09	−1,753,450.00	−109,654.00	−42,882.88	24,002.90	3,323,539.00
flux_measured (m3/d)	254	44,296,870.85	−2,025,790.00	1,205,017.50	31,246,319.02	53,661,902.50	283,484,104.00
flux_measured_tdiff (m3/d)	252	30,273,557.86	−183,476,009.00	−3,049,036.25	131,068.79	1,642,733.75	211,659,304.00
head (m)	1333	3.54	−8.03	−0.7	1.12	2.48	19.95
head_disp (m)	654	1.32	−4.31	−0.64	0.17	0.83	5.73
head_sdiff (m)	68	2.35	−7.05	−0.96	0.35	1.85	5.18
head_spinup (m)	508	3.48	−28.34	−0.54	1.16	2.9	16.97
head_tdiff (m)	1545	1.23	−7.94	−0.43	0.10	0.67	8.59
priority_wells (m)	693	2.33	−12.65	−0.39	0.33	1.24	6.69
priority_wells_disp (m)	430	1.93	−11.89	−0.24	−0.01	0.64	5.16
priority_wells_tdiff (m)	661	0.94	−7.65	−0.28	−0.03	0.26	3.78

^a

Parameter groups summarized in table 2.

Table 8.    

Calibration statistics summary of residuals for the base realization of iteration two for the Cache model. [Right]

Table 8.    Calibration statistics summary of residuals for the base realization of iteration two for the Cache model. [Right]

Absolute residuals
Standard deviation	Minimum	25th percentile	Mean	75th percentile	Maximum
231,360.44	8,507.20	73,984.78	227,810.53	297,297.75	1,855,420.00
719,483.12	21.30	58,126.00	562,786.38	789,971.60	4,218,250.00
411,096.54	133.00	21,550.22	236,860.47	274,625.25	3,323,539.00
44,283,665.83	99.44	1,285,020.00	31,264,957.30	53,661,902.50	283,484,104.00
27,568,749.01	1,370.00	665,427.38	12,483,938.48	12,079,324.75	211,659,304.00
2.83	0	0.67	2.41	3.07	19.95
0.91	0.01	0.34	0.98	1.33	5.73
1.57	0.03	0.47	1.78	2.49	7.05
2.7	0	0.61	2.48	3.48	28.34
0.92	0	0.22	0.83	1.12	8.59
1.74	0	0.45	1.58	2.2	12.65
1.65	0	0.24	1.01	1.01	11.89
0.76	0	0.09	0.54	0.66	7.65

Table 9.    

Calibration statistics summary of residuals for the base realization of iteration one of the Grand Prairie model. [Left]

[Residual, observed value minus simulated value; Absolute residual, absolute value of observed value minus simulated value]

Table 9.    Calibration statistics summary of residuals for the base realization of iteration one of the Grand Prairie model. [Left]

Parameter groupa [units]	Count	Residual
		Standard deviation	Minimum	25th percentile	Mean	75th percentile	Maximum
flux_estimated (m3/d)	659	1,166,288.68	−13,217,000.00	−190,450.00	−202,138.72	170,814.00	3,787,800.00
flux_estimated_sdiff (m3/d)	146	5,998,759.30	−12,131,000.00	−409,800.00	2,804,665.42	4,408,250.00	39,597,000.00
flux_estimated_tdiff (m3/d)	659	1,319,580.32	−8,606,000.00	−303,790.35	−112,286.49	93,965.50	9,460,000.00
grandprairie_mrva_priority_wells (m)	799	3.34	−12.56	−3.63	−1.35	0.49	18.44
grandprairie_mrva_priority_wells_disp (m)	465	1.26	−4.63	−0.09	0.69	1.16	7.52
grandprairie_mrva_priority_wells_tdiff (m)	770	0.9	−11.26	−0.1	0.05	0.16	8.07
grandprairie_mrva_spinup_wells (m)	451	4.37	−12.48	−4.6	−1.52	0.84	23.82
grandprairie_mrva_wells (m)	2,960	4.97	−24.24	−1.4	0.78	3.07	29.27
grandprairie_mrva_wells_disp (m)	1,291	1.29	−8.15	−0.92	−0.38	0.28	6.08
grandprairie_mrva_wells_tdiff (m)	2,714	1.35	−13.73	−0.32	0.09	0.33	12.94
grandprairie_sparta_priority_wells (m)	697	8.13	−51.82	−13.59	−9.95	−5.51	25.26
grandprairie_sparta_priority_wells_disp (m)	309	5.56	−37.16	−1.79	−1.08	1.53	8.23
grandprairie_sparta_priority_wells_tdiff (m)	650	5.44	−40.2	−1.85	−0.22	1.69	31.76
grandprairie_sparta_wells (m)	285	9.79	−39.25	−17.17	−12.93	−8.25	18.43
grandprairie_sparta_wells_disp (m)	95	2.71	−6.92	−0.36	0.69	0.68	15.02
grandprairie_sparta_wells_tdiff (m)	251	3.78	−17.99	−0.85	−0.04	0.84	17.78
head_spinup (m)	323	14.56	−42.13	−11.06	−4.45	2.17	33.39
priority_wells(m)	263	3.93	-9.67	-2.49	0.4	3.31	9.5

^a

Parameter groups summarized in table 3.

Table 9.    

Calibration statistics summary of residuals for the base realization of iteration one of the Grand Prairie model. [Right]

Table 9.    Calibration statistics summary of residuals for the base realization of iteration one of the Grand Prairie model. [Right]

Absolute residual
Standard deviation	Minimum	25th percentile	Mean	75th percentile	Maximum
1,071,073.52	410.44	63,570.31	503,554.96	481,650.24	13,217,000.00
5,264,379.01	5,100.00	824,750.00	4,010,102.41	4,748,000.00	39,597,000.00
1,170,115.54	19.80	59,398.95	619,810.79	603,679.50	9,460,000.00
2.38	0	0.78	2.71	3.71	18.44
1.06	0	0.25	0.97	1.23	7.52
0.82	0	0.05	0.39	0.37	11.26
2.6	0.02	1.87	3.82	5.21	23.82
3.61	0	0.9	3.5	4.79	29.27
0.98	0	0.32	0.93	1.16	8.15
1.14	0	0.12	0.73	0.85	13.73
6.99	0.1	6.4	10.78	13.69	51.82
4.81	0.01	0.56	3	3.44	37.16
4.38	0	0.68	3.23	4.05	40.2
7.7	0.42	10.08	14.27	17.22	39.25
2.46	0.02	0.25	1.32	1.07	15.02
3.14	0	0.34	2.1	2.16	17.99
10.66	0	2.53	10.86	16.65	42.13
2.16	0.01	1.72	3.3	4.36	9.67

The primary streamflow observation groups of flux_measured and flux_estimated exhibited mean residuals of 31,246,319.02 cubic meters per day (m³/d) (12,771.45 ft³/s) and 98,889.08 m³/d (40.4 ft³/s) (table 8). Simulated streamflows at sites with estimated streamflow observations (flux_estimated observation group) exhibited a better fit compared to the measured data (fig. 13A–D). The large mean residual for the flux_measured was likely caused by the lack of measured streamflow data in the Cache model domain and the very low observation weights assigned based on the groundwater-level centric focus of the model. The only site with nonzero weighted and measured streamflows was the Cache River at Patterson, Ark. (USGS station 07077500) (fig. 7A), which meant that PESTPP–IES was not affected by the residuals for those observations and did not attempt to reduce it in the parameter upgrade process. Spatial differences in streamflow (flux_sdiff or flux_estimated_sdiff), an indication of the magnitude of gaining or losing flow along a stream reach, were oversimulated for the Cache River and Bayou De View, with a mean residual of −547,164.94 m³/d (223.65 ft³/s), which indicated the model allowed less leakage of stream water back into the aquifer (table 8). However, the mean residual was about 20 to 50 percent of mean streamflow and therefore the model is a reasonable representation of the natural system dynamics. The locations with estimated observations such as Bayou De View or the L’Anguille River tended to undersimulate streamflow particularly during low flow periods. This undersimulation indicated that stream inflows, estimated by the random forest surface-water model from Dietsch and others (2023), may have had a low bias during low flow periods or there were shallow processes, such as interflow or bank storage, not simulated by this groundwater-flow model that contributed water to streams during low flow periods (Leaf and others, 2023).

The “balanced” Grand Prairie model fit to observations was assessed using the base realization from iteration one of the PESTPP–IES calibration run. Iteration one was selected to avoid an overfit solution given that many of the prior parameter values for aquifer properties were derived from the “MRVA-focused” calibration run. Also, other studies have indicated that a realization chosen from iterations one through three is best for adequate history matching and reasonable parameter fields (White, 2018; Hunt and others, 2021; Leaf and others, 2023; Traylor and others, 2023). The measured/estimated compared with simulated one to one plots and histograms of residuals for each observation group from the calibrated base realization of iteration one for the “balanced” Grand Prairie model are shown in appendix figures 5.1–5.36.

The “balanced” Grand Prairie model adequately simulated streamflow across the model domain. The mean absolute residual for the flux_estimated observation group was −503,554.96 m³/d (−205.82 ft³/s), which is about 1 to 30 percent of estimated mean total flow for major streams including the White River, La Grue Bayou, Bayou Meto, and Wabbaseka Bayou (figs. 7B, 14A–F; table 9). Most streams exhibited a good fit between estimated and simulated streamflow during low flow periods. Spatially, streams with estimated observations in the southeastern region of the model in Arkansas County, such as La Grue Bayou (fig. 14A), exhibited some undersimulation of streamflow during low flow periods at multiple locations. The undersimulation of streamflow for La Grue Bayou and other streams in the southeastern region of the model may have been caused by uncertainty in diversions for irrigation reservoirs that led to more water being diverted in the model than what was actually diverted, or the streambed Kv values were not low enough to keep water from leaking through the streambed; however, this bias was not enough to affect the model’s ability to better understand La Grue Bayou’s stream-aquifer dynamics (fig. 14A).

Overall, the “balanced” Grand Prairie model adequately simulated absolute groundwater levels and spatial and temporal trends in groundwater levels in the MRVA, and oversimulated absolute groundwater levels in the middle Claiborne aquifer; however, spatial and temporal trends simulated in the middle Claiborne aquifer were accurate (fig. 15A–F). The mean absolute residuals for the grandprairie_mrva_priority_wells and grandprairie_sparta_priority_wells observation groups were 2.71 and 10.78 m, respectively, which indicated that the “balanced” model simulated groundwater levels in the MRVA well, but there was substantial bias in the middle Claiborne aquifer (table 9). The mean residual for the grandprairie_mrva_priority_wells and grandprairie_sparta_priority_wells observation groups were −1.35 m and −9.95 m, respectively, which indicated a small oversimulation bias in the MRVA and a large oversimulation bias in the middle Claiborne aquifer (table 9). Spatial trends in simulated groundwater levels for the CGWA of the MRVA exhibited minor bias for the grandprairie_mrva_priority_wells, whereas groundwater levels outside of the CGWA (grandprairie_mrva_wells_residuals observation group) exhibited some undersimulation bias (mean residual of about 4 m) in the region between Bayou Meto and the Arkansas River and some oversimulation bias (mean residual of about −5 m) east of the White River (fig. 12B). Temporal trends in simulated groundwater levels, particularly for the focus calibration period from January 1, 2010, through December 31, 2018 (grandprairie_mrva_priority_wells_disp observation group), within the CGWA of the MRVA exhibited minor under simulation bias with a mean residual of 0.69 m, and adequately simulated pumping signals and overall changes in groundwater levels across the focus simulation period (fig. 12B). Spatial trends in groundwater levels outside the CGWA of the MRVA also exhibited minimal bias with a mean residual of −0.38 m for the grandprairie_mrva_wells_disp observation group (table 9, fig. 15B). Additional measured compared with simulated plots of groundwater levels for the calibrated base realization of iteration one for the Grand Prairie model are included in appendix figures 5.1–5.36 and 6.1–6.2.

Spatial bias in the simulated middle Claiborne aquifer groundwater levels varied mostly from north to south across the Grand Prairie region. The “balanced” simulation exhibited the smallest oversimulation bias in Tier 1 and the central part of Tier 2 CGWAs with a mean residual in groundwater levels of about −4 m; the largest oversimulation bias was at the southern edges of Tier 2 and Tier 3 with mean residuals of about −11 m (fig. 15B). This increase in residuals farther south in the CGWA corresponded to the presence of the Vicksburg-Jackson confining unit, which subcrops between the MRVA and middle Claiborne aquifer. The spatial bias in simulated groundwater levels for the middle Claiborne aquifer indicated that the “balanced” Grand Prairie model better simulated the unconfined system where the MRVA and middle Claiborne aquifer are hydraulically well connected north of the Vicksburg-Jackson subcrop in southern Lonoke and Prairie counties compared to the confined middle Claiborne aquifer beneath the Vicksburg-Jackson confining unit in the central and southern region of the model domain in Jefferson and Arkansas counties (fig. 12B). Overall, the oversimulation and spatial bias in middle Claiborne aquifer groundwater levels can be attributed to uncertainty in agricultural pumping locations and volumes in the middle Claiborne aquifer across the Grand Prairie region.

Discussion of Calibrated Parameter Values

The Cache and Grand Prairie models were calibrated to produce reasonable parameter fields and values in addition to an acceptable fit to observations. Reasonable parameter fields should be relatively smooth with few extreme values except where conceptually appropriate (Hunt and others, 2019). The calibrated parameter fields for both models included aquifer properties of Kh, Kvani, SS, and SY; streambed Kv; areal recharge; pumping from wells; and general head and river conductance.

For the Cache model, the base realization of iteration two produced reasonable parameter fields and good fit to the observations (fit was discussed in the “Calibrated Model Fit” section of this report). Calibrated Kh for the MRVA model layers (1–8) ranged from 0.03 to 500 m/d (fig. 16). The largest MRVA Kh values (500 m/d) were in a region of higher resistivity in model layer 5. The smallest mean Kh for the MRVA was in a low resistivity region of layer 1 (18.86 m/d; fig. 16; table 10). The mean Kh for layers 1–4 and 5–8 was about 28 and 45 m/d, respectively, and similar to calibrated Kh values from previous studies (table 10; Reed, 2003; Clark and Hart, 2009; Rashid and others, 2015). The larger mean Kh values for deeper MRVA layers are consistent with the conceptual understanding that the coarser grained sediments are present deeper in the alluvial system (Reed, 2003). Model layers 9 through 15, which represent the Tertiary units, exhibited varying ranges of Kh. The Kh values for model layers 9–15 ranged from 0.13 to 500 m/d with mean Kh values from about 28 to 75 m/d (fig. 16). The largest Kh values were in the middle Claiborne aquifer and the smallest Kh values were in the Midway confining unit, which is consistent with the conceptual model of these units as a sand aquifer and clay dominated confining layer (fig. 16; Minsley and others, 2021). The Midway confining unit is most prominently visible in the AEM resistivity data from Minsley and others (2021; fig. 3A,D) and the Kh fields in figure 16 where it subcrops beneath the MRVA with an eastward dip for the entire Cache model. The shallowest subcrop of the Midway confining unit is in model layer 9 at the western end of the model (dips east) and is visible as the region of Kh less than 0.02 m/d (layer 9 in fig. 16). The middle Claiborne aquifer is the region of larger Kh to the east of the Midway confining unit and directly beneath the MRVA. The middle Claiborne aquifer Kh ranged from 1 to 500 m/d, similar to the MRVA; therefore, the MRVA and middle Claiborne aquifer are hydrologically connected in this region. The Kh values in the transition region between the Midway confining unit and the middle Claiborne aquifer represent the Wilcox aquifer system (lower and middle Wilcox aquifers); the transition between these two aquifers is difficult to distinguish in this region of the MERAS (Minsley and others, 2021). Model layers 16–18, which were defined using the framework layers from Hart and others (2008), represented the bottom of the lower Claiborne and Wilcox aquifers. The Kh ranges from about 0.02 to 100 m/d (fig. 16). Layer 16 was calibrated with the lowest mean Kh values of 0.2 m/d, which is consistent with interbedded clays in the lower Claiborne aquifer. Mean Kh for layers 17 and 18 was about 66 and 10 m/d for the Wilcox aquifer system, respectively. The pilot points multipliers for Kh by model layer increased by a mean of about 7 percent (mean multiplier of 1.068 for all model layers); locally, Kh pilot point values exhibited a range of 0.297–3.264. A summary of Kh pilot point multiplier values by model layer for the calibrated Cache model is provided in table 7.1. Precise values of Kh that cannot be distinguished from figure 16 are available in the model archive files associated with this report (Traylor and Weisser, 2024).

Table 10.    

Summary of calibrated aquifer properties for the base realization of iteration two of the Cache model.

[m/d, meter per day; Kh, horizontal hydraulic conductivity, Kvani; vertical anisotropy of hydraulic conductivity; Kv, vertical hydraulic conductivity; SY, specific yield]

Table 10.    Summary of calibrated aquifer properties for the base realization of iteration two of the Cache model.

Layer	Kha (m/d)	Kvanib (Kh/Kv)	SYc (unitless)
1	18.86	10.18	0.23
2	25.24	4.7	0.25
3	33.79	3.65	0.24
4	35.65	9.61	0.29
5	51.58	5.19	0.28
6	29.87	5.17	0.29
7	52.62	7.28	0.28
8	47.19	6.93	0.26
9	28.56	5.76	0.27
10	29.88	6.34	0.26
11	48.51	5.11	0.27
12	51.04	7.12	0.24
13	43.14	4.79	0.25
14	34.63	5.58	0.23
15	75.37	4.65	0.27
16	0.02	9.88	0.2
17	66.18	16.94	0.14
18	9.94	6.58	0.14

^a

Mean for layers 1–4: 28.39; layers 5–8: 45.32; layers 9–18: 38.73.

^b

Mean for layers 1–8: 6.59; layers 9–18: 7.28.

^c

Mean for layers 1–9: 0.27; layers 9–18: 0.23.

Table 11.    

Summary of calibrated Streamflow Routing package and River package properties for the base realization of iteration two of the Cache model.

[rhk, streambed hydraulic conductivity used in Streamflow Routing package; SY, specific yield; SFR, Streamflow Routing package; RIV, River package]

Table 11.    Summary of calibrated Streamflow Routing package and River package properties for the base realization of iteration two of the Cache model.

Stream	Minimum	Mean	Maximum
Vertical conductivity (rhk) used in SFR package, in meters per day
Village Creek	0.02	0.25	1.07
Cache River	0.02	0.3	2.89
Bayou De View	0.01	0.14	0.7
L'Anguille River	0.03	0.18	0.84
Riverbed conductance used in RIV package, in square meters per day
Black River	27,686	70,324	122,702
White River	28,678	93,338	271,915

Cache model Kvani ranged from about 1 to 33 (dimensionless) across all model layers. Mean Kvani for the MRVA and deeper Tertiary layers was about 6.59 and 7.28, respectively (table 10). Kvani ratios less than 10 indicated that vertical groundwater flow across model layers was not substantially impeded in the MRVA or the deeper Tertiary units (Fetter, 2001). The primary limiting factor in movement of water across all model layers was the Kh values. The largest Kvani value in the MRVA was about 32, located in model layer 4 (fig. 17). The MRVA Kvani values are consistent with coarse-grained alluvial sediments that facilitate vertical flow in a highly pumped system such as the MRVA in the Cache region. The largest Kvani value in the Tertiary layers 9–15 was about 33, located in one location in model layer 10 (fig. 17). Layer 17 exhibited a larger Kvani (about 17) than layers 16 or 18 but was still in the range where vertical flow between the layers was not substantially impeded. The Kvani values for the fine-grained Midway confining unit are smaller than previously published values in Clark and Hart (2009); however, the values in this study did not negatively affect the ability of the Cache model to accurately simulate past conditions, and parameter compensation was not evident in the other parameters based on reasonable parameter values and parameter fields for this model. The pilot point multipliers for Kvani by model layer increased by a mean of about 7 percent (mean multiplier of 1.066 for all model layers); locally, Kvani pilot point values exhibited a range of 0.275–3.725. A summary of Kvani pilot point multiplier values by model layer for the calibrated Cache model is available in table 7.1. Precise values of Kvani that cannot be distinguished from figure 17 are available in the model archive files associated with this report (Traylor and Weisser, 2024).

Cache model SY ranged from 0.08 to 0.35 (dimensionless) and exhibited a mean SY of about 0.27 in the MRVA, which were consistent with previously published values (Reed, 2003; Clark and Hart, 2009; Rashid and others, 2015). The Cache model in this report provided a higher resolution of storage and aquifer properties than all previous models that simulated the Cache region whereby those previous studies specified single constant values or large zones for the MRVA. Calibrated SY fields in the heavily pumped layers of the MRVA (layers 4–7) exhibited multiple regions of larger SY values of 0.3 to 0.35 (fig. 18A; table 10). The Tertiary units (layers 9–18) exhibited a mean SY of 0.23. The SS values for all model layers ranged from 5.0x10⁻⁷ to 1x10⁻⁶ 1/m and were consistent with values from Reed (2003) (fig. 18B). The pilot point multipliers for SY by model layer increased from prior values by a mean of about 9 percent (mean multiplier of 1.079 for all model layers); locally, SY pilot point values exhibited a range of 0.337–4.152 (fig. 18). A summary of SY pilot point multiplier values by model layer for the calibrated Cache model is provided in table 7.1. Precise values of SY that cannot be distinguished from figure 18 are available in the model archive files associated with this report (Traylor and Weisser, 2024).

Areal recharge multipliers assigned by surface connectivity zone in the Cache model (fig. 8A) exhibited minimal to moderate deviations from their prior values of 1.0. The calibrated multipliers ranged from 0.57 to 1.49 and a mean of 1.08. The mean value by physiographic region (MRVA, Crowley’s Ridge, and Ozark Plateau) of the Cache model domain was 1.01 for the MRVA, 1.16 for Crowley’s Ridge, and 1.07 for the Ozark Plateau. The calibrated multiplier of 1.01 for the MRVA indicated that recharge from the SWB model deviated less from the spatial trends of the surface connectivity in the MRVA region than Crowley’s Ridge, likely because there were more observations and more focus on recharge in the MRVA, but the Crowley’s Ridge recharge also did not deviate substantially from the SWB model spatial trends (fig. 19A). Temporally, the stress period multipliers ranged from 0.07 to 0.87 (fig. 19B) with a mean of 0.22 and a standard deviation of 0.11; the prior value was 0.2. This 13- to 97-percent reduction in the SWB model estimated recharge (water leaving the root zone) indicated that the unsaturated zone substantially attenuated the downward movement of water toward the water table and multipliers of this magnitude were necessary in the absence of an unsaturated zone in the model simulation (fig. 19B).

Streambed Kv in the SFR package for the Cache model ranged from 0.01 to 2.89 m/d with a mean value of 0.25 m/d, similar to those in Reed (2003), and also varied substantially by reach (fig. 19A). Based on streambed Kv fields, streams in the Cache model were well connected to the MRVA in a few areas of short reaches, with most other reaches exhibiting small streambed Kv values that limited the stream’s connection with the MRVA. Mean streambed Kv was largest for the Cache River at 0.3 m/d in comparison among Bayou De View, the L’Anguille River, and Village Creek. The largest streambed Kv values (2.45 and 2.89 m/d) were on the northern reaches of the Cache River; Bushy Creek Ditch, a tributary to the L’Anguille River also had several reaches with streambed Kv of 2.16 m/d (fig. 19A). The smallest streambed Kv values were 0.01 m/d located on a northern reach of Bayou De View (fig. 19A). Riverbed conductance for the Black River and White River was 27,686 to 271,915 m²/d. The Black River was less connected to the MRVA with a mean conductance of 70,324 m²/d, and the White River exhibited a higher degree of connection to the MRVA with a mean conductance of 93,338 m²/d (fig. 19A, table 11). The upper reaches of the White River, below the confluence with the Black River and the White River near the outflow of the model domain, exhibited the largest riverbed conductance values. Like the streambed Kv in the SFR package, these riverbed conductances varied spatially from reach to reach, with large conductance reaches often adjacent to smaller conductance reaches (fig. 19A). The least conductive reaches were on the Black River in the northwest part of the model (fig. 19A).

Calibrated well pumping multipliers by crop type and stress period (wel_croptype_mult), and by data source (wel_datasource_mult) did not deviate substantially from their prior values of 1.2. These minimal grid scale adjustments to pumping stresses during the calibration process indicated that each data source provided pumping rates across the entire simulation period that were conceptually accurate after increasing the prior multipliers to 1.2. The calibrated combined crop type and stress period multipliers had a mean value of 1.21 to 1.25 for each of the crop types across all 148 stress periods (table 12). The multipliers by data source (SWUDS, MERAS 2, and Irrigation Water Use Model [IWUM]) exhibited a range from 1.0 to 1.33 (Traylor and Weisser, 2024). The largest deviation from the prior was in agricultural pumping (IWUM) multiplier, which was reduced by 0.2 to 1.0 (Traylor and Weisser, 2024). The multiplier for the spin-up stress periods (stress periods 2–7), which adjusted pumping from the MERAS 2 data source, increased by 0.13 to 1.33. The municipal pumping from the SWUDS data source decreased to 1.17. The calibrated water use multipliers by crop type for each stress period in the Cache model are shown in figure 8.1. Precise values of SY that cannot be distinguished from figure 18 are available in the model archive files associated with this report (Traylor and Weisser, 2024).

Table 12.    

Summary of the calibrated parameters by parameter group for the base realization of iteration two for the Cache model.

[--, unitless; posterior, post calibration]

Table 12.    Summary of the calibrated parameters by parameter group for the base realization of iteration two for the Cache model.

Calibration parameter group name	Calibration parameter group description	Units	Mean posterior value	Lower posterior value	Upper posterior value
k_zonea	Horizontal hydraulic conductivity by resistivity zone	Meters per day	41.38	0.016	256.19
k_pp_layer0	Horizontal hydraulic conductivity for pilot points for model layer 1	--	1.00	0.382	2.83
k_pp_layer1	Horizontal hydraulic conductivity for pilot points for model layer 2	--	1.11	0.382	2.67
k_pp_layer2	Horizontal hydraulic conductivity for pilot points for model layer 3	--	1.10	0.416	3.26
k_pp_layer3	Horizontal hydraulic conductivity for pilot points for model layer 4	--	1.05	0.457	2.51
k_pp_layer4	Horizontal hydraulic conductivity for pilot points for model layer 5	--	1.10	0.397	3.18
k_pp_layer5	Horizontal hydraulic conductivity for pilot points for model layer 6	--	1.06	0.329	2.22
k_pp_layer6	Horizontal hydraulic conductivity for pilot points for model layer 7	--	1.08	0.387	2.25
k_pp_layer7	Horizontal hydraulic conductivity for pilot points for model layer 8	--	1.07	0.374	2.55
k_pp_layer8	Horizontal hydraulic conductivity for pilot points for model layer 9	--	1.01	0.430	2.34
k_pp_layer9	Horizontal hydraulic conductivity for pilot points for model layer 10	--	1.04	0.297	2.34
k_pp_layer10	Horizontal hydraulic conductivity for pilot points for model layer 11	--	1.09	0.376	2.64
k_pp_layer11	Horizontal hydraulic conductivity for pilot points for model layer 12	--	1.08	0.318	2.42
k_pp_layer12	Horizontal hydraulic conductivity for pilot points for model layer 13	--	1.04	0.406	2.23
k_pp_layer13	Horizontal hydraulic conductivity for pilot points for model layer 14	--	1.04	0.415	3.12
k_pp_layer14	Horizontal hydraulic conductivity for pilot points for model layer 15	--	1.13	0.438	2.75
k_pp_layer15	Horizontal hydraulic conductivity for pilot points for model layer 16	--	1.07	0.369	1.91
k_pp_layer16	Horizontal hydraulic conductivity for pilot points for model layer 17	--	1.04	0.457	2.16
k_pp_layer17	Horizontal hydraulic conductivity for pilot points for model layer 18	--	1.09	0.479	1.98
kvani_zonea	Vertical anisotropy by resistivity zone	Meters per day	7.01	1.450	24.92
kvani_pp_layer0	Vertical anisotropy for pilot points for model layer 1	--	1.10	0.441	2.41
kvani_pp_layer1	Vertical anisotropy for pilot points for model layer 2	--	1.06	0.364	2.64
kvani_pp_layer2	Vertical anisotropy for pilot points for model layer 3	--	1.08	0.411	3.73
kvani_pp_layer3	Vertical anisotropy for pilot points for model layer 4	--	1.07	0.507	2.31
kvani_pp_layer4	Vertical anisotropy for pilot points for model layer 5	--	1.05	0.403	2.91
kvani_pp_layer5	Vertical anisotropy for pilot points for model layer 6	--	1.04	0.409	3.03
kvani_pp_layer6	Vertical anisotropy for pilot points for model layer 7	--	1.08	0.399	2.22
kvani_pp_layer7	Vertical anisotropy for pilot points for model layer 8	--	1.09	0.275	2.53
kvani_pp_layer8	Vertical anisotropy for pilot points for model layer 9	--	1.06	0.360	3.10
kvani_pp_layer9	Vertical anisotropy for pilot points for model layer 10	--	1.05	0.378	3.33
kvani_pp_layer10	Vertical anisotropy for pilot points for model layer 11	--	1.08	0.393	3.17
kvani_pp_layer11	Vertical anisotropy for pilot points for model layer 12	--	1.08	0.411	2.26
kvani_pp_layer12	Vertical anisotropy for pilot points for model layer 13	--	1.07	0.288	2.29
kvani_pp_layer13	Vertical anisotropy for pilot points for model layer 14	--	1.03	0.360	2.62
kvani_pp_layer14	Vertical anisotropy for pilot points for model layer 15	--	1.04	0.313	2.22
kvani_pp_layer15	Vertical anisotropy for pilot points for model layer 16	--	1.09	0.378	3.36
kvani_pp_layer16	Vertical anisotropy for pilot points for model layer 17	--	1.09	0.316	2.21
kvani_pp_layer17	Vertical anisotropy for pilot points for model layer 18	--	1.04	0.393	2.32
ss_zonea	specific storage by resistivity zone	1 per meter	0.000001	0.000	0.00
sy_zonea	Specific yield by resistivity zone	--	0.24	0.118	0.34
sy_pp_layer0	Specific yield for pilot points for model layer 1	--	1.12	0.396	2.71
sy_pp_layer1	Specific yield for pilot points for model layer 2	--	1.04	0.379	2.10
sy_pp_layer2	Specific yield for pilot points for model layer 3	--	1.06	0.355	2.88
sy_pp_layer3	Specific yield for pilot points for model layer 4	--	1.07	0.435	2.64
sy_pp_layer4	Specific yield for pilot points for model layer 5	--	1.07	0.379	2.82
sy_pp_layer5	Specific yield for pilot points for model layer 6	--	1.07	0.339	2.57
sy_pp_layer6	Specific yield for pilot points for model layer 7	--	1.09	0.337	4.15
sy_pp_layer7	Specific yield for pilot points for model layer 8	--	1.08	0.368	2.44
sy_pp_layer8	Specific yield for pilot points for model layer 9	--	1.08	0.452	2.45
sy_pp_layer9	Specific yield for pilot points for model layer 10	--	1.11	0.434	3.75
sy_pp_layer10	Specific yield for pilot points for model layer 11	--	1.09	0.365	2.82
riv_cond_mult	River package conductance multiplier	--	1.42	1.346	1.49
ghb_cond_mult	General head boundary multiplier for the north Ozark and south Ozark interface	--	1.31	1.053	1.56
wel_datasource_mult	WEL package multiplier by data source	--	1.17	0.998	1.33
wel_swirr_mult	WEL package multiplier by surface-water irrigated area	--	0.86	0.573	1.00
rzonemulta	Recharge multipliers by shallow resistivity zone	--	1.08	0.570	1.49
rch_pp_mult	Recharge pilot point multipliers	--	1.04	0.613	1.83
rchspmult	Stress period recharge multipliers	--	0.22	0.068	0.87
wel_croptype_mult	WEL package multiplier by crop type	--	1.22	0.629	2.00
sfr_inflow_mult	Streamflow Routing package inflow multiplier	--	1.00	0.902	1.08
sfr_runoff_mult	Streamflow Routing package runoff multiplier	--	1.01	0.660	1.50
sfr_kv_mult	Streamflow Routing package vertical streambed hydraulic conductivity multiplier	--	0.27	0.008	2.36

^a

Mean value for multiple zones.

For the Grand Prairie model, the base realization of iteration one produced acceptable parameter fields and fit to the observations (fit was discussed in the “Calibrated Model Fit” section of this report). Calibrated Kh for the MRVA model layers (1–8) ranged from 1.0 to 300 m/d (fig. 20). Each MRVA layer exhibited areas with Kh values of 300 m/d, with layer 3 exhibiting the most widespread areas with 300 m/d Kh values (fig. 20; Traylor and Weisser, 2024). For all MRVA layers, the areas with Kh values of 300 m/d were correlated to areas of higher resistivity from the AEM data (fig. 20). The mean Kh for layers 1–4 and 4–8 was about 78 and 48 m/d, respectively, and within the range of calibrated Kh values from previous studies (table 13; Reed, 2003; Clark and Hart, 2009). The larger mean Kh values for shallower MRVA layers are consistent with the AEM resistivity facies zones, which exhibited higher resistivity zones in the shallower layers but is not consistent with the overall concept that coarser grained sediments are present deeper in the alluvial system that was in the Cache model (Reed, 2003). Model layers 9 through 15, which represented the Tertiary-age units, including the Vicksburg-Jackson confining unit, middle Claiborne aquifer, and the lower Claiborne confining unit, exhibited varying ranges of Kh based on their representative geologic unit. The Kh values for model layers 9–15 ranged from 1.0 to 263 m/d with a mean value of 7.9 m/d (fig. 20; table 13). The largest Kh values were in small areas in layers 9 and 15 that corresponded to high resistivity areas (Minsley and others, 2021). The Vicksburg-Jackson confining unit was most prevalent in the southern half of the Grand Prairie model and was prominent in the Kh fields for layers 9–13 where it is represented by smaller Kh values (fig. 20). The middle Claiborne aquifer is unconfined and in hydraulic connection with the MRVA in the northern half of the model domain where the Vicksburg-Jackson confining unit is absent. In the southern half of the model domain, where the Vicksburg-Jackson confining unit is present, the middle Claiborne aquifer is confined and exhibited minimal connection with the overlaying MRVA. Model layers 16–18, which were defined using the framework layers from Hart and others (2008), represented the bottom surface of the lower Claiborne and Wilcox aquifers, like in the Cache model. In Grand Prairie model layers 16–19, Kh ranges from about 1.0 to 8.8 m/d and mean value of 2.1 m/d, which were like values from McKee and Clark (2003) and Clark and Hart (2009) (fig. 20). Layer 16 calibrated with the lowest mean Kh values of 0.2 m/d, which is consistent with interbedded clays in the lower Claiborne aquifer. On average, the pilot point multipliers for Kh by model layer increased by about 16 percent (mean multiplier of 1.16 for all model layers); locally, Kh pilot point values exhibited a range of 0.171–5.567 (table 7.2). A summary table of Kh pilot point multiplier values by model layer for the calibrated Grand Prairie model is in table 7.2. Precise values of SY that cannot be distinguished from figure 20 are available in the model archive files associated with this report (Traylor and Weisser, 2024).

Table 13.    

Summary of calibrated parameters by parameter group for the base realization of iteration one of the Grand Prairie model.

[m/d, meter per day; Kh, horizontal hydraulic conductivity, Kvani, vertical anisotropy of hydraulic conductivity, calculated as Kh/Kv; Kv, vertical hydraulic conductivity]

Table 13.    Summary of calibrated parameters by parameter group for the base realization of iteration one of the Grand Prairie model.

Layer	Kha (m/d)	Kvanib (Kh/Kv)	SYc (unitless)	SSd (1 per meter)
1	74.76	168.46	0.27	1.48E-06
2	60.87	236.72	0.30	1.43E-06
3	129.64	198.82	0.30	1.00E-06
4	45.02	82.39	0.30	1.80E-06
5	69.14	104.99	0.30	1.10E-06
6	44.04	158.41	0.28	1.05E-06
7	38.68	675.77	0.28	1.18E-06
8	38.17	89.50	0.26	9.10E-07
9	7.35	117.71	0.30	6.31E-07
10	10.55	533.39	0.27	8.44E-07
11	12.11	99.88	0.26	9.56E-07
12	8.03	191.36	0.29	1.14E-06
13	11.87	126.60	0.32	1.47E-06
14	3.62	264.44	0.30	2.00E-06
15	2.09	265.04	0.30	8.05E-07
16	1.67	600.14	0.13	6.11E-07
17	2.18	413.61	0.16	1.80E-06
18	3.4	359.67	0.21	3.94E-06
19	1.16	86.76	0.12	5.62E-07

^a

Mean for layers 1–4: 77.57; layers 1–8: 62.54; layers 9–15: 7.95.

^b

Mean for layers 1–8: 214.38; layers 9–19: 297.18.

^c

Mean for layers 1–8: 0.29; layers 9–19: 0.25.

^d

Mean for layers 1–8: 1.24E-06; layers 9–19: 1.34E-06.

Table 14.    

Summary of calibrated Streamflow Routing and River package parameters for the base realization of iteration one of the Grand Prairie model.

[rhk, streambed hydraulic conductivity used in Streamflow Routing package; SY, specific yield; SFR, Streamflow Routing package; RIV, River package]

Table 14.    Summary of calibrated Streamflow Routing and River package parameters for the base realization of iteration one of the Grand Prairie model.

Stream	Minimum	Mean	Maximum
Vertical conductivity (rhk) used in SFR package, in meters per day
Bayou Meto	0	6.80E-03	3
Bayou Two Prairie	0	1.03E-07	1.60E-06
Cache River	0	2.36E-07	6.50E-06
White River	0	1.8E-07	6.50E-06
Wattensaw Bayou	0	2.53E-07	4.72E-06
Riverbed conductance used in RIV package, in square meters per day
Arkansas River	64,890	237,368	557,150

Grand Prairie model Kvani (fig. 21) ranged from about 43 to 5,127 (dimensionless) across all model layers and overall was substantially higher than in the Cache model (fig. 17). Mean Kvani for the MRVA (layers 1–8) and deeper Tertiary layers (layers 9–19) was about 214 and 297, respectively (table 13). Kvani ratios greater than 100 indicated that vertical groundwater flow across model layers was substantially impeded (Fetter, 2001). The primary limiting factor in the vertical movement of water across all model layers was the aquifer characteristics indicated by the anisotropy ratios. The largest Kvani value in the MRVA was about 5,127, located in lower AEM resistivity zones of model layer 7 (figs. 4B, 8B, 21; Minsley and others, 2021). Large anisotropy values were also present in shallow MRVA layers, which is consistent with a shallow clay lens in the Grand Prairie model domain (figs. 4B, 20) that limits areal recharge. The largest Kvani value in the Tertiary layers 9–15 was about 4,654 and located in low resistivity zones of model layer 10 (Minsley and others, 2021) and spatially consistent with the Vicksburg-Jackson confining unit (approximately layers 9–12; figs. 4B, 21; table 15). Layers 16–19 exhibited a mean Kvani from 87 to 600 and exhibited some heterogeneity (fig. 21). The pilot point multipliers for Kvani by model layer increased by a mean of about 39 percent (mean multiplier of 1.39 for all model layers); locally, Kvani pilot point values exhibited a range of 0.140–9.365. A summary of Kvani pilot point multiplier values by model layer for the calibrated Grand Prairie model is provided in table 7.2.

Grand Prairie model SY ranged from 0.05 to 0.35 (dimensionless) and exhibited a mean SY of about 0.29 in the MRVA (layers 1–8), which were consistent with previously published values (fig. 22A; Reed, 2003; Clark and Hart, 2009). The Grand Prairie model, like the Cache model, provided a higher resolution of storage and aquifer properties for the MRVA than all previous models that simulated the Grand Prairie region with single constant values or large zones. Calibrated SY fields in all MRVA layers exhibited multiple regions of larger SY values of 0.05 to 0.35, which indicated the preference of the calibration toward higher storage values (fig. 22A). The Tertiary units (layers 9–18) exhibited a mean SY of 0.25, which is about 25 percent higher than in the middle Claiborne aquifer focused model from McKee and Clark (2003). The SS values for all model layers (fig. 22B) ranged from 5.6x10⁻⁷ to 3.9x10⁻⁶ 1/m and were consistent with values from Reed (2003) and McKee and Clark (2003) (table 13). The pilot point multipliers for SY by model layer increased by a mean of about 14 percent (mean multiplier of 1.138 for model layers 1–11); locally, SY pilot point multiplier values for the MRVA exhibited a range of 1.108–3.594. A summary of SY pilot point multiplier values by model layer for the calibrated Grand Prairie model is provided in table 7.2.

Areal recharge multipliers by surface connectivity zone exhibited minimal deviations from their prior values of 1.0. The calibrated multipliers by surface connectivity zone ranged from 0.78 to 1.17 and a mean of 1.02 (table 15). The calibrated multiplier of 1.02 for the MRVA indicated that recharge from the SWB model correlated well with the spatial trends of the surface connectivity (figs. 8B, 23A). Temporally, the stress period multipliers ranged from 0.10 to 1.19 with a mean of 0.39 and a standard deviation of 0.24 (fig. 23A). The mean deviation from prior values was −0.002, which indicated that overall, the calibration did not substantially alter the total amount of recharge to better match temporal trends across the simulation period. However, like the Cache model, the range in stress period multipliers (0.1 to 1.19) of the SWB model estimated recharge (water leaving the root zone) indicated that the unsaturated zone substantially attenuated the downward movement of water toward the water table for individual stress periods and locations (fig. 23A,B).

The rhk values in the SFR package were very low and ranged from 0 to 3.9x10⁻⁶ m/d (fig. 23A) with a mean value of 1.2x10⁻⁶ m/d, which were substantially lower than the stream rhk values derived from conductances in Clark and Hart (2009). There was some variance in rhk by reach (fig. 23A). Based on streambed Kv fields, streams in the Grand Prairie model exhibited limited connection to the MRVA in most areas with a few reaches that exhibited larger rhk values and a greater connection to the MRVA. Mean rhk was largest for Bayou Meto at 0.0068 m/d compared to Bayou Two Prairie, the Cache River, White River, and Wattensaw Bayou (fig. 23B; table 14). The largest rhk values (3 m/d) were Bayou Meto (fig. 23B; table 14). Riverbed conductance for the Arkansas River was 64,890 to 557,150 m²/d with a mean value of 237,368 m²/d, which equates to an rhk of 0.95 m/d (3.11 ft/d) (Langevin and others, 2017); similar to the Arkansas River values from Clark and Hart (2009); however, the model from Clark and Hart (2009) simulated fewer streams than the Grand Prairie model in this study. Like the rhk in the SFR package, these riverbed conductances varied spatially from reach to reach, with large conductance reaches often adjacent to smaller conductance reaches (fig. 23A).

Well pumping multipliers by crop type, data source, and stress period did not deviate substantially from their prior values and base multipliers of 1.0 (fig. 8.1). The calibrated crop type multipliers exhibited a range from 0.73 to 1.42 with a mean of 1.03; all crop type multipliers decreased between 0 and 10 percent from their prior values. Like the Cache model data source multipliers, these small grid scale adjustments to pumping stresses during the calibration process indicated that each data source provided pumping rates across the entire simulation period that were conceptually accurate and only minor fine-tuning of the pumping multipliers was needed during the calibration process. The combined data source (agricultural, MERAS 2, SWUDS, and thermoelectric) and stress period multipliers had a mean value of 1.03 for all 148 stress periods (table 15). The largest variance (standard deviation) was for the MERAS 2 pumping values and was in accord with the greater uncertainty in pumping during the first seven stress periods where irrigation practices were less known. The mean of the multiplier for the municipal pumping from the SWUDS data source increased from 1.002 to 1.003.

Table 15.    

Summary of the calibrated parameters by parameter group for the base realization of iteration one of the Grand Prairie model.

[--, unitless; posterior, post calibration]

Table 15.    Summary of the calibrated parameters by parameter group for the base realization of iteration one of the Grand Prairie model.

Calibration parameter group name	Calibration parameter group description	Units	Mean posterior value	Lower posterior value	Upper posterior value
k_zone	Horizontal hydraulic conductivity by resistivity zone	Meters per day	42.41	0.88	300.00
k_pp_layer0	Horizontal hydraulic conductivity for pilot points for model layer 1	--	1.14	0.26	3.88
k_pp_layer1	Horizontal hydraulic conductivity for pilot points for model layer 2	--	1.18	0.22	4.76
k_pp_layer2	Horizontal hydraulic conductivity for pilot points for model layer 3	--	1.13	0.21	3.35
k_pp_layer3	Horizontal hydraulic conductivity for pilot points for model layer 4	--	1.14	0.20	4.09
k_pp_layer4	Horizontal hydraulic conductivity for pilot points for model layer 5	--	1.16	0.20	4.71
k_pp_layer5	Horizontal hydraulic conductivity for pilot points for model layer 6	--	1.16	0.21	4.73
k_pp_layer6	Horizontal hydraulic conductivity for pilot points for model layer 7	--	1.07	0.24	4.00
k_pp_layer7	Horizontal hydraulic conductivity for pilot points for model layer 8	--	1.15	0.14	5.43
k_pp_layer8	Horizontal hydraulic conductivity for pilot points for model layer 9	--	1.15	0.19	4.77
k_pp_layer9	Horizontal hydraulic conductivity for pilot points for model layer 10	--	1.19	0.17	4.20
k_pp_layer10	Horizontal hydraulic conductivity for pilot points for model layer 11	--	1.14	0.25	4.27
k_pp_layer11	Horizontal hydraulic conductivity for pilot points for model layer 12	--	1.13	0.27	5.95
k_pp_layer12	Horizontal hydraulic conductivity for pilot points for model layer 13	--	1.17	0.22	3.55
k_pp_layer13	Horizontal hydraulic conductivity for pilot points for model layer 14	--	1.15	0.25	5.18
k_pp_layer14	Horizontal hydraulic conductivity for pilot points for model layer 15	--	1.12	0.24	5.57
k_pp_layer15	horizontal hydraulic conductivity for pilot points for model layer 16	--	1.24	0.32	4.78
k_pp_layer16	Horizontal hydraulic conductivity for pilot points for model layer 17	--	1.15	0.19	3.12
k_pp_layer17	Horizontal hydraulic conductivity for pilot points for model layer 18	--	1.23	0.23	3.49
k_pp_layer18	Horizontal hydraulic conductivity for pilot points for model layer 19	--	1.16	0.20	4.83
kvani_zone	Vertical anisotropy by resistivity zone	Meters per day	196.70	6.67	1100.53
kvani_pp_layer0	Vertical anisotropy for pilot points for model layer 1	--	1.15	0.16	4.66
kvani_pp_layer1	Vertical anisotropy for pilot points for model layer 2	--	1.09	0.22	3.56
kvani_pp_layer2	Vertical anisotropy for pilot points for model layer 3	--	1.12	0.19	5.76
kvani_pp_layer3	Vertical anisotropy for pilot points for model layer 4	--	1.15	0.26	3.85
kvani_pp_layer4	Vertical anisotropy for pilot points for model layer 5	--	1.15	0.23	4.85
kvani_pp_layer5	Vertical anisotropy for pilot points for model layer 6	--	1.14	0.21	6.32
kvani_pp_layer6	Vertical anisotropy for pilot points for model layer 7	--	2.26	0.40	8.60
kvani_pp_layer7	Vertical anisotropy for pilot points for model layer 8	--	2.26	0.52	9.27
kvani_pp_layer8	Vertical anisotropy for pilot points for model layer 9	--	2.34	0.41	7.94
kvani_pp_layer9	Vertical anisotropy for pilot points for model layer 10	--	2.34	0.47	9.37
kvani_pp_layer10	Vertical anisotropy for pilot points for model layer 11	--	1.18	0.14	4.21
kvani_pp_layer11	Vertical anisotropy for pilot points for model layer 12	--	1.17	0.28	5.66
kvani_pp_layer12	Vertical anisotropy for pilot points for model layer 13	--	1.12	0.16	3.85
kvani_pp_layer13	Vertical anisotropy for pilot points for model layer 14	--	1.18	0.21	5.79
kvani_pp_layer14	Vertical anisotropy for pilot points for model layer 15	--	1.19	0.20	3.43
kvani_pp_layer15	Vertical anisotropy for pilot points for model layer 16	--	1.07	0.20	2.89
kvani_pp_layer16	Vertical anisotropy for pilot points for model layer 17	--	1.16	0.24	5.13
kvani_pp_layer17	Vertical anisotropy for pilot points for model layer 18	--	1.19	0.23	3.81
kvani_pp_layer18	Vertical anisotropy for pilot points for model layer 19	--	1.16	0.21	2.64
ss_zone	Specific storage by resistivity zone	1 per meter	0.00	0.00	0.00
sy_zone	Specific yield by resistivity zone	--	0.26	0.09	0.35
sy_pp_layer0	Specific yield for pilot points for model layer 1	--	1.13	0.27	3.35
sy_pp_layer1	Specific yield for pilot points for model layer 2	--	1.16	0.18	3.45
sy_pp_layer2	Specific yield for pilot points for model layer 3	--	1.14	0.16	3.08
sy_pp_layer3	Specific yield for pilot points for model layer 4	--	1.15	0.20	3.06
sy_pp_layer4	Specific yield for pilot points for model layer 5	--	1.16	0.24	3.59
sy_pp_layer5	Specific yield for pilot points for model layer 6	--	1.12	0.22	3.00
sy_pp_layer6	Specific yield for pilot points for model layer 7	--	1.11	0.20	3.10
sy_pp_layer7	Specific yield for pilot points for model layer 8	--	1.12	0.16	3.02
sy_pp_layer8	Specific yield for pilot points for model layer 9	--	1.15	0.28	3.35
sy_pp_layer9	Specific yield for pilot points for model layer 10	--	1.15	0.21	3.05
sy_pp_layer10	Specific yield for pilot points for model layer 11	--	1.13	0.20	3.12
riv_cond_mult	River package conductance multiplier	--	0.30	0.30	0.30
ghb_cond_mult	General head boundary multiplier for the North Ozark and South Ozark interface	--	1.66	1.66	1.66
wel_croptype_mult	WEL package multiplier by crop type	--	1.03	0.73	1.42
wel_layer_aemframework_mult	WEL package multiplier by model layer and resistivity zone	--	1.11	0.24	5.12
wel_swirr_mult	WEL package multiplier by surface-water irrigated area	--	0.89	0.00	1.00
wel_fraction_moved	WEL package multiplier defining fraction of pumping moved from the MRVA to model layer 15 of the Middle Claiborne aquifer	--	0.49	0.49	0.49
rzonemult	Recharge multipliers by shallow resistivity zone	--	1.02	0.78	1.17
rch_pp_mult	Recharge pilot point multipliers	--	0.94	0.22	1.20
rchspmult	Stress period recharge multipliers	--	0.39	0.10	1.19
wel_datasource_mult	WEL package multiplier by data source	--	1.03	0.49	1.79
sfr_inflow_mult	Streamflow Routing package inflow multiplier	--	1.00	0.86	1.18
sfr_runoff_mult	Streamflow Routing package runoff multiplier	--	1.03	0.39	2.00
sfr_kv_mult	Streamflow Routing package vertical streambed hydraulic conductivity multiplier	--	0.02	0.00	0.32

The parameter fraction_move_aiwum_mrva2sparta increased about 29 percent from 0.378 to 0.486, which indicated that PESTPP–IES preferred to increase the amount of pumping in the middle Claiborne aquifer to improve the calibration. As discussed in the “Calibration Parameters” section of this report, this parameter moved pumping from MRVA model layers 5 through 8 to layer 15. This parameter likely increased to compensate for the large oversimulation bias in groundwater-level residuals in the middle Claiborne aquifer throughout the simulation. Given the uncertainty in agricultural pumping from the middle Claiborne aquifer in the Grand Prairie region, this was an acceptable parameter adjustment that illuminates the need for better water use information in the middle Claiborne aquifer within the Grand Prairie region.

In the Grand Prairie model, well pumping multipliers by the coupled model layer and AEM framework zone (wel_layer_aemframework_mult) exhibited moderate deviations from their prior values, which illuminated PESTPP–IES’ preference for adjustments to pumping by layer in the Grand Prairie model (Traylor and Weisser, 2024). The mean multipliers ranged from about 0.455 to 5.125 and exhibited a mean value across all layers of 1.11 (table 15). PESTPP–IES increased or decreased multipliers for layers 1 through 15 from 9 to 63.8 percent. Pumping multipliers for layers 16 and 19 were increased by 346.2 and 412.5 percent, respectively, which indicated that pumping in the lower Claiborne aquifer and middle and lower Wilcox aquifers may have been underrepresented in the model prior to calibration. The calibrated well pumping multipliers by crop type for each stress period in the Cache model are shown in figure 8.2.

Posterior Parameter and Simulated Observation Distributions

The distribution of parameter values from the prior (or initial) PESTPP–IES ensemble compared with the posterior (or final) chosen iteration (iteration two for Cache and iteration one for Grand Prairie) ensemble provides an indication of parameter uncertainty and how it changed through calibration. Smaller posterior parameter distributions compared to their initial distributions indicated a reduction in the uncertainty. Conversely, a wider or similar distribution for posterior parameters indicates that the calibration did not reduce the uncertainty of that parameter or group of parameters; perhaps because its parameter identifiability was too low, it was not sensitive to the observation dataset or it already had an optimal distribution (Doherty and Hunt, 2009). The ensemble parameter distribution produces a range of simulated equivalent observations that can also be used to understand the uncertainty of outputs such as simulated groundwater levels. A summary of several salient parameter groups and simulated groundwater levels at priority wells for the calibration process for each model follows.

The Cache model’s posterior ensemble distribution was generally smaller and had reduced parameter uncertainty compared to the prior ensemble. Groups of parameters that exhibited the most reduction in their distribution included Kvani (Kvani_zone parameter group) for zone 9, several of the recharge zone multipliers (rzonemult parameter group), SFR package Kv multipliers by zone (sfr_mult_zone parameter group), and pumping multipliers by data source (wel_datasource parameter group) for all three pumping data sources, which are IWUM (analogous to AIWUM here), MERAS2, and SWUDS (fig. 24A–G). The agricultural pumping data source (IWUM) posterior parameter ensemble exhibited the largest reduction in its distribution, which indicated that the calibration gained the most information from the observations of that parameter and that the prior ensemble distribution was very wide and a reflection of how calibration can reduce uncertainty in important parameters. The calibration did not substantially change the distribution of most Kh and Kvani by zone and pumping by crop type parameters, which indicated that those parameters were likely less identifiable (figs. 9.1, 9.2; Doherty and Hunt, 2009).

The Grand Prairie model exhibited a marginal reduction in posterior ensemble parameter distributions for most parameters because the prior ensemble for many groups of parameters was based on a preliminary PESTPP–IES calibration ensemble from the “MRVA-focused” model instead of previously published values, as in the Cache model. The calibration did not substantially reduce the distribution of most Kh and Kvani by zone parameters, which indicated that those parameters were likely at reasonable values prior to the automated calibration (figs. 9.3, 9.4). Parameters that exhibited the largest reduction from prior to posterior parameter distributions were the well pumping multiplier parameters for the “other” crop type (parameter wel_croptype_mult) and for the fraction of pumping moved from the MRVA to model layer 15 of the middle Claiborne aquifer (parameter fraction_move_aiwum_mrva2sparta; fig. 25A, B). The uncertainty in pumping from the middle Claiborne aquifer within the Grand Prairie region was reflected in the model’s preference to increase pumping in the deeper aquifer whereby the posterior distribution of the fraction_move_aiwum_mrva2sparta parameter became more normally distributed rather than bimodal (fig. 25B).

The reduction in parameter uncertainty from the prior to posterior parameter ensembles for both models resulted in a reduction in the posterior distribution of simulated groundwater levels. In the Cache model, the calibration process reduced the 2007 through 2018 distribution of groundwater levels in observation well 352505090565301, located in the southern section of the Cache CGWA, from about 40 m to about 10 m (a 75 percent reduction) (figs. 26B, 7A). Cache model observation well 352726090523101, also located in the southern section of the Cache CGWA, exhibited a reduction from about 35 m to about 10 m (a 71 percent reduction) (figs. 26A, 7A). A similar reduction was present in the observation ensemble for groundwater levels in the Grand Prairie model within the northern section of the Grand Prairie CGWA where observation wells 344543091510601 and 345129091455801 exhibited a reduction from about 60 m to about 22 m (a 63 percent reduction) (figs. 26C, D, 7B).

Cache Model Groundwater Budgets

The groundwater budgets for each area of the Cache model were characterized by irrigation and nonirrigation season trends in water use whereby the largest outflows were to pumping wells and occurred during the peak growing season (May through September). Across the entire model domain, the primary source of water to the wells was releases from unconfined groundwater storage (fig. 27A; tables 16 and 17). Mean 2007 through 2018 releases from unconfined aquifer storage were 50,886,613 m³, which was about 48 percent of all inflows (table 16). Stream leakage (from streams simulated with the SFR package) and areal recharge were other substantial inflows, constituting about 33 and 9 percent of total inflows, respectively (table 16). The large outflows to wells each irrigation season are due to pumping from agricultural wells for irrigation and aquaculture use; municipal water use did not contribute a substantial amount to total annual outflows. Mean 2007 through 2018 aquaculture, irrigation, and municipal pumping were 54,657,872 m³, constituting about 52 percent of all outflows (table 17).

The groundwater budget for the Cache CGWA Tier 3 was characterized by large outflows to irrigation wells during the late spring and summer months and supported primarily by releases from groundwater storage, which was similar to the entire Cache model domain (fig. 27B; tables 18 and 9). From 2007 through 2018, mean annual outflows to aquaculture, irrigation, and municipal wells were about 64 percent of total outflows (table 19). Releases from unconfined storage were about 49 percent of total inflows. Other substantial inflows included leakage from streams simulated with the SFR package (about 15 percent), areal recharge (about 8 percent), and inflows from the surrounding aquifer (about 28 percent) (table 18). Other substantial outflows included replenishment to unconfined storage (about 34 percent). Outflows to the surrounding aquifer were minimal at about 1 percent (table 19).

The groundwater budgets for the entire Cache model domain and the CGWA Tier 3 indicated that there was more surface water and groundwater connection outside of the Tier 3 region than within (tables 16 and 18). This was evident visually whereby simulated water table altitudes for dry (September 1, 2014) and wet (May 1, 2009) conditions exhibit higher water table altitudes along the Cache River and Village Creek, which indicated losing stream conditions. Water table altitudes along the White River were lower than the surrounding aquifer during dry conditions, which indicated a gaining stream (fig. 28A,B). During wet conditions, the White River exhibited a losing stream owing to the higher water table altitudes than surrounding aquifer (fig. 28C,D). This behavior along the White River indicated that river stage was important in determining the level of connection with the underlying aquifer. Simulated stream leakage during dry conditions exhibited a number of stream reaches drying out primarily along the central and southern reaches of Bayou De View and tributaries of the L’Anguille River, whereas dry stream reaches during wet conditions were simulated mainly for small tributaries of Bayou De View and the L’Anguille River (fig. 28B,D).

Grand Prairie Model Groundwater Budgets

The groundwater budgets for each area of the Grand Prairie model, including the MRVA and the Tertiary-age aquifers, are characterized by irrigation and nonirrigation seasonal trends in water use whereby the largest outflows were to pumping wells that operated during the peak growing season (May through September) to irrigate crops. Across the entire model domain, the primary source of water to the wells was the release of groundwater from unconfined storage (fig. 29A; tables 20 and 21). The total mean 2007 through 2018 release from unconfined aquifer storage was 46,704,367 m³, which was about 48 percent of all inflows (table 20). Connection with surface water was minimal in the Grand Prairie model throughout the entire domain. The mean annual stream leakage (from the SFR package) was about 218 m³, which is less than 1 percent of total inflows to the groundwater. Areal recharge was the primary source of water and accounted for about 35 percent of total inflows (table 20). The large outflows to wells during each irrigation season were due to pumping from agricultural wells for irrigation and aquaculture use, whereas municipal and thermoelectric water use did not contribute a substantial amount to total annual outflows. Mean outflow from 2007 through 2018 for aquaculture, irrigation, and municipal pumping was 52,053,330 m³, about 54 percent of all outflows (table 21).

The groundwater budget for the Grand Prairie CGWA Tier 3 was similar to the entire Grand Prairie model domain and was characterized by large outflows to irrigation wells during the late spring and summer months and supported primarily by releases from groundwater storage and limited connection to surface water (fig. 29B; tables 22 and 23). From 2007 through 2018, the mean annual outflow to aquaculture, irrigation, and municipal wells was about 56 percent of total outflows (table 23). Releases from unconfined storage were about 48 percent of total inflows (table 22). Other substantial inflows included areal recharge (about 28 percent) and inflows from the surrounding aquifer (about 24 percent) (table 22). Other substantial outflows included replenishment to unconfined storage (about 35 percent) and outflows to the surrounding aquifer (about 8 percent) (table 23).

Simulated groundwater budgets for the Grand Prairie model indicated minimal connection between surface water and groundwater for reaches simulated with the SFR package, particularly in the Grand Prairie region (between the Arkansas River and White River) (figs. 2B, 30A,B). There was more substantial connection simulated between the Arkansas River, which used the RIV package, and the groundwater system. This surface water and groundwater connection during dry (September 1, 2014) and wet (May 1, 2009) periods was evident in the higher simulated water table along the Arkansas River (fig. 30A,C). The primary cone of depression under the CGWA Tier 3 and between the Arkansas River and White River region was also similar for dry and wet periods. There were more dry stream reaches simulated during dry conditions than in wet conditions (fig. 30B,D); however, the dry conditions were in upper reaches of smaller streams with flows that generally originated within the model domain. During dry conditions, the simulated water table exhibited more drawdown across the model domain from pumping that was visible in the simulated water table altitudes (fig. 30A,C). Surface-water irrigation effects on the water table were also more prominent during wet conditions where water stored in irrigation reservoirs supplemented pumping in the region with the most drawdown (fig. 30A,B).

Pumping in the Tertiary-Age Aquifers

Both models simulated pumping in the Tertiary-age aquifers below the MRVA. For the Grand Prairie model, the magnitude of pumping was similar in the MRVA (model layers 1–8) and Tertiary-age aquifers (model layers 9–19), including the middle Claiborne aquifer (fig. 31A–C). This was a deviation from previous conceptual and numerical models of the Grand Prairie region in Clark and others (2011a,b), which stated that pumping in the middle Claiborne aquifer of northeastern Arkansas was one to two orders of magnitude less than the pumping from the MRVA. The total amount of pumping in the MRVA was about 5 percent less than the pumping from the Tertiary-age aquifers from April 2007 through December 2018. The mean monthly total pumping was about 2,162,345 m³ for the MRVA and about 2,267,726 m³ for the Tertiary-age aquifers. The monthly trend in pumping indicated that peak pumping from both aquifers generally occurred in July, which is in accord with the typical irrigation trends for crops. Mean monthly pumping volumes were largest for wells in the MRVA from June through September, whereas municipal pumping in the Tertiary-age aquifers accounted for larger mean monthly pumping volumes from October through May (fig. 31D). The conceptual change in pumping for the Tertiary-age aquifers, including the middle Claiborne aquifer was due to the movement of pumping from the MRVA to the middle Claiborne aquifer to improve the calibrated model’s fit to observed groundwater levels in the Tertiary-age layers of the model; this approach was deemed reasonable based on the uncertainty in Tertiary-age and middle Claiborne aquifer pumping in the Grand Prairie region, as described in the “Calibration approach” and “Groundwater-Flow Model Calibration Results” sections of this report.

Spatial distributions of calibrated pumping for the Cache model and Grand Prairie model MRVA and Tertiary-age aquifers are shown in figure 32A, B, respectively. In the Cache model, the MRVA (model layers 1–8) pumping primarily occurred in the central and western side of the model domain, and the deeper Tertiary-age aquifer (model layers 9–18) pumping occurred in the central and eastern side of the model domain (fig. 32A). In the Grand Prairie model, pumping from the MRVA and Tertiary-age aquifers was widespread across the model domain (fig. 32B).