Modeling Lake Titicaca's water balance: the dominant roles of precipitation and evaporation

In the face of climate change and increasing anthropogenic pressures, a reliable water balance is crucial for understanding the drivers of water level fluctuations in large lakes. However, in poorly gauged hydrosystems such as Lake Titicaca, most components of the water balance are not measured directly. Previous estimates for this lake have relied on scaling factors to close the water balance, which introduces additional uncertainty. This study presents an integrated modeling framework based on conceptual models to quantify natural hydrological processes and net irrigation consumption. It was implemented in the Water Evaluation and Planning System (WEAP) platform at a daily time step for the period 1982–2016, considering the following terms of the water balance: upstream inflows, direct precipitation and evaporation over the lake, and downstream outflows. To estimate upstream inflows, we evaluated the impact of snow and ice processes and net irrigation withdrawals on predicted streamflow and lake water levels. We also evaluated the role of heat storage change in evaporation from the lake. The results showed that the proposed modeling framework makes it possible to simulate lake water levels ranging from 3808 to 3812 m a.s.l. with good accuracy (RMSE = 0.32 m d⁻¹) over a wide range of long-term hydroclimatic conditions. The estimated water balance of Lake Titicaca shows that upstream inflows account for 56 % (958 mm yr⁻¹) and direct precipitation over the lake for 44 % (744 mm yr⁻¹) of the total inflows, while 93 % (1616 mm yr⁻¹) of the total outflows are due to evaporation and the remaining 7 % (121 mm yr⁻¹) to downstream outflows. The water balance closure has an error of −15 mm yr⁻¹ without applying scaling factors. Snow and ice processes, together with net irrigation withdrawals, had a minimal impact on variations in the lake water level. Thus, Lake Titicaca is primarily driven by variations in precipitation and high evaporation rates. These results will be useful for supporting decision-making in water resource management. We demonstrate that a simple representation of hydrological processes and irrigation enables accurate simulation of water levels. The proposed modeling framework could be replicated in other poorly gauged large lakes because it is relatively easy to implement, requires few data, and is computationally inexpensive.

1 Introduction

1.1 On the need for an integrated water balance in large lakes

Lakes are water reservoirs of vital importance for the development of regions as they provide many ecosystem services, including fisheries, water supply, tourism, and energy generation (Sterner et al., 2020). However, these services can be impacted by fluctuating lake water levels. For instance, Yao et al. (2023) showed that, in the period from 1992 to 2020, there was a significant decrease in water levels in 43 % of natural lakes (457), an increase in 22 % (234), and nonsignificant trends in 35 % (360). Understanding the main drivers of fluctuations in water levels is crucial for effective lake management, which requires a realistic water balance that accounts for both natural processes and anthropogenic pressures (Wurtsbaugh et al., 2017). Several studies on large lakes (> 500 km²) (e.g., Rientjes et al., 2011; Vanderkelen et al., 2018; Wale et al., 2009) have estimated water balance under the assumption that net water withdrawals in the contributing catchments are negligible. However, this assumption may no longer be valid due to changing climate conditions and increased competition for water uses, potentially leading to reduced upstream inflow (Wurtsbaugh et al., 2017). For example, Schulz et al. (2020) demonstrated that net withdrawals for irrigation exacerbate the decline in storage at Lake Urmia, which is also impacted by climate change.

In large lakes, it is essential to adopt integrated water balance modeling, which represents the interactions and feedbacks between natural hydrological processes and water management within a single modeling framework (Niswonger et al., 2014). Unlike traditional decision support systems applied to large lakes (Hassanzadeh et al., 2012), which typically simulate natural flows and irrigation water requirements independently, integrated modeling enables these processes to be simulated in a coupled and dynamic manner. A few studies have attempted to address an integrated water balance in a large lake hydrosystem (e.g., Hosseini-Moghari et al., 2020; Lima-Quispe et al., 2021). These studies tend to focus on some components of the water balance, while others address them superficially. For example, Hosseini-Moghari et al. (2020) focused on estimating the upstream inflow of Lake Urmia, with less emphasis on direct precipitation over the lake and evaporation, which both play very important roles in the water balance of large lakes (Gronewold et al., 2016). This is partly because large lake hydrosystems involve numerous hydrological processes, and there are often insufficient data to represent these processes in detail and evaluate them comprehensively. The issue is further complicated in transboundary lake hydrosystems, where hydrometeorological monitoring is not always coordinated (Gronewold et al., 2018).

1.2 Challenges in estimating the water balance of large lakes

It is widely recognized that large lakes have a major influence on regional climate (Scott and Huff, 1996; Su et al., 2020). For example, it has been observed that direct precipitation over the African Great Lakes is more intense than in their surrounding areas (Anyah et al., 2006; Kizza et al., 2012; Nicholson, 2023; Thiery et al., 2015). According to Scott and Huff (1996), this is due to the differences in heat capacities between the lake surfaces and the surrounding areas as well as the large amount of moisture that lakes provide for the lower atmosphere, which can lead to increased cloudiness and precipitation over lakes. Estimating precipitation based on ground stations, which are mainly located in the surrounding areas, can lead to inaccuracies. Despite the current availability of remotely sensed datasets, they have been shown to still have significant biases (Hong et al., 2022; Satgé et al., 2019). Regarding direct evaporation, this depends not only on meteorological conditions, but also on the size of the lake, water depth, and water clarity, which all influence the energy balance due to changes in water temperature and vertical mixing (Lenters et al., 2005). Thus, estimating lake evaporation based on meteorological data alone can lead to inaccuracies (Bai and Wang, 2023). The energy balance method is considered to be one of the most appropriate and accurate for estimating evaporation (Lenters et al., 2005) but requires large quantities of data, meaning that it is generally difficult to implement. The original Penman formulation, which does not include changes in heat storage, has been used to estimate lake evaporation (e.g., Kebede et al., 2006; Lima-Quispe et al., 2021). However, Blanken et al. (2011) observed a 5-month delay between peak net radiation and evaporation due to heat storage in Lake Superior in North America. One of the limitations of estimating the change in heat storage is the lack of water temperature data at different depths, and so models are used to simulate the thermal stratification dynamics of the water (e.g., Antonopoulos and Gianniou, 2003).

For upstream inflow, ideally, measured streamflow data will be available. However, there are always ungauged catchments contributing to lakes (Wale et al., 2009). Upstream inflow is mostly estimated with hydrological models (e.g., Rientjes et al., 2011; Zhang and Post, 2018). For the ungauged catchments, regionalization methods are applied based on the parameters obtained in the gauged catchments (e.g., Guo et al., 2021). Basic lumped rainfall-runoff simulations may fall short due to the complex interplay of natural hydrological processes and water management. Integrating both aspects is essential for hydrosystems under significant anthropogenic pressure (e.g., Ashraf Vaghefi et al., 2015; Fabre et al., 2015; Hublart et al., 2016). In high mountain catchments, snow and ice processes significantly impact hydrological responses. Estimating melt is challenging because it is difficult to obtain accurate forcing data (e.g., precipitation and temperature) in high-elevation areas where the measurement network is very sparse (Ruelland, 2020) as well as control data (e.g., upstream-area streamflow and glacier mass balance). In this context, temperature-index approaches (Hock, 2003) are more suitable than energy balance approaches and can produce simulations with acceptable accuracy using a reduced number of parameters and forcing data (e.g., Ruelland, 2023). In terms of water management, according to Wu et al. (2022), 60 % of freshwater withdrawals worldwide are made for agricultural irrigation. Catchment-scale irrigation has been addressed using approaches based on soil water deficit (Kannan et al., 2011; Shadkam et al., 2016) and those that additionally consider irrigation scheduling (Githui et al., 2016; McInerney et al., 2018). Regardless of the approach chosen, one of the limitations is the lack of measured irrigation data (McInerney et al., 2018), which hinders the evaluation of irrigation simulations. This evaluation can consequently only be undertaken indirectly by attempting to more realistically reproduce observed outlet discharges and accounting for net consumption for water uses within the catchments (e.g., Fabre et al., 2015; Hublart et al., 2016).

Net groundwater exchanges are neglected in some studies under the assumption that these fluxes are very small (Duan et al., 2018; Lima-Quispe et al., 2021). Lake–groundwater interactions have been addressed using conceptual (Parizi et al., 2022) and physically based models (Vaquero et al., 2021; Xu et al., 2021), chemical and isotopic balances (Bouchez et al., 2016), and the water balance (Chavoshi and Danesh-Yazdi, 2022). The water balance is a fairly easy option, as the net flux is the result of the other components. However, it is crucial to dispose of accurate measurements or estimates of the other water balance terms to avoid propagating uncertainty (Chavoshi and Danesh-Yazdi, 2022). Hydrochemical or isotopic analyses are considered accurate (Bouchez et al., 2016) but can be very costly for a large lake. Modeling approaches are often limited by data availability (Barthel and Banzhaf, 2016), particularly when using models to dynamically simulate surface water–groundwater interactions (Xu et al., 2021). On the other hand, downstream outflow in exorheic lakes can be estimated by direct measurements (Chebud and Melesse, 2009), a rating curve relating lake level to outflow (Sene and Plinston, 1994), and as a residual of the water balance (Duan et al., 2018).

1.3 Placing Lake Titicaca in the context of an integrated water balance

