Next Article in Journal
A New Method to Estimate Reference Crop Evapotranspiration from Geostationary Satellite Imagery: Practical Considerations
Next Article in Special Issue
Pan-European Calculation of Hydrologic Stress Metrics in Rivers: A First Assessment with Potential Connections to Ecological Status
Previous Article in Journal
Spatial and Temporal Variation of Dissolved Heavy Metals in the Mun River, Northeast Thailand
Previous Article in Special Issue
Estimation of Water Budget Components of the Sakarya River Basin by Using the WEAP-PGM Model
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Modeling Riparian Restoration Impacts on the Hydrologic Cycle at the Babacomari Ranch, SE Arizona, USA

1
U.S. Geological Survey, Western Geographic Science Center, 520 N. Park Avenue, Ste. #102K, Tucson, AZ 85719, USA
2
U.S. Geological Survey, Arizona Water Science Center, 520 N. Park Avenue, Tucson, AZ 85719, USA
3
Lacher Hydrological Consulting, Tucson, AZ 85719, USA
4
Hydrology and Atmospheric Sciences, University of Arizona, Tucson, AZ 85719, USA
5
BIO5 Institute, University of Arizona, Tucson, AZ 85719, USA
*
Author to whom correspondence should be addressed.
Water 2019, 11(2), 381; https://doi.org/10.3390/w11020381
Submission received: 21 December 2018 / Revised: 12 February 2019 / Accepted: 18 February 2019 / Published: 22 February 2019
(This article belongs to the Special Issue Hydrologic Modelling for Water Resources and River Basin Management)

Abstract

:
This paper describes coupling field experiments with surface and groundwater modeling to investigate rangelands of SE Arizona, USA using erosion-control structures to augment shallow and deep aquifer recharge. We collected field data to describe the physical and hydrological properties before and after gabions (caged riprap) were installed in an ephemeral channel. The modular finite-difference flow model is applied to simulate the amount of increase needed to raise groundwater levels. We used the average increase in infiltration measured in the field and projected on site, assuming all infiltration becomes recharge, to estimate how many gabions would be needed to increase recharge in the larger watershed. A watershed model was then applied and calibrated with discharge and 3D terrain measurements, to simulate flow volumes. Findings were coupled to extrapolate simulations and quantify long-term impacts of riparian restoration. Projected scenarios demonstrate how erosion-control structures could impact all components of the annual water budget. Results support the potential of watershed-wide gabion installation to increase total aquifer recharge, with models portraying increased subsurface connectivity and accentuated lateral flow contributions.

1. Introduction

Given a growing concern regarding global water shortages, the potential to increase groundwater recharge offers promise for the sustainable management of growing populations [1,2]. Although the concept of managed aquifer recharge (MAR) is relatively new, it is being implemented worldwide, especially in aridlands where shortages are common, and the crises are more imminent. Artificial recharge of groundwater is achieved by facilitating the downward movement of surface water into the soil, first as infiltration and then as recharge when it reaches the water table [3].
Erosion-control structures (ECS; i.e., check dams, gabions, etc.) impede storm runoff and soil erosion, facilitating infiltration [4]. Various researchers have documented increased groundwater levels at check dams in Yemen, Saudi Arabia, Pakistan, Spain [5], and India [6], finding increases in recharge varying from 6% to 40% [2]. ECS cause water to be detained in pools, with volume dependent on size and structure. Maximum depths of 0.77 m were documented in ephemeral shallow pools formed behind ECS in SE Arizona [7]. This is similar to what has been described for beaver ponds, from 0.7 m to 2.0 m [8,9]. The formation of in-channel pools and ponds can increase storage and residence time, altering streamflow. Norman et al. [10,11,12,13] found decreases in flood events, but increases in volume of flow and increases in sediment retention behind ECS in SE Arizona and northern Sonora, Mexico. Vegetation is maintained at sites restored with, and downstream from, ECS—despite drought conditions in aridlands [14,15]. Nichols et al. [16] found increased soil moisture through the channel bank soil profiles at check dams vs. at control sites in SE Arizona.
We posit that synthetic ECS impact flows in a manner similar to beaver dams to support increased recharge rates. Dams created by beavers along the stream channel result in ponds that also raise the water table in the adjacent riparian zone [8,17]. White [9] studied the convective flow patterns beneath beaver dams, where underflow carries convected streamwater beneath the structures. Gurnell [18] describes a correlation between beaver ponds and lateral riparian groundwater levels. These structures increase lateral connectivity, forcing water sideways and creating diverse wetland environments [17]. Puttock et al. [19] hypothesized that beaver-constructed features increase water storage within the landscape, creating a stepped profile channel.
To consider the effects of increased recharge via ECS on the entire water budget, an understanding of response and transport from headwaters to the larger watershed is necessary at different scales. Arnold et al. [20] used the Soil and Water Assessment Tool model (SWAT) to discern recharge, finding some recharged water re-emerges as lagged baseflow, while some is lost to evapotranspiration (ET). Sun and Cornish [21] used SWAT for recharge estimation in the headwaters of the Liverpool Plains in Australia, finding that recharge only occurred in wet years with variations primarily explained by climatic factors (rather than land-use changes). Results of their study demonstrate the need for catchment-scale modeling (not point-scale) to predict catchment-scale recharge and suggest that improvements to SWAT might be needed for dry-climate modeling. Kim [22] combined the outputs of SWAT with the Modular Three-Dimensional Finite-Difference Groundwater Flow Model (MODFLOW; [23]) to simulate groundwater recharge, allowing for the interaction between the stream network and the aquifer in the Musimcheon Basin of Korea.
In this paper, we explore ECS-induced recharge at an aridland ranch in SE Arizona using the quasi-distributed watershed model, SWAT, with the fully distributed groundwater model, MODFLOW (U.S. Geological Survey: Reston, VA, USA). In 2013, Borderlands Restoration, LLC (Limited Liability Company), installed a series of large gabion structures in the Vaughn Canyon sub-watershed, of the Upper Babocomari Watershed (ex. Figure 1) [24,25,26]. Our hypothesis is that hydrologic modeling can be used to combine watershed- and point-scale data to understand and predict the dynamic functions controlling infiltration and recharge at gabions and to construe a catchment-scale understanding of the hydrologic budget. We describe how data collected in the study area were used to build a conceptual model, validate predictions, and extrapolate findings into future scenarios.

Study Area

