In most watershed-modeling studies, flow is calibrated at one monitoring site, usually at the watershed outlet. Like many arid and semi-arid watersheds, the main reach of the Santa Cruz watershed, located on the Arizona-Mexico border, is discontinuous for most of the year except during large flood events, and therefore the flow characteristics at the outlet do not represent the entire watershed. Calibration is required at multiple locations along the Santa Cruz River to improve model reliability. The objective of this study was to best portray surface water flow in this semi-arid watershed and evaluate the effect of multi-gauge calibration on flow predictions. In this study, the Soil and Water Assessment Tool (SWAT) was calibrated at seven monitoring stations, which improved model performance and increased the reliability of flow predictions, in the Santa Cruz watershed. The most sensitive parameters to affect flow were found to be curve number (CN2), soil evaporation and compensation coefficient (ESCO), threshold water depth in shallow aquifer for return flow to occur (GWQMN), base flow alpha factor (ALPHA_BF), and effective hydraulic conductivity of the soil layer (CH_K2). In comparison, when the model was established with a single calibration at the watershed outlet, flow predictions at other monitoring gauges were inaccurate. This study emphasizes the importance of multi-gauge calibration to develop a reliable watershed model in arid and semi-arid environments. The developed model, with further calibration of water quality parameters will be an integral part of the Santa Cruz Watershed Ecosystem portfolio Model (SCWEPM), an online decision support tool, to assess the impacts of climate change and urban growth in the Santa Cruz watershed.
Introduction
Watershed models are used to analyze the effects of land use and climate change on water resources and provide important information for decision making. These models should be calibrated and validated with available data.1 Careful consideration should be taken during model calibration, because results can vary widely based on calibration method. Even more attention in calibration should be given when dealing with watersheds in arid and semi-arid regions where flow is not continuous. Watershed in these regions are characterized by relatively larger extremes for hydrological components including low annual precipitation but high intensity storms, high evaporation, low baseflow but high flash floods, and runoff loss through stream bed.2
The Soil and Water Assessment Tool (SWAT3) was originally developed to predict the impact of management practices on water, sediment and agricultural chemicals in watersheds comprising different soils, land use and management conditions over long time periods. SWAT has been used for studying the impact of land use change4–567 and climate change8–910 as well as their combined/interactive effect11–1213 on water resources. SWAT is also coupled with ecosystem models to study the tradeoffs between water quality and ecosystem services.14–15 SWAT has also been used for modeling arid and semi-arid watersheds throughout the world.16–17181920
This paper presents a careful calibration of SWAT for the Santa Cruz watershed to provide the best representation of surface-water flow and hydrologic processes. In the US-Mexico border region of the desert Southwest, where ecosystem health is dependent on limited water resources, the model will be used to assess contaminant and sediment transport, a step towards identifying risk to water resources for the US Geological Survey (USGS) US Mexico Border Health Initiative (BEHI).21 Products of the research will also be utilized in the USGS Santa Cruz Watershed Ecosystem Portfolio Model (SCWEPM), an online decision support tool, to promote the use of information relevant to water allocation and land management.22
The Santa Cruz watershed, located in southern Arizona and northern Sonora, Mexico, is 9,100 km2 with spatially heterogeneous climate, land use/land cover, soils, and elevation. Because the main reach is discontinuous, except in a few large flood events, the flow characteristics at the outlet do not represent the entire watershed. The objectives of this study were to investigate the effect of multi-gauge calibration on flow prediction in a semi-arid watershed and develop a reliable model to analyze the effects of various future scenarios of land use and climate change in the Santa Cruz watershed.
Study Area
The Santa Cruz River (Fig. 1), a tributary of the Colorado River, originates in southern Arizona and flows south towards Mexico before re-entering the United States. Most of the upper Santa Cruz watershed is located in southern Arizona (85%); the remaining 15% is located in northern Mexico. The length of the main reach from its headwaters in Arizona's San Rafael Valley to the watershed outlet (Fig. 1) is about 212 km. Based on a land-use map depicting the watershed in 1999 (Fig. 2),23 the watershed is dominated by shrub land (~70%) followed by forested land (~15%) and developed land (10%). Elevations in the watershed range from 610 m to 2,884 m (Fig. 1). Soils in the upper Santa Cruz watershed are dominated by loam and sandy loam (Fig. 3). Although some sections of the river are perennial because of point-source discharge from wastewater treatment plants (Fig. 1), the majority of the river's length is ephemeral. The climate of the Santa Cruz watershed is distinguished by mild-winter and high-summer temperatures. Precipitation distribution follows a bimodal pattern in which the majority falls during mid-summer and winter months.23 Based on 48 years of data (1960-2007), the average annual precipitation ranges between 200 mm and 800 mm. The average annual precipitation over the whole watershed is approximately 430 mm. The average annual temperature is 66 °F. While, the average annual maximum temperature is 82 °F, the annual minimum temperature is 49 °F.
Methods
This study is focussed on modeling a semi-arid watershed. Results from this study are applicable to similar systems. Methodology described for this hydrological modeling study cover details on SWAT model set up, sensitivity analysis, model calibration and validation. The prediction performance with a single outlet calibrated model is compared with a spatially calibrated model.
Figure 1.
Map portraying the Santa Cruz watershed, elevation changes, streams, point locations of NOAA weather stations USGS stream-gaging stations, and WWTP point sources.