Lake Titicaca, located on the Altiplano of the tropical Andes of South America, is one of the highest large lakes in the world and an interesting case study for an integrated water balance. This lake is part of a vast endorheic catchment and is connected by the Desaguadero River to Lake Poopó in Bolivia (see Fig. 1) (Lima-Quispe et al., 2021). As a transboundary lake shared by Peru and Bolivia and as a poorly gauged hydrosystem, it faces many of the aforementioned challenges. The region experiences significant interannual climate variability (Garreaud and Aceituno, 2001), which, coupled with complex water management issues (Revollo, 2001), intensifies the difficulties in managing Lake Titicaca. These challenges include extreme hydrological events (droughts and floods), lake releases, and water pollution (Revollo, 2001; Rieckermann et al., 2006). Water levels measured in Puno (Peru) have fluctuated by approximately 6 m over the past century, with the lowest recorded in 1943–1944 and the highest recorded in 1984–1986 (Sulca et al., 2024), causing USD 125 million in flood damage (Revollo, 2001). In response to these challenges, a management plan was developed in the early 1990s for both Lake Titicaca and the Altiplano hydrosystem (Revollo, 2001). A key component of this plan was the construction of an outflow gate to regulate lake releases and the establishment of operating rules. Although the outflow gate was completed in 2001, lake releases remain nearly the same as under natural conditions because the operating rules have not yet been implemented. Addressing these water management challenges requires an accurate integrated water balance allowing better knowledge of the drivers of the lake water level variations.

Figure 1 The main geographical features of the Lake Titicaca hydrosystem and the locations of the main streamflow gauges. The reference year for the limits of the glacier is 2000 (RGI Consortium, 2017). The reference year for the croplands is 2010 (Ministerio de Desarrollo Rural y Tierras, 2011; Ministerio del Ambiente, 2015).

Unlike other large lakes, very few studies have been conducted on Lake Titicaca. The only study modeling the water balance of Lake Titicaca and Lake Poopó was the one by Lima-Quispe et al. (2021) using the Water Evaluation and Planning System (WEAP) platform with a monthly time step for the period 1980–2015. The study aimed to distinguish between the relative contributions of climate and irrigation management to water level fluctuations. However, the modeling approach proposed by the authors has a significant limitation because it is based on a scaling factor for precipitation over the lake to close the water balance, which clearly introduces additional uncertainty. Other methodological shortcomings include (i) the omission of snow and ice processes, which can play a non-negligible role in this high-elevation region; (ii) the estimation of evaporation using the Penman method, without accounting for changes in heat storage; and (iii) the use of historical monthly averages (humidity, wind speed, and incoming solar radiation) to calculate reference evapotranspiration and evaporation, without considering interannual variability.

1.4 Scope and objectives

In addressing the challenges and limitations of representing hydrological processes in poorly gauged large lakes such as Lake Titicaca, we pose the following key question: how can a reliable water balance be estimated? To answer this, we present an integrated modeling framework based on conceptual models to estimate the water balance of Lake Titicaca more reliably. The modeling framework is applied at a daily time step for the period 1982–2016, allowing us to represent water level fluctuations over a wide range of hydroclimatic conditions. The specific questions are the following: to what extent are water level variations sensitive to net irrigation withdrawals and to snow and ice processes? What is the role of heat storage change in evaporation from the lake? To address these questions, new approaches are introduced for (i) predicting upstream inflow, including hydrological sensitivity to net irrigation consumption and snow and ice processes; and (ii) estimating evaporation from the lake using the Penman method while accounting for changes in heat storage.

2 Material

2.1 The Lake Titicaca hydrosystem

Lake Titicaca, located at 3812 m a.s.l. on the Altiplano of the tropical Andes of South America, covers an area of approximately 8340 km². The elevation of the catchments that contribute to the lake ranges between 3812 and 6300 m a.s.l. (average 4200 m a.s.l.) and covers an area of approximately 48 780 km². The lake has an average volume of 958 km³ and a maximum depth of 277 m according to the bathymetry carried out between 2016 and 2019 (Autoridad Binacional del Lago Titicaca, 2021). Lake Titicaca is of regional hydrological importance, and the outflows of Lake Titicaca represent up to 79 % of the inflows of Lake Poopó (Lima-Quispe et al., 2021). The Lake Titicaca Authority (Spanish abbreviation ALT) was also created as an autonomous binational entity with the mission of managing the lake.

2.2 Climate data

Daily precipitation and air temperature (see Fig. 2) were obtained from the data generated in Bolivia (Ministerio de Medio Ambiente y Agua, 2018) with the Gridded Meteorological Ensemble Tool (GMET) (Clark and Slater, 2006; Newman et al., 2015). GMET has a spatial resolution of 0.05° for the period 1980–2016. It is based on a probabilistic method using ground station data, with further details provided by Clark and Slater (2006) and Newman et al. (2015). Lima-Quispe et al. (2021) used the same data. In some catchments in the study area, Satgé et al. (2019) evaluated 12 satellite-based products and found that MSWEP and CHIRPS products were the most promising at the 10 d time step. As a result, for daily precipitation, four datasets, i.e., GMET, MSWEP, CHIRPS, and basic interpolation of ground station data with inverse distance weighting (IDW) (Ruelland, 2020), were initially tested according to daily hydrological sensitivity analyses. The results showed that GMET led to more accurate simulations in most catchments and Lake Titicaca (not shown here for the sake of brevity). Daily data on relative humidity, wind speed, and solar radiation are very scarce in space and over time. Many values were missing in the time series from the weather stations, thus calling their quality and representativeness into question. For this reason, we used reanalysis data from ERA5-Land (Muñoz-Sabater et al., 2021) (see Appendix A).

Figure 2 Spatial distribution of (a) annual precipitation and (b) mean temperature for the hydrological period 1980–2016 according to GMET (Gridded Meteorological Ensemble Tool).

Figure 2 shows the spatial pattern of precipitation (see Fig. 2a) and air temperature (see Fig. 2b) based on GMET for the period 1980–2016. Annual precipitation varies between 440 and 1100 mm (mean 725 mm). The wettest areas are concentrated in the northwest and over Lake Titicaca. The driest areas are to the south. The spatial distribution does not show generalized dependence on elevation, but there are small areas on the eastern and western margins where precipitation increases with elevation. The annual mean air temperature varies between −2 and 11 °C (average 6 °C). The coldest areas are the western and eastern areas, coinciding with the highest mountains. The warmest areas are located over Lake Titicaca. The spatial distribution of the air temperature depends on the elevation.

2.3 Snow and glacierized areas

According to MODIS snow cover (Hall et al., 2002) computed over the period 2000–2016 based on a method described in Ruelland (2020), 80 % of the upstream catchments are completely free of snow throughout the year. Areas where the snow cover persists for more than 20 % of the year are located above 4700 m a.s.l. These high-elevation areas represent less than 0.5 % of the total surface area. Glacierized areas are located above 4600 m a.s.l. and represented 231 km² in the early 2000s (RGI Consortium, 2017), i.e., less than 0.5 % of the total area. The estimated glacier volume is ∼ 12 km³ (Farinotti et al., 2019), which represents ∼ 1.3 % of the mean volume of Lake Titicaca (958 km³).

2.4 Croplands and irrigable area

According to land cover maps (Ministerio de Desarrollo Rural y Tierras, 2011; Ministerio del Ambiente, 2015) (see Fig. 1), cropland covers 8069 km² (17 % of the upstream catchments). Agriculture is largely traditional, rainfed, and constrained by droughts, frosts, and hailstorms (Garcia et al., 2007). The main crops include forage grasses, potatoes, grain barley, and quinoa (INTECSA et al., 1993c). Only 40 % of arable land is cultivated due to crop rotation and agroclimatic constraints (INTECSA et al., 1993c). Potatoes are planted in October–November and harvested after 6 months, while quinoa has a similar cycle. Beans and onions, mostly irrigated, are planted from July to September. The land cover maps (Ministerio de Desarrollo Rural y Tierras, 2011; Ministerio del Ambiente, 2015) do not distinguish between rainfed and irrigated areas (see Fig. 1), limiting the identification of changes in irrigated areas over time. The ESA CCI-LC dataset (ESA, 2017) also does not identify irrigated areas in the Altiplano.

