Hydrologic Analysis of an Intensively Irrigated Area in Southern Peru Using a Crop-Field Scale Framework

Majes is one of the largest agricultural areas in the Arequipa region (southern Peru). Low seasonal precipitation and increasing water demands for agricultural irrigation, industry, and human consumption have made water supply projections a major concern. Agricultural development is becoming more extensive in this dry, sunny climate where crops can be grown year-round. However, because this type of project usually involves significant perturbations to the regional water cycle, understanding the effects of irrigation on local hydrology is crucial. Based on the watershedscale Soil and Water Assessment Tool (SWAT), this investigation focuses on the impacts of intensive irrigation on hydrological responses in the Majes region. This study is unique because we allow for crop-field scale input within our regional-scale model to provide information at this smaller scale, which is important to inform local stakeholders and decision makers. Each hydrologic response unit (HRU) was generated to represent an individual crop field, so that management practices could be applied according to real-world scenarios. The management file of each HRU was modified to include different operation schedules for crop rotation, irrigation, harvest, and tillage. The model was calibrated and validated against monthly observed stream discharge during the 2009–2020 period. Additionally, evapotranspiration, irrigation water volume, and daily stream discharge downstream of the local river (Siguas) were used to verify the model performance. A total of 49 sub-basins and 4222 HRUs were created, with 3000 HRUs designated to represent individual crop fields. The simulation results revealed that infiltration from agricultural activities in Majes represents the majority of annual groundwater return flow, which makes a substantial contribution to streamflow downstream of the Siguas River. Simulations also suggested that groundwater flow processes and the interactions between surface and groundwater have a major impact on the water balance of the study area. Additionally, climate variability had a higher impact on surface runoff than on groundwater return flow, illustrating that the groundwater component in the study area is important for future water resources resiliency under expected climate change scenarios. Finally, there is a need to perform a follow-up implementation to provide a guideline for decision-makers to assess future sustainable water resources management under varying climatic conditions for this arid irrigated system.


