Assessing the Impact of Artiﬁcial Recharge Ponds on Hydrological Fluxes in an Irrigated Stream–Aquifer System

: Artiﬁcial recharge ponds have been used increasingly in recent years to store water in underlying aquifers and modify baseline groundwater gradients or alter natural hydrologic ﬂuxes and state variables in an aquifer system. The number of constructed ponds, their geographic spacing, and the volume of water diverted to each pond can have a signiﬁcant impact on baseline system hydrologic ﬂuxes and state variables such as groundwater head, with the latter sometimes rising to cause waterlogging in cultivated areas. This study seeks to quantify the impact of recharge ponds on groundwater state variables (head, saturated thickness) and associated ﬂuxes within an irrigated stream-aquifer system. We use a numerical modeling approach to assess the impact of a set of 40 recharge ponds in a 246 km 2 region of the South Platte River Basin, Colorado on localized groundwater head, regional groundwater ﬂow patterns, and groundwater interactions with the South Platte River. We then use this information to determine the overall inﬂuence of recharge ponds on the hydrologic system. A linked agroecosystem–groundwater (DayCent-MODFLOW) modeling system is used to simulate irrigation, crop evapotranspiration, deep percolation to the water table, groundwater pumping, seepage from irrigation canals, seepage from recharge ponds, groundwater ﬂow, and groundwater–surface water interactions. The DayCent model simulates the plant–soil-water dynamics in the root zone and soil proﬁle, while MODFLOW simulates the water balance in the aquifer system. After calibration and testing, the model is used in scenario analysis to quantify the hydrologic impact of recharge ponds. Results indicate that recharge ponds can raise groundwater levels by approximately 2.5 m in localized areas, but only 15 cm when averaged over the entire study region. Ponds also increase the rate of total groundwater discharge to the South Platte River by approximately 3%, due to an increase in groundwater hydraulic gradient, which generally offsets stream depletion caused by groundwater pumping. These results can assist with groundwater resource management in the study region, and generally provide valuable information for the interplay between pumping wells and recharge ponds, and their composite effect on groundwater–surface water interactions. In addition, the developed linked DayCent-MODFLOW modeling system presented herein can be used in any region for which recharge rates should be calculated on a per-ﬁeld basis. Author Contributions: Conceptualization, C.D. and R.T.B.; methodology, C.D. and R.T.B.; software, C.D.; validation, C.D. and R.T.B.; formal analysis, C.D.; investigation, C.D. and R.T.B.; resources, C.D. and R.T.B.; data curation, C.D.; writing—original draft preparation, C.D.; writing—review and editing, R.T.B.; visualization, C.D. and R.T.B.; supervision, R.T.B.; project administration, R.T.B.;


Introduction
Artificial recharge ponds are often used in semi-arid and arid regions to store water in underlying aquifers [1][2][3] or alter the baseline groundwater gradients in an aquifer system. In terms of the former, recharge ponds can be part of an overall approach to manage aquifer recharge (MAR), with the benefits of not losing water to evaporation and requiring very little land for regional water storage. In terms of the latter, recharge ponds can be used to augment streamflow, as a means of offsetting stream depletion caused by groundwater pumping located in the alluvium of river corridors [4]. This is needed when groundwater is pumped out of priority in water right systems. For either purpose, water often is diverted Hydrology 2022, 9,91 2 of 17 from nearby streams or rivers and deposited into the recharge pond sites, with the water seeping through the pond bed into the underlying aquifer material.
There is a wide range of methods that can be used to recharge an aquifer, but five make up the vast majority of MAR techniques, including well, shaft and borehole recharge (57%); spreading methods (29%); in-channel modifications (7%); induced bank filtration (6%); and rainwater and runoff harvesting [2]. MODFLOW is the most commonly used groundwater flow model to simulate the impact of artificial recharge volumes and rates on local and regional groundwater levels and gradients [2].
Several studies have used a modeling approach to assess the impact of these different artificial recharge options on groundwater systems. Researchers [5] used a MODFLOW model to examine the impacts of artificial recharge on stream depletions in the Spokane Valley-Rathdrum Prairie aquifer of Idaho and Washington, USA, based on hypothetical recharge scenarios from potential injection wells and infiltration basins throughout the valley. Both MAR methods were found to be successful in increasing groundwater discharge to the Spokane River. Similarly, others [6] found that artificial recharge near the Upper San Pedro River (Arizona) can sustain baseflows and offset stream depletion caused by pumping. In another study [7], MODFLOW was used to forecast the effect of artificial recharge on regional water supplies for rural drinking water. The volume of recharge can also create a groundwater mound that prevents inflow of contaminated groundwater from irrigated fields. While these studies demonstrate the influence of artificial recharge on a groundwater system and associated fluxes, none quantify the influence of existing recharge ponds using historical recharge pond volumes and pumping rates.
To satisfy this gap in recharge pond understanding, the objective of this study is to quantify the influence of seepage from existing recharge ponds on groundwater system-response variables and fluxes in a highly managed irrigated stream-aquifer system. The study is performed for a 246 km 2 irrigated region located northeast of Denver, Colorado, wherein artificial recharge ponds are used in conjunction with water right augmentation plans. The impact of recharge ponds on the groundwater system is assessed using a linked agroecosystem-groundwater modeling system [8], with MODFLOW used to simulate groundwater head, groundwater storage, and the majority of groundwater sources and sinks, and DayCent used to estimate spatio-temporal recharge patterns to the water table. Specifically, the model, once calibrated and tested against field data, is used to explore the impact of recharge ponds on water table elevation and volumetric exchange rates between the aquifer and the South Platte River. The model is also used to determine if the recent implementation of recharge ponds is achieving the intended purpose, i.e., to offset the impact of pumping on streamflow depletion. The methods presented herein can be used to assess the impact of recharge ponds in other agricultural groundwater systems.

Study Area
The study area encompasses a 246 km 2 region located 64 km northeast of Denver, Colorado, within the South Platte River Basin ( Figure 1A) and specifically within the conductive Quaternary alluvium of the basin. The area includes the towns of Gilcrest and LaSalle, with a total population of about 3500. The study area is used mainly for agriculture. The main crops grown are corn, alfalfa, and grass pasture. The irrigation type is approximately 50% flood irrigation and 50% sprinkler irrigation, with the irrigation season from April through October and irrigation water obtained from four irrigation canals (diverting water from the South Platte River) or from the alluvial aquifer via groundwater pumping wells. Figure 1B shows the location of Gilcrest and LaSalle, the South Platte River, the four irrigation canals, and the location of 339 pumping wells (green dots) and 39 monitoring wells, with the latter used to monitor groundwater levels in the region. Many wells have pumping rates higher than 5450 m 3 /day (1000 gal/min). The topography includes a broad fluvial valley along the South Platte River. The land surface has an elevation of 1510 m in the south, lowering to an elevation of 1410 m in the northeast within the South Platte River channel. Highly permeable deposits are found in the central part of the aquifer. Groundwater flow is generally from south to north, following the topography, with groundwater discharging to the South Platte River [9]. Within the past 10-15 years, groundwater levels in the region have risen, leading to flooded basements, waterlogging of cultivated fields, and failure of septic systems.
Hydrology 2022, 9, x FOR PEER REVIEW 3 of includes a broad fluvial valley along the South Platte River. The land surface has an e vation of 1510 m in the south, lowering to an elevation of 1410 m in the northeast with the South Platte River channel. Highly permeable deposits are found in the central part the aquifer. Groundwater flow is generally from south to north, following the topograph with groundwater discharging to the South Platte River [9]. Within the past 10-15 yea groundwater levels in the region have risen, leading to flooded basements, waterloggi of cultivated fields, and failure of septic systems.  Within the appropriation doctrine of Colorado water law, the majority of the 339 pumping wells are junior in water right to water use from the South Platte River, and therefore any streamflow depletion (stream seepage or a decrease in groundwater discharge to the river) induced by pumping must be replaced by other sources of water. This has led to the construction of recharge ponds in the area (blue dots in Figure 1B), with the recharge from these ponds and the resulting rise in groundwater gradient and groundwater discharge to the South Platte River used to offset the pumping-induced streamflow depletion. There are 40 recharge ponds in the study area.
A MODFLOW model of the study area has previously been constructed and tested against measured groundwater levels from the network of 39 monitoring wells [10]. The model is based on water data (pumping, recharge, etc.) from the Colorado Decision Support System, and has been used to assess the causes of high water tables in the region, i.e., the relative influence of irrigation, canal seepage, and pumping on the regional water table [10]. Modeling results have indicated that recharge from surface water irrigation and canal seepage are the main controls on water table elevation, and therefore may be the cause of the high water table. However, the influence of constructed recharge ponds on water table elevation has not been properly assessed, since the majority of recharge ponds in the LaSalle/Gilcrest area have been constructed since 2013 and the simulation period of the original model application terminated in 2012.

Groundwater Numerical Modeling (MODFLOW)
This paper modifies the original MODFLOW model of the region [10] to extend the simulation through 2018 and to include recharge from the DayCent model. MODFLOW is a FORTRAN-written program that numerically solves the three-dimensional ground-water flow equation, using a finite-difference method [11]. The finite difference method is employed by discretizing the ground-water system spatially into a grid of cells. Each cell represents a volumetric portion of the aquifer and contains uniform hydraulic properties (hydraulic conductivity, specific yield, specific storage). A water balance equation is established for each grid cell, and the groundwater storage and associated head is solved at each time step according to groundwater inflows and outflows in three directions and the combined sources or sinks (e.g., pumping, recharge, groundwater-surface water exchange). The following equation shows the water balance for a cell in an unconfined aquifer [12]: where, in Equation (1), x, y, z are the three dimensions; h is groundwater head (L); K is hydraulic conductivity (L/T). S s is specific storage (1/T). ϕ is porosity taken equal to specific yield S y ; F s is the fraction of the cell thickness that is saturated; and f(F) is a function of F s that typically is set to 1.0 [11]. Inflows or outflows across each cell face is estimated using Darcy's law (Equation (2)): where Q is volumetric flow rate [L 3 /T], A is the cross-sectional area perpendicular to the flow [L 2 ], h 1 − h 2 is the head difference across the prism parallel to flow [L], and L is the length of the prism parallel to the flow path [L]. The MODFLOW flow model constructed for the study region has grid cells of 152.4 m × 152.4 m, resulting in 120 rows and 120 columns. The aquifer is discretized vertically by 10 layers, with the top elevation of the first layer extracted from the DEM (see Figure 2A) and the bottom of the lowest layer formed from the base of each borehole. Layer thickness ranges from 1 to 10 m, with the thinnest layers near the top of the model, for more accuracy in recharge exchange between DayCent and MODFLOW. Model active cells are shown in Figure 2B and represent the extent of the alluvial aquifer. 2A) and the bottom of the lowest layer formed from the base of each borehole. Lay ness ranges from 1 to 10 m, with the thinnest layers near the top of the model, accuracy in recharge exchange between DayCent and MODFLOW. Model active shown in Figure 2B and represent the extent of the alluvial aquifer.  In the study area, the South Platte River alluvial aquifer is a heterogeneous geologic unit composed of interbedded gravel, sand, silt and clay underlain by low-permeability bedrock shale. Aquifer thickness varies from 0 to more than 30 m, with most of the area having a thickness of 15-25 m. In the previous modeling study for this region [10], a 3D aquifer material map was developed by interpolating material data from 450 boreholes, using the Kriging method. The values of hydraulic conductivity and specific yield for each material type (clay, clay and silt, silt, silt and sand, sand, sand and gravel, gravel) were estimated via model calibration in the previous study [10], with values assigned to cells within each of the 10 layers of the MODFLOW grid. Figure 3 shows the estimated material map for layers near the top, middle, and bottom portions of the aquifer. The aquifer material is more conductive (sand and gravel) within and near the alluvium corridor of the South Platte River (shown along the north of the study area in black squares). Pockets of clay are present in more abundance in the middle layer, and the bottom of the aquifer has more sand and gravel areas. Figure 3D (lower right-hand map) shows the distribution of clay in the area. Clay layers have accumulated around the mid-south portion and east boundary of the model. Table 1 lists the final values of hydraulic conductivity K, specific yield S y , and specific storage S s for the various material types within the aquifer, as found through model calibration in the previous study [10].
Hydrology 2022, 9, 0 6 of 17 using the Kriging method. The values of hydraulic conductivity and specific yield for each material type (clay, clay and silt, silt, silt and sand, sand, sand and gravel, gravel) were estimated via model calibration in the previous study [10], with values assigned to cells within each of the 10 layers of the MODFLOW grid. Figure 3 shows the estimated material map for layers near the top, middle, and bottom portions of the aquifer. The aquifer material is more conductive (sand and gravel) within and near the alluvium corridor of the South Platte River (shown along the north of the study area in black squares). Pockets of clay are present in more abundance in the middle layer, and the bottom of the aquifer has more sand and gravel areas. Figure 3D (lower right-hand map) shows the distribution of clay in the area. Clay layers have accumulated around the mid-south portion and east boundary of the model. Table 1 lists the final values of hydraulic conductivity K, specific yield S y , and specific storage S s for the various material types within the aquifer, as found through model calibration in the previous study [10].  There are multiple groundwater sources and sinks, including irrigation recharge, precipitation recharge, canal seepage, recharge pond seepage, pumping, and aquifer-river exchange ( Figure 1B). There are four major irrigation canals in the study area: Union Ditch, Evans No 2, Farmers Independent, and Hewes Cook. The rate of canal seepage is calculated based on South Platte River canal diversion data from the Colorado Division of Water Resources (CDWR) and a specified portion (10%) of the canal water that seeps into the aquifer. The monthly volume of groundwater extraction from each pumping well and the monthly volume of water diverted into each of the 40 recharge ponds also are obtained from the CDWR database. We assume that recharge from the ponds is equal to the volume diverted to the ponds, minus evaporation from the pond surface. Daily evaporation depths (m/d) are computed using the Penman method [13], and then multiplied by the pond surface area (m 2 ) to obtain a daily volume of evaporation (m 3 ). Pumping rates are included in the MODFLOW model using the Well package; recharge from irrigation, canal seepage, and recharge pond seepage are included using the Recharge package. The calculation of recharge from irrigation is performed using DayCent, described in the next section. Monthly stress periods are used. Streamflow and stream-aquifer water exchange for the South Platte River (western boundary of the model domain) and Beebe Draw (eastern side of the domain) are simulated using the Stream Flow Routing (SFR) package, using a stream water budget and Manning's equation to compute the stream flow rate and depth. The width of the South Platte River was set to 45 m, and the width of Beebe Draw was set to 6 m. Time series of monthly rates (m 3 /day) for pumping, canal seepage, pond recharge, and groundwater recharge (deep percolation from DayCent) are plotted in Figure 4. Note the seasonal dependency of rates, as cultivation activities during the growing season (April-October) drive the sources and sinks of the aquifer.

Agroecosystem Model (DayCent)
The DayCent agroecosystem field-scale model is used to compute the deep percolation (i.e., recharge) to the water table. The DayCent model [14,15] is a medium complexity agroecosystem model. The major sub-models of DayCent are plant growth, soil water flow, soil organic carbon cycling, soil nutrient cycling, and greenhouse gas emission. Major inputs for the model are daily weather, soil physical properties, plant type, and management practices. DayCent has been widely used for carbon and nitrogen simulations in agroecosystems [16][17][18][19]. The crop growth/production sub-model has been used in simulations of agroecosystem dynamics globally [16,[20][21][22][23][24]. Recently, the DayCent modeling code has been improved in simulating crop canopy development, crop growth, and water use [15,25]. This newest version of DayCent is used in this study for linking with MOD-FLOW. DayCent is written in FORTRAN and C.
The DayCent modeling code includes the main water balance components for a soil profile: infiltration of precipitation and irrigation, surface runoff, deep percolation, evapotranspiration (ET; evaporation and transpiration), and capillary rise of groundwater: where ΔSi is the net change in soil water at the end of day i and i-1. In this equation, P, RO, and DP are precipitation, runoff, and deep percolation on day i, respectively. Inet is the net irrigation on day i. GW is the groundwater contribution if a shallow water table is present. ETc is the actual evapotranspiration on day i. All units are in cm day-1. The soil profile is defined by users and is usually less than three meters in depth. The input soil parameters include soil texture, bulk density, field capacity, wilting point, and saturated hydraulic conductivity. Precipitation and irrigation application are the source of surface water. A portion runs off the field, and the remainder infiltrates into the soil. As water flows in the root zone, a portion is removed by evapotranspiration (ET), and the remainder percolates to the water table. Reference ET is simulated in the model using either the standardized Penman-Monteith method [26] or the Hargreaves method [27], with the latter used when only air temperature is available. Crop coefficients are used in conjunction with reference ET

Agroecosystem Model (DayCent)
The DayCent agroecosystem field-scale model is used to compute the deep percolation (i.e., recharge) to the water table. The DayCent model [14,15] is a medium complexity agroecosystem model. The major sub-models of DayCent are plant growth, soil water flow, soil organic carbon cycling, soil nutrient cycling, and greenhouse gas emission. Major inputs for the model are daily weather, soil physical properties, plant type, and management practices. DayCent has been widely used for carbon and nitrogen simulations in agroecosystems [16][17][18][19]. The crop growth/production sub-model has been used in simulations of agroecosystem dynamics globally [16,[20][21][22][23][24]. Recently, the DayCent modeling code has been improved in simulating crop canopy development, crop growth, and water use [15,25]. This newest version of DayCent is used in this study for linking with MODFLOW. DayCent is written in FORTRAN and C.
The DayCent modeling code includes the main water balance components for a soil profile: infiltration of precipitation and irrigation, surface runoff, deep percolation, evapotranspiration (ET; evaporation and transpiration), and capillary rise of groundwater: where ∆S i is the net change in soil water at the end of day i and i − 1. In this equation, P, RO, and DP are precipitation, runoff, and deep percolation on day i, respectively. I net is the net irrigation on day i. GW is the groundwater contribution if a shallow water table is present. ET c is the actual evapotranspiration on day i. All units are in cm day-1. The soil profile is defined by users and is usually less than three meters in depth. The input soil parameters include soil texture, bulk density, field capacity, wilting point, and saturated hydraulic conductivity. Precipitation and irrigation application are the source of surface water. A portion runs off the field, and the remainder infiltrates into the soil. As water flows in the root zone, a portion is removed by evapotranspiration (ET), and the remainder percolates to the water table. Reference ET is simulated in the model using either the standardized Penman-Monteith method [26] or the Hargreaves method [27], with the latter used when only air temperature is available. Crop coefficients are used in conjunction with reference ET to estimate potential ET for each crop type. The ET is partitioned into potential evaporation and potential transpiration as a function of the green canopy coverage and residue coverage. The green canopy coverage (CC) is calculated from Beer's law [28,29]: where k (dimensionless) is the light extinction coefficient of the vegetation, and GLAI is green leaf area index (m m-1). The GLAI and CC approach was recently added to DayCent [15]. Water uptake by crop roots is limited by soil available water. In this study, DayCent is run for each of the active grid cells (11,372 total, see Figure 2B) in the top layer of the MODFLOW model, to facilitate linkage between DayCent deep percolation and MODFLOW's recharge package. For each DayCent model, the root zone is divided uniformly into 14 layers for a total thickness of 210 cm. The main inputs for DayCent include daily weather data, soil data, and crop type. The weather data (precipitation and temperature) are obtained from the Peckham station, part of CoAgMet (Colorado Agricultural Meteorological Network). The soil physical property data for the top soil containing the fraction of gravel, silt, and clay are collected from SSURGO (Soil Survey Geographic Database). The main crops are alfalfa, corn, grass pasture, wheat, sugar beets. Crop schedule files are created for each crop, including tillage and irrigation schedule information. Yearly crop rotation data are collected from satellite images from the National Agricultural Statistics Service (NASS) ( Figure 5) and mapped to the DayCent models based on the geological locations.
where k (dimensionless) is the light extinction coefficient of the vegetation, and GLAI is green leaf area index (m m-1). The GLAI and CC approach was recently added to DayCent [15]. Water uptake by crop roots is limited by soil available water. In this study, DayCent is run for each of the active grid cells (11,372 total, see Figure  2B) in the top layer of the MODFLOW model, to facilitate linkage between DayCent deep percolation and MODFLOW's recharge package. For each DayCent model, the root zone is divided uniformly into 14 layers for a total thickness of 210 cm. The main inputs for DayCent include daily weather data, soil data, and crop type. The weather data (precipitation and temperature) are obtained from the Peckham station, part of CoAgMet (Colorado Agricultural Meteorological Network). The soil physical property data for the top soil containing the fraction of gravel, silt, and clay are collected from SSURGO (Soil Survey Geographic Database). The main crops are alfalfa, corn, grass pasture, wheat, sugar beets. Crop schedule files are created for each crop, including tillage and irrigation schedule information. Yearly crop rotation data are collected from satellite images from the National Agricultural Statistics Service (NASS) ( Figure 5) and mapped to the DayCent models based on the geological locations.

Model Simulation, Calibration and Testing
The simulation period was 2012-2020. The MODFLOW model is loosely coupled with the DayCent model. The DayCent model runs the entire simulation first and generates deep percolation amounts in the output files. Then, the deep percolation is mapped directly to the corresponding MODFLOW grid cell, to be included in the Recharge package. The model is tested by comparing simulated and measured groundwater heads, with simulated values from grid cells that correspond to locations of monitoring wells (see Figure 1). In DayCent, one of the largest model factors impacting deep percolation is irrigation efficiency (i.e., the fraction of applied irrigation water that is used by the crops). Irrigation efficiency varies with irrigation type, weather, and irrigation management practices. In this study, irrigation efficiency is treated as an uncertain parameter that requires estimation. A range of irrigation efficiency values (25% to 100%) are tested in the DayCent-MODFLOW simulation, to determine which value provides the best fit between simulated

Model Simulation, Calibration and Testing
The simulation period was 2012-2020. The MODFLOW model is loosely coupled with the DayCent model. The DayCent model runs the entire simulation first and generates deep percolation amounts in the output files. Then, the deep percolation is mapped directly to the corresponding MODFLOW grid cell, to be included in the Recharge package. The model is tested by comparing simulated and measured groundwater heads, with simulated values from grid cells that correspond to locations of monitoring wells (see Figure 1). In DayCent, one of the largest model factors impacting deep percolation is irrigation efficiency (i.e., the fraction of applied irrigation water that is used by the crops). Irrigation efficiency varies with irrigation type, weather, and irrigation management practices. In this study, irrigation efficiency is treated as an uncertain parameter that requires estimation. A range of irrigation efficiency values (25% to 100%) are tested in the DayCent-MODFLOW simulation, to determine which value provides the best fit between simulated and observed groundwater head. Simulations were run successively until a best fit, in terms of root mean square error (RMSE), was found.

Estimating the Impact of Recharge Ponds on System Hydrologic Features
After the model is built, calibrated, and tested, the impact of the recharge ponds on groundwater state variables and associated fluxes is quantified. This is accomplished using two sets of simulations:

2.
Removing each recharge pond from the system one-at-a-time (40 simulations).
For both experiments, output groundwater heads are compared with original head values to quantify the extent and magnitude of recharge pond impact on localized and regional groundwater levels. In addition, exchange rates between the aquifer and the South Platte River are compared with original exchange rates, to determine the impact on groundwater return flows to the river.
A third set of simulations are run to determine the interplay between pumping wells and recharge ponds. The intent of constructing the recharge ponds is to offset stream depletions caused by groundwater wells pumping out of priority, to enable growers to still use groundwater for irrigation when they are junior in the Colorado water right system. The first simulation of this set runs the model without either groundwater pumping or recharge ponds, to determine pre-pumping patterns and rates of groundwater discharge to the South Platte River. A second simulation runs the model with groundwater pumping but no recharge ponds, to determine the amount of stream depletion caused by the use of the pumps. A third simulation runs the model with both groundwater pumping and recharge ponds (i.e., the baseline simulation), to determine if the pre-pumping rates of groundwater discharge are achieved, i.e., if the inclusion of the recharge ponds offsets the effect of the pumping wells on stream depletion. Whereas the first two simulation sets provide general insights into the impact of recharge ponds on hydrologic fluxes in the irrigated stream-aquifer system, the third simulation set determines if the recharge ponds are effective in the overall water management of the system. Figure 6A shows the RMSE of comparing the simulated groundwater head and observed head with different irrigation efficiencies. The lowest RMSE is 1.41 m with 50% irrigation efficiency, and therefore this irrigation efficiency % was used for all scenario simulations. This level of irrigation efficiency is indicative of a mixture of flood irrigation and sprinkler center-pivot irrigation systems. The ME and MAE are 0.6 m and 1.3 m. Figure 6B is the 1:1 plot of modeling head vs. observed head from more than 30,000 observations. Time series of the comparison of modeled and observed head are plotted as well for several monitoring well locations ( Figure 6C,D). The modeling results successfully catch the seasonal fluctuation in groundwater level. Figure 7 displays the total volumes (millions of m 3 ) associated with each groundwater source and sink during the 2012-2020 simulation period, for the baseline simulation. Among these, groundwater recharge, recharge, and canal seepage are the biggest components in water balance, corresponding to 749, 383, 345 million m 3 for the entire simulation period. Of the outputs (pumping and groundwater discharge), groundwater discharge to the South Platte River accounts for 94%. Of the inputs, groundwater recharge and canal seepage account for 46% and 38%, respectively. Note that a portion of groundwater pumping returns to the aquifer as irrigation recharge.        Figure 8 shows the decrease in groundwater head with the application of 0% (i.e., baseline simulation) to 80% of the current recharge pond seepage rates. The groundwater head can drop as low as 2.3 m in some areas and 0.15 m over the entire study area when no recharge ponds are applied. The average head decrease ranges from 0.14 to 0.67 m for the Gilcrest town area, a location where high groundwater levels have caused significant infrastructure damage. Figure 9 shows the change in hydrologic fluxes for each of these scenarios. The recharge pond seepage volume (red bars) increases from 0 percent to 100 percent, with the absence of recharge ponds (0% scenario), resulting in the largest decrease (2.8%) in groundwater discharge to the river being 2.8%. Eliminating recharge pond seepage decreases the regional groundwater gradient, decreasing the amount of groundwater flowing towards and discharging to the South Platte River. In general, decreasing recharge pond seepage can have a significant effect on local groundwater levels, lowering the water table sufficiently to decrease the impact on building foundations, residential basements, and wastewater treatment pond lining. However, the decrease in pond seepage results in an overall decrease in groundwater discharge to the South Platte River, which impacts the water rights of the river. Figure 10 shows a time series of percent decrease in groundwater discharge to the South Platte River, for each of the recharge pond seepage scenarios. The average decrease ranges from 0.37% to 3.6%. The decrease can be as high as 10%, around the summer months of 2016. The scenario with no recharge pond seepage (0% pond recharge) shows the largest decrease in groundwater discharge, corresponding to the result shown in Figure 9.

Recharge Pond Impact on Groundwater Head and Groundwater Discharge
(2.8%) in groundwater discharge to the river being 2.8%. Eliminating recharge pond seepage decreases the regional groundwater gradient, decreasing the amount of groundwater flowing towards and discharging to the South Platte River. In general, decreasing recharge pond seepage can have a significant effect on local groundwater levels, lowering the water table sufficiently to decrease the impact on building foundations, residential basements, and wastewater treatment pond lining. However, the decrease in pond seepage results in an overall decrease in groundwater discharge to the South Platte River, which impacts the water rights of the river. Figure 10 shows a time series of percent decrease in groundwater discharge to the South Platte River, for each of the recharge pond seepage scenarios. The average decrease ranges from 0.37% to 3.6%. The decrease can be as high as 10%, around the summer months of 2016. The scenario with no recharge pond seepage (0% pond recharge) shows the largest decrease in groundwater discharge, corresponding to the result shown in Figure 9.    . The interaction between the aquifer and the South Platte River corresponds to these different scenarios. X-axis is the percent of baseline recharge pond seepage. The percent values above the bars are percent decrease in groundwater discharge from the baseline simulation. The furthest right bars are for the baseline simulation. The impact of each individual recharge pond is also assessed. Figure 11A shows the average head decrease over the study region as each recharge pond is removed from the system. Pond #28 has by far the largest impact on groundwater head. Figure 11B shows the head drop map corresponding to a removal of Pond #28 from the groundwater system. The groundwater head decreases by 2 m in the localized area around the recharge pond. Figure 12 shows the impact of each individual recharge pond on groundwater discharge to the South Platte River. Pond #28 has the largest impact on groundwater discharge, with a daily average decrease of 6674 m 3 /day when this pond is removed from the system. The impact of each individual recharge pond is also assessed. Figure 11A shows the average head decrease over the study region as each recharge pond is removed from the system. Pond #28 has by far the largest impact on groundwater head. Figure 11B shows the head drop map corresponding to a removal of Pond #28 from the groundwater system. The groundwater head decreases by 2 m in the localized area around the recharge pond. Figure 12 shows the impact of each individual recharge pond on groundwater discharge to the South Platte River. Pond #28 has the largest impact on groundwater discharge, with a daily average decrease of 6674 m 3 /day when this pond is removed from the system.    Figure 13 summarizes the results from the third set of simulations, to quan interplay between pumping wells and recharge ponds. The diagram shows a cross of the model domain, from the high-elevation bluffs on the east to the South Plat on the west, with the groundwater sloping from west to east and groundwater di ing to the South Platte River. Arrows and values indicate the water balance of the under the baseline 2012-2020 simulation period, with gray arrows representing amounts of groundwater flux and black text indicating fluxes in millions of m 3 . Un scenario of no pumping and no recharge ponds, net groundwater discharge is 6 m 3 . When pumping is included, the discharge rate decreases by 2.9%, to 644 × Including recharge ponds with pumping results in a net discharge of 666 × 10 6 0.5% increase in discharge compared to the no pumping scenario. These results i that the recharge ponds indeed offset the stream depletion caused by pumping, an increase the discharge by 3 × 10 6 m 3 .  Figure 13 summarizes the results from the third set of simulations, to quantify the interplay between pumping wells and recharge ponds. The diagram shows a cross-section of the model domain, from the high-elevation bluffs on the east to the South Platte River on the west, with the groundwater sloping from west to east and groundwater discharging to the South Platte River. Arrows and values indicate the water balance of the region under the baseline 2012-2020 simulation period, with gray arrows representing relative amounts of groundwater flux and black text indicating fluxes in millions of m 3 . Under the scenario of no pumping and no recharge ponds, net groundwater discharge is 663 × 10 6 m 3 . When pumping is included, the discharge rate decreases by 2.9%, to 644 × 10 6 m 3 . Including recharge ponds with pumping results in a net discharge of 666 × 10 6 m 3 , or a 0.5% increase in discharge compared to the no pumping scenario. These results indicate that the recharge ponds indeed offset the stream depletion caused by pumping, and even increase the discharge by 3 × 10 6 m 3 .
of the model domain, from the high-elevation bluffs on the east to the South Platte River on the west, with the groundwater sloping from west to east and groundwater discharging to the South Platte River. Arrows and values indicate the water balance of the region under the baseline 2012-2020 simulation period, with gray arrows representing relative amounts of groundwater flux and black text indicating fluxes in millions of m 3 . Under the scenario of no pumping and no recharge ponds, net groundwater discharge is 663 × 10 6 m 3 . When pumping is included, the discharge rate decreases by 2.9%, to 644 × 10 6 m 3 . Including recharge ponds with pumping results in a net discharge of 666 × 10 6 m 3 , or a 0.5% increase in discharge compared to the no pumping scenario. These results indicate that the recharge ponds indeed offset the stream depletion caused by pumping, and even increase the discharge by 3 × 10 6 m 3 .

Summary and Conclusions
In this study, the impact of recharge ponds on the groundwater system of the Gilcrest/La Salle area in the South Platte River Basin is assessed. This is achieved by introducing a new loosely coupled model that links the MODFLOW groundwater flow model and the DayCent agronomic hydrologic model. Daily deep percolation simulated by field-scale DayCent models is mapped to the MODFLOW grid cells, for use in the Recharge package. In this model application, we assume that water exchange between the soil profile and the water table is uni-directional (downward). If exchange in both directions is desired, the reader is referred to a modeling study [8] that tightly couples DayCent and MODFLOW.
The model is used to quantify the influence of 40 artificial recharge ponds on groundwater head and groundwater-surface water interactions during the 2012-2020 historical period. Several sets of simulations are run to quantify the effect, collectively and individually, of recharge pond seepage on local and regional groundwater levels and groundwater discharge to the South Platte River. Results show that groundwater levels can be decreased by up to 2.3 m in some areas of the study region, which can have a significant effect on historical groundwater flooding. Damage to building foundations, residential basements, and wastewater treatment plant bottom liners, all of which have been impacted by high water tables in recent years, could be mitigated by a decrease in groundwater head of this magnitude.
However, this decrease in groundwater levels results in a decrease in groundwater discharge to the South Platte River by approximately 3%. As the intent of the recharge ponds is to increase groundwater discharge and thereby offset stream depletions caused by groundwater pumping, and simulation results do indicate that pond seepage offset but do not over-offset depletions during the 2012-2020 time period. Mitigating high water table issues in the region can be achieved only by (1) modifying fluxes of other sources and sinks of groundwater, or (2) modifying or relaxing the adjudication of water law, which dictates the need for offsetting pumping-induced stream depletion, in this region.