Development of Strategy for SWAT Hydrologic Modeling in Data-Scarce Regions of Peru

The headwaters of the Chili River draining to the El Frayle reservoir in the Arequipa Department of Peru was selected as the study area (Fig. 1). El Frayle is one of the seven water supply reservoirs located upstream of the city of Arequipa, the second-largest city in Peru (the eighth reservoir diverts water for irrigation to the Siguas River). With a total drainage area of $\mathrm{1,026}{\mathrm{km}}^2$, El Frayle has the largest catchment among four reservoirs without upstream regulation (ANA 2014). El Frayle receives about 400 mm in annual precipitation, and almost 85% occurs December through March (Ministerio del Ambiente 2019). Understanding how much water is available and when it is available in headwater catchments will play an important role in water regulation for the entire region because these catchments are major sources of water for public supply, agriculture, and industrial use (Stensrud 2016).

SWAT performance assessment was conducted based on daily hydrologic records for the El Frayle reservoir that are available from the Autonomous Authority of Majes (AUTODEMA) from January 2009 to the present (AUTODEMA 2019). The monitoring data include reservoir water level, effective volume, and discharge, which is determined by representatives from local and regional water authorities. Evaporation, based on a Class A pan located at the El Frayle climate station (Fig. 1), is also provided each day.

SWAT requires physiographic inputs including topography, soil characteristics, land cover, and climate data for model development. For elevation, the 12.5-m terrain-corrected elevation data provided by the Advanced Land Observation Satellite (ALOS)–Phased Array type L-band Synthetic Aperture Radar (PALSAR) were obtained from the National Aeronautics and Space Administration (NASA) Earthdata (2019a). Radiometrically terrain-corrected (RTC) elevation corrects geometric distortions caused by side-looking radar (UAF-ASF 2019). The ALOS-PALSAR-RTC data set was found to provide more accurate elevation in southern Colombia compared to other global elevation data sets by Correa-Muñoz et al. (2019), and has been found useful for river network extraction in semiarid regions (Niipele and Chen 2019), which makes it suitable for this study.

Daily precipitation and temperature data for the five closest weather stations (hereafter called Station_PT) were obtained from the national meteorology and hydrology service of Peru (Ministerio del Ambiente 2019). No stations collect data inside the El Frayle watershed (Fig. 1), or at the higher elevations ($>\mathrm{5,000}\mathrm{m}$). Therefore, a 1-km-resolution gridded, terrain-corrected precipitation and temperature data set developed by Moraes et al. (2019) based on the SENHAMI station data (hereafter called Grid_PT) was also used as an alternative to the station data [Fig. 2(a)]. The average annual precipitation of the El Frayle watershed ranged from 288 to 437 mm. Therefore, subbasins delineated based on elevation with a precipitation range greater than 44 mm (which is 10% of maximum annual average) were manually split into smaller subbasins with lower ranges, resulting in the creation of 106 subbasins (Fig. 2).

SWAT uses climate data from the climate station nearest to a subbasin’s centroid for its hydrological processes (Neitsch et al. 2011). Therefore, using the Station_PT data set (with five stations), the same precipitation and temperature data sets were assigned to several subbasins even with different spatial climate ranges [Fig. 2(b)], while using Grid_PT, precipitation and temperature time series were obtained for pixels associated with the centroids of each subbasin, which resulted in 106 different precipitation and temperature data sets.

Daily streamflow [$Q_{\mathrm{in}}(t)$] was determined using a water balance of the El Frayle reservoir based on daily discharge from the reservoir [$Q_{\mathrm{out}}(t)$], reservoir water level [$H(t)$], and volume [$V(t)$], precipitation [$P(t)$], and pan evaporation [$E(t)$] available from AUTODEMA (Fig. 3) (AUTODEMA 2019). Plotting AUTODEMA records for volume versus water level showed that there is a strong polynomial relationship between volume and water level: {$V_{est}(t)=0.3\times {[H(t)]}^2-\mathrm{2,387.7}\times H(t)+\mathrm{4,757,334}$}. This relationship for estimating volume [$V_{est}(t)$] was used because some volume values appeared erroneous. Then the surface area [$A_{est}(t)$] was estimated by fitting a line to the ratio of observed volume change to the change in water level

A few zeros and extremely high pan evaporation measurements that appeared to be erroneous were also replaced with the average of measured evaporation on a day before and a day after. Finally, the reservoir inflow was estimated from the water balance as

Soil information available for the El Frayle watershed includes maps of soil taxonomy and suitability and a set of soil profiles developed by the Ministry of Environment of Peru (MINAM 2017). Both maps use the same 59 polygons with map unit code as a common polygon field (Table 1). The soil taxonomy map provides order, suborder, group, and subgroup using the US standard soil taxonomy classification (USDA-NRCS 2015). Six dominant taxonomic subgroups (with more than 5% coverage) were present in the watershed (Table 1). The second map describes the capacidad de uso mayor, translated as soil suitability in English. The soil suitability represents the dynamic status of soil agrological quality based on the analysis of edaphic, climatic, and relief characteristics (SENACE 2009). It addresses suitability for cultivation, pasture, forest, or protection based on limitations due to soil, salt, erosion risk, drainage, flood risk, and weather (Table 1). The two soil maps were overlaid to create a map of taxonomy and suitability combinations, referred to as soil classes [Fig. 4(a)]. Soil classes representing less than 2% of the watershed area were merged into the dominant class in the same taxonomic subgroup.

The third type of soil information is a rich set of soil profiles that were described for 160 locations within the Arequipa Department. Profile descriptions include detailed observation and measured properties such as texture and organic matter content that are needed for SWAT simulation, and also taxonomy and suitability classes (Gobierno Regional Arequipa 2016). The challenge was identifying soil profiles that could be representative of each soil class present in the El Frayle watershed because no profile was available for soil classes in the watershed. To bridge this gap, soil profiles were assigned to polygons based on the combination of soil taxonomy and suitability. Soil classes that did not have an observed soil profile were merged with one with the same taxonomy and the most similar suitability code (Table S1).

When more than one soil profile was available for one soil class, the profiles were combined using the slicewise aggregation method (Beaudette et al. 2013). In this method, layers (slices) of a soil class were defined based on horizons of the associated soil profiles, then sand, silt, clay, organic matter, and rock content of layers from all soil profiles were averaged. For a detailed explanation of processes and examples of aggregated profiles or to access the assigned soil properties for the entire department, see Daneshvar et al. (2020b).

The SWAT model requires many soil characteristics (Arnold et al. 2013), some of which (e.g., soil texture, organic matter, and rock fragment contents) could be obtained directly from the soil profile data. Soil properties not directly available were estimated based on soil texture and other available data, using pedotransfer functions and a methodology developed by Daneshvar et al. (2020b). The final result was a set of soil properties for each soil class. Properties for an example soil class (MM-P3se) as well as the range of estimated properties for all soil profiles are given in Table 2. Estimated soil properties for the remaining eight classes (Table 1) are available in Tables S2–S9. The wide spatial variations in estimated soil properties reveal the necessity for regional soil database development because global soil databases such as the FAO/UNESCO soil map of the world (FAO/UNESCO 2003) provided only one soil type for the entire El Frayle watershed.

The Ministry of Environment of Peru provides a regional land cover map and classification at three levels (Gobierno Regional Arequipa 2016). The land uses representing the largest areas in the watershed are pajonals and tolares, which are shrubs or rangeland, while the rest of the watershed includes other shrubs and grasslands (yaretales and césped de puna) in addition to wetlands, water, and bare areas [Fig. 4(b); Table 3]. The SWAT model requires many land cover properties to simulate daily and seasonal changes in evapotranspiration (Arnold et al. 2013). The SWAT land cover database includes properties for US and global land cover classes, but does not include properties specific to local vegetation around the world (Arnold et al. 2013). Daneshvar et al. (2020a) developed a methodology to establish new classes for local Peruvian land cover (Table 3) in two steps. Developing new land cover definitions rather than modifying existing ones as others have done (i.e., Uribe et al. 2013) will facilitate their use in future modeling in the region. In order to develop a regional land cover database for SWAT, the most similar SWAT land cover was chosen as a base (Table 3) by reviewing images and consulting with experts in the region (Daneshvar et al. 2020a). The new land cover class was populated using the SWAT base class parameters prior to adjustments being made to better represent local land use. Once the base class was established, three parameters related to land cover seasonal cycles were adjusted for the local vegetation: ALAI_MIN, BLAI, and T_OPT (Table 3).