Relevance of Swat model
The SWAT model was developed for application in large watersheds, uses readily available inputs, is computationally efficient, and provides the capability to study long-term impacts of perturbations (eg, climate and land use change). SWAT simulates potential nonpoint source and point source pollutants to track contaminants and can also be used to investigate the costs and benefits, in terms of ecosystem services, of various future scenarios of land use and climate change.21,22
Soil and Water Assessment Tool(Swat)description
SWAT is a continuous-time simulation, semi-distributed, quasi-process-based watershed model.24 The major model inputs are topography, soil properties (such as texture, soil erodibility, hydraulic conductivity, hydrologic soil group, soil depth, organic matter content, available water capacity), land use/cover type, weather/climate, and land management practices. The SWAT model subdivides the watershed into several sub-watersheds, which are further divided into hydrological response units (HRUs) according to topography, land use, and soil. Surface runoff from daily precipitation is estimated using a modification of the Soil Conservation Service (SCS) curve number (CN) method.25,26 Daily runoff can also be estimated using sub-daily precipitation using the Green-Ampt method.26 SCS CN method was used in this study. Runoff from all HRUs in the sub-watershed yields total sub-watershed discharge. Flow in SWAT is routed through channels using Muskingum routing method or variable storage coefficient method.26 The latter was used in this study. Estimates of potential evapotranspiration (PET) are made using the modified Penman Montieth method, Hargreaves method or Priestley-Taylor method.26 Penman Montieth method was used in this study. A Kinematic Storage method is used to predict lateral flow. In SWAT, erosion and sediment yields from each HRU are predicted based on the modified soil loss equation (MUSLE).27 Channel sediment is routed based on the modified Bagnold's sediment transport equation.28 SWAT also has Nitrogen (N) and Phosphorus (P) process representations in which mineralization, decomposition, and immobilization are represented. Inorganic N and P are calculated in surface runoff, lateral flow and groundwater flow as products of volume and average concentration. Organic N and P transport are calculated with sediment and runoff using a loading function.29 The SWAT theoretical manual26 provides a detailed description of the processes that are incorporated in SWAT. There are over 900 published papers on the development and application of SWAT.30
Figure 2.
Land use distribution map of the Santa Cruz watershed in 1999 reclassified from Villarreal et al.23