The available data come from the inventory of irrigation systems. The inventory on the Bolivian side was created in 2012 (Ministerio de Medio Ambiente y Agua, 2013). For the Peruvian side, the rights of use granted by the Autoridad Nacional del Agua (ANA) until 2023 were used (https://snirh.ana.gob.pe, last access: 12 May 2023). These inventories provide information on location (latitude and longitude), irrigable area, and volume granted (Peru) or a reference volume of irrigation (Bolivia). Figure 3 illustrates the irrigation systems in terms of location and irrigable area. The irrigable area is 767 km² (see Table 1), which represents 1.6 % of the upstream catchments that contribute to Lake Titicaca. Only 9.5 % of the croplands are located in “irrigable areas”, i.e., cropland within the area of influence of an irrigation system that can potentially be irrigated. However, not all of the irrigable area is irrigated because irrigation depends on the availability of water in space and over time. Then, we assumed that the irrigable area was constant over the period 1980–2016. Figure 3 also shows that most of the irrigation systems cover an area of less than 5 ha; i.e., small-scale irrigation predominates. Also, irrigation is mostly practiced by smallholder farmers. Furrow irrigation is the most common system, and its efficiency is about 35 % (Autoridad Nacional del Agua, 2009; Instituto Nacional de Recursos Naturales, 2008; Ministerio de Medio Ambiente y Agua, 2013). The main sources of water are rivers and reservoirs (see Fig. 3). The main reservoir is Lagunillas (see Fig. 3), with a capacity of 500×10⁶ m³. The remaining 15 reservoirs have capacities of less than 30×10⁶ m³. Due to a lack of data on dam management, streamflow regulation was not accounted for, assuming that it has minimal impact on natural flows.

Table 1 Main characteristics of the gauged and ungauged catchments that contribute to the lake. Glacier area was estimated using RGI 6.0 glacier outlines (RGI Consortium, 2017). Cropland area (including rainfed and irrigated area) was estimated using the 2010 land cover maps of Peru and Bolivia (Ministerio de Desarrollo Rural y Tierras, 2011; Ministerio del Ambiente, 2015). Irrigable area was estimated using data from the inventories of agricultural land use rights (Peru) (https://snirh.ana.gob.pe, last access: 12 May 2023) and irrigation systems (Bolivia) (Ministerio de Medio Ambiente y Agua, 2013). The reference year for irrigable area in Bolivia is 2012, and in Peru it has been updated to 2023. Elevations were extracted from a digital elevation model (DEM) at 90 m spatial resolution from the Shuttle Radar Topographic Mission (SRTM; Jarvis et al., 2008).

Figure 3 Locations of irrigation systems in terms of irrigable area (Ministerio de Medio Ambiente y Agua, 2013; https://snirh.ana.gob.pe, last access: 12 May 2023).

2.5 Glaciological and hydrological control data

2.5.1 Geodetic mass balance of the glaciers

Within the study area, no observations are available at the scale of small glacierized catchments. Only geodetic mass balance data are computed at the scale of the entire Lake Titicaca hydrosystem (e.g., Dussaillant et al., 2019; Hugonnet et al., 2021a). The Hugonnet et al. (2021a) dataset, which is based on ASTER satellite stereo imagery, is available for the period 2000–2019. The mass balance at a 5-year time step is provided for RGI 6.0 glacier outlines. The error range of the Hugonnet et al. (2021a) dataset is smaller than the Dussaillant et al. (2019) dataset, and the interpolation of glacier elevation changes is based on Gaussian process regressions. Therefore, we used the Hugonnet et al. (2021a) dataset.

2.5.2 Streamflow records

Seven streamflow gauges (see Table 1 and Fig. 1) with daily records were used in this study. The gauged (ungauged) area represents 76 % (24 %) of the total area of the catchments that feed the lake. The quality of the Peruvian gauge data can be considered satisfactory since monthly streamflow gauging is performed to calibrate the rating curves, but on the Bolivian side the quality of the Escoma and Achacachi data is questionable. According to SENAMHI–Bolivia, the river stages measured in Escoma and Achacachi are prone to systematic errors due to erroneous measurements made by observers with limited measurement training (Escoma) and/or changes in the geomorphology of the riverbed (Achacachi). Streamflow gauging is only carried out twice in a hydrological year.

2.5.3 Lake water levels

We had access to data recorded at two water level gauges: Puno and Huatajata (see Fig. 1). The Puno gauge (also known as Muelle ENAFER) is managed by SENAMHI–Peru, while the Huatajata gauge is managed by SENAMHI–Bolivia. The daily historical water levels from Puno are more reliable. In the case of Huatajata, inconsistencies were detected in the records made prior to 1998. Additionally, during a field visit, it was observed that the Huatajata measurement scale is prone to displacement. Therefore, in this study we used data from Puno, which provide continuous daily water levels (m a.s.l.) over the period 1982–2016.

3 Method

3.1 Modeling framework used to quantify the water balance in a high mountain lake hydrosystem

The water balance of Lake Titicaca was modeled at a daily time step for the hydrological period from 1 September 1981 to 31 August 2016 (hereafter 1982–2016) using a lumped model store following the equation

$$\begin{array}{}\text{(1)}& P_{\mathrm{lake}}+Q_{\mathrm{in}}-E_{\mathrm{lake}}-Q_{\mathrm{out}}={\displaystyle \frac{\mathrm{d}h}{\mathrm{d}t}}\pm Q_{\mathrm{gw}}\pm \mathit{\epsilon },\end{array}$$

where P_lake, Q_in, E_lake, and Q_out are, respectively, direct precipitation over the lake, inflow from upstream catchments, evaporation from the lake, and downstream outflow. The term $\mathrm{d}h/\mathrm{d}t$ represents the storage change in the lake over a time window. Q_gw represents the net groundwater exchange, and ε represents the error that cannot be explained by the components of the water balance. The unit of the water balance terms is millimeters per day.

WEAP (Yates et al., 2005) was adapted and used to represent the water balance dynamics. The models in WEAP typically seek a compromise between data availability and the complex representation of hydrological processes. This is essential in the context of poorly gauged regions, where it is not possible to represent all hydrological processes in sufficient detail. Unlike the study by Lima-Quispe et al. (2021), which also uses WEAP, this study uses a daily time step and a different approach to simulating irrigation water allocation. We acknowledge that there is an overlap in the precipitation and air temperature data as well as in the irrigable area on the Bolivian side. However, in order to implement the model at a daily time step, it was necessary to collect new data updated to the required timescale. In addition, new data were available, such as lake bathymetry and irrigable area on the Peruvian side. Figure 4 shows the main processes and parameters of the modeling chain used.

Figure 4 Main processes and parameters (in red) of the modeling chain used to simulate the water balance of Lake Titicaca. P, T, ETo, Pe, P_lake, Q_in, E_lake, Q_gw, Q_out, and h stand for, respectively, precipitation, air temperature, reference evapotranspiration, effective precipitation, direct precipitation over the lake, upstream catchment inflow, evaporation from the lake, net groundwater exchange, downstream outflow, and lake storage. Root zone and deep zone stores were modified based on Yates et al. (2005).

3.1.1 Direct precipitation over the lake

P_lake was extracted from GMET for the outline of the lake. The large area and volume of Lake Titicaca favor absorption of solar radiation and result in higher water temperatures than in the surrounding areas, which, in turn, induces convection and higher precipitation over the center of the lake (Roche et al., 1992). However, the magnitude and spatial distribution of precipitation over the lake are not understood well. GMET included two precipitation gauges located on two different islands in the lake. In this study, the extracted data were used directly without correcting for scaling factors. Significant underestimation of precipitation could lead to significant error.

3.1.2 Upstream inflow

Q_in was estimated using a conceptual modeling approach that combines a degree-day model for simulating snow and ice processes with the Soil Moisture Model (SMM, part of WEAP) (Yates et al., 2005) to simulate the processes contributing to the generation and regulation of water storage and water flow in the catchments, including irrigation. The model was applied using the same 100 m elevation bands in each catchment to account for snow and ice accumulation and melt, and glacierized and non-glacierized areas were distinguished in each elevation band.

For the snow and ice processes, a degree-day model was applied that considered two stores: one for ice and one for snow (see Fig. 4) in a semi-distributed mode with 100 m elevation bands. However, each glacier in each catchment was simulated separately. For snow accumulation, the liquid (P_rain, mm) and solid (P_snow, mm) fractions of total precipitation were estimated from a linear separation between the snow (T_s) and rain (T_l) temperature thresholds according to values (see Table 2) recommended in Ruelland (2023). Potential (maximum) snowmelt was calculated as DDF_snow(T−T_m), where DDF_snow is the degree-day factor (mm d⁻¹ °C⁻¹), T is the air temperature (°C), and T_m is the melting temperature threshold (°C). T_m was calibrated according to two values (T_m,max and T_m,min), where the maximum value occurs in austral summer and the minimum value in winter. In the outer tropical regions, the amplitude of the diurnal range of air temperature is indeed considerable in winter. This means that it is warmer during the daytime, which can increase melting, while cold conditions prevail at night. The seasonal variation of T_m was calculated using the following equation:

$$\begin{array}{}\text{(2)}& T_{\mathrm{m}}={\displaystyle \frac{T_{\mathrm{m},\max }+T_{\mathrm{m},\min }}{\mathrm{2}}}+{\displaystyle \frac{T_{\mathrm{m},\max }-T_{\mathrm{m},\min }}{\mathrm{2}}}\sin \left(\mathrm{2}\mathit{\pi }{\displaystyle \frac{D+\mathrm{81}}{\mathrm{365}}}\right),\end{array}$$

where T_m,max and T_m,min are the maximum and minimum temperature thresholds (°C) and D is the Julian day. The maximum value of T_m was assumed to occur on 10 January and the minimum on 12 July.

Table 2 Parameters of the upstream catchment model of Lake Titicaca and the associated fixed values or ranges tested. The ranges presented were used to generate the random sample of hypercubes in the Monte Carlo approach.

Actual snowmelt (M_snow) was determined as a function of maximum snowmelt and snow accumulation. For ice melt (M_ice), the same approach was used as for potential snowmelt, except that DDF_snow was replaced by an ice degree-day factor (DDF_ice). Ice melts when it is not covered by a snowpack. The daily mass balance (B, mm w.e.) and effective precipitation in glacierized areas (Pe_g, mm w.e.) in each elevation band (j) were computed as follows:

$$\begin{array}{}\text{(3)}& {\displaystyle }{B}_j=P_{\mathrm{snow},j}-M_{\mathrm{snow},j}-M_{\mathrm{ice},j},\text{(4)}& {\displaystyle }{\mathrm{Pe}}_{\mathrm{g},j}=P_{\mathrm{rain},j}+M_{\mathrm{snow},j}+M_{\mathrm{ice},j}.\end{array}$$

The annual mass balance in each evaluation band, B_a,j, was estimated as the sum of the daily mass balance in a hydrological year. The annual mass balance was also calculated for individual glaciers (B_a,g) for comparison with the available glaciological and geodetic mass balance. B_a,g was calculated as

$$\begin{array}{}\text{(5)}& B_{\mathrm{a},\mathrm{g}}={\displaystyle \frac{\sum _{j=\mathrm{1}}^n(B_{a,j}\times A_{\mathrm{g},j})}{A_{\mathrm{g}}}},\end{array}$$

where A_g,j and A_g are the glacier area in the elevation band j and the total glacier area, respectively.

The glacierized surface area (RGI Consortium, 2017) was fixed for the period simulated. The area provided for the year 2000 was considered an intermediate value for the period 1982–2016. Ice thickness was also assumed to be infinite. The air temperature in each elevation band (T_j) was estimated as $T_{\mathrm{GMET},j}+\mathrm{\Gamma }(Z_{\mathrm{GMET}}-Z_j$), where T_GMET,j is the air temperature derived from GMET for each elevation band, Z_GMET is the mean areal elevation signal from GMET in the elevation band j, Z_j is the mean elevation of the elevation band, and Γ is a constant temperature lapse rate that was set to the value calculated from GMET (i.e., 5.8°C km⁻¹). The precipitation extracted from GMET for each elevation band was used directly with no modification.

Regarding SMM, it is a one-dimensional model based on two stores (see Fig. 4). The first store represents the root zone and the second the deep zone (Yates et al., 2005). The model without irrigation has seven free parameters as shown in Fig. 4, of which the crop coefficient (K_c) can be set using reference values from the literature. In addition, there are two parameters associated with the initial states of the two stores called z₁ and z₂. An additional parameter was included for runoff routing using the Muskingum equation. SMM is driven by precipitation and reference evapotranspiration estimated by the modified Penman–Monteith method (Maidment, 1993) for a grass crop 0.12 m in height and with a surface resistance of 69 s m⁻¹. The climate input data are detailed in Sect. 2.2. The effective precipitation in the elevation band j of both the non-glacierized and glacierized fractions is given as

$$\begin{array}{}\text{(6)}& {\mathrm{Pe}}_j=(P_{\mathrm{rain},\mathrm{ng},j}+M_{\mathrm{snow},\mathrm{ng},j})A_{\mathrm{ng},j}+{\mathrm{Pe}}_{\mathrm{g},j}(\mathrm{1}-A_{\mathrm{ng},j}),\end{array}$$

where P_rain,ng and M_snow,ng refer to rainfall and snowmelt in the non-glacierized fraction (mm). The term A_ng is the relative area of the non-glacierized fraction.

In SMM, water requirements (WRs) for irrigation are determined by crop evapotranspiration (from seasonal crop coefficients (K_c) and reference evapotranspiration) and the depletion of available water in the root zone store (see Fig. 4). K_c adjusts the reference evapotranspiration to reflect crop-specific characteristics (Allen et al., 1998), such as phenology. It was derived using cropping calendar and crop type data (Autoridad Nacional del Agua, 2009, 2010; Instituto Nacional de Estadística, 2015). The lower and upper irrigation threshold parameters (L_irr and U_irr; see Table 2 and Fig. 4) dictate both the timing and quantity of water used for irrigation (Yates et al., 2005). When the relative soil moisture of the root zone store drops below the lower threshold, a water requirement is triggered and irrigation is supposed to be applied up to the upper threshold (Yates et al., 2005). The irrigation use of runoff (IUR) method was used to allocate water. This method consists of setting or calibrating a percentage (IUR) of a catchment's runoff (before the runoff reaches the main river) that can be used for internal irrigation. IUR focuses on water allocation at the catchment scale, especially when hundreds of irrigation systems are to be represented together. Simulating each of the irrigation systems shown in Fig. 3 individually would not be feasible at the scale of this study, and the IUR approach is thus better suited for this. The irrigation net consumption, IRR_net, was calculated as

$$\begin{array}{}\text{(7)}& {\mathrm{IRR}}_{\mathrm{net}}=min\left(Q_{\mathrm{wi}}\times \mathrm{IUR},\mathrm{WR}\right)\times (\mathrm{1}-{\mathrm{IRR}}_{\mathrm{rf}}),\end{array}$$

where Q_wi is runoff without irrigation (mm), IUR is a calibration parameter expressed as a percentage, WR is the irrigation water requirement (mm), and IRR_rf is the irrigation runoff fraction expressed as a percentage. The term $min\left(Q_{\mathrm{wi}}\times \mathrm{IUR},\mathrm{WR}\right)$ is the water withdrawn for irrigation. IRR_rf is calculated as follows: (i) in the first iteration, SMM simulates Q_wi; (ii) in the second iteration, runoff is simulated assuming that the full WR is supplied; and, (iii) finally, IRR_rf is estimated based on how much runoff would flow due to irrigation alone.

3.1.3 Evaporation from the lake

E_lake was estimated using the Penman method for open water (Penman, 1948). This method is justified because it requires fewer data than an energy balance approach but is not as simple as a temperature-based approach. The Penman method also attempts to incorporate the energy balance in a simplified manner and includes mass transfer. The equation is given as

$$\begin{array}{}\text{(8)}& E_{\mathrm{lake}}={\displaystyle \frac{\mathrm{\Delta }}{\mathrm{\Delta }+\mathit{\gamma }}}{\displaystyle \frac{R_{\mathrm{n}}-G}{\mathit{\lambda }}}+{\displaystyle \frac{\mathit{\gamma }}{\mathrm{\Delta }+\mathit{\gamma }}}f\left(U_{\mathrm{2}}\right)(e_{\mathrm{s}}-e_{\mathrm{a}}),\end{array}$$

where E_lake is the evaporation (mm d⁻¹), Δ is the slope of the vapor pressure curve (kPa °C⁻¹), λ is the latent heat vaporization set at 2.45 MJ kg⁻¹, and γ is the psychrometric constant (kPa °C⁻¹). The terms R_n and G are net radiation at the water surface and heat storage changes (MJ m⁻² d⁻¹), respectively. f(U₂) is the function of wind speed measured at 2 m above the lake surface that is equal to c(a+bU₂), where the constants a=10, b=5.4, and c=0.26 for Lake Titicaca were taken from Delclaux et al. (2007). Also, e_s is the vapor pressure at the evaporating surface (kPa), and e_a is the vapor pressure at 2 m above the lake surface (kPa). Δ is given as $(e_{\mathrm{s}}-e_{\mathrm{a}})/(T_w-T)$, where T_w and T are the evaporating surface temperature and air temperature (°C), respectively.

R_n is the sum of net shortwave radiation (K) and net longwave radiation (L). K is given as K_in(1−α), where K_in is the incident solar radiation (MJ m⁻² d⁻¹) and α is the albedo of the water surface. The L component is the difference between the incident flux (L_in) emitted by the atmosphere and clouds and outgoing radiation (L_out) from the evaporating surface. L_in and L_out can be estimated with Eqs. (9) and (10), respectively. For L_in, we used the equation calibrated by Sicart et al. (2010) on the Zongo glacier, which is located at a distance of about 100 km from Lake Titicaca. The authors suggest that the calibration can be used in the tropical Andes.

$$\begin{array}{}\text{(9)}& {\displaystyle }{L}_{\mathrm{in}}=C{\left({\displaystyle \frac{e_{\mathrm{a}}}{T+\mathrm{273.15}}}\right)}^{\mathrm{1}/\mathrm{7}}(\mathrm{1.67}-{\mathit{\tau }}_{\mathrm{atm}}.83)\mathit{\sigma }(T+\mathrm{273.15}{)}^{\mathrm{4}},\text{(10)}& {\displaystyle }{L}_{out}={\mathit{\epsilon }}_{\mathrm{w}}\mathit{\sigma }(T_w+\mathrm{273.15}{)}^{\mathrm{4}},\end{array}$$

where, for a daily time step, C is equal to 1.24; e_a is the vapor pressure (hPa); σ is the Stefan–Boltzmann constant (set at $\mathrm{4.90}\times {\mathrm{10}}^{-\mathrm{9}}$ MJ m⁻² K⁻⁴ d⁻¹); τ_atm is the atmospheric transmissivity (–) that can be approximated as $K_{\mathrm{in}}/S_{\mathrm{extra}}$; S_extra is the theoretical shortwave irradiance (MJ m⁻² d⁻¹) at the top of the atmosphere; and ε_w is the water emissivity set at 0.98 (–).

The term G (heat storage change) in Eq. (8) was estimated using the approach and assumptions of the study conducted by Pillco Zolá et al. (2019) on Lake Titicaca.

$$\begin{array}{}\text{(11)}& G=c_{\mathrm{w}}{\mathit{\rho }}_{\mathrm{w}}{\displaystyle \frac{V_{\mathrm{mix}}}{A_{\mathrm{lake}}}}{\displaystyle \frac{\mathrm{d}{T}_w}{\mathrm{d}t}},\end{array}$$

where c_w is the specific heat of water ($\mathrm{4.18}\times {\mathrm{10}}^{-\mathrm{3}}$ MJ kg⁻¹ °C⁻¹), ρ_w is the water density (1000 kg m⁻³), V_mix is the volume above the mixing depth (m³), and A_lake is the surface area of the lake (m²). d$T_{\mathrm{w}}/\mathrm{d}t$ is the change in water temperature (°C) over the time interval (d).

Air temperature was obtained from GMET, while all the other meteorological variables were obtained from ERA5-Land (Sect. 2.2). Since there are no long-term measurements of lake surface water temperature (LSWT) and the remotely sensed datasets do not cover the entire study period, the Air2Water model (Piccolroaz et al., 2013; Toffolon et al., 2014) was used to simulate LSWT. Calibration and evaluation were performed against ARC-Lake V3 remotely sensed data (MacCallum and Merchant, 2012) (see Appendix B).

3.1.4 Downstream outflow

Q_out was simulated using the rating curve shown in Fig. 5a. This curve was established 30 years ago based on a hydrodynamic simulation of the Desaguadero River (INTECSA et al., 1993b). The elevation corresponds to the vertical datum of Peru. The rating curve was used for the entire study period and implemented in the model in combination with the lake bathymetry (see Fig. 5b) carried out between 2016 and 2019 by ALT.

Figure 5 Information used to simulate Q_out for (a) the rating curve and (b) the bathymetry. The rating curve of the lake outlet to the Desaguadero River was obtained from the master plan (INTECSA et al., 1993b). In both figures the elevation is referenced to the Peruvian vertical datum. The bathymetry (i.e., the relationship between lake water level and storage volume) was recorded by ALT between 2016 and 2019.

3.1.5 Storage change

$\mathrm{d}h/\mathrm{d}t$ was calculated directly from the water levels measured at the Puno gauge. Storage change is basically the difference between the water levels of the current time step and the previous time step.

3.1.6 Net groundwater exchange

Q_gw was considered negligible. According to INTECSA et al. (1993a), the leakage from Lake Titicaca to the aquifers is very limited, and the lake can be considered an almost completely closed surface system. This is because the lake bed is composed of sediments with very low permeability. In this case, the only areas of high permeability would be limited to alluvial deposits saturated by water that mostly flows towards the lake. According to the same study, the inputs from alluvial deposits were 0.56 m³ s⁻¹. Therefore, omitting Q_gw from the water balance is justified.

3.2 Evaluation of the modeling framework

The performance of the lake water balance model was evaluated using both the error term of Eq. (1) and the RMSE computed between observed and simulated water levels. Since the Q_in and E_lake were modeled, intermediate calibration and evaluation were necessary. Evaporation measurements are not available in the study area to evaluate the performance of the E_lake. Therefore, a formal calibration and evaluation procedure was only implemented for the Q_in estimates. The procedure was applied sequentially, first to obtain the model parameters simulating snow and ice processes (see Appendix C) and then to calibrate and evaluate the upstream catchment model using the Nash–Sutcliffe efficiency index (Nash and Sutcliffe, 1970) calculated on the root-mean-square-transformed streamflow (NSE_sqrt) in order to provide an intermediate fit between high and low flows (Oudin et al., 2006).

3.2.1 Upstream catchment model

The upstream catchment model shown in Fig. 4 has 15 parameters. However, seven of the parameters were set to reduce the number of free parameters. Four parameters (T_m,max, T_m,min, DDF_snow, and DDF_ice) related to snow and ice store were set to values obtained in the Zongo catchment (see Appendix C). The simulated mass balance of all glaciers of the Titicaca hydrosystem was compared with the geodetic mass balance (Hugonnet et al., 2021a) for the period 2000–2009. Similarly, two parameters of the irrigation module, L_irr and U_irr, were set to 80 %. Winter et al. (2017) used a threshold of 100 % for furrow irrigation in California. However, a value of 80 % is reasonable for our study area because it is irrigated under conditions of limited water availability. X (routing store) was set to a default value (0.2) in WEAP. A total of eight free parameters were kept, as shown in Table 2. The procedure used to obtain the set of parameters with the best performance and the subsequent evaluation consisted of four steps.

First (Step 1 in Fig. 6), the model was run for the period 1980–2016 (of which the first 2 years were used as a warmup period) with 10 000 parameter sets generated from a random sample of hypercubes from the Monte Carlo approach within the parameter intervals tested (see Table 2). Second (Step 2 in Fig. 6), the best-performing parameter sets in terms of NSE_sqrt were identified along with subperiods consisting of (i) 5 continuous years (i.e., seven subperiods between 1982 and 2016) and (ii) 5 discontinuous years identified as coldest, warmest, driest, and wettest (i.e., four discontinuous subperiods between 1982 and 2016). The 11 best-performing parameter sets were selected using different subperiods. Third (Step 3 in Fig. 6), the mean of 11 streamflow simulations generated with the selected parameter sets was calculated. Fourth (Step 4 in Fig. 6), the performance of the mean of the streamflow simulations over the 11 subperiods was evaluated using the NSE_sqrt criterion.

Figure 6 Procedure used to evaluate the upstream catchment model to estimate the Q_in in each gauged catchment.

Two objectives justify the use of seven continuous and four discontinuous subperiods. The first objective was to evaluate the transferability of the model parameters to nonstationary conditions within the period 1982–2016, including particularly contrasted subperiods in terms of precipitation and temperature. In addition, continuous subperiods were suitable for assessing the transferability of the water allocation parameter (IUR). If there had been a significant and sustained increase in the irrigation withdrawals, the IUR parameter would not have been transferable over time and would therefore have varied over the period 1982–2016. However, it is important to note that irrigable area, crop type, and K_c remain under the assumption of stationarity. The second objective was to consider the parameterization uncertainty through uncertainty envelopes and confidence intervals in the simulations resulting from the 11 parameter sets used over the whole 1982–2016 period.

3.2.2 Sensitivity assessment of model predictability to irrigation and snow and ice processes

To evaluate the sensitivity of the modeling chain to net irrigation consumption and snow and ice processes, different processes were progressively excluded from the modeling structure. Three model structures were tested in the following configuration: (i) BasicModel+IRR+SNOW/ICE, which represents the full reference model structure (as shown in Fig. 4); (ii) BasicModel+IRR, where the processes associated with snow and ice (accumulation and ablation) are excluded but irrigation is maintained; and (iii) BasicModel, where snow and ice processes and irrigation consumption are excluded. The objective was to evaluate whether a simpler structure in terms of the hydrological processes considered performed better or worse than a more complex structure in simulations of catchment streamflow and lake water levels. This is justified because, a priori, snowfall and the proportion of glacierized area and irrigable area are very limited in the catchments that contribute to the lake (see Table 1), and their impact on the streamflow and water level prediction may be negligible. In addition, the data used to represent and control these processes are very limited, which can lead to inaccuracies that could worsen the simulations of streamflow and lake water levels instead of improving them. Table 2 shows the active and inactive parameters in the three model structures.

3.2.3 Transferring parameters to the ungauged catchments

The approach used to transfer the parameters to the ungauged catchments consisted of the following steps: (i) the median of the parameter sets obtained for the seven gauged catchments was calculated for each subperiod, thus generating 11 parameter sets; (ii) the upstream catchment model was run for the 11 parameter sets; and (iii) the mean of the 11 streamflow simulations was calculated, including the confidence interval.

4 Results

4.1 Modeling chain performance

4.1.1 Performance of the geodetic mass balance simulated by the model for Titicaca glaciers

Figure 7 shows the comparison between the simulated and geodetic mass balances. The scatterplots reveal significant variability in model performance, with some glaciers (represented by each point) close to the identity line and others deviating significantly. The model simulates a more negative glacier mass balance compared to the geodetic glacier mass balance. Figure 7a displays glaciers according to their surface area, while Fig. 7b shows them based on their mean elevation. The model performed more effectively in the catchments that concentrate 92 % of the glacierized area (Achacachi, ungauged catchments, Tambillo, and Escoma). The model performed very poorly in Ramis (see Table 3), but the proportion of glacierized area in that catchment is very small (0.1 %). Therefore, the proportion is expected to have a limited impact on streamflow predictions. Additionally, the model performs much better in the case of large glaciers than small glaciers (see Table 3). For example, for large glaciers representing 84 % of the glacierized area, the weighted simulated mass balance for all the catchments has a bias of 11 %, which is within the error range of the geodetic mass balance (see Table 3). Much more obvious is the dependence on elevation (see Fig. 7b). The model performed very poorly on glaciers with mean elevations below 5100 m a.s.l., which represent only 10 % of the simulated glaciers. At elevations above 5200 m a.s.l., i.e., 68 % of the glaciers, the points are distributed around the identity line. The notable variability in model performance (see Fig. 7b) could be attributed to inaccuracies in precipitation data for some catchments, as estimating precipitation in high-elevation remote areas remains a complex challenge (Ruelland, 2020). At the catchment scale, the underestimation and overestimation of glacier-wide mass balance are compensated for, and the biases are relatively small (see Table 3).

Figure 7 Scatterplots comparing simulated and geodetic glacier mass balances for 2000–2009, based on the remotely sensed observation of Hugonnet et al. (2021a). The dot size represents the (a) glacier area and (b) mean elevation. The dashed line indicates the identity line, while the gray line represents the error in the geodetic glacier mass balance.

Table 3 Comparison of simulated and geodetic glacier mass balances (MBs) in upstream catchments for the period 2000–2009. The catchment-scale mass balance was calculated as an area-weighted average of each glacier. The error corresponds to the geodetic glacier mass balance and was calculated as a weighted average for each catchment.

4.1.2 Performance and sensitivity of the upstream catchment model to irrigation and snow and ice processes

Figure 8 shows the distribution of the performance of the three modeling options (BasicModel+IRR+SNOW/ICE, BasicModel+IRR, and BasicModel) in the evaluation made for each catchment. A striking feature is that the performance in the catchments on the Peruvian side (Ramis, Ilave, Coata Unocolla, and Huancane) was significantly better than in the catchments on the Bolivian side (Escoma, Tambillo, and Achacachi). The performance of the three models is distributed symmetrically in the boxplots in each catchment. This suggests that there are no significant performance differences between the three models. Therefore, Q_in is not very sensitive to net irrigation consumption or snow and ice processes. This reflects the fact that the upstream catchments of Lake Titicaca are dominated by a pluvial monsoonal climate regime and that contributions from glacierized areas and snow have little influence on the Q_in prediction. However, in terms of mean and median NSE_sqrt, BasicModel+IRR+SNOW/ICE and BasicModel+IRR generally performed slightly better than BasicModel, but the differences are marginal: in most cases the difference is between 1 % and 3 %. However, in Achacachi, the NSE_sqrt obtained with BasicModel+IRR+SNOW/ICE is 6 % higher than the value obtained with the BasicModel+IRR and BasicModel models.

Figure 8 Performance distributions of the three modeling options (BasicModel+IRR+SNOW/ICE, BasicModel+IRR, and BasicModel) in each gauged catchment in the evaluation. The size of the sample in each boxplot is 11 (resulting from the procedure presented in Fig. 6).

4.1.3 Performance of the modeling chain with respect to lake water levels

Figure 9 shows the simulation of daily water levels in Lake Titicaca over the 1982–2016 period. Three different estimates of Q_in were used, one from BasicModel+IRR+SNOW/ICE, one from BasicModel+IRR, and one from BasicModel. The results show that the models are able to simulate the amplitude and frequency of annual, interannual, and decadal water level fluctuations reasonably well. Visually, the performances of the three Q_in estimates appear to be relatively similar. Based on the ε term, BasicModel+IRR performed better than BasicModel+IRR+SNOW/ICE and BasicModel. The difference in the ε term between BasicModel+IRR+SNOW/ICE and BasicModel+IRR was 0.01 mm. When snow and ice processes and net irrigation consumption were excluded, the error increased by 0.03 mm d⁻¹. However, the differences in the error were marginal. Figure 9 shows that, in the mostly dry years of the 1990s, BasicModel+IRR+SNOW/ICE and BasicModel+IRR simulated daily water levels better. However, the performance measured by the RMSE differed in the error term. In that case, BasicModel performed best. This is because the RMSE was calculated using daily water level data (cumulative change in storage over time), whereas the ε term was computed directly from the water balance at each time step. Inaccuracies in some time steps were then propagated to later time steps due to the slow response of the lake. The RMSE obtained was very small compared to the average water level (3809.7 m a.s.l.).

Figure 9 Performance of the modeling chain compared with the water levels of Lake Titicaca measured in Puno. Modeling was evaluated using the error (ε) term of the lake water balance and RMSE. Three simulations of water levels are presented because the models produced three different estimates of Q_in.

A striking feature of Fig. 9 is the systematic underestimation of daily water levels between 2001 and 2010. This is related to inaccuracies in the estimation of water balance terms. The hydrological response of Lake Titicaca is relatively slow, and it was possible to verify that significant errors were present between 2001 and 2004 but were not underestimated over the whole decade (2001–2010) (see Fig. D1). In 2001, the outflow gate to regulate outflows into the Desaguadero River was completed, and dredging of the Desaguadero River began. This could mean fewer outflows into the Desaguadero River and consequently more storage in the lake. However, even assuming that there were no downstream outflows in those years would not compensate for the underestimation. Looking for other sources of error, discussions with ALT revealed that the lake sometimes received inflows from the Desaguadero River, especially in wet years. This could be the case, since 2001 and 2004 were wet hydrological years, but other than verbal communication there are no records to support such a claim.

BasicModel+IRR+SNOW/ICE produced a reasonably realistic simulation with an even smaller error, with marginal differences, than BasicModel. Consequently, BasicModel+IRR+SNOW/ICE can be considered a reference structure for providing an estimate of the water balance of the Lake Titicaca hydrosystem.

4.2 Simulated water balance

4.2.1 Simulated water balance in upstream catchments

Table 4 shows the mean annual water balance for the hydrological period 1982–2016 simulated by BasicModel+IRR+SNOW/ICE. Some terms (AE_ni, IRR, and Q_in) were calculated as the average of 11 ensemble members. At the scale of all the upstream catchments (see Table 4), the annual precipitation is 723 mm, of which 6 % is estimated to be snowfall. Thus, despite the elevations at which the study area is located (> 3810 m a.s.l.), the precipitation regime is clearly dominated by rainfall. This regime has a unimodal pattern with a strong seasonal cycle between the rainy and dry seasons (see Fig. 10). This is very characteristic of the outer tropical areas, where precipitation mainly occurs in the summer and the dry season in winter. In the upstream catchments, snowmelt accounts for 6 % of the total input ($P_{\mathrm{rain}}+M_{\mathrm{snow}}+M_{\mathrm{ice}}$). Snow remains on the ground for only a very short time, as there is practically no time lag between snowfall and snowmelt (see Fig. 10). Snowfall and snowmelt have a less accentuated seasonal cycle, as the rare precipitation that occurs in the dry season is mostly snow. The contribution of ice melt was simulated to be 6 mm yr⁻¹ on average, which represents only 1 % of the total input in upstream catchments. However, the contribution of ice melt is slightly higher in the Achacachi and Escoma catchments, where the simulated value is 7 %. Most ice melt occurs between October and November (see Fig. 10) because the temperature is above the melting threshold and glaciers are mostly free of snow cover. Δ_ice is negative, indicating a loss in ice stock.

Table 4 Mean annual water balance (mm) in the upstream catchments for the hydrological period 1982–2016 simulated with BasicModel+IRR+SNOW/ICE. P, P_rain, P_snow, M_snow, M_ice, AE_ni, IRR, Δ_ice, Δ_snow, Δ_SM, and Q_in represent total precipitation, rainfall, snowfall, snowmelt, ice melt, actual evapotranspiration in non-irrigated areas, net consumption of irrigation water, variation in ice storage, variation in snow storage, variation in soil moisture storage, and streamflow in upstream catchments, respectively. The water balance follows the equation $Q_{\mathrm{in}}=P-{\mathrm{AE}}_{\mathrm{ni}}-\mathrm{IRR}-{\mathrm{\Delta }}_{\mathrm{ice}}-{\mathrm{\Delta }}_{\mathrm{snow}}-{\mathrm{\Delta }}_{\mathrm{SM}}$.

Figure 10 Seasonal cycle (monthly average for the period 1982–2016) of the water balance in the upstream catchments of Lake Titicaca simulated with BasicModel+IRR+SNOW/ICE. P, P_snow, M_snow, M_ice, AE, IRR, and Q_in represent total precipitation, snowfall, snowmelt, ice melt, actual evapotranspiration, irrigation net water consumption, and streamflow in the upstream catchments. For some terms (AE, IRR, and Q_in), the gray bars were estimated as the mean of the 11 ensemble members resulting from the procedure shown in Fig. 6. The predictive uncertainty is presented for both the entire prediction range (i.e., all predictive uncertainty) and the 95 % confidence interval. “All predictive uncertainty” was estimated as the maximum and minimum values of the ensemble members. The terms associated with snow and ice are not subject to predictive uncertainty because fixed parameters were used (see Table 2).

Annual actual evapotranspiration (AE) was simulated at 570 mm, of which 2 % corresponds to net irrigation water consumption (IRR). AE is highest between January and March and lowest between July and September (see Fig. 10). IRR is concentrated in the transition season (see Fig. 10). The simulated streamflow is about 153 mm, which represents 21 % of the total outflow in upstream catchments (AE+Q_in), the remaining 79 % being actual evapotranspiration. IRR represents 7 % of Q_in. The seasonal cycle of Q_in shows that the peak is in February, while Q_in is low from June to November. The predictive uncertainty associated with the upstream catchment model parameters is shown in Fig. 10 for the terms AE, IRR, and Q_in. The prediction range of AE is very narrow (see Fig. 10) compared to that of IRR and Q_in, indicating that evapotranspiration is less sensitive to the model parameters. The range of prediction of IRR is quite wide (see Fig. 10), indicating that the simulations are very sensitive to SMM parameters. The range of the prediction of Q_in is also wide (see Fig. 10). However, the average predictions fitted the measured streamflow and model performance reasonably well (see Fig. 8).

4.2.2 Simulated water balance in the lake

Figure 11a shows the annual evolution of the water balance, and Table 5 shows the long-term average values. Over the period 1982–2016, the average annual precipitation over the lake was 744 mm (σ=144 mm) and the inflow from the upstream catchments 958 mm (σ= 392 mm). This means that 44 % of the inflows come from direct precipitation over the lake, while the remaining portion (56 %) comes from upstream catchments. Regarding outflows, the annual evaporation from the lake is 1616 mm (σ= 28 mm) and the downstream outflow is 121 mm (σ= 191 mm). Thus, 93 % of the losses are due to evaporation and 7 % due to downstream outflow. The measured storage change for the period 1982–2016 was −50 mm, indicating a drop in the water level. The simulated change in storage was −35 mm, which indicates an overestimation. Therefore, the water balance closure has an error of about −15 mm. Compared to evaporation, P_lake and Q_in showed significant interannual variability (see Fig. 11a). E_lake was subject to less pronounced interannual variability but showed an increasing trend over the period 1982–2016 due to the increase in temperature (about +0.1 °C per decade). This means that the interannual variability of water levels depends to a large extent on variations in precipitation over the lake and in the upstream catchments. The highest values of both P_lake and Q_in occurred between 1985 and 1986 and caused large floods around the lake and in the Desaguadero River, where discharges reached about 900 mm yr⁻¹. In the 1990s, inflows were lowest in the period studied, resulting in a mostly negative change in storage. Substantial inflows to the lake in the early 2000s led to a significant positive change in storage, although this was followed by another dry period that lasted until 2012.

Table 5 Lake Titicaca water balance components simulated for the period 1982–2016. The lake water balance follows the equation $P_{\mathrm{lake}}+Q_{\mathrm{in}}-E_{\mathrm{lake}}-Q_{\mathrm{out}}-\mathrm{d}h/\mathrm{d}t+\mathit{\epsilon }=\mathrm{0}$.

Figure 11 Water balance of Lake Titicaca for the period 1982–2016 in terms of (a) interannual variability and (b) seasonal cycle. The values in parentheses correspond to the mean annual or monthly values for the period 1982–2016. The lake water balance follows the equation $P_{\mathrm{lake}}+Q_{\mathrm{in}}-E_{\mathrm{lake}}-Q_{\mathrm{out}}-\mathrm{d}h/\mathrm{d}t+\mathit{\epsilon }=\mathrm{0}$.

Figure 11b shows a very marked seasonal cycle for precipitation and upstream inflow. For example, the monthly peak in upstream inflow is about 230 mm (in February), while in the dry season the values are very close to 15 mm (in September). One of the features is the lag of upstream inflow with respect to direct precipitation over the lake, which is evidence of the relatively slow hydrological response in the upstream catchments because of the size of the catchment area. The seasonal cycle of evaporation is less marked than that of precipitation and upstream inflow. Although air temperature shows strong seasonality, evaporation is also influenced by other meteorological variables, and heat storage plays a critical role in the seasonal cycle. The peak of the mean monthly evaporation is around 170 mm (in January), while the minimum value is around 95 mm (in August). In Fig. 11b, it is also interesting to observe the seasonal cycle of the error term, which is mostly positive in the rainy season and mostly negative in the dry season. This could indicate that the lake receives net groundwater inflow during the rainy season and is subject to net outflow during the dry season. However, the magnitudes cannot be attributed directly to net groundwater flow because the error term includes both the uncertainty associated with estimating the other terms of the water balance (P_lake, Q_in, E_lake, and Q_out) and assumptions concerning the groundwater flow.

5 Discussion and conclusions

5.1 Main findings

This study presents three main novelties. First, our integrated modeling framework accurately simulates the daily water balance of Lake Titicaca without requiring scaling factors. Consequently, the propagation of uncertainty in estimating components of the water balance is significantly reduced. For instance, Fig. 12 illustrates how omitting the calibrated precipitation scaling factor used by Lima-Quispe et al. (2021) leads to unrealistic simulations of lake water levels. Our modeling approach also benefits from (i) a rigorous calibration and evaluation procedure for simulating upstream inflows (see Fig. 6), (ii) the estimation of evaporation from the lake using the Penman method while accounting for LSWT and heat storage, and (iii) estimates of reference evapotranspiration and lake evaporation that account for climate variability, using the ERA5-Land dataset for humidity, wind speed, and solar radiation. This study provides a realistic water balance that estimates most hydrological processes (see Tables 4 and 5), although the role of groundwater remains a major unknown. Its magnitude is expected to be a small component of the total water balance.

Figure 12 Comparison of simulated water levels by BasicModel+IRR+SNOW/ICE (aggregated to the monthly time step) with those obtained by Lima-Quispe et al. (2021) without applying the scaling factor to precipitation over the lake.

Second, through the hydrological sensitivity analysis, we demonstrate that net irrigation withdrawals and snowmelt and ice melt have a minimal impact on lake level fluctuations, indicating that this is primarily driven by rainfall and evaporation variability. However, this does not diminish the importance of glaciers. In fact, glaciers are significant at the scale of the headwater catchments, particularly for supplying water to large cities such as El Alto and La Paz (Buytaert et al., 2017; Soruco et al., 2015), maintaining wetlands like the bofedales (Herrera et al., 2015), and supporting irrigation (Buytaert et al., 2017). In most gauged catchments, incorporating irrigation resulted in only slight improvements in modeling performance. Nonetheless, this approach made it possible to estimate the net consumption due to irrigation at the scale of the catchments that contribute to the lake. Although this consumption is currently low, it is expected to increase significantly due the climatic and anthropogenic changes in the study area. It should also be noted that this process was based not only on soil water deficit, but also on local knowledge of farmers' water allocation practices.

Third, we disentangle the role of the change in heat storage in estimating E_lake. Annual evaporation (1616 mm yr⁻¹) is comparable to the evaporation (∼ 1600) of other low-latitude lakes (Wang et al., 2018) and aligns with previous studies of Lake Titicaca (∼ 1600 mm yr⁻¹) (Delclaux et al., 2007; Pillco Zolá et al., 2019). Despite the low air temperature over Lake Titicaca due to its high altitude, the evaporation rate is quite high. This is largely due to net radiation, although humidity, wind speed, and changes in heat storage also play a significant role in the seasonal variation. Regarding heat storage changes, the lake reaches maximum heat gain in October and maximum heat loss in May (see Fig. 13a). Neglecting changes in heat storage leads to overestimation of evaporation during the lake's heating period and underestimation of evaporation during the cooling period (see Fig. 13b). A similar finding was reported by Bai and Wang (2023) for Lake Taihu in China. Although several studies have investigated evaporation from Lake Titicaca using various methods (Carmouze, 1992; Delclaux et al., 2007; Pillco Zolá et al., 2019), our estimates are innovative because they are based on long-term data, including recent periods. Additionally, the accuracy of our estimates was underpinned by a water balance with a small error term, which enhances the reliability of our findings.

Figure 13 (a) Seasonal variations in heat storage and (b) the role of heat storage in the seasonal variation of evaporation from the lake. The figures are based on long-term (1982–2016) average values.

The periods of rising and falling water levels are closely linked to direct precipitation over the lake and upstream inflows (see Fig. 11a), which is mostly influenced by interannual precipitation variability. Understanding the effects of climate oscillations on precipitation variability is therefore crucial for understanding water level changes. Some authors (e.g., Garreaud and Aceituno, 2001; Jonaitis et al., 2021) noted that the interannual variability of precipitation in the region is mainly driven by the El Niño–Southern Oscillation (ENSO). During its warm phase, conditions are typically dry, while during the cold phase conditions are usually wet, although this relationship is not always consistent (Garreaud et al., 2003). For instance, Jonaitis et al. (2021) observed negative precipitation anomalies during La Niña phases and positive precipitation anomalies during El Niño phases in the Lake Titicaca region, though these anomalies were not statistically significant. Segura et al. (2016) argued that El Niño plays an important role in interannual precipitation variability and that decadal and interdecadal variations are influenced by sea surface temperature (SST) anomalies in the central–western Pacific. Therefore, variations in water levels cannot be attributed to ENSO alone. Sulca et al. (2024) found that interannual variations in water levels are related to SST anomalies in the southern South Atlantic and that interdecadal and multidecadal variability can be explained by Pacific and Atlantic SST anomalies. Additionally, they noted that multidecadal variations are linked to North Atlantic SST anomalies and southern South Atlantic SST anomalies.

5.2 Limitations of the modeling framework

Model forcing and evaluation data are the main sources of uncertainty. Daily precipitation and air temperature from GMET rely heavily on the spatial representativeness of ground-based measurements, which are very sparse at elevations above 4000 m a.s.l. As a result, the estimated snow and ice processes above that elevation may be affected by significant inaccuracies. In part, this could explain the poorly simulated mass balances at some glaciers compared to the geodetic observations (which are also uncertain). As shown by Ruelland (2020), the lack of high-altitude stations can have serious implications for correct estimations of precipitation and temperature lapse rates and thus the difficulty in realistically representing the accumulation and ablation processes. The precipitation estimated by GMET on the lake can also be questioned because it is only based on a few stations located on certain islands, which does not guarantee that they will correctly represent local convective phenomena (Gu et al., 2016; Nicholson, 2023), due to the very large surface area of the lake. Consequently, the underestimation of lake water levels in the early 2000s may be associated with underestimation of precipitation.

The simulation of snow and ice processes also needs discussing. The degree-day approach used in the present study might not be perfectly suited to simulating the melting of snow and ice at daily time steps in tropical glaciers, where melting energy is not always correlated with air temperature (Sicart et al., 2008; Rabatel et al., 2013). It may be more appropriate to use models based on energy balance, although they require large quantities of data. For instance, Frans et al. (2015) applied the DHSVM driven by reanalysis data fitted with local measurements in the Zongo catchment. For streamflow prediction over the period 1992–2010, while the performance they obtained at the monthly time step was satisfactory, this was not the case when the evaluation was conducted at a daily time step. The authors argue that the poor performance was due to inaccuracies in the forcing data. This suggests that using a more sophisticated model does not necessarily lead to more realistic streamflow simulations and that uncertainty in the input data may be more important than the structure of the model used. Another limitation is the direct transfer of snow and ice model parameters obtained in the Zongo catchment to the whole Titicaca hydrosystem, although temperature-index-based model parameter transfer has been shown to work relatively well, i.e., with only a small reduction in model performance in the southern Swiss Alps (Carenzo et al., 2009) and western Canada (Shea et al., 2009). The unsatisfactory performance of mass balance simulations on some glaciers may rather be due to inaccurate precipitation and air temperature data. On the other hand, the estimated glacier mass balances did not consider changes in glacier area, thickness, and volume over time. The variation in ice stock over the period 1982–2016 is negative (−5 mm yr⁻¹), which is consistent with the effects of global warming and in agreement with in situ observations (Rabatel et al., 2013). As the surface area of the glaciers in our model was obtained in 2000 and glaciers are shrinking worldwide, melting before the year 2000 may be underestimated, and after the year 2000 it may be overestimated. However, the biases are limited to some extent by the choice of an intermediate glacier area for the modeled period. If it is intended to simulate future changes in glaciers, it may be beneficial to include morphometric glacier changes in the model, drawing inspiration from simple approaches in the literature (e.g., Seibert et al., 2018). To initialize the model, global glacier thickness datasets could be used (e.g., Farinotti et al., 2019; Millan et al., 2022)

Estimates of lake evaporation and reference evapotranspiration also have certain limitations that are worth mentioning. Reanalysis data were used for some forcing data (humidity, solar radiation, and wind speed). Wind speed data were adjusted using the bias found at one meteorological station, which is not necessarily representative of the entire spatial domain. However, we believe that the impact is not very significant because the aerodynamic component only accounts for about 20 % of the reference evapotranspiration and evaporation. This is consistent with previous research in the study area, which indicated that, rather than the aerodynamic component, it is the radiative component that contributes significantly to the reference evapotranspiration (Garcia et al., 2004) and evaporation from the lake (Delclaux et al., 2007). Concerning evaporation, the estimated change in heat storage is based on certain assumptions that could be questionable. Instead of water temperature at different depths, we used LSWT, as suggested by Pillco Zolá et al. (2019), and the magnitude and seasonal variation in our evaporation estimates are in agreement with those of Carmouze (1992) based on temperature measurements taken at different depths.

Concerning irrigation, the main limitations are linked to the stationarity of the irrigable area and crop types. Some authors (e.g., Geerts et al., 2006) reported an increase in quinoa production in the Altiplano. However, quinoa is usually a rainfed crop and consequently requires limited irrigation. Information was also missing concerning changes in the crops cultivated over time. K_c derived from satellite images (Pôças et al., 2020) could be an interesting way of filling this gap. However, small irrigation systems predominate in our study area, which could have implications for the accuracy of the estimates. In addition, the locations of the inventoried irrigation systems are referential, since there are no maps showing delimited areas and it would consequently first be necessary to delimit the irrigable area and then use spectral indices to derive the crop coefficients. Despite these limitations, it should be noted that the reference evapotranspiration was driven by time-varying meteorological data. The water requirement is thus not entirely stationary.

The contributions of groundwater in the catchments were only simulated with the deep store of SMM. This is an important simplification because it neglects the dynamic interaction between rivers and aquifers. The use of pumped groundwater for irrigation was also not taken into consideration because irrigation inventories (Ministerio de Medio Ambiente y Agua, 2013) indicated that the proportion of groundwater in the supply was very small. Therefore, we would expect the impacts on streamflow prediction and water supply results to be negligible. Groundwater–lake interaction was also neglected. Figure 11b shows that the error term exhibits a seasonal variation, being positive during the rainy season and predominantly negative during the dry season. Linking the error term to net groundwater flow suggests that groundwater–lake interactions are seasonally variable. Lake water levels fluctuate by an average of 67 cm over the hydrological year, reaching a maximum in April and a minimum in December. In Fig. 11b, the error term has the highest positive values between December (26 mm) and January (54 mm), indicating a gain in net groundwater flow to the lake. When the lake reaches high water levels, losses to groundwater tend to dominate. This dynamic suggests that there could be a reversal of the hydraulic gradient throughout the year depending on the water level of the lake and the groundwater. However, it is important to note that the error term reflects not only the net groundwater flow, but also the uncertainty in estimating the other components of the water balance.

5.3 Prospects

Despite the aforementioned limitations, our modeling framework and its associated results on the water balance components will be useful for supporting decision-making in water resource management in Lake Titicaca because they represent climatic, glacio-snow-hydrological, and water allocation components. Contrary to the perceptions of some stakeholders, who often attribute the lake's water level variations to water withdrawals for irrigation or glacier retreat, this study demonstrates that Lake Titicaca's variations are primarily driven by rainfall and evaporation variability. In the next stage, water management scenarios could be evaluated in the context of climate change. Some of the management scenarios that could be explored include increasing the irrigable area, the efficiency of irrigation systems, and lake releases. Designing these scenarios would require collaborative work with stakeholders, in particular with the Lake Titicaca Authority. Exploring future scenarios will be crucial for identifying and planning intervention actions to ensure the sustainability of Lake Titicaca. Such experiments could also serve as a replicable model for other poorly gauged large lakes around the world. The conceptual models within the modeling framework are easy to apply, require minimal data, and are computationally inexpensive. Several of these models are part of WEAP, which is openly accessible to developing countries (for academic purposes and public institutions). Additionally, we provide in the current article detailed equations for models that are not included in this platform.

Appendix A: Performance of the ERA5-Land dataset

The ERA5-Land (Muñoz-Sabater et al., 2021) dataset is available from 1950 to the present at a spatial resolution of 0.1° and an hourly time resolution. Data were aggregated at a daily time step. The quality of the ERA5-Land dataset was evaluated using data recorded at a station located at El Alto airport as it was the only one with reliable humidity and wind speed data and relatively long time series. For relative humidity, the performance was acceptable, as ERA5-Land was able to adequately reproduce both the annual magnitudes and the seasonal cycle. For wind speed, ERA5-Land satisfactorily represented the seasonal cycle albeit with a significant systematic bias (see Fig. A1a). The bias evaluated at El Alto airport was used to correct the wind speed (see Fig. A1b) in the spatial domain of the study.

Figure A1 Wind speed performance of ERA5-Land evaluated at El Alto gauge in terms of (a) seasonal cycle and (b) interannual variation for the period 2003–2020.

Appendix B: Simulation of LSWT

A conceptual lumped model called Air2Water was thus used to estimate LSWT from air temperature (Piccolroaz et al., 2013; Toffolon et al., 2014). The Air2wateR 2.0.0 version available in R was used. ARC-Lake V3 data (MacCallum and Merchant, 2012) available for the period 1995–2012 were used to calibrate and evaluate the Air2wateR model. The spatial resolution of ARC-Lake V3 is 0.05°, and the temporal resolution is daily. The model was calibrated for the period 1995–2003 and evaluated using the period 2004–2012 (see Fig. B1). The performance of Air2Water in the independent control period was acceptable and captured seasonal and interannual variations very well. The LSWT of ARC-Lake V3 had previously been compared with data measured at a buoy located in Lake Titicaca that recorded LSWT (Lazzaro et al., 2021) between 2019 and 2023. The comparison was made between the cell intersecting the measurement point and compared the range of fluctuations and seasonality, since there is no temporal link between the two sources of data. Measured LSWT and LSWT obtained from ARC-Lake V3 fluctuated in a similar range and showed the same seasonal cycle (see Fig. B2).

Figure B1 LSWT simulated with the Air2wateR model (Piccolroaz et al., 2013; Toffolon et al., 2014) for Lake Titicaca and its calibration and evaluation performances.

Figure B2 LSWT obtained from the buoy site (16.25° S, 68.68° W) in Arc-Lake V3 (MacCallum and Merchant, 2012) and OPLT (Lazzaro et al., 2021).

Appendix C: Parameters of the snow and ice model

On-site observations are only available for the Zongo catchment, which is the glacier-dominated catchment adjacent to our study area (see Fig. 1). Streamflows are measured at the Tubo gauge for a catchment of 3.5 km², of which 52 % is glacierized area. The annual mass balance estimated by the glaciological method (based on direct sampling of changes in the mass of the glacier) was also available for the hydrological period 1991–2016 (see Fig. C1b). The approach used was to obtain the snow and ice model parameters (see the parameters in Fig. 4 and Table C1) in the Zongo catchment and then transfer them directly to BasicModel+IRR+SNOW/ICE. Firstly, the catchment model parameters were calibrated against the streamflow of the Tubo gauge for the hydrological period 2000–2010. The previous 2 hydrological years were used as a warmup period. The irrigation module was disabled because the non-glacierized area was completely dominated by rocks and no crops are cultivated at the high elevations of the Zongo catchment. The parameters with the best performance against the observed streamflow were obtained from 10 000 simulations generated with random hypercube sampling of the Monte Carlo approach. The range of parameters tested is shown in Table C1. NSE_sqrt was used as the objective function to select the best-performing parameter sets for the calibration period. Secondly, the model was evaluated against the mass balance observed using the glaciological method (see Sect. 2.5.1) for the period 1992–2016 and using the RMSE criterion. The evaluation was performed with other internal variables of the model, because the goal was to obtain the set of parameters associated with the snow and ice stores. Good internal consistency of the glacier model was thus desirable. In addition, the 25-year time window made it possible to evaluate the ability of the model to simulate the interannual and decadal variation in the mass balance of the Zongo glacier. Figure C1 shows the model's calibration (see Fig. C1a) and evaluation (see Fig. C1b) performances. Figure C1a shows that the model reproduced streamflow reasonably well. Mass balance was also simulated reasonably well. The RMSE obtained (349 mm) is in the range of the error of the mass balance based on the glaciological method (400 mm) (Sicart et al., 2007).

Figure C1 Performance of the catchment model in the Zongo catchment versus (a) the observed streamflow at the Tubo gauge and (b) the annual mass balance estimated using the glaciological method. The mass balance estimated using the glaciological method has an error of 400 mm yr⁻¹ (Sicart et al., 2007).

Table C1 Parameters of the model implemented in the Zongo catchment and their associated fixed values or ranges tested in order to obtain the parameters of the snow and ice model.

Appendix D: Simulated and observed water levels

Figure D1 Performance of the modeling chain compared with the water levels of Lake Titicaca measured in Puno. Three simulations of water levels are presented because the models produced three different estimates of Q_in.