Minimum leaf area index (LAI) for plants during dormant period (ALAI_MIN), and maximum leaf area index (BLAI) were estimated based on the remotely sensed leaf area index of pixels representing each local vegetation. These values were extracted from 4-day, 500-m-resolution leaf area index time series from 2010 to 2017 that were obtained from the Moderate Resolution Imaging Spectroradiometer (MODIS) Visible Infrared Imaging Radiometer Suite (VIIRS) subsets (Myneni et al. 2015; ORNL DAAC 2018a). LAI time series of locations covered by a majority of each vegetation type were used to estimate average annual LAI cycle for that type. Finally, ALAI_MIN and BLAI were obtained from the average annual LAI cycle.

The optimal temperature for plant growth (T_OPT) was updated based on the average annual temperature of locations where each vegetation occurred in the watershed (Moraes et al. 2019). For more details on the process and examples of parameters estimation see Daneshvar et al. (2020a).

Estimated ALAI_MIN and BLAI range from 0.1 to $0.2{\mathrm{m}}^2/{\mathrm{m}}^2$ and 0.2 to $1.5{\mathrm{m}}^2/{\mathrm{m}}^2$, respectively, while T_OPT varies from 0.35°C to 3.85°C (Table 3). This shows the importance of regional adjustment because SWAT default values for the original classes (shrubland and grassland) were $0{\mathrm{m}}^2/{\mathrm{m}}^2$, $2–2.5{\mathrm{m}}^2/{\mathrm{m}}^2$, and 25°C, respectively. Using the default values could result in very different values for maximum transpiration rates and seasonal limits for plant growth.

The SWAT simulation was conducted from 2009 to 2017, with the first year of simulation used as warmup and not included in the analysis of results. First, the model was calibrated (from 2010 to 2013) and validated (from 2014 to 2017) for predicted daily streamflow into the reservoir. Then predicted evapotranspiration from water and land were compared with pan evaporation measurement of the reservoir and MODIS measurements for two dominant vegetation types, respectively. Finally, the uncertainty of simulated water balance ratios to estimated inputs were evaluated.

Uncertainty of SWAT outputs were quantified to assess the significance of choices made when developing the model input weather, soil, and land cover data sets. For weather data, outputs of SWAT calibrated with Station_PT and Grid_PT (calibration parameters listed in Table 4) were compared, while for the assessment of soil and land cover data, outputs were used only from SWAT calibrated with Station_PT weather data. Model outputs from the two weather data sources showed that Grid_PT generated 7 mm more streamflow (82 mm instead of 75 mm) as a result of 6 mm more precipitation (402 mm instead of 396 mm), but less surface runoff (2.2 mm instead of 6.6 mm). This resulted in a 67% lower surface runoff/precipitation ratio for Grid_PT relative to Station_PT (Table 6). SWAT simulations were also run swapping the CN_2 values for Station_PT and Grid_PT, and Grid_PT still resulted in lower surface runoff. As a result of the decreased surface runoff, the Grid_PT scenario resulted in slight (less than 8%) increases in base flow, recharge to deep aquifer, and ET (Table 6). Given the minimal differences in annual average precipitation between the weather data products, the differences in surface runoff, base flow, and recharge are likely due to differences in how that precipitation is distributed across the watershed.

Uncertainty assessment of estimated soil properties included the widest range in inputs and resulted in the widest uncertainty range for water balance components (Table 6). Broad ranges of soil texture and organic matter content in soil profiles directly affected percolation rates, which resulted in wide uncertainty ranges for base flow (up to 46%) and recharge to the deep aquifer (up to 36%). Less variation in simulated plant available water resulted in lower uncertainty (up to 16%) for ET estimation. The regression equations that were used to estimate soil hydrologic properties also have uncertainties (Saxton and Rawls 2006), which is not explored here, but would increase the overall uncertainty of model outputs.