Data compilation
Input data required for SWAT were collected (Table 1) and formatted as required by SWAT. The digital elevation model (DEM), land use/land cover, and soil map for the US and Mexico were joined to make a single input layer for SWAT. A 30-m land use/land cover map was developed for the Santa Cruz watershed.23 The US General Soil Map (STATSGO2) were used to represent soils in the US and joined with the Food and Agriculture Organization of the United Nations (FAO) soils dataset that described soils south of the border. Higher resolution Soil Survey Geographic (SSURGO) data was available for the US portion of the watershed, but not adapted to this study. Although higher resolution soil data impacts the distribution of hydrological response units (HRUs) in the SWAT model, it does not lead to obvious improvements of stream flow simulation in large watersheds.37 Observed stream discharge data from 7 USGS stream-gaging stations (Table 2) was compiled to facilitate model calibration and validation.
Model Set up
SWAT model for the Santa Cruz watershed was set up using the GIS interface of SWAT (ArcSWAT). The latest version of ArcSWAT, ie, ArcSWAT 2009.93.7a31 was used for model parameterization in this study. The model development process sequentially involves watershed delineation, land use/soil/slope reclassification, HRUs creation, weather data incorporation, addition of WWTP discharge, and writing input files and model simulation.
Watershed delineation
The watershed outlet was defined at Cortaro gauge (Fig. 1) where the Santa Cruz River exits the city limits of Tucson. The DEM was used for generating flow paths and flow direction in a Geographic Information System (GIS). The DEM was masked with the predefined watershed boundary22 to delineate a stream network in the Santa Cruz watershed. A 50 km2 threshold was used to define the origin of a stream. Locations of stream-gaging stations, water quality sampling points, and the waste water treatment plants (WWTPS) obtained from the USGS, Arizona Department of Environmental Quality (ADEQ) and Arizona Department of Water Resources (ADWR) were documented in a separate database table (.dbf), to be integrated as “outlet location coordinates” for generating specific sub-watersheds for later analysis. A total of 131 sub-watersheds were created for the Santa Cruz watershed (Fig. 1).
Table 1
Sources of data used in setting up the SWAT model for the Santa Cruz watershed.

Land Use/soil/slope classification
Land use, soil, and slope inputs are required to be reclassified by SWAT into database codes, so that their properties can be used later by SWAT. Land use conditions from year 199923 was reclassified into 13 land use/land cover classes using NLCD 2001 lookup table incorporated in SWAT, (Table 3). The 1999 land-use map was selected because the observed-flow data available was from the same time period (Fig. 2). Likewise, the soil map was reclassified into 23 different SWAT soil classes (Fig. 3; Tables 4 and 5). The watershed was reclassified into three slope classes (0%-5%, 5%-10% and >10%; Table 6).
Table 2
Seven USGS stream-gaging stations in the Santa Cruz watershed used in this study.

Hydrologic Response Units(Hrus)creation
Next, the land use, soil, and slope layers were used to create a unique combination or HRU within each sub-watershed. HRUs are lumped land areas within the sub-watershed that comprise a unique combination of these data. A total of 7,702 HRUs were created for the Santa Cruz watershed. A HRU threshold of 1% for land use, 5% for soil and 10% for slope was then used to reduce the final number of HRUs to 2,583. A threshold is used at a sub-watershed scale so a sub-watershed that contains a land cover, soil or slope smaller than the threshold will be incorporated as the dominant land cover, soil or slope closest to it. For most of the surface hydrological processes, SWAT analyzes each HRU individually. The number of HRUs in a sub-watershed ranged from 2 to 53 based on the area and heterogeneity of each subwatershed. The area of HRUs ranges from 0.002 km2 to 18.8 km2.
Weather Data incorporation
Climate data (daily values of precipitation, minimum and maximum temperature) from the National Oceanic and Atmospheric Administration (NOAA) weather stations35 within and up to 50 km outside the watershed boundary (Fig. 1) were obtained and analyzed. Missing records were estimated by interpolation from neighboring weather stations and the Parameter-elevation Regressions on Independent Slopes Model (PRISM).38,39 A daily dataset from 1960 to 2007 was used as input for both precipitation and temperature. In SWAT, each sub-watershed uses data from the nearest station as input.
Table 3
SWAT land use classes.

