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

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.


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.
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.7m 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 subwatershed, 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].
Water 2018, 10, x FOR PEER REVIEW 3 of 20

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 m 3 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 m 3 /s by 2030.Lacher [32] also simulated up to 0.03 m 3 /s of pumping-induced declines in baseflow in the upper half of the Babocomari over the 21st century.This tributary contributes about 7.4M m 3 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 m 3 (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 m 3 /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  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.

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.
Water 2018, 10, x FOR PEER REVIEW 6 of 20 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.

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.

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.

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 (Figures 5 and 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.
Water 2018, 10, x FOR PEER REVIEW 7 of 20

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.

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 (Figures 5 and 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.

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.

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.
. Figure 6.Photo taken 5 min after flow onset.Guidelines were superimposed later to assist in estimating ponding extent by C. Fandel.

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.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/m 3 ).

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.

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 m 3 /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.

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 km 2 ; 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.
(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 km 2 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.

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.

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].
Water 2018, 10, x FOR PEER REVIEW 11 of 20 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.

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.

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 pregabion 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].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].

Erosion and Deposition
The LiDAR survey area of interest (Figure 5) depicted erosion (528.2 m 2 ) and deposition (497.1 m 2 ) 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.

MODFLOW Change Prediction
Lacher [25] simulated change in head, or pressure, at our study site to simulate 500× natural recharge (936 m 3 /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 m 3 /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.

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 (r 2 = 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].

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).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).

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-

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 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 m 3 , which is ~123 metric tons over the period, around the structure itself.In the channel only, sedimentation was ~15.25 m 3 , 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.

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.

Figure 2 .
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 .
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 3 .
Figure 3. Map depicting the Upper Babocomari Watershed, hydrology, dams, rain gages, and gaging station mapped for this study.

Figure 3 .
Figure 3. Map depicting the Upper Babocomari Watershed, hydrology, dams, rain gages, and gaging station mapped for this study.

Figure 4 .
Figure 4. Flow chart describing phases of research and hydrologic model applications.Figure 4. Flow chart describing phases of research and hydrologic model applications.

Figure 4 .
Figure 4. Flow chart describing phases of research and hydrologic model applications.Figure 4. Flow chart describing phases of research and hydrologic model applications.

Figure 6 .
Figure 6.Photo taken 5 min after flow onset.Guidelines were superimposed later to assist in estimating ponding extent by C. Fandel.

Figure 7 .
Figure 7. Image depicting the LiDAR-derived data collected at gabion #2 in Vaughn Canyon (Image by B. Forbes).

Figure 6 .
Figure 6.Photo taken 5 min after flow onset.Guidelines were superimposed later to assist in estimating ponding extent by C. Fandel.

Figure 7 .
Image depicting the LiDAR-derived data collected at gabion #2 in Vaughn Canyon (Image by B. Forbes).

Figure 7 .
Figure 7. Image depicting the LiDAR-derived data collected at gabion #2 in Vaughn Canyon (Image by B. Forbes).

Figure 8 .
Figure 8. Selected high, intermediate, and low infiltration scenarios with and without a gabion present.

Figure 8 .
Figure 8. Selected high, intermediate, and low infiltration scenarios with and without a gabion present.

Water 2018 ,
10, x FOR PEER REVIEW 14 of 20

Figure 9 .
Figure 9. Chart portraying change predictions for water budget in 2050 if gabions are installed throughout the Babocomari Watershed.

Figure 9 .
Figure 9. Chart portraying change predictions for water budget in 2050 if gabions are installed throughout the Babocomari Watershed.

Table 2 .
Parameters used to calibrate the Upper Babocomari Watershed Model, where variant ranges of default values depend on input.