Changes in hydrology due to land cover parameterization was assessed by replacing the three estimated plant growth parameters (ALAI_MIN, BLAI, and OPT) with values for each of the four primary land cover types in turn. This resulted in low uncertainty ranges (up to 6%) in estimated water balance ratios (Table 6). The three scenarios based on estimated properties for pajonal, tolares, and yaretales (with $\mathrm{ALAI}\text{\_}\mathrm{MIN}=0.1{\mathrm{m}}^2/{\mathrm{m}}^2$, BLAI ranging from $0.2{\mathrm{m}}^2/{\mathrm{m}}^2$ to $0.75{\mathrm{m}}^2/{\mathrm{m}}^2$, and T_OPT ranging from 0.35°C to 3.85°C) resulted in no change in the water balance components. The fourth scenario based on estimated properties for césped de puna ($\mathrm{ALAI}\text{\_}\mathrm{MIN}=0.2{\mathrm{m}}^2/{\mathrm{m}}^2$, $\mathrm{BLAI}=1.5{\mathrm{m}}^2/{\mathrm{m}}^2$, and $\mathrm{T}\text{\_}\mathrm{OPT}=3.6°\mathrm{C}$) reduced surface runoff by 3% and base flow and recharge to deep aquifer by 6%, while increasing ET by 3%. This suggests that the differences in parameters between pajonal, tolares, and yaretales were not hydrologically significant for this watershed, but the higher values for ALAI_MIN and BLAI for césped de puna were. Results showed that land cover uncertainty had a bigger impact on ET than weather uncertainty, but less on the partitioning of outflow into surface runoff and base flow. Changes in water use by land cover also affected recharge rates, but magnitude was similar to that of changing the weather data. Despite the wide range of surface runoff for two weather data, soil uncertainty had the biggest impact on water balance ratios (Fig. 8). Regional soil sampling and property measurement will significantly reduce these uncertainty ranges; however, the calibrated model provided satisfactory hydrologic predictions proving that developed methodologies for soil and land cover databases can be applied to similar watersheds in data-scarce regions.

The average annual water balance for an 8-year (2010–2017) simulation revealed that El Frayle watershed is a base-flow-driven watershed, where only 8% (2.4%–9.7%) of average annual streamflow (82 mm) is generated from surface runoff [Fig. 9(a)]. This is also in agreement with the observed data and the fact that no precipitation occurred for almost 7 months [Fig. 9(b)], but naturalized observed flow showed that there is a consistent low flow throughout the year [Fig. 5(a)]. Between 57% and 67% (base scenario 62%) of annual precipitation (396 or 402 mm) is removed by ET (227–267 mm; base scenario 246 mm) and only 10% of it reaches streamflow through surface runoff (2–7 mm; base scenario 7 mm) or lateral flow (19–45 mm; base scenario 35 mm). The remaining (95–133 mm; base scenario 104 mm) percolated to groundwater, some of which reaches the stream through return flow, while some percolates to deeper layers.

Results presented here are limited to those from the default simulation (with Station_PT weather data), so users are directed to the results of the uncertainty analysis (Table 6; Fig. 8) presented previously to put these results in context. Surface runoff had the smallest uncertainty range (6.1–7.3 mm), though the low surface runoff, especially that associated with Grid_PT (2.2 mm), served as a constraint. Base flow ranged from 56 to 91 mm, recharge to deep aquifer could vary from 27 to 37 mm, and ET had the widest variation and ranged from 227 to 267 mm.

Precipitation in the watershed is mostly in the form of rainfall (with only about 1% snowfall) that primarily occurs between December and April (with more than 92% of annual rainfall), and the region is mostly dry for the rest of the year with less than 5 mm of precipitation in May through August [Fig. 9(b)]. This wet period corresponds to warmer conditions in the high-altitude regions. Reservoir recharge from streamflow occurs from December to April and typically peaks in February [Fig. 9(b)]. Monthly streamflow averages $0.5\mathrm{mm}/\mathrm{month}$ from the end of the rainy season (in April) until early December. Simulated annual maxima streamflow into the reservoir ranges from $18{\mathrm{m}}^3/\mathrm{s}$ (in 2010) to $65{\mathrm{m}}^3/\mathrm{s}$ (in 2012) (with coefficient of variation of 0.4). On average, surface runoff is generated on only 38 days per year, all of which is associated with heavy precipitation days, where at least one of five stations recorded more than 8 mm of precipitation and average precipitation of the watershed was more than 1 mm. Daily streamflow into the reservoir dropped to less than $0.02{\mathrm{m}}^3/\mathrm{s}$ for 39 days during the dry year of 2014 (230 mm of precipitation), making this an intermittent stream despite a contributing area of more than $\mathrm{1,000}{\mathrm{km}}^2$. The limited duration of the rainy season (December to April) and large variability in monthly maximum flow rates between wet years (e.g., 75 mm in February 2012) and dry years (e.g., 19 mm in February 2010) highlights the significance of reservoir storage and groundwater supply for sustainable water supply in semiarid regions of the Peruvian Andes.

The calibrated and validated model of El Frayle watershed describes a system where infiltration capacity is relatively high compared to typical precipitation intensities, so that over land flow generation is only associated with limited heavy precipitation days of the year (e.g., only 17 days in 2014).

Water moves through the root zone of the silty loam soils relatively rapidly, and due to the high land surface slopes appears to move quickly through near-surface layers to contribute to local streamflow, although we cannot say with certainty whether a perched saturated zone develops or whether it travels as interflow in the unsaturated zone. A physically significant portion (6.9%–9%) of the annual water balance is lost to the local watershed in the form of deep recharge to regional aquifers.

Vegetation is water limited for much of the year, with simulated ET rates well below potential even during the wet season. The model was not sensitive to changes in transpiration potential, as evaluated through changes to LAI, because actual ET is instead controlled by water availability. However, sensitivity to some of the controls that limit transpiration during periods of low humidity was not evaluated and may be important here. In these experiments, ET estimation is most sensitive to the soil parameters that control how much water is held in the root zone as plant available water versus how much percolates through the root zone to the shallow aquifer.

Precipitation product choice is a challenge in this data-sparse region. The low-elevation station network likely underestimates watershed average precipitation by approximately 1.5% (6 mm), but their distribution across the watershed resulted in a 67% (5 mm) increase in surface runoff compared to the Grid_PT data set. Grid_PT predicted more precipitation than Station_PT for lower-elevation river valley areas to the north and east of the reservoir, but less precipitation on steep slopes in the north and south of the watershed (Figs. 1 and 2). The higher precipitation on steeper slopes, in particular, can help explain the higher surface runoff generation of the model that used Station_PT. Precipitation falling on steeper slopes is more likely to result in surface runoff, while precipitation falling on the flatter river valley is more likely to infiltrate. There was minimal temporal variation between the two data sets because Grid_PT was developed based on the Station_PT data set, but because Grid_PT is highly correlated with the elevation, there are some differences in snow versus rain (Fig. S4). Because the model was recalibrated, the choice of precipitation data set did not impact the estimate of available water, merely the flow path, with the estimate of groundwater recharge being the most sensitive.

Groundwater resources are largely undeveloped in this region and there is high uncertainty regarding their yields, but this study simulates an average annual groundwater recharge of $29\mathrm{mm}/\mathrm{year}$ (or $290{\mathrm{m}}^3/\mathrm{ha}$). Annual water balance analysis found that only 20% of annual precipitation resulted in local streamflow. Several recent global assessments have considered a basin to experience water scarcity when total withdrawals for human uses approach 40% of streamflow (e.g., Kummu et al. 2016). For this basin, that implies that only 8% of average annual precipitation is available for allocation. Such thresholds should be considered in water management planning with respect to future climate, urban growth, and irrigation projections. In terms of agricultural production, total annual water yield with respect to environmental water need in the region is $320{\mathrm{m}}^3/\mathrm{ha}$, meaning that a drainage area the size of the El Frayle watershed (102,600 ha, which covers 8% of the Quilca watershed), located in the elevation range of 3,000–5,000 m, can support water needed for up to 5,259 ha of garlic or 14,149 ha of peas (which are two common crops in the region) (Mekonnen and Hoekstra 2011; FAO 2019).