Addition of Discharge from Waste Water Treatment Plants(Wwtps)
The Santa Cruz watershed has sewage treatment plant discharges of water, sediment and nutrients, to be incorporated into the model. In order to account for the loadings from these point sources, average daily loading data (discharge volume in cubic meters) were added for the WWTPS at Nogales, Ina Road, and Roger Road to the point-source location database file in respective sub-basin locations. Continuous discharge from these WWTPS is responsible for perennial flow in several segments in the Santa Cruz watershed and SWAT routes these through the channel network and integrates them with the loadings generated from the land areas.26
Model simulation
Based on continuity of observed data available for calibration and validation, a 26-year period (1982 to 2007) was selected for simulation. The model was executed at a daily time step but analyzed at a monthly time scale. The first five years of the simulation were used to condition the model in order to minimize uncertainties due to unknown initial conditions, particularly antecedent soil moisture conditions. A conditioning period of at least 1 year and up to 5 years for the SWAT model is recommended to initialize and develop reasonable starting values for model variables.40,41
Table 4
Soil based on texture classes.

Simulation results
The simulated monthly average flows were compared with the observed data from 7 gauges inside the watershed (Fig. 1) using Percent Bias (PBIAS) as a quantitative measure. PBIAS measures the average tendency of the simulated data to be greater or less than the observed flow data. The optimal value of PBIAS is 0, with values near zero representing model simulations with minimal mass balance errors. Positive values indicate model overestimation of discharge, and negative values indicate model underestimation of discharge (Equation 1).42
where O and S are observed and simulated values, respectively.The simulated and observed flow deviated considerably with the model portraying an extreme overestimation in the uncalibrated state (Table 7). A sensitivity analysis was then done to find the most sensitive parameters for model calibration.
Sensitivity analysis
A sensitivity analysis was done at the stream-gaging station closest to the outlet (Cortaro). This analysis identified the most influential model parameters for simulating the observed data. The method incorporated in the ArcSWAT interface combines the Latin Hypercube (LH) and One-factor-at-a-Time (OAT) sampling.43,44 This method was used for sensitivity analysis in this study. With this sensitivity analysis method, SWAT runs (p + 1)*m times, where p is the number of parameter and m the number of loops. Then, for each loop, a parameter set is selected such that a unique area of the parameter space is sampled, and used for running baseline simulation. Next, selecting one at a time, each parameter is randomly selected and its value is changed by a defined percentage and a simulation is executed. Once all the parameters have been varied, the algorithm then locates a new sampling area by changing all the parameters.45,46 In this way, using 26 parameters and the default of 10 loops, with parameters changing by default 5%, the model was run 270 times to identify important parameters (Table 8).
Table 5
SWAT soil classes.

Model Calibration and validation
The model was calibrated using a manual calibration approach for flow calibration. A “split sampling” method was adapted for calibration/validation where roughly one half of the observed flow data for the available period was used for calibration, and the remaining half was used for validation.24 In this study, Cortaro, the stream gaging station closest to the watershed outlet, is referred to as the watershed outlet. The model was first calibrated and validated at Cortaro (Fig. 1). Then, predicted flow was validated at all other locations to make sure the model was representing the watershed as a whole. Although the model was validated at Dodge, which was relatively closer to the outlet (Fig. 1), the model failed to validate at all of the other gauge stations, which implied that the model was not representing the entire watershed.
Table 6
SWAT slope classes.

Table 7
Table comparing simulated vs. observed flow before model calibration, using PBIAS%.

Table 8
Flow parameters used for sensitivity analysis, their definitions, rank, and mean change in objective function (SSE) due to 5% change in value.

Thus the parameters were further adjusted in a multi-gauge calibration at 6 other gauges inside the watershed. Calibration was initiated from the upstream gauge, starting at Nogales and then sequentially calibrating at Tubac, Continental, and Tucson (Fig. 1). Within the Rillito Creek watershed, the model was first calibrated at Vail and then at Dodge (Fig. 1). Finally, the flow was re-calibrated at Cortaro. Parameters were adjusted (Tables 9 and 10) in the areas not yet calibrated, to develop a single model for the entire Santa Cruz watershed. It should be noted that some parameters which were found to be less sensitive yet are preferred over more sensitive ones, because they improve the shape of the hydrograph, and runoff and baseflow ratios.
The performance of the calibration was evaluated qualitatively based on visual comparison of time series plots and flow duration curves (FDCs) and quantitatively based on the coefficient of determination (R2), Nash-Sutcliffe Efficiency (NSE), and Percent Bias (PBIAS), three common model evaluation measures.42 NSE indicates how well the plot of observed versus simulated data fits the 1:1 line. NSE is computed as:
where, O and S are observed and simulated values, respectively.Ō is the mean of observed values. NSE values range between -∞ and 1 with NSE of 1 indicating a perfect simulation. Simulation results are often considered to be satisfactory when NSE is greater than 0.5.42 The coefficient of determination (R2) is a measure of collinearity between observed and simulated data, and ranges between 0 and 1 (Equation 3).19
Although R2 > 0.5 is acceptable for modeling, a higher value is considered better.42
Results and Discussions
Sensitivity analysis
The sensitivity analysis performed at Cortaro provided an understanding of which processes and their associated parameters in SWAT were primarily controlling the hydrology of the Santa Cruz watershed, particularly with respect to flow at the outlet. The curve number (CN2) was found to be the most sensitive parameter followed by soil evaporation and compensation coefficient (ESCO), threshold water depth in shallow aquifer for return flow to occur (GWQMN), Baseflow alpha factor (ALPHA_BF), and effective hydraulic conductivity of the soil layer (CH_K2). Runoff was depicted to be controlled by curve number (CN2) because the SCS curve number method was used reflecting that the characteristics of land and soil control runoff generation. Curve number is a function of soil, land use and antecedent moisture content, and is an important parameter for calibrating flow in any region and in any model that applies this runoff generation method.
Table 9
Adjusted parameters and their respective default values.

The soil evaporation compensation factor (ESCO) controls the evapotranspiration from soil which allows for modifying the depth distribution used to meet the soil evaporative demand, to account for the effect of capillary action, crusting and cracks. The baseflow recession constant (ALPHA_BF) is a direct index of groundwater flow response to changes in recharge; it describes the rate at which stream-flow decreases when the stream channel is recharged by groundwater. The baseflow recession constant is very important in a semi-arid watershed, where the streamflow recession is quick. Effective channel hydraulic conductivity controls the loss of water through the streambed; ie, the transmission loss from the stream bed is a controlling mechanism for flow in this system. Effective hydraulic conductivity of the stream is set to 0 by default in the SWAT model, meaning no loss from stream bed is expected, which needs to be increased in semi-arid basins. These parameters were also found to be sensitive in Walnut Gulch watershed, a subwatershed of the neighboring San Pedro watershed in southern Arizona.47 Research showed that Walnut Gulch, which has the lowest mean monthly discharge by an order of magnitude compared with other watersheds used in that study, has a wide range of variation in hydrograph uncertainty, particularly for high flows (-18.2% to +86.7%). This emphasizes the need for careful consideration, and perhaps model adjustments, when evaluating a watershed with long-term low flows and flashy storm responses. The adjustment of identified parameters in the Santa Cruz watershed likewise helped improve the flow prediction during calibration.
Calibration and validation
At both of the headwater gauges of the Santa Cruz watershed, Nogales and Vail (Fig. 1), streamflow is more persistent than at downstream gauges. Despite good model performance during the calibration period (PBIAS = -4%; Table 9), peak flows were consistently underestimated during the validation period (PBIAS = -43%; Table 11) at Vail (Fig. 4E). FDCs at Vail showed that the model underestimated high flows and low flows and overestimated intermediate flows (Fig. 4E). At Nogales, the model was able to predict peak flow and low flow conditions (Fig. 4A), but under-predicted flow by 18% during the overall calibration period and over-predicted by 26% during the validation period (Table 11). In particular, model-predicted and observed flow did not match for December 2007 at Nogales (Fig. 4A), when the model simulated high flow but there was no observed flow. The model mostly underestimated low flows (Fig. 5A).
Table 10
Adjusted parameters and their respective adjusted values for model calibration.

