Physiography

The entire study area lies within the Great Plains physiographic province (Fenneman, 1931). Several subdivisions of the Great Plains province are present within the study area and can be differentiated by topographic regions that include the Bluffs and Escarpments, Dissected Plains, Plains, Sandhills, and Valleys (fig. 1B; Conservation and Survey Division, 2019). The Dissected Plains (also known as the Dissected Loess Plains) is present in portions of Buffalo, Custer, Dawson, Frontier, Gosper, Lincoln, and Logan Counties (fig. 1B) and are characterized by level to rolling hills and loess bluffs that grade into the Platte River valley to the south and east (Weaver and Bruner, 1948). Areas of considerable relief with hills of 100 ft or more were created by wind and stream erosion and produced the dissected nature of the loess deposits (Weaver and Bruner, 1948). The Dissected Plains in Frontier and Gosper counties feature deeply incised canyons and flat uplands known as “tablelands.” The dominant soil type is silt that includes the Holdrege Silt Loam and Colby Silt Loam (Weaver and Bruner, 1948; University of Nebraska-Lincoln, 2018). The Plains region located in portions of Adams, Clay, Gosper, Hall, Hamilton, Howard, Kearney, Phelps, Polk, and York counties is characterized by gently rolling hills with less relief than the Dissected Plains. The Sandhills region is present primarily in outlying portions of the main Sandhills region in portions of Buffalo, Hall, Howard, Kearney, Lincoln, Merrick, and Phelps counties (fig. 1B). The Valleys region is primarily located along the Platte River valley lowlands in portions of Buffalo, Hall, Howard, Kearney, Lincoln, Merrick, and Phelps counties (fig. 1B). The Valley lowlands are typically flood-plain areas and have minimal topographic relief with an eastward slope of 6 to 8 ft per mile (Peterson, 2009; fig. 1A–B).

Climate

The climate in the study area transitions from subhumid continental classification in the east, specified in Köppen (1936), to dry subhumid continental classification in the western part of the study area (Conservation and Survey Division, 1998); that is, annual precipitation decreases westward throughout the study area. The climate is characterized by warm and humid summers and cold and windy winters. The 30-year (1981 to 2010) average annual climate normal for precipitation is 29.12 inches at Columbus; 25.23 inches at Kearney; and 23.71 inches at Gothenburg (fig. 1A) (National Climatic Data Center, 2019). The 30-year (1981 to 2010) average climate normal for summer minimum and maximum temperatures is 63.5 and 74.3 degrees Fahrenheit (°F) at Columbus; 60.3 and 72.4 °F at Kearney; and 60.4 and 72.9 °F at Gothenburg (National Climatic Data Center, 2019). The 30-year (1981 to 2010) average climate normal for winter minimum and maximum temperature is 15.7 and 25.1 °F at Columbus; 14.9 and 26.2 °F at Kearney; and 16.0 and 27.9 °F at Gothenburg (National Climatic Data Center, 2019).

Reference ET (ET_ref) was measured at University of Nebraska-Lincoln Extension field sites at Clay Center, Nebr. (not shown) located in the center of Clay County and in North Platte, Nebr. (not shown), about 25 miles west of Brady from 1983 through 2003 (fig. 1A; Irmak and Skaggs, 2011). Average annual ET_ref at Clay Center was 61.7 inches and 62.4 inches for North Platte, which is substantially higher in the study area than average annual precipitation (24.3 inches) and is typical for the subhumid climate classifications. Average monthly ET_ref rates exceed precipitation rates throughout the year measured at locations near Clay Center and North Platte (Irmak and Skaggs, 2011; National Climatic Data Center, 2019). Average annual AET rates for a period between 2000 and 2009 are about 24.8 inches in the study area (Szilágyi and Kovacs, 2010). However, local AET values generally exceed precipitation by as much as 130 percent for areas with widespread irrigation of crops, whereas areas with natural vegetation generally exhibit less AET compared to precipitation (Szilágyi and Kovacs, 2010). Further, AET is generally lower in the eastern part of the study area than the western area for natural vegetation while irrigated crops exhibit similar AET values across the study area (Szilágyi and Kovacs, 2010).

Land Use, Crop Coefficients, and Water Use

Land use and water use in the study area are linked and are important characteristics of the supply and demand driven hydrologic system. Prior to the Civil War (1861 to 1865), the primary land uses in Nebraska were range, pasture, and grass (Hiller and others, 2009). The adoption of the Homestead Act in 1862, the end of the Civil War in 1865, and Nebraska’s statehood in 1867 encouraged settlers to move westward into the study area. In 1895, total cropland area was about 90 percent of 2005 total cropland area (Hiller and others, 2009). Since the mid-1960s, crop diversity in Nebraska decreased, and cropland is now dominated by corn and soybeans (Hiller and others, 2009).

The primary land use since the 1890s has been cultivated crops such as corn, soybeans, and winter wheat. By 1895, canals were developed to deliver surface water for crop irrigation to areas in Buffalo, Dawson, Phelps, and Kearney counties (Hiller and others, 2009; Peterson, 2009). The first irrigation wells were drilled and began pumping groundwater from the alluvial aquifer along the Platte River around 1900 (Nebraska Department of Natural Resources, 2017). Most irrigated land prior to 1940 was irrigated with surface water from canals diverted onto fields using flood irrigation techniques, but after 1940, with improvements in well-drilling technology and later the invention of the more efficient center-pivot irrigation systems, groundwater-irrigated land increased as dryland crops were converted to irrigated cropland. Since 1940, about 29,000 irrigation wells have been drilled in the study area, with most drilled between 1955 and 1990 (fig. 3; Nebraska Department of Natural Resources, 2017). In 2016, the density of irrigation wells in the study area was 3.8 wells per square mile; because of this high density, the CPNRD does not allow the drilling of new irrigation wells or development of new irrigated acres unless other irrigation wells or acres are retired (Central Platte Natural Resources District, 2019). Some irrigators have surface-water rights and groundwater wells for irrigation; the irrigated land that receives water from surface water and groundwater is described as “commingled.” Irrigators with commingled land primarily use their surface-water right and may supplement with groundwater from wells when necessary.

Within the study area, long-term land use data were available from the Cooperative Hydrology Study (2017). These data indicate that by 1950, land use in the study area was about 53 percent pasture (2,632,858 acres), 29 percent dry cropland (1,658,659 acres), 13 percent other uses (406,679 acres), and 5 percent irrigated cropland (227,884 acres; Cooperative Hydrology Study, 2017). The other land uses include open water, urban, riparian forest and wetlands, and roads. Between 1950 and 2016, dry cropland and pasture were converted to irrigated land supplied by groundwater wells. In 2016, the distribution was 47 percent irrigated land with 43 percent irrigated by groundwater (2,289,189 acres), 39 percent pasture (1,926,600 acres), 10 percent dry cropland (495,197 acres), and 4 percent other uses (215,094 acres; Cooperative Hydrology Study, 2017; fig. 4A, B). The development of irrigated cropland from 1950 to 2005 has taken place primarily from Buffalo County in the central part of the study area to Polk and York Counties in the eastern part of the study area (fig. 4B).

Crop coefficients (Kc), which are properties ascribed to plants and used to estimate AET, vary depending on land use or crop type. Allen and others (1998) reported Kc values for the two most common land uses, rangeland and corn, that vary from 0.15 to 1.05 and 0.15 to 1.20 for different parts of the growing season cycle. Crops typically are irrigated between May and September each year based on planting dates and harvesting dates from the U.S. Department of Agriculture (1997 and 2016). Further, the largest amounts of irrigation typically occur in July and August when there is the most difference between ET_ref and precipitation. The fraction of AET that comes from plant transpiration also varies across a single year, much like Kc values. For example, AET during the nongrowing season is predominantly from the evaporation because plants are absent or dormant. Alternately, AET is predominantly transpiration during the middle of the growing season when leaf area index is largest, and more plant surface area is available to transpire.

The primary water use has evolved since the late 1800s from surface-water irrigation to groundwater irrigation. In the late 1890s, surface water for irrigation was the primary water use by way of several canals (Peterson, 2009; fig. 1A). Groundwater irrigation development prior to 1940 was limited to about 760 wells (fig. 3) that irrigated about 158,000 acres. Based on 5-year water use survey data, groundwater irrigation became the primary water use after 1950, and from 1985 to 2015 the total groundwater withdrawals for irrigation were 480 to 840 million gallons per day (Mgal/d) or about 75 to 93 percent of the total water use in the Buffalo, Dawson, Hall, and Merrick counties (U.S. Geological Survey, 2017; fig. 5). Prior to the development of center-pivot irrigation systems in the early 1950s, groundwater pumped for irrigation was applied using less efficient flood irrigation methods. Flood irrigation efficiency prior to 1940 was assumed to be 50 percent (Irmak and others, 2011).

Groundwater use for irrigation within the study area is affected by precipitation; irrigation withdrawals are generally higher during years of lower-than-average precipitation and less in years of higher-than-average precipitation. For example, the reduction in groundwater use for irrigation in 2010 can be attributed to an increase in precipitation compared to other years. Surface water used for irrigation accounted for about 32 to 132 Mgal/d or about 2 to 20 percent of total use (fig. 5). All public supply is from groundwater, and those withdrawals accounted for about 19 to 31 Mgal/d or about 3 to 5 percent of total withdrawals (fig. 5). In 2016, the primary water use in the study area was groundwater for irrigation and the secondary water use was groundwater used for public supply (Dieter and others, 2018).

Surface Water

The surface-water network consists of man-made irrigation canals and natural streams generally flowing west to east for five subregional watersheds: the Big Blue, Little Blue, Loup, Platte, and Republican Rivers (fig. 1C). The Platte River watershed is the primary watershed that constitutes the central region and includes the major streams such as the Platte River, Wood River, and Prairie Creek (fig. 1C). The Platte River flows eastward through the study area from Brady, Nebr., in the west through Kearney, Nebr., and Grand Island, Nebr., in the central region, before flowing out of the study area a few miles east of Duncan, Nebr. (not shown) (fig. 1A). The Platte River has the largest average annual flows in the study area. The mean annual streamflow at the Platte River near Duncan, Nebr. (USGS streamgage 06774000) is 1,420 cubic feet per second (ft³/s) for the period of record from 1895 to 2016 (U.S. Geological Survey, 2017). The Platte River is a braided stream often with two or three main channels for much if its path through the study area (Alexander and others, 2013). The South Loup and Loup Rivers flow along the northern boundary of the study area and have more tributaries draining from the north than in the study area (fig. 1C). The CPNRD boundary is approximately coincident with the Platte River watershed between Gothenburg, Nebr., and Columbus, Nebr. (fig. 1A). The Big Blue watershed is in the eastern part of the study area and includes Big Blue River and minor streams such as the Lincoln Creek and the West Fork of the Big Blue River (fig. 1C). The southern part of the study area includes the Little Blue River watershed in the east, which is drained by minor streams such as Cottonwood Creek (not shown) and Big Sandy Creek, and the Republican River watershed in the southwest, which includes minor streams such as Deer Creek and Muddy Creek (fig. 1C). Streams that flow out of the study area include the Platte River near Columbus, Nebr.; the Big Blue River in eastern Polk County; and Lincoln Creek and Beaver Creek in eastern York County (fig. 1C). Muddy Creek and Deer Creek flow out of the study area in Frontier County (fig. 1C). The spatial location of each stream was derived from the National Hydrography Dataset (McKay and others, 2012).

Some reaches of streams leak stream water into the groundwater system, primarily in the central and eastern region (Peterson and Carney, 2002). Stream leakage depends on stream physical properties such as vertical hydraulic conductivity of the streambed, streambed thickness, channel width, and the hydraulic gradient between the stream stage and the groundwater. Calibrated vertical streambed hydraulic conductivity from groundwater-flow models developed in Peterson (2009) and Peterson and others (2016) ranged from 0.1 to 10 feet per day (ft/d). Data were unavailable to define streambed thickness; therefore, streambed thickness was assumed to be a uniform constant value of 3 ft for all streams. Stream channel widths were defined using recent areal satellite imagery from Google Earth (Google, 2018). Streambed hydraulic conductivity was also assumed to be related to the predominant soil type in that location because streams in the study area are shallow and generally do not cut into bedrock, except in the southwestern part of the study area (fig. 1B).

A network of eight canals divert surface water from the Platte River for irrigation. The earliest canals began diverting water in 1895 and include Cozad Canal, Dawson County Canal, Gothenburg Canal, Kearney Canal, Orchard-Alfalfa Canal, and Sixmile Canal (fig. 1C). Elm Creek Canal (not shown) and Thirtymile Canal began diverting water from the Platte River in 1932 and 1928, respectively (fig. 1C). By 1940, the eight canals (Cozad Canal, Dawson County Canal, Orchard-Alfalfa Canal, Gothenburg Canal, Sixmile Canal, Kearney Canal, Thirtymile Canal, and Elm Creek Canal) were operating in the CPNRD with the water rights to divert an estimated 200,000 acre-feet per year (acre-ft/yr) (Peckenpaugh and others, 1987; Nebraska Department of Natural Resources, 2019), which is similar to the estimated annual diversion of 193,000 acre-feet in recent years; diversion amounts can vary slightly from year to year (Peterson, 2009). After 1940, the Central Nebraska Public Power and Irrigation District (CNPPID) constructed the Tri-County and Phelps Canals, which diverted surface water from the Platte River west of the study area but had about 130,000 acre-feet leakage through the canal and lateral beds within the study area each year (Peterson, 2009). An estimated 40 percent of the diverted water to the CPNRD canals, about 80,000 acre-ft/yr, leaks through the canal and lateral beds and recharges the aquifer (Peterson, 2009; Peterson and others, 2016). The annual leakage has contributed to increases in base flow of the Platte River and some tributaries in the area (Peckenpaugh and others, 1987). Two reservoirs, Johnson Lake and Elwood Reservoir, built in 1941 and 1974, respectively, exhibit stages above the water table in the area and leak water to the underlying aquifer (Central Nebraska Public Power and Irrigation District, 2019). Hydraulic head values for Johnson Lake and Elwood Reservoir were determined using the mean lake/reservoir stage for each stress period after construction (Central Nebraska Public Power and Irrigation District, 2019).

Hydrogeology and Groundwater

The Northern High Plains aquifer is a part of the High Plains aquifer (Peterson and others, 2016) and constitutes the primary groundwater aquifer in the study area. The geologic units in the study area consist of Quaternary-age valley-fill deposits, dune sand, loess, and alluvium, and Tertiary-age Ogallala Formation silt and sandstone (Gutentag and others, 1984). The most recent units are the Holocene-age valley-fill deposits located in stream valleys that consist of coarser materials such as sands and gravels, and the wind-blown dune sand deposits present in isolated parts of the study area (fig. 6). The wind-blown Pleistocene loess deposits are present throughout the study area and contain silt and fine-grained sands and clays. Pleistocene alluvial deposits, commonly referred to as “Paleo-channels,” also are present throughout the study area as stream deposits that consist of sands and gravels. The geologic units with similar hydraulic properties such as water storage and permeability were grouped into hydrostratigraphic units (HUs) and described in Cannia and others (2006) for incorporation into groundwater-flow models developed in Luckey and Cannia (2006), Carney (2008), and Peterson (2009). The spatial distribution of HUs in the study area is presented as a fence diagram in figure 36 of Cannia and others (2006).

The Ogallala Formation is the principal geologic unit that forms the Northern High Plains aquifer (Gutentag and others, 1984). The Northern High Plains aquifer is the primary aquifer in the study area, is hydrologically connected to the Quaternary-age alluvial aquifers, and includes the Quaternary-age deposits and Tertiary-age Ogallala Formation (Cannia and others, 2006). The study area contains HUs 1–6 from Cannia and others (2006), where HUs 1 and 2 consist of Quaternary-age valley-fill, loess deposits, and alluvial aquifers; HUs 3 and 4 consist of the finer grained and less permeable Quaternary-age loess deposits and Upper-Tertiary-age portions of the Ogallala Formation; and HUs 5 and 6 consist of the sands, sandstones, silts, and gravels of the Ogallala Formation. Although HUs 3 and 4 are less permeable deposits, they are not confining units between coarser HUs 1 and 2 and HUs 5 and 6; therefore, each of the six HUs are hydrologically connected within the study area. The base of the Northern High Plains aquifer in the study area is the Cretaceous-age Pierre Shale, which is HU 10 (fig. 6). The Nebraska portion of the Northern High Plains aquifer has exhibited a loss of 6 million acre-feet of recoverable storage from predevelopment to 2015 (McGuire, 2017). The Northern High Plains aquifer serves as the source for all groundwater irrigation and public supply wells in the study area.

The groundwater is characterized by west to east regional groundwater flow (Peterson and others, 2016). Within each watershed in the study area, consistent with the regional flow system, groundwater typically flows west to east and either flows out of the study area on the eastern edge or discharges as base flow to streams. Locally, groundwater flows toward streams in the western portion of the study area where steams are gaining flow from groundwater, and groundwater flows away from streams in the central and eastern portions of the study area where streams are losing flow to groundwater. Groundwater discharge to streams, referred to as “base flow,” is a component of flow in most streams but not all reaches. Base flow is affected by pumping of wells near streams (Kollet and Zlotnik, 2003).

Aquifer properties, particularly horizontal hydraulic conductivity (Kh) and specific yield (Sy) have been evaluated in previous studies. Houston and others (2013) estimated Kh and Sy at test holes in the Northern High Plains aquifer, using lithologic logs to interpret Kh and Sy for vertical intervals of the aquifer. Peterson (2009) included the calibrated Kh for HUs 1 through 6; HUs 1 and 2 had Kh values of about 10 and 155 ft/d, respectively (table 1; table 3 in Peterson, 2009). The intervals representing HUs 3 and 4 consisted of predominantly silt. These two HUs had similar Kh of about 8 ft/d and acted as a semiconfining unit for most of the study area (Peterson, 2009). The average Kh for the interval representing HUs 5 and 6 was about 33 and 10 ft/d, respectively (Peterson, 2009).

Mean Kh estimates for test holes were highest for the interval representing HUs 1 and 2 at 66.6 ft/d (table 1; Houston and others, 2013). Mean Kh for the intervals representing HUs 3–6 were similar; however, the range in Kh for the upper of these intervals was from less than 1 to 325 ft/d and for the lower interval from about 5.6 to 67.7 ft/d, indicating that the lower interval sediments are much more homogeneous than the interval representing HUs 3 and 4 (table 1). Sy values were similar to the average Sy for the Northern High Plains aquifer of about 0.15 (table 1; Houston and others, 2013; McGuire, 2017).

Additionally, recharge rates to groundwater through the unsaturated zone ranged from 0.2 to 10 inches per year (in/yr) based on values measured across the CPNRD on irrigated land, dryland, and rangeland (Steele and others, 2014). Calibrated recharge simulated in Peterson (2009) had a range of about 1 to 7 in/yr (see fig. 16 and table 9 from Peterson, 2009). Calibrated recharge simulated in Peterson and others (2016) had a range of about 2 to 5 in/yr (see fig. 20B from Peterson and others, 2016).

Climate and Landscape Components

The landscape in the study area includes land-use characteristics such as crop type, ET characteristics for each crop type, soil types, and surface-water characteristics. The landscape is characterized by an extensive surface-water network that includes major streams as described in the “Surface Water” section of this report. The climate and landscape water subsystems include inflows from precipitation; root uptake; inflows from surface-water deliveries or groundwater used for irrigation, runoff, and deep percolation past the root zone; and outflows of ET (broken down by source of water).

Precipitation is the largest inflow to the landscape, although locally the inflows from irrigation can exceed precipitation, particularly in drought conditions. Most surface-water deliveries occur to meet irrigation demands between May and September (table 1.1). Outflows from the landscape subsystem include ETp, ET of irrigation water, ETg that passes through the landscape, deep percolation past the root zone, and runoff of precipitation to streams. ET of irrigation water is generally the largest outflow from the landscape subsystem. Measured annual AET rates varied based on land cover and precipitation. Rangeland had the highest rates of AET that range from 22 to 36 in/yr depending on precipitation and location (Irmak, 2014). Regional irrigated and dry cropland annual AET rates were about 21.6 and 16.8 in/yr, respectively (Irmak, 2014). The AET of irrigated corn, the predominant crop in the study area, was 25.9 in/yr based on nearby measurements (Kranz and others, 2008). Based on average AET rates and land use reported by the Cooperative Hydrology Study (2017), the estimated annual volume of total AET of precipitation and irrigation was about −8,600,000 acre-feet for pre-1940 and about −9,500,000 acre-feet for the recent development period (2011–16; table 2). The uncertainty associated with these estimates was difficult to quantify owing to the lack of AET rates available for all land uses within the study area and the lack of an uncertainty assessment in the mentioned studies. The range of AET for the recent development period (2011–16), was about −8,800,000 to −11,000,000 acre-ft/yr based on the range of irrigated corn AET and rangeland AET in Kranz and others (2008), Irmak (2014), and Cooperative Hydrology Study (2017). The ET of irrigation water is conceptually difficult to estimate because precipitation and irrigation occur on the same area across similar time periods; therefore, the landscape ET component in table 2 combines ETp and ET of irrigation water sources.

Deep percolation past the root zone (often referred to as recharge) in the study area has been estimated at 0.2 to 2.1 in/yr below rangeland (that is, pasture), 2.7 to 6.3 in/yr below groundwater-irrigated cropland, 10 in/yr below surface-water irrigated cropland, and 0.5 to 2.5 in/yr below dry cropland (Steele and others, 2014). Additionally, recharge in the study area that occurs as leakage from the unlined irrigation canals can exceed 10 in/yr (Peterson and others, 2016). The highest recharge rates in areas without canal influence were about 4 to 6 in/yr in the flat, densely groundwater-irrigated cropland in the Platte River valley (Peterson and others, 2016).

Surface-Water Components

Surface-water inflows include streamflow that enters the study area, groundwater that discharges through the streambed as base flow, and runoff. Outflows include leakage of stream water into the underlying aquifer and outflows of streamflow where streams exit the study area. Runoff from precipitation was estimated for the study area using base flow separation techniques in the USGS Groundwater Toolbox (table 4; Barlow and others, 2014; 2017) from streamgage data for the Platte River, Prairie Creek, Silver Creek, and Wood River (table 4; fig. 1A). Streamgages were present at locations of major stream inflows and outflows in the study area for most or all of the development period, which provided adequate data to estimate the gaged flows in and out of the surface-water system. Minor streams that flowed into the study area that were ungaged along the Loup River boundary streams were estimated with flows of about 5 ft³/s for the entire period of interest based on stream width from aerial imagery and an assumed depth of 1 ft. The largest gaged inflows are at the North Loup River near St. Paul, Nebr. (USGS streamgage 06790500; fig. 1A), which has an average discharge of 853 ft³/s for the period of record 1928 to 2016 (U.S. Geological Survey, 2017). However, the Loup River system is a boundary stream; the largest gaged inflows for a nonboundary stream are at the Platte River at Brady, Nebr. (USGS and Nebraska Department of Natural Resources [NeDNR] streamgage 06766000; fig. 1A) with an average total inflow of about 746 ft³/s for the period from 1939 to 2016 (table 4). The largest outflows are at the Platte River about 5 miles downstream from the Platte River near Duncan, Nebr. (USGS streamgage 06774000: fig. 1A), which has an average discharge of 1,742 ft³/s for the period from 1939 to 2016 (U.S. Geological Survey, 2017). For the pregroundwater irrigation development period (pre-1940), total flow into the study area at the Platte River at Brady, Nebr. (USGS streamgage 06766000) was about 1,541 ft³/s. The outflow at the Platte River near Duncan, Nebr. (USGS streamgage 06774000) for the pre-1940 period was about 1,620 ft³/s; the overland flow portion of total streamflow was estimated to be an outflow of about 497 ft³/s (about 360,000 acre-ft/yr or 1 inch) within the Platte River watershed, based on values obtained from Cooperative Hydrology Study (2017). For the recent period (2011 to 2016), total streamflow and runoff at the Platte River at Brady, Nebr. (USGS streamgage 06766000) were about 1,102 ft³/s and 427 ft³/s, respectively. Recent period (2011–16) total streamflow and runoff at the Platte River near Duncan, Nebr. (USGS streamgage 06774000) were about 2,495 and 884 ft³/s, respectively (table 4). Runoff was estimated to be an outflow of about 460,000 acre-ft/yr (about 635 ft³/s or 1.11 inches) within the Platte River watershed, which contributes to the runoff at the Platte River near Duncan, Nebr. (USGS streamgage 06774000) based on values obtained from Cooperative Hydrology Study (2017; fig. 1A).

Groundwater Components

Inflows to the groundwater-flow system include recharge from the landscape subsystem deep percolation component, recharge as canal leakage from the CPNRD and CNPPID canal systems, stream leakage, and cross-boundary flow into the active model area. Recharge rates are described in the “Climate and Landscape Components” section of this report. Outflows from the groundwater-flow system include cross-boundary flow out of the active model area, discharge to streams as base flow, withdrawals from irrigation wells, withdrawals from municipal wells, and ETg. Transient changes in groundwater inflows and outflows are balanced by increases or decreases of groundwater in storage. When inflows are greater than outflows, aquifer storage increases (rise in groundwater levels), and when inflows are less than outflows, aquifer storage decreases (decline in groundwater levels).

Groundwater-flow subsystem components as annual inflows on the left-hand side and annual outflows on the right-hand side are presented in equation 1:

Groundwater generally flowed into the active model area at the northwestern and southwestern boundary and flowed out along portions of the southwest and eastern boundaries during pre-irrigation development period of 1895 and throughout the irrigation development period after 1895 (fig. 7C). Data were not available at observation wells near the study area boundaries to estimate the amount of groundwater that flowed across the study area boundaries during the pre-irrigation period. Therefore, estimates of cross-boundary flow were derived using Darcy’s Law (Fetter, 2001), where hydraulic gradients were obtained from decadal groundwater levels simulated in Peterson and others (2016) and aquifer transmissivities (Kh and aquifer thickness) from Houston and others (2013) and Peterson and others (2016). Hydraulic conductivities across these boundaries were 16 to 158 ft/d, hydraulic gradients were 0.001 to 0.004, and aquifer thicknesses were about 142 to 588 ft (Houston and others 2013, Peterson and others 2016). These estimates of cross-boundary groundwater flow indicated that groundwater inflow for the pre-irrigation development period (about 1940) was about 50,000 acre-ft/yr and outflow was about 180,000 acre-ft/yr. These estimates of cross-boundary groundwater flow were assumed to be reasonable for the period prior to widespread groundwater irrigation development (about 1940) because withdrawals from the few active irrigation wells prior to 1940 were not enough to cause substantial declines in the water table at the boundaries of the study area and the decadal groundwater-level contours from Peterson and others (2016) showed little change across these boundaries for the April 30, 1940, and 1940s periods (see fig. 15 from Peterson and others, 2016). Estimates of groundwater-flow inflow to the study area for the recent period from 2011 to 2016 were about 50,000 acre-ft/yr and outflows were about 140,000 acre-ft/yr. The decrease in cross-boundary outflows is likely associated with the widespread increase in outflows to irrigation wells after 1940, which is indicated by the migration of groundwater-level contours eastward in some areas shown in the 1940s and 2000s decadal groundwater-level contours in Peterson and others (2016).

The groundwater-flow system is connected to the surface-water system through stream-aquifer interaction. Water-table contour maps indicate groundwater levels were above adjacent stream surfaces in the western portion of the active study area, which results in a net outflow from the groundwater to streams (Summerside and others, 2001; fig. 7C). In the eastern portion of the study area, stream surfaces are generally higher than the adjacent groundwater levels, which results in a net inflow to the groundwater from streams (Summerside and others, 2001); for example, along the Platte River approximately between Cozad and Grand Island, Nebr. (fig. 1A).

The Platte River, the primary stream in the study area, is connected to the groundwater-flow system for the entire length that it flows through the study area (fig. 1A). Along some reaches, the groundwater-flow system discharges base flow to the Platte River, and along other reaches the groundwater-flow system receives inflows from the Platte River as stream leakage (Peterson and Carney, 2002). A groundwater discharge analysis in Peterson and Carney (2002) estimated base flow at several streamgages on the Platte River during the fall low-flow season: the average groundwater discharge as base flow was from −3.0 to 4.0 ft³/s per mile, where negative values indicate a losing reach of the stream. Estimated base flow at Platte River at Brady, Nebr., streamgage (USGS streamgage 06766000) was about 675 ft³/s for the recent development period and 570 ft³/s for the pre-1940 period (table 3). Estimated base flow at Platte River near Duncan, Nebr. (USGS streamgage 06774000) was estimated to be about 1,611 ft³/s for the recent development period and 649 ft³/s for the pre-1940 period (table 3).

Recharge (analogous to deep percolation from the landscape subsystem) is the primary inflow to the groundwater-flow system for the pre-1940 and recent development periods. Recharge rates are highest on irrigated lands and lowest on rangeland (see deep percolation values described in the “Climate and Landscape Components” section of this report and based on measured values from Steele and others, 2014). Stream leakage is also estimated to be a substantial inflow to the groundwater along “losing” stream reaches of the Platte River between Overton (not shown) and Grand Island, Nebr. (fig. 1A). Although groundwater irrigation is the primary groundwater outflow and is typically associated with groundwater-level declines, localized increases in groundwater storage changes for the recent development period from inflows such as canal leakage and decreased outflows to groundwater ET have contributed to a small increase in storage prior to 1940 (McGuire, 2017). The increase in groundwater levels because of canal leakage has been monitored in the study area and has a local influence on the water table that does not reflect groundwater levels in other areas that do not have canal leakage inflows each year (U.S. Geological Survey, 2017).

The increase in irrigated cropland and irrigation well development has resulted in a substantial increase in outflows from groundwater irrigation wells since the 1940s (Peterson, 2009; Peterson and others, 2016). Based on a consumptive use deficit (when available precipitation is less than the amount required for a crop to fully transpire) of about 10 inches per growing season for corn, the dominant crop type, an estimated 20 inches of water was required to be pumped during an average climate year to irrigate about 158,000 acres of cropland (pre-1940 land use) for an estimated annual volume pumped of 263,000 acre-feet (table 3; Irmak and others, 2011). Estimated ranges of groundwater irrigation in table 3 are based on the range of efficiency appropriate for the time period based on values in Irmak and others (2011). Prior to the development of center pivots in the 1950s, flood or furrow irrigation was the common technique to apply groundwater pumped for irrigation, and these efficiencies had a range of 50 to 65 percent. By 2016, the number of irrigation wells was about 32,000 (fig. 3; Nebraska Department of Natural Resources, 2017) and irrigation efficiency improved with development of center-pivot technology; 80 percent of groundwater irrigation is delivered by sprinklers in Nebraska, of which most are more efficient center pivots (Irmak and others, 2011; Johnson and others, 2011). The average rate of groundwater irrigation application decreased to about 10 to 11 in/yr, but with the increase in acreage to about 2,169,000 acres of cropland, the annual estimated volume pumped increased to about 1,990,000 acre-feet (table 3). With the development of irrigation technology, such as the Low Energy Precise Application, the efficiency of groundwater pumping improved to as much as 90 percent (Irmak and others, 2011). Metered data for irrigation pumping were unavailable for more than 99 percent of irrigation wells in the CPNRD, and irrigation technology data such as traditional sprinkler or Low Energy Precise Application were also not available. Lack of these data also contributed to the uncertainty of the conceptual estimates of groundwater irrigation.

Boundary Conditions

A combination of groundwater divides, surface-water features, and available input data extents were used to define the “active” domain and the boundary conditions simulated by the CPIHM around the CPNRD focus area. Simulated groundwater contours from the Northern High Plains groundwater-flow model (Peterson and others, 2016) were the principal source of information on the location of the western and southern hydrologic boundaries surrounding the CPNRD (fig. 7B). The western areal boundary of the study area was selected as the region extending north from central Frontier County through Brady, Nebr. (fig. 1A) to the South Loup River in southeastern Logan County about 10 miles upstream from South Loup River at Arnold, Nebr. (USGS streamgage 06781600, fig. 1A). Based on average decadal groundwater contours from Peterson and others (2016), groundwater flows into the study area across the northern and southern portions of the western boundary for this study, except for a region just north of the Platte River near Brady, Nebr., where the groundwater flows south into the Platte River, and the region in Frontier County where groundwater flows approximately southeast and creates a no-flow groundwater boundary (fig. 15 from Peterson and others, 2016; fig. 7C).

The South Loup River was chosen as the northern boundary of the CPIHM, west of the confluence of the Loup River (fig. 7C). The Loup River forms the northern boundary of the study area to its confluence with the Platte River in Platte County (fig. 7C). Hydrologically, the northern boundary that includes the South Loup and Loup Rivers is a mapped groundwater divide and creates a no-flow groundwater boundary (Peterson, 2009). Groundwater flows across much of the eastern boundary of the study area in Clay, York, and Polk Counties (fig. 7C). The northern and eastern boundaries are coincident with the boundaries from models in Peterson (2009) and Cooperative Hydrology Study (2017).

The southern boundary extended from the groundwater divide west of Deer Creek in Frontier County through central Gosper, Phelps, and Kearney Counties, then along the Big Sandy Creek in Clay County (fig. 7C). The southern boundary of the model approximately follows the groundwater divide created by a groundwater mound under the CNPPID canal system in central Gosper, Phelps, and Kearney Counties that is a result of persistent leakage recharge from the canals (fig. 7C). The southern boundary contains areas where groundwater flows east or southeast out of the study area or flows parallel to the model boundary and creates a no-flow groundwater boundary (fig. 7C). Vertical boundaries of the active domain of the CPIHM were the land surface as the upper boundary and top of the Pierre Shale as a no-flow lower boundary (HU 10 from Cannia and others, 2006; Peterson, 2009; Peterson and others, 2016).

Layering Scheme

The vertical discretization of the CPIHM, between the vertical boundaries (land surface and Pierre Shale), was defined by the HUs described in Cannia and others (2006) where present in the study area. The six HUs (HUs 1–6) present in the study area were combined into three hydrologically similar groups based on permeability, geologic age, and spatial relation to reduce the number of vertical layers that would need to be represented in the CPIHM (fig. 8). HUs 1 and 2 were combined into numerical model layer 1, HUs 3 and 4 were combined into numerical model layer 2, and HUs 5 and 6 were combined into numerical model layer 3 (fig. 8). Hydrologic connection remained between all three layers after grouping (figs. 6 and 8). Model layer 1 is the Quaternary alluvium and loess from HUs 1 and 2. Group 1 is present everywhere in the study area with an average thickness of 163 ft and a range from less than 10 ft in parts of the Platte River valley to greater than 500 ft in parts of Dawson and Custer counties (figs. 7B, 8). The alluvium and loess typically consisted of gravels, sands, and silts. Model layer 2 combined the lower Quaternary and upper Tertiary Silt identified in Cannia and others (2006) as HUs 3 and 4. There were a few areas in the central and eastern portion of the study area where these units were absent (Cannia and others, 2006). Group 2 is thicker toward the east with an average thickness of 57 ft and a range from less than 10 ft in the west and central areas to 300 ft in the eastern portion of the study area.

Model layer 3 is the Ogallala Formation identified in Cannia and others (2006) as HUs 5 and 6. Group 3 also incorporated airborne electromagnetic data from Cannia and others (2017) that were used to refine the aquifer base. The Ogallala Formation consists of sands, silts, and clays at various intervals. The Ogallala Formation is present in the western and central portions of the study area and thins eastward with an average thickness of 242 ft and ranges from less than 10 ft in the eastern portions to 500 ft in the western portions of the study area (figs. 7B, 8).

Spatial and Temporal Discretization

The CPIHM was spatially discretized into a grid of orthogonal blocks, called cells, horizontally and vertically. The orthogonal grid consisted of 163 rows and 327 columns with horizontal cell sides of 2,640 ft by 2,640 ft (0.5 mile by 0.5 mile or 160 acres) (table 5). The CPIHM was vertically discretized into three layers of varying spatial extents and thickness according to the hydrogeology and aquifer characteristics described in the “Study Area Description” and “Layering Scheme” sections of this report. The CPIHM includes 159,903 total cells (fig. 7A) where model layers 1, 2, and 3 include 30,788 active cells; 30,009 active cells; and 17,712 active cells, respectively. Based on the hydrostratigraphic data, layer 1 (Quaternary-age alluvial and loess deposits) was present, and therefore active, throughout the study area (fig. 7A). Layer 2 was active throughout most of the study area but was absent in some parts of the central and eastern model domain. The absence of layer 2 where layers 1 and 3 were present appeared as “holes” in the model grid and will be referred to as such throughout the report. Layer 3 was present in the west and central model domain and thinned out to the east (fig. 7A). MODFLOW–NWT cannot simulate discontinuous layers such as the absence of layer 2 at the “holes”; inactive cells are simulated as no-flow boundaries. Most of the layer 2 holes have layers 1 and 3 present and hydrologically connected. To simulate this hydrologic connection in the CPIHM, cells in the layer 2 holes were activated and given an artificial layer thickness of 5 ft with a hydraulic conductivity equal to the average of layers 1 and 3, and anisotropy was set to a 1:1 ratio to permit vertical flow in the artificial layer 2 (shown as “layer 2 holes” on fig. 7B).

The predevelopment model was run as a 500,000,000-day transient stress period to ensure that an approximate groundwater storage equilibrium was reached (table 6). The groundwater levels generated by the predevelopment model were used as the initial groundwater levels for the development model. The predevelopment model approximated average steady-state groundwater conditions prior to the construction of canals for surface-water irrigation and wells for groundwater irrigation.

The CPIHM development period was temporally discretized into 610 stress periods to simulate transient conditions from May 1, 1895, to December 31, 2016 (table 6). The transient simulation start date was selected as May 1, 1895, because it was coincident with the first documented surface-water diversions for irrigation in the CPIHM active domain. Also, two previous groundwater-flow models that included the CPIHM domain began transient simulations in 1895 to account for the onset of surface-water diversions for irrigation (Peterson, 2009; Peterson and others, 2016). The first 170 stress periods of the development model were discretized into periods of time that reflected the typical irrigation season from May 1 to September 30 and nonirrigation season from October 1 to April 30, with 2-week time steps and were meant to provide a reasonable starting point for the period of interest (April 1, 1982, to December 31, 2016; table 6). Stress periods 171 to 610, from May 1, 1980, to December 31, 2016, were discretized into monthly time periods to improve the temporal resolution of outputs generated for the period of interest, after April 1, 1982, again with 2-week time steps (table 6). Temporal resolution was increased after May 1980 to monthly stress periods to meet the objectives of the study and to facilitate comparison of simulated groundwater levels from different months within an irrigation season with the CPNRD’s baseline groundwater levels from the spring of 1982. The CPIHM scenario period was temporally discretized into 396 stress periods to simulate transient conditions from January 1, 2017, to December 31, 2049 (table 6). Stress periods were discretized into monthly time periods with 2-week time steps to maintain the temporal resolution with the development period model (table 6).

Landscape Inputs and Configuration (Farm Process)

Landscape inputs are necessary for simulation of landscape processes within the MF–OWHM, as outlined in the “MODFLOW–One-Water Hydrologic Model Theory and Approach” section of this report (fig. 2). Landscape properties specified within the FMP for this study included land surface altitude, soil types, ET_ref, precipitation, land use and crop types, Kc, irrigation flags, root depths and pressures, fractions of transpiration (FTRs), fractions of evaporation from irrigation (FEI), fractions of inefficient losses to surface water from precipitation (FIESWP), and fractions of inefficient losses to surface water from irrigation (FIESWI). Surface-water supply sources and routing, the irrigation supply well information, and the water accounting units were specified within the FMP. The water accounting units are referred to hereafter as “water-balance subregions” (WBSs) (Hanson and others, 2014a; Boyce and others, 2020). The WBSs were delineated for the CPIHM based on subregional watersheds, surface-water irrigation zones, groundwater irrigation zones, CPNRD GWMAs, and other natural resource district boundaries. In total, there are 212 WBSs in the CPIHM (Traylor, 2023). The 212 WBSs were combined into 24 “supergroup” zones to facilitate analysis of the water budgets for areas of interest within the CPIHM (fig. 9). The 24 GWMAs within the CPNRD were used to define the supergroups (fig. 9). Hereafter, the term “GWMA” will be used to reference results pertaining to the supergroups.

The land surface altitude dataset was created using a combination of 1- and 2-meter (m) light detection and ranging (lidar) digital elevation model datasets from the Nebraska Department of Natural Resources (2018) and a 10-m digital elevation model dataset from the U.S. Geological Survey (2018) and resampled to CPIHM model cell size. The precipitation (PPT) and ET_ref datasets used in the CPIHM were developed for each cell and stress period, generated from interpolation of climate data at as many as 77 weather stations in the COHYST area (Cooperative Hydrology Study (2017, p. 5–9), then clipped to the study area. The PPT and ET_ref datasets were constructed using the same weather stations and methods described in Cooperative Hydrology Study (2017), which included data from 1950 to 2013. Outside of the 1950 to 2013 period, climate datasets were obtained from the High Plains Regional Climate Center (2018) for each weather station included in the Cooperative Hydrology Study (2017) (fig. 10A, B).

The Cooperative Hydrology Study (2017) land-use dataset specified multiple crop types per model cell in acres, and these values were converted to fractions of model cell area for use in the CPIHM (Traylor, 2023). For example, if crop type “groundwater-irrigated corn” in a specific cell was 64 acres, then the fraction used in the CPIHM cell dataset was 0.4 (for 160-acre cells). The distribution of land uses for 1982 and 2005 is shown in figure 1.1 (Dappen and others, 2007; also see figure 4-D-1 from Cooperative Hydrology Study, 2017). Soil types were grouped into sand, sandy loam, and silty clay based on classifications from a geographic information system soil dataset in University of Nebraska-Lincoln (2018) (fig. 7B). The land-use datasets were derived from those used in Cooperative Hydrology Study (2017) and included Cooperative Hydrology Study (2017) irrigated and dryland crop types for 1950 to 2013 partitioned into commingled irrigated, dryland uses, groundwater irrigated cropland, and surface-water irrigated cropland (table 7). The Cooperative Hydrology Study (2017) crop type “other” was then split into urban, roads, open water, and riparian forest and wetlands using the National Land Cover Database datasets to identify locations of those separate land uses (table 7) (Yang and others, 2018). The 1950 land-use dataset was used for 1895 to 1950, which was an acceptable approximation based on Hiller and others (2009), and the 2013 dataset was repeated from 2014 to 2016.

The Kc values were specified for each crop type and stress period and interpolated and scaled to the various stress period lengths. When possible, daily Kc values were calculated following the methods used in the Cooperative Hydrology Study (2017). Additionally, the out-of-season value for cropland used in Cooperative Hydrology Study (2017) was used to represent fallow land and open water during the winter season, which was specified as November 15 to March 15. Monthly riparian forest Kc values from Hall and Rus (2013) were used for forest values. Published Kc values, planting dates, and crop growth period lengths were taken for all other land-use types from Allen and others (1998). Daily Kc values were calculated for all land uses by interpolating between the planting dates and growth stage dates. Monthly Kc values were then obtained by averaging the daily rates for the given month. Additional sources for Kc values from the Nebraska Agricultural Water Management Network (2018) were compared to the previously mentioned monthly Kc values and, where possible, were adjusted as needed to align with published values. Kc values assigned to each land use and stress period are available in the model archive associated with this report (Traylor, 2023)

Each land use was flagged for irrigation or no irrigation (dryland) for each stress period. Additionally, this “irrigation flag” included the irrigation efficiency type for the land use referenced from the separate “on-farm efficiency” (OFE) input table that specifies the irrigation efficiency as a fraction of the irrigation water pumped that is used by the crops (Boyce and others, 2020). Irrigation flags were designated for each commingled, groundwater, and surface-water crop type for potential irrigation stress periods (that is, May through September for irrigated crop types except winter wheat, which was flagged for potential irrigation from March to October). Irrigation efficiencies were specified in the development period model as 0.7, 0.55, and 0.8 for commingled, surface-water irrigated, and groundwater irrigated crops in the bi-annual irrigation stress periods from 1895 to 1980, respectively; and 0.72, 0.65, and 0.9 for commingled, surface-water irrigated, and groundwater irrigated crops in the monthly irrigation stress periods from May 1980 to December 2016, respectively, based on values from Irmak and others (2011).

Fractions of transpiration (FTRs) and FEI were assigned for each crop type. FTRs and FEI for common crop types like corn were derived from models that predicted corn transpiration and ET (Kimball and others, 2016). FTRs and FEI for winter wheat were derived from a study in China (Kang and others, 2003) and Nebraska (Irmak and others, 2016). Root depths for each land use ranged from 1 to 15 ft and were obtained from Irmak and Rudnick (2014) and U.S. Department of Agriculture (2016). Root pressure values as depths were obtained from table 7 in Hanson and others (2014b). The FIESWP and FIESWI values were chosen based on values of similar parameters within Cooperative Hydrology Study (2017). The runoff calculated by FIESWP or FIESWI was routed to a nearby stream cell specified in the semirouted returns feature of FMP.

The development model simulated the proliferation of groundwater irrigation wells in the study area. The irrigation wells that supply the CIR for each groundwater-irrigated crop were specified in the “farm wells” block within the FMP. Each irrigation well was assigned a pumping location in the model grid (layer, row, and column designation; fig. 11) and a WBS to which the well may supply water in addition to a maximum pumping rate. Farm wells that intersected stream cells were removed to improve model stability. Irrigation well start and end dates were based on construction dates and decommission dates in the well database (Nebraska Department of Natural Resources, 2017). Surface-water deliveries from canal diversions were specified as an available source of water for 36 WBSs (156–189, 211, and 212) because they were coincident with the surface-water irrigated acres. The amount of surface deliveries available each stress period depended on the diversion amount from the streamflow routing network to the canals (Nebraska Department of Natural Resources, 2019).

MODFLOW Inputs and Configuration

MF–OWHM utilization of the Newton solver for the CPIHM required the use of the Upstream Weighting package (Niswonger and others, 2011), including hydraulic properties such as Kh, anisotropy, Sy, and specific storage (Ss) for each model layer. The precalibration initial values for the hydraulic properties are described in the “Hydrogeology and Groundwater” section. In an unconfined system, Sy is a more influential property than Ss because it quantifies 3 or 4 orders of magnitude more water. Further, assignment of a single Ss value across all active cells for each layer was a technique taken from the Northern High Plains aquifer model calibration from Peterson and others (2016). The aquifer storage (Ss and Sy), Kh, and anisotropy values were adjusted during the calibration process.

The General-Head Boundary package (GHB; Harbaugh and others, 2000) simulated inflow and outflow of groundwater across boundaries of the active model (640 cells) and the interaction of Johnson Lake and Elwood Reservoir with the groundwater-flow system (31 cells; fig. 1A; table 5). The GHB is a head-dependent flux boundary that requires a specified hydraulic head and conductance for each GHB designated cell. The flux between a GHB and adjacent cell is the product of the GHB cell conductance and the difference in hydraulic heads between the GHB cell and the adjacent cell. The conceptual model of groundwater flow (in the “Boundary Conditions” section of the report) indicated that there was lateral groundwater flow into the model domain along some portions of the western model boundary and lateral outflow along some southern and eastern portions of the model based on simulated groundwater-level contours from Peterson and others (2016) (fig. 7B). The GHB cells were assigned for these inflow and outflow edges of the model (fig. 12). The hydraulic head values were determined using a combination of rasterized decadal average groundwater levels from Peterson and others (2016) and interpolated heads from measured groundwater-level data at observation wells coincident or near boundary cells. Initial GHB heads were extracted from the rasterized grids of the decadal groundwater-level contours from Peterson and others (2016) and checked against measured groundwater levels when available (U.S. Geological Survey, 2017). For differences greater than 5 ft, GHB head values were specified by extracting values from raster grids interpolated from the measured data near boundary cells. Hydraulic conductance along the boundaries were estimated based on local aquifer properties and layer thicknesses (Peterson, 2009; Houston and others, 2013). Lake and reservoir stage data were used to specify the GHB heads for both reservoirs. Lake GHB conductance was manually adjusted such that the simulated groundwater levels below the lakes were similar to the measured groundwater levels from a nearby observation well (USGS 404046099504501) located on the southern shore of Johnson Lake (U.S. Geological Survey, 2017). Because of the proximity of Elwood Reservoir to Johnson Lake, the GHB conductance used for Elwood Reservoir was the same as that defined for Johnson Lake.

The Streamflow Routing package (SFR; Niswonger and Prudic, 2005) simulated the flow and routing of surface water in streams and canals and their interaction with the landscape and groundwater-flow system within the study area and along the northern boundary. The SFR is a head-dependent flux boundary that interacts with the groundwater-flow system of each designated SFR cell based on the relation among the stage of the stream, streambed conductance, and the hydraulic head of the aquifer. Flow of water between the SFR and the groundwater-flow system is dictated by the conductance, which is the product of the hydraulic conductivity of the streambed and the area of the stream channel, divided by the streambed thickness. The SFR network included 4,318 cells, or reaches, and 219 segments, which are a group of reaches with uniform hydraulic properties, to represent major streams such as the Platte River, South Loup River, Middle Loup River (not shown), Loup River, Wood River, Prairie Creek, Big Blue River, and the perennial reaches of their perennial tributaries (fig. 12). The upper perennial reaches of the tributaries included in the SFR network were determined based on density of streams using geographic information system techniques and assessment of satellite imagery to identify segments of streams that interact with the groundwater based on identification of a wet channel and identification of vegetation along the stream that is assumed to be supported by groundwater. Flows and stage are calculated with the SFR as the difference among inflows from upstream reaches; runoff; tributaries; groundwater discharge to the reach as base flow; and outflows to downstream reaches, canal diversions, and leakage into the aquifer. The SFR simulated stream depth using Manning’s equation and assumed a rectangular channel (Niswonger and Prudic, 2005). Therefore, physical stream properties such as vertical hydraulic conductivity of the streambed, stream width, streambed roughness coefficient, and streambed thickness were specified for each segment. The Gage package (Merritt and Konikow, 2000) was used to extract time series of simulated streamflow for some SFR reaches for calibration to measured flow data.

Canals within the CPNRD focus area, and their associated diversions from the Platte River, were simulated with the SFR (fig. 12). The SFR simulated the routing of diverted streamflow into Cozad, Dawson, Gothenburg, Kearney, Orchard-Alfalfa, Sixmile, and Thirtymile canals (Nebraska Department of Natural Resources, 2019). Stress period averaged canal diversion values were created from daily diversions data available from the Nebraska Department of Natural Resources (2019) and were input to the SFR (table 1.1). Note that the CNPPID canal system located south of the Platte River was not simulated with the SFR. The leakage of diverted water through the canal beds of the CNPPID canals was simulated with the Recharge package (Harbaugh and others, 2000), wherein these recharge rates were estimated using change in groundwater-level values from measured groundwater-level hydrographs near the canals and mean Sy values of the aquifer. Additionally, the Recharge package was used to simulate the leakage from B-1 Reservoir in Dawson County (fig. 12); recharge rates were calculated based on the difference between reservoir inflow and evaporation data to estimate leakage through the lakebed of about 85 percent of inflows (Nebraska Department of Natural Resources, 2019).

Groundwater withdrawals for municipal supply (fig. 11) were simulated using the Well package (Harbaugh and others, 2000). The municipal well spatial and temporal information were obtained from the NeDNR Registered Well Database (Nebraska Department of Natural Resources, 2018). Well depths were used to determine the model layer from which the wells were pumped, and construction dates were used to determine the first active stress period of each well. Withdrawal rates were estimated for wells in communities or time periods without withdrawal data using withdrawal data from municipalities with similar populations. A total of 198 municipal wells were simulated (fig. 11). Municipal wells are hereafter referred to as production wells.

Calibration Parameters

A total of 740 parameters were specified for the predevelopment and development CPIHM within six parameter groups, and 435 of those parameters were adjustable during the automated calibration process (table 8). Parameterization (the specification of model inputs as calibration variables, or parameters that PEST may adjust) was based on prior data availability, local knowledge, and uncertainty in their values. In environmental modeling, there is usually a lack of available data to quantify model inputs for every model cell and stress period. The parameterization scheme (that is, the structure and number of parameters used to represent the unknown natural-system input values of a model) is an important aspect of modeling because it can have a substantial effect on the ability of a model to gain information from the calibration dataset and appropriately simulate the natural system (Fienen and others, 2010). The parameter value and its imposed lower and upper bounds represent the reasonable range of values or the precalibration (prior) uncertainty of the model input, based on prior knowledge of the system (Fienen and others, 2010). Excluding a model input from the adjustable parameter set, either as a “fixed” parameter or excluding it entirely from the PEST framework, equates to “hard coding” the model input and is similar to saying that model input is known with 100-percent certainty, which is highly unlikely in complex environmental models. However, the model may be insensitive to some model inputs, which cause far smaller effects on model outputs than other inputs. Although some parameters may be insensitive to the calibration process because they may not gain information from the calibration dataset, they can still cause changes in scenario results, which is an important consideration when specifying parameters. Additionally, some model inputs can be well constrained (less uncertainty) because they may have ample data to specify an acceptable value that produces acceptable model results, whereas others may be poorly constrained (more uncertainty) and require more freedom during calibration to find optimal values.

Calibration parameters were combined into groups defined as aquifer storage properties (aqprops), streamflow routing properties (sfrprops), farm or landscape properties (farmprops), pilot points for Kh (pilotpts), climate properties (climprops), and pilot points for aquifer anisotropy (apilotpts; table 8). As will be described in the “Two-Phase Calibration Approach” section, certain parameters were fixed to improve model stability, prevent spurious parameter values, and better simulate ET. After several calibration attempts, 305 of the 740 parameters were converted to nonadjustable or “fixed” parameters within the PEST framework so the final calibration used 435 adjustable parameters. The groups farmprops, aqprops, and climprops were all fixed for the final calibration, whereas the calibration treated the sfrprops, pilotpts, and apilotpts as adjustable.

Adjustable Parameters for Final Calibration

The streamflow routing properties (sfrprops) group of parameters included 11 adjustable parameters for vertical hydraulic conductivity of streambed (Ksb). Ksb parameters were defined for the lower and upper portion of each SFR reach and “tied” in PEST so that the Ksb remained the same for each SFR segment and did not change for reaches within a segment. Unique Ksb parameters were defined for natural streams according to soil class used in the CPIHM and for each canal simulated with the SFR (soil classification in fig. 1B; parameters in fig. 13). Initial values for calibration of Ksb for natural streams were set to less than 2 ft/d and canal Ksb values were set between 1 and 10 ft/d based on manual trial-and-error testing prior to automated calibration through which model stability was improved. The lower and upper bound were set to 0.001 and 15 ft/d, respectively, which is the acceptable range for the study area based on values from Peterson and others (2016).

Vertical hydraulic conductivity and Kh for each model cell was specified using interpolation between pilot points (see Doherty and others, 2010b). Vertical hydraulic conductivities were set using anisotropy ratios, also interpolated between pilot points. The pilot points for Kh (pilotpts) were specified as adjustable parameters to adequately represent the spatial heterogeneity of Kh across each model layer without needing to adjust Kh for each active model cell. A total of 124 pilot points were used: 49 for layer 1, 48 for layer 2, and 27 for layer 3 (fig. 13). In this study, pilot points were generated using the FloPy python module (Bakker and others, 2016), interpolation factors were generated using the “ppk2fac” PEST Groundwater Data Utility (Doherty 2018), and model input arrays generated from interpolation between the pilot points used the PyEmu module (White and others, 2016). The pilot point distribution was based on the presence of an active extent of each model layer. Pilot point density was chosen as one pilot point for every 25 active model cells. Pilot points within a few cells of a stream or general head boundary were manually moved away from the boundary to avoid the influence of boundary conditions. For this study, the vertical interval estimates of Kh from Houston and others (2013) were assigned to aquifer intervals that aligned with the HUs from Cannia and others (2006) and that correspond to the model layers (table 1). Initial values for individual pilot points were selected by extracting Kh values from two-dimensional arrays generated through interpolation of test hole data in Houston and others (2013) and summarized in table 1 (table 1 is a summary of the test hole data, not the interpolated grids used to extract the values).

The pilot points for aquifer anisotropy (apilotpts) were specified for the same locations as the Kh pilot points (fig. 13). Initial values for model layer 1 anisotropy pilot points were specified as 10:1 ratio (Kh:Kv) in accord with published values in Domenico and Schwartz (1990) for unconsolidated sands, silts, and clays of alluvial aquifers (Cannia and others, 2006). Layer 2 anisotropy at pilot points was specified as 20:1 because of the finer grained and lower Kh values of silts reported for that layer (Cannia and others, 2006). Layer 3 anisotropy at pilot points were specified as 10:1 ratio for sandstones (Domenico and Schwartz, 1990; Houston and others, 2013). Horizontal anisotropy was not simulated in the CPIHM.

Calibration Targets

Calibration targets are measured or estimated values assumed to represent conditions of the hydrologic system at the time of their measurement, with a degree of uncertainty that is based principally on measurement or estimation error. Calibration targets can be single measurements in time or a time series that represents the temporal behavior of the hydrologic system. PEST allows calibration targets to be clustered into groups where the contribution of each group to Φ can be assessed. The CPIHM had the following observation groups specified in PEST (table 9): development period average streamflows per stress period processed from daily average measured streamflows from USGS National Water Information System (NWIS) streamgages (nwisflow), development period average streamflows per stress period processed from daily average measured streamflows from the NeDNR streamgages (dnrflow), average streamflow differences between streamgages (sflowdiff), predevelopment model average streamflows per stress period processed from daily average measured streamflows from the USGS NWIS streamgages (ssflow), Predevelopment model average streamflow differences per stress period processed from daily average measured streamflows from the USGS NWIS streamgages (ssflowdiff), predevelopment model average groundwater levels processed from measured pre-1950 groundwater-levels in NWIS (sswls), development period measured groundwater-levels from USGS NWIS (realwls), development period differences in measured groundwater levels between stress periods (realdif), and development period measured groundwater-levels across transects of the Platte River, collected by the CPNRD (realwlst). Two sets of observation points were included that did not involve actual measurements but were used in the analysis of the calibrated model: scenario period “potential” observations that can be used for dataworth analysis (potheads), and scenario period “future” observations used in the predictive uncertainty analysis, located at the geometric center of each Groundwater Management Area (futheads) that were incorporated for the predictive uncertainty analyses (table 9). The calibration was not affected by any observation targets in the potheads or futheads groups because they were zero weighted. The CPIHM contained 40,711 weighted calibration targets.

Groundwater Levels

The CPIHM was calibrated to groundwater levels for the predevelopment and development models. The predevelopment model was calibrated to 189 estimated groundwater levels at 189 observation wells as the average of measured groundwater levels prior to January 1, 1950, and assigned to the sswls observation group (table 9; fig. 13). The development model was calibrated to 19,662 measured groundwater-level observations at 963 observation wells obtained from the NWIS database (U.S. Geological Survey, 2017) and 9,848 groundwater-level differences, which were second-order observations calculated as the difference between a measured groundwater level at a single location from one stress period to the next, to ensure the model matched changes in groundwater levels in addition to the absolute groundwater level (table 9). A total of 16,840 groundwater-level observations at 725 wells completed in layer 1; 1,065 groundwater-level observations at 125 wells completed in layer 2; and 1,759 groundwater-level observations at 113 wells completed in layer 3 (table 10). Although wells were distributed across the model area, a total of 859 wells had less than 10 observations and 68 wells had 100 or more observations (fig. 13). Western regions of the active model area in Custer, Frontier, and Lincoln Counties had very few observations (fig. 13). About 95 percent of the groundwater-level observations were measured after 1980 (fig. 14); therefore, the calibration incurred a bias toward the post-1980 period. This bias was desirable because one objective of the study was to use the model results to update the CPNRD’s GMP, which compared current management and MADs to 1982 groundwater levels. Consequently, the natural bias built into the model calibration (because most observations were measured after 1980) allowed PEST to focus on adjusting parameters to match observations in the post-1980 period. All groundwater-level calibration targets and locations are available in the model archive associated with this report (Traylor, 2023).

Two-Phase Calibration Approach

Calibration of the CPIHM involved two phases: a manual adjustment of parameters, followed by the automated calibration using BeoPEST. The two-phase calibration approach was adopted after several attempts to calibrate the CPIHM in a fully automated fashion because many of the adjustable parameters produced unrealistic results and resulted in significant model instability. Landscape parameters and some groundwater parameters were fixed to improve model stability during calibration, prevent spurious calibration results, and better match the conceptual estimates of landscape ET. The conceptual landscape ET estimates were based on measured data or calibrated outputs of previously developed models, which provided a blueprint for acceptable ranges of parameter values and outputs of ET and groundwater irrigation, for the CPIHM.

An assessment of parameter sensitivities to observations using a matrix of sensitivities, known as a Jacobian matrix, from a preliminary automated calibration attempt with all parameters set as adjustable, showed that the growing season ET_ref scale factors and Sy multipliers were some of the most sensitive (table 1.3). Preliminary automated calibration runs revealed that PEST preferred to adjust the Sy multipliers beyond acceptable ranges, even with heavy regularization weights applied, and as a result the multipliers were converted to fixed parameters before subsequent calibration. Further, the storage parameters were fixed in the calibration of layers 2–5 for the Eastern Model Unit groundwater-flow model from Peterson (2009) and the single layer Northern High Plains aquifer model from Peterson and others (2016); both models included the model domain for this study and used the same input test hole data derived from Houston and others (2013). This assessment of parameter sensitivities and methods from previous studies led to the implementation of the two-phase calibration approach that included the fixing of landscape and aquifer storage parameters to reasonable values based on available data.

In the initial manual adjustment phase, parameters were adjusted within the PEST framework to improve the fit between calibration targets and their simulated equivalent values. This phase focused on adjusting the growing season ET_ref scale factors to improve the CPIHM’s simulation of landscape ET and irrigation pumping throughout the development period model to values that were similar to the conceptual ET values and irrigation pumping from table 2. Model outputs were postprocessed to produce annual values of those ET and irrigation fluxes, and the recent 5 years of the simulated outputs (2011–16) were compared to the conceptual values from table 3. As mentioned in the “Fixed Parameters for the Final Calibration” section of this report, the ET_ref scale factors required adjustment from their 1.0 values that were standard for the other scale manually adjusted factors. The unadjusted average annual inputs of ET_ref shown in figure 10B were about 57 inches, which was about 7 to 8 percent less than the measured values from Irmak and Skaggs (2011). Consequently, initial manual trial and error testing indicated that simulated total AET (which is directly controlled by the ET_ref, and therefore the ET_ref scale factors) for all crop types was too low compared to the conceptual model water-budget values. Additionally, the undersimulated AET was responsible for the undersimulation of groundwater withdrawals for irrigation compared to available data and conceptual model values described in the “Conceptual Model of the Hydrologic System” section of this report, because lower ET corresponded to lower CIR. Initial automated calibration runs, with ET_ref scale factor parameters set as adjustable, showed that they were highly sensitive parameters as PEST preferentially adjusted them and they hit their upper bounds when the upper limit was set to even an unrealistic 3.0 (3.0 equates to a 3 times multiplier on ET_ref). Trial-and-error manual adjustment was used to determine the best initial values of each ET_ref scale factor to combat overfitting and maintain agreement with the conceptual model landscape ET.

After manual calibration, when the simulated recent water-budget fluxes were in accord with the conceptual model values, the final set of input parameter values obtained from the manual calibration were set as the initial parameter values for the automated PEST calibration. Therefore, initial parameter values at the beginning of the automated calibration phase were derived from a combination of local knowledge of the study area, information accrued from other scientific studies of the study area or similar hydrologic settings, and information gained about the CPIHM during the manual calibration phase. To honor these initial parameter values and their information and to keep PEST from deviating too far from the initial values, Tikhonov regularization was applied during the automated calibration, also described in “Regularization of Parameters” section of this report (Doherty and Hunt, 2010; Doherty, 201531).

Comparison of Calibration Targets to Simulated Equivalent Values

The comparison of calibration targets to their simulated equivalent values is an assessment of the fit between measured or observed data such as groundwater levels and streamflows with corresponding simulated values produced by the CPIHM. During the calibration process, PEST compared the calibration targets to corresponding simulated values and calculated the difference between the observation targets and simulated equivalent values, which for the remainder of the report will be referred to as the “residual.” A positive residual indicates an underestimation of simulated values by the CPIHM and a negative residual indicates an overestimation of simulated values by the CPIHM. The Gage package extracted the simulated streamflows from SFR reaches coincident with the location of the 53 streamgages for comparison to the streamflow calibration targets. The mod2obs PEST utility was used to extract the corresponding simulated groundwater levels for comparison to the calibration target groundwater levels (Doherty, 2018).

Groundwater Levels

The calibrated predevelopment CPIHM groundwater levels displayed an adequate fit to the calibration targets for the sswls observation group with an average absolute residual of 8.2 ft (table 11). Therefore, the predevelopment CPIHM provided acceptable initial conditions for the transient development-period CPIHM. Spatial trends in the calibrated sswls included an overestimation north of the Platte River and an underestimation south and east of the Platte River and minimal bias along the Platte River and Wood River from west to east (fig. 15). The calibrated development CPIHM groundwater levels displayed an adequate fit to the calibration targets for the realwls and realwlst observation groups, with average absolute residuals of 9.6 and 2.8 ft, respectively (figs. 15, 16A, and 2.3; table 11). The combined realwls and realwlst average absolute groundwater-level residuals for model layers 1, 2, and 3 were 6.1, 12.4, and 7.4 ft, respectively (table 11) for the whole development-phase model, which indicated that the CPIHM adequately simulated the groundwater levels for each layer. The fit between the measured and simulated groundwater levels in the realwls and realwlst groups was somewhat better within the CPNRD area after May 1, 1980, because the calibration of the CPIHM focused on the minimization of groundwater-level residuals for this area and period. This improvement in fit occurred mostly in layers 1 and 2 for the CPIHM and within the CPNRD; layer 3 exhibited an absolute residual that was larger (7.6 ft for the CPIHM and 7.5 ft for the CPNRD) after May 1, 1980, compared to the absolute residual for all time periods (7.4 ft for the CPIHM and 7.3 ft for the CPNRD); however, the difference was about 0.1 ft (table 12). The combined realwls and realwlst CPNRD average absolute groundwater-level residuals after May 1, 1980, for model layers 1, 2, and 3 were 4.0, 10.7, and 7.5 ft, respectively (table 12). Further, the calibrated development CPIHM groundwater levels displayed an adequate fit to the calibration targets for the realdif observation group with average absolute residuals of 1.2 ft and a standard deviation of 2.8 ft (table 11). These residuals indicated that the development-period CPIHM captured the transient nature of the groundwater-subsystem as groundwater-levels changed.

Table 11.    

Calibration results statistics for the groundwater-level and streamflow observation groups of the Central Platte Integrated Hydrologic Model.

[<, less than; >, greater than; --, could not calculate a value]

Table 11.    Calibration results statistics for the groundwater-level and streamflow observation groups of the Central Platte Integrated Hydrologic Model.

Observation group	Average	Standard deviation	Minimum	25 percent	50 percent	75 percent	Maximum	Average absolute residual (feet)
Calibrated groundwater-level residual statistics, in feet
realwls	−1.4	14.2	−95.9	−5.9	0.4	5.9	137.0	9.6
realwlst	0.1	4.1	−20.6	−1.7	−0.2	2.0	21.2	2.8
sswls	2.6	11.7	−29.4	−4.1	0.7	7.2	56.8	8.2
realdif	<0.1	2.8	−80.3	−0.4	>−0.1	0.5	43.1	1.2
Calibrated streamflow statistics, in cubic feet per second
dnrflow	15.1	1,060.7	−15,012.1	−53.7	−1.0	28.0	11,831.9	356.4
nwisflow	−191.8	833.0	−15,226.4	−191.4	−10.8	21.5	7,890.3	464.6
ssflow	2,230.0	1,518.2	1,156.4	1,693.2	2,230.0	2,766.7	3,303.5	1,073.5
sflowdiff	171.0	1,088.0	−4,230.5	−268.3	76.7	506.5	5,845.0	695.2
ssflowdif	−2,147.0	--	−2,147.0	−2,147.0	−2,147.0	−2,147.0	−2,147.0	0.0

Table 12.    

Average simulated groundwater levels and their residuals by layer in the Central Platte Integrated Hydrologic Model domain and the Central Platte Natural Resources District domain for the entire development model period (May 1, 1895, to December 31, 2016) and post-May 1, 1980.

[CPIHM, Central Platte Integrated Hydrologic Model; CPNRD, Central; Platte Natural Resources District]

Table 12.    Average simulated groundwater levels and their residuals by layer in the Central Platte Integrated Hydrologic Model domain and the Central Platte Natural Resources District domain for the entire development model period (May 1, 1895, to December 31, 2016) and post-May 1, 1980.

Model layer	Count		CPIHM domain			CPIHM domain post-May 1, 1980		CPNRD			CPNRD domain post-May 1, 1980
	Observ-ations	Observation wells	Average simulated groundwater-level, in feet above mean sea level	Average residual, in feet	Average absolute residual, in feet	Average residual, in feet	Average absolute residual, in feet	Average simulated groundwater-level in the CPNRD, in feet above mean sea level	Average residual, in feet	Average absolute residual, in feet	Average residual, in feet	Average absolute residual, in feet
1	16,840	725	2,052.06	−0.7	6.1	−0.4	6.0	2,010.06	2.3	4.0	2.3	4.0
2	1,065	125	2,013.42	−2.7	12.4	−3.6	11.5	1,959.01	6.9	10.3	7.4	10.7
3	1,759	113	2,211.60	1.5	7.4	1.4	7.6	2,182.39	1.1	7.3	1.1	7.5

The calibrated groundwater-level residuals after May 1, 1980, exhibited spatial bias throughout the CPIHM domain. The average groundwater-level residuals at observation wells south of the Platte River were mostly negative, which indicated that the CPIHM overestimated the groundwater levels in this region; some residuals in this region exhibited a large bias of greater than 25 ft (fig. 15), including well 159 (fig. 17A). The average groundwater-level residuals at observation wells north of the Platte River were mostly positive, which indicated that the CPIHM underestimated the groundwater levels in this region (fig. 15), including well 352 (fig. 17B). Some of these wells were within the CPNRD boundary, particularly in Dawson, Buffalo, Hall, Howard, and Nance counties (fig. 15). The residuals along the Platte River were generally less than 5 ft (figs. 15 and 17C–F). Most of the residuals for observation wells that were overestimated by the CPIHM were located to the south and east of the Platte River outside the CPNRD boundary, particularly in Frontier, Gosper, Phelps, Kearney, Adams, Clay, Hamilton, and York Counties; because of their locations outside the CPNRD, they were weighted lower than the observations inside the CPNRD boundary, which caused them to have less effect on the Φ during the calibration. Therefore, PEST did not try to improve the fit for calibration targets outside the CPNRD boundary as much as the targets inside the CPNRD. Greater variability in groundwater levels outside of the Platte River Valley topographic region was also expected because of more dynamic terrain and depth to water in the Valley-side Slopes, Dissected Plains, and Sandhills topographic regions compared to the Valley region of the Platte River (fig. 1B). Additional comparisons of measured and simulated groundwater altitudes at other locations used in the model calibration process are available in the USGS data release associated with this report (Traylor, 2023).

Streamflows

The calibrated predevelopment CPIHM undersimulated streamflow by 2,230 ft³/s for the ssflow observation group (table 11). However, given that (1) the ssflow observations were more uncertain than measured streamflows because they were estimated based on measured flows between 1930 and 1940, (2) the simulated predevelopment groundwater levels exhibited an adequate fit to the sswls observations, and (3) streamflow is a more dynamic observation than groundwater levels, and (4) major inflows to the Platte River are a transient input for the development CPIHM, the marginal fit to the ssflow observations still resulted in an adequate calibration for the development-period model. Most of the calibrated development-period CPIHM streamflows displayed a passable fit to the calibration targets (fig. 16E–H). Calibration targets from the nwisflow group exhibited a coefficient of determination of 0.776 with their simulated equivalent values. The nwisflow group exhibited a minor positive bias in their residuals, which indicated that the CPIHM oversimulated streamflows. The CPIHM accurately simulated the flows along the Platte River, including the outflow point of the Platte River near Duncan, Nebr. (fig. 18), streamgage (USGS streamgage 06774000; fig. 1A), which was important for the efficacy of the water balance during the simulation. The largest bias in flow along the Platte River occurred at splits in the complex braided network of channels such as a north and south channel at the Platte River near Overton North Channel, Nebr. (USGS streamgage 06767998) and the Platte River near Overton South Channel, Nebr. (USGS streamgage 06767999): residuals were −961 and −108 ft³/s, respectively (fig. 15). However, the single channel Platte River near Overton, Nebr. (USGS streamgage 06768000) residual was only −31 ft³/s (fig. 15). These residuals were expected because the model was not designed to simulate all the local complexities of a braided channel system such as Platte River; rather it was designed to simulate the primary channels of the Platte River and other streams in the study area. Further, the dnrflow residuals for the development-period CPIHM exhibited a large overestimation of flows on the Platte River streamgages, which was expected because these observations were weighted low in the calibration and consequently the parameters in PEST were not as sensitive to these observations compared to the nwisflow group. Additional comparisons of measured and simulated streamflow at 53 streamgages used in the model calibration process are available in the USGS data release associated with this report (Traylor, 2023).

The final phi distribution by observation group is presented in figure 19. The observation group with the largest contribution to Φ was from the nwisflow group (26.2 percent), followed by the realdif (24.2 percent) and realwls (19.9 percent) groups (fig. 19). The primary calibration focused on reducing the residuals of the realwls, realdif, and nwisflow groups more than other groups; therefore, the final Φ distribution, where the top three were nwisflow, realdif, and realwls, reflects the purpose of the study objectives.

Calibrated Parameters

This section of the report describes the final calibrated values for the automated values of the 435 adjustable parameters in the six parameter groups for the CPIHM. The following sections will provide a summary of the final calibrated parameter values and a comparison to the initial parameter values to illustrate how the calibration changed the parameter values to improve the fit between the calibration targets and their simulated equivalent values. The summaries will focus on the hydraulic parameters contained in the core MODFLOW–NWT functionality of the CPIHM because the landscape and climate parameters contained in the FMP functionality of the CPIHM were fixed before the automated calibration process. A sensitivity analysis also summarizes which parameters had the largest effect on the model outputs, and parameter identifiability indicates how well the calibration target dataset informed the parameters, which is discussed in the next section of the report.

Calibrated Groundwater-Flow Parameters

Calibrated Kh estimated at pilot points (fig. 13) and interpolated to the model grid kept the general trends and magnitude of the initially estimated hydraulic conductivity described in Houston and others (2013). The calibrated average Kh for layers 1, 2, and 3 were about 70, 32, and 35 ft/d, respectively (table 13). These numbers are all slightly higher (within 5 ft/d) than the initial starting values (table 1). Calibrated layer 1 Kh was generally less than 50 ft/d in the western region and more than 50 ft/d in the eastern region of the CPIHM (fig. 20A). Calibrated layer 2 Kh was generally less than 50 ft/d in almost all regions of the CPIHM except for the northeast corner where some Kh values were greater than 200 ft/d (fig. 20B); calibrated values greater than 200 ft/d were also simulated in the coincident COHYST model region in Peterson (2009). Calibrated layer 3 Kh was generally less than 40 ft/d in almost all regions of the CPIHM except for a west-central region where some Kh values were greater than 100 ft/d (fig. 20C). The highest Kh values for each layer were distributed around pilot points that had high initial values and were not outside the range of reasonable values for that layer; the high Kh values agreed with the test hole data from Houston and others (2013).

Calibrated anisotropy, the ratio of horizontal to vertical hydraulic conductivity (Kh:Kv) estimated at pilot points (fig. 13) and interpolated to the model grid remained in the general range of acceptable values for each layer. The average vertical anisotropy for layers 1, 2, and 3 were 9.73, 9.86, and 9.99, respectively, which were equated to vertical hydraulic conductivity values for all layers (between 0.33 and 31.70 ft/d), and the calibrated average vertical hydraulic conductivity for layers 1, 2, and 3 were 7.86, 3.41, and 3.60 ft/d, respectively (table 13). Calibrated layer 1 vertical hydraulic conductivity was between less than 7 ft/d in most of the western CPIHM domain and primarily between 7 and 20 ft/d in the southern and eastern domain (fig. 20D). Layer 2 vertical hydraulic conductivity was less than 10 ft/d for most of the CPIHM domain (fig. 20E). Calibrated layer 2 vertical hydraulic conductivity increased from the initial values in much of the central and southern CPIHM domain where it was present and generally decreased in the northern and eastern regions (fig. 20E). Calibrated layer 3 vertical hydraulic conductivity remained near initial values for most of the CPIHM domain except for two areas in the west where values were greater than about 10 ft/d (fig. 20F). The Sy, which was fixed during calibration, is shown in figure 20G–I.

Parameter Sensitivity and Identifiability

The automated calibration process includes the calculation of the sensitivity of simulated observations to model parameters, where sensitivity is the change in simulated observation divided by the change in a model parameter; PEST records sensitivities of simulated observations to each parameter in the Jacobian matrix (Doherty, 2005). Parameter sensitivities were extracted from the Jacobian matrix file written by PEST and composite parameter sensitivities were analyzed from the sensitivity file, also written by PEST, and are available in the model archive (Traylor, 2023). Parameter sensitivities can provide some insight into parameters that were most impactful during the calibration process. The parameter sensitivities were assessed by parameter group and the most sensitive group was the anisotropy pilot points (apilotpts) with an average composite sensitivity of 0.351. The Kh pilot points (pilotpts) group had a slightly smaller average sensitivity of 0.342. The streambed hydraulic conductivity group (sfrprops) was the least sensitive parameter group. The apilopts group and, to a slightly less extent, pilotpts were the most sensitive because the anisotropy parameters affected the magnitude of vertical flow in the development model. The present gradient between the layers where average simulated groundwater levels were lowest in layer 2 (table 12) coupled with the large amount of groundwater irrigation wells that pumped in the development model (more than 20,000 by 1980; fig. 3), affected the vertical flow dynamics and the groundwater levels in the transient system.

Parameter sensitivities for the final Jacobian matrix used for calibration, with fixed landscape parameters, were low compared to an assessment of parameter sensitivity with all parameters set as adjustable (fig. 21A and table 1.3). In the Jacobian matrix, calculated with all parameters set to adjustable, July and August development period ET_ref scale factor parameters were the most sensitive and, during preliminary calibration runs, exhibited large deviations from initial values that affected AET, deep percolation (recharge to the water table), and irrigation pumping (table 1.3). The Sy multiplier for layer 1 was also one of the most sensitive parameters (12 of 740) that deviated from reasonable values during preliminary calibration runs (table 1.3). The parameter sensitivities under the fixed scheme influenced the forecast uncertainties and are described in the “Forecast Uncertainty Analysis” section of this report.

Parameter sensitivities are one metric that can be useful in understanding the parameter behavior during calibration. In complex environmental models such as the CPIHM, parameters are often correlated because there are far fewer observations than parameters, which can limit the ability for PEST to fully resolve some or all parameters with a unique solution during calibration. Quantification of how PEST is able to resolve unique parameter values is known as “parameter identifiability.” Parameter identifiability was assessed and provided insight into the information gained by each parameter from the calibration targets dataset during the calibration process and how well each parameter was uniquely estimated and correlated with other parameters (Doherty and Hunt, 2009; Doherty, 2015). Parameter identifiability was analyzed using pyEMU, an open-source Python module that offers a suite of tools to interrogate PEST files and analyze uncertainty (White and others, 2016). Identifiability values ranged from zero to one, where a value of one represented a parameter that could be uniquely resolved based on the information in the calibration target dataset; any error associated with the parameter was a result of noise or uncertainty in the calibration targets. Conversely, an identifiability of zero represented a parameter that gained no information from the calibration target dataset and did not cause a change to the residuals that were calculated at the calibration targets. Identifiability values between zero and one indicate that the parameter gained some information from the calibration targets during the calibration, but it also shares some information gained by the calibration target dataset with other parameters and cannot be uniquely solved (Doherty, 2015).

Overall, the pilot points for Kh parameters were the most identifiable as a group given the spatial distribution of the pilot points throughout the model domain and the number of observations in the calibration targets dataset. Individually, the most identifiable parameters in the CPIHM were the hcond_a2 (SFR vertical hydraulic conductivity of the streambed for soil type A), l2hkpp78 (layer 2 Kh at pilot point 78), and l3hkpp23 (layer 3 Kh at pilot point 23; fig. 21B), which were located downgradient of many groundwater-level observations that informed the groundwater-flow direction and reaction to stresses throughout the simulation period. Parameter hcond_a2 exhibited a low sensitivity of about 0.74 and was not in the top 50 most sensitive parameters shown in figure 21A. Whereas it could be the most uniquely solved in the calibration (highest identifiability, shown in figure 21B), the low sensitivity indicated that it did not have a substantial impact on the calibration. The high identifiability of parameter hcond_a2 is likely because it was one of the parameters that spanned much of the model domain and could gain information from observations at many of the 53 streamgage observation locations also distributed across the model domain. Conversely, parameter l2hkpp78 exhibited a high sensitivity of about 1.54 (second most sensitive parameter), which, coupled with its high identifiability, indicated that it was more uniquely solvable and impactful to the calibration than other parameters except hcond_a2. Parameter l3hkpp23 exhibited a very low sensitivity of about 0.10 (332nd most sensitive parameters out of the 435 parameters); however, its identifiability indicated that it was more uniquely solvable than most parameters, but it was not particularly impactful to the calibration (fig. 21A–B).

Landscape Water Budgets for the Central Platte Integrated Hydrological Model Domain

The calibrated predevelopment CPIHM produced landscape water-budget inflows and outflows that were meant to approximate the average conditions prior to surface-water and groundwater irrigation development and to provide reasonable initial conditions for the development CPIHM. Precipitation was the only inflow to the landscape subsystem with an average annual flux of 10,937,662 acre-ft/yr (26.6 inches). Total outflows were distributed among Ep, Tp, Ei, Ti, deep percolation, and runoff. The primary outflow was Ep with an average annual flux of about 6,362,644 acre-ft/yr (58 percent of the outflows). The sum of Ep and Tp, evapotranspiration of precipitation (ETp), accounted for 97 percent of the outflows from the landscape subsystem with an average annual flux of 10,577,561 acre-ft/yr (25.8 inches). The average annual flux for deep percolation and runoff were 91,687 acre-ft/yr (0.22 inches) and 268,414 acre-ft/yr (0.65 inches), respectively (Traylor, 2023).

The calibrated development CPIHM produced landscape water-budget inflows and outflows that maintained the general trends and magnitudes similar to the conceptual understanding of the landscape water subsystem (fig. 22; table 14). Total inflows to the landscape, prior to groundwater irrigation development, were largely from precipitation, with minor amounts from surface-water deliveries for irrigation. The average annual flux of precipitation was 9,978,276 acre-ft/yr (24.3 inches if spread equally across the model domain), or about 93 percent of the inflows to the landscape for the development period (May 1, 1895, to December 31, 2016) (table 14). The average annual total outflows from the landscape during the development period were 10,709,293 acre-ft/yr and the primary outflows were to Tp (4,184,162 acre-ft/yr) and Ep (3,751,102 acre-ft/yr); ETp accounted for about 80 percent of the outflows from the landscape at 7,935,263 acre-ft/yr (table 14 shows Ep and Tp). Compared to the estimates for evapotranspiration of irrigation plus precipitation for the CPIHM in table 2, which ranged from 8,600,000 to 11,000,000 acre-ft/yr, the sum of all evaporation and transpiration columns in table 14 total 8,554,930 acre-ft/yr. Deep percolation (recharge to the water table) and runoff were also substantial outflows from the landscape across the development period (1895–2016) with average annual fluxes of 1,122,257 (about 2.7 inches) and 1,032,106 (about 2.5 inches) acre-ft/yr, respectively (table 14). Further, the deep percolation values on irrigated land and dry cropland simulated by the CPIHM were in accord with recharge measurements from Steele and others (2014) and the overall average annual deep percolation of about 2.7 inches is within the ranges simulated by the models in Peterson (2009) and Peterson and others (2016) for this study area.

Temporal trends in the landscape budget were marked by the onset of widespread groundwater irrigation that was responsible for the increase in total inflows to the landscape during the development period; average annual total inflows for 1895–2016 and 2011–16 were 10,709,293 and 12,084,641 acre-ft/yr (table 14). The increase in total inflows led to an increase in total outflows primarily through an increase in deep percolation and ET of irrigation water (sum of Ei and Ti) (fig. 22). The ETp declined after about 1960; however, total landscape-derived ET, the ETp and irrigation water (sum of Ep, Ei, Tp, and Ti), exhibited less than a 2-percent increase in the recent development period (2011 to 2016) compared to the pre-1940 development period (fig. 23). The landscape water budgets for the development CPIHM were similar from 1895 to about 1960 whereby the budget was primarily affected by inflows from annual precipitation and outflows to ET of that precipitation. After about 1960, the increase in groundwater irrigation for crops added additional outflows as ET of the irrigation water (fig. 22, 23). The large increase in irrigation after 1960 was the primary cause of the increase in deep percolation from about 660,000 acre-ft/yr (1895–1960 average) to about 1,700,000 acre-ft/yr (1961–2016 average; fig. 22). Another factor that contributed to the substantial increase in deep percolation was a 7-percent increase in average precipitation after 1960.

In the last 5 years of the development model simulation (2011–16), the average annual volume of total landscape ET, which included ET of precipitation and irrigation water sources, was 8,420,099 acre-ft/yr, about 11 percent less than the conceptual model estimate of 9,500,000 acre-ft/yr for the same time period but only 4.5 percent less than the estimated range of uncertainty for ET from the conceptual estimates in table 2 (table 14). The average deep percolation and runoff volumes simulated for the last 5 years of the development period (2011–16) were 2,056,936 and 1,607,607 acre-ft/yr, respectively. Monthly trends from 1980–2016 indicated that the largest rates of deep percolation occurred in the spring (April and May) when precipitation was highest, and ETp was relatively low compared to peak growing season (July) ET (fig. 24).

The manually calibrated development period average annual ET_ref increased from 57.4 to 74.5 inches, a 30-percent increase from the precalibrated ET_ref (fig. 25). The reduction in ET_ref after 1980 is a result of different scale factors applied to the bi-annual and monthly stress periods to account for monthly dynamics in ET; however, this reduction was not present in the annual simulated total landscape ET, which increased by less than 2 percent for the development period owing to an increase in ET from irrigation water, as stated previously in the discussion of temporal trends (figs. 22, 23). The average monthly ET_ref increased 13-percent from the precalibrated ET_ref, for May through September (fig. 26, 10B). The increase in ET_ref was necessary because it led to an increase in the CWD for all crop types and increased the CIR for irrigated crops to values that were similar to conceptual estimates of total landscape ET and groundwater irrigation pumping.

Although some simulated landscape flow components such as irrigation pumping and ET were a little different than conceptual estimates, they still represented the conceptual model within expected uncertainty. The landscape budget simulated by the CPIHM can be considered a more precise representation of the landscape water budget than the conceptual flux estimates from table 2, because the conceptual estimates are uncertain, but also because the simulated budgets are calibrated and reconciled with simulated groundwater flow.

Landscape Water Budgets for Groundwater Management Areas and Other Domains

Simulated landscape water budgets were assessed for each GWMA and the CPNRD (table 14). The budgets for each GWMA and the CPNRD were similar to the CPIHM budget and were characterized by primary inflows from precipitation and large outflows to ET. Deep percolation was affected by three factors: precipitation; the amount of irrigation water from irrigation wells or surface-water deliveries; and ETg, which was linked to the depth of the water table. Like the CPIHM domain, GWMAs with large irrigation volumes per acre, such as GWMA 3 and 14 (fig. 9), also exhibited large volumes of recharge per acre (table 14, table 3.1). Similarly, GWMAs with large fluxes of precipitation (in the eastern region of the CPIHM) and shallow water tables with large fluxes of ETg, such as GWMAs 19, 22, and 23, also exhibited the large recharge volumes per acre (table 14, table 3.1). GWMAs 1, 6, and 7 exhibited the smallest deep percolation (recharge) volumes per acre because they were in the western region of the CPIHM with less precipitation and less irrigation (fig. 9; table 14, table 3.1). Larger amounts of deep percolation occurred on irrigated lands compared to dryland in the CPNRD (fig. 27). The average annual rates of deep percolation on irrigated land after 1980 (about 5.1 inches) were double the amount on nonirrigated land (about 2.4 inches), and both were within the ranges measured in Steele and others (2014). Further, the relation among precipitation, irrigation, and deep percolation was evident during years of extreme climate conditions. For example, deep percolation on irrigated and nonirrigated cropland was greater than the amount of irrigation pumpage in 1993 because it was the wettest year between 1981 and 2016 (fig. 27). As a result, the CIR was minimal and there was surplus precipitation that went to deep percolation. Conversely, 2002 and 2012 exhibited large irrigation pumpage and little recharge because these were severe drought years. CIR for the severe drought years was very high and there was minimal precipitation available for recharge (fig. 27).

Groundwater-Flow Budgets for the Central Platte Integrated Hydrologic Model Domain

The calibrated predevelopment CPIHM produced groundwater-budget inflows and outflows that were meant to approximate the average conditions prior to surface-water and groundwater irrigation development and to provide reasonable initial conditions for the development CPIHM. Total inflows and outflows were about 193,669,000 and 194,105,766 acre-ft/yr, respectively. Approximate steady-state conditions were achieved whereby the change in groundwater storage across the 500,000,000-day transient stress period was −572 acre-ft/yr. Stream leakage was the primary inflow to the groundwater subsystem, with a flux of 144,885,246 acre-ft/yr (75 percent of total inflows). Inflows to the study area across model boundaries accounted for 18 percent of total inflows and recharge from deep percolation accounted for 7 percent of total inflows. The primary outflow was to ETg, which constituted 53 percent of total outflows. Discharge of groundwater to streams as base flow constituted 33 percent of total outflows, and outflows across model boundaries accounted for 14 percent of total outflows (Traylor, 2023).

The calibrated development CPIHM produced groundwater-flow budget inflows and outflows that reflected the transient conditions simulated in the development period model (table 15, 16). The largest inflow component was recharge (analogous to deep percolation from the landscape), with an average development period (May 1, 1895, to December 31, 2016) annual volume of 1,122,257 acre-ft/yr (table 15). Stream leakage (602,318 acre-ft/yr) was another major annual inflow, which indicated that the stream system in the CPIHM was an important contributor to the groundwater-flow system (table 15). The largest outflows were to irrigation wells (693,171 acre-ft/yr). Changes in groundwater storage were substantial during the development period; the replenishment of groundwater storage was −985,912 acre-ft/yr and releases from groundwater storage were 863,519 acre-ft/yr. The large releases from storage (groundwater-level declines) began to make a substantial impact on the groundwater-flow budget after about 1960 owing to the increase in groundwater pumping for irrigation (fig. 28A). The net storage changes were driven primarily by the relation between outflows to irrigation wells and farm net recharge each year; farm net recharge is the difference between recharge from deep percolation and ETg (fig. 28B). Generally, groundwater levels increased when outflows to irrigation wells were less than farm net recharge (fig. 28B). The replenishment of storage came from a combination of groundwater flowing into the study area, recharge (deep percolation), and stream leakage (fig. 28C, shown as positive values of net streams). Prior to widespread irrigation development (about 1940), recharge was the primary source of water to replenish storage for most average and wet years. During dry periods, stream leakage became the dominant source of water to replenish storage (fig. 28C). After 1980, the average annual releases from groundwater storage were similar to those from recharge at about 1,957,000 acre-ft/yr and 1,941,000 acre-ft/yr, respectively. On average, the increase in outflows to irrigation wells were generally supported by increases in recharge throughout the development period except in dry years (for example, 2002, 2012) when farm net recharge was much less than outflows to irrigation wells. After 1980, even in dry years inflows from farm net recharge were generally greater than net streams (fig. 28D). Net replenishment to storage (negative net storage; groundwater-level rise) generally occurred during wet years or when the absolute magnitude of farm net recharge was greater than the absolute magnitude of irrigation wells (fig. 28B). Also, after 1980, the average annual depth (where depth represents a water-budget component expressed in terms of volume per area) of groundwater pumped by irrigation wells was 10.8 inches (1,698,000 acre-ft/yr), which was similar to the conceptual estimated depth of 10 inches (1,800,000 acre-ft/yr) from table 3 (fig. 27).

The monthly trend in the simulated groundwater-flow budgets was typical of systems with large amounts of groundwater-irrigated cropland and characterized by large irrigation demands during the growing season. Monthly trends were assessed for the 1981 to 2016 period when the development model stress periods lengths changed to monthly (table 6). Recharge was the largest average monthly inflow from October to June (fig. 29). The largest volumes of recharge occurred in April, May, and October when groundwater pumping to irrigation wells was minimal or absent and precipitation was relatively high (fig. 29). The highest average monthly precipitation in the CPIHM was in June, but there was less recharge compared to April and May because the high landscape ET (1,306,656 acre-ft/yr) removed much of the precipitation that would have otherwise become recharge to the water table. The lowest volumes of recharge occurred in driest months of December, January, and February (fig. 29). The largest inflows to the groundwater-flow system during the 2 months with the highest irrigation (July and August) were caused by the release of groundwater from storage that supported the groundwater pumping during that time; the replenishment to storage occurred primarily during April, May, and October when recharge was highest (fig. 29). The average volume of water pumped for irrigation during July (814,815 acre-ft/yr) was more than twice as large as any other outflows for any month (replenishment to groundwater storage in May: −428,771 acre-ft/yr) (fig. 29).

Groundwater-Flow Budgets for the Groundwater Management Areas and Other Domains

Individual simulated groundwater-flow budgets were assessed for each GWMA and the CPNRD for the development CPIHM. The budgets for each area were similar to the CPIHM budget, which is characterized by the largest inflows of recharge, largest outflows to irrigation wells, and large fluxes to and from groundwater storage (table 15). Other major inflows to many GWMAs and the CPNRD were the groundwater flows from adjacent zones upgradient and stream leakage (table 15). The spatial orientation of the CPNRD and the 24 GWMAs within the CPNRD along the Platte River corridor, perpendicular to the regional west to east direction of groundwater flow, caused most of the GWMAs to exhibit substantial inflows of groundwater from adjacent zones upgradient and outflows to adjacent downgradient zones (table 15). In the CPNRD, prior to 1960, the primary inflow was from streams and after 1960 the primary inflow was farm net recharge (recharge from deep percolation minus ETg) (fig. 28D). Like the CPIHM, the relation between outflows to irrigation wells and farm net recharge was the primary driver of net changes to storage within the CPNRD later in the development period; net replenishment to storage (groundwater-level rise) generally occurred when the difference between the absolute value of irrigation wells and farm net recharge was small (less than 270,000 acre-ft/yr). These results indicated that the groundwater was under more stress when irrigation pumping was high or recharge was low (fig. 28E). Net flow to and from the CPNRD (net zone flow) was a net outflow for 106 of the 122 years simulated in the development period and indicated that the loss of groundwater within the CPNRD may be attributed to a combination ETg, base flow, and outflow to irrigation wells (fig. 28D).

For the recent development period (2011–16), groundwater wells were the largest outflow for 21 of the 24 GWMAs, indicating that irrigation pumping is the dominant stress on most regions of the CPNRD. GWMAs 5, 19, and 23 were dominated by outflows of net zone flow, farm net recharge, and net streams, respectively, rather than irrigation pumping (fig. 28F, table 3.6). On average, changes in storage were driven by outflows to irrigation wells and farm net recharge for most GWMAs (fig. 28F, outflows are negative net values). Generally, GWMAs that exhibited net release to storage (groundwater-level decline) from 2011 to 2016 had larger imbalances between the magnitude of irrigation and farm net recharge, which indicated that recharge increases owing to pumping in these areas were less substantial and not enough to support those levels of pumping during that time period (fig. 28F, G). Further, the GWMAs that had smaller outflows to irrigation wells, coupled with low amounts of recharge, but exhibited net releases from storage (groundwater-level declines) generally had larger net outflows to adjacent zones, which indicated that net zone flow was an important contributor to the stability of the groundwater-flow system for these GWMAs (for example, GWMAs 1, 5, 6, 7; fig. 28G). Inflows from adjacent zones were larger than recharge (deep percolation) in 15 of the 24 GWMAs; however, net zone flow was less than zero (a net outflow) for 14 of 24 GWMAs, which indicated that more groundwater leaves a zone than enters (fig. 28F, table 15). For example, GWMAs 20 and 21 exhibited net releases from storage (positive net storage) despite having more net inflows from farm net recharge than irrigation pumping because they had substantial net outflows to adjacent zones (as negative net zones) (fig. 28F). The source of inflows to a GWMA was important in understanding the stability of the groundwater-flow system. The stability of the groundwater-flow system for each GWMA is influenced by activities and stresses within that area, such as pumping and recharge, but also the activities in GWMAs upgradient. Declines in an upgradient GWMA can cause less inflow of groundwater to the adjacent downgradient GWMAs even if the activities and stresses in the downgradient GWMA remain constant. If a GWMA receives more inflows from adjacent zones than other inflows, the groundwater-flow system of the GWMA is more influenced on activities upgradient. GWMAs 4, 8, 21, 23, and 24 each have greater than double the volume of inflows from adjacent zones compared to recharge (deep percolation) and are examples of areas that are more influenced by activities upgradient (table 15).