Introduction
Intensive agricultural irrigation projects bring many benefits to rapidly developing regions worldwide. These include, but are not limited to, providing significant contribution to both local and global food supply, and activating the economy through job creation. Unfortunately, these essential projects can also result in unintended yet significant local issues such as surface and groundwater scarcity through changes in the water balance and natural hydrology, especially in arid and semi-arid regions [1][2][3][4][5]. Depending on factors such as climate and geology, intensive irrigation projects can additionally result in debris flows [6].
Countless studies around the world have focused on the effects of agricultural projects on water resources with a broad range of understanding and conclusions. These studies commonly use hydrological models as the primary evaluating tool to determine the impact of irrigation activities on local water cycles under arid and semi-arid climates [7][8][9][10][11][12][13][14]. Some previous representative rigorous watershed irrigation studies are discussed below to provide context for this work. Liu [15] used a hydrological model at the regional scale with an embedded irrigation module and newly developed annual irrigation maps to quantify the recharge, transpiration, streamflow, and groundwater responses from irrigation in the Republican River Basin (central United States), a watershed that covers about 65,000 km 2 . According to the author's simulations, irrigation affects groundwater recharge throughout the year, while its influence on surface runoff is concentrated in the growing season. In addition, agricultural activities affect recharge and surface runoff more during dry years, and agricultural return flow volumes change depending on the season. In Saudi Arabia, on the other hand, Mahmoud and Alazba [16] evaluated the hydrological responses of land cover changes induced by human activities in an arid region (Al-Baha District) using 10 years of remote-sensing data, concluding that land-use changes such as deforestation and agricultural irrigation strongly affect the availability of water resources. After several years of water shortages, Li et al. [17] simulated the hydrological effects of agricultural irrigation in the Yellow River Basin (China) using 46 years of daily streamflow data, daily precipitation and temperature data from 77 meteorological stations, and data pertaining to cropland and large reservoirs, concluding that the establishment of irrigation projects has been the main variable contributing to streamflow decreases within the 795,000 km 2 watershed. However, to best inform stakeholders and decision makers, model output should be provided in the smallest spatial unit possible, such as the specific crop-field scale; but such smaller scale modeling that can inform both local and regional-scale water budgets have not been conducted, other than for one study in the semi-arid western U.S. [18]. In addition, very few data-based modeling studies addressing high irrigation impacts on local and regional hydrology have been conducted in true arid regions outside of China [17,19], or even for semi-arid regions outside the U.S. and China [3,20], and particularly not for areas with year-round growing seasons.
In contrast to the above, it is additionally believed that a most accurate analyses to evaluate the effects of agricultural irrigation on local hydrology can be accomplished at smaller scales, as demonstrated by Zheng et al. [21], Liu et al. [22], Shi et al. [23], Awan et al. [24], and Wen et al. [25], among many other researchers. Most of those studies used the Soil Water Assessment Tool (SWAT) model, being a universally-accepted mathematical model for representation of reality in watersheds with agricultural activities [18,[26][27][28][29][30][31]. SWAT uses a robust approach on the soil-water balance and allows for relatively complete agricultural management practices (e.g., crop growth, irrigation, fertilization, harvest, and tillage) in irrigated areas. For example, using an event-based irrigation schedule within the SWAT model in the 600 km 2 Barr Creek catchment (Australia), McInerney et al. [32] evaluated the impact of irrigation schedules on evapotranspiration (ET), streamflow, and potential aquifer recharge to improve agricultural productivity and ecosystem sustainability. Additionally, Daneshvar et al. [33] developed a new land-cover database based on local land-cover and satellite data to evaluate the performance of SWAT in the headwaters of the Chili River (southern Peru), with limited input data. Moreover, Ahn et al. [3] applied the SWAT model to quantify changes in horizontal water transfers and vertical water budget under a drought-adaptive management strategy in the semi-arid Rincon Valley (New Mexico, NM, USA), a 2466 km 2 irrigated area. This study considered detailed hydrological components over time and space, and demonstrated the important role of surface water infiltration during the irrigation season. However, in most cases, detailed crop growing/rotation and canal irrigation processes that control hydrological water budget in agricultural regions were not taken into consideration.
SWAT has been approved worldwide as the state-of-the-art-model for agricultural water evaluations. However, its application in arid intensive irrigation arid in South America is still challenging due to the lack of available flow and management data for model calibration and validation. Contributing to the global efforts to understand the smaller-scale effects of intensive irrigation on local water cycles, this study focuses on Majes (Arequipa Region), which is one of the most important agricultural developments in southern Peru. The study area is located on a flat debritic plateau within an extremely arid desert, and the Majes agricultural site in particular is developed above a narrow valley carved by the Siguas River. The Majes irrigation project was established in 1971 to activate the local economy and to supply the region's agricultural needs [34]. However, due to long-term intensive irrigation practices, problems have arisen in this region from an increase in landslide activity along the deeply incised valley [6,35]. Additionally, an advantageous local climate additionally allows crops to be grown year-round. However, concern over future availability of water in the entire Arequipa Region as demands for the resource continue to increase to supply the growing agricultural sector, among other activities [36].
There have been a few local studies focused on water-related effects of agriculture in the area, including the work of Cruz [37], who studied the effects of soil tillage on soil moisture; Parra [38], who evaluated the effects of soil tillage on soil water losses; Sanchez [39], who determined the effects of polyacrylamide on irrigation water availability; Alata [40], who worked on crop conversion models; and Amado [41], who characterized the bacterial population of Majes' irrigation water. However, the sustainability of current agricultural water usage, as well as the effects of irrigation on the area's hydrological behavior and response is unknown. Furthermore, the vast majority of studies done in Majes are largely restricted to either specific governmental regulations or studies aimed at field experiments and theoretical investigations [42,43], which is time consuming, expensive, and cannot be brought to a larger scale [44]. To the best of our knowledge, no mathematical modeling work has been undertaken in this region, which could provide detailed information on hydrological components, contribute to adaptation measures for future land management scenarios that involve sustainable agricultural water use, and offer insights for other similar irrigated areas within this arid region and similar ones elsewhere.
Our study presents a unique modeling framework to evaluate the impacts of agricultural irrigation on water resources in the arid region. Specifically, in an intensively managed agricultural watershed, smaller scale input is needed for accurate hydrologic predictions of agricultural landscapes with diverse crops. In addition, more detailed hydrological water budgets are needed for decision makers that can account for field-level activities. In our approach, individual field management such as detailed fertilizer application, crop growing/rotation, tillage for each crop, and canal irrigation processes are represented at a crop-field scale to make the model more accurate. The smallest spatial unit of SWAT is defined as hydrologic response unit (HRU), which is not an efficient way to discretize physically meaningful boundaries at the field scale. Therefore, we applied an approach to define HRUs by cultivated field boundaries. The ability of the model to capture the relative influence of artificial management features on the system water cycles is a key aspect of this approach.
Considering the above information, the overall objective of this study was to incorporate smaller scale input, more detailed model input, which is required to develop an accurate local-or regional-scale hydrologic model for agricultural lands with numerous and diverse crop types. In addition, the method can allow for more smaller-scale (i.e., fieldscale) detailed information for agricultural water management, although this latter feature will be used in forthcoming publications. This more detailed, smaller-scale framework will enable a full evaluation of the impact of intense irrigation on the local and regional water cycle under extreme aridity in the Majes region and Siguas River watershed, which is a critical region for Peruvian economy. This objective was accomplished via modification and the application of the SWAT modeling approach to capture streamflow, infiltration, and groundwater return flows from local and regional agricultural activities.

Study Area
Majes is one of the largest agricultural areas in southern Peru. However, the region is vulnerable to a number of irrigation-related challenges, as previously mentioned, due to long-term excessive water usage (i.e., little attention is paid to efficiency), deteriorated irrigation infrastructure, and inadequate drainage systems [42,45,46]. This irrigation project is located on a high plateau that lies on the west bank of the Siguas River (one of the main tributaries of the Chili-Vitor River system), which eventually leads into the Pacific Ocean. This river system (Siguas-Chili, Vítor) is the main source of industrial, agricultural, and urban water used in the Arequipa Region ( Figure 1) [43]. The Siguas River watershed is located between longitudes 15 • 44 W and 16 • 38 W, and latitudes 71 • 40 S and 72 • 29 S, with a total drainage area of 3357 km 2 and a total channel length of 380 km. Monthly mean temperature and precipitation from four climate monitoring stations (Figure 1, green dots) are plotted in Figure 2. Temperature ranges within each location is quite stable, but mean temperatures vary due to elevation changes within the drainage area. Annual average temperature in Majes is 18 • C, approximately 8 • C higher than the northern portion of the study region. In Majes, however, daily temperature varies from 9.8 • C to 25.6 • C, with a mean annual temperature of 18.1 • C, which makes agricultural production possible all the year round, as previously mentioned. Climate in the entire study area is considered to be arid to semi-arid, with mean annual rainfall varying from 17 mm (in Majes) to 473 mm (in Madrigal). The majority of annual rainfall occurs between December and April (summer months), primarily during brief but intense convective storms. and diverse crop types. In addition, the method can allow for more smaller-scale (i.e., field-scale) detailed information for agricultural water management, although this latter feature will be used in forthcoming publications. This more detailed, smaller-scale framework will enable a full evaluation of the impact of intense irrigation on the local and regional water cycle under extreme aridity in the Majes region and Siguas River watershed, which is a critical region for Peruvian economy. This objective was accomplished via modification and the application of the SWAT modeling approach to capture streamflow, infiltration, and groundwater return flows from local and regional agricultural activities.

Study Area
Majes is one of the largest agricultural areas in southern Peru. However, the region is vulnerable to a number of irrigation-related challenges, as previously mentioned, due to long-term excessive water usage (i.e., little attention is paid to efficiency), deteriorated irrigation infrastructure, and inadequate drainage systems [42,45,46]. This irrigation project is located on a high plateau that lies on the west bank of the Siguas River (one of the main tributaries of the Chili-Vitor River system), which eventually leads into the Pacific Ocean. This river system (Siguas-Chili, Vítor) is the main source of industrial, agricultural, and urban water used in the Arequipa Region ( Figure 1) [43]. The Siguas River watershed is located between longitudes 15°44′ W and 16°38′ W, and latitudes 71°40′ S and 72°29′ S, with a total drainage area of 3357 km 2 and a total channel length of 380 km. Monthly mean temperature and precipitation from four climate monitoring stations (Figure 1, green dots) are plotted in Figure 2. Temperature ranges within each location is quite stable, but mean temperatures vary due to elevation changes within the drainage area. Annual average temperature in Majes is 18 °C, approximately 8 °C higher than the northern portion of the study region. In Majes, however, daily temperature varies from 9.8 °C to 25.6 °C, with a mean annual temperature of 18.1 °C, which makes agricultural production possible all the year round, as previously mentioned. Climate in the entire study area is considered to be arid to semi-arid, with mean annual rainfall varying from 17 mm (in Majes) to 473 mm (in Madrigal). The majority of annual rainfall occurs between December and April (summer months), primarily during brief but intense convective storms.   Crop water requirements in this region are always in excess of local precipitation (i.e., negative water balance), because this is a desert-type climate with extreme aridity. Therefore, the Majes project started to divert water from the Colca River in 1983 to irrigate an agricultural area of 160 km 2 . The Colca River, located in the neighboring Camaná Basin, is released from the Condoroma Reservoir, then captured and diverted by the "Tuti" diversion structure. From there, the flow is transferred to the Siguas watershed towards an 88 km tunnel and, finally, a 12 km open channel. Eventually, the released water in the Siguas River is captured by the "Pitay" diversion structure (see Figure 1, orange cross) and conducted to Majes through irrigation canals.

Model Description
SWAT is a physically based and spatially semi-distributed model, operating on a daily time step [49]. This model was designed to predict agriculture-related hydrological processes, including surface and soil water flow, nutrient and sediment transport, crop growth, land use, and water management, in complex watersheds under different management scenarios [50][51][52][53]. The watershed is delineated into sub-basins, based on river network and user-defined threshold drainage area, which is further divided into "hydrologic response units" (HRUs) with unique land-use type, soil attributes, and management operations. Simulated results from each HRU are summed in each sub-basin and then routed through river network until the watershed's outlet. SWAT simulations are based on the water balance equation for the land phase of the hydrologic cycle (Equation (1)), where hydrological responses in SWAT are based on the water balance within the soil Crop water requirements in this region are always in excess of local precipitation (i.e., negative water balance), because this is a desert-type climate with extreme aridity. Therefore, the Majes project started to divert water from the Colca River in 1983 to irrigate an agricultural area of 160 km 2 . The Colca River, located in the neighboring Camaná Basin, is released from the Condoroma Reservoir, then captured and diverted by the "Tuti" diversion structure. From there, the flow is transferred to the Siguas watershed towards an 88 km tunnel and, finally, a 12 km open channel. Eventually, the released water in the Siguas River is captured by the "Pitay" diversion structure (see Figure 1, orange cross) and conducted to Majes through irrigation canals.

Model Description
SWAT is a physically based and spatially semi-distributed model, operating on a daily time step [49]. This model was designed to predict agriculture-related hydrological processes, including surface and soil water flow, nutrient and sediment transport, crop growth, land use, and water management, in complex watersheds under different management scenarios [50][51][52][53]. The watershed is delineated into sub-basins, based on river network and user-defined threshold drainage area, which is further divided into "hydrologic response units" (HRUs) with unique land-use type, soil attributes, and management operations. Simulated results from each HRU are summed in each sub-basin and then routed through river network until the watershed's outlet. SWAT simulations are based on the water balance equation for the land phase of the hydrologic cycle (Equation (1)), where hydrological responses in SWAT are based on the water balance within the soil profile, including ET, surface runoff, lateral flow, percolation, and groundwater return flow.
The terms in Equation (1) are defined as total volume of water per unit area (i.e., depth of water), where S t is the final soil moisture content (mm); S o is the initial soil moisture content on day i (mm); P i is the precipitation on day i (mm); ET i is the evapotranspiration on day i (mm); Q i,sur f is the surface runoff on day i (mm); Q i,seep is the percolation entering the vadose zone through the soil profile on day i (mm); Q i,gw is the groundwater return flow on day i (mm); and t is the total time (days).
Surface runoff from daily precipitation was simulated using the Soil Conservation Service curve number method (SCS-CN) [54]. At each time step, CN (curve number) values related to surface runoff were adjusted according to the soil moisture present [49], which were calculated based on land use/land cover, soil type, and management conditions [55]. The process of infiltration water percolated through each soil profile into the root zone and below was represented based on a storage routing technique [56]. The Muskingum routing method [57] was used to route streamflow through the sub-basins, while potential evapotranspiration (PET) was calculated using the Pennman-Monteith approach [56]. SWAT has been implemented worldwide for thousands of watersheds under different environmental conditions [3,9,58,59], as previously mentioned. However, neither SWAT nor other integrated watershed models have been implemented enough in Peru, especially for watersheds with sparse data.

Soil Water Assessment Tool (SWAT) Model Inputs
Major input data for SWAT include digital elevation models (DEM), stream networks, soil properties, spatially variable land-use designations, and historical climatic information including precipitation, maximum and minimum average daily air temperature, solar radiation, relative humidity, and wind speed. The source and applied data are summarized in Table 1. ArcSWAT 2012 was used in this investigation for performing SWAT (Revision 664) simulations [60]. Topographic parameters such as drainage network, number of subbasins, and slope degree and lengths were defined from the DEM using the Advanced Land Observing Satellite (ALOS) Phased Array L-band Synthetic Aperture Radar (PALSAR) Radiometrically Terrain-Corrected 12.5-m pixel size, provided by the Alaska Satellite Facility Distributed Active Archive Data Center (ASF DAAC) (see Figure 3A). A highresolution stream geographic information system (GIS) layer from Arequipa's National Water Authority (Autoridad Nacional del Agua-ANA) was overlaid on the DEM map to improve the hydrographic segmentation and sub-basin boundary delineation. Daily discharge from the Colca River provided by Majes Autonomous Authority (Autoridad Autónoma de Majes-AUTODEMA) was imposed in the Siguas watershed through the positive point source specified near its head (Figure 1, yellow triangle). The diversion volume was provided as a daily varying time series, with an average daily discharge of 12.5 m 3 /s during the past 14 years (2007-2020). A 1:600,000-scale land-cover map obtained from the Peruvian Ministry of the Environment (Ministerio del Ambiente) ( Figure 3B) and a 10 km grid size soil map from the Food and Agricultural Organization (FAO) were applied to create each HRU [61].  Climatic information was obtained from the Peruvian National Service of Meteorology and Hydrology (Servicio Nacional de Meteorología e Hidrología-SENAMHI), which provided 14 years (2007-2020) of weather input data, including daily precipitation and daily air temperature (minimum and maximum) documented from all weather stations shown in Figure 1. As can be seen from the figure, three weather stations in high altitude areas (green dots) are located outside of the study boundary, and only the Majes station is inside the irrigation district. Considering that the majority of streamflow in the Siguas River is diverted from the Colca River (as previously mentioned), precipitation from the air temperature (minimum and maximum) documented from all weather stations shown in Figure 1. As can be seen from the figure, three weather stations in high altitude areas (green dots) are located outside of the study boundary, and only the Majes station is inside the irrigation district. Considering that the majority of streamflow in the Siguas River is diverted from the Colca River (as previously mentioned), precipitation from the upstream region should not have substantial influence on the streamflow and will not further affect the irrigation activities. Moreover, Moraes et al. [62] recently developed a terrain-corrected, precipitation and temperature dataset to show the spatiotemperal distribution of climatic information. The figure in Appendix A shows a comparison conducted between the SWAT-interpolated climate distributions and Moraes's climate model. The pattern of both precipitation and temperature are well represented.
Information related to relative humidity, solar radiation, and wind speed was obtained from the SWAT Global Weather Database, although they were only available until 2014. The SWAT's weather generator tool (a stochastic simulator) [63] was used to determine missing data during the weather definition processes using measured long-term precipitation and temperature data [56]. Furthermore, the streamflow gauging station located along the Siguas River above the city of Majes ( Figure 1, pink star), and the continuous daily irrigation diversion data from the Pitay flow station (Figure 1, orange crossed circle), which were collected from AUTODEMA, were used to evaluate SWAT's simulated streamflow and irrigation diversion, respectively, during the period between January 2009 and February 2020.

Model Configuration
HRU boundaries delineated by SWAT would not typically coincide with the cultivated field boundaries. In this study, a method similar to that developed by Wei et al. [18] was applied to modify the soil layer, resulting in HRUs in line with the boundaries of the cultivated fields, so that crop rotations and management operations within each unit could be assessed. The map of cultivated field boundaries was provided by AUTODEMA and delineated into shapefiles by hand. This method considers each HRU-represented cultivated field as an independent calculation, with only one type of crop growing during a specified period of time. The major steps of this approach were to (1) intersect cultivated field polygons with the soil layer in order to add the boundary of each field into the original soil map; (2) assign a new unique soil ID to each field, with the original soil properties retained; and (3) overlay land use layers, modified soil/cultivated fields layers, and slope in ArcSWAT to generate the HRUs map ( Figure 3). The soil modification should be processed before setting up the model in the ArcSWAT interface, and can be achieved by adding lines with new soil IDs and properties in the user's soil database (the data table belongs to the geodatabase of SWAT2012.mdb). The designated crop cultivation method was later specified in each HRU's management input file. Using this method, the Siguas River watershed was delineated into 49 sub-basins and 4222 HRUs. Out of them, 3000 HRUs are designated to represent cultivated fields.
Parameters of current management practices such as planting date, crop rotations, irrigation and fertilization initiation, harvest, and tillage operations, were defined separately for each cultivated crop and are provided as inputs to the model (through management input files (.mgt)). The planting, harvest and kill, and tillage dates for each crop type is listed in Table  2; these values were obtained from the reports of Huanca [48] and Quiroz [64], as well as from the opinion of the Research Manager in charge of the Universidad Nacional de San Agustin de Arequipa (UNSA)'s Research Station at Majes (Julio Flores, personal communication). Although avocado trees are the dominant fruit species cropped in the area, avocado does not have information in SWAT's crop database. Therefore, it was assumed that avocado trees were already mature and that they kept growing as usual after harvesting each April. The index values of avocado tree's growth, bio-energy ratio, and harvest were adjusted according to the field knowledge provided by Julio Flores, as well as average biomass numbers found in the literature [65,66]. SWAT's auto-fertilization algorithm was initiated at the beginning of each year so that no nutrient stress was generated. Harvest operation is scheduled by using "harvest only" (permanent crops in Table 2) or "harvest and kill" (normal crops) function in SWAT. The "harvest only" function removes plant biomass solely and allows the plant to continue growing, while the "harvest and kill" function removes the plant yield portion, then converts the remaining biomass to residue on the soil surface. If the "harvest and kill" function was applied, the model assumes that the crop has to be replanted before it can start growing. The crop rotation strategy specified for this study was based on local cropping practices, with rotations fixed throughout the simulation period to evaluate the effects on local hydrological components. An example of the SWAT's management input file (.mgt) with two-year pepper/potato rotation practices is provided in Table 3.  The auto-irrigation algorithm was implemented in SWAT, whereby irrigation for each HRU (field) is triggered based on a plant's water demand, which is controlled by the water stress threshold {AUTO_WSTR} (use of this specific bracket denotes the SWAT model input parameter as written in the SWAT's user guide). When plant stress falls below this threshold fraction due to water stress, the model automatically applies the available water from the source to the HRU up to the maximum amount {IRR_MX} specified by the user [56].
In SWAT, the irrigation source can be specified as a river reach, reservoir, shallow aquifer, deep aquifer, or a source outside the watershed. In this study, irrigation water is simulated by a river reach (surface canals) from the Siguas River at the Pitay diversion point shown in Figure 1 (orange cross). The diversion amount is limited by water availability in the river. The maximum amount of irrigation water {IRR_MX} applied for auto-irrigation events was set to 60 mm, based on the maximum daily irrigation diversion divided by the total cultivated field area. User-defined overall water loss from the system can be assigned by the irrigation efficiency parameter {IRR_EFF}, while the surface runoff parameter {IRR_ASQ} was used to implement the surface runoff ratio return back to the system. These values for different irrigation techniques such as sprinkler, drip, or flooding were taken from the FAO's database [67] and specified to vary on each field. The overall procedures to modify and apply the SWAT model to the intensively irrigated plateau is shown in Figure 4. can be assigned by the irrigation efficiency parameter {IRR_EFF}, while the surface runoff parameter {IRR_ASQ} was used to implement the surface runoff ratio return back to the system. These values for different irrigation techniques such as sprinkler, drip, or flooding were taken from the FAO's database [67] and specified to vary on each field. The overall procedures to modify and apply the SWAT model to the intensively irrigated plateau is shown in Figure 4.

SWAT Sensitivity Analysis, Calibration, and Validation
Model sensitivity, calibration, and validation were performed using the SWAT's Calibration and Uncertainty Procedures software (SWAT-CUP) [68]. One-at-a-time sensitivity analysis was implemented manually to identify dominant parameters governing streamflow, with selected parameter ranges determined by literature results from other modelling studies applied to irrigated regions [3,[69][70][71]72]. Then, a global sensitivity analysis using the p-value and t-stats for prioritization was regressed against the objective function of the Nash-Sutcliffe coefficient (NSE), which is the commonly used method to quantify overall performance for integrated watershed climate-hydrology-water quality models [73][74][75].
The simulation period for this study was from January of 2007 to February of 2020, with the first two years considered as a "warm-up" period [76]. Calibration was processed by comparing monthly simulations to observed streamflow at the stream gauge ( Figure  1, pink star) from January 2007 to December 2015. Then, the calibrated model was validated from January 2016 to February 2020. Five hundred model simulations were performed at each iteration for a total of 2500 simulations, using the sensitive parameters listed in Table 4. The new parameter ranges were adjusted after each iteration to obtain the lowest simulated and observed residuals based on the objective function. For further

SWAT Sensitivity Analysis, Calibration, and Validation
Model sensitivity, calibration, and validation were performed using the SWAT's Calibration and Uncertainty Procedures software (SWAT-CUP) [68]. One-at-a-time sensitivity analysis was implemented manually to identify dominant parameters governing streamflow, with selected parameter ranges determined by literature results from other modelling studies applied to irrigated regions [3,[69][70][71][72]. Then, a global sensitivity analysis using the p-value and t-stats for prioritization was regressed against the objective function of the Nash-Sutcliffe coefficient (NSE), which is the commonly used method to quantify overall performance for integrated watershed climate-hydrology-water quality models [73][74][75].
The simulation period for this study was from January of 2007 to February of 2020, with the first two years considered as a "warm-up" period [76]. Calibration was processed by comparing monthly simulations to observed streamflow at the stream gauge (Figure 1, pink star) from January 2007 to December 2015. Then, the calibrated model was validated from January 2016 to February 2020. Five hundred model simulations were performed at each iteration for a total of 2500 simulations, using the sensitive parameters listed in Table 4. The new parameter ranges were adjusted after each iteration to obtain the lowest simulated and observed residuals based on the objective function. For further verification, the monitoring canal diversion record from 2009 to 2020 was converted to irrigation water and compared to the SWAT's simulated irrigation amount, to verify the model's performance. [77] Next, ET values for the Majes irrigation area were verified as well. Large-scale in situ ET measurements are expensive, difficult, and time-consuming [78]. Thus, the Moderate Resolution Imaging Spectroradiometer (MODIS) data carried by the National Aeronautics and Space Administration's (NASA) Terra and Aqua satellites were used as observed data to verify the SWAT's simulated monthly ET outputs (for HRUs). The MODIS 8-day cumulative ET data were averaged to daily ET and re-calculated to be monthly ET values based on the actual observation date, providing a reliable source to analyze crop growing processes within the cultivated area. Spatiotemporal ET maps were downloaded from 2009 to 2020, with a spatial resolution of 500 m. Moreover, a set of streamflow field measurements was collected from the Siguas River (downstream from Majes), between September 2019 and February 2020. Samples were collected from the Siguas River (Figure 1, red box) during 5 sampling campaigns. Difficult site access made automated monitoring and frequent field visits impossible. These datasets were collected to further test model results downstream from the Siguas River. Regarding model performance, the Nash-Sutcliffe Efficiency (NSE), percent bias (PBIAS), and coefficient of determination (R 2 ) were used to quantify the calibration and validation results at the gauging station. These indices were calculated using 86 s (2)-(4), as follows: where Q s is the SWAT's simulated streamflow; Q o is the observed streamflow; Q o is the arithmetic mean of the observed streamflow; and n is the total number of observations. These measurements are widely used in hydrological modeling, thus serving as a benchmark for performance evaluation. NSE ranges from −∞ to 1, with values close to 1 denoting better model performance. The PBIAS values describe whether the simulation results under or overestimate the observed values, ranging between −100 and ∞, with values closer to 0 representing better model performance. Finally, R 2 is an index used to investigate the degree of linear relationship between observed and simulated data, with values varying from 0 to 1, with values close to 1 indicating a good fit.

Model Uncertainty
The model's uncertainty was analyzed using the Sequential Uncertainty Fitting Ver. 2 (SUFI-2) algorithm, which accounts for all sources of uncertainties such as measured data (e.g., observed flowrate), driving variables (e.g., climate data), input parameters, and model structure. The propagations of the sources' uncertainties led to uncertainties in the simulated output variables [79]. SUFI-2 plots all uncertainties on the parameters, expressed as the 95% probability distributions, by using the Latin hypercube sampling method [68]. This algorithm tries to capture most of the measured data within the 95% prediction uncertainty (95PPU) band. The p-factor and r-factor are used to compare the band representing measured data (plus its error) with 95PPU for model simulation. The p-factor is a fraction of measured data (plus its error) bracketed by the 95% prediction uncertainty (95PPU) band, which is calculated at 2.5% and 97.5% levels of the cumulative distribution of an output variable obtained through Latin hypercube sampling [68]; it ranges from 0 to 1, where values approaching 1 indicate better model simulation, considering the uncertainty. The r-factor, on the other hand, is another index reflecting the average thickness of the 95PPU band, with values approaching 0 (i.e., narrower band) being ideal. For reference, p-factors greater than 0.7 and r-factors less than 1.5 were recommended by Abbaspour et al. [79] as satisfactory values.

Parameter Sensitivity Analysis
The global sensitivity results are listed in Table 4 and the final calibrated values are shown in Table 5. Many of these 21 parameters have also been identified as among the most relevant ones for SWAT's calibration by many other studies performed around the world [80][81][82]. Parameters related to surface water dynamics such as SCS-CN for antecedent soil moisture condition {CN2} and the evaporation soil compensation factor {ESCO} demonstrated very high sensitivities, with the former being the most sensitive in this study. The SCS-CN impacts the model simulation of the watershed's runoff, which is a function of land use, antecedent soil moisture conditions, and soil permeability. As shown on Table 5, {CN2} was reduced by 36%, revealing that surface runoff was overestimated by its initial value. It might also indicate that soil properties (e.g., infiltration capacity) were not properly reflected by the initial {CN2} value [80]. The {ESCO} parameter (the fifth ranked parameter) influences baseflow periods and can be used for evaporation demands to account for the effect of capillarity of unsaturated soils [83]. As the {ESCO} value increases, the model defines a lower capacity to retain water in the soil, leading to more lateral flow and total runoff to streams. The irrigation-related parameters were also highly sensitive in the calibration process; irrigation efficiency {IRR_EFF} and the maximum depth of irrigation water {IRR_MX}, which can be applied on the HRU each time irrigation is trigged, were ranked as second and third, respectively. The high sensitivity shown by {IRR_MX} and {IRR_EFF} is consistent with the findings by Ouessar et al. [69] and Ahn et al. [3], who revealed that intense irrigation activities have a substantial impact on the spatial variation of the watershed's hydrology. Similarly, parameters related to groundwater processes also had an important impact on streamflow, being the lag time that it takes for water to leave the bottom of the soil profile and enter into shallow aquifer {GW_DELAY} the most sensitive parameter among them; this reflects the hydrogeological properties of the aquifer system and unsaturated zone, as well as the depth to the groundwater table. Another sensitive groundwater-related parameter was the threshold depth of water in the aquifer specified for return flow to occur ({GWQMN}), which has great impact on groundwater discharge to streams and reflects the sensitivity of the shallow groundwater aquifer properties and groundwater-surface water interactions. The sensitivity associated with these two parameters suggests that, during infiltration from irrigation/precipitation, groundwater makes substantial contribution to streamflow in the Siguas River. Additionally, the fourth ranked parameter (available water capacity of the soil profile, {SOL_AWC}) and the sixth ranked parameter (saturated hydraulic conductivity of the infiltration zone, {SOL_K}) control infiltration water, which were also found to be highly sensitive in other arid to semi-arid watersheds [3,84,85]. The increased values of {SOL_AWC}, {SOL_K}, and decreased value of {CN2}, represent greater soil water storage of infiltrated water, reflecting a reduction of surface runoff.

Model Calibration and Validation
The time series of simulated and observed monthly streamflow comparison for Pitay station is shown in Figure 5, with open area and gray-shaded area indicating the calibration and validation period, respectively. The precipitation shown in Figure 5 is averaged over the entire Siguas River watershed. The streamflow in this figure reflects the general seasonality of the wet and dry season experienced in the headwaters of the Colca River. With 5 iterations, model optimization resulted in an "acceptable" uncertainty, with the modeled p-factor to be 0.98 (calibration) and 0.82 (validation), and r-factor being 0.73 (calibration) and 0.63 (validation), indicating that 98% and 82% of monitored streamflow is enveloped within a really narrow band during the calibration and validation period, respectively. The statistics for the calibration period (NSE = 0.88 and R 2 = 0.89) and validation period (NSE = 0.76 and R 2 = 0.78) indicated a "very good" match between simulated and observed flow in such a way that the model was able to capture the timing and magnitude of streamflow changes [86]. The PBIAS results (PBIAS = −5.3% for calibration and PBIAS = 4.3% for validation, respectively) showed a "satisfactory" performance [86], with the model under-estimating streamflow by 5.3% during the calibration period and over-estimating it by 4.3% during the validation period. Overall, calibration/validation results agree with the reported values documented in the literature [71,87,88].

Model Calibration and Validation
The time series of simulated and observed monthly streamflow comparison for Pitay station is shown in Figure 5, with open area and gray-shaded area indicating the calibration and validation period, respectively. The precipitation shown in Figure 5 is averaged over the entire Siguas River watershed. The streamflow in this figure reflects the general seasonality of the wet and dry season experienced in the headwaters of the Colca River. With 5 iterations, model optimization resulted in an "acceptable" uncertainty, with the modeled p-factor to be 0.98 (calibration) and 0.82 (validation), and r-factor being 0.73 (calibration) and 0.63 (validation), indicating that 98% and 82% of monitored streamflow is enveloped within a really narrow band during the calibration and validation period, respectively. The statistics for the calibration period (NSE = 0.88 and R 2 = 0.89) and validation period (NSE = 0.76 and R 2 = 0.78) indicated a "very good" match between simulated and observed flow in such a way that the model was able to capture the timing and magnitude of streamflow changes [86]. The PBIAS results (PBIAS = −5.3% for calibration and PBIAS = 4.3% for validation, respectively) showed a "satisfactory" performance [86], with the model under-estimating streamflow by 5.3% during the calibration period and over-estimating it by 4.3% during the validation period. Overall, calibration/validation results agree with the reported values documented in the literature [71,87,88]. Results show a good representation of the baseflow during the entire simulation period. However, the model had lower accuracy to predict peak flows. For example, the model overestimated peak flows in 2012 and underestimated the events that occurred between 2015 and 2017. It seems that this issue is typical not only in SWAT simulations [89], but also generally in integrated watershed model calibrations [74,75]. For this study, the authors were more interested in accurate simulation of baseflow and low-flows because these are strongly impacted by irrigation, and analyses of agricultural irrigation as it relates to water management in arid regions was a primary purpose of this research. Given this goal, the model calibration for our intended purpose is actually better than suggested by the overall metrics stated above. Results show a good representation of the baseflow during the entire simulation period. However, the model had lower accuracy to predict peak flows. For example, the model overestimated peak flows in 2012 and underestimated the events that occurred between 2015 and 2017. It seems that this issue is typical not only in SWAT simulations [89], but also generally in integrated watershed model calibrations [74,75]. For this study, the authors were more interested in accurate simulation of baseflow and low-flows because these are strongly impacted by irrigation, and analyses of agricultural irrigation as it relates to water management in arid regions was a primary purpose of this research. Given this goal, the model calibration for our intended purpose is actually better than suggested by the overall metrics stated above.
The relatively poor simulation of peak flows likely happened during intense storms, primarily during the summer months. The lack of a sufficient number and well-distributed climatic gauging stations within the watershed is possibly one reason for under or overestimating peak flows, since rainfall recording locations in the model were attributed to an area much larger than that actually affected by the events, and inadequate meteorological input data cannot completely represent, for example, El Niño and La Niña impacts on local precipitation. However, the monitoring system adopted by Peru's "Autoridad Nacional del Agua" (National Water Authority, ANA) are not yet sufficient to compose a complete database for agro-hydrological studies.

Model Verification
The structure at the Pitay station (Figure 1, orange cross) diverted the majority of Siguas streamflow to irrigate croplands in the Majes project, which is a significant part of the water balance and an important system to sustain irrigated agricultural in the Siguas River watershed. Therefore, verifying the accuracy of model-simulated diversions at the Siguas River is important. The diversion volume associated with agricultural irrigation can be controlled via the following parameters: (1) heat units and (2) irrigation efficiency factor ({IRR_EFF}, i.e., the overall water lost including evaporative loss and conveyance loss). The calibration results show that 25% water losses in transit from the source of supply to the pointed HRU's soil profile. Thus, 75% of diversion was used for crop irrigation during the simulation period. Based on this approach, the monthly verification comparison between calculated irrigation water and irrigated volume simulated by SWAT (considering 14 years) is shown in Figure 6A. As can be seen from the figure, irrigation volumes are very high in summer seasons and low during winter months. This result is explained by significantly less evapotranspiration resulting from reduced incoming solar radiation during winter seasons, together with lower average daily temperatures in that same period (see Figure 2). Following the evaluation criteria defined by Moriasi et al. [86], comparison results indicated a "good" performance for R 2 (0.78) and PBIAS (5%), while an NSE value of 0.51 can be considered as "satisfactory". Model results indicated that the average irrigation volume simulated by SWAT (2.83 × 10 7 m 3 /month) was similar to the average irrigation volume calculated from observed records (2.88 × 10 7 m 3 /month). Therefore, the calibrated SWAT model can provide accurate diversion rates and irrigation amounts for the study region.
In addition, matching the pattern of irrigation return flow from hundreds of crop fields could also affect the model's ability to simulate overall hydrological responses within the watershed. The irregular and random management operations carried out locally by different farmers are difficult to account for, leading to unpredictable water yields from precipitation and irrigation events. This statement is also reported by several studies that showed a discrepancy between simulated and observed streamflow in irrigated regions [27,90]. Moreover, the one-dimensional analytical lagging routine of SWAT model is not able to provide physical representation of aquifer heterogeneity properties and regional groundwater hydraulic gradients, which can finally cause the mistiming of groundwater return flows to the river. This issue is relevant for most integrated watershed models, which typically have relatively poor representation of complex aquifers systems [72,91,92], but could be improved in future investigations by linking the SWAT model with a more rigorous groundwater hydrology model such as MODFLOW (modular three-dimensional finite-difference groundwater model), to represent groundwater flow processes in spatiallyvarying aquifer systems, as was performed by Bailey et al. [93] and Wei and Bailey [72]. However, implementation of a more complex groundwater model is extremely challenging in arid regions of developing countries, which typically have little or no available aquifer data [94].
Monthly comparison of ET between MODIS and SWAT was performed only for the Majes irrigated area, from 2009 to 2020 ( Figure 6B). The graph not only indicated the monthly variation of ET, but also showed the yearly ET pattern. As can be seen from the figure, the overall dynamics and spread of the simulated and estimated ET match each other well. The statistical execution of ET results fits a "good" performance criteria, with an NSE value of 0.76 and an R 2 value of 0.77. PBIAS can be categorized as "very good", with the obtained 3% value [86]. In general, the mean value of estimated ET from MODIS was 58.3 mm/month, slightly less than the 59.6 mm/month simulated by SWAT, with the highest ET occurring in February (summer) and the lowest in July (winter). Additionally, ET estimated by SWAT was higher in some winter periods and lower in some summer seasons, compared to the MODIS values. One explanation for this difference is that the ET estimated by MODIS has some uncertainties because the meteorological properties are not unique within a 500 m MODIS pixel to induce ET variations for certain periods [95,96]. ET depths varied from 36.5 mm to 101.8 mm during wet seasons (summers), whereas during dry months (no rainfall) it ranged from 29.2 mm to 31.2 mm. Moreover, as shown in Figure 6A,B, variations of irrigation volume and ET have similar patterns, since evapotranspiration in agricultural areas is affected not only by precipitation, but also it determines how often irrigation events occur. These results could provide useful information to identify seasonal events and improve understanding of the interactions between agricultural activities and climate variability. In addition, matching the pattern of irrigation return flow from hundreds of crop fields could also affect the model's ability to simulate overall hydrological responses within the watershed. The irregular and random management operations carried out locally by different farmers are difficult to account for, leading to unpredictable water yields from precipitation and irrigation events. This statement is also reported by several studies that showed a discrepancy between simulated and observed streamflow in irrigated regions [27,90]. Moreover, the one-dimensional analytical lagging routine of SWAT model is not able to provide physical representation of aquifer heterogeneity properties and regional groundwater hydraulic gradients, which can finally cause the mistiming of Considering the extremely high flows, limited access and the brief but intense storms, it was unsafe to collect peak-flow data during spring and summer seasons. Additionally, COVID-19 restrictions further prevented streamflow monitoring between March and November of 2020. Therefore, only five observed values were available during the low-flow periods (from September 2019 to February 2020), for the estimated daily flow downstream from Majes (Figure 1, red box). The daily time series plot is shown in Figure 6C, with the black dot line representing the simulated values and red dots representing the observed streamflow. SWAT simulated results varied from 0.3 m 3 /s to 39 m 3 /s. The model was able to estimate the magnitude of streamflow during the low-flow period, comparing them to observed data, highlighting the 3 February 2020 event, which caught the recession phase of the hydrograph, in addition to peak and baseflow. Groundwater has substantial impact on the water balance of the study region, so this particular result that matches the low and minimum flows is important for the overall goals of this study, to evaluate the impact of irrigation on groundwater return flows. A better understanding of the daily streamflow behavior and monthly or seasonal trends would improve future water quality predictions and better inform agricultural best management practices, but a longer time of observation data would be needed.

Water Balance Analysis by Agricultural Irrigations
Using the calibrated modeling results, a water-balance analysis was conducted from 2009 to 2020 to evaluate the relative flow components in the Siguas River watershed. SWAT results showed a negative average balance of −1755 mm/year between water input and PET, suggesting that the system is under a strong water deficit (i.e., the differences between annual water gains and losses are 1755 mm, being the latter a lot larger than the former). Hydrological deficits of 850.2 mm/year were also reported by Rivas-Tabares et al. [82] in a study performed in a semi-arid irrigation region located in Spain. For the entire study area of the present investigation, the annual average depth of water inputs (556 mm) was consumed by the following components: ET (60%), percolation (31%), total runoff at the watershed outlet (34%), and re-evaporation from the shallow aquifer (3%). Percolation (i.e., the seepage water vertically flowing out of the soil profile) is controlled by soil's water holding capacity, infiltration rate, and evapotranspiration uptake. In particular for the total runoff, 4% and 8% of the water input was converted to surface and lateral flow, respectively, while 22% was represented by groundwater flow. This large percent of percolation and groundwater flow can reveal that the river water dynamics are strongly affected by the groundwater aquifer's hydraulic properties [97]. These findings indicate that the surface-aquifer-stream interactions are important for water dynamics in the study region.
Annual averaged hydrological elements (mm) comparisons between the irrigated area in Majes and the non-irrigated area of the Siguas River watershed are shown in Figure 7. Results indicate that the annual precipitation occurring upstream of the study area is higher than that at Majes. Considering such difference, water requirements in Majes were supplemented with water diverted outside from the Colca watershed [45]. Land use changes for agricultural activities in the irrigated area produced 49% more annual ET, due to the availability of moisture in the soil. As expected, surface runoff values in non-irrigated areas were found to be in the lower range because a large portion of the terrain is covered by herbaceous plants with higher roughness, which can slow down surface runoff and allow more infiltration into shallow soils [98]. The downward percolation is further minimized by horizontal subsurface heterogeneities, which could cause lateral discharge to the cliff's walls. Moreover, average soil percolation and groundwater discharge values at Majes were 148.7 mm and 126.4 mm, respectively, significantly higher than the non-irrigated area due to the increased water volumes applied to the surface in the irrigated areas. and allow more infiltration into shallow soils [98]. The downward percolation is further minimized by horizontal subsurface heterogeneities, which could cause lateral discharge to the cliff's walls. Moreover, average soil percolation and groundwater discharge values at Majes were 148.7 mm and 126.4 mm, respectively, significantly higher than the nonirrigated area due to the increased water volumes applied to the surface in the irrigated areas. Soil type is another important driver that affects the local hydrologic cycle, especially soil processes and groundwater flows. The predominant soil type in the irrigated area is loam, which has higher moisture holding capacity, and has infiltration rates high enough to drain irrigation water into the aquifer system, as compared to other less-permeable soil types [99]. As a result, a large fraction of irrigation is infiltrated into groundwater through soil percolation, suggesting that irrigation return flow is a significant feature of water availability. Additionally, the processes of percolation and aquifer recharge affected by infiltration can increase perched water table levels and minimize particle cohesivity within the soil profile due to saturation [100]. It can further lead to slope instability on the cliff wall and result in landslides along the edge of the cliff in this intensively irrigated arid region [101]. As a consequence, geophysical studies are recommended to further investigate shallow and deeper subsurface heterogeneities and porous materials to enable a better understanding of aquifer heterogeneity and groundwater flow processes. Figure 8 illustrates the distribution of monthly water yield results (i.e., surface flow, lateral flow, and groundwater flow) during the 14-year simulation period in the irrigated Soil type is another important driver that affects the local hydrologic cycle, especially soil processes and groundwater flows. The predominant soil type in the irrigated area is loam, which has higher moisture holding capacity, and has infiltration rates high enough to drain irrigation water into the aquifer system, as compared to other less-permeable soil types [99]. As a result, a large fraction of irrigation is infiltrated into groundwater through soil percolation, suggesting that irrigation return flow is a significant feature of water availability. Additionally, the processes of percolation and aquifer recharge affected by infiltration can increase perched water table levels and minimize particle cohesivity within the soil profile due to saturation [100]. It can further lead to slope instability on the cliff wall and result in landslides along the edge of the cliff in this intensively irrigated arid region [101]. As a consequence, geophysical studies are recommended to further investigate shallow and deeper subsurface heterogeneities and porous materials to enable a better understanding of aquifer heterogeneity and groundwater flow processes. Figure 8 illustrates the distribution of monthly water yield results (i.e., surface flow, lateral flow, and groundwater flow) during the 14-year simulation period in the irrigated area. In general, most of the water entering the Siguas River comes from the aquifer (87%), followed by surface runoff (11%) and lateral flow (2%), which means that groundwater makes the main contribution to the streamflow downstream from Majes. As can be seen from Figure 8, simulated monthly surface runoff had a small seasonal periodicity. Monthly discharges are characterized for having higher depths (1.56-1.6 mm) between January and March (summer season) and lower depths (1.13-1.17 mm) between June and August (winter season), which follows the natural annual oscillations (Figure 3). Regarding monthly lateral flows, the median discharges maintained relatively low values with wide ranges in the first four months, an indication of a high inter-annual variability during summer seasons, but becoming more stable during the rest of the year. However, monthly groundwater flows showed higher average depths from 11.8 mm to 12.1 mm between May and August, with a 4-month lag time for the high-flow season, compared to surface flow, which goes directly into the river. The lag time implies a delay between water percolating out of the soil and groundwater flowing to the stream, due to the longer traveling time of return flows traveling through the vadose zone and the aquifer system, until reaching the river [15]. Moreover, the steady upper and lower quartile ranges of monthly groundwater flow data revealed a low inter-annual variability. Irrigation is applied all-year-round and, therefore, the impact of climate variation on surface runoff is more substantial than on groundwater return flows. This is important because it suggests that the groundwater component is relevant for future water resources resiliency under expected but uncertain future climate change. 21, 13, x FOR PEER REVIEW 20 of 29

Limitations of Data and Methods
The developed SWAT model can be a useful tool for crop planning and decision making to secure water resources for future generations. However, there was no stream gauge located at the outlet of the watershed to completely evaluate the model's behavior of the entire hydrologic cycle, leading to a large uncertainty of model's performance in the irrigation area. As a result, this study provided additional components to verify hydrological variables. However, the acceptable match between observed and simulated irrigation water can help confirm that the automatic irrigation algorithm trigged the proper amount of irrigation water for crops; the comparison between MODIS-and SWAT-simulated ET data provided confidence about the plant's growing processes; and the manually monitored streamflow during low-flow periods downstream from the study site indicated that the applied model was able to capture the recession phase of the hydrograph, which is affected by subsurface water responses. These good matches between simulated and observed data improved the reliability of the model. Additionally, as mentioned above, due

Limitations of Data and Methods
The developed SWAT model can be a useful tool for crop planning and decision making to secure water resources for future generations. However, there was no stream gauge located at the outlet of the watershed to completely evaluate the model's behavior of the entire hydrologic cycle, leading to a large uncertainty of model's performance in the irrigation area. As a result, this study provided additional components to verify hydrological variables. However, the acceptable match between observed and simulated irrigation water can help confirm that the automatic irrigation algorithm trigged the proper amount of irrigation water for crops; the comparison between MODIS-and SWATsimulated ET data provided confidence about the plant's growing processes; and the manually monitored streamflow during low-flow periods downstream from the study site indicated that the applied model was able to capture the recession phase of the hydrograph, which is affected by subsurface water responses. These good matches between simulated and observed data improved the reliability of the model. Additionally, as mentioned above, due to the difficult access and safety consideration, no peak-flow samples were collected during the 2019 and 2020 summer seasons. Moreover, coronavirus disease 2019 (COVID-19) restrictions further prohibited additional field samples after February of 2020. Therefore, the streamflow sampling campaigns were performed during lowflow periods, which offered limited values in validating daily streamflow simulation. Although peak-flow points reflect direct rainfall-runoff processes and low-flow points tend to show groundwater return flows from irrigation events (the overall goal of this study), a longer time of observation data will be needed to better inform agricultural best management practices.

Management Recommendations for Future Applications
A number of studies have successfully illustrated that the spatially explicit SWAT hydrological model could be widely used as a practical tool to estimate the impact of management strategies prior to execution, in order to provide science-based evaluations for policy-and decision-making [102,103]. One of the most common applications of SWAT as a decision-making tool is to use optimal irrigation technologies. Currently, irrigation is mainly applied using traditional flood systems, which can be easily replaced by more efficient drip irrigation that is characterized by having lower evaporation losses [104]. According to the simulation results, system ET, soil percolation, and groundwater return flow are typically directly related to irrigation activities at Majes. Drip irrigation technologies maintain higher levels of moisture in the root zone, meeting the crop's growth requirements with lower volumes of water stored in the soil surface which, in turn, is able to slow down soil moisture evaporation and reduce the rate of seepage processes. However, the benefit of flood irrigation is that it can be used to help manage salt build-up in the soil, while this kind of salt amelioration would not be achieved by drip irrigation [105]. The application of mulch (i.e., a thick layer of organic material that covers the soil surface, in order to decrease evaporation rates) is another strategy to increase irrigation efficiency [106]. Following the same principle, decreasing the intensity of tillage and leaving residue in place result in more residues remaining on the soil surface from a previous crop, which can help to reduce surface runoff and soil evaporation while increasing the efficiency of irrigation. Excessive soil percolation can also be minimized by the improvement of conveyance efficiency (e.g., the use of water retaining polymers and/or biochar to minimize seepage volumes, evaporation rates, and available water for plants to use [104]). Polymers as a soil ameliorant may be more efficient at water retention and, therefore, more appropriate at Majes, while biochar brings the added benefit soil health while contributing to the mitigation of climate change through carbon sequestration. In this regard, explicit applications of these management alternatives have a key role to increase irrigation efficiency; thus, the amount of water applied to irrigated fields can be decreased, which subsequently lead to the reduction of groundwater recharge from deep percolation, decreasing also return flows and, as previously mentioned, mitigating the potential for agricultural pollution loading into the Siguas River. Moreover, pollutants such as heavy metals, which can be naturally found in the region's waters, can be treated with low-cost local organic materials, as concluded by Garcia-Chevesich et al. [107], resulting in healthier agricultural products. Similarly, less seepage water from efficient irrigation practices results in reduced risk of agricultural-induced landslides [6], a problem already present in the region [35].
In addition, surface water resources are fed by glaciers and precipitation occurring in the high Andes mountains in southern Peru. Climate change projections by general circulation models (GCM) strongly support the scenario in which spatial and temporal precipitation and temperature patterns will be shifting, resulting in changes in the timing of precipitation, reduced snowpack, and more evaporation in Peru [108]. Similarly, a large number of glaciers in the Peruvian Andes have already disappeared and water supply is expected to decrease in the coming years, which may eventually affect hydrological processes and exert pressure on farmers [109]. Thus, future climate projections should be part of water management planning in this dry region. For example, the increase in temperature could lead to soil moisture deficits, increase irrigation requirement, and generate risks of crop reductions, due to raising evapotranspiration and decreasing soil moisture within cultivated areas. As a consequence, strategic and operational aspects should be combined for the development of sustainable water management policies that consider climate change [110]. The SWAT model developed in this investigation can be easily linked to regional climate models (under different temperature and water availability scenarios), together with information on current and projected agricultural practices and regional crop types to evaluate future water extraction for irrigation projects, hydrological budget, and system sustainability, in order to determine whether it will be possible to keep current cultivation patterns in Majes or whether significant agricultural management changes will be necessary to avoid water scarcity and crop losses.
Due to the increasing consumption of irrigation water in Majes, on the other hand, Siguas River streamflow downstream from Majes has experienced significant decreases in recent decades, with its baseflow always below 0.3 m 3 /s during the 2019 dry season (see Figure 6C). Maintaining relatively high river flow levels is a critical element in monitoring river conditions to support sustainability of the aquatic ecosystem. Similarly, as cultivated areas increase every year, more water will need to be diverted from the Colca River located in the adjacent watershed, in order to fulfill reductions in water supply, which will affect the ecohydrology of the Colca River Basin as well. In order to reduce the degradation of riparian ecosystems in both rivers, a variety of management alternatives can be applied to decrease diversion water and maintain high river flows. A promising solution is an efficiency strategy that deliberately use less irrigation water but still aims to satisfy the crop yield at certain times in the growing period [111]; for example, applying less water during early vegetative and grain-filling stages that are less sensitive to drought stress. Crops such as cotton, maize, potato, sunflower, and sugar beet can provide feasible and acceptable irrigation options for minimal yield reductions with limited supplies of water. In order to ensure successful deficit of irrigation in the study region, it is necessary to consider the soil's capacity to retain water and understand crop yield responses to water stress. With the assistance of our model, management practice modifications such as adopting flexible planting dates, applying less water for different crops, decreasing plant population, and dynamic uses of land can all be evaluated. Additionally, the applied model's results suggest the resiliency of groundwater under the impact of climate variation. Applying groundwater pumping as a supplementary source instead of surface water diversion for irrigation could also contribute to irrigation demands. Advancements such as conjunctive use of surface and groundwater sources would significantly enhance the flexibility and applicability of irrigation water to meet agricultural needs and maintain higher flows in the Siguas River. However, a groundwater quality analysis need to be conducted due to the possibility of high concentrations of salinity, metals, and nutrients in the groundwater [2]. Management operations must fully account for surface and groundwater interactions, which need further efforts to link the current SWAT model with the physically based groundwater model MODFLOW to address the role of groundwater irrigation in the hydrological cycle due to the heterogeneity of the aquifer system.
Considering all of the above, this is a water-scarce region. In other words, unless efficient irrigation practices are strongly applied and sustainable importation projects introduced that transport water from other sources such as neighboring aquifers or even the coast (desal), the entire region is destined to suffer the consequences. Luckily, local agencies and the Universidad Nacional de San Agustín de Arequipa are already working on how to ensure the sustainable use of water resources in the long term, as the region's demands continue to increase every year. In the meantime, using models to simulate environmental processes from agriculture at the Siguas River watershed can enhance local heterogeneities in terms of climate, land processes, crop growth, and management to be taken into account in order to discover cost-effective combinations of measures that could facilitate future decision-making. The main benefit of our model is the fact that each cultivated crop type has been considered as an HRU; thus, a particular agricultural management strategy can be defined independently for each crop type and can be transferred locally with fewer political and practical difficulties. On the other hand, our model is able to handle the implementation of a small number of strategically allocated best management practices on the targeted crops, in order to quantify the interactions among such practices (for each crop type), improving research capacity on how water management influences the hydrologic cycle of the cultivated area. However, to implement strategies identified with this approach, high operational complexity and large efforts are required, because the strategies need to be changed at the HRU level (i.e., crop by crop).

Summary and Conclusions
The purpose of this study was to set up a physically-based, spatially-distributed SWAT model to evaluate the influence of irrigation activities on agricultural water management, contributing with valuable insights for a better understanding of the local water cycle responses from the intense irrigation project in the important Majes region and in the Siguas River watershed, a semi-arid to arid plateau located in Southern Peru. In this study, an approach was applied to define HRUs by cultivated field boundaries, and individual field management such as detailed irrigation activities, management operations for each crop, and canal irrigation processes are represented at the field scale to best inform the stakeholders and decision makers. The model was used to characterize hydrological processes, including a detailed description of surface runoff, evapotranspiration, soil percolation, lateral flow, and groundwater return flow, which helped to investigate the spatial variability of hydrological fluxes during different seasons. The model outputs were tested against observed stream discharge, calculated irrigation water, and MODIS estimated evapotranspiration.
Results indicated that the calibrated model is a useful tool to provide satisfactory performance to represent natural flows and agricultural management operations, describing hydrological responses and improving our understanding of the interactions between agricultural activities and climate variability in irrigated areas located in dry regions. In general, most of the water entering the Siguas River (downstream from the irrigated area) comes from Majes. By comparing the water budget between the irrigated and non-irrigated areas, current irrigation practices increase groundwater recharge and return flow in the Majes region. In addition, the impact of climate variation on surface flow is more substantial than on the groundwater return flow, illustrating that the groundwater component in this area is important for future water resources resiliency under the uncertainty of current climate trends. Due to the aridity of the area, the effect of agricultural activities on the system's water balance depends on complex local factors including soil type, precipitation, and irrigation efficiency.
The crop-field scale results provided by this study can be used to perform a follow-up implementation to investigate the effects of management alternatives on water quantity from irrigated agricultural areas to the receiving water bodies, ensuring that such management occurs optimally in terms of quantity, time, and location, providing options for policy makers to regulate water consumption and sustainable usage of the resource. Additionally, the spatial-temporal analysis can be extended to predict other modelling parameters such as water quality, but a longer time of observation data and additional water quality processes-related modules should be included. This model can further be used to predict how changes in precipitation and temperature are likely to affect water yields and soil water dynamics in the next decades, with the ultimate goal of exploring options to mitigate the influence of climate change and human consumption. One limitation of SWAT, however, is that the model was unable to represent satisfactory streamflow patterns, especially the observed peak flows. In this study area, groundwater has a substantial impact on the system's hydrological budget and, therefore, it should be explored and improved in future investigations to obtain better estimations of groundwater flow processes in heterogeneous aquifer systems by linking the SWAT model to the groundwater model MODFLOW.
Finally, a clear conclusion is that, based on the results obtained, irrigation efficiency in Majes is far from appropriate. Moreover, significant changes to minimize seepage volumes must be applied urgently, thus reducing the risk of future water shortages and geological hazards.