The river is perennial at Tubac because of the continuous supply of treated effluent from the Nogales International WWTP. Here, the model underestimated the flow by 10% during the calibration period and overestimated by 7% during the validation period (Table 11). Similar to Nogales, there were inconsistencies in December 2007 (Fig. 4B) at this gauge, where the model matched peak flows, but overestimated high flows and underestimated medium to low flows (Fig. 5B).
The river becomes ephemeral between Tubac and Continental. While the model performed well in the calibration period (Fig. 4C), it underestimated high flows during the validation period by 19% (Table 11). The model mostly underestimated high flows and overestimated low flows (Fig. 5C). Flow characteristics at Tucson (Fig. 4D and Fig 5D) are similar to Continental (Fig. 4C and Fig. 5C), demonstrating that these two gauges are connected during high flows. The model mostly underestimated medium to high flows except for the highest and overestimated medium to low flows (Fig. 5D). Thus while the model performed better during the calibration period (Table 11), it underestimated the flow by 51% during the validation period.
The model performed well at Dodge during both the calibration and validation periods (Fig. 4F). However, there were three cases (August 1992, February 1994 and 1997), where the model predicted high flows in response to precipitation input, contradicting the negligible observed flow (Fig. 4F) and thereby overestimated flow in the calibration period (Table 11). FDCs indicated that the model overestimated all medium and low flows (Fig. 5F).
The river becomes perennial again at Cortaro because of discharge from two nearby wastewater treatment plants (Roger and Ina). The model predicted the flow extremely well from 1996 until 2007 (Fig. 4G). Inconsistencies between observed and simulated flow were remarkable from 1991 to 1995. Inconsistency in the results was attributed to the uncertainty in precipitation input data. This supposition was supported by better validation performance for this gauge (Fig. 5G and Table 11). While the model predicted the peak flow well, it underestimated all other flow conditions (Fig. 5G). The model overestimated the flow by 72% during the calibration period and underestimated the flow by 1% during the validation period (Table 11).
Model Predictions and Possible Sources of uncertainties
The SWAT model was calibrated well at all of the gaging stations based on all model evaluation parameters. However, the model overestimated the total flow during the calibration period at Dodge and Cortaro. It was not possible to correct this using available precipitation data. An auto calibration algorithm (Parasol43) included in the ArcSWAT interface was used with sum of squared errors (SSE) as an objective function to see if its application would improve the model prediction, but results were inferior to manual calibration. The model did validate at all stations, including Dodge and Cortaro, which indicates relative confidence in the model. Performance during validation periods was not as good as the performance in the calibration periods for 5 of the gauge stations (Nogales, Tubac, Continental, Tucson and Vail), which is a common result in watershed modeling studies. While the performance at Dodge during the validation period was as good as during the calibration period, the model performed better in the validation period at Cortaro. At Vail, the model consistently underestimated medium and high flows throughout the validation period. On the other hand, at Cortaro, the model predicted significant flows in several months during the calibration period in contrast to no or very low observed flow. Aside from the uncertainty attributed to the precipitation input, the SCS curve number method, which works on daily rainfall depths, does not consider the duration and intensity of precipitation.17 Representing this precipitation characteristic is necessary for semi-arid watersheds like the Santa Cruz, where high-intensity short duration precipitation occurs.17
Table 11
Model performance statistics.

Figure 5.
FDCs of simulated and observed mean monthly flows at (A) Nogales, (B) Tubac, (C) Continental, (D) Tucson, (E) Vail, (F) Dodge, and (G) Cortaro.

Although no obvious differences in streamflow simulation in large watersheds with different soil data resolution has been found,37 the low-resolution soil data used might not be able to represent the proper soil conditions at the sub-watershed level. Improvement in streamflow prediction using higher-resolution data has been observed;48 however this approach also requires more effort in preparing and calibrating the model with finer resolution. The differences can also be attributed to the model's inability to perfectly represent the antecedent moisture content of the soil before each event. Moreover, rivers in semi-arid regions are more extreme and less predictable than those of rivers in humid regions because of their spatial and temporal flow variations.49,50
At Cortaro, it was noticed that the model generated flow corresponding to significant precipitation events based on input data, but the observed ‘no flow’ indicated that either the precipitation event was not large enough to generate flow, or the model overestimated the antecedent moisture content which resulted in higher flows. This result could be due to underestimation of transmission losses but more likely is due to the uncertainty in precipitation data, since the model responded well in the validation period. Overall, the performance statistics (NSE and R2) and graphs are deemed acceptable42 for an arid watershed.
The model underestimated the low flow conditions in downstream relatively wetter reaches, yet overestimated the low flows in upstream gauges. The data were analyzed after ignoring data corresponding to no flow conditions in the observed data. In dry reaches, the SWAT model generally predicted higher baseflow and in comparatively wet reaches, underestimated baseflow. Although the river remains longitudinally disconnected much of the time, it is connected during big events, as illustrated from time series plots, and the model simulates most of the peak flows throughout the watershed during both the calibration and validation periods.
Findings from this research demonstrate that it is important to determine the parameter sets that hold true for the entire watershed when calibrating at this scale in arid/semi-arid regions, then identify and adjust only a few parameters, based on the observed data at the monitoring points. Since there are only a few months with observed flow in the stream reaches, a longer calibration period is required in arid regions compared to humid regions for better model calibration, otherwise the model can be biased to specific flow conditions and may not work well as for later validation periods, which might be dominated by different flow conditions. Due to high stream bed water losses in semi-arid and arid streams, precipitation at the sub-watersheds closest to the outlet has the most impact on the outlets flow, particularly during low- and medium-storm events.
For initiating the hydrological process at each subwatershed, SWAT assigns precipitation input from the station closest to the centroid of the sub-watershed. Therefore, having only a few precipitation gauges in a large watershed will increase uncertainty in the precipitation input. This result was reflected in many cases in the Santa Cruz watershed, where the model predicted high flow contrary to no or very-low flow in the observed data and vice versa. More precipitation gauges would reduce this uncertainty effect and improve model calibration and validation in large semi-arid/arid watersheds. The flow prediction accuracy would also increase with more precipitation gauge stations spatially distributed throughout the Santa Cruz watershed. Many other studies51–5253 have highlighted the importance of improved precipitation input for obtaining better results using SWAT. The availability of climate data plays a vital role in model performance and accuracy. In particular, the spatial variability of precipitation data represents one of the major limitations in hydrologic modeling of large watersheds.3,41
Effect of Multi-Gauge calibration
The prediction ability of spatially calibrated SWAT model was found to be better compared to the model calibrated at single watershed outlet. With outlet calibration, although the model was able to accurately predict the flow at Dodge, which was relatively closer to the outlet (Fig. 1), the model completely failed to simulate the flow at gauges including Nogales (NSE = -0.54), Tubac (NSE = -0.1), and Tucson (NSE = 0.02), which implied that the model was not representing the entire watershed. The performance of the model improved with multi-gauge calibration when compared to single outlet calibration (Table 12). The improvement was noteworthy at Nogales (NSE = 0.54), Tubac (NSE = 0.41), and Tucson (NSE = 0.52) (Table 12). The performance for the validation period was also improved for the outlet (Table 12). The results from this study suggested that single outlet calibration of the watershed in arid and semi-arid regions can be misleading and thus requires spatial calibration for capturing the spatial heterogeneity and discontinuities in the watershed. This result is particularly important if the developed model, as in this case, is used for assessing the impacts resulting from climate change and urban growth on water resources and identifying sensitive areas for implementing management practices and other water related decision making.
Table 12
Validation performance with only outlet calibration vs. multi-gauge calibration.