The Upper Babocomari sub-watershed is a tributary of the San Pedro River that originates near Cananea, Sonora, Mexico, flowing north into southeastern Arizona, United States. The San Pedro River is ~338 km long, ~58 km of which were designated as the San Pedro Riparian National Conservation Area (SPRNCA) by Congress in 1988 to protect and enhance the desert riparian ecosystem in Southeast Arizona (Figure 2). Water availability in the SPRNCA is dependent on flood events and groundwater discharge to the river (baseflow). This groundwater largely originates as precipitation in the Huachuca Mountains, west of the river [27] and generally flows from the mountain fronts toward the center of the valley [28]. Multiple regional groundwater models have been developed to evaluate groundwater conditions and model future water-use scenarios and predict how groundwater withdrawals will impact the Upper San Pedro regional aquifer [28,29,30,31,32,33,34]. Pumping in the Fort Huachuca, Sierra Vista, and Huachuca City areas has created a cone of depression in the groundwater table [31,32,35,36,37,38,39].
The Babocomari River is one of the main tributaries to the San Pedro River. It is classified as a predominantly ephemeral stream, but has perennial stretches supported by the regional groundwater system. The Babocomari River lies mostly within the boundaries of the Babacomari Ranch (spelling varies from river to ranch) just north of the Fort Huachuca Military Reservation. The Brophy family bought the Babacomari Ranch in 1935 and practice conservation ranching. Groundwater flow patterns in the headwaters area of the Babocomari River are poorly defined due to a lack of monitoring wells in the area, but the Arizona Department of Water Resources (ADWR) well registry depicts levels in the Babocomari Watershed between 1166m and 1816 m [24]. Pools of standing water imply the existence of shallow water tables along the Babocomari River, but a Fort Huachuca test well depicted average declines of 11.1 cm/yr [36], and the ADWR [40] measured declines of up to 1.2 m in wells near the far west end of the Babacomari Ranch over the 2001–2006 period.
Groundwater discharge to the Babocomari River was estimated at 2.5M m3 per year [36]. Sharma et al. [41] analyzed stream flow and groundwater data collected by the US Bureau of Land Management (BLM) on the San Pedro and Babocomari rivers to establish a more efficient monitoring program for the SPRNCA and reported measurements in the Babocomari ranging from no flow to 0.04 cubic meters per second (m3/s) from intermittent gaging between 1990 and 1995. Corell et al. [31] predicted that water levels in upper reaches of Babocomari could decline up to 4.6 m and Babocomari River flows decline from 0.03 to 0.01 m3/s by 2030. Lacher [32] also simulated up to 0.03 m3/s of pumping-induced declines in baseflow in the upper half of the Babocomari over the 21st century. This tributary contributes about 7.4M m3 of water from both groundwater and surface-water sources, per year to the San Pedro River (Arizona Department of Water Resources 2005). Evapotranspiration (ET) losses on the Babocomari River were estimated at ~0.74M m3 (10% flow) [31].
A large structure was installed in the 1940s as a federal Works Progress Administration (WPA) project for erosion control and as an irrigation reservoir, with concrete spillway, some cement dykes, and ~91 m earthen dam. There is a drainage pipe through the dam that is not regularly maintained and an irrigation pipe occasionally (i.e., once/year) used to water an agricultural field to the east. The dam protects the Babocomari Cienega, one of the last intact and largest southwestern wetlands left in the region [42]. Cienegas have water tables at or near the surface, highly unusual and sought-after in aridlands. Soils in the wetland are high in clay and organic matter and the plant community is a mixture of Sporobolus sp. (Sacaton) and grass-like plant species of wetlands (sedges and rushes) with several rare plants including one endangered species, a rare orchid, Canelo Hills Ladies Tresses (Spiranthes delitescens) [43]. In addition, a high density of phreatophyte vegetation inhabits the downstream edge including a cattail (Typha sp.) invasion which extends downstream at least 1 km from the lower pond (wetland).
Beaver were extirpated from the San Pedro River by 1894, but were reintroduced in 1999 by the BLM, who brought 15 beavers to the SPRNCA; these populations are expanding [44,45]. In 2016, we identified new beavers in the Babocomari River, with 2 beaver dams and 1 lodge; the larger beaver pond had a ~50 m-wide dam. As noted, beaver dams can affect the hydraulic flow, sediment loading, nutrient transport, and biological functions of a stream, as well as development or degradation of the stream channel itself [8,18,46,47].
The Upper Babocomari gaging station (# 09471380) has operated by the US Geological Survey (USGS) in cooperation with the BLM since 2000. It is located approximately 29 km upstream of the confluence with the San Pedro River and approximately 6.4 km east of the ranch headquarters (Figure 3). Lacher [25] attributes some of the perennial flow here to the groundwater levels close to the surface of the streambed (indicating a water table reliably intersecting the stream channel), as well as slow seepage from the pond at the ranch headquarters. Average discharge at the gage is ~ 0.06 m3/s (Table 1) with high flows occurring once or twice per year in response to rainfall during the summer monsoon.
Rainfall is spatially distributed in the watershed, ranging from 33 cm to > 75 cm annually [39], with bimodal temporal distribution. Summer months (July–October) receive ~ 71% of the rainfall (approximately 26% of which occurs in August) and 17% occurs in winter months (December–January).
The geology of bedrock has a potentially strong influence on surface- and groundwater flow and storage properties in the watershed. The mountains around the Babocomari River are composed of Mesozoic and Paleozoic sedimentary rocks including limestone, overlying older Proterozoic granites, and cut by thrust faults that leave east-west trending ridges and steep northwest strikes with vertical displacement [48,49]. The Babocomari River is aligned with one of these faults [43]. The impacts of these structures depend on their depths and permeabilities. Limestone can be especially permeable, and cavernous rock supporting perennial springs and streams has been documented in the Huachuca Mountains [50].
The combination of mountains with higher precipitation, cavernous-limestone conduits, subsurface discharge from the mountain block, and shallow bedrock is likely a significant factor in maintaining groundwater levels and streamflow in the upper Babocomari River area, even in the presence of the pumping center and cone of depression to the east. The Babocomari watershed overlies two structural basins, Sonoita in the west and San Pedro in the east. The Sonoita Basin is shared among three watersheds. These basins contain basin fill and alluvial deposits of varying degrees of consolidation and poorly constrained thickness and groundwater depth at many locations. Thus, groundwater flow direction and the subsurface flow divide are not well known.
The Babocomari watershed boundary to the west and northwest is a low ridge of alluvium (basin fill) in the Sonoita grasslands [43,49]. Upland soils are clay and loam and produce luxuriant plains grassland vegetation [51]. These soils are mostly underlain by conglomerate and fanglomerate bedrock [48,49]. Soils found in the riparian grasslands are formed in younger (Holocene) alluvium and include soil with sandy loam and silt loam textures with some clay (i.e., Pima Series with clayloam and silty clayloam textures) [43,52]. Soils found along rivers and in stream channels are sandier and gravellier in texture [36,43,52].
Three types of vegetation, largely determined by aspect and elevation, exist in the watershed: (i) the Mountains include oak woodland and chaparral, with Douglas Fir and Ponderosa Pine; (ii) the Basin floor is desert grassland; and (iii) the Riparian zone contains shrubs, grasses, and herbs, as well as trees (Fremont cottonwood) that thrive in areas with water tables within 1–2 m of land surface [36]. Robinett and Kennedy [43] established riparian monitoring stations in 2010 and revisited them annually at the Babocomari for 4 years, finding stable or improving vegetative trends, good ecological conditions in riparian plant communities, and positive geomorphic conditions.

2. Methods

In the upper Babocomari watershed, a series of gabions were installed in Vaughn Canyon, an ephemeral tributary, by Borderlands Restoration in 2014–2015. This paper describes how research conducted in the field is used in an iterative process with groundwater and watershed models to define a study area, extrapolate recharge estimates, quantify the volume of water that structures could detain, and estimate the larger influence on components of the water budget (Figure 4). Soil moisture and infiltration metrics were used as proxies for groundwater availability and recharge. Subsurface drainage response is explored in relationship to land management and recharge. We present the benefits of using a model to examine and predict hydrologic response.

2.1. Data Collection and Analysis

Several data types were collected between 2014 and 2016, to determine infiltration rates in Vaughn Canyon stream channels and estimate the permeability of the unsaturated zone between the land surface and the aquifer. Parameters derived from the data collection and analyses were used to further understanding of a conceptual model and validate input to hydrologic models.

2.1.1. Infiltration

Fandel et al. [53,54,55] tested low-cost methods for quantifying the total infiltration induced by gabion construction over the course of a single flow event. Time-lapse cameras, buried temperature loggers, and pressure transducers recorded data from June 2015 to February 2016 at Gabion #2 (Figure 5 and Figure 6). All instruments were synchronized to collect a data point every five minutes for the duration of the study period. Soil properties were measured and averaged at 5 different locations in Vaughn Canyon. At the beginning of the study, a borehole permeameter was used to determine in situ saturated hydraulic conductivity (K), and soil samples were taken to quantify dry bulk density and porosity at the structure [54].
Point infiltration fluxes were derived using established fluid and heat transport equations [56] and applied to temperature time series and soil properties, assuming a unit hydraulic gradient. These point fluxes were upscaled to determine the infiltration flow rate (Qinf) in the gabion’s area of influence, based on the surface area of ponded water at each timestep, as inferred from camera footage [54]. Integrating the flow rate over the event duration yielded the total infiltration volume upstream of Gabion #2 for that flow event. To quantify the gabion’s impact, this was compared to the estimated infiltration volume before gabion construction. To account for uncertainty associated with visual estimates of ponded surface area, a range of possible areas was calculated based on measured gabion dimensions superimposed on the photos (Figure 6), observed channel morphology, and measuring the channel floor from satellite imagery to calculate the total wetted surface area between gabions.

2.1.2. Erosion and Deposition

Leica MS60 (Leica Geosystems: Balgach, Switzerland) and GS14 GPS (Leica Geosystems: Balgach, Switzerland)) equipment was used to measure Gabion # 2 (Figure 7). The Leica MS60 is a terrestrial- Light Detection and Ranging (LiDAR) scanner that can interface with RTK GPS equipment that provides all scan data with accurate ground control. Four reference marks were driven into the banks and act as ground control points (GCP) for the project. These GCP were surveyed with RTK GPS equipment and corrected to the NAD 83 datum and are used to correct all scanned points to the survey datum and used for subsequent scans to measure land-surface change. The two point-cloud scans have an accuracy of 0.15 cm and a total of ~5.3 million points.
In CloudCompare, an open-source 3D point-cloud and mesh-processing software, we ran the Cloth Simulation Filter [57] on the two-point clouds to generate a rough bare-ground point cloud, generated a minimum height model using the “ground points”, and exported the datasets as rasters. The channel is well defined, but the sides of the arroyo are completely obscured by Sporobolus sp. (Sacaton) and other vegetation higher up that decreases the accuracy of a digital elevation model (DEM) away from the sandy channel. Using a shapefile digitized to include only the sediment channel and avoid areas along the banks with incomplete scanning (which included grass heights), we clipped the extent of the raster and subtracted the baseline dataset from the original to create a DEM of difference (DoD).
To derive the area and volume calculations from the rasters, we reclassified our DoD and simplified values to either positive or negative values and used cell counts to estimate area of deposition and erosion. The zonal statistics tool was used to calculate the mean depth of each zone (positive and negative zones) and net volume and sediment changes. This was converted to mass by multiplying the specific gravity of sediment (average 2 metric tons/m3).

