SUSTAINABILITY OF THE HIGH PLAINS AQUIFER : FROM DELIBERATE IMPACTS TO UNINTENDED CONSEQUENCES By Erin Marie King Haacker A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Environmental Geosciences Doctor of Philosophy 201 8 ABSTRACT SUSTAINABILITY OF THE HIGH PLAINS AQUIFER : FROM DELIBERATE IMPACTS TO UNINTENDED CONSEQUENCES By Erin Marie King Haacker The High Plains Aquifer (HPA) is a vital water resource for the agricultural production of the United States. Most of the aquifer is being used at an unsustainable rate, creating an irrigation economy that relies on resource exhaustion. In the very long term, much of the aquifer will not be usable for irrigated agricul ture, but this leaves a senescence look like? What are the primary drivers of decline, and how are patterns of decline changing? What are the implications for land use change and management? This project uses two approaches to address these questions. The first is statistical: through spatial statistics and re gression analysis, I investigate the patterns and relationships in the water table data. The second approach is process - based modeling, which uses the physical framework of the water cycle and the landscape to examine the distribution of water in the aquif er. These approaches are highly complementary. Statistical methods offer a n empirical estimate of the water table in space and time. These methods suggest relationships between parameters used in modeling, such as the lag between precipitation and water table change, that might be difficult to incorporate a priori into a process - based model due to the difficulty in observing the relevant processes. Process - based models are grounded in the physics and geography of the system, and so they are able to accommodate scenarios that do not strictly follow the historical course of events . This counterfactual modeling is helpful in teasing apart the roles of the many processes that make up the water cycle, and in understanding the likely results of different management scenarios. Both process - based and statistical treatments demonstrate th at the aquifer is going dry in many places. The Southern High Plains (SHP) region in Texas and New Mexico is at greatest risk . In th e SHP , recharge is low, and groundwater is virtually the sole source for irrigation . The rate of groundwater decline is link ed with the climate and with crop prices, but these factors only account for 25 - 40 % of the total rate of depletion or recovery across the aquifer. In the SHP, conversion from dryland farming to perennial grassland, as through the Conservation Reserve Prog ram (CRP), does not confer protection to the water table. Reductions in irrigation water can increase the lifespan of the aquifer . M odel results indicate that in the most vulnerable part of the SHP, a 50% reduction in irrigation could have slowed the incre ase in the depleted area of the HPA by a factor of six between 2001 and 2014, from about 3.3% to about 0.6% of land area drying out in that timeframe. This holds whether the decrease is at the intensive margin (less irrigation per acre) or extensive margin ( fewer acres irrigated). Meanwhile, with no irrigation during that 15 - year period, the water table would have recovered by about 0.93 m, enough saturated thickness to provide water for about two irrigation seasons . iv This effort is dedicated t o five great teachers and mentors who set me on this path : Gail Wintczak, Denise Senor, Dave Trexler, Marc Hendrix, and Bill Woessner. Thank you. v ACKNOWLEDGEMENTS Throughout this degree program, I received support from the National Science Foundation through NSF WSC Award 1039180, and from a University Distinguished Fellowship awarded by the Michigan State University Graduate School. The office staff at the Environm ental Science and Policy Program and the Department of Earth and Environmental Sciences do the hard work that keep s everything running: Marcy Heberer, Karessa Weir, Jackie Bennet, Heidi Lynd, Pam Robinson, Brittany Walter , Judi Smelser, and Dallas Coryell . A mong my graduate school cohort, I would particularly like to thank my fellow musketeers, Autumn Parish and Alex Kuhl . My family is also terrifically supportive (and/or still very tiny ). I would like to thank my committee, six people who are Very Busy an d Important, for always finding time to offer valuable feedback. I n the last leg of my research, I accepted a research position in the Nebraska Water Center at the University of Nebraska - Lincoln . Thank you to Chittaranjan Ray and the Ogallala Water Coordin ated Agriculture Project (USDA NIFA Award 2015 - 09812) for the support and additional insights that were integral in finishing my degree program at Michigan State. 8 7. An y mistakes and opinions originate from the author and do not reflect an official position of any institution. vi PREFACE If the High Plains Aquifer were a lake, it would have been the 7th largest in the world by volume - t, following decades of decline, it finally would have lost 7th place to Lake Huron, dipping precipitously from about 3560 km 3 to 3535 km 3 of water in storage . The High Plains stretches across parts of eight states, in an area the size of Sweden. Its wate r resources are distributed very unevenly. The greatest volume of water, nearly 80% of the total water in the aquifer, is in the Northern High Plains, primarily in Nebraska. The hotter southern parts of the aquifer, lacking large surface water features, re ly heavily on the remaining 20%. Agriculture is the basis for the economy, and irrigation water is the basis for agriculture across most of the region (Smidt et al. 2015; Fenichel et al. 2016). Primary crops include corn, soybeans, winter wheat, sorghum, c otton, and alfalfa. It is likely that nearly everyone in the United States has eaten, worn, or burned (as ethanol) products fed by High Plains Aquifer water. The HPA exists due to geologic happenstance. It is an unconfined aquifer whose matrix was washed across an uneven surface of bedrock, so that the surface is much smoother than the aquifer base. This means that the paleovalleys and uplands were a major determinant of saturated thickness prior to development of the resource. The distribution of sediment thickness, combined with weather and irrigation patterns, has led some areas to decline in saturated thickness more rapidly than other areas. Historically, the aquifer was in equilibrium from the time its sediments were vii deposited (10 - 12 million years ago ; Gutentag et al. 1984) until the first irrigation wells were drilled at the end of the 19th century. Early settlers to the High Plains used windmills to pump enough water for a few acres of vegetables, and watered cattle at springs, without much impact on water resources. Eventually, with improvements in well technology and access to credit, demand for irrigation water increased (Green stewardship began to catch up to resource use, and groundwater management areas were initiated across the High Plains. However, many of these districts were created to serve the needs of the irrigators. Extending the life of the aquifer is not always seen as As a res ult of its importance - it is the largest freshwater aquifer in the United States (Dennehy 2000) - the HPA has been studied extensively through the years. Regional studies such as the Cooperative Hydrology Study (COHYST) in the Platte River Corridor, the G roundwater Availability Models (GAMs) in Texas, and modeling of the Republican River Basin, have added to extensive United States Geological Survey (USGS) reports starting with the High Plains Regional Aquifer Study. As computers, data availability, and mo deling technologies have improved, there are now unprecedented opportunities to create complex, linked models, both statistical and physical, to further characterize the aquifer in space and time. This dissertation relies almost equally on statistical and process - based descriptions of the HPA. These approaches are very different, but complement each other. Spatial statistics are mostly empirical; they are used to look for relationships, viii correlations, and patterns, without many assumptions of cause and effec t relationships . Process - based, deterministic models can be used to test scenarios, since the physical relationships between the inputs are built into the models. Modeling methods can ainst the kind of observation datasets that can be built and tested with statistics. Empirical and process - based methods are not always juxtaposed, and each approach can be used to inform the other. In this dissertation, statistically - derived water table elevations are used to describe changes through time, but also as a comparison dataset for groundwater modeling. Model boundaries are informed by saturated thickness estimates, and estimates of change in water table elevation through time. The process - base d models incorporate a large number of inputs that affect the water cycle, while the statistical models look at just one or a few datasets for empirical relationships. This creates two lines of evidence for results - not completely independent, but methodo logically distinct, to allow for a stronger sense of the system. For example, when I first began designing the groundwater model presented in Chapters 3 and 4, my model boundaries extended to the escarpments of the Southern High Plains. As I worked on the initial calibration for the groundwater model, I kept running into problems with certain areas going dry. A literature search indicated a wide range of hydraulic conductivity estimates for the region, which disagreed by two to three orders of magnitude. F inally, as I was examining the statistical estimates for saturated thickness, I realized that the same areas that were going dry in the process - based model were in places that I had statistically determined to have no appreciable ix saturated thickness above the bedrock. This caused me to radically rethink my groundwater model boundaries. Ultimately, I shrank my groundwater model domain to include only the most important irrigated subregion of the SHP, with no - flow boundaries where the aquifer thins out. This led almost immediately to a very stable model with steady - state conditions that were similar to predevelopment estimates. Statistical and process - based models each work in their own way , but both methods build our understanding in the same direction, a little closer to the reality of the HPA. This dissertation can be thought of as being divided into two section s, each of which includes two chapters . The first section is a series of stat istical treatments of the aquifer through time. It is guided by well observations, which are the most direct measure of the changes that have occurred since the development of High Plains water resources. The statistics mostly describe the existing data, a nd therefore allow the aquifer to tell its story fairly directly. The relative simplicity of the analyses, and the need for as many observations as possible, encourages the study of the aquifer as a whole. The price for this directness and scope is that th e analysis does not account for natural processes of the water cycle. estimates of water table elevation . This builds on work from the USGS, such as the High Plains Regional Aquifer - System Analysis (RASA) (Gutentag et al. 1984), Luckey et al. (1981), and McGuire (2003). Through the use of sophisticated spatial statistics, this x study allows more data to be in corporated for the water table estimates, for an annual time series that stretches back to 1935. This time series is used to forecast declines using cell - by - cell linear regression, providing an estimate for the pumping horizon if trends continue. A version of this study was published in the journal Groundwater , and is cited in other chapters as Haacker et al. (201 6 ). The second study further investigates the pattern of water table change shown by well hydrographs. It uses a statistical technique called segmented regression. Segmented regression is a form of linear regression, except that two lines are fit to a data ser ies, on either side of a breakpoint - an inflection where the pattern changes. For well hydrographs, this might mean that the water table begins to recover, or starts declining faster. The timing and direction of breakpoints are then compared to the implem entation of water management areas, to show that these management districts have a complex interaction with the pattern of water table decline. For most groundwater management areas, their initiation coincides with mostly negative breakpoints, when the aqu ifer starts declining faster; five to ten years later, there is a rise in positive breakpoints, indicating recovery, or at least slower depletion. However, cause and effect are difficult to sort out with purely statistical, rather than process - based, inves tigations. The next section of this dissertation focuses on process - based modeling of the Southern High Plains in parts of New Mexico and Texas. An integrated surface water and groundwater water model was created using the Landscape Hydrology Model (LHM). The model area extends across the southern part of the aquifer, up to the xi Canadian River Valley. Within the broader LHM domain, fully saturated groundwater flow is simulated for the most heavily irrigated part of this region. This Southern High Plains mo del was used to answer two sets of questions, one focusing on the surface water component, and the other on groundwater. The first uses spatial Conservation Reserve Program (CRP) data to examine the changes that might occur from land cover transition, in t his case from agriculture to shortgrass prairie. Grassland has a higher ratio of transpiration to evaporation, and less runoff than agricultural land. As a result, dryland farms in the Southern High Plains have notably higher recharge rates than the unbrok en prairie (Stanton et al. 2011; Scanlon et al. 2012). When land is transitioned back from dryland farming to mixed grasses, this change in recharge is reversed, suggesting that CRP may have the greatest positive effect for water resources in areas where a griculture has led to rising water tables and salinization. In the fourth study, four irrigation scenarios are run to estimate the relative effects of irrigation reduction at the intensive (less water applied per acre) and extensive ( fewer acres watered) margins. It is shown that limitations similar to those in intensive management areas, such as the Kansas Local Enhanced Management Area, decrease the rate of decline substantially, but are not enough to make that part of the aquifer sustainable. A reductio n in irrigated acreage, without retiring land from agriculture, can slow groundwater decline, but even a complete cessation of irrigation does not cause a rapid recovery in water tables. xii TABLE OF CONTENTS LIST OF TABLES ................................ ................................ ................................ ....................... xiv LIST OF FIGURES ................................ ................................ ................................ ...................... xv KEY TO ABBREVIATIONS ................................ ................................ ................................ ...... xx Chapter 1 The High Plains Aquifer From Predevelopment to Senescence ...................... 1 Abstract ................................ ................................ ................................ ......................... 1 Introduction and Background ................................ ................................ .................... 1 Methods ................................ ................................ ................................ ......................... 5 The High Plains Aquifer ................................ ................................ ......................... 5 Water Level Interpolation ................................ ................................ ....................... 6 Characterizing Water Levels and the Development Timeframe ...................... 9 Annual Water Level Changes ................................ ................................ .............. 15 Water Level Decline Forecast ................................ ................................ ............... 17 Results and Discussion ................................ ................................ .............................. 19 Predev elopment Map and Aquifer Storage ................................ ....................... 19 Annual Water Level and Storage Changes ................................ ........................ 19 Recent Rates of Decline and Depletion Forecast ................................ ............... 28 Conclusions ................................ ................................ ................................ ................. 30 Chapter 2 Effects of Management Areas, Drought, and Commodity Prices on Groundwater Decline Patterns across the H igh Plains Aquifer ................................ ......... 33 Abstract ................................ ................................ ................................ ....................... 33 Introduction and Background ................................ ................................ .................. 34 Groundwater Depletion and Management ................................ ........................ 34 Study Site ................................ ................................ ................................ ................. 36 Groundwater Management Areas ................................ ................................ ....... 40 Factors Affecting Irrigation Demand ................................ ................................ .. 42 Methods ................................ ................................ ................................ ....................... 43 Analysis Objectives and Procedure ................................ ................................ ..... 43 Spatial and Temporal Datasets ................................ ................................ ............ 43 Well Hydrographs ................................ ................................ ................................ . 45 Breakpoint Analysis ................................ ................................ ............................... 46 Cross - Correlation and Multiple Linear Regression Analysis .......................... 49 Results and Discussion ................................ ................................ .............................. 52 Regional Patterns ................................ ................................ ................................ ... 52 Effects of Groundwater Management Areas ................................ ...................... 54 Saturated Thickness ................................ ................................ ............................... 59 xiii Cross - Correlation and Multiple Linear Regression Results ............................ 60 Conclusions ................................ ................................ ................................ ................. 66 Chapter 3 Effec ts of Land Use Change on Water Resources of the Southern High Plains ................................ ................................ ................................ ................................ ....................... 70 Abstract ................................ ................................ ................................ ....................... 70 Introduction and Background ................................ ................................ .................. 71 Study Area ................................ ................................ ................................ ............... 71 The Hydrological L andscape ................................ ................................ ............... 72 Land Retirement Programs ................................ ................................ ................... 73 Purpose of Study ................................ ................................ ................................ .... 74 Methods ................................ ................................ ................................ ....................... 75 Models and Input Data ................................ ................................ ......................... 75 Scenar ios ................................ ................................ ................................ .................. 80 Results and Discussion ................................ ................................ .............................. 83 Baseline Model Results ................................ ................................ .......................... 83 CRP Transition Scenarios ................................ ................................ ...................... 85 Conclusions ................................ ................................ ................................ ................. 86 Chapter 4 Groundwater Decline and Water Resource Management in the Texas Panhandle ................................ ................................ ................................ ................................ .... 89 Abstract ................................ ................................ ................................ ....................... 89 Introduction and Background ................................ ................................ .................. 90 Study Area ................................ ................................ ................................ ............... 90 Water for Irrigati on in the Southern High Plains ................................ .............. 91 Methods ................................ ................................ ................................ ....................... 94 Model Description ................................ ................................ ................................ . 94 Hydraulic Conductivity Calibration ................................ ................................ ... 98 Irrigation Scenarios ................................ ................................ .............................. 100 Results and Discussion ................................ ................................ ............................ 103 Water Tables and Saturated Thickness ................................ ............................. 103 Irrigable Area ................................ ................................ ................................ ........ 106 Conclusions ................................ ................................ ................................ ............... 108 REFERENCES ................................ ................................ ................................ ........................... 112 xiv LIST OF TABLES Table 1.1. A comparison of storage changes from recent publications. ............................. 25 Table 2.1. Factors in the linear regression models shown in Figure 11. Not all cross - correlates shown in Figure 10 were significant in the multiple linear regression models. negative numbers indicate that the variable is measured prior to the annual break fraction, e.g. a lag of - the coefficient of the variable with regards to the net annual break fraction; positive estimates mean that high values of the variable are associated with positive breakpoints, and the relationship between the variable and the annual break fraction is due to chance. NHP: Northern High Plains; CHP: Central High Plains; SHP: Southern High Plains. ... 51 Table 3.1. Data sources for LHM. ................................ ................................ ............................ 76 Table 3.2. CRP land was primarily identified as grassland in the NASS Cropland Data Layer (CDL) for the available years in the study area (Myneni et al. 2014). ..................... 81 Table 4.1. Data sources for the MODFLOW groundwater model. ................................ ..... 96 xv LIST OF FIGURES Figure 1.1. Map showing the location of the High Plains Aquif er and its regions, with U.S. state boundaries and postal abbreviations. Major rivers of the High Plains are labeled in blue. Aquifer boundary from the study by Qi (2010). Projection is North ongitude. ................... 5 Figure 1.2. Graph showing the number of measurements used in water table interpolation, and th e number of measurements that were removed for quality control purposes. ................................ ................................ ................................ ................................ ....... 8 recorded for wells in the predevelopment dataset for (A) the USGS and (B) this project. ................................ ................................ ................................ ................................ ....................... 11 Figure 1.4. (A) The Northern High Plains, shown with well measurements (orange) and points from the National Hydrography Dataset (NHD) (blue). T hese points were used to create a predevelopment map. (B, C) A conceptual diagram showing the effect of NHD point inclusion. (B) Without NHD points, the interpolation technique ignores valleys where groundwater intersects the land surface, potentially over estimating groundwater levels. (C) With NHD points, the water level is constrained to the land surface. Note that NHD points were only used in areas where the initial water level interpolation was above the land surface; thus this technique does not cause the ............... 1 2 Figure 1.5. Specific yield estimates from McGuire et al. 2012. ................................ ............ 17 Figure 1.6. The water level elevation (piezometric surface) during predevelopment, in meters. ................................ ................................ ................................ ................................ .......... 18 Figure 1.7. (A) Water level change, from the predevelopment su rface to 2012. (B) Hydrographs of individual wells from eachof the regions of the High Plains (USGS site IDs: SHP, 341420101441602; CHP, 362746102364102; NHP, 404345098560001). (C) Average water level change across the three regions of the High Plains. ......................... 20 Figure 1.8. Mean standard deviation of water table elevation krigi ng estimates. ........... 21 Figure 1.9. The percent remaining saturated thic kness in the SHP in 2012 compared to 1943 (the year of maximum saturated thickness for the SHP). Note that areas shown in dark gray are locations where the saturated thickness was interpolated to be zero at xvi predevelopment. ................................ ................................ ................................ ........................ 22 Figure 1.10. (A) The total volume of water in the High Plains Aquifer from 1935 to 2012. (B) Water remaining in storage from 1935 to 2012, as a percentage of the maximum water in storage. ................................ ................................ ................................ ......................... 24 from 1992 to 2012. ................................ ................................ ................................ ...................... 27 Figure 1.12. Percent of aquifer area with less than 9 m saturated thickness in future decades, assuming a linear pattern of decline. ................................ ................................ ...... 28 Figure 1.13. Projected depletion timeframe for the entire High Plains Aquifer, using a 9 m threshold saturat ed thickness for aquifer utility and assuming a linear pattern of decline. Gray areas are in little danger of depletion in the foreseeable future. ................ 29 Figure 2.1. State - level groundwater doctrines for the High Plains (Joshi 2005). The three regions of the High Plains Aquifer Northern, Central, and Southern a re overlayed in gray. Some states use a modified system as described in Smidt et al. (2016). .................. 37 Figure 2.2. (A) Mean annual precipitation in millimeters, from 1930 - 2015 (PRISM); (B) mean annual temperature in degrees Celsius, from 1930 - 2015 (PRISM); (C) groundwater level decline from predevelopment through 2016, estimate d using methods from Haacker et al. (2016); (D) irrigated/non - irrigated area from MODIS - MIrAD satellite data, 2002, 2007, and 2012. Groundwater management areas with governance below state level are outlined in black. ................................ ................................ ............................... 38 Figure 2.3. Establishment year for groundwater management areas that are predominantly on the High P lains. Select establishment years are shown on the map for reference. States are labeled with postal abbreviations. ................................ ....................... 39 Figure 2.4. Example well hydrographs with segmented regressions for well locations displayed in the map (right). GMA establishment dates are noted for each example well. R 2 values for segmented regressi ons are reported for wells A and C. Wells B and D did not have statistically significant breaks in slope; R 2 for those wells are shown for linear regression. GMA = Groundwater Management Area. Water table elevations are shown in meters above sea level. ................................ ................................ ................................ ......... 48 Figure 2.5. Annual break fraction through time for the (A) Nor thern, (B) Central, and (C) Southern HPA. ................................ ................................ ................................ ..................... 53 Figure 2.6. Positive and negative annual break fractions over time in select Nebraska Natural Resources Districts (NRDs). NRD outlines are shown at bottom. ....................... 55 xvii Figure 2.7. Positive and negative annual break fractions over time in Groundwater Management Districts (GMDs), Kansas. GMD outlines are shown at bottom. ................ 56 Figure 2.8. Positive and negative annual break fractions over time in select Underground Water Conservation Districts (UWCDs), Texas. UWCD outlines are shown on the map at bottom. ................................ ................................ ................................ ... 58 Figure 2.9. Multiple linear regression models of net annual break fraction through time. Net annual break fractions are the sum of positive and negative breakpoints in a given year, divided by the number of measurements in that year. The gray shading indicates the disparity between modeled and observed net annual break fractions. Model adjusted R 2 values for the Northern, Central, and Southern High Plains Aquifer are 0.40, 0.25, and 0.25, respectively. ................................ ................................ ................................ ...... 62 Figure 2.10. Net number of breakpoints normalized to the number of wells measured in a yea r, cross - correlated with factors affecting irrigation water demand. Cross - correlation graphs are shown for the factors that are independently most strongly correlated with net break rate in each region. The most significant weather and commodity price vari years before the breakpoint occurs. A lag of zero is the correlation in the same year as the net annual break fraction. Dashed horizontal lines indicate the threshold for correlation to be sign ificant at p<0.05. National crop prices are in 2016 dollars. ............. 63 Fi gure 2.11. Time series of observed, expected, and residual net breakpoints with respect to the year of establishment for groundwater management areas. Year 0 is the year in which the management area was established; negative and positive numbers on the x - a xis represent the number of years pre - and post - establishment. A. The expected number of net breakpoints, based on crop prices and moisture conditions, is shown in light gray, compared with the observed net breakpoints in dark gray. B. The residual net br eakpoints the difference between expected and observed is shown in red (net negative residual) or blue (net positive residual). ................................ ................................ . 65 Figure 3.1. Location, showing Southern High Plains boundary cut off at the Canadian River, with county outlines. ................................ ................................ ................................ ..... 71 Figure 3.2. Playa lakes shown as points. Inset shows actual boundaries of playas. ........ 72 Figure 3.3. (A) Average enrollment in CRP for all counties located at least 85% within the H igh Plains, and counties in eight High Plains states, but with 15% or less area on the High Plains. (B) Map of percentage of land in CRP in the Southern High Plains in 2016. ................................ ................................ ................................ ................................ .............. 74 xviii Figure 3.4. Examples of LHM inputs. (A) Percent sand in the top layer of soil, derived from SSURGO data; note some discrepancies at county bound aries. Urban areas within model boundaries are blank. (B) Elevation in meters. The land surface of the High Plains slopes from west to east. (C) Irrigated area from MIrAD - US for 2012. .................. 77 Figure 3.5. Model boundaries. Light blue is the LHM model boundary; darker blue is the MODFLOW groundwater model, which overlaps the L HM surface model area. A 5x5 km block of 100 cells is shown for reference to model scale. ................................ ........ 78 Figure 3.6. Results for baseline analysis, averaged over the model period, 2001 - 2014. (A) Runo n to cells, in meters. Most of the model area is internally drained. (B) Recharge, i.e. percolation below the root zone. ................................ ................................ ............................. 80 Figure 3.7. Land use in the SHP in 2016, CDL. Land was reclassed in LHM as nine land types, including agriculture, grassland, shrubland, and barren. ................................ ........ 82 Figure 3.8. (A) Simulated versus observed groundwater levels for the baseline run, which did not incorporate information from the spatial CRP dataset. (B) The baseline run has groundwater elevation estimate residuals within a margin of error of observations (two standard deviations around the mean, shown). ................................ ... 84 Figure 3.9. Average deep percolation, cm, for the same CRP land area in three scenarios: baseline (CRP is grassland), mo st likely scenario if CRP were ended (1/3 of CRP area reassigned as agriculture), and all CRP land area converted to agriculture. Recharge is consistently higher in farmland. (A) Average deep percolation during the growing season (Jun - Sep), land transition scenarios minus baseline. (B) Average deep percolation during the non - growing season (Nov - Apr). ................................ ................................ .......... 85 Figure 4.1. Boundary of the Southern High Plains, showing groundwater conservation districts in the Texas Panhandle. The groundwater model overlies about half of the High Plains Underground Water Conservation District No. 1, which was instituted in 1951 and adopted its first Groundwater Management Plan in 1998. ................................ . 91 Figure 4.2. Thickness of Ogallala Formation sediments in the Southern High Plains. The black interior border is the extent of the groundwater model in this study. .................... 92 Figure 4.3. Irrigated area in the Southern High Plains in 2007, whic h is the coverage used as a baseline in the LHM - MODFLOW model. Most of the irrigated area overlies locations where the estimated saturated thickness in 2001 is greater than 9 m. The black interior border is the extent of the groundwater model in this study . .............................. 93 Figure 4.4. Closeup of the MODFLOW groundwater model boundar y, showing specified head cells (dark blue) and 50 meter contours for water table elevations in 2001, xix in meters above sea level. The High Plains Underground Water Conservation District No. 1 is shown in gray. ................................ ................................ ................................ .............. 94 Figure 4.5. (A) Saturated hydraulic conductivity from the Texas Groundwater Availability Model of the Ogallala Aquifer (Deeds et al. 2015). (B) Hydraulic conductivity used in the MODFLOW model for this study, modified from USGS estimates (Cederstrand and Becker 1998b). ................................ ................................ ........... 98 Figure 4.6. Observed versus simulated head for the steady - state model used to develop boundary conditions and hydraulic conductivity parameters for MODFLOW. .............. 99 Figure 4.7. Irrigation scenarios used in this study. Left: irrig ation extent in 2007, used in the baseline and per - hectare application limited scenarios. Right: closeups showing irrigated extent in the three irrigated scenarios. The fourth scenario had no irrigation. ................................ ................................ ................................ ................................ ..................... 101 Figure 4.8. Water table observations versus simulated water tables, 2001 - 2014, for the baseline scenario (no limit s to extent of irrigated land or intensity of water applied). Residuals are two standard deviations around the mean. ................................ ................. 103 Figure 4.9. (A) Saturated thickness for the baseline scenario, with full irrigation extent and intensity. (B) Saturated thickness in the no - irrigation scenario minus the baseline. (C) Saturated thickness when irrigated area is limited minus the baseline. (D) Satu rated thickness when irrigation intensity is limited minus the baseline. Saturated thickness was greater than the baseline in all irrigation limited scenarios, but the spatial distribution of saturated thickness was different depending on the density of ir rigation in a locality. All saturated thicknesses are for 2014, at the end of the model runs. ....... 104 Figure 4.10. Mean monthly irrigation water applied by scenario. ................................ .... 105 Figure 4.11. (A) Head difference for the limited - extent scenario and the limited - intensity scenario. Blue indicates that the limited - extent scenario outperformed the limited - intensity scenario, while red indicates that the limited - intensity scenario performed better in terms of water table elevation at the end of the model period. (B) Average water applied per year in the baseline (unrestricted irrigation) scenario. ...................... 105 Figure 4.12. Irrigable area over time, based on a 9 m threshold for saturated thickness. ................................ ................................ ................................ ................................ ..................... 107 xx KEY TO ABBREVIATIONS CDL Cropland Data Layer CHP Central High Plains CLU Common Land Units COHYST Cooperative Hydrology Study CREP Conservation Reserve Enhancement Program CRP Conservation Reserve Program EPA Environmental Protection Agency ET Evapotranspiration DNR Department of Natural Resources GAM Groundwater Availability Model GMA Groundwater Management Area gpm gallons per minute GWCD Ground Water Conse rvation District HPA High Plains Aquifer hr - hour km - kilometers LEMA Local Enhanced Management Area LHM Landscape Hydrology Model m - meters MIrAD MODIS Irrigated Agriculture Dataset for the United States xxi MODIS Moderate Resolution Imaging Spectroradiometer MODFLOW Modular Three - Dimensional Finite - Difference Groundwater Flow Model NASS National Agricultural Statistics Service NED National Elevation Dataset NHD National Hydrography Dataset NHP Northern High Plains NIFA National In stitute of Food and Agriculture NRD Natural Resources District NSF National Science Foundation PDSI Palmer Drought Severity Index PRISM - Parameter Elevation Regression on Independent Slopes Model RASA Regional Aquifer - System Analysis SHP Souther n High Plains SSURGO Soil Survey Geographic Database USDA United States Department of Agriculture USGS United States Geological Survey UWCD Underground Water Conservation Districts WSC Water Sustainability and Climate yr year 1 Chapter 1 The High Plains Aquifer From Predevelopment to Senescence A version of t his chapter appeared in the journal Groundwater ( 54 : 2, pp. 231 - 242). Abstract A large imbalance between recharge and water withdrawal has caused vital regions of the High Plains Aquifer (HPA) to experience significant declines in storage. A new predevelopment map coupled with a synthesis of annual water levels demonstrates that aquifer storage has declined by approximately 410 km 3 since the us estimates. If current rates of decline continue, much of the Southern High Plains and parts of the Central High Plains will have insufficient water for irrigation within the next twenty to thirty years, whereas most of the Northern High Plains will expe rience little change in storage. In the western parts of the Central and northern part of the Southern High Plains, saturated thickness has locally declined by more than 50%, and is currently declining at rates of 10 - 20% of initial thickness per decade. Th e most agriculturally productive portions of the High Plains will not support irrigated production within a matter of decades without significant changes in management. Introduction and Background The High Plains Aquifer supplies nearly a third of all grou ndwater for irrigation 2 in the United States (Dennehy et al. 2002). Groundwater pumping for irrigation on the High Plains creates a large deficit between discharge and recharge (Butler et al. 2013; Scanlon et al. 2012). Water level decline across the High P lains Aquifer causes impacts to natural ecosystems, such as reduced stream flow (Scanlon et al. 2012). These declines also have significant economic impacts. According to an analysis by Torell et al. (1990), depth to water is one of the primary determinant s of the price of irrigable farmland on the High Plains, though reductions in saturated thickness have a greater impact on pumping rates (Hendricks and Peterson 2012). Water table declines have already led some to shift from irrigated to dryland farming, a s the saturated thickness of the aquifer no longer supports pumping in some areas (Terrell et al. 2002). Degradation of the aquifer began long before the issue of water sustainability was addressed, though problems were recognized in some areas a few decad es following development (White et al. 1946; Gaum 1953). Stewardship of the remaining available water in the aquifer will depend on a thorough understanding of its history to provide a sound basis for future projections. While the United States Geological Survey (USGS) has published a number of state - of - the - aquifer reports (e.g. Luckey et al. 1981; McGuire 2004, 2009, 2012), these do not contain comprehensive annual estimates of water level changes. This paper extends previous work by the USGS to provide gr eater coverage of water level change in space and time, using a new method of estimation. A spatiotemporal archive of water levels is necessary to understand how the aquifer has been affected by the last century of resource use. 3 A thorough assessment of w ater level change requires understanding the state of the aquifer prior to development. While the predevelopment condition of the aquifer development (Bredehoeft 2002), it does give a baseline from which to compare the effects of prolonged impact. For example, when Luckey et al. (1981) published a map of water level changes since predevelopment, they suggested that development occurred beginning in the 1930 s in Texas, grad ually progressing northward during subsequent of 1980. McGuire (2004) introduced the simplifying assumption that development had quently been considered as a development date for the entire aquifer (e.g., Stanton et al. 2011; Scanlon et al. 2012; McGuire 2012). While 1950 does predate widespread installation of high - capacity pumps (Luckey and Becker 1999), making it a good candidate for a generalized aquifer - wide predevelopment date, such an assumption can mask important local - to regional - scale changes resulting from earlier groundwater pumping, as indicated in McGuire (2012). Aquifer management can be guided by projections of futur e resource depletion. To address current unsustainable water withdrawals, groundwater management districts in Texas and Kansas , and Natural Resources Districts in Nebraska, are working to introduce conservation and management goals to increase the lifespan of the aquifer. goal that seeks to maintain at least 50% of the saturated thickness as of Jan. 1, 2010 by 4 2050 (High Plains UWCD 2014 ). Additionally, in Kansas, a Local Enha nced Management Area (LEMA) is being used to emplace conservation programs that can locally increase the usable lifespan of portions of the aquifer (Kansas State Senate 2012; Whittemore et al. 2014). Efforts in these two areas may serve as a model for futu re conservation measures applied more broadly across the HPA. Meeting conservation targets requires both quantification of current water level decline trajectories and projections of future decline. Scanlon et al. (2012) suggest an overall depletion tim eframe of 630 years for the whole aquifer, including the areas in the NHP that are not being depleted at a significant rate. They also conclude that the lifespan of the aquifer varies greatly among and within subregions of the High Plains. Wilson (2007), i n a report to the Kansas Legislature, projected that much of the Central High Plains in western Kansas would be unable to provide water at a sufficient rate for irrigated agriculture within 25 years. Any plan for the future of the aquifer must take into a ccount its past and present. A history of High Plains water levels is presented here that demonstrates the spatial and temporal complexity of changes using annual water level maps. I present an approach to identify and map the changes in water level that h ave followed a complex history of aquifer development. Finally, I extend this historical perspective with a spatially explicit projection of depletion that can help guide management of the aquifer at local to regional scales. 5 Methods The High Plains Aquifer The High Plains is a physiographic region in the central United States, which can be subdivided (Stanton et al. 2011) into the N orthern High Plains (NHP), Central High Plains (CHP), and Southern High Plains (SHP; Figure 1 . 1 ). The significance of the High Plains Aquifer to agriculture in the Un ited States has fostered an extensive literature, including a review of the physical and hydrogeological aspects of the High Plains (Gutentag et al. 1984). The aquifer was formed by the eroded Figure 1 . 1 . Map showing the location of the High Plains Aquifer and its regions, with U.S. state bound aries and postal abbreviations. Major rivers of the High Plains are labeled in blue. Aquifer boundary from the study by Qi (2010). Projection is North America Albers Equal Area ° longitude. 6 detritus of the growing Rocky Mountains over the course of the Miocene and early Pliocene epochs (24 to 2 million years ago), which generated a large a ssemblage of alluvial fans that were subsequently cut off from the sediment source by south - flowing streams. A large proportion of the aquifer was also deposited as wi ndblown sediment. This changing environment of deposition led to variations in hydraulic conductivity in the modern aquifer. Most of the regional aquifer is unconfined, with discontinuous clay, silt, sand, and gravel layers that range from highly consolida ted (Gutentag et al. 1984) to unconsolidated. The water - bearing units, particularly in Nebraska, include low hydraulic conductivity layers (Dlubac et al. 2013); portions of the Republican River Basin have thus been modeled previously as a leaky confined aq uifer ( X. Chen et al. 2005). Groundwater generally flows from west to east, following the trend in the topography of the High Plains land surface. Water Level Interpolation Water levels are measured at distinct locations, as depth to water where a well in tersects the saturated zone or where groundwater intersects the land surface (streams). The water levels must be interpolated between these locations. Here I apply wi th explicit estimates of error (Isaaks and Srivastava 1989). The interpolated variable was change in water level, as opposed to water level elevation or depth to water. Interpolating changes in water level, rather than interpolating water levels directly, has the advantage that changes in areas far from data points tend to revert to a number 7 close to zero. While the absence of a monitoring well does not indicate a lack of water level change, this is a conservative assumption, and most areas of profound chan ge have early and frequent well measurements. This minimizes the effects of areas where kriging error might overwhelm the signal from water level changes, as in areas where very little water withdrawal occurs. In previous studies (e.g., McGuire 2012), wat er level change has been estimated by subtracting measurements taken in one year from measurements taken in the same well in a previous year, and interpolating this difference. I modify this method in several ways. First, I calculate the difference between the current year and prior two from the interpolated surfaces, allowing many more measurements to be included in an infrequently sampled dataset, since any well measure d in a given year can be used in years. An initial water level map had to serve as the basis for computing subsequent water level changes, which in this case is provided by the predevelopment water level map described below. I also employ a different interpolation technique and additional post - processing steps. Water level measurements from the High Plains Aquifer, spanning from 1876 to 2012, were obtained from the United St ates Geological Survey (USGS). These measurements were combined with static water levels from wells listed as points of diversion from the Nebraska Department of Natural Resources (DNR). For the interpolation period in this study, 1935 - 2012, total annual n umbers of measurements 8 used from these sources ranged between 1,972 and 21,013 ( Figure 1 . 2 ). These measurements were then filtered for quality control in a three step process. First, aquifer base ele vations and land surface elevations were assigned to the wells and used to remove measurements falling outside these limits (i.e., negative depth to water or water level below the aquifer base). Second, a temporal filter was applied, removing measurements collected during the active pumping season (defined as May through September inclusive). Third, measurements beyond three standard deviations from the mean for a given well were also removed, since a few measurements existed that were orders of magnitud e away from the mean of measurements from that well; these Figure 1 . 2 . Graph showing the number of measurements used in water table interpolation, and the number of measurements that were removed for quality control purposes. 9 outliers were attributed to measurement or recording errors rather than behavior of the groundwater system. All land surface elevations were derived from the National Elevation Dataset (NED) (Gesc h 2007; Gesch et al., 2002), at 1 arcsecond resolution. The National Hydrography Dataset (NHD) (USGS 2012 ) was used to describe stream locations. The bottom of the aquifer was interpolated from a contour map from Cederstrand and Becker (1998 a ) using the To po to Raster tool in ArcGIS. Rasterizing contoured datasets attempts to minimize terracing, but artifacts are still present in some areas of the aquifer. The base of the aqui fer is one of the boundaries that defines the saturated thickness of the aquifer. Thus, the accuracy of the aquifer base is an important constraint on the accuracy of the estimates, which is particularly significant given that saturated thickness largely d etermines the rate at which water can be pumped from an aquifer. All raster surfaces, including the aquifer base, water levels, water level change, and land surface elevation, were created or resampled to a grid at 250 m x 250 m resolution, to ensure unif ormity among datasets. Characterizing Water Levels and the Development Timeframe Water levels under predevelopment conditions, prior to significant anthropogenic effects, can be represented as a water storage balance in which mean recharge equals mean di scharge over years to decades. After development, pumping 10 from wells significantly increases net aquifer discharge. Recharge is altered by land use changes that impact evapotranspiration and runoff (e.g., change from grassland to irrigated or dryland monoc ulture) (Scanlon et al. 2007) as well as return flow from irrigation water and increased infiltration. Creating a predevelopment map from water level measurements is a paradoxical process. Where the water table does not intersect the land surface, water l evel measurements are taken from wells. Wells are typically drilled to extract water, thus development has prevailed across the approximate 450,000 km 2 aquifer. As a result, a predevelopment water levels map can best be created as a composite of predevelopment states across the aquifer, rather than using a particular point in time, when few measurements are available. To characterize the predevelopment state, a dataset was c reated representing the earliest measurements across the aquifer. Starting with the oldest measurement in each available well, a series of circular buffers were created ranging from 2 to 10 km in radius. The 8 km radius was selected as described below. A m easurement was retained only if it was the oldest water level measurement within its 8 km radius buffer (neighborhood), regardless of the number of wells present in its surrounding neighborhood. Declustering was not applied, despite the heterogeneous spaci ng of wells across the aquifer, in part because this neighborhood method generated a relatively regular spacing of measurements within the dataset. The resulting wells, each associated with the oldest water level measurement in an 8 - km radius, were collate d 11 into a dataset of 7958 measurements across the High Plains. Another dataset of 20852 (USGS High Plains Water Level Monitoring Study, 2006) was compared to the predevelopmen t dataset generated in this study ( Figure 1 . 3 ). Additional declustering did not improve the accuracy of the estimated surfaces, and was not used in the final analysis. The root mean square error of the kriged surface Figure 1 . 3 . asurement recorded for wells in the predevelopment dataset for (A) the USGS and (B) this p roject. 12 was approximately 10 meters. The kriged surface is not constrained to lie between the land surface elevation and the aquifer base; thus it was post - processed to enforce these limits. This post - processing involved: 1) a second interpo lation step in which additional points were Figure 1 . 4 . (A) The Northern High Plains, shown with well measurements (orange) and points from the National Hydrography Dataset (NHD) (blue). These points were used to create a predevelopment map. (B, C) A conceptual diagram showing the effect of NHD point incl usion. (B) Without NHD points, the interpolation technique ignores valleys where groundwater intersects the land surface, potentially overestimating groundwater levels. (C) With NHD points, the water level is constrained to the land surface. Note that NHD points were only used in areas where the initial water level interpolation was above the land surface; thus this technique does not cause the interpolated surface of the water level to scale and would represent considerable vertical exaggeration. 13 added from surface water feature elevations and 2) an explicit constraint to aquifer top and bottom elevations. A reduced sample of surface water feature elevations were added under the assumption that kriged are as above the land surface represented locations where groundwater and surface water intersect ( Figure 1 . 4 ) . To accomplish this, areas where the interpolated surface was above the land surface were isolated, and stream courses in these areas were selected from the National Hydrography Dataset (NHD). First, these stream segments were converted into several million point features created from the vertices of the NHD flow lines. As the inclusion of all of these points in building a revised kriging model would overwhelm the water level measurements from the wells, a 13 km x 13 km grid was created to overlay the aquife r and a maximum of one NHD point was chosen at random in each grid cell. The elevations at the NHD points, retrieved from the NED, were added to the point dataset for a second round of kriging. A total of 1984 stream points were merged with the 7958 predev elopment elevations to create a final predevelopment point dataset. After the second kriging, any remaining areas where the water level elevation lay above the land surface or aquifer base were truncated appropriately. A new kriged predevelopment water le vel surface was created with the 9942 points of the combined well and NHD points. The water level measurements that were used to generate this surface are highly variable in age ( Figure 1 . 3 ) and are not meant to average predevelopment state of the aquifer. The predevelopment water level dataset was developed as a synthesis of the 14 earliest water level measurements recorded in a given area. The size of that area, which will be referred to as the well neighborhood, is a defining factor in the density of points chosen as predevelopment measurements. The 8 km neighborhood radius provide d the greatest ratio of signal to noise (i.e., the smallest ratio between the root mean squared error of the kriging model relative to the measurements used, and the number of points included) compared to other neighborhood sizes tested (2, 5, and 10 km ra dii). Very small neighborhoods result in dense data. In other words, choosing the oldest measurement in any given 2 km radius will cause far more wells to be included than in a dataset comprising the oldest wells in a 10 km radius. However, this leads to the inclusion of wells that were drilled after the aquifer had already been impacted by nearby pumping. Inclusion of additional well measurements increases the local heterogeneity of measurement ages in the predevelopment dataset. A neighborhood larger tha n 10 km would exclude some early measurements and provide less data coverage. Calculations using the Theis equation also indicated that pumping would have little impact on water levels 8 km away after multiple years of pumping. These calculations were don e using transmissivity estimates typical for the Southern High Plains. Specific yield rather than specific storage was used, to simulate late - stage drawdown in an unconfined aquifer (Freeze and Cherry 1979). The first measured wells in a neighborhood tend to be similar in age, indicating a local onset of development (or monitoring). If more than one well from a neighborhood had data in the first year with recorded measurements, both measurements were kept in the database, though this 15 was uncommon. The 8 km radius was used for all measurement dates, and may be a conservative estimate for lower - capacity wells drilled prior to the 1960s when high - capacity wells (>100 gpm) became common (Luckey and Becker 1999). Annual Water Level Changes Maps of annual changes in water level were generated for water years 1935 - 2012. First, current year water level change was calculated by subtracting measured water two - year average was chos en to remove some of the noise within the dataset that can accumulate with the differencing method. For the first two years (1935 and 1936), the predevelopment surface was used instead of the prior two year average. The water level changes were then interpolated using ordinary kriging with an exponential semivariogram model (Isaaks and Srivastava 1989). The geostatistical structure of the water level change data shifts as development progresses, so parameters of the s emivariogram were fit to the data at decadal intervals, as finer intervals of updated variogram parameters did not appear to improve estimates. Annual interpolations were automated in Python to increase speed, prevent human error, improve uniformity, and m ake the process reproducible. The result of the kriging is an interpolated map of the change in water level relative to the average of the past two years. The kriged surface was then added to the average water level elevation from the previous two years to produce an estimate of the current year water level. The partial sill decreased nearly linearly from approximately 21.5 to 10.5 m 2 from 16 1935 - 2012, because the density of measurements increased, and the difference between water level changes in a given ra dius tended to decrease while variability in head increased. The sill (which includes the partial sill plus the nugget) is the amount of variation that exists between points at the distance of the range, where the semivariogram flattens out (Isaaks and Sri vastava 1989). Other semivariogram parameters did not change significantly in the optimization through time, thus the nugget was fixed at 5 m 2 and the major range was held constant at 42000 meters, as it did not change significantly based on the variogram fits through time. Once water level maps were generated from the annual kriging, a second interpolation step, using surface water feature elevations and truncation of areas with estimated water table elevations above the land surface and below the aquifer base, was applied as described above. Water levels were also constrained to change no more than 7 m per year (representing three standard deviations from the average change) within the resultant 250 m x 250 m grid. This occurred over very small portions o f the aquifer in any given year, typically less than 1% of aquifer area. Areas with zero saturated thickness at predevelopment were also assigned the aquifer base elevation for subsequent water years. Volume estimates were obtained by multiplying estimated saturated thickness by the specific yield of the aquifer (McGuire et al. 2012; Figure 1 . 5 ). 17 Water Level Decline Forecast After the final set of water level maps were generated, the trend in water level elevation was calculated for each grid cell from 1993 - 2012. A twenty - year window was long enough to sho w meaningful trends in water level over a period of relatively monotonic water level changes; a linear trend proved to be a good fit to the data over Figure 1 . 5 . Specific yield estimates from McGuire et al. 2012. 18 that timespan for areas of rapid decline. This analysis resulted in a water level trend map which was used to project water levels into the future, should recent rates of groundwater extraction continue. A date o f depletion was then calculated as the year in which the saturated thickness remaining would fall below a threshold of 9 m, likely indicating insufficient aquifer transmissivity to maintain pumping rates for high - Figure 1 . 6 . The water level elevation (piezometric surface) during predevelopment, in meters. 19 capacity irrigation wells in that cell. Re sults and Discussion Predevelopment Map and Aquifer Storage The neighborhood method used in this analysis results in a dataset that is smaller and represents earlier well measurements on average than those included in the USGS predevelopment dataset ( Figure 1 . 3 ). Earlier measurements are less likely to have large well measurements, whi le still maintaining adequate spatial coverage, also reduces the likelihood that the average well in the dataset has been greatly impacted by pumping as of its first available measurement. This new predevelopment water level map ( Figure 1 . 6 ) of the High Plains Aquifer, consistent with previous work, shows a clear gradient from a high water level elevations in the west to a lower water level elevations in the east, simi lar to maps from Cederstrand and Becker (1999) and Luckey et al. (1981). This gradient drives the regional flow of groundwater toward the east, where the eastern margin of the aquifer was a historic area of natural discharge (White et al. 1946). The total volume water in storage in the High Plains Aquifer, under predevelopment conditions, is estimated at 4230 km 3 . Annual Water Level and Storage Changes Water level change has been most profound in the Southern and Central High Plains, consistent with other studies and with individual well hydrographs, as well as 20 remotely sensed estimates ( Figure 1 . 7 ; Breña - Naranjo et al. 2014). The bulk of the NHP has experienced little net change . M uch of this area has sparse data availability, and consequently high mean standard deviations relative to the rest of the aquifer ( Figure 1 . 8 ). The CHP has seen the greatest declines in the western regions of the aquifer, as noted in Wilson (2007). In the SHP, much of th e Figure 1 . 7 . (A) Water level change, from the predevelopment surface to 2012. (B) Hydrographs of individual wells from eachof the regions of the High Plains (USGS site IDs: SHP, 341420101441602; CHP, 362746102364102; NHP, 404345098560001). (C) Average water level change across the three regions of the High Plains. 21 predevelopment saturated thickness has been depleted ( Figure 1 . 9 ), particularly in the n orthern part of the region. When the interpolated water level change was add ed back to the mean water level for the previous two years, areas with no measurements Figure 1 . 8 . Mean standard deviation of water table elevation kriging estimates. 22 Figure 1 . 9 . The percent remaining saturated thickness in the SHP in 2012 compared to 1943 (the year of maximum saturated thickness for the SHP). Note that areas shown in dark gray are locations where the saturated thickness was interpo lated to be zero at predevelopment. 23 showed little change. This meant that areas that were not developed in 1935 stayed relatively consistent , close to the predevelopment levels, until measurements began to be available. While the predevelopment water lev el map in this study was developed using a carefully selected subset of the data made available by the USGS and Nebraska DNR, the annual interpolation method employed here takes advantage of more of the available data than previously published interpolatio ns. To increase temporal consistency, prior reports (e.g. McGuire 2004, 2009, 2012) included only wells present in both the predevelopment and later comparison periods, thus excluding the majority of well data, as a large number of wells are measured only once or a few times. Prior studies were intended to give semi - regular updates of aquifer status, not necessarily for the types of water resource investigations that annual maps can provide. The effort described in this paper also presents methodological e nhancements, including the application of kriging in a two - step interpolation process which enables the addition of surface water features elevations in areas with gaining streams. The inclusion of surface water feature points helped define stream valleys and low water levels generated by loss of groundwater to gaining streams, lowering the interpolated water levels ( Figure 1 . 4 ) and minimizing the need for subsequent t runcation at the land surface. the High Plains Aquifer contained about 4160 km 3 of water. Note that this is 70 km 3 lower than the predevelopment estimate, due likely to the drier overall climate and drought conditions prevalent during the 1930 s, in combination with the selection of the earliest data at each location in the predevelopment map. Of the 4160 km 3 of water in 24 storage, approximately 2948 km 3 was in the North ern High Plains, 887 km 3 was in the Central High Plains, and 326 km 3 was in the Southern High Plains. The annual maps were used to calculate changes in water level and volume (i.e., change from 1935 to 2012 as opposed to change from predevelopment to 2012 ). The predevelopment map represents a composite of water levels measured prior to significant local development, rather than the piezometric surface in a particular year. By 2012, aquifer storage volume had decreased to about 3750 km 3 of water (NHP: 2940 km 3 , CHP: 636 km 3 , SHP: 171 km 3 ), a total decline of about 410 km 3 ( Figure 1 . 10 ). Figure 1 . 10 . (A) The total volume of water in the High Plains Aquifer from 1935 to 2012. (B) Water remaining in storage from 1935 to 2012, as a percentage of the maximum water in storage. 25 Table 1 . 1 . A comparison of storage changes from recent publications. The USGS predevelopment measurement dataset yields an aquifer volume of 4086 km 3 when interpolated using the same method as the predevelopment dataset generated in this study. This volume is 65 km 3 less than our estimate of storage in 1950. McGuire (2012) estimated that net decrease in storage was 304 cubic kilometers from predevelopment to 2011, somewhat less than the estimate of 328 km 3 from predevelopment to 2007 from Stanton et al. (2010) and McGuire (2009 ; Table 1 . 1 ). This is about 15% less than the change calculated in this study for 1950 - 2011 (385 km 3 ), and 30% less than change from predevelopment to 2011 as calculated in this study ( about 460 km 3 ). The disparities among these estimates indicate the degree to which changes in interpolation methods can alter the resulting estimates. While this study employs more complex methods than previous reports, uncertainties remain difficult to quantify when t aking into account all datasets upon which these estimates rely. It is unclear whether differences in total storage volume are within the bounds of uncertainty for the various reports. However, the annual changes in water level agree substantially with pre vious reports, and offer greater potential for certainty in areas with high data availability. The overall pattern of depletion shows a general trend of development followed by periods of rapid depletion ( Figure 1 . 10 ), most of which is monotonic with little Est. Water in Storage, km 3 Est. Decline in Storage, km 3 1950 2007 2011 Estimated in this report 4151 3776 3766 385 (1950 - 2011) Stanton 2011/McGuire 2009 3914 3586 330 (1950 - 2007) McGuire 2012 3950 3650 300 (1950 - 2011) 26 recovery following the onset of development. Declines in water level and volume became much more rapid in the 1950 s, which coincides with the installation of many high - capacity wells (Luckey and Becker 1999). Each subregion (Northern, Central, and Southern) has a distinct history of exploitation and a unique pattern of response. The NHP, which includes the thickest part of the aquifer (the Nebraska Sand Hills), ex perienced the greatest delay in irrigation development. This region is less susceptible to groundwater depletion than other areas of the aquifer due to both higher precipitation and recharge. In addition, withdrawals are lower over much of the NHP, and the onset of groundwater irrigated agriculture was significantly delayed in the region relative to the SHP and CHP ( Figure 1 . 3 ). There is no clear signature of overall aq uifer development based on depletion for this region, although a few counties have experienced water - level declines. This is in broad agreement with previously published maps of water level change in Nebraska (e.g. Young et al. 2012). The vast majority of the decline in storage has occurred in the SHP and CHP, where total storage has always been much smaller than the volume of water in the NHP. The CHP region has a heterogeneous response to development, with some areas experiencing water level reductions o n the order of 50 meters ( Figure 1 . 7 a ). This is true in some areas where the total saturated thickness is 80 meters or less. Other areas experienced very little chang e, including some areas with very little saturated thickness despite being within the bounds of the geologic formations of the aquifer, which affected the distribution of irrigated agriculture. Overall declines in the CHP have 27 followed a fairly linear pattern, with a decrease in saturated thickness on the order of two meters per decade since the mid - 1950 s, with a recent increase in the rate of decline beginning in 2010 due to a widespread drought. The SHP has seen the greatest percentage declines in saturated thickness. In the northern third of the SHP, average saturated thickness has declined by approximately 75% ( Figure 1 . 9 ), and much of the remaining water is inaccessible for high capacity Figure 1 . 11 . Rate of decline in percent per decade, using a from 1992 to 2012. 28 pumping due to the thinning of the saturated zone. Water levels in the SHP declined most rapidly prior to the 1970's ( Figure 1 . 7 c ), at which point depletion slowed, the cause of which is under further investigation using complemen tary approaches including both statistical and process - based modeling. Recent Rates of Decline and Depletion Forecast Figure 1 . 11 shows the current rate (calculated from 1993 - 2012) of water level decline as a percentage per decade of initial saturated thickness. Some of the most agriculturally productive areas of the High Plains are experiencing rapid depletion, including much of the SHP, the western CHP, and the southwestern NHP. Sustained rapid drawdown due to irrigation wells will lead to a continued increase in the percent of the High Plains area where the aquifer cannot supply adequate water for irrigation (defined in this analy saturated thickness). If drawdown continues at the present rate ( Figure 1 . 12 ), the Figure 1 . 12 . Percent of aquifer area with less than 9 m saturated thickness in future decades, assuming a linear pattern of decline. 29 portion of the aquifer unsuitable for irrigation will proceed from less than 25% (including areas that have never supported irrigated agriculture) as of 2012 to approximately 40% by 2100. Estimates from the US Census of Agriculture and satellite imagery ind icate that roughly 13% of the HPA was irrigated as of 2002 (Stanton et al. 2011). Thus, the aquifer underlying a substantial fraction of currently irrigated land, Figure 1 . 13 . Projected depletion timeframe for the entire High Plains Aquifer, using a 9 m threshold saturated thickness for aquifer utility and assuming a linear pattern of de cline. Gray areas are in little danger of depletion in the foreseeable future. 30 including most of the SHP and CHP, would become depleted by 2100. Figure 1 . 13 illustrates the dramatic increase in depleted area across the SHP and CHP, defined as a location with less than 9 m of saturated thickness. The western parts of the SHP and CHP are at particular risk of depletion in the next few decades. Parts of the aquifer that have been heavily exploited in the past will continue to be drawn down by 10 - 20% per decade for the remainder of their lifespan unle ss drastic changes in management are implemented. In contrast, many areas of the aquifer are in little danger of depletion in the immediate future, including parts of the High Plains that have never been suitable for irrigated agriculture along with those in more humid portions of the aquifer such as the eastern CHP and NHP. Much of the Southern High Plains has already been depleted, and more is likely to be depleted in the next few years. Depletion slows as more water is lost from storage, causing farmer s to switch from irrigated to dryland farming. Nevertheless, some of the most heavily exploited parts of the aquifer are likely to become unusable for irrigation in the next few decades. Conclusions Available water storage in the aquifer has declined by approximately 410 km 3 since 1935, with approximately 15% more storage decline up to 2011 than estimated by prior studies (McGuire, 2012; Stanton et al. 2011). Although the total depletion is a little less than 10% of initial aquifer storage, water level de clines have not been uniformly distributed throughout the aquifer. Most of the water storage in the aquifer is in the 31 NHP, where little net decline has occurred. In contrast, large portions of the SHP and CHP are experiencing more than 10% decline per deca de. This analysis indicates that based on current trends, the lifespan of some of the most productive portions of the aquifer can be measured on the order of decades, with substantial depletion of currently irrigated areas of the SHP and CHP by 2100. The current approach to resource use is managed mining, with the understanding that, given the low rates of recharge prevalent in the SHP and CHP, the aquifer is effectively a nonrenewable resource in many locations. Recharge rates have not been sufficient to maintain aquifer levels in any large area of the High Plains without significant reduction in discharge to surface water. Much of the SHP and CHP subregions constitute a finite reservoir that can be exploited more efficiently with improved stewardship, wh ich will involve decreased pumping, perhaps coupled with enhanced recharge with water from streams and rivers when sufficient water is available, or from water diversions. Improvements in technology and data availability are making it possible to correlate water level changes in wells with gravimetric changes measured over broad scales, using the GRACE satellite (e.g., Breña - Naranjo et al. 2014). The predevelopment water level is not a goal for future restoration of water levels, nor does it alone provide a great deal of insight into the expected lifespan of the aquifer. However, it does give an idea of the starting point of the aquifer, along with a realistic baseline for modeling studies. A steady - state predevelopment groundwater model for the entire aqui fer would provide reciprocal verification. Any study of the past encourages projection of future trends. In the case of the 32 High Plains Aquifer, the continuation of past management practices is likely to lead to wholesale depletion across much of the cent ral and southern aquifer regions, forcing major changes in land use and production. However, alteration of management practices is already a primary focus for many farmers and supporting organizations. This aquifer depletion, while unlikely to cause a sud den crisis (as its trajectory differs across space), may ultimately reduce agricultural production to the levels seen on the edges of the SHP and western margin of the CHP, where only a minimal aquifer thickness has been present, and which by comparison ha ve never been major areas of production. food security, economy, and environment. Continued water level decline is inevitable while current management policies continue, but th e sustainability of this resource the largest freshwater aquifer in the country is a vital concern. More concerted stewardship efforts must be put into place if we are to manage the depletion of the High Plains Aquifer in an economically and socially r esponsible way. 33 Chapter 2 Effects of Management Areas, Drought, and Commodity Prices on Groundwater Decline Patterns across the High Plains Aquifer Abstract I use an 82 - year record of water table data from the High Plains Aquifer to introduce a new application of segmented regression to hydrogeology, and evaluate the effects of droughts, crop prices, and local groundwater management on groundwater level trajectories. Across the High Plains, there are discernable regional cycles of faster and slower water table declines. About 25 - 40% of the variation in decline rates is attributable to the effects of climate and crop prices on the demand for irrigation water. Both droughts and high crop prices contribute to cycles of faster decline, with variable lags for the on set of recovery. The Southern High Plains takes longer to begin recovery than the Central High Plains; the Northern High Plains responds much more quickly, especially to changes in climatic conditions. T he establishment of groundwater management areas is o ften accompanied by faster water table declines for 10 - 15 years at a higher rate than predicted based on climate and commodity prices alone. Segmented regression can provide a tool for groundwater managers to evaluate the effectiveness of management strate gies on groundwater storage and decline, using readily available water table data. 34 Introduction and Background Groundwater Depletion and Management Groundwater in major areas of the High Plains Aquifer is being rapidly depleted (McGuire 2012). For the past few decades, attention has focused on the likely pathways and timescales to aquifer depletion (e.g., Scanlon et al. 2012; Steward and Allen 2016; Haacker et al. 2016). Irrigation is becoming more and more dependent on groundwater. Over - allocation of s urface water resources, shifting demands on water supplies, changing climate and variable precipitation, and growing awareness of the ecological steady source of water. Me anwhile, agricultural intensification and increased demand for animal feed are placing greater stresses on existing water supplies. Despite the critical role that groundwater plays in providing water for irrigation, groundwater management in the United Sta tes and elsewhere, is in its early stages relative to surface water management (Megdal et al. 2014). For example, in California, formal statewide management of groundwater began in 2014 with the Sustainable Gro undwater Management Act (Babbit et al. 2018), users of surface water (Reisner 1986). Groundwater laws in the United States have primarily been enacted at the state level (Kirchhoff and Dilling 2016). Despite most states having policies in place, the elevation of the water table continues to decline in most semi - arid areas with groundwater irrigation (e.g., Smidt et al. 2016). Many states have established smaller, localized groundwater management areas (GMAs) to distribute control of resources 35 based on unique local conditions, including water level declines. These GMA s were frequently created to serve the needs of local irrigators, rather than the needs of local ecosystems. However, some groundwater management areas have evolved in function over time in response to changing local priorities and governance. Megdal et al. (2014) found that the main concerns for groundwater governance at the state level include groundwater contamination, user conflicts, and water level declines; despite the physical nature of these concerns, GMA s have generally been evaluated based on social metrics rather than environmental change (Kirchhoff and Dilling 2016). About two decades after local GMA s were established in the High Plains, a survey in Colorado and Kansas found that two thirds of irrigators felt that their local groundwater management areas should work to reduce groundwater depletion (White and Kromm 1995). Only half felt that their GMA s should attempt to maintain the resource over a specific time horizon, although th is may be changing as declines have become more severe. The physical effects of GMA s have generally not been included in studies that estimate the lifespan of aquifers such as the High Plains Aquifer, in the central United States. Previous statistical mod els of aquifer lifespan have assumed that decline is linear (Haacker et al. 201 6 , Scanlon et al. 2012) or logistic (Steward and Allen 2016), or otherwise excluded GMA s from consideration as predictive variables for groundwater levels (Sahoo et al. 2017). T hese assumptions of the trajectory of water table change preclude any consideration of the changing role of management districts in statistical or machine learning models. Moreover, well records do not support the assumption of 36 monotonic decline that is bu ilt into linear or logarithmic fitting. It is vitally important to recognize that the lifespan of the aquifer can be extended by slowing or reversing declines; on the other hand, an acceleration of decline suggests a failure of management to work toward gr oundwater sustainability. Here, I demonstrate that these cycles of faster and slower decline have occurred throughout the High Plains Aquifer, as a result of variable climate, commodity prices, and localized groundwater governance. I quantify changes in s lope (breakpoints) of groundwater levels through time and correlate those breaks in slope with drought and fluctuations in crop prices. I compare the net number of annual breaks with those expected before and after the establishment of GMA s. To better unde rstand the role of GMA s in altering the trajectory of water level changes across the High Plains Aquifer, I employ a method rarely used in water resources studies, segmented regression (Muggeo 2008; Shao and Campbell 2002); I demonstrate its utility for assessing the influence of groundwater management areas on water levels. Study Site The High Plains Aquifer is the source for nearly one third of all groundwater used for irrigation in the United States. 98% of water extracted from the High Plains Aquifer is used for irrigation (Dennehy et al. 2002). Two features make this aquifer a particularly important resource: it is very large, and much the overlying arable land receives insufficient precipitation for reliable rainfed agric ulture. With 450,000 km 2 of land area over parts of eight states, the High Plains Aquifer contained approximately 37 3,500 km 3 of water in 2016, although it is very unevenly distributed. Much of the High Pla ins region has a saturated thickness that is barely sufficient for irrigation (Haacker et al. 2016), depending on the hydraulic conductivity and other factors (Hecox et al. 2002). However, the saturated thickness varies widely, from effectively zero in muc h of New Mexico, to hundreds of meters in the Nebraska Sand Hills. Much more water is withdrawn than returns to the aquifer each year (Scanlon et al. 2012; Haacker et al. 2016). Water from the High Plains Aquifer contributes significantly to state economie s. Figure 2 . 1 . State - level groundwater doctrines for the High Plains (Joshi 2005). The three r egions of the High Plains Aquifer Northern, Central, and Southern are overlayed in gray. Some states use a modified system as described in Smidt et al. (2016). 38 From 1996 - 2005, Kansas farmers annually withdrew about $110 million USD worth of groundwater (Fenichel et al. 2016). The combination of hydrogeologic and economic factors make groundwater management from the High Plains region both difficult and necessa ry. Based on c limate (potential evapotranspiration, ET) and hydrogeologic gradients across the High Plains, the aquifer is commonly divided into three subregions: Northern, Central, and Southern (Stanton et al. 2011; Figure 2 . 1 ). Mean annual temperatures are higher to the south, and precipitation roughly doubles from west to east (Mahmood and Hubbard 2002; Figure 2 . 2 a,b). The Northern High Plains, which has the lowest average evaporation demand, has greater mean saturated aquifer thickness Figure 2 . 2 . (A) Mean annual precipitation in millimeters, from 1930 - 2015 (PRISM); (B) mean annual temperature in degrees Celsius, from 1930 - 2015 (PRISM); (C) g roundwater level decline from predevelop ment through 2016, estimated using methods from Haacker et al. (201 6 ); (D) irrigated/non - irrigated area from MODIS - MIrAD satellite data, 2002, 2007, and 2012. Groundwater management areas with governance below state level are outlined in black. 39 and recharge than the Central or Southern High Plains, where the water - bearing geologic units are thinner and there is little natural recharge. Across most of the Central and Southern High Plains, pumping for irrigation has caused large declines in groundwater levels ( Figure 2 . 2 c). Although many early irrigators in Texas believed that their groundwater was an inexhaustible resource Figure 2 . 3 . Es tablishment year for groundwater management areas that are predominantly on the High Plains. Select establishment years are shown on the map for reference. States are labeled with postal abbreviations. 40 (Thei s, 1937; White et al., 1946). However, increased understanding and awareness has not forestalled the intensive growth of regional farm economies over the past eight decades economies that rely on unsustainable groundwater resources ( Figure 2 . 2 d). Groundwater Management Areas Nebraska, Colorado, Kansas, and Texas together cover about 80% of the area overlying the High Plains Aquifer. The se four states have comparable G MA s. Key features of these GMA s are that they retain some jurisdiction over local natural resources, including groundwater, and can raise funds and create rules and obligations for irrigators within their boundaries. The 41 GMA s in this study are predomina ntly within the boundaries of the High Plains Aquifer ( Figure 2 . 3 ) and were typically instituted a few decades after local aquifer development. Beyond these parameters, GMA s are diverse in their operation and functions. s Districts, which are primarily funded through property resources at a semiautonomous, local level. G MA s in Colorado and Texas, on the other hand, were organized largely to support the needs of local irrigators and maintain profits (White and Kromm 1995). In Texas the comparable GMA structures are Underground Water Conservation Districts, not to be confused with state - designated Groundwater Management Areas, which play a very different, top - down governance part of their name) were created specifically to thwart the adoption of state - level 41 groundwater use restrictions (Green 1981). Kansas has taken an intermediate approach between Nebraska and Texas to create management districts, ensuring a conservation agenda for groundwater districts that are established, and yet mandating that districts be created only at the request of irrigators (Bossert 1993). Despite their roles in supporting regional irrigation, some GMA s in Colorado, Texas, and Kansas have increasingly adopted the role of conservation entities; a prime example is the creation of the Sheridan 6 Local Enhanced Management Area in Kansas Groundwater Management District 4. The Local Enhanced Management Areas are not specifically examined in this paper because the original Sheridan 6 district covers a small area (six sections of land, 256 km 2 ), about one third the size of the smallest GMA co nsidered in this study (750 km 2 ) and entirely within an existing GMA . Regulations within GMA s generally include guidelines for well spacing, monitoring, and capping of abandoned wells (White and Kromm 1995). These may be binding or voluntary, and may come with subsidies or penalties. The institution of pumping restrictions is much more recent, and highly variable across the aquifer. This paper focuses on the initiation of management, since this event is comparable across all management areas and thus yield s the largest dataset for statistical analysis. The motivation for this study was to identify any changes in groundwater level trajectory that occur when a localized GMA is established on the High Plains. Following the establishment of GMA s, positive infle ctions in water levels as declines slow or even reverse might indicate conservation of water resources. The mechanisms of 42 conservation within GMA s can take numerous forms including: taxation, education, regulation, establishing social norms, monitoring, an d incentivizing more efficient technologies. Factors Affecting Irrigation Demand Aside from factors under the control of GMA s, a host of conditions can affect the demand for irrigation water in the High Plains. Changes in agricultural technology and mana gement may be driven by broad market and non - policy factors, including shifts in crop prices, changes in crop genetics, or interstate river basin compacts. There are regulatory, legal, and policy actions that may affect water use at scales greater than GMA s (i.e. state, interstate, or federal). The region is very susceptible to drought and groundwater depletion. Isolating every factor that may have influenced past water level trajectories is beyond the scope of this paper. Rather, I seek to control for the most important market - based and physical forces that might have affected water use, and for which data are available across the entire area for the decades since aquifer development. Climate and crop prices are two such factors. Groundwater serves as a r eservoir in times of drought, and as an additional supplement for high yields when commodity prices are high. Agricultural drought occurs at a regional scale, whereas crop prices are national. Both of these patterns stand in contrast to other important det erminants of water demand, such as extreme climatic events and local economic conditions. Therefore, these factors reflect overarching conditions for irrigators, which are 43 applicable throughout and beyond the scale of individual management areas. Methods Analysis Objectives and Procedure Our analysis takes four parts, as follows. 1) Algorithmically identify breaks in breakpoints and aquifer depletion, climate, and commodity pric es, within GMA s and aquifer regions. 3) Develop a model of the influence of the most significant crop price and climate parameters on changes in water level trajectory. Finally, 4) using these models to control for the influence of climate and prices, eval uate patterns in net annual breaks relative to the establishment of GMA s. Scripts are hosted at github.com/ehaacker. Spatial and Temporal Datasets A spatial database of GMA boundaries was created for Texas, Kansas, Nebraska, and Colorado, which included t he year of formation of each GMA ( Figure 2 . 3 ). Texas has both the earliest - (1951) and latest - formed (1998) groundwater management areas. Most GMA s were formed in the late 1960 s to mid - 1970 s. Cross - correlation and multiple linear regression analyses required regional data on Palmer Drought Severity Index (PDSI; Palmer, 1965; NOAA 2017); precipitation from the PRISM reanalysis product (PRISM 2004, updated 2017); and c ommodity prices (NASS 2017 a ). PDSI, precipitation, and commodity prices are available throughout the region and time period considered in this study. 44 Annual PDSI data were available (NOAA 2017) at low resolution (~ 2 degrees square of latitude and longitu de, which is ~250,000 square kilometers and roughly the scale of the Southern High Plains). PDSI incorporates precipitation and temperature to derive a more holistic measure of agricultural drought than precipitation alone. Negative PDSI numbers indicate d rought, while positive numbers signify unusually moist conditions. The PDSI magnitude indicates the severity of the deviation from the regional mean climate. The PDSI incorporates both precipitation and temperature, and can thus provide a more holistic rep resentation of plant stress than precipitation alone. Annual PRISM precipitation data at a 4 km resolution (PRISM Climate Group 2004) were used directly in addition to PDSI. PRISM is a reanalysis product that uses precipitation gauge measurements and othe r available data to create an estimate for average annual regional precipitation. The correlation coefficient between PDSI and precipitation is about 0.8. National prices for corn, sorghum, and cotton were obtained from the National Agricultural Statistics Service (USDA 2017). Records used for this analysis were for the commodity total price received for the crop year, adjusted to 2016 dollars. Prices for corn, sorghum, and cotton are strongly correlated with one another (>0.95). An estimate of aquifer bas e elevation was needed to calculate saturated thickness. The aquifer base elevation data was amalgamated from four sources: the Texas Groundwater Availability Model for Texas, New Mexico, and Oklahoma (Deeds and Jigmond 2015); a contour map provided by the Kansas Geological Survey (kgs.ku.edu) for the Kansas High Plains; borehole estimates from USGS Data Series 777 45 (Houston et al. 2011) in the Northern High Plains; and an earlier contour map from the USGS (Cederstrand and Becker 1998 a ) for parts of Colorado that were not covered by other datasets. Point - scale datasets and contour maps were converted from feet to meters and interpolated to create continuous rasters of aquifer base elevation, and these rasters were combined using a raster catalog function in A rcGIS. Well Hydrographs This analysis used USGS depth - to - water data for the High Plains Aquifer (USGS 2017; Qi 2010). Prior to analysis of the depth - to - water data, I performed a series of quality control procedures. Wells with a depth to water measurement of zero were omitted, along with well measurements that were flagged in the USGS database to indicate that levels were not static. Measurements taken during the growing season (May to September inclusive) were excluded to minimize seasonal effects from irr igation pumping. Water level measurements for each well were then averaged over the October 1 to September 30 water year. For example, measurements used for the 1980 water year were taken between October 1, 1979, and April 30, 1980, thereby excluding the i rrigation season. Water table elevations were calculated by subtracting the depth to water from the ground surface elevation, determined at the well location from the 1/3 arcsecond elevation grid from the USGS 3D Elevation Program (USGS 2016). To analyze changes in water table trajectories, I restricted the dataset to wells that had measurements for ten or more (not necessarily consecutive) water years. Most wells had gaps in their record with no measurements for one or more water years; 46 specifically, many wells are measured every other year. Thus, the ten measurements for a given well may be spread out over fifteen to twenty years. The breakpoint analysis required interpolating data for missing years, thus any wells with a measurement gap greater than two consecutive years were excluded. Following these quality control steps, an initial dataset of 106,716 wells with at least one valid measurement was reduced to a database with 12,441 individual well hydrographs. Valid wells for breakpoint analysis have reco rds as early as 1929, with the maximum number of records occurring in 1970, when there were 2132 valid wells. Additionally, saturated thickness of the aquifer was calculated for each well, by subtracting the aquifer base elevation from the water table elev ation. Since the aquifer base elevation is constant, the slope of saturated thickness over time is the same as the slope of the water table over time. For each year, an average saturated thickness value was calculated for all wells with and without breakpo ints in the Northern, Central, and Southern High Plains. Breakpoint Analysis The first part of our analysis used segmented regression (Shao and Campbell 2002; Muggeo 2008), a broken - line fitting technique, to identify statistically significant changes in slope, or breakpoints, in water table elevation over time (well hydrographs; Figure 2 . 4 (Muggeo 200 8), which uses a modified linear regression to algorithmically identify a slope break (or breaks) in a data series (illustrated in Figure 2 . 4 47 pa ckage I fit a linear model to each well water level time series, allowing for one possible breakpoint defined as the decimal year at which a break in slope occurs. Allowing for only a single breakpoint does simplify the often complex drawdown history of ea ch well, but importantly, this restriction avoids overfitting, which increases confidence in the physical and statistical significance of each breakpoint. After breakpoints were determined, I calculated the significance of each breakpoint using the Score test (Muggeo 2016). If the Score test calculated p - value was less than 0.05, the breakpoint for that well was considered significant, and was included in further analysis. Since the initial breakpoint estimate is initialized using a random seed rather than being a user - input variable, the segmented regression was run 100 times for each well hydrograph, and the mode of the resulting breakpoint estimates was used for further analysis. The mode breakpoint year differed from the mean in less than 2% of well hyd rographs, suggesting low sensitivity to the initial breakpoint estimate. Of the 12,022 wells (97% of total wells) with identified significant breaks in slope, more than 96% had a significant breakpoint located in all 100 runs. For wells with a significant break in slope in all runs, I then quantified the slope change across the break as either positive or negative deflections. A positive deflection across the breakpoint indicates slowing decline, recovery of the water table (e.g., Figure 2 . 4 a), or accelerating increases in the water table (in < 4% of positive breakpoints). Likewise, negative deflections indicate the start of water table decline, acceleration of decline (e.g., Figure 2 . 4 c), or slowing of water table increase (in < 3% of negative breakpoints). 48 I then defined an aggregate measure of breakpoint behavior within desired s ummar where is the net annual break fraction calculated as the difference of the positive and negative break fractions; is the year of each time series, and represent the total count of positive and negative breakpoints in a given year and area, and; is the total count of well records for the year in the summary area. The total positive ( ) and Figure 2 . 4 . Example well hydrographs with segmented regressions for well locations displayed in the map (right). GMA establishment dates are noted for each example well. R 2 values for segmented regressions are reported for wells A and C. Wells B and D did not have statistically s ignificant breaks in slope; R 2 for those wells are shown for linear regression. GMA = Groundwater Management Area. Water table elevations are shown in meters above sea level. 49 negative ( ) breakpoints w ere calculated individually as the annual sum of positive and negative breakpoints within one of two spatial frameworks: the 41 GMA s across the four states (TX, KS, CO, and NE), and the three High Plains Aquifer regions (Northern, Central, and Southern). T he number of breaks was divided by the number of measurements to normalize the effects of years with only a small number of wells, particularly toward the beginning of the time series. Aside from net break fraction, two additional single - valued spatially - s ummarized breakpoint time series were investigated: the mean annual magnitude of slope change and net breakpoint count alone. I found the net breakpoint count and mean annual magnitude of slope change were susceptible to changing numbers of well observatio ns through time, and were therefore less meaningful overall. Furthermore, the net break fraction was more stationary than these other assessment variables according to an augmented Dickey Fuller test used to ensure that the mean of a data series did not ch ange over time (a requirement for further statistical analysis). Cross - Correlation and Multiple Linear Regression Analysis To examine the influence of market - based and physical factors on well hydrograph breakpoints, I compared the annual time series of b reak fractions with annual time series of the regional precipitation, PDSI, and national prices for corn, sorghum, and cotton, using cross - correlation functions. Additionally, I compared the break fractions with time series of annual differences ( ) in climate and commodity prices. The cross - correlation functions test the magnitude and direction of 50 simple correlation between pairs of variables, both for the present - year and successively increasing annual lags, to determine whether correlation ma gnitudes are higher when the two time series are offset by one or more years. This occurs, for example, when changes in crop prices or drought indices are associated with positive or negative break fractions a few years later. I next built a multiple line ar regression model describing net break fractions for each of the three major High Plains Aquifer regions, starting with independent variables at all lags with statistically significant ( p <0.05) correlation according to the CCF analysis. Regression models were built using data from 1950 - 2016, to avoid the models were used to quantify how much variation in net break fraction could be attributed to climate and economic deman d for crops. Independent variables (climate and crop price at each significant lag) were eliminated from the models if the regression coefficients were not significant to p <0.10 or better. Most of the lags for these parameters were not significant in the multiple linear regression model, even though they were significant in the CCF analysis, because of high autocorrelation within each independent variable. For example, the price of agricultural products is much less volatile than precipitation; thus, this a decent predictor for prices next year. Because of this autocorrelation, the inclusion of a moving average or sum of years (for example, to represent prolonged high prices or drought) did not add information to the model, because it led to increased 51 autoco rrelation. If a parameter was significant to the linear regression model at multiple lag years ( p <0.10), then all significant lag years were retained ( Table 2 . 1 ). After creating a multiple linear regression model for each region of the High Plains Aquifer, the models were used to estimate the expected net number of wells with Variable Lag Estimate St. Dev. Error Significance NHP Corn Price - 4 0.0017 0.00084 2.0 0.051 . NHP Palmer Drought Index - 1 - 0.00067 0.00021 - 3.2 0.0022 ** NHP Palmer Drought Index - 2 - 0.00043 0.00023 - 1.9 0.067 . NHP Palmer Drought Index - 3 - 0.00049 0.00022 - 2.2 0.034 * CHP Precipitation - 2 - 0.000084 0.000032 - 2.6 0.011 * CHP Precipitation - 3 - 0.000071 0.000033 - 2.2 0.032 * CHP Change in Sorghum Price 0 - 0.0033 0.0013 - 2.6 0.011 * SHP Palmer Drought Index - 3 - 0.00039 0.00020 - 1.9 0.057 . SHP Palmer Drought Index - 4 - 0.00042 0.00022 - 1.9 0.059 . SHP Change in Drought Index 0 0.00052 0.00018 2.9 0.0047 ** Table 2 . 1 . Factors in the linear regression models shown in Figure 11. Not all cross - correlates shown in Figure 10 were significant in the multiple linear regression number of years between the variable and the net annual break fraction; negative numbers indicate that the variable is measured prior to the annual break fraction, e.g. a lag of - 4 is four years prior to the net annual break fraction. fficient of the variable with regards to the net annual break fraction; positive estimates mean that high values of the variable are associated with deviation and error of the es estimate of likelihood that the relationship between the variable and the annual break fraction is due to chance. NHP: Northern High Plains; CHP: Central High Plains; SHP: Southern High Plains. 52 breakpoints in a given year, based on the expected response of water table elevation change s to the time series of commodity prices and climate conditions. The models yielded an annual dataset of expected net break fraction, which was multiplied by the total number of well measurements in that region in that year, to estimate the expected net nu mber of breakpoints. For each GMA , the expected net number of breakpoints was determined for the years before and after the GMA was established. This expected net number of well breakpoints was then compared with the net observed breakpoints for the years before and after GMA s were established. The difference between net observed breakpoints and net expected breakpoints was calculated to generate an annual dataset of net residual breakpoints. This residual value represents the net number of breakpoints that are not explained by commodity prices and climate conditions, and can thus be attributed to some other factor. Results and Discussion Regional Patterns The Northern, Central, and Southern High Plains each have distinct cycles of positive and negative defl ections in water table trajectory over time ( Figure 2 . 5 ). Within a region, years with many positive breakpoints tend to have very few wells with negative breakpoints, and vice versa. This was unexpected, considering the variation in well design and management strategies that might otherwise lead to random patterns of annual break fractions. Despite this consistency, strong positive or negative signals seen in one region are frequently muted or absent in other regions. For example, in the 53 ins experienced a significant acceleration in aquifer depletion, evidenced by dominantly negative breakpoints, while the Central High Plains witnessed some slowing of declines. For the period of analysis, the earliest breakpoint cycles in the Northern Hig h Figure 2 . 5 . Annual break fraction through time for the (A) Northern, (B) Central, a nd (C) Southern HPA. 54 Plains are mostly positive, indicating slower declines or shifts from decline to recovery. These positive breakpoint cycles may be influenced by irrigation from the Platte River, representing increased recharge from surface flows. A large positive break occurs in the negative breakpoints peaking in the year 2000. This general pattern also holds true for the Central High Plains, but this corresponds with a period of relative r ecovery in the Southern High Plains. A large negative deflection occurred in the Central High Plains in 2010, when nearly 10% of measured wells experienced a negative breakpoint. This signal may reflect the onset of the drought that became severe in 2012, during which all regions of the High Plains Aquifer experienced negative breakpoints despite a simultaneous reduction of irrigated area (Deines et al. 2017). Patterns in the Southern High Plains are unique in exhibiting several cycles of positive breakpo breakpoints. This suggests that positive breakpoints represent different areas that begin recovery at different points in time, rather than the cycles of decline and recovery that are seen elsew here. The spatial distribution of breakpoints confirms this interpretation. Effects of Groundwater Management Areas I then examined the detailed effects of specific GMA s in three states: Nebraska (Natural Resource s Districts), Kansas (Groundwater Managemen t Districts), and Texas (Underground Water Conservation Districts). Colorado, which has comparable 55 irrigation districts, lacked sufficient data for independent spatial analysis, but is included in regional assessment of the Northern and Central High Plains. In Nebraska, the establishment of Natural Resource s Districts was followed mostly by positive breakpoints for the first 3 - 10 years ( Figure 2 . 6 ). The cycle of negative and positive breakpoints continued after about fifteen years of well recovery. NRD s in the western part of the state are more difficult to assess, since fewer long - term water Figure 2 . 6 . Positive and negative annual break fractions over time in select Nebraska Natural Resources Districts (NRDs). NRD outlines are shown at bottom. 56 table records are available. Within the NRDs , there are distinct areas with larger or smaller changes in slope, which may be useful to identify for targeted management prog rams. The breakpoint cycles in the Kansas High Plains are readily identifiable ( Figure 2 . 7 ), but there are more years with simultaneous positive and negative breakpoints than Figure 2 . 7 . Positive and negative annual break fractions over time in Groundwater Management Districts (GMDs), Kansas. GMD outlines are shown at bottom. 57 the other two regions. By comparison, in other regions, breaks within a single year tend to be either positive or negative, so that the net break fraction is almost equal to the total break fraction. This may indicate more heterogeneity within the Central High Plains than in the Northern and Southern High Plains. In Kansas, the establishment of Groundwater Management Districts was associated with a strong signal of negative breakpoints in the west, although the central Groundwater Management Districts showed less change in slope during this period. This may indicate that the Groundwater Management Districts were formed as a response to strong declines in water levels, which is supported by the recent history of Local Enhanced Management Areas (Golden and Guerrero 2017). In western Kansas, the imposition of management districts was followed by a long period of slowing declines. However, the three western Groundwater Management Districts have a few wells with faster declines even in years with overall recovery. Western Kansas also has many wells with large changes in slope at their breakpoints, particularly Groundwater Management District 3. Central Kansas Groundwater Management Districts 2 and 5 are slightly more volatile, with a shorter cycle of decline and recovery, higher peak annual break fractions, and very few wells with breakpoint deflections that are opposite the regional trend for that year. Overall, central and western Kansas tend to behave independently in terms of their breakpoint cycles, likely reflecting differences in geology. Central Kansas has younger, less compact sediments, facilitating recharge. Groundwater Management Districts 1, 3, and 4 in western Kansas ar e also in the driest part of the state. The UWCDs of Texas show a strong regional signal ( Figure 2 . 8 ). Unlike the 58 Nebraska NRDs and the Kansas GMDs , the Texas districts were established over several decades. Figure 2 . 8 indicates that the periods of rapid decline and recovery in Texas tend to be similar in areas with and without a conservation district, a lthough the creation of management districts is generally followed by a period of recovery within Figure 2 . 8 . Positive and negative annual break fractions over time in select Underground Water Conservation Districts (UWCDs), Texas. UWCD outlines are shown on the map at bottom. 59 ten to fifteen years. High Plains Underground Water Conservation District No. 1 has several periods with positive breakpoints, without interceding negative br eakpoints, which may indicate depletion rather than recovery. Changes in slope following breakpoints tend to be largest in the Southern High Plains. Saturated Thickness It would be natural for wells to have positive breakpoints when they, or nearby wells, are no longer yielding enough water for efficient irrigation, and pumping stops. The effect of well yield decline could not be assessed directly using this statistical design. The operative factor for well yield decline is saturated thickness. As a statis tical timeseries, annual saturated thickness is equivalent to annual water table elevation. This is because the saturated thickness is equal to water table elevation minus the aquifer base elevation, which is a constant at a given well location. Saturate d thickness was assessed indirectly by comparing the annual average saturated thickness in a region to the average saturated thickness in wells with breakpoints in that year. Positive breakpoints tend to occur in wells whose saturated thickness is less tha n the regional average, while wells that have negative breakpoints (faster decline) tend to have greater than average saturated thickness. On average, wells with negative breakpoints have 9.5, 4.8, and 4.8 meters greater than average saturated thickness in the year of their breakpoint in the Northern, Central, and Southern High Plains, respectively. For positive breakpoints, the average differences for the Northern, Central, and Southern High Plains are - 5.9 m, - 2.4 m, and - 2.5 m, respectively. This is 60 like ly a reflection that some wells experience positive breakpoints because the aquifer is depleted, causing abandonment of irrigation wells. Areas with greater saturated thickness will support more intensive pumping when irrigation demands increase. Saturated thickness is not a perfect proxy for water yield at any particular well, despite being informative at a regional scale. Many wells are not screened to the aquifer base, particularly older wells and those in the Northern High Plains (Perrone and Jasechko 2 017), and well yield depends on other factors such as the well diameter and pump location (Hecox et al. 2002; Foster et al. 2017). Thus, this does not indicate a threshold for saturated thickness, but rather is a qualitative finding consistent with the hyp othesis that parts of the aquifer are no longer being pumped because they are going dry. Cross - Correlation and Multiple Linear Regression Results A cross - correlation analysis shows that the Palmer Drought Severity Index (PDSI), precipitation, and crop pri ces, as well as yearly changes in these conditions, can have a significant regional effect on net breakpoints. These factors generally have lags between 0 and 7 years prior to changes in net annual break fractions. Figure 2 . 10 shows the correlation between the net annual break fraction, regional drought (Northern and Southern High Plains) or precipitation (Central High Plains), and national crop prices (Northern and Southern High Plains) or change in crop prices (Central High P lains) in 2016 dollars. In the Northern High Plains, the Palmer Drought Severity Index (PDSI) is the climate parameter most strongly correlated with net annual break fractions. There is a 61 strong signal of positive breakpoints following periods of drought. Corn prices are positively correlated with net breakpoints in the Northern High Plains. Water tables tend to recover following periods of high commodity prices. In the Central High Plains, precipitation is more strongly correlated than PDSI with net ann ual break fraction. Well recovery tends to follow wet periods, and increased declines follow dry periods. Increases in commodity prices tend to be followed by a higher fraction of negative breakpoints in the following one to six years. In the Southern Hig h Plains, there is a clear relationship between drought and net annual break fraction. Dry periods are associated with increased rates of decline, and wet periods are associated with recovery, with a peak lag of about four years. This is a longer lag than observed for the Central or Northern High Plains, each of which has a lag of two to three years. In all cases, however, the correlation coefficient is negative, meaning that drought is followed by positive breakpoints (recovery). The results suggest that i n the High Plains Aquifer, drought is the driving factor in water level trends, with relatively little change resulting from wet periods. The multiple linear regression model for each region was based on the crop (corn, sorghum, or cotton) and climatic con ditions (PDSI or precipitation) time series that was found to be most correlated (either positively or negatively) with the net positive and negative breakpoints (Figure 8; Table 1). These regression models had adjusted R 2 values of 0.40, 0.25, and 0.25 fo r the Northern, Central, and Southern High Plains, respectively. 62 Table 1 shows that the net annual break fraction in the Northern High Plains was positively correlated with the price of corn, and negatively correlated with the preceding three years of drough t. Drought is followed by positive breakpoints, Figure 2 . 9 . Multiple linear regression models of net annual break fraction through time. Net annual break fractions are the sum of positiv e and negative breakpoints in a given year, divided by the number of measurements in that year. The gray shading indicates the disparity between modeled and observed net annual break fractions. Model adjusted R 2 values for the Northern, Central, and Southe rn High Plains Aquifer are 0.40, 0.25, and 0.25, respectively. 63 indicating rapid recovery. High crop prices boost positive breakpoints four years later, which suggests that farmers are using more water while crop prices are high, and Figure 2 . 10 . N et number of breakpoints normalized to the number of wells measured in a year , cross - correlated with factors affecting irrigation water deman d. Cross - correlation graphs are shown for the factors that are independently most strongly correlated with net break rate in each region. The most significant weather and commodity price variables for each region are shown . Lag refers to the number of years before the breakpoint occurs . A lag of zero is the correlation in the same year as the net annual break fraction. Dashed horizontal lines indicate the threshold for correlation to be significant at p<0.05. National crop price s are in 2016 dollars. 64 decrease their water use when they are not expecting higher - than - usual compensation. Climate and crop prices explained less of the variation in the Central and Southern High Plains than in the Northern High Plains, but still accounted for about 25% of the break fraction. In the Central High Plains, precipitation and changes in crop price were more significant than PDSI or crop prices. Decreases in precipitation were followed by positive breakpoints after 2 - 3 years, suggesting a slightly longer timeframe fo r recovery after drought, and potentially a lag in the effects of drought on recharge. Declines in crop price are associated with negative breakpoints in the same year, suggesting that farmers may use more water to maximize their yields when their expected return per unit of yield is low. This stands in contrast to the Northern High Plains, where farmers increase their water application when prices are high. In the Southern High Plains, crop prices were not significant in the multiple linear regression mod el. Drought (low PDSI) was followed by positive breakpoints 3 - 4 years later, suggesting a longer lag for recovery than in the Central or Northern High Plains, including greater delayed effects of drought such as low soil moisture. The most significant fact or for the net annual break fraction in the Southern High Plains was the change in drought index. When conditions were wetter than the previous year, there was a high fraction of positive breaks in the same year, which strongly suggests that farmers are wo rking to maximize the benefits of precipitation by pumping less in wet years. The multiple linear regression models for the three regions of the High Plains were used to create an annual estimate of expected net annual break fractions based on 65 climate an d commodity prices. These estimates were translated into expected net annual break fractions before and after the establishment of the 41 GMA s in the study ( Figure 2 . 11 ). Positive deflections were observed in 78% of GMAs (32 out of 41) within ten years of the implementation of the GMA , typically exceeding the net positive break fraction predicted by the multiple linear regression models. Interestingly, based on climatic conditions and crop prices, the linear regression models for the High P lains would predict net negative breakpoints 10 - 15 years after GMA establishment, but this is not observed. Instead the net breakpoints were positive, though small in magnitude. However, there is a marked tendency for water tables to have negative inflecti ons just before and for a few years following enactment of the GMA , also beyond what is expected based on climate and crop prices. There are multiple possible reasons Figure 2 . 11 . Time series of observed, expected, and residual net breakpoints with respect to the year of establishment for grou ndwater management areas. Year 0 is the year in which the management area was established; negative and positive numbers on the x - axis represent the number of years pre - and post - establishment. A. The expected number of net breakpoints, based on crop pric es and moisture conditions, is shown in light gray, compared with the observed net breakpoints in dark gray. B. The residual net breakpoints the difference between expected and observed is shown in red (net negative residual) or blue (net positive residual ). 66 for this association between GMA establishment and faster rates of water table decline. I t could simply be that GMA s tend to be founded during crises, where water tables are declining rapidly. This is supported by the fact that the expected net number of breakpoints based on climate and crop prices is negative (more negative breakpoints than p ositive breakpoints expected) in the years shorty before and after GMA establishment, although the actual number of negative breakpoints is far higher than would be expected from climate and crop prices alone. It is also possible that the announcement of G MA implementation causes farmers to begin using more water. When restrictions are written as proportional or absolute reduction based on current restrictions to be e nacted in the near future. G MA s also generally oversee water use monitoring programs, which incentivize many farmers to report that they are using their full water rights rights that may be taken away in some regions if it can be proven that they are not being used. Conclusions Previous attempts to characterize water level declines have assumed that changes in water level were linear (e.g. Haacker et al. 2016) or logistic (Steward and Allen 2016). This analysis shows that rates of change have significant variations with time, due to climate, crop prices, and other factors that determine the range of management options available to irrigators. Despite the many factors that complicate interpretation of breakpoints, and multiple management strategies used by farmers, 67 cycles of rapid decline and relative recovery occur at the regional level. This is good news, because if the water table did decline monotonically, it would suggest that depletion is inevitable, and that extending the life of the aquifer would be difficult or impossible. Fluctuation in the rate of water table decline suggests an opportunity to enhance the lifespan of the resource through improved management. Segmented regression of well hydrographs provides a new method to evaluate groundwater re sources over time. Breakpoint analysis is particularly appropriate when there is an expectation that a specific change, such as the imposition of management, should lead directly to a recovery of the water table, or at least a decrease in the rate of decli ne. Breakpoint analysis could also provide a rubric for groundwater managers, as an alternative to saturated thickness targets. Rather than striving to limit water table declines, for example, a GMA could adopt a goal of reducing the rate of water table decline. Well hydrographs are relatively easy and inexpensive to obtain compared with water use data, and changes in the water table incorporate the effects of all management and surface processes, in cluding changes in pumping and recharge. Process - based modeling would be necessary to interrogate the causative relationships in the system, and the sensitivity of hydrograph slopes to various conditions. Despite the difficulty in establishing causal rela tionships, climate and crop prices both appear to influence the rates of decline, and therefore saturated thickness. It is remarkable that periods of faster and slower declines follow a strong regional pattern, despite pumping decisions being made at the f arm level. This suggests that a regional scale of management is likely to be effective, since most wells in a region have the same 68 trends, regardless of differences between individual farms. Technology and cropping decisions clearly play a significant role in aquifer declines. It is likely that some breakpoints are related to the adoption of more efficient irrigation systems that both encourage water conservation within a field, and simultaneously encourage irrigation of larger areas. It is also likely that the three regions respond differently to periods with increased rainfall. For example, higher rainfall in Nebraska likely becomes recharge, while in Texas it might improve yields and boost profit enough to encourage farmers to take advantage of the opport unity to drill more wells. It is also notable that the latest cycle of negative breakpoints, starting in the mid - which biofuel use became widespread and corn prices rose. For the most part, regional patterns are consiste nt within a given year, exhibiting either accelerating or slowing declines, with few wells showing breaks opposite the trend of nearby wells. This suggests that the rapidity of decline is largely determined at the regional scale, rather than the farm scale , which in turn suggests that regional management may be effective. Managers may also be interested in examining differences in farming practices between wells with and without breakpoints in a given cycle. With the exception of the Nebraska Natural Resou rce s Districts, most GMA s are formed by irrigators to facilitate groundwater use, rather than limit it (White and Kromm 1995). They invariably perform useful services such as capping wells, and may enforce existing regulations such as water rights priority for senior users and spacing between wells, which may have a small effect on water level declines. However, only a 69 few areas to date have seen large, mandated and enforced decreases in withdrawals. In the future, it is possible that breakpoint analysis co uld help evaluate the effectiveness of programs such as the Sheridan 6 Local Enhanced Management Area in Kansas, in which irrigators voluntarily imposed pumping limits. 70 Chapter 3 Effects of Land Use Change on Water Resources of the Southern High Plains Abstract In 2016, the Conservation Reserve Program (CRP) encompassed about 9,000 square kilometers in the Southern High Plains (SHP). CRP contract acres, which are converted from row cropping to prairie grasslands, are thus a major land cover, representing about 10 % of all land in the SHP. Land cover plays an integral role in water resource management, particularly in the semi - arid Great Plains. In the SHP, plant communities and management play significant roles in the balance between evapotranspiration, infiltratio n, and runoff; these in turn affect the rate of aquifer replenishment or decline. This study uses process - based modeling of the water cycle in conjunction with spatially explicit CRP data to quantify the potential effects of changing land use on deep perco lation in the SHP from 2001 - 2014. Three scenarios were developed: a baseline scenario in which CRP land is retained in shortgrass prairie, a scenario with one third of CRP fields are shifted to agriculture, and an endmember case with all CRP land reclassif ied as agriculture. On average, the model predicts that when existing CRP lands are reassigned as agriculture, average recharge increases by 78% during the growing season (Jun - Sep) and 72% during the non - growing season (Nov - Apr). This translated to a chan ge in the annual average recharge over the whole domain from 1.9 cm to 2.0 cm. Additional recharge in the scenario in which all CRP 71 land was reclassified as agriculture ranged from 1.2 x 10 7 m 3 in model year 2001 to 2.4 x 10 6 m 3 in model year 2005. Introduction and Background Study Area The Southern High Plains (SHP) is a unique natural and agricultural ecosystem in the semi - arid Great Plains of the United States. It includes much of northwest Texas and the eastern margin of New Mexico ( Figure 3 . 1 ) . The major water - bearing unit is the Ogallala Formation. For the purposes of this study, the SHP is extended northward to the valley of the Canadian River, where the Ogallala Formation is de eply incised. Figure 3 . 1 . Location, showing Southern High Plains boundary cut off at the Canadian River, with county outlines. 72 The water cycle in the SHP exhibits a unique balance between recharge, evaporation (ET), and surface runof f to playa lakes. In this semi - arid area precipitation is distributed unevenly in space and time the average annual precipitation was 46 cm/year for the 2001 - 2014 model period, but there is a strong east - west gradient in precipitation. The east received an average of 62 cm/year of precipitation, compared to 29 cm/year in the western part of the study region (PRISM 2004 ). In general there is a bimodal summer/fall precipitation pattern, with a large portion of precipitation from high - intensity events, leadi ng to large interannual variability in rainfall. The Hydrological Landscape Rather than draining to streams that flow away from the study area, runoff in the SHP is concentrated in playa lakes ( Figure 3 . 2 ), ephemeral features which build up Figure 3 . 2 . Playa lakes shown as points. I nset show s actual boundaries of playas. 73 sufficient head to percolate through the unsaturated zone. Approximately 95% of surface water drains internally within the region (Wood and Osterkamp 1987; Osterkamp and Wood 1987). Conversion from shortgrass prairie to dryland agriculture has been observed to increase the percolation rate in interplaya areas (Stanton et al. 2011; Voldseth et al. 2007). Deep percolation from playa lakes depends on local run off; thus, plant communities are a primary mediator of recharge in the SHP. The extensive root system of the native shortgrass prairie is very efficient at intercepting soil water (Scanlon et al. 2005). By contrast, annual field crops such as wheat, cotton , corn, and soybeans are less effective at drawing soil water due to their relative lack of root mass and higher internal root pressure (i.e. higher wilt point; Jackson et al. 1996). Land Retirement Programs In 1986, the Land Retirement Program, part of th e Soil Conservation Service, was replaced by the Conservation Reserve Program (CRP). CRP was codified in the 1985 Farm Bill (U.S. Cong. 1985) with a stated purpose of environmental conservation for the ation programs, CRP requires and subsidizes the replacement of row crops with a mix of grasses similar to the native plant community, rather than simply fallowing. In 1992, the program adopted an Environmental Benefits Index (EBS) to help determine the ren t per acre that the government could offer for a given parcel of land. The EBS prioritized water quality, likelihood of wind erosion, slope, and habitat, among other factors (Hellerstein 2006); 74 however, water resource quantity was not an explicit factor in the EBS. Since 1987, about 10% of the land in the SHP has been enrolled in CRP in any given year (NASS 2017 b ; Figure 3 . 3 ). Purpose of Study In a previous study examining the effects of CRP near the study region (Texas County in the Oklahoma panhandle), Rao and Yang (2010) postulated that CRP could increase recharge compared with dryland crops. Their study describes a correlation of 0.81 between the percent of land enrolled in CRP and average water table recovery. Their results were within a 90% confidence interval, lower than the usual standard of 95% confidence for results to be considered significant. T heir proposed mechanism for this water table increase was that the deep roots of grasses would facilitate, rather than Figure 3 . 3 . (A) Average enrollment in CRP for all counties located at least 85% within the High Plains, and counties in eight High Plains states, but with 15% or less area on the High Plains. (B) Map of percentage of land in CRP in the Southern Hi gh Plains in 2016. 75 intercept, recharge, which is contrary to findings from studies of land transition from undisturbed prairie to dryland agriculture (Scanl on et al. 2005; Stanton et al. 2011; Voldseth et al. 2007). This study examines the effects of the Conservation Reserve Program on the Southern High Plains. I hypothesize that by transitioning crops such as cotton and corn to grasses that are able to inte rcept a larger proportion of precipitation, CRP has reduced the amount of water that runs off to playa lakes as well as deep percolation between playas, and thus reduced recharge to the aquifer relative to the likely amount of recharge if that land were wh olly or partially farmed. This study uses an integrated surface - and groundwater modeling approach to estimate the magnitude of change in runoff, recharge, and evapotranspiration resulting from transitions between grassland and row crops in the Southern H igh Plains. Methods Models and Input Data Spatial CRP data for this study was provided by the National Agricultural Statistics Service (NASS). This data was used to inform scenarios for the Landscape Hydrology Model (LHM) (Kendall 2009, Hyndman et al. 200 7), with integrated MODFLOW (USGS Modular Three - Dimensional Finite - Diffe rence Ground - Water Flow Model; Harbaugh et al. 2000). Table 3 . 1 lists data inputs and source in formation. LHM is a deterministic hydrology model that incorporates a large number of physical surface processes. It estimates the flux of water and energy, incorporating 76 Dataset Name Dataset Source Land use Cropland Data Layer (CDL) www.nass.usda.gov/Research_and_Science/ Cropland/SARS1a.php USDA, National Agricultural Statistics Service, 2016 Cropland Data Layer Surface hydrology National Hydrography Dataset National Wetlands Inventory nhd.usgs.gov/wbd.html USDA, USGS, EPA collaboration, 2014 www.fws.gov/wetlands USFWS, 2014 Weather National Land Data Assimilation System ldas.gsfc.nasa.gov/nldas/NLDAS2forcing.php Cosgrove et al. 2003 Soils Soil Survey Geographic Database sdmdataaccess.sc.egov.usda.gov Soil Survey Staff, NRCS, 2014 Leaf Area Index Moderate Resolu tion Imaging Spectroradiometer modis.gsfc.nasa.gov/data/dataprod/mod15.php Myneni et al. 2014 Irrigation MODIS Irrigated Agriculture Dataset for the United States earlywarning.usgs.gov/USirrigation Brown and Pervez 2014 Elevation National Elevation Dataset Digital Elevation Map nationalmap.gov/elevation.html Gesch et al. 2002 Aquifer boundary United States Geological Survey co.water.usgs.gov/nawqa/hpgw/html/GIS.html Qi 2010 Starting heads United States Geological Survey Groundwater Historical Data waterdata.usgs.gov/nwis/uv/ ?referred_module=gw Haacker et al. 201 6 Observations United States Geological Survey Groundwater Historical Data waterdata.usgs.gov/nwis/uv/ ?referred_module=gw Table 3 . 1 . Data sources for LHM. 77 climate, plant growth, irrigation, soil texture, and other inputs ( Table 3 . 1 ; Figure 3 . 4 ). LHM simulates the full surface water cycle, including deep percolation, evaporation, tr anspiration, and runoff. The MODFLOW groundwater modeling framework is run as a component of LHM. MODFLOW simulates groundwater movement through an aquifer. Its primary inputs include aquifer dimensions such as thickness and specific yield, boundary condi tions, hydraulic conductivity, recharge, and starting heads (see Chapter 4). Hydraulic conductivity was fitted using initial estimates from Cederstrand and Becker (1998 b ), as described in Chapter 4, and starting heads were generated using statistical inter polation from well data, as described in Haacker et al. (201 6 ). Aquifer base elevation was obtained from the Texas Groundwater Availability Model for the High Plains Aquifer System (Deeds and Jigmond 2015). The MODFLOW groundwater model Figure 3 . 4 . Examples of LHM inputs. (A) Percent sand in the top layer of soil, derived from SSURGO data; note some discrepancies at county boundaries. Urban areas within model boundaries are blank. (B) Elevation in meters. The land surface of the High Plains slopes from west to east. (C) Irrigated area from MIrAD - US for 2012. 78 included a smaller subset of the SHP ( Figure 3 . 5 ) to exclude areas that are only weakly connected to one another in the post - development aquifer system. The one - la yer MODFLOW model was initially built as a steady - state model with an approximation of natural recharge from LHM, and predevelopment starting head (Haacker et al. 201 6 ). Estimates of average hydraulic conductivity disagree by several orders of magnitude ( Cederstrand and Becker 1998 b ; Luo and Pederson 2012; Deeds and Jigmond 2015). The steady state model provided reasonable agreement with predevelopment groundwater head observations using a 50% reduction of estimated hydraulic conductivity from Cederstrand and Becker (1998 b ), which in turn is based on a digitization of Gutentag et al. (1984) estimate derived from aquifer matrix lithology. The base of the Ogallala Formation was used as the model bottom, because the Ogallala Formation has a much higher hydraul ic conductivity than the underlying formations Figure 3 . 5 . Model boundaries. Light blue is the LHM model boundary; darker blue is the MODFLOW groundwater model, which overlaps the LHM surface model area. A 5x5 km block of 100 cells is shown for reference to model scale. 79 (Deeds et al. 2015, Deeds and Jigmond 2015) and is thus a reasonable no - flow boundary (Anderson et al. 2015). See Chapter 4 for more details regarding groundwater modeling methods. The linked LHM - MODFLOW mode l was run from 2001 through 2014, the time period for which all input data was available; specifically, MODIS Leaf Area Index data does not begin until 2001, and NLDAS weather data was available through 2014 at the time of model construction. A 500 m grid was used for both the LHM and MODFLOW models. This grid size is approximately the scale of an agricultural field in this region. MODFLOW is called as a module within LHM. The surface model apportions precipitation and irrigation between runoff, ET, and re charge, based on climate and landscape parameters. LHM saves deep percolation data in a format readable to MODFLOW, incorporating a delay in recharge arrival at the water table of up to six months based on a one - dimensional empirical model (Kendall 2009). MODFLOW simulates groundwater head change and discharge to drains, and total water balance. Spatial data for CRP was obtained through a research agreement with USDA. The spatial CRP data for this study was only complete within the model area for the year 2016, thus scenarios were based on 2016 land use and CRP coverage. These data were provided as polygons that occupy portions or all of individual Common Land Units (CLU). Spatially explicit CLU and CRP data cannot be shared freely, as stipulated by the 20 08 Farm Bill (U.S. Cong. 2008). All information directly related to CRP lands are presented as county - level or higher aggregate data, preventing identification of enrollment status of individual fields or landowners. 80 Scenarios Three land - use scenarios were developed according to the land use class assigned to land enrolled in CRP in 2016. A baseline scenario was run using the 2016 Cropland Data Layer ( USDA 2016 ) with CRP assigned as grassland ( Figure 3 . 6 ; Table 3 . 2 ). The two alternative scenarios were created by changing the land use class of fields enrolled in CRP. O involved an intermediate level of land use transition, with 30% of fields in CRP reclassified as agricult ure. This corresponds with the estimated rate of conversion to agriculture upon the expiration of CRP contracts (Morefield et al. 2016), and is similar to what would be expected if CRP had terminated in 2001, at the beginning of the model Figure 3 . 6 . Results for baseline analysis, averaged over the model period, 2001 - 2014. ( A) Runon to cells, in meters. Most of the model area is internally drained. (B) Recharge, i.e. percolation below the root zone. 8 1 Table 3 . 2 . CRP land was primarily identified as grassland in the NASS Cropland Data Layer (CDL) for the available years in the study area (Myneni et al. 2014). period. No irrigation was applied to 2016 CRP land regardless of the scenario, since CRP payment rates are based on rent for dryland farming, and land that transitions from CRP to agriculture is unlikely to be irrigated except at the beginning of the contr act as plant cover is being established. 70% of land identified as CRP was grass according to the 2016 Cropland Data Layer (CDL). The rest was primarily shrubland ( Table 3 . 2 ). LHM includes nine land cover categories: urban, agriculture, grassland, shrubland, deciduous, coniferous, water, wetland, and barren. Land use in LHM is used to determine root distribution, canopy height, and canopy conductance. For the three scenarios, all 2016 CRP land was reclassified as either grassland or agriculture. Root distribution for LHM was modified from Jackson et al. (1996). LHM does not incorporate specific root mass, but rather uses a representative proportional distribution with depth for each land use type. This root distribution, which uses an exponential model to estimate the proportion of root ma ss at a given depth (Jackson et al. 1996), was adjusted so that grass roots would extend more deeply into the soil, since 2014 CDL 2015 CDL 2016 CDL Average Grassland or Pasture 72.27% 70.83% 71.41% 71.50% Shrubland 11.22% 13.70% 11.73% 12.22% Misattributed as Cropland 6.13% 5.00% 5.90% 5.68% Hay or Alfalfa 3.31% 3.75% 4.16% 3.74% Fallow or Idle Cropland 2.75% 2.57% 2.61% 2.64% Developed Open Space 2.61% 2.56% 2.14% 2.44% Water or Wetland 1.10% 0.97% 1.15% 1.07% Other 0.98% 1.05% 1.35% 1.13% 82 the predominant crops in the study area ( Figure 3 . 7 ) have shallower rooting systems than native grasses. Minimum canopy height was 0.5 m for grass and 0.0 for crops, while both had a maximum canopy height of 1 m, scaled according to seasonal LAI. Canopy conductance, a measure of how much water a plant uses relative to the vapor pressure deficit, was also lower for grassland (0.007) than for agriculture (0.0116; J.M. Chen et al. 2005). The result is that potential evapotranspiration (PET), calculated hourly by a modified Penman - Monteith equation (Kendall 2009), is higher in cropland versus grasses for equivalent LAI. However, grasslands green up earlier, shifting the seasonality of ET and extending the grassland growing season. Also, grasslands are Figure 3 . 7 . Land use in the SHP in 2016, CDL. Land was reclassed in LHM as nine land types, including agriculture, grassland, shrubland, and barren. 83 capable of withdrawing water from the soil at a higher rate than crops, given by their greater root mass ( which translates to reduced root resistance). Typically, permanent wilting point is defined as the volumetric soil moisture at which plants can no longer draw water from the soil, and thus permanently wilt. Measured as a matric potential in meters of water , typical wilt points for native rangeland are approximately - 600 m to - 800 m, while the wilt point for crops is about - 150 m, so crops are far less effective at removing soil moisture in dryer conditions (Scanlon et al. 2005). Wilt point was modified from standard LHM default values in order to account fo r lower wilt points in grasses. Results and Discussion Baseline Model Results In the baseline scenario, deep percolation was roughly consistent with estimates for average recharge in the Southern High Plai ns (Stanton et al. 2011). Average wetland runon was 0.11 m/yr, meaning that on average, wetlands received a total of 11 cm of water (equivalent total depth) over the course of a year. Wetlands contributed 3.4 times as much recharge in proportion to upland areas, although wetlands make up less than 5% of the model area. Average recharge flux was 1.9 cm/yr over the entire model area in the baseline scenario. This is similar to predevelopment estimates of 1.5 cm/yr (Wood and Sanford 1995). This is lower than e stimates from Stanton et al. 2011, which ranged from 0.25 - 3.6 cm/yr for non - playa areas, and 9.1 cm/yr in playa lakes, for an overall average of 2.8 cm/yr. In the baseline scenario, where land enrolled in CRP was 84 average recharge in CRP areas was 1.2 cm/yr, about 40% less than in the scenario with all CRP fields reclassified as agriculture. Simulated groundwater levels agreed with observed water table elevations within a margin of error of two standard deviations ( Figure 3 . 8 ). Simulated water tables are generally less than five meters in error compared with measurements. A starting head bias is evident, likely because starting heads are derived through statistical interpolation of water table measurements (Haacker et al. 2016), and suggests that the model should be run for a few years to generate the initial water table. Surface water flows were not available for model validation, because there are no channels in the study area with sufficient flow to warrant a stream gage. Figure 3 . 8 . ( A) Simulated versus observed groundwater levels for the baseline run, which did not incorporate informatio n from the spatial CRP dataset. ( B) The baseline run has groundwater elevation estimate residuals within a margin of error of observations (two standard deviations around the mean, shown). 85 CRP Transition Scenarios In land use class, overall average annual recharge was about 1 mm less than in the the most - likely scenario for termination of the Conservation Reserve Program (30% of Figure 3 . 9 . Average deep percolation, cm, for the same CRP land area in three scenarios: baseline (CRP is grassland), most likely scenario if CRP were ended (1/3 of CRP area reassigned as agriculture), and all CRP land area converted to agricult ure. Recharge is consistently higher in farmland. ( A) Average deep percolation during the growing season (Jun - Sep), land transition scenarios minus baseline. ( B) Average deep percolation during the non - growing season (Nov - Apr). 86 CRP fields transitioned to agriculture), average recharge was intermediate between the baseline and full - reclassification scenarios. Upland recharge, from deep percolation between playa lakes, was consistently higher when CRP areas were reclassified as agriculture ( Figure 3 . 9 ). The effect of land use transitions varied from year to year, although recharge was consistently higher in scenarios where CRP fields were reclassified as croplands. In the non - growing seas on (Nov - Apr; Figure 3 . 9 a), recharge in the scenario with all CRP land reclassified as agricultural use was 78% higher than recharge in the same area when CRP was class ified as grassland. During the growing season (Jun - Sep; Figure 3 . 9 b), recharge for land reclassified as agriculture was 72% greater. In all scenarios, recharge was co ncentrated during wet years such as 2005, 2007, and 2011, particularly during the growing season. When CRP land was reclassified as agricultural, the difference between dry and wet - year deep percolation increased relative to the other scenarios. Conclusio ns Owing to the lack of external drainage, the water cycle of the SHP is uniquely sensitive to changing land cover. Transitioning from grassland to dryland agriculture has the potential to increase recharge. When sustainability efforts are focused on a sin gle metric environmental health, there is eventually a tradeoff with other aspects of the ecosystem (Hoekstra and Wiedmann 2014). In the case of CRP adoption, rehabilitating the landscape may cause recharge to decrease to something closer to predevelopment levels, which may hamper recovery of the water table outside of the 87 irrigation season. Model experiments suggest that land use change has the greatest effect during wet years ( Figure 3 . 9 ). Absolute differences in deep percolation were less pronounced during the non - growing season, likely due to the dormancy of grasses in the winter. Unlike watershed - based surface hydrology models, LHM estimates the flux of water on a c ell - by - cell basis, making it particularly suitable to the Texas Panhandle due to the predominance of internally drained areas. Results were consistent with observations for cropland recharge and water table elevation. Irrigated land is unlikely to be enro lled in the Conservation Reserve Program. Federal rents are based on the rates for dryland rental. However, some states have invested in the complementary Conservation Reserve Enhancement Program (CREP), a state - federal partnership. CREP focuses on retirem ent of irrigated farmland, and pays commensurately higher rents. While it provides for limited enrollment compared to CRP, CREP has the potential to improve water table recovery rates more effectively than CRP. The pumping rate exceeds recharge by approxim ately an order of magnitude in the SHP, so cessation of pumping will more than offset a decrease in recharge, as long as it is not offset by increased pumping elsewhere. However, CREP is used mostly in areas where pumping affects surface water flows, such as near the Platte River in Nebraska. Neither Texas nor New Mexico had a CREP agreement with USDA during the model period. Another study on the HPA, in Finney County, Kansas, found that farmers were much more likely to enroll in CRP in areas with little s aturated thickness, even though 88 fields enrolled in CRP are dryland, not irrigated (Nellis et al. 1996, Egbert et al. 1998). This suggests that farmers and managers may hope to facilitate water table recovery through CRP enrollment; otherwise, there would b e no difference between rents for dryland above the High Plains Aquifer and dryland in nearby counties, or dryland with a lot or a little underlying saturated thickness. Because grassland enrolled in CRP exhibits less recharge than dryland agriculture, the program may in fact may be more useful in areas where arable land is becoming salinized due to waterlogging from rising water tables. CRP involves a large federal allotment. $15.4 billion was paid directly to farmers in the High Plains states (Texas, New Mexico, Oklahoma, Kansas, Colorado, Nebraska, Wyoming, and South Dakota) from 1986 - 2016 (NASS 2017 b ). County rents on the High Plains averaged $48.25 per acre, while counties that are outside the boundaries of the HPA but within the same eight states avera ged $43.13/acre. CRP involves a number of improvements in ecosystem services, but this may occur at the expense of recharge without additional interventions designed to increase deep percolation, such as wetland resto ration. 89 Chapter 4 Groundwater Decline and Wate r Resource Management in the Texas Panhandle Abstract Groundwater management on the Southern High Plains of Texas and New Mexico presents multiple scientific and political challenges. The agricultural economy of the region is anchored by the High Plains Aquifer (HPA), which is experiencing significant groundw ater level declines due to decades of groundwater mining. The region lacks alternative water sources, thus water use planning is necessary to extend the aquifer lifespan. This study uses an integrated surface hydrology and groundwater modeling approach to investigate the effects of water conservation in the most heavily irrigated portion of the Texas Panhandle. I contrast reduction in irrigated acreage (termed the extensive margin in economics) with reduction in application intensity (intensive margin) usin g process - based hydrologic modeling. If pumping had ceased in 2001, our 14 - year simulations show that water levels in the Southern High Plains would have recovered by an average of 0.062 m/y, whereas with no pumping limits, the aquifer declined by about 0.46 m/y. This suggests that the aquifer has about a 7:1 ratio of depletion to recovery including return flow. Under the business - as - usual scenario, approximately 30% of irrigable land was lost in 14 years. Two alternative management scenarios were run, ea ch of which decreased irrigation withdrawals by about 50%. Decreasing irrigated land by 50% or imposing application 90 limits of 25 cm on existing irrigated land reduced the loss of potentially irrigable area to 15% and 15.5%, respectively. Limits to water ap plication were most effective in areas with a high density of irrigated land, whereas reductions in the extent of irrigated land were more effective in areas that already had a low density of irrigated acreage. Limits to the extent of irrigated land improv ed the ability of remaining irrigated areas to adjust to drought. Introduction and Background Study Area Irrigation on the Southern High Plains (SHP) is both essential to the regional economy and unsustainable in its current form. Prolonging the life of the aquifer can increase the total crop yield from irrigation over the long term, along with total profit (Foster et al. 2014). The study area is predominantly in Texas, which is the only state in his doctrine asserts that groundwater is private property and not subject to most limits for use. Texas was also the first state in the High Plains region to create water conservation districts, starting with High Plains Underground Water District Number 1 ( Figure 4 . 1 ), which created a basis for water planning; its name was later changed to High Plains Underground Water Conservation District No. 1. High Plains UWCD 1 has adopted four ten - year management plans, starting in 1998 (High Plains UWCD 1998). Since 2005, these plans were modified and adopted in conjunction with the designated Groundwater Management Areas (GMAs) of Texas, 91 each of which encompass es multiple Und erground Water Conservation Districts. The GMAs are administered at the state level, with the help of the Texas Water Development Board (Mace et al. 2008). Groundwater management plans must use the estimates of available groundwater from the Groundwater Av ailability Models (GAMs) for the local aquifers (Deeds 2015; Deeds et al. 2015). Each UWCD creates its own groundwater plan, but there was no effective mechanism of enforcement during the model period. Water for Irrigation in the Southern High Plains This study focuses on the impacts of varying irrigation scenarios on the Figure 4 . 1 . Boundary of the Southern High Plains, showing groundwater conservation districts in the Texas Panhandle. The groundwater model overlies about half of the High Plains Underground Water Conservation District No. 1, which was instituted in 1951 and adopted its first Groundwater Management Plan in 1998. 92 groundwater within a portion the Southern High Plains ( Figure 4 . 1 ). This area is one of the most heavily irrigated subregions of the High Plains; although it encompasses less than 5% of the total aquifer area, it accounts for almost 10% of the irrigated land in the HPA based on MODIS data (Myneni et al. 2015; Qi 2010). Unlike other parts of the High Plains, this area possesses no perennial streams. Virtually all irrigation water in this region is from the High Plains Aquifer. The groundwater model area overlies a paleovalley, a bedrock low that existed prior to deposition of the Ogallala Formation, which resulted in thicker water - bearing sediments ( Figure 4 . 2 ). The groundwater model is bounde d by an escarpment to the east, with historically flowing springs (White et al. 1946) that have mostly stopped flowing. Less saturated thickness, and less irrigation, occurs to the north, south, and west of the groundwater model area ( Figure 4 . 3 ). The predominant direction of groundwater flow is to the southeast ( Figure 4 . 4 ). The m odel area is almost entirely Figure 4 . 2 . Thickness of Ogallala Formation sediments in the Southern High Plains. The black interior border is the extent of the groundwater model in this study. 93 within Groundwater Management Area 2 and the High Plains Underground Water Conservation District No. 1. This area is among the least sustainable in the HPA (Scanlon et al. 2012; Haacker et al. 201 6 ), and is one of the parts of the HPA that could most benefit from conservation measures, based on average saturated thickness and estimated hydraulic conductivity (Foster et al. 2017). This area of the Texas Panhandle is vital to the regional economy, ma king water planning and conservation essential. Irrigation water can be conserved in two ways: by Figure 4 . 3 . Irrigated area in the Southern High Plains in 2007, which is the coverage used as a baseline in the LHM - MODFLOW model. Most of the irrigated area overlies location s where the estimated saturated thickness in 2001 is greater than 9 m. The black interior border is the extent of the groundwater model in this study. 94 limiting the amount of water applied to a given field (referred to in economics as the intensive margin), or limiting the amount of irrigated land (similarly , the extensive margin). Ultimately, all interventions to reduce the use of groundwater are effective in one or both of these dimensions. For example, adopting more efficient irrigation technology will limit water applied per hectare, while the enforcement of minimum well spacing requirements can limit the number of irrigated hectares. Methods Model Description Four scenarios were designed to investigate irrigation water conservation techniq ues using the Landscape Hydrology Model (LHM; Hyndman et al. 2007; Kendall 2009; Luscz et al. 2017, Wiley et al. 2010) along with its integrated groundwater module Figure 4 . 4 . Closeup of the MODFLOW groundwater model boundary, showing specified head cells (dark blue) and 50 meter conto urs for water table elevations in 2001, in meters above sea level. The High Plains Underground Water Conservation District No. 1 is shown in gray. 95 (MODFLOW). The model was built with a 500 - m grid spacing for both the surface and groundwate r components and run hourly from 2001 - 2014. The root zone of the surface component extends down to 2.5 meters depth and was divided into seven layers, increasing by factors of two starting with a top layer that is 2.5 cm thick. Beneath the root zone is a d eep unsaturated zone that extends to the water table. Water that percolates below the root zone becomes groundwater recharge after a delay corresponding to the thickness of the unsaturated layer. The fully saturated groundwater system was simulated using a single layer, and weekly stress periods. LHM requires seven main types of input data: soils, hydrography, topography, land use, weather, leaf area index (LAI), and groundwater model inputs. Soil data was derived from the SSURGO dataset. Hydrography came from the National Hydrography Dataset (NHD) and the corresponding w atershed b oundary d ataset. Topography data 30 meters in this area. Land use data came from the National Land Cover Datasets (NLCD), for 2001, 2007, and 2011, and from the MiRAD dataset (Myneni et al. 2014). LHM was driven using hourly climate data from the NLDAS - 2a forcing dataset, a model/data reanalysis product at 8 km resolution. Additionally, MODIS LAI data were used to drive vegetation change acr oss the landscape. See Chapter 3 for a table of LHM inputs. Groundwater model inputs are summarized in Table 4 . 1 , and included specific yield from the USGS (Cederstra nd and Becker 1998b), Ogallala Formation base elevation map from the Texas Groundwater Availability Model for the High Plains 96 Aquifer (Deeds and Jigmond 2015), and starting heads for the transient model estimated using methods d escribed in Haacker et al. (2016 ). A more detailed description of the hydraulic conductivity calibration is described below. The MODFLOW model uses a combination of constant - head and no - flow boundary conditions ( Figure 4 . 4 ), which were informed by the direction of groundwater flow. No - flow boundary conditions are located adjacent to areas with little or no saturated thickness, and along flowlines that are normal to th e equipotential lines shown in Figure 4 . 4 . Constant head was assigned along the northwestern boundary. The starting head conditions in 2001 were developed by interpola ting measured water Dataset Name Dataset Source Land Surface Elevation National Elevation Dataset Digital Elevation Map nationalmap.gov/elevation. html Gesch et al. (2002) Accessed from NED 2011 Hydraulic Conductivity Digital map of hydraulic conductivity Modified from Cederstrand and Becker 1998b water.usgs.gov/GIS/metadata/ usgswrd/XML/ofr98 - 548.xml Base Elevation Texas Groundwater Availability Model Deeds and Jigmond (2015) Observations Steady State - Predevelopment dataset Transient Water table observations 2001 - 2014 Steady State Methods from Haacker et al. 2016 Transient Data from USGS Groundwater Data for the Nation, waterdata.usgs.gov/nwis/gw Specific Yield Digital specific yield map McGuire et al. 2012 water.usgs.gov/GIS/metadata/ usgswrd/XML/sir12 - 5177_hp_sp_yield.xml Starting Heads Steady State - Predevelopment water table Transient Water table observations 2001 - 2014 Methods from Haacker et al. 2016 Table 4 . 1 . Data sources for the MODFLOW groundwater model. 97 table (Haacker et al. 201 6 ). To delineate the constant head boundary, the absolute magnitude of annual change was summed for the available annual datasets, 1935 - 2016. The constant - head boundary in the northwest was placed along an area where fluctuations in saturated thickness estimates have been minimal relative to nearby areas in the seven decades of irrigation prior to the model run. Outside of the model boundary to the northwest, there is less measurement data available due to low sa turated thickness, while further inside the model area, there have been large declines. Future work will incorporate a model spinup, running the model for several years prior to the test period, to generate physically based starting heads, and to enable ad justment of the northwest constant head boundary and incorporation of more of the visible paleovalley in Figure 4 . 2 . The extent of active cells was decided based on saturated thickness and irrigated extent ( Figure 4 . 3 ). The southern edge of the boundary is an area with very little saturated thickness, and the east is bounded by the escarpment that forms the edge of the High Plains, dissected with drains that were defined using flow accumulation estimates from the land surface elevation. When MODFLOW is run as a component of LHM, it uses the drains from the LHM surface model. While more active model area could have been included to the north and south, there would be limited hydrogeologic connection outside the paleovall ey. The selected model area represents the majority of irrigated land in the SHP and is the area of the HPA that is most heavily irrigated with groundwater, both in terms of application rate and density of groundwater irrigated acres. It is likely the most vital area of the HPA for future enhanced groundwater 98 management. Hydraulic Conductivity Calibration The MODFLOW model was built using MODFLOW - 2000 (Harbaugh et al. 2000), initially as a st eady - state, predevelopment model, and then incorporated into LHM for transient simulations from 2001 - 2014. Hydraulic conductivity ( Figure 4 . 5 ) was an important parame ter in the model, but estimates range over several orders of magnitude for the region, and distributions are not consistent among estimation techniques (Cederstrand and Becker 1998 b ; Deeds and Jigmond 2015; Luo and Pederson 2012). A modification of the USGS estimate of hydraulic conductivity, which is based on a lithologic assessment from the High Plains Regional Aquifer Study (Gutentag 1984), provided a reasonable fit between simulated and measured predevelopment heads for the steady - state model ( Figure 4 . 6 ; Haacker et al. 201 6 ). The USGS estimated hydraulic Figure 4 . 5 . ( A) Saturated hydraulic conductivity from the Texas Groundwater Availability Model of the Ogallala Aquifer (Deeds et al. 2015). ( B) Hydrau lic conductivity used in the MODFLOW model for this study, modified from USGS estimates (Cederstrand and Becker 1998b). 99 conductivity across the region based on the lithology of the bedrock at a location. The midpoint of the range in each location was used as the initial estimate of hydraulic conductivity calibration. When that estimate proved too high for steady - state model convergence, other modifications of hydraulic conductivity were attempted including estimates from Luo and P ederson (2012) and calibrated hydraulic conductivity values from the Texas GAM (Jigmond et al. 2015). The model with the low end of the range provided by the USGS for each lithology also failed to converge. However, a 50% reduction in hydraulic conductivi ty from the USGS lithologic estimates offered good agreement with water table observations ( Figure 4 . 6 ). It was less prone to flooding than steady state solutions with lower estimates of hydraulic conductivity, and did not fail to converge due to dry cells, which occurred at higher estimates o f hydraulic conductivity. Simulated Head, m Observed Head, m Figure 4 . 6 . Observed versus simulated head for the steady - state model used to develop boundary conditions and hydraulic conductivity parameters for MODFLOW. 100 Irrigation Scenarios Four scenarios were run. The baseline scenario adheres to current practices and land use, with no strict irrigation limit, and with irrigated area assigned using the 2007 MODIS - MIrAD dataset (Table 1). The 20 07 MIrAD dataset was selected because it has the greatest irrigated extent out of the three available MIrAD years (2002, 2007, and 2012), as the effects of drought in 2002 and 2012 likely reduced irrigated area (Deines et al. 2017). Two intermediate scenar ios were designed to investigate the effects of irrigation reduction through decreased irrigated land and decreased water applied per hectare. The irrigation intensity - lim i ted scenario kept the irrigated acreage from 2007 MIrAD - US ( Figure 4 . 3 ), but set a limit of 25 cm for irrigation water applied. This is similar to irrigation limits that have been set in the Local Enhanced Management Area (LEMA) Sheridan County SD - 6 within Kansas Groundwater Management District 4 (Kansas Dept. of Agriculture 2013) and some Natural Resource Districts in Nebraska, which have capped pumping at 140 cm (55 in) over 5 years, thus averaging no more than 28 cm (11 in) per year. 101 The irrigation extent - limited scenario allowed the same maximum per - hectare irrigation total (0.7 meters, higher than the seasonal water demand for almost all cells), and reduced the irrigated area by 50 %. This was accomplished by removing half the irrigated area, by overlaying a 1x1 km grid in which every other cell was removed from Figure 4 . 7 ). This was a large 500 - m LHM - MODFLOW grid, but also small enough that the overall distribution of irrigation across the area remained sim ilar. This simulates the effect of reduced irrigated area, perhaps due to greater well spacing, which might occur due to decreased well Figure 4 . 7 . Irrigation scenarios used in this study. Left: irrigation extent in 2007, used in the baseline and per - hectare application limited scenarios. Right: closeups showing irrigated extent in the three irrigate d scenarios. The fourth scenario had no irrigation. 102 yields or to top - down management. A final scenario was run with no irrigation, to give a sense of how much recovery woul d be expected in the absence of pumping. All agriculture was dryland; land use classes were not changed in LHM. LHM automatically scheduled irrigation based on simulated soil moisture. An irrigation event was triggered when a soil had less than 33% of its plant available water remaining in the root zone, based on soil properties and water content. For each irrigation event, 2.54 cm of water was applied to the entire cell. When a soil has 33% of its plant available water, the plant can still take up water, b ut the soil is fairly dry. This was considered a reasonable proxy for farmer decision - making. Well yield was set at 136 m 3 /hr (600 pm ) based on typical well yields for the prevalent hydraulic conductivity and saturated thickness in the SHP (Hecox et al. 20 02). Agricultural planting and harvesting dates were selected to match the typical day of year for cotton and corn (planting: day 133, harvesting: day 290). No irrigation occurred outside of the growing season, or when leaf area index was below a leaf area /ground area threshold of 0.0005. During a LHM simulation, the surface components create inputs for a subsequent transient run of the MODFLOW groundwater module. MODFLOW acquires inputs for recharge and pumping (as negative recharge) for all stress period s using the surface hydrology estimated in LHM. Only water for irrigation was withdrawn from the aquifer, as water use for municipalities and industries is only a small part of the water budget (<5%) in this area, and is largely not consumptive. 103 Results and Discussion Water Tables and Saturated Thickness Groundwater levels for the baseline scenario were within two standard deviations (about 5 meters) of observed water table el evations. Residuals get closer to zero as the scenario progresses, suggesting that starting heads might be biased ( Figure 4 . 8 ). In addition, to improve model stability , initial saturated thickness was set at a minimum of five meters. Since the aquifer is very thin in some places, this may have been enough to bias model outcomes. The water table declined more slowly in both limited irrigation scenarios relative to the ba seline scenario ( Figure 4.9 ); however, this reduction did not occur equally across the study area. In the areas with the highest Simulated Head (m) Observed Head (m) Year Figure 4 . 8 . Water table observations versus simulated water tables, 2001 - 2014, for the baseline scenario (no limits to extent of irrigated land or intensity of water applied). Residuals are two standard deviations around the mean. 104 Figure 4 . 9 . (A) Saturated thickness for the baseline scenario, with full irrigation extent and intensity. (B) Saturated thickness in the no - irrigation scenario minus the baseline. (C) Saturated thickness when irrigated area is limited minus the baseline. (D) Saturated thickness when irrigation intensity is limited minus the bas eline. Saturated thickness was greater than the baseline in all irrigation limited scenarios, but the spatial distribution of saturated thickness was different depending on the density of irrigation in a locality. All saturated thicknesses are for 2014, at the end of the model runs. density of irrigated land, limited - extent was not as effective for water table recovery as the limited - intensity scenario; Figure 4 . 11 ). I n areas where irrigated fields are widely spaced, further reductions to irrigated land have a greater effect than water application limits. Conversely, where most of the landscape was irrigated, water tables remained higher when application limits were p ut into place. This is due in part to the seasonal 105 effects of pumping for the limited - extent and limited - intensity scenarios. When irrigated land was restricted (limited - extent), farmers watered throughout the growing season, while limit ed application (limited - intensity) led to farmers hitting their irrigation quotas fairly early in the season ( Figure 4 . 10 ). Future work will control for this, since wa ter - restricted farmers are more likely to schedule irrigation according Figure 4 . 11 . ( A) Head difference for the limited - extent scenario and the limited - intensity scenario. Blue indicates that the limited - extent scenario outperformed the limited - intensity scenario, while red indicates that the limited - intensity scenario performed better in terms of water table elevation at the end of the model period. ( B) Average water applied per year in the baseline (unrestricted irrigation) scenario. Figure 4 . 10 . Mean monthly irrigation water applied by scenario. 106 to plant growth stages. More sparsely irrigated areas also tend to occur on marginal agricultural soils, some of which are poorly drained; reduction in irrigated land (limited - extent) had a greater effect in these areas. Shifts in individual crops were not considered due to the complexity of such scenarios. LHM currently uses the same parameters (e.g. root distribution, canopy properties, and albedo) for all agricultural land uses, re gardless of cultivar. In the No Irrigation run, the average water table recovery over the course of the simulation was 0.93 m, but this varied by location, as shown in Figure 4.9 b. Some areas continued to see slight groundwater declines without irrigation , and some areas saw large increases in water levels. Some areas of the model dewatered almost immediately, suggesting a need for a groundwater model spinup period. Where the aquifer was already too thin to support irrigation at the start of the model run, comparatively little increase in saturated thickness occurred. In the most heavily irrigated parts of the study area, water tables in this scenario rebounded as much as 20 meters over the course of the model run. This assumes that the aquifer material has not compacted enough to reduce effective porosity. Irrigable Area Although limiting irrigated land and limiting water application intensity were both effective in maintaining total irrigable area compared with the baseline scenario, they exhibited differe nt patterns in time and space. In drought years, when water application per hectare was limited, farmers were unable to adjust water use. While 107 limits to irrigated land led to slightly more usable area throughout the scenario, drawdown was greater in drought years such as 2011 - 12 when application intensity was not limited. Overall pumping was realistic, based on metere d and estimated irrigation data from the Texas Water Development Board for crops grown in the Southern High Plains (Turner et al. 2011). In both reduced irrigation scenarios less irrigated land and limited water application intensity pumping was reduc ed by ~50% relative to the baseline, although the seasonal pattern was different between the two limited irrigation scenarios ( Figure 4 . 10 ). This seasonal shift is re flected in the difference between saturated thicknesses in the baseline and the limited irrigation scenarios ( Figure 4.9 ). As shown in Figure 4 . 12 , these reductions in pumping would be somewhat effective in prolonging the useable period of the aquifer, but they are far more profound than any politically or economically feasible reduction in the near term. In the final year of the model run (2014), 2 7% of the model area was no longer usable for irrigation, based on a threshold Figure 4 . 1 2 . Irrigable area over time, based on a 9 m threshold for saturated thickness. 108 of 9 meters of saturated thickness for sufficient well yield (Hecox et al. 2002). In the limited irrigated land and application intensity scenarios, about 15% of land was no lon ger usable for irrigation, and more land was becoming usable each year relative to the baseline. In the baseline scenario with business as usual pumping, 4.6% of the model area (900 km 2 ) had insufficient saturated thickness after 14 years (0 - 9 m saturated thickness in Figure 4.9 ). By comparison, a large reduction in irrigation led to less dry aquifer area, to similar degrees in the extent - limited scenario (1.8% dry) or intensity - limited scenario (1.9% dry). In the scenario without irrigation, irrigable are a increased slightly, from 87% in 2001 to 90% in 2014 ( Figure 4 . 12 ), and saturated thickness was greater than zero everywhere by the end of the model period. The aver age water table decline for the baseline scenario was 0.46 m/yr. The average decline was equivalent for extensive and intensive marginal irrigation reduction, at about 0.2 m/yr. Conclusions The Texas Panhandle is taking strides toward more sustainable wate r management by adopting management goals that are both economically acceptable and environmentally feasible, known locally as Desired Future Conditions (Mace et al. 2008). But at a fundamental level, its groundwater resources will not suffice to support i rrigated agriculture on a broad scale using current practices. There is too much imbalance between recharge and crop water demand (Scanlon et al. 2006). The usable period of the aquifer will be extended by decreasing irrigation, either through irrigated 109 a rea or amount of water applied per acre. It is environmentally and economically important that this occur via careful resource management, rather than failure of wells as the aquifer is exhausted (Foster et al. 2017). Both the limited intensity and limited acreage scenarios were unrealistically strict, but each had aquifer areas that were no longer usable for irrigation by the end of the 15 - year model period. This illustrates the fact that even large reductions in irrigation will not be sufficient for aquif er recovery. It is very significant that the areas with the greatest proportion of irrigated land cover had the largest saturated thickness response to application limits, whereas areas with sparse irrigation had more response to further decreases in irrig ated land. These forms of management are not mutually exclusive, and each could be targeted at the areas where they will likely be most beneficial. Texas created the first groundwater management conservation districts in the High Plains (see Chapter 2), an d it is possible that relatively restrictive LEMA - like management areas may become a popular and effective style of resource stewardship despite the statewide Rule of Capture. As discussed in Hornbeck and Keskin (2014), adjustments to irrigation occur diff erently on the intensive and extensive margin, depending on the timescales involved. When crop water demand or water availability change, the intensive margin is the first to adjust crops are watered more or less. Over time, if climate, aquifer, and mark et conditions continue, the pattern of irrigated land will evolve in response. While Hornbeck and Keskin (2014) studied the historical effects of increased access to water in the HPA, similar adjustments are likely to occur as the aquifer is depleted. An 110 u nderstanding of these factors will allow targeted management to prolong the life of the resource and offer maximum relief to farmers, who often have large sums invested in irrigation equipment. These results could provide a scientific basis to set improve d desired future conditions for the Texas Underground Water Conservation Districts. Even an average decline of 0.2 m/yr across the model area is clearly not sustainable. If the most vulnerable areas of the aquifer are left unprotected, then maintaining eve n 50% of saturated thickness is not sufficient. The easiest places to maintain saturated thickness are those that have already been depleted below a threshold in which large volume pumps are functional. As irrigation technology improves, it tends to chang e the balance between intensive and extensive irrigation (Li and Zhao 2016). As these scenarios show, substantial changes to either of these margins can cause large changes in the total water balance. This suggests that intensive versus extensive water use may have a greater effect on other environmental indicators, such as habitat availability and farm effluent. Hendricks and Peterson (2014) point out that the marginal value of water decreases with increasing application, which suggests that reductions on the intensive margin are more cost - effective. Currently, groundwater management plans are evaluated by the Texas State Water Board using a MODFLOW model (Deeds and Jigmond 2015). The inclusion of a landscape model such as LHM could provide improved cons ervation opportunities, by giving a more complete picture of the water cycle and the apportionment between 111 recharge, runoff, and ET; this could help provide integrated conservation plans including maintaining playa wetlands. The water plans adopted by Grou ndwater Management Area 2, starting in 1998, have generally included the provision that 50% of saturated thickness will remain in the aquifer in 50 years, known as the 50/50 Management Standard (twdb.texas.gov). However, as of 2016 when the last standard w as adopted, the desired future condition for GMA 2 is seven to eight meters (23 - 28 ft) of average drawdown from 2012 to 2070. While the results of this study are meaningful for the water balance of the SHP, future work should consider the optimal spatial d istribution of reductions in water to minimize decreased well yields, and investigate the combination of intensive and extensive reductions in water use for yields as well as water table recovery. Further parameter estimation and sensitivity analysis is al so necessary for predictive modeling, particularly for hydraulic conductivity and specific yield. This will enable closer comparison between the baseline model scenario and the actual declines observed during the model period. Future work will include scen arios on the 50 - year timescale used for groundwater planning in Texas, using irrigation limits that are feasible for management targets. While decline appears inevitable as long as the SHP is used for irrigation, the pattern of decline, and methods for man aging drawdown, have large implications for the agricultural economy of the region. 112 REFERENCES 113 REFERENCES Anderson, M.P., W.W. Woessner, and R.J. Hunt, 2015. Applied Groundwater Modeling: Simulation of Flow and Advective Transport . Academic Press, San Diego, CA . 564p. Zellmer, 2018. The Future of Groundwater in California: Lessons in Sustainable Management from Across the West. Report from the Environmental Defense Fund and Daugherty Water for Food Global Institut e at the University of Nebraska, Lincoln, NE. 124p. Bossert, W.A., 1993. Overview of Kansas Groundwater Management Districts: Their Duties, Authorities and Expectations. Water Organizations in a Changing West , Natural Resources Law Center, University of Colorado School of Law. Bre dehoeft, J., 2002 . The water budget myth revisited: why hydrogeologists model. Ground Water 40 (4) : 340 - 345. Breña - Naranjo, J.A., A.D. Kendall, and D.W. Hyndman, 2014 . Improved methods for satellite - based groundwater storage estimates: A decade of monitoring the High Plains Aquifer from space and ground observations. Geophysical Research Letters 41 (17) : 6167 - 6173 ; doi:10.1002/2014GL0 61213. Brown, J.F. and M.S. Pervez, 2014. Merging remote sensing data and national agricultural statistics to model change in irrigated agriculture . Agr icultural Systems 127 : 28 - 40 ; doi:10.1016/j.agsy.2014.01.004. Butler, J., R. Stotler, D. Whittemore, a nd E. Reboulet, 2013. Interpretation of water - level changes in the High Plains aquifer in western Kansas. Groundwater 51(2): 180 - 190; doi:10.1111/j.1745 - 6584.2012.00988x. Cederstrand, J.R. and M.F. Becker, 1998a . Digital m ap of b ase of a quifer for High Pl ains a quifer in parts of Colorado, Kansas, Nebraska, New Mexico, Oklahoma, South Dakota, Texas, and Wyoming . USGS Open - File Report 98 - 393 , Reston, VA . Cederstrand, J. R. and M.F . Becker, 1998b. Digital map of hydraulic conductivity for the High Plains aqu ifer in parts of Colorado, Kansas, Nebraska, New Mexico, Oklahoma, South Dakota, Texas, and Wyoming . USGS Open - File Report 98 - 548, Reston, VA. Cede rstrand, J.R. and M.F. Becker, 1999. Digital map of predevelopment water levels for the High Plains Aquifer in parts of Colorado, Kansas, Nebraska, New Mexico, Oklahoma, 114 South Dakota, Texas, and Wyoming . USGS Open - File Report 99 - 264 , Reston, VA . Chen, J.M., X. Chen, W. Ju, and X. Geng , 2005. Distributed hydrological model for mapping evapotranspiration using re mote sensing inputs. Journal of Hydrology 305(1) : 15 - 39. Chen, X., Y. Yin, J.W. Goeke, and R.F. Diffendal, 2005. Vertical movement of water in a High Plains aquifer induced by a pumping well. Environmental Geology 47 (7) : 931 - 941. Cosgrove, B.A., D. Lohmann, K.E. Mitchell, P.R. Houser, E.F. Wood, J.C. Schaake, A. Robock, C. Marshall, J. Sheffield, Q. Duan, and L. Luo, 2003. Real time and retrospective forcing in the North American Land Data Assimilation System (NLDAS) project. Journal of Geophysical R esearch: Atmospheres 108(D22) ; doi:10.1029/2002JD003118 . Deeds, N.E. and M . Jigmond, 2015. Numerical Model Report for the High Plains Aquifer System Groundwater Availability Model . Prepared for the Texas Water Development Board, Austin, TX. Deeds, N. E., J . J. Harding, T . L. Jones, A. Singh, S. Hamlin, and R . C. Reedy (eds.), 2015. Final Conceptual Model Report for the High Plains Aquifer System Groundwater Availability Model . Prepared for the Texas Water Development Board, Austin, TX. Deines, J.M., A.D. Kendall, and D.W. Hyndman, 2017. Annual i rrigation d ynamics in the US Northern High Plains d erived from Landsat s atellite d ata. Geophysical Research Letters 44(18): 9350 - 9360; doi:10.1002/2017GL074071 . Dennehy, K. F. , D.W. Litke, and P.B. McMahon, 2002. T he High Plains Aquifer, USA: groundwater development and sustainability. Special Publications of the Geological Society of London 193 (1) : 99 - 119 ; doi:10.1144/GSL.SP.2002.193.01.09. Dlubac, K., R. Knight, Y. Song, N. Bachman, B. Gr au, J. Cannia, and J. Williams, 2013. Use of NMR logging to obtain estimates of hydraulic conductivity in the High Plains aquifer, Nebraska, USA. Water Resources Research 49 (4) : 1871 - 1886 ; doi:10.1002/wrcr.20151. Egbert, S.L., R.Y. Lee, K.P. Price, R. Boy ce, and M.D. Nellis, 1998. Mapping conservation reserve program (CRP) grasslands using multi seasonal thematic mapper imagery. Geocarto International 13(4) : 17 - 24. Fenichel, E.P., J.K. Abbott, J. Bayham, W. Boone, E.M.K. Haacker, and L. Pfeiffer, 2016. M easuring the value of groundwater and other forms of natural capital. 115 Proceedings of the National Academy of Sciences 113(9) : 2382 - 2387. Foster, T., N. and A.P. Butler, 2014. Modeling irrigation behavior in groundwater systems. Water R esources R esearch 50(8) : 6370 - 6389 ; doi:10.1002/2014WR015620 . Foster, T., N. Brozovi , and A. P. Butler, 2017. Effects of initial aquifer conditions on economic benefits from groundwater conservation. Water Resources Research 53(1) : 744 - 762; doi:10.1002/2016WR019365 . Freeze, R.A. and J.A. Cherry, 1979. Groundwater . Prentice Hall, Englewood Cliffs, NJ. 604 pp. Gaum, C. H., 1953. High Plains, or Llano Estacado, Texas - New Mexico. In The physical and economic foundation of natural resources, IV, Subsurface facilities o f water management and patterns of supply Type area studies , Interior and Insular Affairs Committee, House of Representatives, United States Congress, pp. 94 - 104. Gesch, D., M. Oimoen, S. Greenlee, C. Nelson, M. Steuck, and D. Tyler , 2002. The National El evation Dataset. Photogrammetric Engineering and Remote Sensing 68 (1) : 5 - 11. Gesch, D.B., 2007. The National Elevation Dataset. In Digital Elevation Model American Society for Photogrammetry and Remote Sensing, Bethesda, MD. Golden, B. and B. Guerrero, 2017. The economics of Local Enhanced Management Areas in southwest Kansas. Journal of Contemporary Water Research & Education 162(1) : 100 - 111; doi :10.1111/j.1936 - 704X.2017.03262 .x . Green, D.E., 1981. Land of the Underground Rain: Irrigation on the Texas High Plains, 1910 - 1970 . University of Texas Press , Austin, TX . 295 pp. Gutentag, E., F. Heimes, N. K rothe, R. Luckey, and J. Weeks, 1984. Geohydrology of the High Plains aquifer in parts of Colorado, Kansas, Nebraska, New Mexico, Oklahoma, South Dakota, Texas, and Wyoming: High Plains RASA . USGS P rofessional P aper 1400 - B. Haacker, E.M., A.D. Kendall, a nd D.W. Hyndman, 2 01 6 . Water level declines in the High Plains Aquifer: p redev elopment to resource senescence. Groundwater 54(2) : 231 - 242; doi:10.1111/gwat.12350 . Hansen, Z.K. and G.D. Libecap , 2004. Small farms, externalities, and the Dust Bowl of 116 the 1930s. Journal of Political Economy 112(3) : 665 - 694. Harbaugh, A.W., E.R. Banta, M.C. Hill, and M.G. McDonald, 2000. MODFLOW - 2000, the U. S. Geological Survey Modular Ground - Water Model User Guide to Modularization Concepts and the Ground - Water Flow Process . USGS Open - F ile Report 00 - 92, Reston, VA . Hecox, G.R., P.A. Macfarlane, a nd B.B. Wilson, 2 002. Calculation of y ield for High Plains w ells: r elationship between saturated thickness and well yield . Kansas Geological Survey Open - File Report 2002 - 25C, Lawrence, KS . Hellerstein, D., 2006. USDA Land Retirement Programs: Status, Hi story, and Trends. In Agricultural Resources and Environmental Indicators , Economic Information Bulletin No. 16 , Wiebe , K. and N. Gollehon , eds., USDA Economic Research Service, Washington, D.C . Hendricks, N.P. and J.M. Peterson, 2012. Fixed effects estim ation of the intensive and extensive margins of irrigation water demand. Journal of Agricultural and Resource Economics 37 (1) :1 19. High Plains UWCD, Underground Water Conservation District No. 1, 1998. The High Plains Underground Water Conservation Distr ict No. 1 Management Plan, 1998 - 2008 . http://www.twdb.texas.gov/groundwater/docs/GCD/hpuwcd1/hpuwcd1_m gmt_plan1998.pdf High Plains UWCD, Underground Water Conservation District No. 1, 2014. Water District Management Plan 2014 - 2024 . http://static.squaresp ace.com/static/53286fe5e4b0bbf6a4535d75/t/544fde09e4b 07b58025e6e3a/1414520329694/ManagementPlan.pdf. Hoekstra, A.Y. and T.O. footprint. Science 344(6188) : 1114 - 1117. Hornbeck, R. and P. Keskin, 2 014. The evolving impact of the Ogallala Aquifer: a gricultural adaptation to groundwater and climate. American Economic Journal: Applied Economics 6(1) : 1 90 - 219. Housto n, N.A., S.L. Gonzales - Bradford, and A.T. Flynn, 2011. Locations of wells and well log records used to develop the base of the Northern High Plains Aquifer underlying parts of Colorado, Kansas, Nebraska, South Dakota, and Wyoming . USGS Data Series 777, Reston, VA. Hyndman, D.W., A.D. Kendall, a nd N.R. Welty, 2007. Evaluating temporal and sp atial 117 variations in recharge and streamflow using the Integrated Landscape Hydrology Model (ILHM). In Subsurface H ydrology: D ata I ntegration for P roperties and P rocesses , D.W. Hyndman, F.D. Day - Lewis, K. Singha, eds., Geophysical Monograph Series vol. 171, American Geophysical Union, Washington, D.C . Isaaks, E.H. and R.M. Srivastava, 1989. An Introduction to Applied Geostatistics . Oxford University Press, UK. 561 pp. Jackson, R.B., J. Canadell , J.R. Ehleringer, H.A. Mooney, O.E. Sala, and E.D. Schulze, 1 996. A global analysis of root distributions for terrestrial biomes. Oecologia 108(3) : 389 - 411. Joshi, S .R. , 2005. Comparison of groundwater rights in the United States: l essons for Texas . Thesis submitted to Texas Tech University. Kansas Department of Agriculture, 2013. Order of Designation Approving the Sheridan 6 Local Enhanced Management Area within Groundwater Management District No. 4, April 17, 2013 . 12 WATER 8366. Kansas State Senate, 2012. Establishing local enhanced manage ment areas . Senate Bill 310. http://www.kslegislature.org/li_2012/b2011_12/measures/documents/sb310_ enrolled.pdf. Kendall, A.D., 2009. Predicting the impacts of land use and climate on regional - scale hydrologic fluxes . Thesis submitted to Michigan State U niversity. Kirchhoff, C.J. and L. Dilling, 2016. The role of US states in facilitating effective water governance under stress and change. Water Resources Research 52(4 ): 2951 - 2964; doi: 10.1002/2015WR018431 . Li, H. and J. Zhao, 2016. Rebound effects of new irrigation technologies: The role of water rights . Working paper, Michigan State University. Luckey, R. R. and M. Becker , 1999. Hydrogeology, w ater u se, and s imulation of f low in the High Plains a quifer in northwestern Oklahoma, southeastern Colorado, southwestern Kansas, northeastern New Mexico, and northwestern Texas . USGS Water - Resources Investigations Report 99 - 4104 , Reston, VA . Luckey, R.R. , E.D. Gutentag, and J.B. Weeks, 1981 . Water - level and saturated thickness changes, predevelopment to 1980, in the High Plains Aquifer in parts of Colorado, Kansas, Nebraska, New Mexico, Oklahoma, South Dakota, Texas, and Wyoming . USGS Hydrologic Investigations Atlas HA - 652 , Reston, VA . 118 Luo, W. and D.T. Pederson, 2012. Hydraulic conductivity of the High Plains A quifer re evaluated using surface drainage patterns. Geophysical Research Letters 39(2) : L02402 ; doi:10.1029/2011GL050200 . Luscz, E.C., A.D. Kendall, and D.W. Hyndman, 2017. A spatially explicit statistical model to quantify nutrient sources, pathways, an d delivery at the regional scale. Biogeochemistry 133(1) : 37 - 57 ; doi:10.1007/s10533 - 017 - 0305 - 1 . Mace, R.E., R. Petrossian, R. Bradley, W.F. Mullican, and L.A.N.C.E. Christian, 2008. A streetcar named desired future conditions: t he new groundwater availability for Texas (Revised). In The Changing Face Of Water Rights In Texas , 2008 C onference P roceedings . Mace, R.E., C. Ridgeway, and J.M. Sharp Jr, 2004. Groundwater is n o l onger s ecret and o ccult a h istorical and h ydrogeologic a nalysis of the East case. In 100 Years of Rule of Capture: From East to Groundwater Management , Mullican, W.F. III and S. Schwartz, S., eds . , Texas Water Development Board Report 361, Austin, TX . Mahmood, R. and K.G. Hubbard, 2002. Anthropogenic land - use change in the North American tall grass - short grass transition and modification of near - surface hydrologic cycle. Climate Research 21(1 ): 83 - 90. McGuire, V.L., M.R. Johnson, R.L. Schieffer, J.S. Stanton, S.K. Sebree, and I.M. Verstraeten, 2003. Water in s torage and approaches to ground - water management, High Plains aquifer, 2000 . USGS Open - File Report 2003 - 1243 , Reston, VA . McGuire, V.L., 2004. Water - Level c hanges in the High Plains a quifer, p redevelopment to 2002, 1980 to 2002, and 2001 to 2002 . USGS Fac t Sheet 2004 - 3026 , Reston, VA . McGuire, V.L., 2009. Water - level changes in the High Plains aquifer, predevelopment to 2007, 2005 - 06, and 2006 - 07 . USGS Scientific Investigations Report 2009 - 5019 , Reston, VA . McGuire, V.L., 2012 . Water - l evel and s torage c hanges in the High Plains a quifer, p redevelopment to 2011 and 2009 - 11 . USGS Scientific Investigations Report 2012 - 5291 , Reston, VA . McGuire, V. L. , K.D. Lund, and B.K. Densmore, 2012. Specific yield, High Plains aquifer. USGS Scientific Investigations Re port 2012 - 5177 , Reston, VA . Megdal, S.B., A.K. Gerlak, R.G. Varady, a nd L.Y. Huang, 2015. Groundwater governance in the United States: c ommon priorities and challenges. Groundwater 53(5): 677 - 684. 119 Muggeo, V.M. R , 2008. Segmented: an R package to fit regre ssion models with broken - line relationships. R news 8(1) : 20 - 25. Muggeo, V.M.R. , 2016 . Testing with a nuisance parameter present only under the alternative: a score - based approach with application to segmented modelling. Journal of Statist ical Computation and Simulation 86 (15): 3059 3067 ; doi: 10.1080/00949655.2016.1149855 . Myneni, R., Y. Knyazikhin, and T. Park, 2015. MOD15A2H MODIS Leaf Area Index/FPAR 8 - Day L4 Global 500m SIN Grid V006 . NASA EOSDIS Land Processes DAAC , Washington, D.C . doi:10.5067/MODIS /MOD15A2H.006 NASS, National Agricultural Statistics Service, 2017 a . Quick stats . https://quickstats.nass.usda.gov/. Accessed 4/ 2017. NASS, National Agricultural Statistics Service, 2017 b . Conservation Reserve Program Statistics: CRP Enrollment and Rental Payments by County, 1986 - 2016 (xls). https://www.fsa.usda.gov/programs - and - services/conservation - programs/reports - and - statistics/conservation - reserve - program - statistics/index. A ccessed 5/2017. Nellis, M.D., K. Price, S. Egbert, and J. Wu, 1996. Nat ural resource capability of CRP lands as grasslands in southwest Kansas: a remote sensing and GIS perspective. Geocarto International 11(3): 23 - 28. NOAA, National Oceanic and Atmospheric Administration, 2017. Historical Palmer Drought Indices . https://www .ncdc.noaa.gov/temp - and - precip/drought/historical - palmers/. Accessed 4/2017. Osterkamp, W.R. and W.W. Wood, 1987. Playa - lake basins on the Southern High Plains of Texas and New Mexico: Part I. Hydrologic, geomorphic, and geologic evidence for their develo pment. Geological Society of America Bulletin 99(2) : 215 - 223 Palmer, W.C., 1965. Me teorological drought (Vol. 30). Weather Bureau, United States Department of Commerce, Washington, DC. Perrone, D . and S . Jasechko, 2017. Dry groundwater wells in the Unite d States. Environmental Research Letters 12 ( 104002 ); doi: 10.1088/1748 - 9326/aa8ac0 . PRISM Climate Group, 2004. Annual precipitation reanalysis data product . Oregon State University, http://prism.oregonstate.edu, created 4 Feb 2004, accessed 4/ 2017. 120 Qi, S. L., 2010. Digital map of the aquifer boundary of the High Plains aquifer in parts of Colorado, Kansas, Nebraska, New Mexico, Oklahoma, South Dakota, Texas, and Wyoming . USGS Data Series 543 , Reston, VA . Rao, M.N. and Z. Yang, 2010. Groundwater impacts due to conservation reserve program in Texas County, Oklahoma. Applied Geography 30(3) : 317 - 328. Reeves Jr, C.C., 1970. Origin, classification, and geologic history of caliche on the southern High Plains, Texas and eastern New Mexico. The Journal of Geology 78(3) : 352 - 362. Reisner, M., 1986. Cadillac Desert . Penguin Books, New York, New York. 672p. Rosenberg, N.J., D.J. Epstein, D. Wang, L. Vail, R. Srinivasan, and J.G. Arnold, 1999. Possible impacts of global warming on t he hydrology of the Ogallala aquifer region. Climatic change 42(4) : 677 - 692. Sahoo, S., T.A. Russo, J. Elliott, and I. Foster, 2017. Machine learning algorithms for modeling groundwater level changes in agricultural regions of the US. Water Resources Rese arch 53(5): 3878 - 3895; doi:10.1002/2016WR019933. Scanlon, B.R., C.C. Faunt, L. Longuevergne, R.C. Reedy, W.M. Alley, V.L. McGuire, and P.B. McMahon, 2012. Groundwater depletion and sustainability of irrigation in the US High Plains and Central Valley. Pro ceedings of the National Academy of Sciences 109 (24) : 9320 - 9325 ; doi:10.1073/pnas.1200311109. Scanlon, B.R. and R.S. Goldsmith, 1997. Field study of spatial variability in unsaturated flow beneath and adjacent to playas. Water Resources Research 33(10 ): 2 239 - 2252. Scanlon, B.R., I. Jolly , M. Sophocleous, and L. Zhang, 2007. Global impacts of conversions from natural to agricultural ecosystems on water resources: q uantity versus quality. Water Resources Research 43(3); doi:10.1029/2006WR005486. Scanlon, B .R., K.E. Keese, A.L. Flint, L.E. Flint, C.B. Gaye, W.M. Edmunds, and I. Simmers, 2006. Global synthesis of groundwater recharge in semiarid and arid regions. Hydrological processes 20(15): 3335 - 3370. Scanlon, B.R., R.C. Reedy, D.A. Stonestrom, D.E. Prudic, and K.F. Dennehy, 2005. Impact of land use and land cover change on groundwater recharge and quality in the southwestern US. Global Change Biology 11(10) : 1577 - 1593. Shao, Q. and Campbell, N.A., 2002. Applications: Modelling trends in groundwater levels by segmented regression with constraints. Australian & New Zealand Journal 121 of Statistics 44(2) : 129 - 141. Smidt, S.J., E.M. Haacker, A.D. Kendall, J.M. Deines, L. Pei, K.A. Cotterman, H. Li, X. Liu, B. Basso, and D.W. Hyndman, 2016. Complex water ma nagement in modern agriculture: t rends in the water - energy - food nexus over the High Plains Aquifer. Science of The Total Environment 566: 988 - 1001; doi: 10.1016/j.scitotenv.2016.05.127 . Stanton, J., S. Qi, D. Ryter, S. Falk, N. Houston, S. Peterson, S. Wes tenbroek, and S. Christenson , 2011. Selected a pproaches to e stimate w ater - b udget c omponents of the High Plains, 1940 through 1949 and 2000 through 2009 . USGS Scientific Investigations Report 5183 , Reston, VA . Steward, D.R. and A.J. Allen, 2016. Peak groundwater depletion in the High Plains Aquifer, projections from 1930 to 2110. Agricultural Water Management 170 : 36 - 48. Soil Survey Staff, 2014. Soil Survey Geographic (SSURGO) Database. USDA Natural Resources Conservation Service, Washington, D.C . Terrell, B.L., P.N. Johnson, and E. Segarra , 2002. Ogallala aquifer depletion: economic impact on the Texas high plains. Water Policy 4 (1) : 33 - 46. Theis, C.V., 1937. Amount of ground water recharge in the southern High Plains. Eos 18(2) : 564 - 568. Torell, L.A., J.D. Libbin, and M.D. Miller , 1990. The market value of water in the Ogallala Aquifer. Land Economics 66 (2) : 163 - 175. Torres, J . M., D . W. Litke, and S . L. Qi, 1999. Bedrock formations underlying the High Plains Aquifer . USGS Open - F ile Report 99 - 214 , Reston, VA . Turner , C . G . , K. McAfee , S. Pandey, and A. Sunley , 2011. Irrigation metering and water use estimates: a comparative analysis, 1999 2007 . Texas Water Development Board Report 378 , Austin, TX . US Congress 1985. Food Security Act of 1985, P.L. 99 - 198, 99 Stat. 1504, Dec. 23, 1985, As Amended Through P.L. 112 - 240, Effective Jan. 2, 2013. 111 pp. US Congress 2008. Food, Energy, and Security Act of 2008, H.R. 2419. 628 pp. USDA, United States Department of Agriculture, 2016. C ro pland Data Layer, USDA National Agricultural Statistics Service, Washington, D.C. A ccessed 1/2017. USDA , United States Department of Agriculture Natural Resources Conservation 122 Service, USGS, United States Geological Survey, and EPA, Environmental Protecti on Agency, 2014. Watershed Boundary Dataset for High Plains States . Accessed 1/2017. USGS, United States Geological Survey, 2006. High Plains Water - Level Monitoring Study: Predevelopment. http://ne.water.usgs.gov/ogw/hpwlms/data.html. Accessed 11/2013. USGS, United States Geological Survey , in cooperation with U.S. Environmental Protection Agency, USDA Forest Service, and other Federal, State, and Local partners, 2012. National Hydrography Dataset , nhd.usgs.gov. Accessed 5/2012. USGS, United States Geol ogical Survey, 2016. The National Map, 3DEP products and services . The National Map , 3D Elevation Program , http://nationalmap.gov/3DEP/3dep_prodserv.html . Accessed 1/2017. US GS, United States Geological Survey, 2017. Gro undwater levels for the Nation . nwi s.waterdata.usgs.gov/usa/nwis/gwlevels. Accessed 1/2017. Voldseth, R.A., W.C. Johnson, T. Gilmanov, G.R. Guntenspergen, G.R. , and B.V. Millett, 2007. Model estimation of land use effects on water levels of northern prairie wetlands. Ecological Application s 17(2) : 527 - 540. White, S.E. and D.E. Kromm, 1995. Local groundwater management effectiveness in the Colorado and Kansas Ogallala region. Nat ural Resources J ournal 35 : 275 - 307 . White, W.N., W.L. Broadhurst, and J.W. Lang, 1946. Ground water in the High Plains of Texas. USGS Water - Supply Paper 889 - F, Reston, VA. Whittemore, D.O., J .J. Butler Jr, and B.B. Wilson, 201 4 . Assessing the major drivers of water - level declines: new insights into the future of heavily stressed aquifers. Hydrological Sciences Jour nal 61(1): 134 - 145; doi:10.1080/02626667.2014.959958 . Wilson, B., 2007. Ground - w ater l evels in Kansas: a b riefing to the Kansas Legislature House Agriculture and Natural Resources Committee . Kansas Geological Survey Open - file Report 2007 - 1 , Lawrence, KS . Wood, W.W. and W.R. Osterkamp, 1987. Playa - lake basins on the Southern High Plains of Texas and New Mexico: p art II. A hydrologic model and mass - balance arguments for their development. Geological Society of America Bulletin 99(2 ): 2 24 - 230. Wood, W.W. , K.A. Rainwater, and D.B. Thompson, 1997. Quantifying m acropore 123 r echarge: e xamples from a s emi a rid a rea. Ground Water 35(6) : 1097 - 1106. Wood, W.W. and W.E. Sanford, 1995. Chemical and isotopic methods for quantifying ground water recharge in a regional, semiarid environment. Ground Water 33(3) : 458 - 468. Young, A.R., M.E. Burbach, J.T. Korus, and L.M. Howard , 2012 . Nebraska Statewide Groundwater - Level Monitoring Report. Nebraska Water Survey Paper Number 80 . Institute of Agriculture and Natural Resources , University of Nebraska - Lincoln , Lincoln, NE .