Santa Cruz Watershed Ecosystem Portfolio Model(Scwepm)
SCWEPM is an ecosystem services online decisions support tool ( http://lcat.usgs.gov/santacruz/) to assess the impacts of climate change and urban growth, and promote the use of information relevant to water allocation and land management in the Santa Cruz watershed.22 The resulting SWAT model from this study is an integral part of SCWEPM. With further calibration for water quality parameters (sediment and nutrients), the model will be used to examine possible changes in the discharge, evapotranspiration, percolation, surface runoff, transmission losses, water yield, sediment yield and nutrient yield resulting from predicted land use change and climate change scenarios in the Santa Cruz watershed, which will allow bi-national managers to identify problem areas where management activities can be focussed.22 Other impact on ecosystem services being considered include: water quality (nutrient loading), flood control, and erosion potential, which will be visualized and compared using this tool.
Conclusions
This study investigated the effect of multi-gauge calibration for developing a reliable watershed model in the semi-arid Santa Cruz watershed using SWAT. In addition, this paper presents a detailed methodology for developing a watershed model in arid/semi-arid regions. Although SWAT has been applied to many watersheds throughout the world, including watersheds neighboring the Santa Cruz, a detailed calibration, analyses, and recommendations for such an application have not been well documented. This paper is intended to generate a discussion of methodology, and identify standards for calibrating watershed models in arid environments as well as presenting hydrologic results.
Curve number (CN2), soil evaporation and compensation coefficient (ESCO), threshold water depth in shallow aquifer for return flow to occur (GWQMN), Baseflow alpha factor (ALPHA_BF), and effective hydraulic conductivity of the soil layer (CH_K2) were found to be sensitive parameters in calibrating SWAT in semi-arid watersheds like the Santa Cruz. The developed model was calibrated well at all the gauge stations based on all model evaluation parameters and acceptably validated as well. The performance of the model is considered good for an arid watershed, where the response of flow to precipitation events is less predictable than in wetter conditions. The flow prediction accuracy would likely increase with more precipitation gauge stations spatially distributed throughout the Santa Cruz watershed. The performance of the model, as indicated by NSE, much improved with multi-gauge calibration when compared to single outlet calibration at Cortaro. Results from this study accentuate the importance of multi-gauge calibration for developing a reliable watershed model in semi-arid watersheds.
As a next step, the SWAT model will be further calibrated and validated for sediment and nutrient transport. The resulting model is an integral part of the SCWEPM and will be used to examine the tradeoffs, both monetary and non-monetary, of implementing land-management practices on ecosystem services, with consideration for climate change, within the Santa Cruz watershed. Examples of binational ecosystem services being impacted include: water quality (nutrient loading), flood control, and erosion potential.
Author Contributions
Model Calibration/Analysed the data: RN. Wrote the first draft of the manuscript: RN. Contributed to the writing of the manuscript: RN, TM, LMN, JBC. Agree with manuscript results and conclusions: RN, TM, LMN, JBC. Made critical revisions and approved final version: RN, TM, LMN, JBC. All authors reviewed and approved of the final manuscript.
Disclosures and Ethics
As a requirement of publication author(s) have provided to the publisher signed confirmation of compliance with legal and ethical obligations including but not limited to the following: authorship and contributorship, conflicts of interest, privacy and confidentiality and (where applicable) protection of human and animal research subjects. The authors have read and confirmed their agreement with the ICMJE authorship and conflict of interest criteria. The authors have also confirmed that this article is unique and not under consideration or published in any other publication, and that they have permission from rights holders to reproduce any copyrighted material. Any disclosures are made in this section. The external blind peer reviewers report no conflicts of interest.
Acknowledgements
The authors wish to thank the US Geological Survey's (USGS) Geographic Analysis and Monitoring (GAM) Program and US-Mexico Border Environmental Health Initiative (BEHI) as well as the US Environmental Protection Agency's (EPA) Ecosystem Services Research Program (ESRP) and Southwest Ecosystem Services Project (SwESP) for investing in this research. We also want to thank Yongping Yuan and Wenming Nie from the EPA's Landscape Ecology Branch in Las Vegas, NV for their help in acquiring weather and gauge data at the beginning of this study and also Miguel Villarreal for the land use/land cover dataset.