2.2. Models

We examined how coupled field and modeling approaches can be used, both separately and together, to estimate impacts of restoration, finding in general the potential to validate each other and better inform our results.

2.2.1. MODFLOW

Pool and Dickinson [28] released a 5-layer MODFLOW model for the entire Upper San Pedro Basin (USPB), including the portion in Mexico, simulating conditions from 1902 to 2003. Using the Pool and Dickinson model, Lacher [32] updated pumping values across the basin to 2009-10 values, and projected changes in pumping based on US Census data and population estimates for the period 2010 to 2100. Lacher [25,26] used the updated model to evaluate various restoration sites for groundwater response to recharge and for the amount of recharge required to produce significant (>0.03 m3/s) response in baseflow in the Babocomari River within 16 years. The groundwater model used herein was a preexisting regional model of the Upper San Pedro Basin. Details of the model structure, calibration, and error are provided in the original modeling [28]. The 250 m × 250 m grid size of that 5-layer regional model, as well as the lack of shallow groundwater data in the study area, limit the model’s utility for predicting the impacts of small surface-water structures (e.g., check dams) on localized recharge. However, in simulating various hypothetical recharge rates near the study site, the model provided a mechanism for understanding: (1) the potential impacts of local hydrogeology on the ability of surface recharge to manifest as baseflows in the nearby Babocomari River; (2) the effect that recharge location may have on hydrologic responses to recharge; (3) the scale of long-term recharge required to achieve a groundwater driven increase in baseflow in the Babocomari River within a specific time frame.

2.2.2. SWAT

The SWAT is a code that simulates watershed processes and water balance using information on climate, topography, soil and land use [58,59,60]. SWAT divides soils into layers (depth) and the water balance equation is solved for each, depending on saturated conductivity and soil-water content of the layer. When rainfall occurs, runoff and infiltration volumes are estimated. Some water moves vertically through soil layers when water content exceeds capacity, eventually becoming recharge to shallow groundwater. Some of the groundwater contributions can be accessed in SWAT by capillary activity during dry periods [21]. Some water is divided into multiple layers for routing and some moves through the soil via “lateral flow” or unsaturated flow. The combination of surface runoff, lateral flow, and return flow comprises the estimated streamflow. SWAT calculates lateral flow with water redistribution in the soil profile (0–2 m) using the kinematic-storage model to predict storage and flow in each layer [61]. The kinematic-storage model accounts for variations in conductivity, slope and soil-water content and allows for flow upward to an adjacent layer or to the subsurface [60]. This model assumes that the vertical input rate to the saturated zone is a function of the volume of water stored in the unsaturated zone. It cannot simulate the water content gradients in the unsaturated zone.
Calibration is suggested to improve model performance and increase the reliability of flow predictions, especially to analyze the effect of land-use change impacts [62,63]. We modeled the Upper Babocomari Watershed (412 km2; Figure 3) using associated precipitation and discharge gages to calibrate the model for volume estimates and using documented local parameters.

(a) Watershed Delineation, Monitoring Points, and Reservoirs

