Water Resources Evaluation and Sustainability Considering Climate Change and Future Anthropic Demands in the Arequipa Region of Southern Peru

: Climate change and increases in human activities are threatening water availability in the Arequipa Region (southern Peru). However, to date, there has not been a comprehensive inventory of surface water data or an investigation of current surface water conditions or forecasted future conditions resulting from increased anthropic demand or stresses from climate change. This study evaluates surface water resources management including storage, diversions, and conveyance in the Arequipa Region, while creating a tool for the evaluation of future scenarios in the ﬁve main watersheds of this arid region of southern Peru. State-of-the art, open-source modeling software was used. Water uses for each watershed were evaluated against predicted reservoir inﬂows and streamﬂows for different periods. In addition to the above, 12 climate change models and different shared socioeconomic pathways (SSP) were ensembled for the ﬁve watersheds. A semi-distributed approach and an innovative simulation splitting approach was used for each watershed, which allowed for different starting dates for the simulations using all available data obtained from different sources (government and private). Results indicate that the region is expected to have increased ﬂows during the wet season and no signiﬁcant changes during the dry season. Reservoir inﬂows are expected to increase up to 42 and 216% for the lowest and highest SSP evaluated, respectively. Similarly, the model projected streamﬂow increases up to 295 and 704%, respectively. Regarding yearly water availability and considering current and future demands for the watersheds under study, water deﬁcits are not expected in the future if current reservoir storage can be maintained, though it is expected that reservoirs won’t be able to store predicted higher ﬂows, so important volumes of water could be lost during the wet season to the ocean by natural drainage. Given the uncertainty of climate change projections, if future water sustainability is desired, storage and irrigation efﬁciencies should be improved and reservoir sedimentation should be evaluated.


Introduction
Peru is a country highly threatened by climate change as a result of its geography, inequalities, and natural resource-dependent economy [1,2], with the main effects being floods, droughts, landslides, and glacial melting; the severity of these effects is expected to be amplified by El Niño Southern Oscillation (ENSO) [3,4].Studies have been focused on different Peruvian watersheds to evaluate what would be the possible local effects of global warming.However, there is not a single answer for climate change effects on the entire Peruvian territory since results vary by location, final objective, and climate change models being used.For example, the Santa, Mantaro, and Rimac rivers were assessed and a reduction in peak flows was calculated for those three important Peruvian rivers [5].Moreover, some rivers in the Amazonian regions of Peru are expected to increase their streamflow [6], while the Lake Titicaca basin will have more intense, frequent, and prolonged droughts [7].However, the complex human development dynamics make it harder to predict future water availability under unpredictable climatic changes [8].While climate change analysis is usually focused on evaluating precipitation and temperature trends (e.g., [9]), authors around the world are now including human activities when predicting water resources availability (e.g., [10][11][12]).
Considering the above, no significant research has been performed in terms of water resources sustainability involving both climate change and human activities in the Arequipa Region of southern Peru (the focus of this study), besides a few studies in which hydrological models were applied, but without a focus on future climate and water demands.Mendoza et al. [13], for example, used a Soil and Water Assessment Tool (SWAT) model in the Tambo watershed (excluding the drainage area for Pasto Grande, the watershed's main reservoir), while Coloma [14] applied a combination of hydrologic and hydraulic modeling, focusing on the lowest stream gage in that same watershed, using Hydrologic Engineering Center-Hydrologic Modeling System (HEC-HMS) and HEC's River Analysis System (HEC-RAS) models to analyze return periods.Similarly, the United States Agency for International Development (USAID) [15] applied a Water Evaluation and Planning System (WEAP) model to the Quilca watershed, including most of the individual water uses and "naturalizing" all of the reservoirs in the area.Also, in the Quilca watershed, Wei et al. [16] developed a SWAT model to assess the connection between intensive irrigation and groundwater return flows to an adjacent river in the Majes district, concluding that most of the applied irrigated water is lost through percolation, while Daneshvar et al. [17] applied a SWAT model in the El Frayle reservoir (Quilca watershed) using the streamflow naturalization approach.However, while water consumption in the Arequipa Region continues to increase, the combined effects of both climate change and human activities on water resources availability is currently unknown.Furthermore, water governance in the Region is deficient and not well known, and while several government and private entities collect hydrological information, such information is spaced and hard to analyze, while the diversions, storage, and transportation of water in the Region are unclear, resulting in a weakness for sustainable water management.
This study was designed to (1) characterize for the first time the future flows and reservoir operations in the five main watersheds of the Arequipa Region; and (2) evaluate surface water resources sustainability in the Region considering both climate change and human activities.Our study differs from what has been carried out before in the Region since it starts with a detailed description of how surface water is managed in each watershed (with the purpose of building the appropriate models), followed by accurate watershed delineation based on available data, model setup and calibration, and ending with climate change scenarios for multiple models and SSPs, in addition to water consumption projections for each study site, all to evaluate sustainability of surface water resources.The analysis was performed for each of the five watersheds and it evaluates changes based on all of the observed available data, and future climate and use projections for each drainage area.Furthermore, yearly water availability and water deficits were evaluated.Special, atypical emphasis was placed on the reservoir network as it regulates the Region's flows throughout the year.Finally, considering the novelty of this investigation, this manuscript was written in more detail than the average modeling publication, with the intention of having it as a guide for other Peruvian scientists and organizations to replicate the analysis in other watersheds within the country.

Study Area and Methods
The extension of the watersheds of interest and the political boundaries of the Arequipa Region are illustrated in Figure 1, all beginning in the Andes Mountains and draining into the Pacific Ocean to the southwest.Geographic, climatic, and population characteristics for each watershed are listed in Table 1.The climate in those watersheds varies with altitude, with annual precipitation ranging from 0 mm near the coast to 800 mm in the high Andes.The wet season is concentrated during summer months (December through March) [18].Figure 1 also shows the location of climatic and fluviometric stations, reservoirs, infrastructure flow gages, and water diversions that were used to build the hydrological models.Different activities are carried out for the socioeconomic development of the population that live in those watersheds, allowing the generation of income through agriculture, fishery, commerce, mining, and tourism.Finally, to facilitate an easy replication of this study in other areas of Peru, a complete and detailed explanation of the procedures used to develop the hydrological models is found in Supplemental Information (SI) 1.

Results and Discussion
Considering the nature of this investigation, results and their discussion are displayed combining both water management characteristics and hydrological modeling.A discussion of important results is included in each Section, for clarity and better comprehension of the reader.Additionally, even though some paragraphs herein might correspond to Methods, they were included in this Section since they were a result from the investigation (procedures applied after performing literature review, data gathering, and  Due to the nature and complexity of this study, the methods described include the hydrological analysis of the five main watersheds within the Arequipa Region of southern Peru, though small portions of other administrative Departments (Moquegua, Puno, Cusco, Apurimac, and Ayacucho) were also considered since they were included within the watersheds' limits (see Figure 1).All of the tools applied were open source, which allows further collaboration and replicability among any interested party [29]; important since an analysis such as this has never been performed in Peru and it is highly needed in many regions of the country.
Based on available information and communications with local government and private water management entities, each watershed was characterized in terms of where and when water was diverted, stored, and transported (within and between watersheds), all to ensure proper model development.The chosen model was Precipitation-Runoff Modeling System (PRMS), developed by the United States Geological Survey (USGS) (see details in [30]).PRMS was selected since: (1) it has an open source code, making the analysis easier to be replicated in Peru; (2) it is continuously being updated; (3) It has a variety of modules for computing different variables and distributing precipitation and temperatures; (4) it allows the activation/deactivation of Hydrologic Response Units (HRU) of the simulations; (5) it can be easily scripted, which improves its customization; (6) it can simulate all kinds of model complexity; and (7) it can be easily integrated with more advanced codes such as the US Geological Survey modular finite-difference flow model (MODFLOW).
Finally, to facilitate an easy replication of this study in other areas of Peru, a complete and detailed explanation of the procedures used to develop the hydrological models is found in Supplemental Information (SI) 1.

Results and Discussion
Considering the nature of this investigation, results and their discussion are displayed combining both water management characteristics and hydrological modeling.A discussion of important results is included in each Section, for clarity and better comprehension of the reader.Additionally, even though some paragraphs herein might correspond to Methods, they were included in this Section since they were a result from the investigation (procedures applied after performing literature review, data gathering, and surface water management characterization).

Data Gathering and Management
After a remarkable effort, obtained hydroclimatological information needed to run the models was gathered from government and private organizations and organized in a hydroshare folder [31], to make it available to the general public.For hydraulic infrastructure, the main datasets collected were reservoir volumes, precipitation, minimum and maximum temperature, and evaporation rates.Reservoirs were considered in three out of five watersheds (three reservoirs in the Camana, five in the Quilca, and one in the Tambo watersheds).For the Quilca watershed, the Uzuña reservoir was not included for three main reasons: (1) the storage volume is negligible; (2) it was built in the year 2015; and (3) it is disconnected from the main water supply and hydroelectric energy generation systems.Similarly, diversions were considered in two out of five watersheds, with two main diversions in the Camana and Quilca watersheds; both of them start in the Camana watershed, meaning that there is a big dependence of the Quilca watershed on Camana.The first diversion allows water to come from two upstream reservoirs (Pañe and Bamputañe reservoirs) to the top of the Quilca watershed, while the second one was built for diverting flows from the middle of the Camana basin to Quilca; its function is to supply water for the Majes irrigation district, which is located at the edge of the Camana and Quilca watersheds.Details about the reservoirs can be found in Table 2 and a plane view of all of the data gathered is illustrated in Figure 1 (Section 2).The datasets were available starting from the year 1981, but the simulations had different periods based on available stream and reservoir data.All of them had an ending date of 31 December 2016.The Yauca and Ocoña watersheds did not change their starting periods since they only contain streams.The Camana and Quilca watersheds, however, started as soon as 1 January 2006, which is the first date available in the datasets provided by Arequipa's Electricity Generating Enterprise (EGASA).Priority was given to the EGASA (private) datasets instead of the Majes Autonomus Authority (AUTODEMA, government) ones, since the latter had incomplete records, which included missing values, different measurements, and missing variables such as spillways for some reservoirs.

Reservoir Characterization
Since the study area is hydrologically unknown, it is important to explain some important details for each reservoir: -There is only one reservoir (Pasto Grande) in the Tambo watershed, which has two outlets.The first outlet is a channel that diverts the water to the Moquegua watershed, which is located outside of the basins of interest.The second outlet diverts the water to the Tambo watershed, which receives the spillway flows and an annual volume of 8.2 Hm 3 between September and December [33].-There are four reservoirs in the Camana watershed, with three of them regulating water supply for the upper Quilca watershed through the "Pañe-Colca" channel, while one regulates water to be diverted downstream to the middle Quilca watershed for irrigation purposes (Majes Irrigation District).
The Pañe reservoir has natural inflows and two regulated outflows.The first outflow drifts the water toward the "Pañe-Colca" channel, while the second one diverts the spillway to the Camana watershed through the natural stream network.
The Bamputañe reservoir has natural inflows and contributes all its flows to the "Pañe-Colca" channel.The Dique reservoir has three main inflows: (1) the first one is through the "Pañe-Colca" channel, which brings the water from the Pañe and Bamputañe outflows; (2) the second one corresponds to a water intake in the "Antasalla" channel, and (3) the third one has natural inflows coming from the upstream rivers of the Camana watershed.Moreover, it has two outlets, one that diverts the water to the "Zamacola" channel (continuation of the "Pañe-Colca" channel that diverts the water towards the Quilca watershed), and one that routes the spillways to the Camana watershed through what would be the natural routing without infrastructure.The Condorma reservoir is the biggest reservoir concerning storage in the Camana watershed.It has a natural inflow which is the contribution of its natural drainage area, the spillways from the Pañe reservoir, and the spillways of the Dique reservoir.There is only one outflow for scheduled releases and spillways, and routes the outflow for where it would be the natural routing without infrastructure.
-For the Quilca watershed, there are also four reservoirs, which have a similar distribution compared to the Camana watershed; three reservoirs are located upstream, while the one downstream gathers the flows of the other three.
The Challhuanca reservoir has one natural inflow and only one outflow from where spillway and regulated flows are routed.
The Pillones reservoir has two inflows and one outflow.The outflow is through its dam and the inflows come from a tunnel diversion that benefits from the flows coming through the "Zamacola" channel and natural drainage from its basin and the Sumbay one.
The "El Frayle" reservoir is similar to the Challhuanca setup, but with a bigger drainage area.
The "Aguada Blanca" reservoir receives the outflows from its drainage area, the residual flows left downstream of the Pillones tunnel, and the three aforementioned reservoirs.Regarding its outflows, it has one for regulated outflow and spillway through its dam, and one inner tunnel that diverts flows to a hydroelectric plant (Charcani V), which is located not so far downstream of it.
Regarding the Quilca watershed, there is a diversion not so far downstream where the "Zamacola" channel diverts into the Quilca basin, which diverts water into the "Pillones" reservoir.The problem with this diversion is that the remaining downstream flows are not measured and, as a result, the "Aguada Blanca" reservoir calibration will be affected since there is an unmeasured regulated flow (see Figure 2).
For model buildup, the most complicated setup for the reservoirs was the one that connects the "El Pañe", "Bamputañe", "Dique los Españoles", and "Pillones" reservoirs (Camana and Quilca watersheds), mainly due to the fact that there are channels in-between and different inflows and outflows from them. Figure 2 shows a detailed distribution and flow directions for these reservoirs, not only highlighting where data inflows and outflows are located, but also the ones without records.There are many not-mentioned unregulated diversions throughout the extension of the "Pañe-Colca" channel that affected the inflow modeling for the Condoroma reservoir.The same happens for the "Dique los Españoles" reservoir, which has inflows coming from the "Antasalla" and the "Pañe-Colca" channels.In the case of the "Antasalla" channel, all of the inflows coming through this channel and the remaining flows of its upstream river are unknown, which affects the mass balance for the "Dique los Españoles" reservoir and the inflows going into the "Condoroma" reservoir.Regarding the "Pañe-Colca" channel diversion into the "Dique los Españoles" reservoir, there is water being derived into the reservoir by the "Jancolacaya" water intake, which only measures the flows continuing through the channel and not what is going into the reservoir.Luckily, this problem can be taken as a local problem for the "Dique los Españoles" reservoir since flows are measured downstream of the "Pañe-Colca" channel when it changes its name to "Zamacola" channel, which is the channel that takes the water into the Quilca watershed.The effect of all of these unknowns in the inflow volumes for the "Condoroma" reservoir could be negligible since it has the highest storage capacity from its watershed.
Sustainability 2023, 15, x FOR PEER REVIEW 7 of 32 in the inflow volumes for the "Condoroma" reservoir could be negligible since it has the highest storage capacity from its watershed.

Watershed Delineation and Land Use and Soil Classification
The watershed and land-use delineation required processing of data, so this task is presented in Results, rather than in Methods.The main variables required for the model were elevation, soil type, coverage type, and land use for each HRU, which were defined from the land use data geometries.Soil type was not considered for delineating the HRUs since the datasets found had numerous uncertainties and there were no related comprehensive studies.An example of the distribution of some input variables for simulating one of the watersheds is illustrated in Figure 3.

Watershed Delineation and Land Use and Soil Classification
The watershed and land-use delineation required processing of data, so this task is presented in Results, rather than in Methods.The main variables required for the model were elevation, soil type, coverage type, and land use for each HRU, which were defined from the land use data geometries.Soil type was not considered for delineating the HRUs since the datasets found had numerous uncertainties and there were no related comprehensive studies.An example of the distribution of some input variables for simulating one of the watersheds is illustrated in Figure 3.The delineation of basins and subbasins considered drainage points, where stream gages, diversions, or reservoirs were present.Most of the watersheds had missing areas in the lowest parts of the basins since there are no gages in those places except for the Camana and Ocoña watersheds, which have main cities surrounded by agriculture in those areas.There were at least nine subbasins in each watershed (Figure 4), where the basins with "Cuenca" or "basin" as a prefix in the legend are subbasins found in the National Water Authority (ANA) repositories.
It is important to mention that the way the Quilca watershed was modeled included two drainage points since downstream sub-basins without stream gage data available were not considered in the simulation, being the Pitay watershed the smallest (green-colored watershed in Figure 4).Moreover, the number of HRUs varied among watersheds and the amount has a direct relationship with the running time of each simulation, where Camana required longest run time (see Table 3).The delineation of basins and subbasins considered drainage points, where stream gages, diversions, or reservoirs were present.Most of the watersheds had missing areas in the lowest parts of the basins since there are no gages in those places except for the Camana and Ocoña watersheds, which have main cities surrounded by agriculture in those areas.There were at least nine subbasins in each watershed (Figure 4), where the basins with "Cuenca" or "basin" as a prefix in the legend are subbasins found in the National Water Authority (ANA) repositories.
ANA delineations, and areas for the five watersheds.

Model Set-Up
For all five watersheds, simulations were performed based on available data (model file descriptions are provided in SI 2, which includes instructions for model files).For the case of watersheds with reservoirs, a multi-basin approach was performed, while for watersheds without them, only one simulation was needed.The number of subbasins per simulation varied from one (when there was just a stream gage) to multiple subbasins.In the case of reservoirs, there were at least two subbasins, from which one corresponded to the reservoir geometry by itself.This approach has the advantage of reduced running times during calibration, as HRUs of subbasins can be inactivated if needed.
1.When analyzing the simulation setup, it is important to differentiate between drainage areas dominated by streamflows, diversions, and reservoirs (see Figure 5).In the case of drainage areas with reservoirs, a simulation was needed for each, though PRMS allows the simulation of up to four reservoirs at once.One subbasin was needed for each diversion since water transfers were implemented to the entire subbasin and taken out of the resulting subbasin outflow.For the case of streamflow, all could be simulated at once and there was no need to split the simulation based on the number of stream gages, as illustrated in Figure 5.It is important to mention that the way the Quilca watershed was modeled included two drainage points since downstream sub-basins without stream gage data available were not considered in the simulation, being the Pitay watershed the smallest (green-colored watershed in Figure 4).Moreover, the number of HRUs varied among watersheds and the amount has a direct relationship with the running time of each simulation, where Camana required longest run time (see Table 3).

Model Set-Up
For all five watersheds, simulations were performed based on available data (model file descriptions are provided in SI 2, which includes instructions for model files).For the case of watersheds with reservoirs, a multi-basin approach was performed, while for watersheds without them, only one simulation was needed.The number of subbasins per simulation varied from one (when there was just a stream gage) to multiple subbasins.In the case of reservoirs, there were at least two subbasins, from which one corresponded to the reservoir geometry by itself.This approach has the advantage of reduced running times during calibration, as HRUs of subbasins can be inactivated if needed.

1.
When analyzing the simulation setup, it is important to differentiate between drainage areas dominated by streamflows, diversions, and reservoirs (see Figure 5).In the case of drainage areas with reservoirs, a simulation was needed for each, though PRMS allows the simulation of up to four reservoirs at once.One subbasin was needed for each diversion since water transfers were implemented to the entire subbasin and taken out of the resulting subbasin outflow.For the case of streamflow, all could be simulated at once and there was no need to split the simulation based on the number of stream gages, as illustrated in Figure 5.
the reservoir geometry by itself.This approach has the advantage of reduced running times during calibration, as HRUs of subbasins can be inactivated if needed.
1.When analyzing the simulation setup, it is important to differentiate between drainage areas dominated by streamflows, diversions, and reservoirs (see Figure 5).In the case of drainage areas with reservoirs, a simulation was needed for each, though PRMS allows the simulation of up to four reservoirs at once.One subbasin was needed for each diversion since water transfers were implemented to the entire subbasin and taken out of the resulting subbasin outflow.For the case of streamflow, all could be simulated at once and there was no need to split the simulation based on the number of stream gages, as illustrated in Figure 5. 2. To speed up the simulations, climate files for each one of the five watersheds were created.Three main variables are dominant in the inverse distance-elevation weighting (IDE) method of PRMS, from which two main variables were adjusted.
The first variable ("dist_exp") is the most important since it provides more weight in the spatial distribution for a certain weather station, while the second one ("prcp_wght_dist") is a ratio to adjust more weight to distance or elevation, from which 75% percent of weight was given to distance and 25% to elevation.The problem with elevation was the drastic changes in the characteristic topography of the Peruvian Andes, which caused a lot of "noise" in the distribution.To test the first variable, a gridded HRU approach was applied to smooth the distribution of the climate variables.The gridded approach allows seeing a proper distribution of the climate variables compared to the semi-distributed one that is used for the models.The "dist_exp" variable is highly sensitive (i.e., the higher its value, the more influence

2.
To speed up the simulations, climate files for each one of the five watersheds were created.Three main variables are dominant in the inverse distance-elevation weighting (IDE) method of PRMS, from which two main variables were adjusted.The first variable ("dist_exp") is the most important since it provides more weight in the spatial distribution for a certain weather station, while the second one ("prcp_wght_dist") is a ratio to adjust more weight to distance or elevation, from which 75% percent of weight was given to distance and 25% to elevation.The problem with elevation was the drastic changes in the characteristic topography of the Peruvian Andes, which caused a lot of "noise" in the distribution.To test the first variable, a gridded HRU approach was applied to smooth the distribution of the climate variables.The gridded approach allows seeing a proper distribution of the climate variables compared to the semi-distributed one that is used for the models.The "dist_exp" variable is highly sensitive (i.e., the higher its value, the more influence the weather stations will have in nearby cells).Many values were tested for this parameter and a value of 3.5 was chosen since it smooths the distributions best (see Figure 6).The same values were applied to the five watersheds and Climate by HRU (CBH) files were created.
Considering surface outflows (and not reservoir calibration), in the Camana watershed the only area with a different simulation starting date than listed above was the Bamputañe reservoir, which was built in 2009.Therefore, there are three years in which the Bamputañe reservoir would be providing flows to the Condoroma reservoir but are not included in the simulation, which is estimated to be negligible since the storage in the latter is significantly higher.In the case of the Quilca watershed, the starting period of the Aguada Blanca reservoir (located downstream of the other three reservoirs) is affected by the construction year of the upstream-located reservoirs Pillones (in 2007) and Challhuanca (in 2009), but this does not affect the downstream gages from it since there were data available for its outflows since 2006.The Tambo watershed was the one that had the longest simulation period, while the Yauca watershed had only three years of data (gathered by an irrigation association), so the distribution of the data was poor.Furthermore, the Ocoña watershed also had a bad overage for the simulated periods.the weather stations will have in nearby cells).Many values were tested for this parameter and a value of 3.5 was chosen since it smooths the distributions best (see Figure 6).The same values were applied to the five watersheds and Climate by HRU (CBH) files were created.Considering surface outflows (and not reservoir calibration), in the Camana watershed the only area with a different simulation starting date than listed above was the Bamputañe reservoir, which was built in 2009.Therefore, there are three years in which the Bamputañe reservoir would be providing flows to the Condoroma reservoir but are not included in the simulation, which is estimated to be negligible since the storage in the latter is significantly higher.In the case of the Quilca watershed, the starting period of the Aguada Blanca reservoir (located downstream of the other three reservoirs) is affected by the construction year of the upstream-located reservoirs Pillones (in 2007) and Challhuanca (in 2009), but this does not affect the downstream gages from it since there were data available for its outflows since 2006.The Tambo watershed was the one that had the longest simulation period, while the Yauca watershed had only three years of data Before mentioning the calibration results, it is relevant to define which data were most important for certain watersheds, based on our literature search and information gathering.For instance, since the Yauca watershed has only one stream gage, this gage is extremely relevant for this watershed.The situation is different for the other watersheds; it depends on whether there are reservoirs present or just stream gages.In the case of streamflows, the flows vary depending on the season (dry or wet), so both seasons are analyzed, defining a dry period from May to November, and a wet period from December to April.In the case of reservoirs' surface outflows, only scheduled releases were evaluated.Additional relevant information for the Ocoña, Camana, Quilca, and Tambo watersheds is described below: - The Ocoña watershed has only two stream gages.The Ocoña stream gage is more important since it is located in the lowest elevation, representing a larger portion of the watershed.The Salamanca stream gage contributes around 11 and 4.5% for the dry and wet seasons, respectively.-For the Camana watershed, the Pañe, Bamputañe, and Dique reservoirs only contribute their spillway flows to the Camana watershed, which is negligible for a given year since it only happens sometimes during the wet season.Condoroma outflows, on the other hand, are highly tied to streamflow measurements in Tuti; Condoroma outflows are 74 and 22% of the Tuti flow for the dry and wet seasons, respectively.Since the Tuti stream gage diverts a big portion of its flows to the Quilca watershed, it contributes more water to the Huatipa station during the wet season (50%) than during the dry one (2%).Additionally, there is a flow loss of around 29% from the Huatipa to the Hacienda Pampata stream gages during the dry season, which could due to a shortage of data in the Hacienda Pampata time series.Huatipa, on the other hand, contributes 73% of Hacienda Pampata flows during the wet season.- In the Quilca watershed, some external inputs come altogether through the Zamacola channel.The Pañe, Bamputañe, and Dique los Españoles reservoirs contribute 41%, 18%, and 10.35%, respectively.It is important to acknowledge that there are unknown flows from some intakes along the Pañe-Colca channel.The Zamacola channel contributes to the Sumbay stream gage with 82% and 69% for dry and wet seasons, respectively.Moving downstream to the Aguada Blanca reservoir, 53%, 6%, 12%, and 23% of the regular discharges from the reservoir are contributed by the Sumbay, Challhuanca, Pillones, and Frayle diversion/reservoirs, respectively.The last two stream gages in this watershed are the Tingo Grande and Puente del Diablo; it is hard to compare the flows downstream on Puente del Diablo with the Aguada Blanca releases since there are some water uses associated with this stream.There is a flow loss of 30% during the dry season, while for the wet season a one-to-one relationship was found.Meanwhile, for the Tingo Grande flow recording station, there is a contribution of 2% and 8% for the dry and wet seasons, respectively.-For the Tambo watershed, 1% and 3% are contributed for the dry and wet seasons, respectively, between the Pasto Grande reservoir and the Santa Rosa stream gage.
From the calculated relationships, it can be seen that: (1) for the Ocoña watershed, the Salamanca stream gage does not have an important contribution and the main gage is the Ocoña station; (2) for the Camana watershed, the three upstream reservoirs are only contributing spillway flows to the Condoroma reservoir, and it is mainly used to supply flows for the Tuti diversion.Moreover, a big gap between the Tuti and Huatipa stations was detected, from where the Condoroma reservoir and the Tuti diversion characterize the overall diversion/reservoir infrastructure of the regulated part of this watershed, and the Huatipa and Hacienda Pampata stations represent the unregulated part of the watershed, being the Huatipa station the most important since it catches the flows of the drainage areas that supply most of the water; and (3) for the Quilca watershed, the Pañe reservoir from the Camana watershed supplies most of the water through the Zamacola channel.Similarly, this channel supplies most of the water to the Aguada Blanca reservoir compared to the other three reservoirs that this watershed has (El Frayle, Pillones, and Challhuanca).The Aguada Blanca reservoir is the most important one (in terms of surface outflows) since it regulates the Arequipa city water supply.Moreover, the Puente del Diablo station is located downstream from this reservoir and almost at the southernmost point of the city.For this station, the most important contribution is given by the Aguada Blanca reservoir outflows.

Model Calibration
The main runoff-related variable adjusted was "carea_max" (see SI 3), which specifies the percentage of the HRU that produces runoff.The other main variables represent the region between the soil zone and groundwater, which were the most complex.The soil zone included a capillary reservoir (amount of water held in soil moisture), gravity (mobile in the soil zone), and preferential flow (PFR) storage (not included).Hence, after getting familiar with each watershed and how water is managed on them, a conceptual map (see Figure 7) for these storages was prepared to guide the calibration.
The main runoff-related variable adjusted was "carea_max" (see SI 3), which specifies the percentage of the HRU that produces runoff.The other main variables represent the region between the soil zone and groundwater, which were the most complex.The soil zone included a capillary reservoir (amount of water held in soil moisture), gravity (mobile in the soil zone), and preferential flow (PFR) storage (not included).Hence, after getting familiar with each watershed and how water is managed on them, a conceptual map (see Figure 7) for these storages was prepared to guide the calibration.The sensitivity analysis varied depending on the target observation and subbasin.Evaporation rates were not analyzed for sensitivity since one variable was available to adjust.Moreover, a critical factor for the sensitivity results was the boundaries applied to each sensitivity analysis.Overall, the model output was most sensitive groundwater parameters.Figure 8 shows an example of a sensitivity analysis where the metric to be evaluated was the total sensitivity index, with the model output being most sensitive to "gwflow_coef" and "gwstor_min" (the first parameter adjusts the baseflow rates and the second one adjusts the minimum groundwater storage to be preserved in the subbasin).The sensitivity analysis varied depending on the target observation and subbasin.Evaporation rates were not analyzed for sensitivity since one variable was available to adjust.Moreover, a critical factor for the sensitivity results was the boundaries applied to each sensitivity analysis.Overall, the model output was most sensitive groundwater parameters.Figure 8 shows an example of a sensitivity analysis where the metric to be evaluated was the total sensitivity index, with the model output being most sensitive to "gwflow_coef" and "gwstor_min" (the first parameter adjusts the baseflow rates and the second one adjusts the minimum groundwater storage to be preserved in the subbasin).The time series for each observed dataset was calibrated.Figures 9-11 represent examples of calibrated time series for evaporation, reservoir volumes, and streamflow, respectively.Reservoir volumes and streamflows are in English units since it is the default model output for PRMS, as previously mentioned.Moreover, evaporation calibrations got Nash-Sutcliffe Efficiency (NSE) greater than one and very low Root Mean Squared Error (RMSE) since the calibration was performed based on monthly means and there was only a need to adjust one evaporation parameter for each month.A complete view of all calibrated datasets for reservoir evaporation, reservoir volumes, and streamflows are shown in SI 4, SI 5, and SI 6, respectively.
RMSE in the calibrations has many magnitudes and it changes between simulations depending on the minimum and maximum observed records.This statistic was helpful to reduce errors when NSE could not be improved.Starting dates of simulations, the percentage of data available, and the calibration statistics can be seen in Table 4.While most sub-basin calibrations had great fits, it was not possible to calibrate some of them since there was occasionally missing or insufficient data.Furthermore, for some reservoirs such as Challhuanca and Aguada Blanca, the spillways were not used since they dried the reservoirs when simulated.The Tambo reservoir did not have a precipitation gage, which means that the calibration only relied on the Peruvian Interpolated data of the SENAMHI's Climatological and hydrological Observations (PISCO) dataset.Three  While most sub-basin calibrations had great fits, it was not possible to calibrate some of them since there was occasionally missing or insufficient data.Furthermore, for some reservoirs such as Challhuanca and Aguada Blanca, the spillways were not used since they dried the reservoirs when simulated.The Tambo reservoir did not have a precipitation gage, which means that the calibration only relied on the Peruvian Interpolated data of the SENAMHI's Climatological and hydrological Observations (PISCO) dataset.Three RMSE in the calibrations has many magnitudes and it changes between simulations depending on the minimum and maximum observed records.This statistic was helpful to reduce errors when NSE could not be improved.Starting dates of simulations, the percentage of data available, and the calibration statistics can be seen in Table 4.While most sub-basin calibrations had great fits, it was not possible to calibrate some of them since there was occasionally missing or insufficient data.Furthermore, for some reservoirs such as Challhuanca and Aguada Blanca, the spillways were not used since they dried the reservoirs when simulated.The Tambo reservoir did not have a precipitation gage, which means that the calibration only relied on the Peruvian Interpolated data of the SENAMHI's Climatological and hydrological Observations (PISCO) dataset.Three calibrations were considered unsuccessful (detailed below): - The Salamanca subbasin, located in the Ocoña watershed, from which an NSE of 0.35 was obtained, though only a small percentage of importance was associated with the main stream gage.There was not enough data to catch the dry seasons and, more importantly, different baseflow magnitudes during the dry seasons.The hypothesis associated with this is that this stream gage is located downstream of a glacier (see Figure 12), which by having different melting periods is changing this dynamic.Moreover, some previous studies (e.g., [34]) have tried to calibrate this stream gage but were not successful and lesser performance (NSE = 0.20).A study developed by ANA had better results but the simulation was developed on a monthly basis [22].-For the Dique los Españoles reservoir, an NSE of 0.15 was achieved.The problem faced in this reservoir was missing inflows from the hydraulic infrastructure surrounding the structure, and probably high groundwater flows since there is a nearby lagoon.Also, there are missing records from the Jancolacoya and Antasalla water intakes (see Figure 2).From all of the available General Circulation Model (GCM) datasets, 13 models selected for both historic (1940-2014) and future (2015-2100) datasets (see Section 3. details).Three main errors were found in these downscaled GCMs: (1) the minimum peratures were greater than the maximum ones for certain points near the shoreline o study area, which is something that also happened with the PISCO dataset (see SI 1 - The Aguada Blanca reservoir had the worst calibration statistics.Poor calibration was since the downstream flows coming from the Pillones Tunnel were not recorded, so mass balance errors were obtained when trying to do a mass balance between the flows of the Sumbay stream gage (which is located not so far upstream from the Pillones Tunnel diversion) and the flows being diverted to Pillones.The second and most important reason is the dataset for spillways, from which AUTODEMA and EGASA had different records.Furthermore, when activating the spillway flows, the reservoir volume decreased to zero.Some hypotheses for this could be as follows: (1) higher flows are coming from the Pillones reservoir; (2) the precipitation distribution needs more stations in the drainage area contributing to the reservoir inflows; and (3) the spillways are wrongly measured or summarized, since giving one uniform flow for the spillways for an entire day dried the reservoir, and it is not the same to have a high discharge for 24 h and just for six hours.The year 2012 was the most problematic, where a high precipitation event triggered the spillways and the errors (see Figure 13).From all of the available General Circulation Model (GCM) datasets, 13 models were selected for both historic (1940-2014) and future (2015-2100) datasets (see Section 3.1 for details).Three main errors were found in these downscaled GCMs: (1) the minimum temperatures were greater than the maximum ones for certain points near the shoreline of the study area, which is something that also happened with the PISCO dataset (see SI 1); (2)  From all of the available General Circulation Model (GCM) datasets, 13 models were selected for both historic (1940-2014) and future (2015-2100) datasets (see Section 3.1 for details).Three main errors were found in these downscaled GCMs: (1) the minimum temperatures were greater than the maximum ones for certain points near the shoreline of the study area, which is something that also happened with the PISCO dataset (see SI 1); (2) the downscaled points were overestimating precipitation in lower parts of the watersheds, where precipitation is scarce; and (3) for some GCM models, there wasn't a smooth transition of temperatures when going from historical to future datasets.
A bias correction was performed for each precipitation data value for each dataset (see Figure 14 for an example of one data point).This process is straightforward and allowed us to mitigate over-or under-estimations of precipitation and temperatures, while buffering the dispersion of the parameters.
From the bias-corrected time series of temperatures, most watersheds followed the same increase patterns between SSPs and future years (see Table S1 in SI 7).Temperature changes were similar for each watershed.Hence, an overall projection for the five watersheds was successfully performed, in which mean increases of 0.73 ± 0.06, 1.20 ± 0.41, 2.25 ± 1.12, and 2.73 ± 1.25 • C could be expected for SSPs 126, 245, 370, and 585, respectively.
Weighting precipitation for an entire watershed is complicated since the spatial distribution of precipitation varies, particularly with the extreme topography changes in Peru, where low or zero precipitation is found in the lowest parts of the watersheds and most comes from higher elevations.The same procedure as temperatures were used and precipitation values accounted for an increase between 0 and 0.98 mm/day (see Table S2 in SI 7).The larger value is not negligible when calculating the total flow produced by its drainage area.Of the five watersheds, Camana received highest increment in mean precipitation, while Yauca had the lowest, which is expected due to its topography does not have the high elevations compared to the other watersheds.Summarizing each SSP mean daily precipitation increase, 0.10 ± 0.07, 0.23 ± 0.11, 0.19 ± 0.17, and 0.34 ± 0.27 mm/day could be expected for SSPs 126, 245, 370, and 585, respectively.lowed us to mitigate over-or under-estimations of precipitation and temperatures, while buffering the dispersion of the parameters.
From the bias-corrected time series of temperatures, most watersheds followed the same increase patterns between SSPs and future years (see Table S1 in SI 7).Temperature changes were similar for each watershed.Hence, an overall projection for the five watersheds was successfully performed, in which mean increases of 0.73 ± 0.06, 1.20 ± 0.41, 2.25 ± 1.12, and 2.73 ± 1.25 °C could be expected for SSPs 126, 245, 370, and 585, respectively.Weighting precipitation for an entire watershed is complicated since the spatial distribution of precipitation varies, particularly with the extreme topography changes in Peru, where low or zero precipitation is found in the lowest parts of the watersheds and most comes from higher elevations.The same procedure as temperatures were used and precipitation values accounted for an increase between 0 and 0.98 mm/day (see Table S2 in SI 7).The larger value is not negligible when calculating the total flow produced by its drainage area.Of the five watersheds, Camana received highest increment in mean precipitation, while Yauca had the lowest, which is expected due to its topography does not have the high elevations compared to the other watersheds.Summarizing each SSP mean daily precipitation increase, 0.10 ± 0.07, 0.23 ± 0.11, 0.19 ± 0.17, and 0.34 ± 0.27 mm/day could be expected for SSPs 126, 245, 370, and 585, respectively.In summary, the models predict an overall increase in precipitation and temperature within the study area, with this increase being directly related to the severity of SSPs, agreeing with the findings by other authors (e.g., [35]).With this preliminary analysis in mind, each previously simulated set of subbasins was computed using the future climate change datasets.As previously mentioned, three different analyses were performed, and similar characteristics were found through all simulations for each evaluation.Examples for each analysis were developed for monthly means and yearly cumulative rates, where lines denote the mean of the ensembles represented by bands for each SSP (see .The means were used as a reference since there is a lot of uncertainty among climate change models.In summary, the models predict an overall increase in precipitation and temperature within the study area, with this increase being directly related to the severity of SSPs, agreeing with the findings by other authors (e.g., [35]).With this preliminary analysis in mind, each previously simulated set of subbasins was computed using the future climate change datasets.As previously mentioned, three different analyses were performed, and similar characteristics were found through all simulations for each evaluation.Examples for each analysis were developed for monthly means and yearly cumulative rates, where lines denote the mean of the ensembles represented by bands for each SSP (see .The means were used as a reference since there is a lot of uncertainty among climate change models. The following discussion uses the Condoroma Reservoir as a typical example of all of the reservoirs.Monthly reservoir evaporation means demonstrate the positive increase in evaporation moving from SSP 126 to 585 (see Figure 15).Moreover, yearly cumulative rates show the same behavior as monthly ones.The mean trend is conservative for years before 2050, but after that the discrepancies between SSPs are evident (see Figure 16).There is also a discrepancy in the cumulative time series between historic and future datasets, which might be due to the fact that not enough historical data were available.When analyzing reservoir inflow mean changes (see Figure 17), high flows are identified during the wet season and no substantial changes are simulated for the dry season.Moreover, the change in the mean inflows is not as pronounced as the overall change for the ensembles.When analyzing yearly cumulative inflows (see Figure 18), the same tendency as for evaporation is found, where the total volumes per year for each reservoir increase proportionally to SSP scenario.However, here the ensembles show a much higher uncertainty when dealing with higher SSP levels.The following discussion for stream flow uses the Tuti stream gage as a typical example.Streamflows exhibit similar behavior as reservoir inflows for both monthly means (see Figure 19) and annual cumulative flows (see Figure 20), with the only difference being related to the stream gage altitude.The monthly means could have a longer dry period compared to the one shown in reservoirs, which are normally located in the highest elevations and where it rains the most.Since stream gages are all connected when no reservoirs are present, this could modify the downstream stream gages' dry periods since flows could be produced in high-elevation areas.increase proportionally to SSP scenario.However, here the ensembles show a much higher uncertainty when dealing with higher SSP levels.The following discussion for stream flow uses the Tuti stream gage as a typical example.Streamflows exhibit similar behavior as reservoir inflows for both monthly means (see Figure 19) and annual cumulative flows (see Figure 20), with the only difference being related to the stream gage altitude.The monthly means could have a longer dry period compared to the one shown in reservoirs, which are normally located in the highest elevations and where it rains the most.Since stream gages are all connected when no reservoirs are present, this could modify the downstream stream gages' dry periods since flows could be produced in high-elevation areas.In all yearly cumulative time series, it is easy to see that the trends of climate change vary significantly, but this change is pronounced towards extreme values; minimum values are buffered and the tendency for the four analyzed SSPs is similar.Future evaporation projections were conducted for all reservoirs except Pasto Grande (Tambo watershed) since no observed records were available, as previously mentioned.The trends for each reservoir are similar; the increase is directly proportional to the previously analyzed gradient of future temperatures.The maximum increases for the year 2100 (with respect to the year 2025) range between 1-3, 1-10, 12-19, and 17-22% for the SSP 126, 245, 370, and 585, respectively (see Table S3 in SI 7).Reservoirs located in the Quilca watershed have the highest increases in evaporation but, as previously mentioned, the increased ranges do not vary a lot from one dataset to another, so the increases in evaporation could be more associated with the latitude and longitude from where reservoirs are located; the highest altitudes in the watersheds have the lowest temperatures which can mitigate evaporation.Moreover, the Aguada Blanca reservoir, which is considered one of the most important reservoirs for the Quilca watershed, has the highest evaporation rate increases, and it is the most downstream reservoir.
Moving forward to reservoir inflows, the increase for the year 2100 (with respect to the year 2025) are 4-42, −1-43, 75-152, 78-216% for the SSP 126, 245, 370, and 585, respectively (see Table S4 in SI 7).From these values, the lowest relative increase is for the Tambo reservoir and the highest is at the Condoroma reservoir.Similar behaviors are presented when analyzing the total volumes for each reservoir.The Condoroma reservoir is still the  In all yearly cumulative time series, it is easy to see that the trends of climate change vary significantly, but this change is pronounced towards extreme values; minimum values are buffered and the tendency for the four analyzed SSPs is similar.Future evaporation projections were conducted for all reservoirs except Pasto Grande (Tambo watershed) since no observed records were available, as previously mentioned.The trends for each reservoir are similar; the increase is directly proportional to the previously analyzed gradient of future temperatures.The maximum increases for the year 2100 (with respect to the year 2025) range between 1-3, 1-10, 12-19, and 17-22% for the SSP 126, 245, 370, and 585, respectively (see Table S3 in SI 7).Reservoirs located in the Quilca watershed have the highest increases in evaporation but, as previously mentioned, the increased ranges do not vary a lot from one dataset to another, so the increases in evaporation could be more associated with the latitude and longitude from where reservoirs are located; the highest altitudes in the watersheds have the lowest temperatures which can mitigate evaporation.Moreover, the Aguada Blanca reservoir, which is considered one of the most important reservoirs for the Quilca watershed, has the highest evaporation rate increases, and it is the most downstream reservoir.
Moving forward to reservoir inflows, the increase for the year 2100 (with respect to the year 2025) are 4-42, −1-43, 75-152, 78-216% for the SSP 126, 245, 370, and 585, respectively (see Table S4 in SI 7).From these values, the lowest relative increase is for the Tambo reservoir and the highest is at the Condoroma reservoir.Similar behaviors are presented when analyzing the total volumes for each reservoir.The Condoroma reservoir is still the one that would have the biggest changes in the future, reaching maximum increases of The following discussion uses the Condoroma Reservoir as a typical example of all of the reservoirs.Monthly reservoir evaporation means demonstrate the positive increase in evaporation moving from SSP 126 to 585 (see Figure 15).Moreover, yearly cumulative rates show the same behavior as monthly ones.The mean trend is conservative for years before 2050, but after that the discrepancies between SSPs are evident (see Figure 16).There is also a discrepancy in the cumulative time series between historic and future datasets, which might be due to the fact that not enough historical data were available.
When analyzing reservoir inflow mean changes (see Figure 17), high flows are identified during the wet season and no substantial changes are simulated for the dry season.Moreover, the change in the mean inflows is not as pronounced as the overall change for the ensembles.When analyzing yearly cumulative inflows (see Figure 18), the same tendency as for evaporation is found, where the total volumes per year for each reservoir increase proportionally to SSP scenario.However, here the ensembles show a much higher uncertainty when dealing with higher SSP levels.
The following discussion for stream flow uses the Tuti stream gage as a typical example.Streamflows exhibit similar behavior as reservoir inflows for both monthly means (see Figure 19) and annual cumulative flows (see Figure 20), with the only difference being related to the stream gage altitude.The monthly means could have a longer dry period compared to the one shown in reservoirs, which are normally located in the highest elevations and where it rains the most.Since stream gages are all connected when no reservoirs are present, this could modify the downstream stream gages' dry periods since flows could be produced in high-elevation areas.
In all yearly cumulative time series, it is easy to see that the trends of climate change vary significantly, but this change is pronounced towards extreme values; minimum values are buffered and the tendency for the four analyzed SSPs is similar.Future evaporation projections were conducted for all reservoirs except Pasto Grande (Tambo watershed) since no observed records were available, as previously mentioned.The trends for each reservoir are similar; the increase is directly proportional to the previously analyzed gradient of future temperatures.The maximum increases for the year 2100 (with respect to the year 2025) range between 1-3, 1-10, 12-19, and 17-22% for the SSP 126, 245, 370, and 585, respectively (see Table S3 in SI 7).Reservoirs located in the Quilca watershed have the highest increases in evaporation but, as previously mentioned, the increased ranges do not vary a lot from one dataset to another, so the increases in evaporation could be more associated with the latitude and longitude from where reservoirs are located; the highest altitudes in the watersheds have the lowest temperatures which can mitigate evaporation.Moreover, the Aguada Blanca reservoir, which is considered one of the most important reservoirs for the Quilca watershed, has the highest evaporation rate increases, and it is the most downstream reservoir.
Moving forward to reservoir inflows, the increase for the year 2100 (with respect to the year 2025) are 4-42, −1-43, 75-152, 78-216% for the SSP 126, 245, 370, and 585, respectively (see Table S4 in SI 7).From these values, the lowest relative increase is for the Tambo reservoir and the highest is at the Condoroma reservoir.Similar behaviors are presented when analyzing the total volumes for each reservoir.The Condoroma reservoir is still the one that would have the biggest changes in the future, reaching maximum increases of 2217 hm 3 for the year 2100.Pasto Grande is not the one with the lowest volumes but, as previously mentioned, the relative change is the lowest for this reservoir.Moreover, besides Condoroma, reservoirs such as Pasto Grande, Pañe, Aguada Blanca, and El Frayle are the ones with the biggest volume changes.There are multiple possible answers for these changes and specific reservoirs would require further analysis.For instance, the Aguada Blanca reservoir has low storage compared to the other reservoirs, but its drainage area is large.The Pañe reservoir, on the other hand, has a small drainage area but its flows would increase considerably.Some relationships could also be performed based on their designed storage, from where Condoroma, Pasto Grande, and El Frayle stood out.
Regarding streamflows (see Table S5 in SI 7), the increase with respect to 2025 is 17-29, 20-59, 118-457, and 144-704% for SSPs 126, 245, 370, and 585, respectively.It is seen that the increases are greater than in the previous analysis, but this increase is not uniform throughout the watersheds.The biggest increases are found in the Puente del Diablo and Tingo Grande subbasins, from where Tingo Grande is inside the drainage area of the Puente del Diablo streamgage.Hence, these increases could be due to errors regarding the climate change datasets for the Tingo Grande drainage area since there are no known parameters or conditions that could be generating this huge change apart from the natural groundwater influence seen in this subbasin [21].Thus, Tingo Grande is thus highly uncertain and needs further assessment.
Furthermore, the Pitay and Sumbay stream gages, which are the two subbasins with diversions, have the lowest flow increases through the four analyzed SSPs, meaning that in the future these subbasins would still rely on diversions.While methods such as unitary flow are available to compare flows in different sub-basins, it is complicated to compare predicted subbasins' streamflows since these are highly tied to reservoirs' operations, which is the main reason by which the analysis was performed per sub-basin or set of sub-basins; moreover, it is not possible to know future reservoir operations.Hence, the flow magnitudes presented in Table S5 in SI 7 should be carefully considered for management purposes since it is shown that simulated results regarding streamflows are more uncertain than previous analyses.It would be useful to analyze the surface-groundwater interactions for each subbasin to reduce this uncertainty.It may also be useful to include evaluating the influence of glaciers in these basins.For the results presented, the Ocoña watershed will be the one with the highest flow increments, followed by the Camana, Tambo, Yauca, and Quilca watersheds.This assumes that the reservoirs would be able to store all of the documented flows, which may not be true knowing the actual reservoirs' capacities.Hence, in case of reservoir failure, Quilca could move to the third position.Additionally, it is known that some of the reservoirs in the region, Pasto Grande for example, is decreasing in volume due to sediment infilling [36,37].To date, no plan exists to remediate this loss of reservoir volume.Thus, reservoirs may not be able to store these increased flows.
There is a linear relationship between SSP and future years.In other words, it could be expected to have similar flows for the year 2100 in SSP 370 compared to the year 2075 and the SSP 585.Moreover, lower SSP scenarios such as SSP 126 tend to preserve maximum flows over time.

Future Scenarios: Water Usage
Current water usages vary among watersheds [20][21][22][23][24][25][26][27][28], with agriculture exhibiting the highest demand.The other usages (potable water for populations, mining, and industry) are considerably lower (see Figure 21, note logarithmic scale).It is difficult to apportion future water usages accurately in each watershed, as discussed below, but by employing some simplifying assumptions, a useful analysis can be presented.Based on migration studies, substantial population increases are expected only in certain areas: the urban areas surrounding the city of Arequipa, and areas associated with agriculture, particularly large agricultural projects such as Majes [38].The rest of the regions had negative or very low migration rates for the recent 2017 census, so population increases are not expected for most watersheds.Population increases are expected mainly in the Quilca watershed [38].However, future population increases cannot be accurately evaluated.For example, the largest cities in each watershed might demonstrate increases in population.Peru is concerned about population migration to the urban areas, so the country is likely to implement some solutions that could mitigate this migration.
Increases in population could have serious impacts on water availability.Thus, four different scenarios based on the total demand were applied to the current water uses shown in Figure 21.The scenarios consisted of increases in the overall water demand per watershed, from where increases equivalent to 25%, 50%, 100%, and 200% were applied.The only watershed with expected problems would be Yauca; if it's demand triples, water needs would exceed supply (assuming all of the additional summer-season flows could be captured and stored).However, this large increase in demand is unlikely due to the fact that, currently, the area experiences negative migration rates [38]; it is a remote watershed that is not expected to see many increases in population, and no major agricultural projects exist in this area (i.e., small-scale or family farming dominates).
Some water shortages might be expected in the watersheds if demand increases excessively.For example, if demand increases by 300% in the region, the demand will exceed supply, causing a shortage.However, excessive increase, while not expected, could be balanced with new ways to improve water storage capacity.For example, storage capacity (reservoirs) should be maintained, i.e., sediment accumulation in reservoirs needs to be evaluated (something that has not been carried out, as previously mentioned) and controlled, if needed.Indeed, sediment accumulation is the biggest problem that reservoirs Based on migration studies, substantial population increases are expected only in certain areas: the urban areas surrounding the city of Arequipa, and areas associated with agriculture, particularly large agricultural projects such as Majes [38].The rest of the regions had negative or very low migration rates for the recent 2017 census, so population increases are not expected for most watersheds.Population increases are expected mainly in the Quilca watershed [38].However, future population increases cannot be accurately evaluated.For example, the largest cities in each watershed might demonstrate increases in population.Peru is concerned about population migration to the urban areas, so the country is likely to implement some solutions that could mitigate this migration.
Increases in population could have serious impacts on water availability.Thus, four different scenarios based on the total demand were applied to the current water uses shown in Figure 21.The scenarios consisted of increases in the overall water demand per watershed, from where increases equivalent to 25%, 50%, 100%, and 200% were applied.The only watershed with expected problems would be Yauca; if it's demand triples, water needs would exceed supply (assuming all of the additional summer-season flows could be captured and stored).However, this large increase in demand is unlikely due to the fact that, currently, the area experiences negative migration rates [38]; it is a remote watershed that is not expected to see many increases in population, and no major agricultural projects exist in this area (i.e., small-scale or family farming dominates).
Some water shortages might be expected in the watersheds if demand increases excessively.For example, if demand increases by 300% in the region, the demand will exceed supply, causing a shortage.However, excessive increase, while not expected, could be balanced with new ways to improve water storage capacity.For example, storage capacity (reservoirs) should be maintained, i.e., sediment accumulation in reservoirs needs to be evaluated (something that has not been carried out, as previously mentioned) and controlled, if needed.Indeed, sediment accumulation is the biggest problem that reservoirs face around the world, resulting in a worldwide loss of around 20% storage capacity [39], also decreasing the lifespan of the structures.There are ways to prevent this from happening in Arequipa (see [40]).In general, the best way to avoid sediment problems in reservoirs is to build the dam in the right place (location).A small reservoir on an extremely muddy river will rapidly lose its capacity, but a large reservoir on a very clear river may take centuries to lose an appreciable amount of storage.The key factor relies on choosing watersheds that naturally release the least quantity of sediments, i.e., watersheds in good hydrological conditions.For the Arequipa case, the lack of vegetation in the arid environment means that significantly more sediments are being transferred into the reservoirs during storm events.Despite the above, if the dam was already built and sediment deposition needs to be evacuated (as is most likely the case for the Arequipa region), there are effective strategies and techniques that can be applied.The problem of sediment accumulation in reservoirs can be handled in three different locations: (1) in the watershed, using known techniques such as slope stabilization and waterway sediment control, a particularly expensive solution for the study area, considering site conditions; (2) in the reservoir, the most appropriate approach for Arequipa, represented by dredging, flushing, hydrosuction, and the avoidance of settling of fine sediments.Though the approach is normally expensive, it could evacuate a significant portion of deposited sediments downstream; and (3) at the dam using techniques such as dam heightening, sluicing, and turbidity current venting.Considering all of the above, a technical-economic evaluation will be required to decide which approach (or combination of them) should be used to ensure enough storage capacity of reservoirs, based on future conditions.
In addition to the above, the models suggest high flow generation during wet seasons, significantly higher than present day flows, not only in the Department of Arequipa, but also in most of the lower sub-watersheds draining to the Pacific Ocean.The phenomenon is locally known as "huaycos", which are a constant threat to water availability, water quality, and the population's welfare [3].While beyond the scope of this water resources evaluation, huaycos need to be addressed with the necessary structural applications to prevent future disasters.Furthermore, our data analysis indicates that once the capacity of the reservoirs is reached (which can occur early in the wet season), the remaining entering water is simply evacuated downstream, which explains why local rivers still increase their flow rates for months during Peruvian summers, even though there are reservoirs built upstream.Even if a small portion of that water was stored, that would increase the availability of the resource for other uses.Storage capacity alternatives are a function of social, ecological, and economic variables, and are represented mainly by the construction of new reservoirs and the implementation of aquifer recharge facilities.However, both options require further scientific evaluation.This improved storage will be especially important during the dry season; according to model simulations, flows do not increase during those periods, but evaporation increases, reducing reservoir storage, which means increases in the future demands during this season could result in water shortages, particularly for unregulated sub-watersheds (i.e., watersheds without reservoirs).The scenario is likely due to the fact that agriculture occurs year-round in southern Peru and domestic water usage does not decrease during dry seasons.Hence, there is a critical need to increase storage in the unregulated watersheds as well as the regulated ones; specifically, storage that can capture the free flow of water during the wet season, wherein the water supply is most affected by climate change.
Given the uncertainty of future demands, improving irrigation efficiency would be prudent, as crop's evapotranspiration will increase, according to model simulations.Irrigation efficiency is defined as the amount of water used by crops, compared to total water applied through irrigation (e.g., [41]).The highest irrigation efficiency reported in the five watersheds is 0.77 for the sprinkling irrigation system of Majes [16].However, the dominant irrigation in the five watersheds is traditional gravity (flooding) methods, which have lower efficiencies (e.g., [42]).If irrigation efficiencies similar to the ones at Majes should be reached, an increase of the efficiency would be needed; this reduced use would likely never exceed demand due to the fact that agricultural water use dominates the domestic water used by the population; population demands are just 0.42%, 0.54%, 0.54%, 13.19%, and 1.03% of the agricultural demands for the Yauca, Ocoña, Camana, Quilca, and Tambo watersheds, respectively.Moreover, high irrigation efficiencies could only be reached with technologies such as CIMS (California Irrigation Management System), which provides water to meet actual crop demand based on a dense network of weather stations.Our model calculations indicate that such savings would account for a total of 90, 191, 368, 211, and 163 MCM for the Yauca, Ocoña, Camana, Quilca, and Tambo watersheds, respectively.For the Camana watershed, the savings would account for storage greater than its Condoroma reservoir.Similarly, water-retaining polymers could dramatically decrease irrigation demands if applied at a large scale.For instance, it has been found that polymers could increase water retention by 60% in arid and semiarid regions [43].
The previous discussion promotes a bright outlook for sustainable water resources in the Arequipa region, though the discussion relies on the assumption that the additional water could be stored, or water consumption could be decreased with efficient practices.
Water uses such as diversions or reservoirs outflows were not considered when analyzing climate change scenarios and water uses against climate change scenarios, since these flows are controlled by the operators.This means that the analysis reflects the possible changes inside its drainage area (sub-basin).Moreover, groundwater was not considered because of the lack of available records (most uses are based on surface water flows and storage).For example, in the Quilca watershed, agricultural water demands associated with diversions such as the Majes, Pampa de Siguas, and others inside the drainage area of the Pitay stream gage were not considered since those heavily depend on the diversions coming from the Camana watershed.Any specific demands with calculations associated with groundwater were subtracted or not considered in the analysis if not found.
Our approach also has limitations since water demands have spatial variability.For example, mining in some watersheds is concentrated in high-altitude areas, while agriculture could be associated with middle-lower terrains.Similarly, it is important to highlight that our analysis did not include environmental water needs, as there was no information available for this important topic.Future studies should follow the methods applied by other authors that have considered the ecological factor as part of their analyses (e.g., [44]).

Conclusions and Recommendations
For the first time, the Arequipa Region of southern Peru has been characterized and evaluated in terms of water management, after extensive analysis consisting of gathering information from private and public entities to build and calibrate a hydrological model, providing an important tool for future management decisions, thus ensuring sustainability of surface water resources.
After carefully analyzing climate change predictions with the ensembles for each SSP, all reservoirs and stream gages in the studied watersheds show the same trend: a clear increase in flows during wet seasons and no substantial changes during the dry seasons.
Evaporation generally increases due to higher temperatures but there is significant uncertainty in the rates derived from the future climate data.The increase is directly proportional to the SSP level under consideration and has similar magnitudes for all months, which is a different behavior than that predicted for flows, and the resulting loss of stored water could trigger significant deficits during dry seasons if not managed correctly.Therefore, reservoir storage (i.e., sediment deposition) would need to be further assessed since the entire Arequipa Region relies on reservoirs to supply water for all uses in the future.
Most importantly, considering both climate change and hypothetical future demands for the five watersheds, there will be water available to satisfy annual demands, but since most of the flows come during the wet season (when reservoirs currently overflow during peak flows) an increase in watershed storage capacity and irrigation efficiency is strongly recommended to ensure future water sustainability, especially during the dry seasons.Similarly, future studies should consider environmental water demands to ensure healthy aquatic ecosystems in the region, following the analysis by Li et al. [44].
Additionally, this study revealed several important deficiencies that need to be addressed to improve water resources management (and sustainability) evaluations in the Arequipa Region: -There is a clear need to have a unified database for all weather-related variables, streamflows, and reservoir data.It has been found that entities such as ANA, National Meteorology and Hydrology Service (SENAMHI), AUTODEMA, and EGASA, among others, have different datasets and the data are not compiled into one master repository.This causes many problems such as not accounting for weather stations from reservoirs included in the PISCO dataset and discrepancies between databases from one entity to another.Similarly, it is important to have the data publicly available in a quality assured and organized format, to ease the difficulty of future analyses.Moreover, information gathered by private entities should also be included and made available to the public.These improvements in data management would be a step forward toward achieving sustainable development not only in the Arequipa Region but also in the entire country.-For modeling purposes, it is important to include solar radiation measurements in weather stations, especially in areas such as reservoirs and irrigation fields where evaporation and evapotranspiration are a critical piece of the water storage and use evaluation.

-
The inclusion of a weather station in the Pasto Grande reservoir would highly improve conceptualization and modeling of this reservoir.-Based on the climate change analysis, reservoirs such as Condoroma, which buffers all inflows coming from the upstream parts of the Camana watershed, may not be able to mitigate the high flows predicted for SSPs higher than 126.In other words, if the Arequipa Region does not follow the conservative SSP, flooding-related (such as fluvial floods and "huaycos") problems downstream of the reservoirs could worsen, which is something that has been threatening many communities within sub-watersheds draining to the Pacific Ocean.-It is encouraged that the PISCO dataset has continuous development.This dataset may not be as accurate as hard data from weather stations but is a product with a good spatial resolution for Peru, derived from raw weather stations data.Moreover, the inclusion of an integrated database, as previously mentioned, would improve its accuracy.-Given the strong topographical characteristics of Southern Peru, it is recommended to improve the weather stations' density and data recording, specially in areas where high water storage is present (i.e., reservoirs and their drainage area).-During model calibration, spillway flows produced unrealistic behaviors for some reservoirs.Thus, it is recommended to record hourly time series of spillway flows since the high magnitude of these flows is not usually constant for an entire day, which made some reservoirs go dry as soon as they were activated.-Irrigation could be included in the PRMS setup.Land uses associated with irrigation have an HRU identifier.Hence, its inclusion would not be complicated.However, to accomplish this, it would be necessary to monitor irrigation schedules and parameters such as soil moisture and aquifer levels below irrigated areas, which is not currently carried out.Moreover, if an aquifer is generated due to intensive irrigation as in the Majes irrigation district (see [16] for more details), a network of piezometers would be needed to include the groundwater flow direction into the simulation, or more advanced models such as MODFLOW [45] could be linked.-It would be useful to include local water consumption for usages such as drinking water in the city of Arequipa, but at least a monthly time series would be needed for each included water use.Currently, estimates for use lack sufficient detail for detailed analyses (for example, to enable a focus of future water supply versus demand monthly during the dry season).-If it is desired to improve reservoir dynamics, for instance, to use water levels in the reservoirs for model calibration), which could significantly improve reservoir response and dynamics, PRMS source code would need to be modified to include rating curves of volume versus elevation.This approach would be beneficial if there was a strong interaction between reservoir and surrounding aquifer or if there were seasonal changes in the surrounding aquifer's water levels, which would allow the model to analyze the elevation at which water in the reservoirs starts to seep into groundwater storage.This approach would require a network of piezometers to monitor the water table surrounding the reservoir.-After calibrating the five watersheds, the scarcity of data throughout the study areas was evident.A significant amount of data is missing in the reservoir system of the Quilca watershed.There are many water intakes and channels with unknown flows in the Dique los Españoles reservoir and the Pañe-Colca channel.Moreover, only the flows going towards the Pillones tunnel are known and not the ones going downstream to the Aguada Blanca reservoir.Furthermore, implementing a stream gage upstream of the Aguada Blanca would improve calibration and water accounting since the drainage area for this location has been considered to generate high flows.
There is a need for a uniform dataset in all reservoirs and diversions since databases had discrepancies.For instance, the Aguada Blanca reservoir showed a difference of around 50 m 3 /s (in reservoir spillway flow) between two source institutions that documented the same information.For the rest of the watersheds, there are no specific locations but there is a need for more stream gages and stations, and it would be important to have these near cities or irrigation projects to properly evaluate future impacts.-It is extremely important to evaluate reservoir sediment deposition against high flows expected with climate change.The actual capacity of the reservoirs is currently unknown, so an evaluation of sediment accumulation is crucial to ensure future water sustainability in the Region.Hydraulic analysis should be performed to evaluate the possible damages to infrastructure and to downstream inhabitants during the increased "huayco" flows.-A more comprehensive analysis of available groundwaterin the five watersheds would be highly useful to evaluate watershed storage capacity, since it could serve as an additional water supply in the future.Groundwater may represent an important way to store water in the future.-A holistic evaluation of future water supply and storage options is recommended.Such analysis would include evaluating the construction of new reservoirs, the implementation of aquifer recharge facilities, or other solutions such as natural infrastructure (e.g., constructed wetlands), water treatment options such as recycled water and desalination.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/su152316270/s1,SI 1: Details of methods and criteria used to develop the hydrological models (continuation from Section 2. Study area and methods).References  are cited in SI 1; SI 2: Simulation file description; SI 3: Description of the selected parameters for calibration process; SI 4: Calibrated reservoir evaporation charts; SI 5: Calibrated reservoir volume charts; SI 6: Calibrated streamflow charts; SI 7: SSP scenarios for all variables.

Sustainability 2023 , 32 Figure 1 .
Figure 1.Geographical boundries of the five watersheds within the Department of Arequipa, with the location of climatic and fluviometric stations, reservoirs, infrastructure flow gages, and water diversions (Source: own work with data from ANA, n.d.).

Figure 1 .
Figure 1.Geographical boundries of the five watersheds within the Department of Arequipa, with the location of climatic and fluviometric stations, reservoirs, infrastructure flow gages, and water diversions (Source: own work with data from ANA, n.d.).

Figure 2 .
Figure 2. Distribution of flows for the "Condoroma", "Pillones", "Bamputañe", "El Pañe", and "Dique los Españoles" reservoirs, with blue arrows representing flows related to hydraulic infrastructure and blue lines denote flows produced by natural drainage.The regulated label means scheduled releases by the operators; spillway label refers to spillway flows; and regulated labels in red represent scheduled flows not recorded by the operators (without time series).

Figure 2 .
Figure 2. Distribution of flows for the "Condoroma", "Pillones", "Bamputañe", "El Pañe", and "Dique los Españoles" reservoirs, with blue arrows representing flows related to hydraulic infrastructure and blue lines denote flows produced by natural drainage.The regulated label means scheduled releases by the operators; spillway label refers to spillway flows; and regulated labels in red represent scheduled flows not recorded by the operators (without time series).

Figure 3 .
Figure 3. Example of elevation, soil type, coverage type, and land uses for the Camana watershed.

Figure 3 .
Figure 3. Example of elevation, soil type, coverage type, and land uses for the Camana watershed.

Figure 4 .
Figure 4. Subbasins generated for each watershed and watershed boundaries established by ANA (black boundary).

Figure 4 .
Figure 4. Subbasins generated for each watershed and watershed boundaries established by ANA (black boundary).

Figure 5 .
Figure 5. Example of modeling setup for the Camana watershed.Model file descriptions are provided in SI 2.

Figure 5 .
Figure 5. Example of modeling setup for the Camana watershed.Model file descriptions are provided in SI 2.

Figure 6 .
Figure 6.Example(Camana watershed)  illustrating difference between precipitation distribution exponentials for the PRMS's "IDE_DIST" method using a gridded approach for the HRUs.Model file descriptions are provided in SI 2.

Figure 6 .
Figure 6.Example(Camana watershed)  illustrating difference between precipitation distribution exponentials for the PRMS's "IDE_DIST" method using a gridded approach for the HRUs.Model file descriptions are provided in SI 2.

Figure 7 .
Figure 7. Calibration approach for each simulated watershed.

Figure 7 .
Figure 7. Calibration approach for each simulated watershed.

Sustainability 2023 , 32 Figure 8 .
Figure 8. Example of sensitivity analysis for the Yauca watershed.The time series for each observed dataset was calibrated.Figures 9-11 represent examples of calibrated time series for evaporation, reservoir volumes, and streamflow, respectively.Reservoir volumes and streamflows are in English units since it is the default model output for PRMS, as previously mentioned.Moreover, evaporation calibrations got Nash-Sutcliffe Efficiency (NSE) greater than one and very low Root Mean Squared Error (RMSE) since the calibration was performed based on monthly means and there was only

Figure 8 .
Figure 8. Example of sensitivity analysis for the Yauca watershed.The time series for each observed dataset was calibrated.Figures 9-11 represent examples of calibrated time series for evaporation, reservoir volumes, and streamflow,

Figure 9 .
Figure 9. Example of calibrated monthly mean evaporation for the Condoroma reservoir.Calibrated evaporation charts corresponding to other reservoirs are in SI 4.

Figure 9 . 32 Figure 10 .
Figure 9. Example of calibrated monthly mean evaporation for the Condoroma reservoir.Calibrated evaporation charts corresponding to other reservoirs are in SI 4. Sustainability 2023, 15, x FOR PEER REVIEW 15 of 32

Figure 11 .
Figure 11.Example of calibrated streamflows for the Huatipa stream gage.Calibrated streamflow charts corresponding to other stream gauges are in SI 6.

Figure 10 . 32 Figure 10 .
Figure 10.Example of calibrated reservoir volumes for the Condoroma reservoir.Calibrated volume charts corresponding to other reservoirs are in SI 5.

Figure 11 .
Figure 11.Example of calibrated streamflows for the Huatipa stream gage.Calibrated streamflow charts corresponding to other stream gauges are in SI 6.

Figure 11 .
Figure 11.Example of calibrated streamflows for the Huatipa stream gage.Calibrated streamflow charts corresponding to other stream gauges are in SI 6.

Figure 12 .
Figure 12.Glacier land uses for the Salamanca subbasin.

Figure 12 .
Figure 12.Glacier land uses for the Salamanca subbasin.

Figure 12 .
Figure 12.Glacier land uses for the Salamanca subbasin.

Figure 13 .
Figure 13.Example of unsuccessful calibration attempt for the Aguada Blanca reservoir volumes.

Figure 13 .
Figure 13.Example of unsuccessful calibration attempt for the Aguada Blanca reservoir volumes.

Figure 14 .
Figure 14.Precipitation bias correction for one grid cell location of the SSP 126, with (A) showing the time series, (B) representing the annual precipitation, and (C) illustrating a box plot for each dataset used and corrected.

Figure 14 .
Figure 14.Precipitation bias correction for one grid cell location of the SSP 126, with (A) showing the time series, (B) representing the annual precipitation, and (C) illustrating a box plot for each dataset used and corrected.

Figure 15 .
Figure 15.Example of monthly mean evaporation for the Condoroma reservoir considering different SSPs.

Figure 15 .
Figure 15.Example of monthly mean evaporation for the Condoroma reservoir considering different SSPs.

Figure 15 .
Figure 15.Example of monthly mean evaporation for the Condoroma reservoir considering different SSPs.

Figure 16 .
Figure 16.Example of yearly cumulative evaporation for the Condoroma reservoir.

Figure 16 .
Figure 16.Example of yearly cumulative evaporation for the Condoroma reservoir.

Figure 17 .
Figure 17.Example of monthly mean volumes for the Condoroma reservoir.

Figure 18 .
Figure 18.Example of yearly cumulative volume for the Condoroma reservoir.

Figure 17 .
Figure 17.Example of monthly mean volumes for the Condoroma reservoir.

Figure 17 .
Figure 17.Example of monthly mean volumes for the Condoroma reservoir.

Figure 18 .
Figure 18.Example of yearly cumulative volume for the Condoroma reservoir.

Figure 18 .
Figure 18.Example of yearly cumulative volume for the Condoroma reservoir.

Figure 19 .
Figure 19.Example of monthly mean flow for the Tuti stream gage, under different SSPs.

Figure 20 .
Figure 20.Yearly cumulative flow for the Tuti stream gage.

Figure 19 . 32 Figure 19 .
Figure 19.Example of monthly mean flow for the Tuti stream gage, under different SSPs.

Figure 20 .
Figure 20.Yearly cumulative flow for the Tuti stream gage.

Figure 20 .
Figure 20.Yearly cumulative flow for the Tuti stream gage.

Table 3 .
Number of subbasins, HRUs, percentage of area covered by the delineation compared with ANA delineations, and areas for the five watersheds.

Table 4 .
Calibration statistics, simulation starting dates, and the percentage of data available for calibration for each subbasin.