Comparative Analysis of Climate Change Impacts on Meteorological, Hydrological, and Agricultural Droughts in the Lake Titicaca Basin

Abstract

The impact of climate change on droughts in the Lake Titicaca, Desaguadero River, and Lake Poopo basins (TDPS system) within the Altiplano region was evaluated by comparing projected 2034–2064 and observed 1984–2014 hydroclimate time series. The study used bias-corrected monthly climate projections from the fifth phase of the Coupled Model Intercomparison Project (CMIP5), under the Representative Concentration Pathway 8.5 (RCP8.5) emission scenarios. Meteorological, agricultural, and hydrological droughts were analyzed from the standardized precipitation, standardized soil moisture, and standardized runoff indices, respectively, the latter two estimated from a hydrological model. Under scenarios of mean temperature increases up to 3 °C and spatially diverse precipitation changes, our results indicate that meteorological, agricultural, and hydrological droughts will become more intense, frequent, and prolonged in most of the TDPS. A significant increase in the frequency of short-term agricultural and hydrological droughts (duration of 1–2 months) is also projected. The expected decline in annual rainfall and the larger evapotranspiration increase in the southern TDPS combine to yield larger projected rises in the frequency and intensity of agricultural and hydrological droughts in this region.

1. Introduction

Drought is a transitory climatological anomaly characterized by below-normal precipitation during an extended period [1,2]. Drought affects the precipitation regime, hydrological cycle, and socioeconomic systems associated with hydrological conditions [3]. These impacts can be differentiated by their time properties. Below-normal precipitation lasting from weeks to years can be called a meteorological drought. A hydrological drought, characterized by a decrease in streamflow or reduced aquifer recharge [4,5,6], could follow a meteorological drought through the hydrological cycle. A meteorological drought over an extended period of time can also develop into an agricultural drought, which is characterized by a decrease in soil humidity over a prolonged time interval [7,8,9]. Each drought type can affect the population and productive activities (irrigation, rain-fed agriculture, livestock farming, water storage in lakes or reservoirs, etc.), with different intensity, duration, and frequency causing scarcity and increasing water demand. Dry spells can also affect rain-fed agriculture even during the wet season in the Andes [10].

Western South America is frequently affected by drought events because of its complex topographic features and high spatio-temporal hydro-climatic variability [11,12,13,14]. This is the case in the Altiplano region shared by Peru and Bolivia (Figure 1), where droughts affected 460,000 farmers during the 1982/83 period [15]. Economic losses in Peru and Bolivia were estimated at $128 million and $88 million for the 1982/83 and 1989/90 periods, respectively [16]. Twenty-five meteorological drought events were also registered between 2000 and 2010 in Peru, causing crop loss, livestock mortality, and the proliferation of pests and diseases, affecting small agricultural producers, urban populations, and electricity production [17].

Hydrological and agricultural droughts could have multiple causes in the Bolivian Altiplano, where nearly 50% of the active population is dedicated to farming [18]. Diversion of water, from the Desaguadero River, the main tributary of Lake Poopó, for irrigation and increasing evaporation related to local and global warming were identified as the main factors behind the recent disappearance of Lake Poopó in 2015 and 2016 [19]. The drying up of Lake Poopó in 1994 and 1995 was associated with strong negative rainfall anomalies [19,20].

Scientific studies aimed at evaluating the impacts of climate change on droughts have mainly focused on the analysis of climate variables such as precipitation, soil moisture, and runoff by using global circulation models (GCM) [4,6]. GCM–based projections of climate change during the 21st century consistently predict an increase in surface air temperature over the central Andes [21,22,23]. Significant long-term warming has been detected for the Altiplano, with the positive trend in daily maximum temperature being stronger than that of the daily minimum [24,25]. Both trends and GCM projections indicate that evapotranspiration is increasing.

GCM results by [26] projected a decrease in summer precipitation over the central Andes south of 12° S. Regarding climate change, [23] predicted consistent drying for the Bolivian Altiplano between a 1997–2008 baseline and 2050 GCM scenarios. A multi-model ensemble analysis of Coupled Model Intercomparison Project phase 3 (CMIP3) by [21] found an almost year-round increase in westerly flow at mid- and upper levels over the central Andes, resulting in decreased moisture transport from the Amazon and the Atlantic Ocean toward the Altiplano during summer, which could reduce summer precipitation by 10% to 30% from 1970–1999 to 2070–2099. In contrast, a recent study [27] observed that the relationship between easterly winds and austral summer precipitation over the Altiplano has weakened in recent years, leading to more uncertainty in future projections.

Climate change and its impacts on precipitation and evapotranspiration over the Altiplano region therefore require further study. In particular, extreme events such as droughts must be addressed, as changes in intensity, duration, and frequency are highly relevant given the implications of global warming scenarios for the Altiplano [18]. Lake Titicaca is considered a wetland of international importance by its significance in terms of ecology, botany, zoology, limnology, and hydrology. It is therefore important to explicitly examine the potential impacts of future climate change across the Altiplano and estimate the sensitivity of each drought type to different scenarios. Climate change projections along the Andes suggest the need for approaches that contribute to the development of capabilities and knowledge for adaptation [28].

The aim of this study is to investigate the impacts of climate change on the characteristics of meteorological, hydrological, and agricultural drought in the Lake Titicaca, Desaguadero River, and Lake Poopó basins (TDPS system). Frequency, duration, and intensity properties are analyzed for each drought type, using multiple GCM, emission scenarios, and hydrological modeling. Drought indices based on historical frequency distributions of variables such as precipitation, soil moisture, and runoff [4,6] are used to investigate the potential impacts of future climate change on different domains of the hydrological cycle in the TDPS system.

2. Study Area

The northern and central parts of the Altiplano are characterized by the presence of lakes Titicaca and Poopó, two large bodies of water connected by the Desaguadero River. The TDPS (Titicaca, Desaguadero, Poopo, Coipasa Salt Pan) system, as this region is known, extends over an area of 145,000 km² (14.2°–20.04° S; 66.4–71.07° W) in Peru and Bolivia, varying in altitude from 3600 to 6550 masl. Lake Titicaca (8500 km²) is the highest (3809 masl.) navigable lake in the world (Figure 1b). The TDPS system is located between the Pacific coastal deserts of Peru and northern Chile and the humid Bolivian Amazon. It is characterized by a semiarid climate, with mean annual precipitation ranging from 1000 mm in the north to 200 mm in the southwestern part [19]. The region exhibits a marked seasonality, characterized by local boundary layer moisture and zonal winds during the wet season (December–March). During the rest of the year, the Altiplano climate develops from the mid-level westerly flow that brings dry air from the Pacific Ocean region [29]. At interannual time scales, negative and positive rainfall anomalies during El Niño and La Niña are observed in the Altiplano, respectively [30,31,32]. At decadal and interdecadal time scales, [33,34] found that during warming decades in the central-western Pacific Ocean, westerlies are intensified at 200 hPa above the Central Andes, producing decadal periods of a hydrological deficit over the northern Altiplano. This indicates that the interannual variability of precipitation between December and February is likely controlled by the same mechanisms over the Altiplano, consistently with the results obtained by [27].

3. Materials and Methods

3.1. Data

Hydroclimatological data from the National Meteorology and Hydrology Services (SENAMHI) of Peru [35] and Bolivia [36] were collected. The observed datasets included monthly precipitation and maximum and minimum temperature from 1964 to 2014 at 138 meteorological stations along the TDPS (Figure 1). To ensure the greatest monthly data availability, a common period (1984–2014) was selected. To assess data quality, the resulting database, made up of 121 rainfall stations, was submitted to the regional vector method (RVM). The RVM assumed the annual rainfall in the same climate zone was proportional between stations, with little random annual variation as a result of rainfall distribution [37,38].

Monthly mean streamflows from nine hydrological stations, based on their basin size (covering 60% of the TDPS area), were also collected. Streamflow data management was performed using HYDRACCESS software [39]. The first part of the 1984–2014 study period comprised the final years of a wet period in the TDPS system, which was followed by a longer and drier subperiod that started in the 1990s. This interdecadal variability had already been identified by several authors e.g., [22,33].

The outputs of 31 GCM from phase 5 of the Coupled Model Intercomparison Project (CMIP5) under the Representative Concentration Pathway 8.5 (RCP8.5) emission scenario was used in this study. The RCP 8.5 scenario was selected because it usually, as it usually yields the largest climate change impacts [40]. The climate projection datasets were collected from the KNMI Climate Explorer, which is a web application (https://climexp.knmi.nl/start.cgi) for statistical analysis of climate data, providing timely climate datasets.

In climate change studies, uncertainty is associated with issues such as emission scenarios [41], climate model structure [6], and downscaling/bias-correction methods [42]. Downscaling using regional climate models (RCMs) is conditioned by climate dynamics on larger (continental) or even global scale [42]. Some RCM studies [22,26] have been conducted along the TDPS system. Previous drought studies, however, did not consider emission scenarios [43].

3.2. Methodology

As input for hydrological modeling, a reference evapotranspiration (ET) dataset was estimated using the Hargreaves and Samani model [44]. This model recommends a simple equation to estimate ET, using maximum and minimum temperature and solar radiation data (Equation (1)).

$$\mathrm{ET}=0.0023 \left(T_m+17.8\right)(T_{max}-T_{min)}{}^{0.5}\times \mathrm{Ra}$$

(1)

where: ET (mm/day), $T_m$ is the mean temperature (°C), $T_{max}$ is the maximum temperature (°C), $T_{min}$ is the minimum temperature (°C), and Ra is extraterrestrial radiation (MJ/m²/d). This method yields results similar to those obtained by the FAO Penman-Monteith equation for the Altiplano [45,46]. Those authors and [47] found that the Penman-Monteith and Hargreaves-Samani methods are in close agreement with data from lysimeters for pasture varieties used as reference crops. A recent study developed in the Peruvian Altiplano also found a significant correlation (r = 0.55–0.87, p < 0.05) between potential ET calculated with Hargreaves-Samani and Penman-Monteith methods [45]. The Hargreaves-Samani method has already been used for climatic and hydrological studies in the Andes [18,48,49].

Finally, the evapotranspiration dataset estimated using the Hargreaves Samani method and the observed precipitation were interpolated by the Kriging method (0.20° × 0.20° grid-cell size) [50], using up to 121 climatic stations.

3.2.1. Bias-Corrected Climate Projections and Selection of GCM Models

To compare the GCM outputs and observed datasets, a common assessment platform was previously provided. Gridded downscaling based on interpolation was performed by using the nearest-neighbor interpolation method. This is an algorithm that stores all the available data and classifies the new data based on a similarity measurement [51,52,53]. As the spatial resolution of the GCM output ranged from 1° to 4°, they were interpolated to the resolution of the precipitation and ET datasets (0.20° × 0.20° grid-cell size). The nearest-neighbor interpolation method is suitable for the study region, as the TDPS system is mostly flat, with 80% of its area in the altitude range from 3650 to 4300 masl. Furthermore, the interannual variability of precipitation during the wet season is likely controlled by the same mechanisms over the TDPS [27,33,34]. Finally, a smoothed spatial distribution of rainfall for the TDPS was obtained.

To quantify uncertainty in projections of climate change impacts on droughts, all climate model datasets were submitted to a bias correction method (precipitation and evapotranspiration estimated with 0.20° × 0.20° grid-cell size). This was done using distribution-based bias correction algorithms that preserved relative changes [54] at a monthly time scale. The normalized data were used in a parametric mapping approach for adjusting the monthly variability of GCM data to match the variability of the observed data (1984–2014), removing unrealistic jumps. Specifically, the approach used a multiplicative algorithm that preserved properties at a seasonal or monthly scale.

The performance of GCM models was evaluated by comparing observed precipitation P and calculated ET (Equation (1) with observed temperature data) with every one of the 31 GCM correspondent outputs for the 1984–2014 period. Three statistical criteria were used: the Pearson correlation coefficient (CC), root-mean-square deviation (RMSD), and coefficient of variation (CV). To facilitate the interpretation of the statistical results, Taylor diagrams ([22,55,56] were used. In the Taylor diagram, the reference (black dot) corresponded to CC, RMSD, and CV values of 1, 0, and 0, respectively. Thus, in the Taylor diagrams, the position of the GCMs dots relative to the reference dot was an integrated indicator of GCM efficiency in reproducing observed data. The shorter the distance between the GCM and the reference position, the closer the GCM and observed variable estimates. To select the best GCMs, skill scores for each GCM and each variable were calculated according to [54]. This allowed the final selection of climate projections from five GCMs: EC-EARTH, HadGEM2-ES, IPSL-C.M5A-LR, MIROC5, and MPI-ESM-LR [57,58,59,60,61] (

Table 1. Description of the five global circulation models selected for analysis of droughts in this study.

Model	Country	Institute Name	Reference
EC-EARTH model	Europe	European Community Earth-System Model	[57]
HadGEM2-ES model	United Kingdom	Met Office Hadley Centre MOHC	[58]
IPSL-CM5B-LR model	France	Institut Pierre Simon Laplace	[59]
MIROC5 model	Japan	Center for Climate System Research/National Institute for Environmental Studies	[60]
MPI-ESM-LR model	Germany	Max Planck Institut fur Meteorologie	[61]

).

3.2.2. GR2M Model

Complex physically-based hydrological models usually include a significant number of parameters to represent every single component of the water cycle. As a consequence, they require very detailed and comprehensive information, which may not be available or could be too costly to collect. The cost of setting and running a model also increases with the time scale used. Taking into account the data availability for the Altiplano and the goals of this study, the GR2M lumped hydrological model was chosen. GR2M is a monthly water balance model [62,63]. The GR2M model has yielded good results in hydrological studies conducted in basins or catchments of the Andes-Amazon region [49], Peruvian Andes [64,65,66], Ecuadorian Andes [67], and Colombian Andes [68]. It has also been applied in Asia [69,70], Europe [71,72], and Africa [73,74,75,76], among other locations.

The GR2M model is based on the transformation of precipitation to streamflow [77] through two functions: (a) a production store function organized around a reservoir, called reservoir-soil (S), and (b) a transfer function governed by the second reservoir, called gravitational water reservoir (R), which is used to estimate the underground exchange (Figure 2).

Initially, water storage for the S reservoir (X₁ in mm, the capacity of the reservoir-soil) is estimated using input (P) and output (ET). P₁ refers to the excess precipitation that directly enters the R reservoir. P₂ refers to the amount of water as a subsurface flow that also contributes to the R reservoir to estimate the groundwater exchange. The sum of P₁ and P₂ corresponds to the effective precipitation P₃ that enters the gravitational water reservoir (R) (maximum capacity set at 60 mm). The contribution to R is instantaneous at the beginning of the time step, and the reservoir then empties gradually. The change in water level in this R reservoir determines the flow rate of water that can be released or stored (R₂) (Figure 2). R₂ is associated with the X₂ parameter (the coefficient of groundwater exchange) (Figure 2). Finally, to determine runoff (Q), the gravitational water reservoir releases water according to a quadratic function. Calibration of the GR2M model with observed streamflow data was carried out by using a shuffled complex evolution metropolis algorithm [78,79] to find optimum model parameter sets at 9 streamflow gauging stations over the 1984–2014 control period.

To estimate the spatial distribution of the X₁ and X₂ parameters, first the observed precipitation and estimated ET were used as input to the hydrological model for nine basins of the TDPS system and the control period (1984–2014) (

Table 2. GR2M model performance (coefficient of determination R² and Nash-Sutcliffe coefficient NS) and model parameters (capacity of production store X₁ in mm and water exchange coefficient (X₂), using observed datasets for nine sub-basins controlled at the hydrological stations identified as red dots/numbers in Figure 1.

N	Station	Country	Optimal Parameters
			R2	NS	X1	X2
1	Ramis	Peru	0.86	0.84	5.60	0.98
2	Unocolla	Peru	0.72	0.76	5.56	1.16
3	Ilave	Peru	0.80	0.81	5.54	0.98
4	Huancané	Peru	0.85	0.84	5.54	1.13
5	Calacoto Desaguadero	Bolivia	0.50	0.66	9.03	0.36
6	Calacoto Maure	Bolivia	0.40	0.60	5.90	0.85
7	Chuquiña	Bolivia	0.42	0.66	6.02	0.99
8	Ulloma	Bolivia	0.30	0.78	9.30	0.44
9	Escoma	Bolivia	0.62	0.58	5.76	1.00

). X₁ and X₂ parameters were obtained by calibration with observed streamflows (NS ≈ 0.6–0.84, Nash Sutcliffe coefficient). Then X₁ and X₂ grid maps (PM) generated from hydrological modeling for the 1984–2014 period were also used for hydrological modeling for the 2034–2064 period. For unmonitored regions, parameters of neighboring basins were used.

Each pixel contained monthly time series of the precipitation (P) and ET datasets (1984–2014) and X₁ and X₂ values. Thus, monthly time series of soil moisture (S) and runoff (Q) obtained from hydrological modeling (2034–2064 period) were estimated to analyze the changes between 1984–2014 and 2034–2064 in drought properties (intensity, duration, and frequency) of meteorological, agriculture, and hydrological droughts using P, S, and Q data series, respectively.

3.2.3. Drought Analysis

Climate change impacts were analyzed at three levels—meteorological, agricultural, and hydrological droughts—using monthly precipitation (P), soil moisture (S), and runoff (Q) datasets for the 2034–2064 period, relative to the 1984–2014 period. To characterize meteorological, agricultural, and hydrological droughts, the standardized precipitation index (SPI) [80,81], standardized soil moisture index (SSMI) [82,83], and standardized runoff index (SRI) [84], respectively were used. These indices quantified the observed anomaly as a standardized departure from a selected probability distribution function that modeled the raw data. To analyze meteorological, agricultural, and hydrological droughts (SPI, SSMI, and SRI), every grid cell of P, S, and Q datasets were considered as a temporal series.

In the central Andes, the rainy season occurs in the austral summer and the dry season in winter. The seasonality is even stronger and the wet season is shorter in the Altiplano region [18,21]. The choice of timescale was related to the objectives to be achieved [85]. The one-month time scale was applied to every one of the SPI, SSMI, and SRI indices given that seasonality, the semiarid climate, and the fact that even 1-month long droughts during the rainy onset (September–November) and wet (December–April) seasons have serious impacts on crop yields in the Altiplano [86], thus requiring meteorological and agricultural droughts to be analyzed at a detailed time resolution. Hydrological droughts, which are usually studied at longer time scales, were examined also at the 1-month time scale for comparison.

The procedure can be described by the following four steps: (i) fit a probability density function to the frequency distribution of the P, S, and Q variables for the one-month time scale, (ii) use the fitted distribution to estimate the cumulative probabilities by numerical solutions methods, (iii) transform the cumulative probabilities to the Z-values of a normal distribution with zero mean and unit variance. Finally, (iv) in this study, a drought event was registered when the 1-month drought index (SPI, SSMI, or SRI) fell below the threshold (−1.0), which was the same for all drought indices. This threshold indicated that the deviation from average conditions exceeded one standard deviation. A drought in the range −1.0 to −2.0 was considered moderate and severe below −2.0 [6].

Probabilistic distributions for the historical period were developed to be used for future projections. This assumed that the distribution parameters would not change in the future and could be extended to analyze climate change [6,87,88,89]. Time series of the drought indices were estimated over a gridded map in the TDPS using a multi-model ensemble. The duration of a given drought event was defined as the consecutive and uninterrupted time period (one or more months) having an observed drought index below the threshold. The frequency was calculated as the total number of months below the threshold over a defined period (1984–2014; 2034–2064). In addition, the frequencies of drought events according to their duration were also calculated for the 30-year periods. The drought intensity was defined as the average value of the drought index (SPI, SSMI, or SRI) below the threshold over a defined period (1984–2014; 2034–2064).

SPI

Monthly SPI was estimated from the precipitation series fitted by a gamma distribution [3,80,81]:

$$f_{\left(X\right)}=\frac{1}{{\beta }^{\alpha }{\mathsf{\Gamma }}_{\left(\alpha \right)}}{x}^{\alpha -1}{e}^{\frac{-x}{\beta }}, x>0$$

(2)

where $\alpha$ > 0 is a shape parameter, $\beta >0$ is a scale parameter, $x$ is the amount of precipitation, ${ \mathsf{\Gamma }}_{\left(\alpha \right)}$ is the gamma function.

SSMI

SSMI was obtained from the soil water storage series resulting from hydrological modeling (S) (Figure 2). To estimate the distribution of soil moisture, the beta probability distribution was used, because the magnitude of soil water storage was limited [83]. A representation of the beta distribution was as follows:

$$f_{\left(X\right)}=\frac{x^{\alpha -1}}{B \left(\alpha ,\beta \right)}\left(1-x\right){ }^{\beta -1} 0<x<1$$

(3)

where $\alpha$ > 0 is a shape parameter, $\beta >0$ is a shape parameter, and $x$ $is$ the proportional amount of soil water storage. The integral of the probability density function was approximated to a summation over (0, 1).

SRI

The runoff series (Q) was fitted by a lognormal distribution [84], with a probability density function defined as:

$$f\left(x,\mu ,\sigma \right) = \frac{1}{x\sigma \sqrt{2}} exp\left[\frac{-{\left(\ln x-\mu \right)}^2}{2{\sigma }^2}\right], x>0$$

(4)

where $\mu$is the scale parameter, $\sigma$is the shape parameter, and $x$ is the amount of run.

To calculate the relative changes in drought intensity, the months with a drought index (SPI, SSMI, or SRI) below the threshold (−1.0) were filtered from the two 30-year datasets. Then, the drought intensity for each grid cell was calculated as the average of the index values for the filtered (drought) months. The change in percentage was calculated as the difference in drought intensity between the 2034–2064 and 1984–2014 periods divided by the drought intensity of the 1984–2014 period. Meanwhile, the total number of drought months below the threshold over a defined period (1984–2014; 2034–2064) was used to calculate the relative changes in drought frequency analogously. The extended wet season (onset and proper wet season, September–April) was used for the calculations because the precipitation during the dry season is very low (close to zero).

4. Results

4.1. GCM Uncertainty Assessment

Thirty-one GCM outputs from phase 5 of the CMIP5 under the Representative Concentration Pathway 8.5 (RCP8.5) emission scenarios were evaluated. GCM outputs were combined with the GR2M model to estimate the changes in variables that were relevant for water resources management in the TDPS. To analyze the impact of the uncertainty of climate projections on future streamflow predictions, GCM datasets were compared with observed datasets. For conciseness, assessment of only three subbasins of the Lake Titicaca basin to the north and three subbasins of the Desaguadero River in the south are shown in Figure 3, Figure 4 and Figure 5.

Taylor diagrams for the 1984–2014 control period were plotted using climate model datasets for three sub-basins of the Lake Titicaca basin to the north and three sub-basins of the Desaguadero River to the south (Figure 3 and Figure 4). The diagram shows the performance of the GCM before (red dots) and after (blue dots) correction. A high correlation coefficient between the non-corrected GCM and observed datasets is exhibited (r ~ 0.8, p < 0.05), centered root mean square deviation (RMSD) from 40% to 100% and coefficients of variation mostly above 40%, evidenced a large dispersion of the models’ performance (Figure 3).

The bias-corrected datasets (blue dots in Figure 3) showed less dispersion, improved correlation coefficients (r ~ 0.90, p < 0.05), and reduced RMSD at about 40% (Figure 3). As expected, the GCM non-corrected ET datasets showed less dispersion, lower RMSD, and higher correlation with observed datasets than the precipitation datasets (Figure 4). However, the bias-corrected ET datasets still showed improved performance, with a correlation coefficient higher than 0.85 (p < 0.05) and RMSD below 20%.

Finally, five global climate models with the highest skill scores were selected (i.e., EC- EARTH, hadGEM2-ES, IPSL-CM5A-LR, MIROC5, MPI-ESM-LR) [57,58,59,60,61] (Table 1). The MIROC-ESM and MPI-ESM-LR models have previously been used to study the impact of climate change in the Bolivian Altiplano [22]. Precipitation and ET monthly mean analysis using the bias-corrected datasets show that the five GCMs represent quite well the observed annual cycle of precipitation in the Altiplano (Figure 5a), characterized by a dry (May–Sep) and a wet (Dec–Mar) season. A slight overestimation of precipitation was observed during the wet season for some basins (Figure 5a). For potential evapotranspiration, the five models also represented quite well the monthly fluctuations of estimated ET in the six basins (Figure 5b).

4.2. GR2M Model Performance

Table 2 shows the Pearson coefficient of determination R², Nash-Sutcliffe (NS) coefficient, and model parameters (X₁, X₂) for nine sub-basins of the TDPS. In most sub-basins, R² and NS were high (NS ~ 0.60–0.84). They were also higher for the Lake Titicaca tributaries (Peru) than for the Desaguadero River ones. To exemplify the simulation of streamflows from observed precipitation, two sub-basins controlled at stations Ramis (NS = 0.84, Peru) and Calacoto Maure (NS = 0.6, Bolivia) are shown in Figure 6 as representatives of the northern and southern parts of the TDPS system. To exemplify the simulation of streamflows from observed precipitation, two sub-basins controlled at stations Ramis (NS = 0.84, Peru) and Calacoto Maure (NS = 0.6, Bolivia) are shown in Figure 6 as representatives of the northern and southern parts of the TDPS system. Ramis was chosen as the largest and most representative northern sub-basin and theCalacoto Maure sub-basin was selected because it is not affected by Lake Titicaca fluctuations.

The X₁ values suggested a greater storage capacity of the reservoir-soil in the Calacoto Maure basin than in the Ramis one. X₂ values suggested a greater groundwater exchange in Ramis than in Calacoto Maure (Table 2). This was probably associated with the relatively steep topography of the Ramis basin and the semi-flat landscape of the Calacoto Maure one (Figure 1).

4.3. Spatial Distribution of Changes in Drought Characteristics

To evaluate the spatial distribution of changes in temperature, precipitation, and evapotranspiration (1984–2014) associated with future climate scenarios (2034–2064) (Figure 7), a multi-model ensemble mean was calculated using the selected models (EC- EARTH, hadGEM2-ES, IPSL-CM5A-LR, MIROC5, MPI-ESM-LR) for the 2034–2064 period along the TDPS system.

Warming in the 0.5–3.5 °C range is observed in Figure 7a, with the largest temperature rise occurring in the southeastern and the smallest one in the northwestern TDPS (0.5−1.5 °C). According to model projections, the northern and central TDPS will experience a slight increase of 3–6% in future annual precipitation (Figure 7b), while a decrease of up to ~3% is expected in the southern TDPS in Bolivia (Figure 7b). The spatial distribution of ET changes was similar to that of temperature changes: from a 5–8% increase in evapotranspiration in the northern TDPS to 8–10% in the southern zone (Figure 7c).

Because of the strong seasonality of rainfall and the semiarid climate of the TDPS, the extended wet season between September and April was also prioritized for analysis. Monthly relative changes (%) in precipitation and evapotranspiration (

Table 3. Precipitation and evapotranspiration (ET) changes (%) during the wet season (September–April) for the 2034–2064 period relative to 1984–2014 for the basins controlled at the Ramis, Ilave (Peru), and Chuquiña (Bolivia) stations.

Precipitation Changes %	September	October	November	December	January	February	March	April
Ramis	−2.4	−1.8	3.3	8.1	7.3	6.4	6.2	5.4
Ilave	−0.1	−1.5	2.9	6.1	8.2	4.5	3.7	4.9
Chuquiña	−7.4	−7.8	−3.7	3.7	6.0	2.0	1.6	−0.2
Evapotranspiration Changes %
Ramis	8.7	7.9	7.9	7.3	6.8	6.8	7.7	7.7
Ilave	8.9	7.9	8.3	8.7	7.5	7.9	8.8	8.5
Chuquiña	10.1	9.6	10.2	9.5	7.7	8.6	8.8	9.2

) and the percentage changes in the annual coefficient of variation for the 2064–2034 and 2014–1984 periods were estimated (

Table 4. Coefficient of annual variation (%) for the control period (1984–2014) and future period (2034–2064), for the basins controlled at the Ramis, Ilave (Peru), and Chuquiña (Bolivia) stations.

Coefficient of Variation (%)	Ramis	Ilave	Chuquiña
Annual P, control period 1984–2014	83.1	93.4	93.2
Annual P, future period 2034–2064	85.1	95.7	97.0
Annual ET—control period 1984–2014	13.9	17.4	20.3
Annual ET—future period 2034–2064	13.7	16.9	20.0

). Finally, three sub-basins controlled at the Ramis, Ilave (north), and Chuquiña (south) stations were selected for additional comparisons (Table 3 and Table 4). Table 3 shows that the largest monthly rise in precipitation would occur in December or January (the wettest month) in the three sub-basins. The largest decrease occurs in September and October in the southern Chuquiña sub-basin. Table 4 shows a projected rise (2–4%) of the coefficient of variation of annual rainfall; this is largest in the Chuquiña sub-basin.

4.4. Changes in Drought Intensity, Duration, and Frequency

Figure 8a–f shows the spatial distribution of changes in drought intensity and frequency for 2034–2064, compared to 1984–2014, according to the model’s ensemble. Significant differences (%) were observed among the results of meteorological, agricultural, and hydrological droughts. Little change (from −5 to +5%) is expected in the intensity of meteorological droughts in the northern and central TDPS regions (Figure 8a). A diminishing frequency (0 to −20%) of meteorological droughts was observed in the Titicaca basin (northern TDPS) (Figure 8b). In contrast, a moderate rise (5–10%) in meteorological drought intensity and a significant (10–40%) increase in frequency is expected for the southern TDPS (Figure 8a,b). At the annual scale, the precipitation did not show a clear pattern of change in the northern TDPS. Although for some regions of the southwestern TDPS, the precipitation predominantly decreases [18]. The sign (negative/decline in the north and positive/rise in the south) and magnitude (moderate) of the changes in intensity and frequency of meteorological droughts are first associated with the predicted moderate changes in rainfall (Figure 7b) for the climate change scenario 2034–2064, as meteorological droughts are defined only by rainfall anomalies (SPI index). The large changes in frequency are also associated with the expected precipitation decline onset of the rainy season (September–November) for the 2034–2064 period (Table 3).

A rise in intensity (0–20%) and frequency (up to 90%) of agricultural droughts (Figure 8c,d) is predicted for almost the entire TDPS, being greater in the southern TDPS (Figure 8c,d). This is to be expected, as the SSMI index (and soil water balance) used to define agricultural droughts is calculated with the GR2M model using both precipitation P and reference evapotranspiration ET as input variables. The larger changes in the southern TDPS are related to both the reduced precipitation (~3%) and the increased evapotranspiration (~10%) predicted by GCM in this region (Figure 7b,c). A greater increase in frequency than intensity suggests that the interaction between the ET rise and interannual rainfall variability plays a more sensitive role in hydrological modeling when soil moisture content and runoff are estimated.

The increase in intensity (up to 20%, Figure 8e) and frequency (up to 90%, Figure 8f) of hydrological droughts follows a similar spatial pattern as agricultural droughts in the TDPS. The greater frequency (~90%) in the central region (Figure 8f) is attributed to a lower underground exchange estimated for basins controlled at the Ulloma and Calacoto Desaguadero stations (Table 2). In contrast, the regions with the highest underground exchange, identified for the basins controlled at the Chuquiña and Calacoto Maure stations (Table 2), have the lowest frequency increments of hydrological droughts (~25%) (Figure 1 and Figure 8f). This is associated with the close relationship between soil moisture and runoff, as both variables were calculated by the GR2M model using the same rainfall and evapotranspiration inputs. In particular, the similarities in the changes in the intensity of hydrological and agricultural droughts in the southern TDPS are due to the expected decrease in rainfall and increase in evapotranspiration in this region (Figure 7b,c).

To analyze the change in the duration of each drought type, the relationship between the number of droughts (frequency) and their length (duration) was plotted on the same graph for six sub-basins controlled at the Ramis, Unocolla, Ilave, Chuquiña, Calacoto Maure, and Calacoto Desaguadero stations. Our results from the five selected GCMs showed no significant changes for meteorological droughts lasting more than 4 months (Figure 9). Nevertheless, meteorological droughts lasting between 1 and 3 months are more frequent under the RCP 8.5 than in the baseline scenario.

Most GCM-based agricultural drought duration curves are above the baseline, particularly for droughts lasting between 2 and 4 months (Figure 10), evidencing a general increase in the frequency of this type of drought for all sub-basins. This is consistent with agricultural drought changes found in intensity and frequency analyses (Figure 8c,d). It should be noted that no agricultural drought lasting more than 5 months has been identified using SSMI, either for the baseline or RCP 8.5 scenarios, except for the EC-EARTH model. Hydrological drought frequency also increased for most climate change scenarios, particularly for droughts lasting between 1 and 4 months, in contrast, hydrological droughts lasting between 5 and 6 months are gradually reduced (Figure 11).

For the 2034–2064 scenario, the frequency of short-duration droughts lasting 1 to 2 months would increase on average by 30%, 76%, and 61% for the meteorological, agricultural, and hydrological types, respectively (Figure 9, Figure 10 and Figure 11). Whereas the frequency of droughts lasting between 3 and 4 months would increase by 1%, 7%, and 11% for the meteorological, agricultural, and hydrological types, respectively (Figure 9, Figure 10 and Figure 11). These results were consistent with the spatial patterns identified in the frequency of meteorological (Figure 8b), agricultural (Figure 8d), and hydrological droughts (Figure 8f). It is worth mentioning that under the RCP 8.5 scenario, a large increase in the frequency of meteorological droughts lasting 1 month was observed (Figure 9). In contrast, the impacts of climate change on the frequency of agricultural and hydrological droughts will be distributed between durations of 2 and 6 months (Figure 9 and Figure 10).

The frequency of specific drought durations using the 3-month SPI, SSMI, and SRI was also evaluated. A similar temporal pattern in the frequencies of droughts according to their duration using the 1- and 3-month indices was found. The corresponding Figures S1–S3 are provided as supplementary information.

A further evaluation was carried out (Figure 12) to show the changes in the frequency of moderate and severe droughts for the Ilave, Ramis, and Desaguadero basins, controlled at the Ilave, Ramis (Peru), and Chuquiña (Bolivia) stations, respectively (black dots 1, 3, and 7 in Figure 12). Meteorological drought frequency increased in most of the five GCMs scenarios for the southern TDPS. Indeed, an increase in severe droughts is predicted for every basin analyzed (Figure 12a–c).

Fewer moderate meteorological droughts were predicted for basins controlled at the Ramis and Ilave stations (northern TDPS) (Figure 12a,b). This can be explained by the slight rainfall increase expected (3–6%) (Figure 7b) to occur during the wet season (December–April) in the central and northern regions of the TDPS (Table 3). In contrast, the rise in the number of severe meteorological droughts was bigger in the southern basin controlled at the Chuquiña station (Figure 12c). This can be explained by precipitation decreases between September and November in the southern regions (Table 3, Figure 7).

The frequency of agricultural droughts, as estimated from the SSMI, is expected to increase from the baseline (1984–2014) to every one of the five GCM 2034–2064 projections (Figure 12d–f). On average, over 100 moderate droughts and 55 severe droughts were projected for the RCP 8.5 scenarios. Dramatic increases in the frequency of severe droughts, as estimated from the SRI are also expected (Figure 12g–i). The rise in the number of moderate and severe hydrological droughts was bigger in the southern basin controlled at the Chuquiña station than in the northern basins controlled at the Ramis and Ilave stations (Figure 12g–i). This can be explained by the bigger projected ET increase between September and April, and the bigger projected precipitation decrease for the southern TDPS between September and April, as precipitation and ET were used to estimate SRI from the GR2M modeling (Table 3).

5. Discussion

The results of this study show an increase in drought severity in most of the Altiplano region: The impact of droughts is more evident between September and November. In particular, the significant rise in the monthly frequency of hydrological and agricultural droughts could be associated with the expected reduction in rainfall in those months (the rainy season onset). An increase of drought frequency during this season is worrying, as the water demand of agricultural systems reaches its peak during this season [19,90]. A recent article [91] documented GCM projected temperature and precipitation changes of the same signal as this study between September and November over the Altiplano. The impact of these changes would be stronger in the southern TDPS. Indeed, trends in the duration and magnitude of drought events for the Bolivian Altiplano in the last three decades show a similar spatial pattern [18].

Meteorological droughts (duration of less than 3 months) propagate nonlinearly into agricultural and hydrological droughts (duration greater than 3 months). This propagation, shown by SRI and SSMI indices when agricultural and hydrological effects, respectively, are analyzed, is associated with the close relationship between soil moisture and runoff. Both variables are calculated from GR2M modeling using as inputs both rainfall and evapotranspiration; the latter is not considered by SPI analysis for analyzing meteorological droughts. An increase in potential evapotranspiration (from 4 to 10%) is explained by rising temperatures (from 0.5°C to 3.5 °C). In addition, because the ET increase is greater during the onset of the wet season (September–December) it should further contribute to more frequent hydrological and agricultural droughts.

The projected precipitation increases in January (the wettest month) and the rise of the coefficient of variation of annual precipitation suggest an intensification of the hydrological cycle. This is consistent with recent studies that suggest rainfall concentration [10,92,93]. Meteorological droughts in the Altiplano have been associated with El Niño–Southern Oscillation (ENSO events) [94,95]. Indeed, rainfall over the Altiplano decreases during El Niño episodes [30,31]. In addition, [96] also found that warming episodes over the Eastern Pacific Ocean also reduce precipitation over the As most GCMs, including the GCMs selected for this study (EC- EARTH, hadGEM2-ES, IPSL-CM5A-LR, MIROC5, MPI-ESM-LR), project an increase of ENSO episodes at the end of the 21st century [97]. The Pacific teleconnection may explain, at least partially, the increase in the frequency of droughts in the TDPS system. Indeed, studies of precipitation variability using observed datasets have documented trends towards drier conditions in the Altiplano [18].

In this study, a potential issue source of uncertainty is the choice of reference datasets. The potential uncertainty ranges of GCMs make a full characterization impossible. For comparison reasons, only five GCMs with the highest skill score were selected [42,54]. The results above confirm that significant uncertainties exist in future drought simulations. This may arise from different sources such as emission scenarios, GCM structure, and downscaling/bias-correction methods [4,5,6]. In this study, every GCM was downscaled and bias-corrected by algorithms that preserve properties [54], which were comprehensively validated at a seasonal time scale.

Uncertainty of the study results is likely to be influenced by the GCM outputs, in spite of the selection carried out previously. The evaluation of the capacity of the GCM to describe the climate variability (particularly precipitation) at interannual and decadal-interdecadal time scales is suggested. This is still not very frequent in climate change assessments and it was beyond the scope of this study, as the Altiplano climate spatio-temporal variability is not sufficiently known [98].

Moreover, limitations on hydrological model performance are likely associated with input data uncertainty, simplified hydrological processes, and the propagation of drought from meteorological to hydrological and agricultural systems [64]. Indeed, the most prominent disagreements are found for the moderate and extreme frequency of agricultural and hydrological droughts characterized by propagation from meteorological droughts [4,5,6]. However, our results show that it is possible to employ the GR2M model for streamflow simulations in the Altiplano.

6. Conclusions

The impact of climate change on meteorological, agricultural, and hydrological droughts has been evaluated using drought indices estimated from global climate models under the RCP 8.5 emission scenario. The 2034–2064 period relative to 1984–2014 was used to identify expected drought changes in the Altiplano region (the Lake Titicaca, Desaguadero River, and Lake Poopo basins-TDPS system, Peru-Bolivia). According to the ensemble of five selected GCMs, greater temperature and evapotranspiration rises are expected in the southern region of the TDPS system than in the northern one for the 2034–2064 scenario. On the other side, a slight increase in annual precipitation is expected in the northern TDPS and a moderate decrease in the southern TDPS.

Drought projected changes show a north-south gradient, becoming more intense and frequent in the central and southern Altiplano. Furthermore, a nonlinear propagation from meteorological droughts to agricultural and hydrological systems results in drought frequency and duration being amplified from meteorological to hydrological and agricultural droughts. This amplification, shown by SRI and SSMI indices when agricultural and hydrological effects, respectively, are analyzed, is associated with the close relationship between soil moisture and runoff. A similar pattern is observed for agricultural and hydrological drought intensity for most of the TDPS system.

The frequency of meteorological droughts lasting between 3 and 4 months would remain basically the same, but a significant rise in the frequency of meteorological droughts lasting between 1 and 2 months is expected for the climate change scenario for the entire TDPS, that would develop between the dry season and the rainy season onset.

The results obtained in this study are someway limited by the modest performance of most global climate models in the Altiplano region. The relatively large dispersion of the GCM outputs and the associated high uncertainty, particularly for precipitation, propagates from meteorological to hydrological and agricultural droughts estimates. It is clear that there is significant room for improvements in GCM simulations for the Altiplano. A better understanding of the Altiplano spatio-temporal variability of precipitation at different scales and its associated causes could contribute to reducing GCM uncertainty.

Climate change would likely increase the frequency, intensity, and duration of drought events in the Altiplano. Potential impacts of climate change on droughts, extreme precipitation events, or climatic parameters, such as dry-day frequency and wet-day frequency, can contribute directly or indirectly to increased environmental problems in the Andean Altiplano (soil erosion, landslides, floods, dry spells). Adaptation strategies should be developed to mitigate future impacts of climate change on the very vulnerable and growing population of the Altiplano and its sensitive natural systems. This study should provide relevant information for the development of more effective drought management plans.

The lack of observations is one of the main limitations for the performance of agro-climatic and hydro-climatic models in high mountain regions [86,87,88,89]. Simulated climatic datasets from high-resolution models could support the lack of datasets for agro-climatic and hydro-climatic modeling in high mountain regions [26,99,100].