Watershed characteristics were simulated for the Upper Babocomari, using the USGS gaging station as the outlet (Figure 3). A high-resolution DEM, derived from aerial LiDAR flown in 2004 (personal comm. John Hays, Santa Cruz County Flood Control District), includes portions of the headwaters of the Upper Babocomari Watershed. Data were processed to create a 2 m bare-earth model based on points and converted to ASCII format. The highest-resolution DEM available for the rest of the watershed (Cochise County) is 1/3 arc-second (10 m) from the National Elevation Dataset (NED). These were mosaicked and resampled to mimic the higher resolution (2 m) for the study area. We divided our merged DEM-derived slope into 3 classes: (1) 0–5%, (2) 5–25%, and (3) >25%.
In SWAT, surface-water bodies (ponds and wetlands) on the main channel can be modeled as "reservoirs" and have a substantial influence on outflow and discharge. The WPA ponds and the 2 beaver dam ponds into separate “reservoirs” (due to their proximity and smaller size), averaging 1 ha (minimum entry) when the reservoir is filled to the emergency spillway. Almendinger et al. [64] found that seepage from reservoirs in SWAT is not included as a component of groundwater recharge, which does not accurately portray reality, but can be accommodated by disallowing any seepage, by setting the K of the pond bottom to zero [64]. We added 3 model-reporting outlet locations at: (a) Gabion #2 in Vaughn Canyon, (b) the WPA dam, and c) the USGS gaging station (Figure 3).
Soils and land-use data were collected and assembled in a geodatabase. Soils were mapped by the US Department of Agriculture (USDA) Natural Resources Conservation Service (NRCS) at a scale of 1:12,000, in the Soil Survey Geographic Database (SSURGO) data. The USGS National Land-Cover Database (NLCD2011) is a 16-class land-cover classification based on the Landsat Enhanced Thematic Mapper+ (ETM+) from 2011 at a spatial resolution of 30 m [65]. All geospatial data were clipped to the study area and converted to a common Universal Transverse Mercator (UTM) projection.
The temporal and spatial variety of precipitation is critical to understanding hydrologic processes [66]. There are seven precipitation gages in the study area. Four have data available from (MesoWest (http://mesowest.utah.edu)): (i.) Canelo RAWS (QCNA3), (ii.) Audubon Ranch Climate Reference Network (ADRA3), (iii.) Pioneer Airfield (KALK), and (iv.) a private weather station west of Sonoita near the watershed (E0586). Three other privately maintained weather stations were accessed from a cooperative rainfall monitoring network across Arizona (Rainlog.org) [67] (Figure 3). The SWAT weather generator (WGEN_US_COOP_1980_2010) was chosen to simulate other necessary data types (e.g., solar radiation, temperature, etc.).

(b) Iteration, Calibration, and Change Analysis

SWAT was run for the 4+ year study (01/01/2013-2/28/17), using the first 2 years as a warm up. The 411.6 km2 watershed was deliniated into 109 hydrologic response units (HRUs) and 25 sub-basins for processing. We ran a daily timestep using the measured precipitation described above, and calibrated volume of discharge in SWAT using data collected from the gaging station. Three quantitative statistics were used to evaluate the model’s performance of simulated data [68]. Once calibrated, we further checked it using the field data collected over the study period. The SWAT Calibration and Uncertainty Programs (SWAT-CUP) can provide a semi-automated approach using both manual and automated calibration and incorporating sensitivity and uncertainty analysis [69]—however, the autocalibration capability does not capture all information from key data types measured in our study and a manual calibration was carried out instead [70].
We used the model to evaluate alternative management scenarios for restoration, based on findings from our field campaign. SWAT has been applied in many scenario analyses of management practices to evaluate the effects of land-use changes on hydrological response in the Upper San Pedro Watershed [71] and the Santa Cruz Watershed [63]. We used the SWAT model to simulate how installing structures throughout the watershed would impact recharge and in turn, other factors in the hydrologic budget over the long term.

3. Results

This section is divided into three sections to relay (1) findings of field data collection and analysis, (2) results of change prediction from MODFLOW simulations, and (3) results of incorporating these into the SWAT model testing and change prediction results.

3.1. Data Collection and Analysis

3.1.1. Infiltration

To evaluate the impact of a single gabion on infiltration volume, a total of 225 scenarios were run using Gabion #2 dimensions, a range of saturated hydraulic conductivity, and a comparison of ponding patterns pre- and post-gabion installation (Figure 8). Post-monsoon hydraulic conductivity estimated by Fandel [54], using temperature methods, was lower upstream (25–54 mm/hr) than downstream (~100–400 mm/hr), likely due to clogging caused by ponding and consequent deposition of clay. Using these ranges and a 4.5-h flow, the gabion’s impact on infiltration ranged from no influence (gabion did not change ponding patterns) to an increase of 225% (assuming minimal pre-gabion ponding), with the most likely scenarios depicting an increase of approximately 10% [53,54,55]. Figure 8 portrays 6 scenarios with the same infiltration flux and ponded length, but different ponded width (based on time-lapse photography). Because no photo data were available for pre-gabion flow events, the no-gabion estimates were based on a linear decrease in ponded width over time and starting widths derived from channel morphology. The impact of the gabion on total infiltrated volume is the area under the curve with the gabion minus the area under the curve without the gabion [53,54,55].
Major sources of uncertainty in the resulting infiltrated volume estimates are the spatial distribution of bed-sediment hydraulic conductivity, the thermal properties of the bed-sediment, and the spatial extent of ponding over time. Uncertainty could be reduced by taking direct measurements of bed-sediment porosity, hydraulic conductivity, and thermal properties at each point-flux measurement location [53,54,55].

3.1.2. Erosion and Deposition

The LiDAR survey area of interest (Figure 5) depicted erosion (528.2 m2) and deposition (497.1 m2) equally. This likely reflects the additional construction and rehabilitation during the study. We noted significant scour downstream, and some deposition upstream. Several of the downstream gabions were undercut as well. As noted above, one likely result of upstream deposition is reduction of hydraulic conductivity. Although it was not evaluated, downstream scouring could also have resulted in altered hydraulic conductivity.

3.2. MODFLOW Change Prediction

Lacher [25] simulated change in head, or pressure, at our study site to simulate 500× natural recharge (936 m3/d) in MODFLOW, applied over 8 model cells (49 ha), for a simulation time period between March 2014–March 2030. The maximum simulated head change under this scenario was 15.5 m. Maximum head change under the river was roughly 2 m. To examine an extreme version of this scenario, we boosted model inputs to consider recharge at 1000× natural recharge and applied that to 4 model cells (24.5 ha), instead of 8, which resulted in the same annual recharge rate of 936 m3/d, and produced a groundwater mound (local water-level rise), with a maximum change in head of 26.5 m. Although the maximum head change was greater in this simulation than in the 8-cell simulation, the maximum head change under the river was only about 1.5 m by March 2030.
Neither of these simulations produced any change in baseflow in the Babocomari River by 2030, with the vadose zone underlying the streambed still slightly more than 2 m thick at the point where the groundwater mound from the recharge intersected near the confluence of Vaughn Canyon and the Babocomari. However, persistent increases in recharge at Vaughn would likely lead to a larger groundwater mound that would eventually intersect the ground surface and generate baseflow. Our ability to test this hypothesis is limited by the cell size (250 m × 250 m). The actual recharge effects from the restoration efforts are likely to be much more localized than what could be simulated with this model without modifying the grid structure.
Vaughn Canyon cells in MODFLOW have a horizontal K (HK) of 0.21 mm/h [28] with a vertical anisotropy of 5, meaning that HK is 5 times higher than vertical conductivity ~0.04 mm/h. The recharge simulations indicate that recharge rates at the study sites would have to be increased by 400- to 1,000-fold over a 60- to 121-acre area to cause increases in baseflow that are clearly distinguished from model “noise” by 2030. Concentrating recharge over smaller areas close to the river steepens the groundwater gradient (slope) toward the river, thereby raising heads under the river more quickly than with more dispersed recharge.
Taking the maximum increase in infiltration projected (~250%) for one gabion [54], and assuming that 100% of the infiltration becomes recharge, approximately 140 gabions would need to be installed to increase recharge by 500×, as in Lacher (25). If each gabion only increases infiltration by 10% [53], over 450 gabions would be needed.

3.3. SWAT

3.3.1. Calibration

In the Babocomari, the uncalibrated model predicts 14% of precipitation being directed to streamflow from the HRUs (total water yield), where 2/3 of that is “lateral flow”. A manual calibration was done, using approaches from other aridland applications, to improve the SWAT watershed model performance by first calibrating the model to reflect the local water balance and then the stream hydrograph.
Our calibration approach was based on a review of current literature that identified parameters commonly adjusted [11,62,63,70,72]. Sensitive surface-runoff and basin parameters were manually altered to mimic the larger San Pedro Watershed, documented by Yuan and Nie [72], consistent with other aridland watersheds [73] (Table 2). As with model calibration in the nearby Santa Cruz Watershed [63,74], the Chiricahua Mountains [11], and in the Upper San Pedro Watershed [72], the lateral flow component was reduced by nearly an order of magnitude from 9% to 1% of precipitation. This increases surface runoff and ET, as well as slightly increasing the total amount of water entering the aquifers (deep/shallow) via percolation past bottom of soil profile.
Part of surface runoff is referred to as baseflow, which comes from deep subsurface and delayed shallow surface flows. Yuan and Nie [72] eliminated baseflow from simulations by reducing the threshold water level in shallow aquifers for re-evaporation (REVAPMN) and enhancing the re-evaporation coefficient. In this study, the groundwater revap coefficient increased from default 0.02 (min) to 0.2 (max). This doubled the amount of water moving from shallow aquifer to plants/soil profile in the watershed during simulation with no other changes. We also reduced the ESCO coefficient to portray more of the evaporative demand from lower levels to its minimum potential (0), in turn reducing ET. To calibrate surface runoff, we edited the CN2 (SCS runoff curve number) to a low value (multiplied all by 0.5) to match observed water balance volumes. Partly based on research by Yuan and Nie [72], who document transmission losses for the major channel of the Upper San Pedro, we altered the effective channel-bed hydraulic conductivity (CH_K2) to better depict this. Secondarily, calibration was then focused to match the observed data’s hydrograph, as volume was the focus of this study (vs. timing). An informal sensitivity analysis was carried out by varying parameters to create a qualitative measure of parameter sensitivities [62,75]. As described by Niraula et al. [75], after testing each possible parameter, one at a time, and comparing output from the hydrograph to identify which of the identified parameters should be adjusted to mimic flow patterns and volumes portrayed by the USGS gage of the Babocomari, the most sensitive parameters were then manually adjusted to calibrate the model. The final calibration consisted of adjusting only the coefficients listed in Table 2.
Using these settings, the Pearson’s correlation coefficient (r = 0.46) and Coefficient of Determination (r2 = 21.33%) show low collinearity between simulated and observed data; however, these measurements might not be appropriate for daily data [68]. Moriasi et al. [68] suggest percent bias (PBIAS) and Nash-Sutcliffe efficiency (NSE) are more appropriate assessments for daily model output in hydrographs. The PBIAS measures the tendency of simulated values to be larger or smaller than observed ones. Low-magnitude values indicate accurate model simulation, where negative values such as found in this study (−4.89%) indicate model underestimation bias [76]. The NSE ~ 0.21, where values (0.0 > 1.0) demonstrate an acceptable level of performance [59,77].

3.3.2. Change Prediction

The calibrated SWAT model was used to generate predictions out until the year 2050, using SWAT’s weather generator (WGEN_US_COOP_1980_2010). To depict our change scenario in the model input, we increase the saturated hydraulic conductivity of all soil types (Sol_K) and effective hydraulic conductivity in the main channel alluvium (CH_K2) based on our field measured values ranging from 10% to 225% [53,54,55]. This implies many structures need to be installed to effectively increase infiltration (reference the MODFLOW simulations suggesting 140-455 gabions installed). It is important to note that the actual mechanism by which gabions increase infiltration is by slowing and retaining flow, increasing vegetative shading/cooling, ponding duration, groundwater storage potential due to increased fines and organic matter, and area, not by increasing K—which would likely decrease based on fine sediment deposition–but we developed this scenario to use the model to portray the increase in initial infiltration.
Given this scenario, the calibrated model predicts the average annual water budget basin values will be greatly altered (Figure 9). Surface runoff could potentially decrease (0–20%) and the lateral soil runoff could increase (9%–108%), pushing water underground. Groundwater (shallow aquifer) is predicted to increase (4%–42%) and groundwater (deep aquifer) could also increase (17%–50%). Total aquifer recharge could increase (4%–41%) and total water yield will increase (8%–92%). Percolation out of the soil will also increase (4%–42%; Figure 9).

4. Discussion

When water percolating the soil exceed the K of a low permeability layer, a perched water table and subsurface lateral flow can result [78]. Subsurface flow (also known as throughflow, subsurface storm flow, subsurface runoff, and interflow) originates when water infiltrates the soil, and tends to move preferentially laterally (due to anisotropy) through the upper soil horizons toward the stream as ephemeral, shallow, perched groundwater above the main groundwater level [60,78]. On the contrary, Smettem et al. [79] and Brouwer and Fitzpatrick [80] found that soil macroporosity and by-pass flow can prevent subsurface lateral flow. Coes and Pool [81] describe permeable soils in which, when pressurized with repeated infiltrating water, increased storage builds up water volume until it eventually reaches the water table. However, the development of lateral flow is poorly understood and is often ignored in field and modeling studies [78]; Jarvis et al. [82] suggest that more experiments need to be done to understand complex subsurface flow. If you increase surface flow, you may increase lateral flow into the bank and streambed, and thereby increase recharge. We compare this to findings around beaver dams that impact lateral and longitudinal connectivity by introducing roughness and heterogeneity elements that fundamentally change the timing, delivery, and storage of water, sediment, nutrients, and organic matter [17].
MODFLOW simulated recharge at the Vaughn Canyon site resulting in no downstream baseflow response by 2030, even at 1000× natural recharge rates over a 29.5-ha area. While Vaughn responds well to recharge because of its higher hydraulic conductivities, it is too far from the river to affect baseflows by 2030 [25]. The more focused the recharge and the closer to the river, the more immediate the baseflow response is predicted.
Using the SWAT model, it is possible to extrapolate the various changes in infiltration rates demonstrated in the field, resulting from ECS installation. An increase at Gabion #2 of hydraulic conductivity (~10 and 225%) was used to simulate our hypothesized change in the larger watershed. We found that projecting these changes to the year 2050 increases lateral flows that could also increase recharge to the aquifer, though it is unlikely that 100% of infiltrated water actually becomes recharge [66,81]. The general finding to drastically reduce the lateral flow component for aridland SWAT modeling is negated where ECS are installed. To further this modeling application, it would be recommended to test SWAT’s internal calibration tools and apply a more focused and systematic sensitivity analysis to improve the performance of the model, if timing of the hydrograph was a desired goal. An interesting finding of this research is that arid land watershed models should be run to mimic more temperate regions if ECS are installed in riparian channels.
We compared field data with SWAT input/output values, including K measurements, infiltration, and sedimentation data. Fandel et al. [53,54,55] measured Ks values across several gabions during the dry season, using a borehole permeameter. These measurements yielded an average Ks value of ~220 mm/hr. During a flow event, immediately upstream from Gabion #2, Ks was estimated in 1DTempPro from temperature data, yielding an average value of 40 mm/hr. The MODFLOW model used by Lacher [25] applied a HK of 0.21 mm/h and a vertical Ks of 0.04 mm/h. The SWAT model used the SSURGO database classification for the Pima series, with saturated Ks ranging 10-100 mm/h. Terrestrial-LiDAR surveys (July 12 and October 5, 2016) depict a gross volume of sediment deposition recorded between scans ~ 61.5 m3, which is ~123 metric tons over the period, around the structure itself. In the channel only, sedimentation was ~ 15.25 m3, which is ~30 metric tons over the period. In SWAT, there are 4 potential rainstorms that could have induced sediment movement in that timeframe: (i.) 8/3, (ii.) 8/8-10, (iii.) 9/6-7, and (iv.) a small event 9/29, during which, the uncalibrated SWAT model predicted a total of ~44 tons of sediment would be transported with water to the gabion, in close approximation of the T-LiDAR survey.
Rivers in semi-arid regions are less predictable than in more humid or temperate regions [72,74]. Our biggest limitation to the research was lack of groundwater data and measurements. All model inputs, such as the land-cover classifications, can impact model sensitivity to various degrees [83]. Further field studies are required to better understand the mechanisms responsible for the development of perched water tables and subsurface lateral flow in texture-contrast soils [78]. Despite data, model, and field experimentation limitations, we present the benefits of using this combination of field data and modeling to examine hydrologic response and allow for exploration and extrapolation.

5. Conclusions

This study helps quantify how Erosion-Control Structures, such as gabions, allow for increased infiltration in aridland environments and suggest associated impacts on the water budget. Because the study focuses primarily on shallow hydrologic processes, the groundwater model provided a predictive tool for investigating the potential hydrologic responses to long-term aggregated recharge at several locations which were being considered for surface-runoff projects. The simulated responses to hypothetical recharge values helped guide the selection of the surface-water project. Groundwater simulations indicate that recharge rates at the study sites would have to be dramatically increased to affect baseflow downstream. Our field data collection portrays increased average total infiltration volume (~+10%) at gabion installations. We altered the saturated hydraulic conductivity parameter in a hydrologic model to relate the change in soil-water flow rate to the hydraulic gradient and extend our findings throughout a watershed.
Model simulations adapting these findings demonstrate decreases in immediate runoff but extend lateral flow contribution to streamflow. Total groundwater recharge (both shallow and deep aquifer) is predicted to increase as will total water yield when restoration occurs. Given our change scenario of installing gabions throughout the watershed, percolation out of the soil will also increase, ultimately supporting volumes of flow downstream. By concentrating recharge over smaller areas close to the river, a strategically placed, dense network of gabions has the potential to steepen groundwater gradients, thereby raising heads more quickly than with more dispersed recharge. In areas where the groundwater table under the river is relatively shallow, the resulting groundwater mounding from this new source of recharge could potentially intersect the streambed surface and increase baseflow.
In the onset of the SWAT model iteration, we adopted standard aridland protocol to minimize lateral flows for calibration. However, our field studies reflect localized increased hydraulic conductivity at gabion structures that effectively reinstate conditions of lateral flow over long-term scenario testing. Terrestrial-LiDAR surveys validate model predictions of approximately equivalent sediments loads deposited at gabions. The potential to increase water flowing laterally within the soil profile, and eventually entering the main channel, is pinnacle to restoration practitioners in aridland environments.
Simulated results indicate that over time (20–30 years), gabion installation can help sustain precious mountain-front recharge, by forcing lateral flows, extending subsurface contribution to streamflow, and recharging shallow aquifers. The hydrologic model applications allow for visualization and planning of long-term impacts on water budgets and the multiple shared benefits of investing in riparian restoration, especially when those zones exist near major population centers. The results of this study could be used to formulate a proactive groundwater management plan for the study area to promote sustainable use of scarce water resources.

Author Contributions

Conceptualization, L.M.N. and J.B.C.; Methodology, L.M.N. and L.L.; Formal Analysis, L.M.N., J.B.C., L.L., and C.F.; Data Curation, L.M.N., J.B.C., L.L., N.R.W., B.T.F., and T.S.; Writing–Original Draft Preparation, L.M.N., J.B.C., N.R.W., and C.F.; Writing–Review & Editing, L.M.N., J.B.C., L.L., N.R.W. and C.F.; Project Administration, L.M.N., J.B.C., and L.L.; Funding Acquisition, L.M.N., L.L., and T.S.

Funding

Funding for the project was provided by the Walton Family Foundation, with support from the Land Change Science (LCS) Program, under the Land Resources Mission Area of the US Geological Survey (USGS). Swetnam was supported by NSF Awards DBI-0735191 and DBI-1265383.

Acknowledgments

We thank the Brophy family at the Babacomari Ranch, partners from Borderlands Restoration, LLC, and Cuenca los Ojos. We want to recognize the valuable assistance for SWAT modeling provided by Rewati Niraula and James Almendinger. We also appreciate the advice and field support provided by our colleagues at the USGS. Colleagues who helped develop and strengthen our research include from the Appleton-Whittell Research Ranch of the National Audubon Society and the University of Arizona. Lastly, we appreciate the careful peer reviews of this manuscript provided by David C. Goodrich and David Seibert. References to commercial vendors of software products or services are provided solely for the convenience of users when obtaining or using USGS software. Such references do not imply any endorsement by the US Government.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Chowdhury, A.; Jha, M.K.; Chowdary, V.M. Delineation of groundwater recharge zones and identification of artificial recharge sites in West Medinipur district, West Bengal, using RS, GIS and MCDM techniques. Environ. Earth Sci. 2010, 59, 1209–1222. [Google Scholar] [CrossRef]
  2. Renganayaki, P.; Elango, L. A review on managed aquifer recharge by check dams: A case study near Chennai, India. Int. J. Res. Eng. Technol. 2013, 2, 416–423. [Google Scholar]
  3. Bouwer, H. Artificial recharge of groundwater: hydrogeology and engineering. Hydrogeol. J. 2002, 10, 121–142. [Google Scholar] [CrossRef] [Green Version]
  4. Beechie, T.J.; Sear, D.A.; Olden, J.D.; Pess, G.R.; Buffington, J.M.; Moir, H.; Roni, P.; Pollock, M.M. Process-based Principles for Restoring River Ecosystems. BioScience 2010, 60, 209–222. [Google Scholar] [CrossRef] [Green Version]
  5. Martín-Rosales, W.; Gisbert, J.; Pulido-Bosch, A.; Vallejos, A.; Fernández-Cortés, A. Estimating groundwater recharge induced by engineering systems in a semiarid area (southeastern Spain). Environ. Geol. 2007, 52, 985–995. [Google Scholar] [CrossRef]
  6. Pandey, D.N.; Gupta, A.K.; Anderson, D.M. Rainwater harvesting as an adaptation to climate change. Curr. Sci. 2003, 85, 46–59. [Google Scholar]
  7. Norman, L.M.; Sankey, J.B.; Dean, D.; Caster, J.; DeLong, S.; DeLong, W.; Pelletier, J.D. Quantifying geomorphic change at ephemeral stream restoration sites using a coupled-model approach. Geomorphology 2017, 283, 1–16. [Google Scholar] [CrossRef]
  8. Naiman, R.J.; Johnston, C.A.; Kelley, J.C. Alteration of North American Streams by Beaver. BioScience 1988, 38, 753–762. [Google Scholar] [CrossRef]
  9. White, D.S. Biological relationships to convective flow patterns within stream beds. Hydrobiologia 1990, 196, 149–158. [Google Scholar] [CrossRef]
  10. Norman, L.M.; Brinkerhoff, F.; Gwilliam, E.; Guertin, D.P.; Callegary, J.; Goodrich, D.C.; Nagler, P.L.; Gray, F. Hydrologic Response of Streams Restored with Check Dams in the Chiricahua Mountains, Arizona. River Res. Appl. 2016, 32, 519–527. [Google Scholar] [CrossRef]
  11. Norman, L.M.; Niraula, R. Model analysis of check dam impacts on long-term sediment and water budgets in Southeast Arizona, USA. Ecohydrol. Hydrobiol. 2016, 16, 125–137. [Google Scholar] [CrossRef]
  12. Norman, L.M.; Huth, H.; Levick, L.; Shea Burns, I.; Phillip Guertin, D.; Lara-Valencia, F.; Semmens, D. Flood hazard awareness and hydrologic modelling at Ambos Nogales, United States–Mexico border. J. Flood Risk Manag. 2010, 3, 151–165. [Google Scholar] [CrossRef]
  13. Norman, L.M.; Levick, L.; Guertin, D.P.; Callegary, J.; Guadarrama, J.Q.; Anaya, C.Z.; Prichard, A.; Gray, F.; Castellanos, E.; Tepezano, E.; et al. Nogales Flood Detention Study; Open-File Rep. No. 2010–1262; US Geological Survey: Reston, VA, USA, 2010; p. 112.
  14. Norman, L.; Villarreal, M.; Pulliam, H.R.; Minckley, R.; Gass, L.; Tolle, C.; Coe, M. Remote sensing analysis of riparian vegetation response to desert marsh restoration in the Mexican Highlands. Ecol. Eng. 2014, 70, 241–254. [Google Scholar] [CrossRef] [Green Version]
  15. Wilson, N.R.; Norman, L.M. Analysis of vegetation recovery surrounding a restored wetland using the normalized difference infrared index (NDII) and normalized difference vegetation index (NDVI). Int. J. Remote Sens. 2018, 39, 3243–3274. [Google Scholar] [CrossRef] [Green Version]
  16. Nichols, M.H.; McReynolds, K.; Reed, C. Short-term soil moisture response to low-tech erosion control structures in a semiarid rangeland. CATENA 2012, 98, 104–109. [Google Scholar] [CrossRef]
  17. Macfarlane, W.W.; Wheaton, J.M.; Bouwes, N.; Jensen, M.L.; Gilbert, J.T.; Hough-Snee, N.; Shivik, J.A. Modeling the capacity of riverscapes to support beaver dams. Geomorphology 2017, 277, 72–99. [Google Scholar] [CrossRef]
  18. Gurnell, A.M. The hydrogeomorphological e•ects of beaver dam-building activity. Prog. Phys. Geogr. 1998, 22, 167–189. [Google Scholar] [CrossRef] [Green Version]
  19. Puttock, A.; Graham, H.A.; Cunliffe, A.M.; Elliott, M.; Brazier, R.E. Eurasian beaver activity increases water storage, attenuates flow and mitigates diffuse pollution from intensively-managed grasslands. Sci. Total Environ. 2017, 576, 430–443. [Google Scholar] [CrossRef] [Green Version]
  20. Arnold, J.G.; Muttiah, R.S.; Srinivasan, R.; Allen, P.M. Regional estimation of base flow and groundwater recharge in the Upper Mississippi river basin. J. Hydrol. 2000, 227, 21–40. [Google Scholar] [CrossRef]
  21. Sun, H.; Cornish, P.S. Estimating shallow groundwater recharge in the headwaters of the Liverpool Plains using SWAT. Hydrol. Process. 2005, 19, 795–807. [Google Scholar] [CrossRef]
  22. Kim, N.W.; Chung, I.M.; Won, Y.S.; Arnold, J.G. Development and application of the integrated SWAT–MODFLOW model. J. Hydrol. 2008, 356, 1–16. [Google Scholar] [CrossRef]
  23. McDonald, M.G.; Harbaugh, A.W. A Modular Three-Dimensional Finite-Difference Ground-Water Flow Model; U.S. Geological Survey: Reston, VA, USA, 1984.
  24. Norman, L.M. Surface Water Rainfall-Runoff Modeling at the Babacomari Watershed, SE Arizona, with Applications in GIS and RS; Unpublished Report; Walton Family Foundation: Washington, DC, USA, 2013. [Google Scholar]
  25. Lacher, L.J. Technical Memorandum Describing the Groundwater Modeling Study for the Babacomari Ranch Study Area; Lacher Hydrologic Consulting: Tucson, AZ, USA, 2013. [Google Scholar]
  26. Norman, L.M.; Lacher, L.; Seibert, D.; Pulliam, H.R.; Hare, T.; Austin, V.; Villarreal, M.L.; Gray, F.; Callegary, J.B. Delineation and Screening of recharge sites for installation of rock detention structures in the Babocomari River, a tributary of the San Pedro River. Presented at the Science on the Sonoita Plain, Appleton-Whittell Research Ranch, Elgin, AZ, USA, 7 June 2014. [Google Scholar]
  27. Baillie, M.N.; Hogan, J.F.; Ekwurzel, B.; Wahi, A.K.; Eastoe, C.J. Quantifying water sources to a semiarid riparian ecosystem, San Pedro River, Arizona. J. Geophys. Res. 2007, 112, 13. [Google Scholar] [CrossRef]
  28. Pool, D.R.; Dickinson, J.E. Ground-Water Flow Model of the Sierra Vista Subwatershed and Sonoran Portions of the Upper San Pedro Basin, Southeastern Arizona, United States, and Northern Sonora, Mexico; SIR-2006-5228; United States Geological Survey: Reston, VA, USA, 2007.
  29. Pool, D.R.; Coes, A.L. Hydrogeologic Investigations of the Sierra Vista Subwatershed of the Upper San Pedro Basin, Cochise County, Southeast Arizona; US Geological Survey: Reston, VA, USA, 1999.
  30. Freethey, G.W. Hydrologic Analysis of the Upper San Pedro Basin from the Mexico—United States International Boundary to Fairbank, Arizona; Open-File Rep. 82–752; United States Geological Survey: Reston, VA, USA, 1982; p. 52.
  31. Corell, S. Groundwater Flow Model Scenarios of Future Groundwater and Surface Water Conditions: Sierra Vista Subwatershed of the Upper San Pedro Basin—Southeastern Arizona; Modeling Report No. 10B; Arizona Department of Water Resources: Phoenix, AZ, USA, 1996.
  32. Lacher, L.J. Simulated Groundwater and Surface Water Conditions in the Upper San Pedro Basin, 1902–2105, Preliminary Baseline Results; Task 1 Report for December 2010 Contract; Lacher Hydrologic Consulting: Tucson, AZ, USA, 2011. [Google Scholar]
  33. Lacher, L.J. Simulated Near-Stream Recharge at Three Sites in the Sierra Vista Subbasin, Arizona: Tucson; Task 2-4 Report for December 2010 Contract; Lacher Hydrologic Consulting: Tucson, AZ, USA, 2012; 62p. [Google Scholar]
  34. Leake, S.A.; Gungle, B. Evaluation of Simulations to Understand Effects of Groundwater Development and Artificial Recharge on the Surface Water and Riparian Vegetation, Sierra Vista Subwatershed, Upper San Pedro Basin, Arizona; Open-File Rep. 2012–1206; United States Geological Survey: Reston, VA, USA, 2012.
  35. Harshbarger and Associates. Appendix 1—Consultant’s report on water development. In Report on Water Supply, Fort Huachuca and Vicinity, Arizona, by the U.S. Army Corps of Engineers; U.S. Army Engineers District: Los Angeles, CA, USA, 1974; pp. 1–33. [Google Scholar]
  36. Schwartzman, P.N. A Hydrogeologic Resource Assessment of the Lower Babocomari Watershed, Arizona; M.S. Hydrology and Water Resources, The University of Arizona: Tucson, AZ, USA, 1990. [Google Scholar]
  37. Gungle, B.; Callegary, J.B.; Paretti, N.V.; Kennedy, J.R.; Eastoe, C.J.; Turner, D.S.; Dickinson, J.E.; Levick, L.R.; Sugg, Z.P. Hydrological Conditions and Evaluation of Sustainable Groundwater Use in the Sierra Vista Subwatershed, Upper San Pedro Basin, Southeastern Arizona; Sci. Investig. Rep. 2016–5114, No. 90; United States Geological Survey: Reston, VA, USA, 2016.
  38. Schmerge, D.; Corkhill, F.; Flora, S. Water-Level Conditions in the Upper San Pedro Basin, Arizona; A.D.W.R. Water Level Change Map Series Report 3; Arizona Department of Water Resources: Phoenix, AZ, USA, 2006.
  39. Callegary, J.B.; Sosa, I.M.; Villaseñor, E.M.; dos Santos, P.; Saavedra, R.M.; Noriega, F.J.; Huth, A.K.; Gray, F.; Scott, C.A.; Megdal, S.; et al. San Pedro River Aquifer Binational Report; International Boundary and Water Commission: El Paso, TX, USA, 2016. [Google Scholar]
  40. Arizona Department of Water Resources. Upper San Pedro Basin Active Management Area Review Report; Arizona Department of Water Resources: Phoenix, AZ, USA, 2005.
  41. Sharma, V.; MacNish, R.; Maddock, T., III. Analysis of Hydrologic Data Collected by U.S. Bureau of Land Management 1987–1995 and Recommendations for Further Monitoring Programs; Fort Huachuca: Sierra Vista, AZ, USA, 1997. [Google Scholar]
  42. Hendrickson, D.A.; Minckley, W.L. Cienegas: Vanishing climax communities of the American Southwest. In Desert Plants USA; FAO: Washington, DC, USA, 1985. [Google Scholar]
  43. Robinett, D.; Kennedy, L. Babacomari River Riparian Protection Project; Final Report Submitted to: Arizona Water Protection Fund Commission 09-164WPF; US Department of Agriculture: Fort Collins, CO, USA, 2014.
  44. Radke, M. Beaver on the San Pedro River. Presented at the Wildlife and Threatened and Endangered Species Education Forum, Tucson, AZ, USA, 10 August 2013. [Google Scholar]
  45. Wick, M.A. Beaver Making an Arizona Comeback. In Eastern Arizona Courier. Available online: https://www.eacourier.com/ (accessed on 1 December 2018).
  46. Woo, M.-K.; Waddington, J.M. Effects of beaver dams on subarctic wetland hydrology. Arctic 1990, 43, 223–230. [Google Scholar] [CrossRef]
  47. Saksa, P. The Hydrology and Sediment Transport of Low-Gradient, Forested Headwater Streams. Master’s Thesis, Louisiana State University, Baton Rouge, LA, USA, 2007. [Google Scholar]
  48. Du Bray, E.A. (Ed.) Mineral Resource Potential and Geology of Coronado National Forest, Southeastern Arizona and Southwestern New Mexico; U.S. Geological Survey: Reston, VA, USA, 1996.
  49. Cook, J.P.; Youberg, A.; Pearthree, P.A.; Onken, J.A.; MacFarlane, B.J.; Haddad, D.E.; Bigio, E.R.; Kowler, A.L. Mapping of Holocene River Alluvium along the San Pedro River, Aravaipa Creek, and Babocomari River, Southeastern Arizona; Arizona Gelogical Survey: Tucson, AZ, USA, 2009.
  50. Brown, S.G.; Davidson, E.S.; Kister, L.R.; Thomsen, B.W. Water Resources of Fort Huachuca Military Reservation, Southeastern Arizona; Water Supply Paper 1819-D; U.S. Geological Survey: Reston, VA, USA, 1966.
  51. Brown, D.E. Biotic communities of the American Southwest: United States and Mexico [Western States (USA); Great Basin and Pacific Slope States]. In Desert Plants USA; FAO: Washington, DC, USA, 1982; Volume 4. [Google Scholar]
  52. Soil Survey Staff, Natural Resources Conservation Service, United States Department of Agriculture, Soil Survey Geographic (SSURGO). Database for Arizona. 2018. Available online: https://websoilsurvey.nrcs.usda.gov/ (accessed on 1 December 2018).
  53. Fandel, C.; Callegary, J.B.; Ferré, T.P.A.; Norman, L.M.; Scott, C.A. Evaluating the effect of gabions on vertical water flux in an ephemeral stream using wildlife cameras and temperature sensors. Presented at the 2015 Annual Conference of Society for Ecological Restoration—Southwest Chapter, Tucson, AZ, USA, 20–22 November 2015. [Google Scholar]
  54. Fandel, C.A. The Effect of Gabion Construction on Infiltration in Ephemeral Streams. Master’s Thesis, The University of Arizona, Tucson, AZ, USA, 2016. [Google Scholar]
  55. Fandel, C.; Callegary, J.B.; Ferré, T.P.A.; Norman, L.M.; Scott, C.A. Infiltration in ephemeral streams: Quantifying the effect of gabions on vertical water flux using wildlife cameras & temperature sensors. Presented at the Water Resources Research Center Annual Conference, Tucson, AZ, USA, 21 March 2016. [Google Scholar]
  56. Constantz, J. Heat as a tracer to determine streambed water exchanges. Water Resour. Res. 2008, 44. [Google Scholar] [CrossRef] [Green Version]
  57. Zhang, W.; Qi, J.; Wan, P.; Wang, H.; Xie, D.; Wang, X.; Yan, G. An Easy-to-Use Airborne LiDAR Data Filtering Method Based on Cloth Simulation. Remote Sens. 2016, 8, 501. [Google Scholar] [CrossRef]
  58. Arnold, J.G.; Srinivasan, R.; Muttiah, R.S.; Williams, J.R. Large Area Hydrologic Modeling and Assessment Part I: Model Development1. JAWRA J. Am. Water Resour. Assoc. 1998, 34, 73–89. [Google Scholar] [CrossRef]
  59. Gassman, P.W.; Reyes, S.H.; Green, C.H.; Arnold, J.G. The Soil and Water Assessment Tool: Historical development, applications, and future research directions. Trans ASABE 2007, 50, 1211–1250. [Google Scholar] [CrossRef]
  60. Neitsch, S.L.; Arnold, J.G.; Kiniry, J.R.; Williams, J.R. Soil and Water Assessment Tool: Theoretical Documentation, Version 2009; Texas Water Resources Institute: College Station, TX, USA, 2011. [Google Scholar]
  61. Sloan, P.G.; Moore, I.D. Modeling subsurface stormflow on steeply sloping forested watersheds. Water Resour. Res. 1984, 20, 1815–1822. [Google Scholar] [CrossRef] [Green Version]
  62. Niraula, R.; Norman, L.M.; Meixner, T.; Callegary, J. Multi-gauge Calibration for modeling the Semi-Arid Santa Cruz Watershed in Arizona-Mexico Border Area Using SWAT. Air Soil Water Res. 2012, 5, ASWR-S9410. [Google Scholar] [CrossRef]
  63. Niraula, R.; Meixner, T.; Norman, L.M. Determining the importance of model calibration for forecasting absolute/relative changes in streamflow from LULC and climate changes. J. Hydrol. 2015, 522, 439–451. [Google Scholar] [CrossRef]
  64. Almendinger, J.E.; Murphy, M.S.; Ulrich, J.S. Use of the Soil and Water Assessment Tool to Scale Sediment Delivery from Field to Watershed in an Agricultural Landscape with Topographic Depressions. J. Environ. Qual. 2014, 43, 9–17. [Google Scholar] [CrossRef] [PubMed]
  65. Wickham, J.; Homer, C.; Vogelmann, J.; McKerrow, A.; Mueller, R.; Herold, N.; Coulston, J. The Multi-Resolution Land Characteristics (MRLC) Consortium—20 years of development and integration of USA national land cover data. Remote Sens. 2014, 6, 7424–7441. [Google Scholar] [CrossRef]
  66. Goodrich, D.C.; Schmugge, T.J.; Jackson, T.J.; Unkrich, C.L.; Keefer, T.O.; Parry, R.; Bach, L.B.; Amer, S.A. Runoff simulation sensitivity to remotely sensed initial soil water content. Water Resour. Res. 1994, 30, 1393–1405. [Google Scholar] [CrossRef]
  67. Woodard, G.C.; Crimmins, M.; Vazquez, R.; Rupprecht, C. QA/QC Issues Related to Data from Volunteer Citizen Scientist Networks. In AGU Fall Meeting Abstracts; American Geophysical Union: Washington, DC, USA, 2007. [Google Scholar]
  68. Moriasi, D.N.; Arnold, J.G.; van Liew, M.W.; Bingner, R.L.; Harmel, R.D.; Veith, T.L. Model Evaluation Guidelines for Systematic Quantification of Accuracy in Watershed Simulations. Trans. ASABE 2007, 50, 885–900. [Google Scholar] [CrossRef]
  69. Abbaspour, K.C. SWAT-CUP 2012: SWAT Calibration and Uncertainty Programs–A User Manual; Swiss Federal Institute of Aquatic Science and Technology: Dübendorf, Switzerland, 2013; Volume 103. [Google Scholar]
  70. Arnold, J.G.; Moriasi, D.N.; Gassman, P.W.; Abbaspour, K.C.; White, M.J.; Srinivasan, R.; Santhi, C.; Harmel, R.D.; Van Griensven, A.; Van Liew, M.W.; et al. SWAT: Model Use, Calibration, and Validation. Trans. ASABE 2012, 55, 1491–1508. [Google Scholar] [CrossRef]
  71. Hernandez, M.; Miller, S.N.; Goodrich, D.C.; Goff, B.F.; Kepner, W.G.; Edmonds, C.M.; Jones, K.B. Modeling Runoff Response to Land Cover and Rainfall Spatial Variability in Semi-Arid Watersheds. In Monitoring Ecological Condition in the Western United States; Sandhu, S.S., Melzian, B.D., Long, E.R., Whitford, W.G., Walton, B.T., Eds.; Springer: Cham, The Netherlands, 2000; pp. 285–298. [Google Scholar]
  72. Yuan, Y.; Nie, W. Problems and Prospects of SWAT Model Application on an Arid/Semiarid Watershed in Arizona. In Proceedings of the 2015 SEDHYD Conference, Reno, NV, USA, 19–23 April 2015. [Google Scholar]
  73. Veith, T.L.; Van Liew, M.W.; Bosch, D.D.; Arnold, J.G. Parameter Sensitivity and Uncertainty in SWAT: A Comparison across Five USDA-ARS Watersheds. Trans. ASABE 2010, 53, 1477–1486. [Google Scholar] [CrossRef]
  74. Niraula, R.; Meixner, T.; Norman, L.M. Hydrological Modeling of a Semi-arid Santa Cruz Basin. Presented at the 2012 4th Annual Santa Cruz River Researchers’ Day, Tucson, AZ, USA, 29 March 2012. [Google Scholar]
  75. Niraula, R.; Kalin, L.; Wang, R.; Srivastava, P. Determining nutrient and sediment critical source areas with swat: Effect of lumped calibration. Trans. ASABE 2011, 55, 137–147. [Google Scholar] [CrossRef]
  76. Gupta, H.V.; Sorooshian, S.; Yapo, P.O. Status of Automatic Calibration for Hydrologic Models: Comparison with Multilevel Expert Calibration. J. Hydrol. Eng. 1999, 4, 135–143. [Google Scholar] [CrossRef]
  77. Nash, J.E.; Sutcliffe, J.V. River flow forecasting through conceptual models part I—A discussion of principles. J. Hydrol. 1970, 10, 282–290. [Google Scholar] [CrossRef]
  78. Hardie, M.A.; Doyle, R.B.; Cotching, W.E.; Lisson, S. Subsurface Lateral Flow in Texture-Contrast (Duplex) Soils and Catchments with Shallow Bedrock. Appl. Environ. Soil Sci. 2012, 2012, 1–10. [Google Scholar] [CrossRef] [Green Version]
  79. Smettem, K.R.J.; Chittleborough, D.J.; Richards, B.G.; Leaney, F.W. The influence of macropores on runoff generation from a hillslope soil with a contrasting textural class. J. Hydrol. 1991, 122, 235–251. [Google Scholar] [CrossRef]
  80. Brouwer, J.; Fitzpatrick, R.W. Restricting layers, flow paths and correlation between duration of soil saturation and soil morphological features along a hillslope with an altered soil water regime in western Victoria. Aust. J. Soil Res. 2002, 40, 927. [Google Scholar] [CrossRef]
  81. Coes, A.L.; Pool, D.R. Ephemeral-Stream Channel and Basin-Floor Infiltration and Recharge in the Sierra Vista Subwatershed of the Upper San Pedro Basin, Southeastern Arizona; Open-File Rep. 2005–1023; US Geological Survey: Reston, VA, USA, 2005; p. 84.
  82. Jarvis, N.; Koestel, J.; Larsbo, M. Understanding Preferential Flow in the Vadose Zone: Recent Advances and Future Prospects. Vadose Zone J. 2016, 15. [Google Scholar] [CrossRef] [Green Version]
  83. Miller, S.N.; Guertin, D.P.; Goodrich, D.C. Hydrologic Modeling Uncertainty Resulting from Land Cover Misclassification1. JAWRA J. Am. Water Resour. Assoc. 2007, 43, 1065–1075. [Google Scholar] [CrossRef]
Figure 1. Photograph of one of the Vaughn Canyon gabions (Photo by L. Norman 9/11/15).
Figure 1. Photograph of one of the Vaughn Canyon gabions (Photo by L. Norman 9/11/15).
Water 11 00381 g001
Figure 2. Location map depicting the San Pedro Watershed, in the State of Arizona and the US–Mexico border, cities, land ownership, and the Babocomari HUC10 sub-watershed.
Figure 2. Location map depicting the San Pedro Watershed, in the State of Arizona and the US–Mexico border, cities, land ownership, and the Babocomari HUC10 sub-watershed.
Water 11 00381 g002
Figure 3. Map depicting the Upper Babocomari Watershed, hydrology, dams, rain gages, and gaging station mapped for this study.
Figure 3. Map depicting the Upper Babocomari Watershed, hydrology, dams, rain gages, and gaging station mapped for this study.
Water 11 00381 g003
Figure 4. Flow chart describing phases of research and hydrologic model applications.
Figure 4. Flow chart describing phases of research and hydrologic model applications.
Water 11 00381 g004
Figure 5. Diagram portraying location of instruments around Gabion #2 (drawing by C. Fandel).
Figure 5. Diagram portraying location of instruments around Gabion #2 (drawing by C. Fandel).
Water 11 00381 g005
Figure 6. Photo taken 5 min after flow onset. Guidelines were superimposed later to assist in estimating ponding extent by C. Fandel.
Figure 6. Photo taken 5 min after flow onset. Guidelines were superimposed later to assist in estimating ponding extent by C. Fandel.
Water 11 00381 g006
Figure 7. Image depicting the LiDAR-derived data collected at gabion #2 in Vaughn Canyon (Image by B. Forbes).
Figure 7. Image depicting the LiDAR-derived data collected at gabion #2 in Vaughn Canyon (Image by B. Forbes).
Water 11 00381 g007
Figure 8. Selected high, intermediate, and low infiltration scenarios with and without a gabion present.
Figure 8. Selected high, intermediate, and low infiltration scenarios with and without a gabion present.
Water 11 00381 g008
Figure 9. Chart portraying change predictions for water budget in 2050 if gabions are installed throughout the Babocomari Watershed.
Figure 9. Chart portraying change predictions for water budget in 2050 if gabions are installed throughout the Babocomari Watershed.
Water 11 00381 g009
Table 1. Daily discharge, cubic meters per second—statistics based on 16 years of record at Upper Babocomari River (near Huachuca City, AZ) discharge gaging station (#09471380; https://nwis.waterdata.usgs.gov/ accessed on 18 April 2017).
Table 1. Daily discharge, cubic meters per second—statistics based on 16 years of record at Upper Babocomari River (near Huachuca City, AZ) discharge gaging station (#09471380; https://nwis.waterdata.usgs.gov/ accessed on 18 April 2017).
Min–200525th PercentileMedianMean75th percentileMax–2001
0.0050.0130.0240.0620.1050.266
Table 2. Parameters used to calibrate the Upper Babocomari Watershed Model, where variant ranges of default values depend on input.
Table 2. Parameters used to calibrate the Upper Babocomari Watershed Model, where variant ranges of default values depend on input.
ParameterDescriptionFileDefault ValueRangeCalibrated Value
CN2SCS runoff curve number for moisture condition II.mgtVaries×0.5
ESCOSoil evaporation compensation factor(.hru) .bsn0.950–10.75
REVAPMNThreshold water level in shallow aquifer for revap.gw7500–10001000
SOL_AWCAvailable water capacity of the soil layer.sol0–10–1×2.1
CH_K2Effective hydraulic conductivity of channel.rte0(-0.01–500)2
GW_RevapRevaporation coefficient.gw0.020.02–0.20.2
GWQMINDeep percolation loss.gw10000–5000100
GW_DELAY (days)Groundwater delay time.gw310–500100
SURLAGSurface-runoff lag coefficient.bsn410

Share and Cite

MDPI and ACS Style

Norman, L.M.; Callegary, J.B.; Lacher, L.; Wilson, N.R.; Fandel, C.; Forbes, B.T.; Swetnam, T. Modeling Riparian Restoration Impacts on the Hydrologic Cycle at the Babacomari Ranch, SE Arizona, USA. Water 2019, 11, 381. https://doi.org/10.3390/w11020381

AMA Style

Norman LM, Callegary JB, Lacher L, Wilson NR, Fandel C, Forbes BT, Swetnam T. Modeling Riparian Restoration Impacts on the Hydrologic Cycle at the Babacomari Ranch, SE Arizona, USA. Water. 2019; 11(2):381. https://doi.org/10.3390/w11020381

Chicago/Turabian Style

Norman, Laura M., James B. Callegary, Laurel Lacher, Natalie R. Wilson, Chloé Fandel, Brandon T. Forbes, and Tyson Swetnam. 2019. "Modeling Riparian Restoration Impacts on the Hydrologic Cycle at the Babacomari Ranch, SE Arizona, USA" Water 11, no. 2: 381. https://doi.org/10.3390/w11020381

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop