Recent ground thermo-hydrological changes in a southern Tibetan endorheic catchment and implications for lake level changes

Climate change modifies the water and energy fluxes between the atmosphere and the surface in mountainous regions such as the Qinghai–Tibet Plateau (QTP), which has shown substantial hydrological changes over the last decades, including rapid lake level variations. The ground across the QTP hosts either permafrost or is seasonally frozen, and, in this environment, the ground thermal regime influences liquid water availability, evaporation and runoff. Consequently, climate-induced changes in the ground thermal regime may contribute to variations in lake levels, but the validity of this hypothesis has yet to be established.

This study focuses on the cryo-hydrology of the catchment of Lake Paiku (southern Tibet) for the 1980–2019 period. We process ERA5 data with downscaling and clustering tools (TopoSCALE, TopoSUB) to account for the spatial variability of the climate in our forcing data (Fiddes and Gruber, 2012, 2014). We use a distributed setup of the CryoGrid community model (version 1.0) to quantify thermo-hydrological changes in the ground during this period. Forcing data and simulation outputs are validated with data from a weather station, surface temperature loggers and observations of lake level variations. Our lake budget reconstruction shows that the main water input to the lake is direct precipitation (310 mm yr⁻¹), followed by glacier runoff (280 mm yr⁻¹) and land runoff (180 mm yr⁻¹). However, altogether these components do not offset evaporation (860 mm yr⁻¹).

Our results show that both seasonal frozen ground and permafrost have warmed (0.17 ^∘C per decade 2 m deep), increasing the availability of liquid water in the ground and the duration of seasonal thaw. Correlations with annual values suggest that both phenomena promote evaporation and runoff. Yet, ground warming drives a strong increase in subsurface runoff so that the runoff$/$(evaporation + runoff) ratio increases over time. This increase likely contributed to stabilizing the lake level decrease after 2010.

Summer evaporation is an important energy sink, and we find active-layer deepening only where evaporation is limited. The presence of permafrost is found to promote evaporation at the expense of runoff, consistently with recent studies suggesting that a shallow active layer maintains higher water contents close to the surface. However, this relationship seems to be climate dependent, and we show that a colder and wetter climate produces the opposite effect. Although the present study was performed at the catchment scale, we suggest that this ambivalent influence of permafrost may help to understand the contrasting lake level variations observed between the south and north of the QTP, opening new perspectives for future investigations.

1 Introduction

Climate change is amplified in mountainous environments, with major consequences for ecosystems, landscapes, hydrology, human communities and infrastructure (IPCC, 2019). Station observations show that global warming is elevation dependent, with the strongest warming rates being observed at high elevations (Pepin et al., 2015; Wang et al., 2014). Over the Qinghai–Tibet Plateau (QTP), a significant increase in surface air temperatures has been recorded since the 1980s, particularly in the north of the plateau (Zhang et al., 2022). This has been accompanied by a decrease in wind speed, humidification of the air and a general increase in precipitation, although with a strong spatial variability (Bibi et al., 2018). Altogether, these changes have affected the surface energy balance of the plateau through a shift of the Bowen ratio towards more latent heat fluxes, limiting the sensible surface warming (K. Yang et al., 2014).

These changes in water and energy fluxes between the atmosphere and the surface have the potential to alter the hydrological cycle of the QTP, which is the headwater region for major Asian rivers. As such, increasing trends of evaporation over land have been measured (3.8 mm per decade since the 1960s) with strong spatial variability in both the absolute values and increase rates (G. Wang et al., 2020). Changes in the seasonality of river discharge (Cao et al., 2006) and groundwater discharge (Niu et al., 2011) were reported for the same period. Overall, glacier shrinkage has also been observed since the 1960s, with a persistent increase in glacier mass loss rates (Bhattacharya et al., 2021; Hugonnet et al., 2021).

The QTP also features more than 1000 lakes larger than 1 km² (Zhang et al., 2017), most of them located in endorheic catchments. Lake volume changes are therefore attributable to climatic and hydrological changes occurring within the lake catchment, such as glacier melt, ground ice melt, precipitation, evaporation or runoff patterns. The majority of these lakes have experienced a pronounced increase in water levels since the 1990s (Lei et al., 2013, 2014), a trend that was suggested to be mainly driven by changes in precipitation and evaporation patterns (Yao et al., 2018) rather than by an increase in glacier mass loss and runoff (Brun et al., 2020; Zhang et al., 2021a). Nevertheless, lake level variations are not uniform across the QTP and exhibit important spatial variability. Whereas the northern and central QTP have recorded lake expansion, the southern parts of the plateau have experienced lake shrinkage (Qiao et al., 2019; G. Zhang et al., 2021a, 2020). Shrinking lakes have received less attention in the literature than rising lakes because they are fewer. For this reason, the drivers of this shrinkage are still unclear. Qiao et al. (2019) reported that recent lake shrinkage over the QTP could be driven by local precipitation decreases and/or evaporation increases (in relation to air temperature increases). G. Zhang et al. (2020) suggest that the divergent trends in lake level variations across the QTP could be linked to the contrasting evolution of moisture transport between the north and south of the plateau. On longer timescales, lake shrinkage over the QTP during the Holocene seems to be related to variations in the intensity of the Asian monsoon (Chen et al., 2013). Overall, such a complex pattern of rising and shrinking lakes challenges our understanding of the hydrological changes occurring in these High Asian watersheds.

In this regard, new insights into hydroclimatic changes over the QTP can emerge from the investigation of the coupled energy and water fluxes between the ground surface or subsurface and the atmospheric boundary layer. These fluxes are driven by the climate and have a major impact on cold-region hydrology (Pomeroy et al., 2007; Gao et al., 2021; Bring et al., 2016). Indeed, hydrological variables (precipitation, evaporation, runoff) affect the soil water content, which changes its thermal properties and the distribution between latent and sensible fluxes and thus substantially influences the ground thermal regime (Bring et al., 2016; Koren et al., 1999; Martin et al., 2019). In turn, the ground thermal regime modifies the relative proportion of frozen and liquid subsurface water, influencing infiltration possibilities and the amount of water available for evaporation and surface or subsurface runoff (Yi et al., 2006; Carey and Woo, 2001).

So far, climate-induced thermo-hydrological changes over the QTP have received limited attention. Large-scale modeling studies reported changes in the seasonal ground freezing cycles, characterized by a reduction in the frost depth and duration of the frozen period since the 1960s (Qin et al., 2018; C. Wang et al., 2020) and notable ground warming trends in summer and winter (Qin et al., 2021). Similar ground warming trends were reported in the regional modeling study by Qin et al. (2017), along with an increasing trend in evaporation and a decrease in the runoff coefficient over time. Plateau-scale surface energy balance modeling from G. Wang et al. (2020) reported that increasing trends in evapotranspiration could be mainly explained by variations in air temperature and net radiation at the surface.

Complementarily to seasonally frozen ground, permafrost is also a distinctive feature of climate–surface interactions in cold regions. Large-scale permafrost modeling suggests that it covers a significant part of the QTP, mainly as continuous permafrost in the north of the plateau and as discontinuous or sporadic permafrost in the south (Obu et al., 2019). Permafrost on the QTP usually has a low ice content due to limited precipitation and strong evaporation (Wu et al., 2005; Yang et al., 2010). Borehole temperature measurements show that it is a relatively warm type of permafrost (Biskaborn et al., 2019; Wu and Zhang, 2008), and its exposure to high solar radiations makes it sensitive to changes in surface conditions and climate change (Yang et al., 2010). Since the 1960s, climate change has driven permafrost warming across the plateau (Ran et al., 2018; Shaoling et al., 2000). Ran et al. (2018) report that most of the plateau exhibits a warming trend of the ground, comprised of between 0.26 and 0.74 ^∘C per decade, and half of the plateau warms at a rate higher than 0.5 ^∘C per decade. This warming is accompanied by upward migration (of around 100 m between the 1960s and 2000s) and shrinkage of permafrost-covered areas (24 % of the permafrost extent lost between the 1960s and the 2000s; Ran et al., 2018).

Permafrost grounds are characterized by a strong interplay between the ground thermal regime and the land hydrology. Seasonal thawing and freezing of the active layer are driven by the surface energy balance, which, in return, influences surface and subsurface runoff (Kurylyk et al., 2014; Walvoord and Kurylyk, 2016; Sjöberg et al., 2021) and evaporation (Gao et al., 2021). In this regard, both large-scale and regional modeling indicate that thawing permafrost enhances evapotranspiration (Qin et al., 2017; G. Wang et al., 2020). Qin et al. (2017) also report that the increase in evaporation is logically concomitant with a decrease in the runoff coefficient. Additionally, permafrost stores water as ground ice, and its thawing can trigger the release of liquid water in the watershed, contributing up to 15 % of the annual river streamflow (Cheng and Jin, 2013; Yang et al., 2019).

These hydrological changes are tied to various interdependent climate-driven physical processes happening at the ground surface and subsurface (e.g., surface energy balance, infiltration, water phase change, heat conduction). Because these processes exhibit a strong spatial variability in high-mountain environments, it is challenging to represent them accurately together on large spatial scales. Therefore, a deeper understanding of the impact of ground thermo-hydrological changes on the High Asia water cycle can be gained through small-scale physical modeling of these processes. Yet, for now, physics-based approaches at the catchment scale aiming to connect the ground thermo-hydrological regime and the observed hydrological changes on the QTP (such as lake level changes) remain scarce. They are, however, a powerful approach to tackle the following question: how much might climate-driven ground thermal changes affect the water cycle in high-mountain headwater regions? In this study, we use physical land surface modeling to quantify the ground thermo-hydrological changes in an endorheic Tibetan catchment over the last 40 years in response to climate change. We show the interplay in the water and energy fluxes occurring between the atmosphere, the surface and the subsurface and discuss their impact on the hydrology of the catchment and their implications regarding lake level variations.

2 Study area: the Paiku catchment

The Paiku catchment is located in southwestern Tibet, China, close to the border with Nepal (28.8^∘ N, 85.6^∘ E; Fig. 1). Its southern edge lies 7 km from the Shishapangma peak (8027 m a.s.l.). The catchment is endorheic and spans over 78 km from north to south and 66 km from east to west and covers 2400 km². The median elevation of the catchment is 4872 m a.s.l., ranging from 7272 m a.s.l. to its lowest point, with Lake Paiku at 4580 m a.s.l.. Geologically, the catchment is mainly located in the Tethys Himalaya, and thus, an important part of the formations underlying the catchment are metamorphized sedimentary series (Fig. B1 in the Appendix). The southern part of the catchment crosses the Southern Tibetan Detachment, and thus, the southern ridges of the massif belong to the High Himalayan metamorphic formations in the west and to the High Himalayan leucogranites of the Shishapangma massif in the east. The north and northeast ridges are formed by granite intrusions surrounded by metamorphic domes. The inner part of the catchment presents Plio-Quaternary formations such as alluvial fans close to the ridges and inclined alluvial plains in its inner parts (Aoya et al., 2005; Searle et al., 1997; Wünnemann et al., 2015).

Figure 1 The Paiku catchment. (A) Topographic and hydrologic map of the catchment with the glaciers in white, the ephemeral rivers in dark blue and the lake in light blue (elevation: SRTM data courtesy of the US Geological Survey). AWS: automatic weather station. L13 and L15 are surface temperature loggers (Sect. 3.1). (B) Localization of the catchment over the QTP. (C) Monthly temperature and precipitation recorded at the AWS between October 2019 and September 2021.

Automatic weather station (AWS) observations (5033 m a.s.l., October 2019–September 2021; Fig. 1) show that the climate in the catchment is characterized by a relatively small temperature amplitude during the year (around 20 ^∘C, JJA being the warmest months and DJF being the coldest) and significant daily amplitude (up to 10 ^∘C during the warm season). The mean annual temperature is −1.5 ^∘C at the AWS, where night freezing can occur until the beginning of June and resume at the beginning of October. The catchment is dry (200–300 mm yr⁻¹), and precipitation mostly falls as rain during the monsoon (JJAS).

Around 5 % of the catchment is covered by glaciers (RGI Consortium, 2017), which are concentrated in its southwestern part. They feed several proglacial lakes that can reach up to 6 km in length. Geodetic glacier mass budgets show that, similarly to other glaciers in the region, glaciers of the Paiku catchment have undergone sustained mass loss at least since the 1970s, with an average mass balance of −0.3 m water equivalent per year ($\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$) until the beginning of the 2000s and around −0.4 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ thereafter (Bhattacharya et al., 2021). There are more than 10 rivers that drain the catchment towards the lake, and most of them only exhibit a seasonal activity during the monsoon months. The three main ones are (Fig. 1) Daqu (glacier-fed, 450 km²), Bulaqu (glacier-fed, 325 km²) and Barixiongqu (non-glacier-fed, 703 km², Lei et al., 2018).

In the northwest of the catchment, Lake Paiku covers approx. 280 km² (11.5 % of the catchment surface area) and spans over 27 km from north to south. It has a mean water depth of 41 m, with a maximum water depth of 73 m (Lei et al., 2018). It receives water from direct precipitation and from land and glacier runoff, which can be routed at the surface via the river systems or at the subsurface via the alluvial formations. Because it is hydrologically closed, the lake mainly loses water through evaporation. Previous studies reported lake level fluctuations over different timescales. It reached 4665 m a.s.l. (85 m higher than the present level) prior to 25 ka and at the onset of the Holocene (11.9–9.5 ka); afterwards, the lake shrank gradually (Wünnemann et al., 2015). More recently, the lake level decreased by 3.7 m between 1972 and 2015, losing 4.2 % of its surface and 8.5 % of its volume. Measurements have been performed since the end of the 1970s and allow us to accurately know the evolution of the lake level until today (Lei et al., 2021, 2018); they are used in this study to validate our hydrological results (Sect. 3.2.1, Figs. 5d and 6b). At the seasonal scale, the lake level cycle has an amplitude of ∼0.4 m. It is marked by a strong increase during the monsoon period (JJAS), supported by direct precipitation, glacier melt and land runoff. From October and until the next monsoon period, evaporation dominates the lake mass budget, and the level decreases rapidly until January and at a slower rate afterwards (Lei et al., 2021).

3 Material and methods

3.1 Field measurements

An AWS was set up in October 2019 in the south of the catchment at an elevation of 5033 m a.s.l. (Fig. 1). It is equipped with various sensors which record air temperature, pressure, relative humidity, wind speed, incoming and outgoing long- and shortwave radiation, and precipitation every 15 min. The meteorological record extends to September 2021 and covers a period of nearly 2 years. We used it to evaluate and correct the distributed downscaled climatic forcing on which we rely in our modeling framework (Sect. 3.2.2).

Two temperature loggers recorded the surface temperature in the vicinity of the AWS location. Logger 15 (L15) is located at 5055 m a.s.l., 6 km west of the AWS. Logger 13 (L13) is located at 5356 m a.s.l., 12 km west of the AWS (Fig. 1). Both loggers were buried 10 to 15 cm below the surface to avoid direct solar radiation on the sensors and recorded surface temperature at a 20 min time step from October 2017 to October 2018. These surface temperature records were used to evaluate the simulations (Sect. 3.2.4).

3.2 Catchment thermo-hydrological modeling

3.2.1 Conceptual hydrological model for the catchment

To understand the level variations of Lake Paiku over the last 40 years (1980–2019 period), we develop an approach at the catchment scale. Because the catchment is hydrologically closed, the lake receives water input via direct precipitation, land surface and subsurface runoff, and glacier runoff. Conversely, it loses mass via evaporation. Because the quantification of water flows between the lake and potential aquifers surrounding it is difficult (Rosenberry et al., 2015), our approach assumes that these flows are negligible. The present study requires quantification of the different terms of the hydrological balance. Under these assumptions, the hydrological balance of the lake is given by the following equation:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}\mathrm{\Delta }{z}_{\text{Lake}}& ={\text{Precipitation}}_{\text{Lake}}+{\text{Runoff}}_{\text{LandSurf}}+{\text{Runoff}}_{\text{LandSub}}\\ & +{\text{Runoff}}_{\text{Glacier}}-{\text{Evaporation}}_{\text{Lake}}.\end{array}\end{array}$$

The production of forcing data for the catchment (including precipitation) is detailed in Sect. 3.2.2. The land hydrology processes are quantified using the CryoGrid community model (version 1.0) (Westermann et al., 2023), as described in Sect. 3.2.3. Distributed 1D simulations are used to quantify land evaporation and runoff. The routing of water in the catchment is not represented, and the runoff computed for a given simulation is directly accounted for as a water input for the lake. The evaporation from the lake is simulated using the CryoGrid3-Flake model (Langer et al., 2016), as described in Sect. 3.2.5. Glacier melt is not modeled but is estimated for the study period (1980–2019) from remote sensing observations. From these observations, glacier yield is calculated as described in Sect. 3.2.6. Our catchment-scale approach to represent the hydrological balance of the lake is summarized in Fig. 2. Based on this approach, we can evaluate the performance of our framework (Sect. 4.1.2) by comparing the simulated lake balance with the one derived from the detailed observations of lake level variations over the study period (Lei et al., 2018, 2021).

Figure 2 Conceptual hydrological framework for the study.

3.2.2 Forcing-data production and validation

In high-mountain environments, topography creates strong spatial variability in temperature and incoming radiation, which impact the surface energy balance (Klok and Oerlemans, 2002) and the ground thermo-hydrological regime (Magnin et al., 2017). Our approach requires forcing data that (i) capture this variability; (ii) include numerous variables such as air temperature, incoming long- and shortwave radiation, wind speed, specific humidity, rain, and snowfall; and (iii) cover the 40-year study period at a sub-daily time step. The TopoSCALE approach (Fiddes and Gruber, 2014) was developed for this purpose and allows us to downscale reanalysis products like ERA5 (Hersbach et al., 2020) at a high resolution (here $\sim \mathrm{100}\mathrm{m}\times \mathrm{100}\mathrm{m}$).

Additionally, because working at a 10⁻² km² spatial resolution over a 2400 km² catchment would require more than 200 000 forcing files and simulations, we rely on the TopoSUB method (Fiddes and Gruber, 2012) to reduce computational costs. This method uses a SRTM30 digital elevation model to explore redundancies in physiographic parameters of the study area such as elevation, aspect, slope and sky-view factor and to identify groups of high-resolution pixels (100 m×100 m) sharing similar values for these parameters. From there, all the high-resolution pixels belonging to such a group are only described as a single TopoSUB point, for which climatic variables can be downscaled to create one single dataset of climatic time series. The degree of similarity required by TopoSUB to identify groups of high-resolution pixels with redundant physiographic parameters can be adjusted by choosing the final number of TopoSUB points (and thus climate datasets) that should be used to cover the area corresponding to one ERA5 pixel. The Paiku catchment intersects eight ERA5 pixels at 30 km resolution, and we chose to use 50 TopoSUB points within each ERA5 pixel to cover the spatial variability created by the topography in small-scale climate. Ultimately, 368 TopoSUB points are used to cover the catchment. The average level of redundancy (i.e., the average number of high-resolution pixels represented by a single TopoSUB point) is 723 ± 745 (1σ, median: 506, min: 1, max: 4347). Figure C1 shows the distribution of the TopoSUB points and a reconstruction of the topography of the catchment based on this approach. The period covered by the forcing datasets starts on 1 January 1980 and ends on 31 August 2020 (40 years and 8 months).

In the TopoSCALE statistical downscaling approach, we do not rely on the AWS data, and thus the downscaled ERA5 data can be biased, as is often the case over Asia (Jiang et al., 2021, 2020; Jiao et al., 2021; Orsolini et al., 2019). Comparison against the available AWS observations (Fig. D1) indeed highlights notable differences in variables such as air temperature and precipitation. From these differences, we derived monthly bias correction factors that we applied systematically to all 368 of the climate forcing datasets. The catchment averages for precipitation and air temperature are shown in Fig. 3. In this figure and across the rest of the study, we use p values to evaluate the significance of linear trends in the temporal evolution of certain variables (temperature, precipitation, evaporation  etc.). This p value tests the null hypothesis which supposes that the value of the slope is equal to zero. The hypothesis is tested using Student's t test by comparing the distance between the estimated slope and 0 relative to the standard error of the slope. We did not report trends when this p value (probability of a null slope) was higher than 0.005.

Figure 3 Climate forcing data for the land and lake modeling. Annual catchment average air temperature (2 m above ground) and annual total precipitation for the study period. Note that the model is also forced by incoming short- and longwave radiation, humidity, wind speed, and air pressure. Details about the spatial and temporal resolutions of the distributed forcing data are presented in Sect. 3.2.2.

3.2.3 The CryoGrid community model (version 1.0)

To simulate the ground thermo-hydrological regime, we use the CryoGrid community model (Westermann et al., 2023). The CryoGrid community model (CG) is a land surface model designed for applications in cold regions where seasonal frozen ground or permafrost may occur. The model implements heat transfer in a 1D soil column, accounting for freeze–thaw processes of soil water using effective heat capacity (Nakano and Brown, 1972). To do so, soil-freezing curves are based on Dall'Amico et al. (2011), as detailed in Westermann et al. (2013). Vertical water movement in the soil column is based on Richards' equation (Richardson, 1922; Richards, 1931). The soil matric potential and hydraulic conductivity follow van Genuchten (1980) and Mualem (1976). Additionally, to represent the obstruction of connected porosity by ice formation, the hydraulic conductivity is reduced by a factor dependent on the local ice content, following Dall'Amico et al. (2011). The model features the snowpack module called CG Crocus described in Zweigel et al. (2021), which adapts the snow physics parameterizations from the CROCUS scheme (Vionnet et al., 2012) to the native snow module of CryoGrid3 (Westermann et al., 2016). At the surface, the model uses a surface energy balance module to calculate the ground surface temperature and water content. The turbulent fluxes of sensible and latent heat are calculated using a Monin–Obukhov approach (Monin and Obukhov, 1954). Evaporation is derived from the latent heat fluxes using the latent heat of evaporation and is adjusted to the available water in the soil. It occurs in the first grid cell only, but water can be drawn upwards due to matric potential differences. Because vegetation is very scarce in the catchment, we do not expect transpiration to have a strong imprint on evapotranspiration, and our calculations do not unravel evaporation from transpiration.

3.2.4 Model setup and validation

The setup of the CryoGrid community model for the land is presented in Fig. 4. To capture the high spatial variability of mountainous climate, our approach relies on the 368 climate forcing datasets to cover the catchment (see Sect. 3.2.2). This approach enables us to perform spatially distributed modeling. All 368 of the simulations are independent and use the same parameterization. In the absence of direct observations of the soil stratigraphy within the catchment, the soil column was designed to agree with field observations in the region (Yuan et al., 2020; Wang et al., 2009; Hu et al., 2020; Luo et al., 2020; S. Yang et al., 2014; Wang et al., 2008), to be consistent with similar modeling approaches across Tibet (Chen et al., 2018; Song et al., 2020) and to be consistent with input datasets (Shangguan et al., 2013, 2017). Thus, the soil stratigraphy is divided into three units: a top soil (0.3 m thick), a bottom soil (1.7 m thick) and a bedrock unit (extending beyond the depth of interest of the study). An overview of the parameters for each unit, their source and the way they are calculated is presented in Table A1 in the Appendix.

Regarding the processes implemented in the model (Sect. 3.2.3), infiltration according to Richards' equation only occurs in the top and bottom soil units. The bedrock unit has a static water content. Unraveling surface from subsurface flow is an ongoing challenge in catchment-scale hydrology (McDonnell, 2013), and this distinction is important in mountain terrains where these two flows can behave differently due to the complex topography (Seibert et al., 2003; Gao et al., 2014; Hu et al., 2020). For this study, we rely on a simple approach that is based on thresholds regarding the soil water content (porosity and field capacity). These kinds of approaches are thus based on soil properties and have often been used in hydrological modeling studies (Vörösmarty et al., 1989; Shaman et al., 2002; Kelleners et al., 2010; Kampf, 2011; Samuel et al., 2008). In detail, we compute surface and subsurface flow as follows.

On the one hand, surface runoff is computed relative to the saturation level of the soil column. When the entire soil column is saturated (WC = porosity), additional water input from precipitation or snow melt is directly counted as surface runoff. On the other hand, subsurface runoff is computed relative to the field capacity of the ground, which is an input parameter of the model. When the water content (WC) of a ground cell exceeds this field capacity (FC), the amount of water corresponding to WC minus FC is available to produce subsurface runoff. We use the lateral boundary condition LAT_WATER_RESERVOIR from the CryoGrid community model (Westermann et al., 2023) to account for this subsurface runoff. The speed at which this available water exits the soil column towards the lake is calculated with Darcy's law, using the hydrological conductivity of the ground and the mean slope of the catchment as hydraulic slope. Because the model couples thermal and hydrological fluxes, all of these changes in the soil water content can be driven by precipitation input and evaporation but also by water phase change in the ground, such as ice melt.

Because we do not have knowledge of the distributed thermal state with depth over the catchment at the beginning of the simulations, we assume that temperature profiles were in equilibrium with the climate of the 5 first years of modeling (1980–1984). To do so, we start our simulations with a 60-year spin-up of these first 5 years (12 repetitions), which is sufficient to establish a stable temperature profile over the first 9 to 80 m, depending on the simulations, extending beyond the hydrologically active part of the ground (the first 2 m).

Figure 4 Modeling framework for the land hydrology. ERA5 data are downscaled using the TopoSUB and TopoSCALE approaches (Fiddes and Gruber, 2014, 2012). The downscaled data are bias-corrected based on the AWS observations. Distributed 1D simulations are performed using the CryoGrid community model (Westermann et al., 2023). The vertical resolution is indicated with the tick marks on the depth axis.

To validate model simulations, the simulated ground surface temperatures (GSTs) are compared to the two temperature logger time series acquired in the vicinity of the AWS (Sect. 3.1). We used this comparison to calibrate the surface roughness used for the surface energy balance calculations in the model.

The following method is used to produce area-averaged evaporation and runoff (in mm, water equivalent) in a zone of interest. For a given TopoSUB point in this zone, the model produces hydrological values (in m³) using the area of a TopoSUB pixel on the catchment map. Then these values are multiplied by the number of pixels in the zone corresponding to this TopoSUB point in particular, and this is for all the relevant TopoSUB points covering the zone (e.g., evaporation in warm permafrost). Then the area of interest is calculated by counting the number of pixels in the zone of interest and multiplying this number by the area of a pixel. Then the total volume is divided by the total surface for the zone of interest to obtain the final value (in mm).

3.2.5 Lake modeling

The lake thermo-hydrological response to the climatic forcing data is simulated using the CryoGrid3-Flake model (Langer et al., 2016). The two models were coupled by Langer et al. (2016) to simulate the thermal regime of thermokarst lakes (including surficial water freezing and melting) and underlying ground. Here, we use the coupled models mainly to quantify evaporation at the lake surface. In the coupled model, the native surface energy balance module of CryoGrid3 (Westermann et al., 2016) was amended to account for processes tied to free water surface energy balance: (i) the dependence of the albedo of a water surface on solar angle (and thus time of the day) and wind speed (and wave formation), (ii) the dependence of the surface roughness length on wind speed (and wave formation), and (iii) the exponential decay of incoming radiation with depth in the water column. Similarly to the land simulations, the lake simulations were forced by the downscaled ERA5 data (with the TopoSUB and TopoSCALE methodology), with the corrections being derived from the AWS data (Sect. 3.2.2). The simulations were initiated with a 20-year spin-up of the 1980–1984 climate. The simulation results corresponding to the four ERA5 tiles covering the lake were then averaged using the respective spatial footprint of each tile on the lake.

3.2.6 Quantification of glacier mass change

Multiple studies quantified the volume change of the glaciers located within the Paiku catchment in the recent past (1970s to 2020). To our knowledge, there are no field-based measurements of glacier mass balance available in this catchment. As a consequence, we rely solely on geodetic mass balance studies (Brun et al., 2017; Maurer et al., 2019; King et al., 2019; Shean et al., 2020; Hugonnet et al., 2021). All these studies estimated glacier volume changes over periods of 20–30 years from satellite-derived DEMs. As a consequence, we can only estimate the average annual glacier mass balance and not the year-to-year variability. Glaciers occupy approximately 113 km² in the Paiku catchment. They have shrunk over the past 50 years at a rate of 0.44 % yr⁻¹ from an area of 132 km² in 1975 to 122 km² around 2000 and to their current extent (King et al., 2019; Bolch et al., 2019). The average mass balances for the period 1975–2000 and 2000–2020 are −3.9 ± 2.1×10¹⁰ and −5.4 ± 2.4×10¹⁰ kg yr⁻¹, respectively (−4.6 ± 2.5×10⁷ and −6.4 ± 2.8×10⁷ m³ with a 850 kg m⁻³ density). These mass balances correspond to specific mass balances of −0.31 ± 0.17 m of water equivalent per year ($\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$) and −0.47 ± 0.21 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$, respectively.

Regarding glacial runoff, it was estimated to be 320 ± 4 mm yr⁻¹ for the 2001–2010 period by Biskop et al. (2016) using a temperature index approach for ice melt. For the 2000–2018 period, X. Zhang et al. (2020) derived a runoff value of 52 ± 12 mm yr⁻¹ (1.24 ± 0.29×10⁸ m³ yr⁻¹, which we scaled to the basin area). The value we derive of 39 ± 13 mm yr⁻¹ is thus found to have good consistency with the latter one (Sect. 4.1.1).

4 Results

4.1 Model validation and hydrological budget of Lake Paiku

4.1.1 Model validation

Simulated daily ground surface temperatures are in good agreement with the observed ones, showing a bias of −0.2 and 0.6 ^∘C and an RMSE of 1.4 and 1.6 ^∘C for loggers 15 and 13, respectively (Fig. 5a and c). Most of this RMSE is explained by a mismatch between the model and observations in the tails of the temperature distribution, whereas intermediate temperatures exhibit the best agreement with observations.

Figure 5 Model validation. (a) and (c) Modeled mean daily ground surface temperatures compared to measured ground surface temperatures for logger 15 and logger 13 (location on Fig. 1). (b) Modeled annual lake evaporation (blue curve) and comparison with values calculated by Lei et al. (2021) in the light-blue zone. The gray curve shows the smoothed lake level relative to August 2019 based on observations from Lei et al. (2018) (gray points) and Lei et al. (2021) (oscillating gray line). (d) Comparison between the runoffs required to reproduce the observed lake variations (red curve, derived from lake level, lake area, forcing data and lake evaporation) and the sum of the glacier and land runoff we derive from remote sensing observations and modeling (Sect. 3.2). Error bars are associated with the glacier values and come from the geodetic results. Runoff values are expressed as heights scaled to the land surface of the Paiku catchment.

Annual lake evaporation mainly ranges between 800 and 900 mm yr⁻¹ (Fig. 5b), with a mean value of 870 ± 23 mm (1σ). Lake evaporation does not exhibit a linear trend of increasing or decreasing and is mostly dominated by year-to-year variability. Though slightly lower, our evaporation results are in good agreement with the values from Lei et al. (2021), which are derived from local and regional meteorological observations and lake budget calculations (Fig. 5b). We used the simulated evaporation together with the lake level data and lake area data from Lei et al. (2018) and Lei et al. (2021) and the precipitation forcing datasets (3.2.2) to derive the total runoff (land + glacier) required as an input to the lake budget to reproduce the lake variations. This required runoff corresponds to the red line of Fig. 5d. The required runoff volumes are scaled to the land area of the catchment to be comparable with the other variables. Figure 5d also presents the runoff values derived from the land cryo-hydrological modeling and from the glacier remote sensing investigations. Annual volumes are expressed (in mm) over the land part of the catchment (excluding the lake). As presented in Sect. 3.2.6, glacier mass balance values are considered to be constant for the 1980–2000 period and the 2000–2019 period and are, respectively, equal to −4.6 ± 2.5×10⁷ and −6.4 ± 2.8×10⁷ m³ yr⁻¹. The addition of annual precipitation to these values to quantify the total glacier runoff introduces year-to-year variability to the glacier runoff. At the catchment scale, the average glacier runoff over the 40 years is 39 ± 13 mm yr⁻¹.

Over the 40 years, the average annual land runoff value (surface + subsurface) we model is 24 ± 8 mm. Summed together, the land and glacier runoff are in partial agreement with the runoff that is required to close the lake water balance. Annual values are compatible within error bars for 28 out of the 40 years of simulations. The glacier and land runoff are slightly too small to close the lake water balance during the first 20 years and slightly too large for the last 20 years of simulation. Over the whole period, the sum of the glaciers + land runoff produces 95 % of the required runoff. Land runoff is further described in Sect. 4.3, and lake results are described in the following section.

4.1.2 Hydrological budget of Lake Paiku

Our observations, climate data, simulations, geodetic data and the lake level data from Lei et al. (2018, 2021) enable us to quantify the different terms of the lake hydrological budget. We present these results in meters for lake level change based on the average slope of the volume = f (level) relationship (Fig. 6). As the unique output term, evaporation dominates the lake budget with an average annual value of 0.86 m (34.6 m per 40 years; Fig. 6a). Direct precipitation in the lake is the dominant input, with an average annual value of 0.31 m (12.3 m per 40 years), followed by glacier runoff (0.28 m yr⁻¹, 11.3 m per 40 years) and land runoff (0.18 m yr⁻¹, 7.0 m per 40 years). When compared with lake volume observations over the 40 years of the simulation, the simulated lake budget is 1.04 m too negative.

Figure 6 Budget and level of Lake Paiku for the simulation period (1980–2019). (a) The different components of the hydrological budget of the lake according to our framework. Results are given in meters for lake change based on the average slope of the volume = f (level) relationship. (b) Lake level data. Points correspond to observations from Lei et al. (2018, 2021) that we smoothed (green curve, also based on observation points older than 1980). The simulated lake level appears in orange.

Based on our results, we also reconstructed lake level variations that we compare with the observed variations (Fig. 6b). Following our framework, our values are presented at an annual time step. They qualitatively reproduce the overall lake level decrease but tend to overestimate this decrease and show an increasing mismatch with the observations from 0 in 1980 to 2 m in 2005. This mismatch is later compensated for by an increasing lake level trend in our simulation from 2005 to 2019. At the end of the simulation period, the mismatch is 1.04 m, which is consistent with the budget values (Fig. 6a) and the fact that our approach provides 95 % of the required runoff to close the lake budget (Sect. 4.1.1). This pattern of a too-strong decrease followed by an increase is consistent with the comparison between simulated and required runoff presented in Fig. 5d.

4.2 Ground thermal results

Based on our temperature results, we define four categories of ground thermal regimes (Fig. 7a). Cold permafrost refers to the areas of the catchment for which the deepest thaw depth did not exceed 1 m over the 40 years of simulation. For cold permafrost, frozen conditions dominate the first meters of the ground most of the year, and surficial thawing during summer can be interrupted by ground freezing from the surface to the top of the permafrost at night. Warm permafrost refers to the areas of the catchment that present permafrost for the whole duration of the simulation and which are not part of the cold permafrost. These areas are characterized by a distinct seasonal pattern of frozen ground in winter and an active layer in summer. Disappearing permafrost refers to the areas of the catchment that present permafrost at the beginning of the simulation and not at the end. No permafrost refers to the areas without permafrost at the onset of the simulation. The geographical characteristics of each ground category are presented in Table 1, and their distributions throughout the catchment are shown in Fig. 7a. These different ground categories are subsequently used to compare their cryo-hydrological behaviors during the simulation (consistent color code).

Table 1 Cryological classification of the catchment based on the modeled ground temperatures.

At the catchment scale, the 2 m depth temperature (Fig. 7b) shows a pronounced warming trend of 0.17 ^∘C per decade ($p=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{6}}$). This trend is mainly supported by the no-permafrost areas, which underwent a slightly stronger warming trend of 0.2 ^∘C per decade ($p=\mathrm{7}\times {\mathrm{10}}^{-\mathrm{8}}$). Areas with disappearing permafrost, warm permafrost and cold permafrost exhibit smaller trends around 0.1 ^∘C per decade, with decreasing p values (respectively, 0.00001, 0.006 and 0.05 – i.e., non-significant for the last two).

From 1980 to 1989, permafrost covers 27 % of the catchment, and the mean active-layer thickness (ALT) is 1.36 ± 0.51 m (1σ, minimum: 0.11 m and maximum: 2.37 m; Fig. 7c). From 2010 to 2019, permafrost covers 22 % of the catchment. At the scale of the initial permafrost area, this change corresponds to a loss of 19 %. The mean ALT is 1.29 ± 0.49 m (1σ, minimum: 0.11 m and maximum: 2.55 m; Fig. 7d) for this period. Permafrost disappearance (gray zones in Fig. 7d) mainly happens for the low-lying permafrost in the south and at the center of the catchment. It occurs, for the most part, on the outer slopes of the permafrost regions and at the bottom of steep glacial valleys.

Figure 7 (a) Different cryological states of the ground throughout the catchment for the 1980–2019 period (see Table 1). (b) Mean annual ground temperature (MAGT) 2 m deep, averaged for the whole catchment and for the different cryological states of the ground. (c) Average active-layer depth over the 1980–1989 period. (d) Average active-layer depth over the 2010–2019 period. Only locations presenting permafrost at the end of the simulation are assigned a color in the maps in (c) and (d). Locations where permafrost has disappeared are shown in gray on (d).

We also present the average duration of seasonal thaw at a depth of 70 cm averaged over the catchment (Fig. 8a). Because at this depth some areas might present 2 (or more) consecutive years without thawing (highest locations) or without freezing (lowest locations), these areas were excluded from the averaging. In the end, the averaged results account for 89 % of the catchment land area (i.e., excluding glaciers and Lake Paiku). The results show an increasing trend in the duration of the seasonal thaw of +4.6 d per decade ($p=\mathrm{3}\times {\mathrm{10}}^{-\mathrm{4}}$; blue line in Fig. 8a). When looking at the average start and stop days of the seasonal thaw (Fig. 8a, gray lines) in the Julian calendar (day 150 is 30 May and day 300 is 27 October), we note that this increase is mainly caused by a later ending date of the thaw season (Stop date in Fig. 8a; +6.1 d per decade, $p=\mathrm{8}\times {\mathrm{10}}^{-\mathrm{10}}$) and not by an earlier starting date (non-significant trend).

Figure 8 (a) Duration of seasonal thaw 70 cm deep, averaged over the catchment. The asterisk indicates that the presented curves average 89 % of the surface of the catchment (Sect. 4.2). The gray curves and the light-blue area are associated with the right axis and indicate the average start and stop days of the seasonal thaw in the Julian calendar. Values higher than 365 indicate that freezing conditions came back after 31 December. (b) Active-layer thickness (ALT) evolution for warm permafrost. The solid line shows the ALT for simulations experiencing an annual evaporation lower than 150 mm when averaged over the 40 years. The dashed line shows the ALT for simulations with annual evaporation higher than 150 mm. (c) Temporal trends for seasonally frozen ground where there is no permafrost. The asterisk indicates that simulations were excluded if one of the simulated years did not present freezing conditions 70 cm deep (persistence of thawed conditions from one year to another). The presented curves thus average 88 % of the total permafrost-free areas of the catchment. (d) Altitudinal distribution of permafrost in 1980 and 2019. This distribution includes both cold and warm permafrost.

Within warm permafrost, we distinguished the ALT for locations experiencing an average evaporation lower or higher than 150 mm yr⁻¹ during the simulations (Fig. 8b). Whereas locations with average evaporation below 150 mm yr⁻¹ record an active-layer deepening trend of 2.7 cm per decade (p=0.001), this is not the case for locations with an average evaporation higher than 150 mm yr⁻¹ (non-significative trend). This threshold value of 150 mm yr⁻¹ is based on further investigations of the relationships between evaporation and ALT provided in Sect. 5.3.1.

In the permafrost-free areas of the catchment, seasonal frozen ground (Fig. 8c) reaches a depth of 1.43 ± 0.15 m on average and shows a decreasing trend of −6.8 cm per decade ($p=\mathrm{6}\times {\mathrm{10}}^{-\mathrm{4}}$). At 70 cm depth, the average duration of seasonally frozen ground is 136 ± 12 d, with a decreasing trend of −5.3 d per decade ($p=\mathrm{4}\times {\mathrm{10}}^{-\mathrm{4}}$). These values average 88 % of the no-permafrost areas since locations showing persistent thawed conditions at this depth from one year to another were excluded (i.e., minimal seasonal freezing depth lower than 70 cm over the 40 years).

When comparing permafrost spatial distribution between 1980 and 2019 (Fig. 8d), our results show that permafrost distribution above 5750 m a.s.l. was not modified during the simulation. Permafrost disappearance mainly occurred between 5000 and 5750 m a.s.l., with the largest loss reaching 2.5 % of the catchment area between 5250 and 5500 m a.s.l.

4.3 Hydrological results for the land

The mean annual evaporation (land area only) over the simulation time is 180 ± 19 mm (1σ; Fig. 9a). Evaporation shows an increasing trend over the 40 years of +1.01 mm per decade ($p=\mathrm{3}\times {\mathrm{10}}^{-\mathrm{7}}$). Average total runoff over the 40 years is 24 ± 8 mm yr⁻¹ (Fig. 9b) and exhibits an increasing trend of +4.8 mm per decade ($p=\mathrm{8}\times {\mathrm{10}}^{-\mathrm{7}}$). Similarly, surface runoff (13 ± 3 mm yr⁻¹) and subsurface runoff (11 ± 6 mm yr⁻¹) show increasing trends of +1.3 and +3.5 mm per decade ($p=\mathrm{6}\times {\mathrm{10}}^{-\mathrm{5}}$ and $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{7}}$), respectively (Fig. 9b). The surface runoff presented in Fig. 9b includes the snow melt that did not infiltrate the ground. These linear trends we report are high compared to the absolute values of the variables, and their extrapolation backward in time would lead to null values in the recent past, which is unrealistic. This suggests a non-linear evolution of these variables over the 20th century.

We also present the catchment average of the runoff$/$(runoff + evaporation) ratio (Fig. 9c), which is equivalent to runoff$/$(rain + snow − snow sublimation) given the negligible contribution of soil storage variations. Hence, it is the proportion of the water input to the ground surface that is converted into runoff. This proportion is 11 % ± 2 % over the simulation time and shows an increasing trend of +1.23 % per decade ($p=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{7}}$). Figure 9c also shows the average theoretical ratio to maintain a steady lake level (of 17.6 %). This ratio was obtained under the following hypotheses:

• The climate forcing data are the same, and hence, the lake evaporation is the same.

• The glacier contribution is (i) considered to be the same for the historical simulation and this scenario and (ii) is taken as the difference between the total land surface runoff and the red curve of required runoff in Fig. 5; therefore, it is independent of remotely sensed estimates.

• Under these conditions, the runoff increase needed to maintain the lake level is only supplied by land runoff (surface and subsurface) by shifting the runoff$/$(runoff + evaporation) ratio.

The ratio from the historical simulation starts significantly below the theoretical steady lake ratio (10.2 %<17.6 %; Fig. 9c) and increases progressively to 16.0 % in 2019.

Finally, Fig. 9d shows the annual proportion of the liquid / total water averaged for the whole catchment. The value was computed based on the daily water content (liquid and frozen) of the first 2 m of the soil column (the hydrologically active part of the column, Sect. 3.2.4) from which annual averages were derived and used to compute a catchment-scale average. The graph shows that the proportion of liquid water in the total water content increases at around +1.41 % per decade ($p=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}$), indicating that water spends more and more time in the ground in a liquid form, thus being increasingly available for hydrological processes such as evaporation or runoff.

Figure 9 Hydrological results. (a) Annual evaporation averaged over the whole catchment (land area only). (b) Annual runoff averaged over the whole catchment (land area only). The blue curve sums the surface and subsurface runoff. (c) Ratio between runoff and evaporation + runoff averaged over the whole catchment (land area only). The brown line indicates the theoretical average ratio needed to maintain a steady lake level when considering an identical glacier contribution to runoff (details in Sect. 4.3). (d) Annual mean of the (liquid water)$/$(total water) ratio over the first 2 m of ground, averaged over the whole catchment (land area only).

4.4 Sensitivity of evaporation and runoff

We conducted a simple sensitivity test on the climatic conditions (i.e., not a full-scale sensitivity test). We ran the same 40 years of simulations (with thermal initialization) for a climate 1 ^∘C cooler and 30 % wetter (more precipitation) than the historical scenario. We call this new scenario colder and wetter (to be compared with the historical scenario, i.e., the results of the present study presented in the rest of Sect. 4). Results of this experiment are presented in Fig. 10 and Table 2. Because of the difference in climate forcing, the colder and wetter scenario produced a greater amount of cold- and warm-permafrost areas than the historical scenario, as presented on Fig. 10a. Figure 10b shows the proportion of the precipitation reaching the surface (rain + snow − snow sublimation) that produces runoff compared to evaporation for the Paiku catchment. Figure 10c aggregates over the whole catchment the distribution of such precipitation input to the ground between runoff and evaporation for both scenarios. In between them, it also includes the distribution associated with the steady-lake-level scenario of Fig. 9c, which is based on the hypotheses listed as bullet points in Sect. 4.3 (climate forcing of the historical scenario, same glacier contribution, only land runoff increases).

Figure 10 Sensitivity of the distribution between runoff and evaporation to climate. (a) Distribution of the different cryological states of the ground for the historical scenario (presented in Sect. 4.1–4.3) and for an alternative scenario where the climate is 1 ^∘C colder and brings 30 % more precipitation. (b) Runoff as a proportion of the precipitation input to the land (rainfall + snowfall − snow sublimation) for the different cryological states of the ground and for the two climatic scenarios. (c) Catchment-scale ratio between runoff and evaporation for (i) the historical scenario, (ii) for a steady lake level with the same glacier contribution (same as Fig. 9c), and (iii) for the colder and wetter scenario.

Table 2 Distribution between runoff and evaporation for the two scenarios.

^∗ Precipitation input is the input to the ground, counted as rainfall + snowfall − snow sublimation.

The historical scenario shows that cold-permafrost areas produce the highest proportion of runoff, which we attribute to the fact that the ground in these areas is, most of the time, frozen, turning a substantial part of the snow melt and rainfall into surface runoff. When considering grounds with a hydrologically active subsurface (warm permafrost, disappearing permafrost and no permafrost) in the historical scenario, the proportion of runoff increases slightly from warm permafrost to no permafrost. Such an evolution then corroborates the idea that the presence of permafrost tends to increase evaporation at the expense of runoff, as modeled by Sjöberg et al. (2021). Yet, for the colder and wetter scenario, runoff shows a regular decrease from cold to no permafrost, with a more pronounced trend than the historical scenario. Several factors can be at play in this transition and most likely involve (i) a different extent and altitudinal distribution for each cryological type of ground, (ii) an overall reduced intensity of evaporation due to cooler surface temperatures, (iii) a higher soil water content driven by higher precipitation and (iv) difference in the seasonal timings. Altogether, these processes substantially change the proportion of water that ends up as runoff water available for the lake, as highlighted by Fig. 10c.

5 Discussion

5.1 Limitations and potential of the approach

5.1.1 Data usage within the conceptual framework and data scarcity

Our approach relies on a variety of data with regard to their scientific focus (glaciers, ground, lake, atmosphere), their type (in situ observations, remotely sensed data, reanalysis data), their characteristics (pointwise data, distributed data, constant or with various time resolutions) and the way they interact with our models (model parameters, forcing data, validation data, result data in the case of the glacier runoff). Such a diversity arises from our goal to quantify both the ground thermo-hydrological regime and the different terms of the lake budget. This variety also makes it challenging to consistently merge these data into a unique framework. For example, our quantification of the glacier mass change reconstruction is made of two constant values for the study period (1975–2000 and 2000–2020), which limits the relevance of the comparison between the observed lake level variations and the simulated ones.

Yet, the lake level variations are the only hydrological observations available to evaluate the robustness of the runoff we compute. Therefore, we had to combine lake level observations with our precipitation forcing data and lake evaporation quantifications in a simple mass conservation calculation to derive the land runoff to the lake required to reproduce the level variations (red curve in Fig. 5d). In this regard, the sum of the glacier and land runoff we derive over the 40 years corresponds to 95 % of the required runoff to the lake, indicating that the magnitude of our reconstruction is correct. Year-to-year comparison is less accurate, and we suggest that this is the consequence of the aforementioned limitations and also of our modeling strategy, as detailed below.

A main limitation regarding our usage of the data is related to the limited amount of available field observations required to provide robust model parameterizing, climate forcing and in-depth validation of the simulations, both hydrologically and thermally. Regarding climatic forcing data, our AWS measurement offers sound observations to evaluate and adjust the ERA5 data processed with TopoSUB and downscaled with TopoSCALE. Yet, a period of observations longer than 2 years would have enabled more robust corrections and could have allowed us to perform a more advanced statistical downscaling approach, e.g., quantile mapping (Themeßl et al., 2011). As such, the spatiotemporal domain of relevance of these corrections is insufficient to correct data for the whole catchment and for the 40 years of simulations. Overall, considering the strong bias we observe in the raw ERA5 data (Fig. D1), these corrections do represent an important first-order improvement.

Additionally, in the absence of borehole data that would allow us to anchor our parameters into observations, we rely on gridded values designed for hydrological and/or land surface modeling (Sect. 3.2.4 and Appendix A). Because these values might be less reliable than field observations, we chose to average them over the catchment to derive some more robust values. Altogether, this scarcity of field observations is likely to bring significant uncertainties to our analysis. Future efforts should focus on acquiring additional data or developing validation methods based on remotely sensed observations.

5.1.2 Modeling strategy

A limitation in our study is that lateral water flows between land simulation units is ignored. By giving access to the timing of water transport across the catchment, water routing would allow us to investigate temporal hydrological patterns at a monthly or seasonal scale. Because we work at annual and decadal timescales, this limitation has limited consequences for our results. The main consequence is to ignore potential storage effects on the land that would delay the arrival of runoff to the lake. We suggest that it is possible that this limitation partly explains the limited match between computed and required runoff at the annual timescale (Fig. 5). Yet, our subdivision of the catchment based on the different cryological states of the ground allows us to work with hydrological units that are smaller than the catchment and thus present shorter hydrological response times to precipitation.

Additionally, our approach regarding the modeling of runoff is relatively simple, i.e., partition between subsurface and surface runoff based on comparison between the soil water content and field capacity and porosity, respectively. More complex approaches split runoff into more sophisticated categories such as Horton overland flow, Dunne overland flow, subsurface storm flow, etc. (e.g., Savenije, 2010; Gao et al., 2014; Mirus and Loague, 2013). However, over the last decade, the relevance of this type of partitioning between different types of runoff has been questioned (McDonnell, 2013; Gao et al., 2023). In the frame of our study, we find it important to distinguish between surface and subsurface runoff because they generate flows with very contrasted speeds. From a general perspective, this significant difference in flow velocities impacts the hydrological system as a whole (e.g., river discharge, evaporation) and has various consequences throughout the catchment, such as for the water availability for vegetation, erosion and sediment transport.

In the particular case of a cryo-hydrological study, separating surface from subsurface runoff is particularly relevant because both flows do not react in the same way to ground temperature changes. As such, we see our approach as a middle way that allows us to make this distinction based on simple hydrological considerations. Yet, we acknowledge that the classification and quantification of the different types of runoff represent a valuable direction for future investigation of catchment-scale cryo-hydrology in Tibet. Another potential improvement in our modeling approach could be to unravel evaporation from transpiration. However, since vegetation is extremely scarce in the Paiku catchment, which is largely dominated by barren lands, we suggest that this would not significantly affect our results. However, this limitation should be explored in future field and modeling studies.

Conversely, our approach also conveys several important advantages regarding our goal to describe and quantify the ground thermo-hydrological regime of the catchment. The use of TopoSUB enables us to produce results at a resolution of 100 m×100 m over an area of nearly 2400 km², with calculation costs 700 times lower than if each 100 m×100 m pixel was treated individually. Yet, thanks to the clustering method used to produce the forcing dataset (Sect. 3.2.2), the strong spatial variability of the physiography and its impact on the climate and incoming radiations are significant in the forcing data and have a major influence on the ground thermo-hydrological results, as exemplified by the strong spatial variability of ground temperatures (Fig. 7). Beyond elevation, other physiographic parameters such as aspect also influence the results. The mean values of 2 m deep temperature and evaporation over the 40 years for north-facing areas (averaged over the whole catchment and over the 40 years) are 1.3 ^∘C and 163 mm, while they reach 2.9 ^∘C and 197 mm for the south-facing ones. This strong dependence of modeled results on physiography highlights the necessity of taking this into account when modeling the thermo-hydrological regime of the ground in high-mountainous environments. Finally, our approach allows us to couple the physical processes governing both energy and water fluxes at the surface and subsurface and to highlight their interplay, as developed in Sect. 5.3.1.

5.2 Trends in the catchment and across the QTP

5.2.1 Lake hydrological budget and level variations

The total lake level change we simulate is a decrease of 4.11 m. This is qualitatively consistent with the overall observed trend. The mismatch with the observations is limited to a 1.04 m excess in the simulated level drop (Fig. 9a). Our reconstruction shows a decrease of 4.66 m from 1980 to 2007, which is an overestimation of the initial drop. Afterwards, while observations indicate a gradual slowdown of the lake level decrease, we simulate a stabilization followed by a slight increase (0.55 up between 2013 and 2019). The reason for the overall mismatch of 1.04 m can arise from bias (i) in the forcing data (and mainly in the precipitation) used for the land and lake simulations, (ii) in the glacier mass balance estimate, and/or (iii) in the quantification of hydrological processes for the land or for the lake (evaporation, runoff). On top of these potential biases, the difference in trends for the end of the simulation time can be influenced by (i) our estimates of glacier mass changes, which are made of two time averages (one for the 1980–2000 period and one for the 2000–2020 period) and therefore produce very smoothed glacier runoff values that cannot capture variations at the scale of a decade or less, and (ii) the absence of water routing that prevents us from accounting for delays in storage effects on the water supply from the land to the lake.

Additionally, our approach ignores potential water fluxes between the lake and a surrounding aquifer. This can be a possible reason for this mismatch. In the context of a decreasing lake level, an aquifer surrounding the lake can create an additional water inflow when the lake level passes below the piezometric level of the aquifer (Yechieli et al., 1995). We suggest that such an inflow could mitigate the lake level decrease and thus explain the missing water in our reconstruction (Fig. 6b). It could also explain the gradual stabilization of the lake level that our model does not reproduce. This flow is not part of our conceptual hydrological framework even though it likely exists in reality, especially since there is no permafrost near the lake (as we simulate it here), allowing for the existence of such an aquifer (Walvoord and Kurylyk, 2016). Groundwater has been identified as a potential contributor to lake level rise in other regions of the QTP (Lei et al., 2022). In the long run, lake–aquifer systems commonly follow oscillations of the net atmospheric flux of water (precipitation–evaporation) and of the runoff that forces its mass balance (Watras et al., 2014). During these oscillations, the lake can pump water from the aquifer or feed it, depending on the relative difference in the piezometric level between them (Almendinger, 1990; Liefert et al., 2018). Yet, this potential effect is difficult to account for, and its magnitude remains unclear. Therefore, the reasons for the mismatch between observed and simulated lake levels could also be connected to other aspects of our methodology, such as bias in the climatic forcing data and other shortcomings arising from the lack of field data, or hydrological processes, as developed in Sects. 5.1.1 and 5.1.2.

Our reconstruction of the lake budget is informative regarding the contributions of the different inputs and outputs. Regarding lake evaporation, our mean value of 870 ± 23 mm is close to the one modeled by Yang et al. (2016) with the Flake model for Lake Nam (832 ± 69 mm) for the period 1980–2014, but we do not report a significant increasing trend in our results. Yet for the same lake (Nam Co) and a similar period (1980–2016), Zhong et al. (2020) reported an average value of 1149 ± 71 mm (along with an increasing temporal trend) using the Penman formula (Penman, 1948), thus highlighting the potential dependence of the results on the methodology. In our results, direct precipitation to the lake represents 40 % of the inputs, followed by glacial runoff (35 %) and land runoff (25 %). Glaciers are therefore a particularly important contributor to the runoff towards the lake (60 % of the total runoff vs. 40 % for land runoff), which is in contrast to the results from Biskop et al. (2016), who calculated that the runoff input to Lake Paiku was dominated by land runoff (70 % vs. 30 % for the glacier contribution). Here again, these differences likely arise from important differences in input data and methodologies to quantify the different hydrological processes (evaporation, runoff, snow and glacier melt). Yao et al. (2018) reported that, at the QTP scale, the balance between precipitation and evaporation (over land and lake) was dominant over glacier melt for understanding both lake storage increases and decreases. Our reconstruction does not give us access to significant temporal variations of the glacier contribution, but the abovementioned proportions in the contributions to the lake (40 %, 35 % and 25 %) show that the glacier contribution does not dominate the input terms. At the catchment scale, these proportions can vary significantly depending on the glacier coverage. For Lake Selin, Zhou et al. (2015) reported that runoff towards the lake, evaporation from the lake and on-lake precipitation altogether explained 90 % of the lake storage variations for the 2003–2012 period. The catchment of Lake Selin has a very limited glacier coverage, corresponding to 0,63 % of its area (Lei et al., 2013), compared to the Paiku (5 %).

5.2.2 Permafrost and ground temperature trends

Our results indicate that permafrost coverage in the Paiku catchment evolved from 27 % to 22 % of the land area during the simulated period. Such a coverage corresponds to sporadic permafrost (10 %–50 % of the area) and is consistent with recent large-scale estimates of permafrost in the Northern Hemisphere (Obu et al., 2019) and across the QTP (Zou et al., 2017; Ran et al., 2018). This decrease corresponds to a 19 % shrinkage of the 1980 permafrost area, which is higher than the 9 % reported by Gao et al. (2018), a value determined by catchment-scale numerical modeling in the upper Heihe catchment (northeastern QTP) over a similar period. It is also slightly higher than the 13 % decrease modeled by Wang and Gao (2022) from 1971 to 2015 for the Qinghai Lake catchment using a similar approach. Yet, it is smaller than the 34 % loss modeled by Qin et al. (2017) from 1981 to 2015 for the Yellow River Source Region (YRSR, northeastern QTP).

Active-layer (AL) evolution is contrasting throughout the catchment, and a deepening signal is only visible for the locations with limited evaporation (<150 mm yr⁻¹). Given the strong drive of summer climate in relation to ALT, this overall lack of a deepening trend highlights how evaporation can act as an energy intake at the surface (K. Yang et al., 2014), limiting the surface and subsurface heat fluxes and thus AL deepening. In this regard, our results fall in line with the conclusions of Fisher et al. (2016) when observing evapotranspiration and ALTs in boreal forests and also confirm the modeling experiments of Zhang et al. (2021b) on permafrost wetting in arid regions of the QTP. Besides, the lack of an overall deepening trend is consistent with observations from Luo et al. (2018) in the YRSR over the last decade and with the modeled AL from Zhang et al. (2019) at the scale of the QTP for the last 40 years. Where evaporation is limited, we report an AL-deepening trend of 2.7 cm per decade, which is smaller than the 4.8 cm per decade trend modeled by Song et al. (2020) for the YRSR for the same period and smaller than the 4.3 cm modeled by Gao et al. (2018) in the upper Heihe catchment. Yet, it is comparable to the 2 cm per decade value modeled by Wang and Gao (2022) for the Qinghai Lake catchment from 1971 to 2015. Connections between AL deepening and evaporation are discussed in Sect. 5.3.1.

In no-permafrost areas, our simulations show that the thickness of seasonally frozen ground shrinks at a rate of 6.8 cm per decade. This rate is faster than the rate of 3.1 cm per decade quantified by Qin et al. (2018) using the Stefan solution for the YRSR (1961–2016) and faster than the 3.2 cm per decade modeled by Gao et al. (2018, Heihe catchment). However, it is similar to the 6 cm per decade rate modeled by Wang and Gao (2022) in the Qinghai Lake catchment from 1971 to 2015 and smaller than the 12 cm per decade modeled by Qin et al. (2017) for the YRSR (1981–2015). All these values fall within the wide range of 3 to 29 cm per decade reported by C. Wang et al. (2020) when studying seasonally frozen ground over the whole QTP with in situ observations. Regarding timing, we report a decreasing trend of 5.3 d for frozen conditions (70 cm deep) per decade. which is consistent with the decrease of 6.7 d per decade reported by C. Wang et al. (2020) just below the surface.

Regarding the timing of seasonal ground thaw, our results highlight that the increase in the duration of the seasonal ground thaw (at 70 cm) is mostly driven by a progressive delay of the end date of the thaw period. This result contrasts with those from Song et al. (2020) for the same period in the YRSR; they also modeled an increase of the seasonal thaw (at a 2 cm depth), although this was driven by an advancing trend of the start date of the seasonal thaw.

Our warming trends at 4 m depth for permafrost areas is 0.1 ^∘C per decade, which is substantially smaller than the 0.43 ^∘C per decade observed at this depth between 1996 and 2006 in permafrost boreholes along the Qinghai–Tibetan Highway in the northeast of the QTP (Wu and Zhang, 2008). Zhang et al. (2019) reported a 0.13 ^∘C warming per decade of the permafrost top during winter, which is consistent with the trend of 0.14 ^∘C per decade we observe at 2 m depth (mean AL between 1.4 and 1.7 m in our simulations) for the months of December, January and February.

5.2.3 Evaporation and runoff trends

Our results are characterized by (i) an increase in both evaporation and runoff (Fig. 9a and b), mainly driven by an increase in precipitation (Fig. 3 bottom); (ii) a runoff$/$(runoff + evaporation) ratio exhibiting an increasing trend as a result of ground warming and permafrost disappearance, both of which enable more subsurface runoff over time (Figs. 9c and 10d); and (iii) an increase in the proportion of liquid water in the ground compared to ice (Fig. 9d). Regarding all these points, our results show good consistency with the evolution reported by Gao et al. (2018) for the upper Heihe catchment (northeastern QTP) using a similar approach for a comparable period (1971–2013). The increasing trends in evaporation and runoff they report for the thawing season (dominant period for both processes) are comparable with the annual values we report: +10.0 mm per decade for evaporation (our study: +10.1 mm per decade) and +3.3 mm per decade for runoff (our study: +4.8 mm per decade). Similar evolutions are also reported by Wang and Gao (2022) for the Qinghai Lake catchment and by Qin et al. (2017) for the YRSR (1981–2015). These increases in runoff (especially surface runoff) are likely to have an influence on sediment transport. For instance, Li et al. (2021) showed that the current precipitation increase over High Mountain Asia is driving a runoff increase, which contributes to a significant rise in fluvial sediment fluxes. Regarding differences, Qin et al. (2017) modeled a stronger evaporation increase (14.3 mm per decade) linked to a decreasing runoff coefficient. Similarly to Li et al. (2019), we see that an important part of snow melt (49 %) infiltrates the ground and later contributes to runoff and evaporation.

5.3 Cryo-hydrological couplings at the catchment scale and implications for lake level variations

5.3.1 Interdependence of thermal and hydrological variables

Our simulation results enable us to explore the interplay between the fluxes of energy and water at the surface and subsurface. In this regard, we tested the correlation of evaporation with the proportions of liquid water and total water in the ground for cold and warm permafrost, as well as the correlation between evaporation and the duration of seasonal thaw at 70 cm depth (Fig. 11a and b). For permafrost areas (cold permafrost and warm permafrost), evaporation shows a strong correlation with the seasonal distribution between liquid and frozen water, similarly to previous modeling works for the region (Cuo et al., 2015). As such, this correlation suggests that the intensity of seasonal ground thaw plays a role in enabling higher or lower evaporative fluxes. This is likely due to cold surface temperatures strongly reducing water loss from the surface and because moisture delivery to the surface is inhibited when the ground is frozen. We suggest that this dependence is particularly important in the Paiku catchment because evaporation is strong (88 % of the precipitation input to the surface evaporates on average) and because frozen water is the dominant form of water in the ground in permafrost areas (Fig. 11a; the calculation includes the first 2 m below the surface).

Figure 11 Thermo-hydrological couplings. (a) Annual evaporation vs. annual mean of the liquid$/$total water ratio over the first 2 m of ground, averaged for simulations corresponding to cold permafrost and warm permafrost (one dot per year for each permafrost category). (b) Annual evaporation vs. duration of seasonal thaw at a 70 cm depth averaged for simulations corresponding to locations without permafrost (one dot per year). (c) Active-layer deepening trend vs. mean evaporation over the 40 years for each simulation corresponding to warm permafrost (here, one dot corresponds to one TopoSUB point). The color of the dots shows the precipitations averaged over the 40 years for each simulation. The linear regression excludes simulations exhibiting annual evaporation higher than 240 mm. (d) Annual subsurface runoff vs. annual 2 m deep temperature averaged for simulations corresponding to locations with disappearing permafrost (one dot per year). The color of the dot indicates the year of the simulation.

Similarly, evaporation in no-permafrost areas shows a significant correlation with the duration of the seasonal thaw (Fig. 11b). We suggest that this result arises from the fact that frozen ground limits the evaporative fluxes, and thus, years during which the subsurface seasonal thaw is shorter are associated with reduced evaporative fluxes. We also tested the relationship between the linear trend of active-layer deepening and the mean evaporation (over the 40 years of simulation) for warm-permafrost areas (Fig. 11c). Thus, this graph does not present annual values, and one point corresponds to one of the 92 TopoSUB points classified as warm permafrost (values averaged over the 40 years). The graph highlights that TopoSUB points showing an active-layer (AL) deepening trend are associated with low evaporation and precipitation. From there, TopoSUB points with stronger evaporation show no deepening trend or even a shrinkage of the AL. This relationship is contradicted by the highest level of evaporation (>240 mm yr⁻¹) observed for warm permafrost, for which AL deepening is observed again (dark-blue points of the graph). These TopoSUB points with the highest levels of evaporation also correspond to those receiving the largest amount of precipitation.

Runoff also shows a strong connection with the ground thermal regime (Fig. 11d). At the beginning of the simulation, years with an average 2 m deep temperature below 0 ^∘C are associated with limited subsurface runoff (<5 mm yr⁻¹). Over the years, as the ground warms up and permafrost disappears, subsurface runoff increases and can reach 20 to 45 mm yr⁻¹. This result is consistent with the increased subsurface connectivity expected when permafrost thaws (Kurylyk et al., 2014; Gao et al., 2021), which has been both observed (Niu et al., 2016) and modeled (Lamontagne-Hallé et al., 2018; Huang et al., 2020; Gao et al., 2018). We suggest that these substantial changes in subsurface runoff, associated with changes in the ground temperature in Fig. 11d, support the hypothesis of a modification in the hydrological pathways as permafrost thaws.

Altogether, these results suggest a dependence of key variables quantifying the catchment hydrological balance (evaporation, runoff) on the seasonal characteristics and interannual trends of the ground thermal regime (temperature, liquid vs. frozen water content). Similarly to previous studies (Ding et al., 2020; Wang and Gao, 2022), we think these results advocate for the necessity of coupling thermal and hydrological modeling to improve our ability to understand and quantify changes in the hydrological balance of high-mountain catchments. To our best knowledge, along with Gao et al. (2022), our study represents to date the most complete effort to include the variety of coupled climatological, surface and subsurface processes characterizing the climate, hydrology and ground thermal regime of high-mountain catchments in Tibet at a small scale with a high spatial resolution.

5.3.2 Influence of the ground thermal regime on the distribution between runoff and evaporation

Our results indicate that evaporation is particularly strong in the Paiku catchment. Over the 40 years of simulation, 10 % of the total precipitation is converted into runoff, and the rest of the water is either directly returned to the atmosphere from the snowpack via snow sublimation or from the ground surface via evaporation. Comparatively, Gao et al. (2018) observed and modeled a ratio of around 35 % for the Heihe catchment; Qin et al. (2017) reported an average ratio of 33 % for the YRSR, and Li et al. (2014) reported a ratio of 83 % for the Qugaqie catchment (central QTP) but modeling hydrological fluxes only.

Our sensitivity test on evaporation and runoff for slightly different climates (Sect. 4.4) highlights the fact that the role of permafrost regarding the runoff–evaporation distribution is a complex question, as has already been discussed in the literature (e.g., Bring et al., 2016). Some studies have suggested that landscape-scale permafrost thaw would trigger more evaporation (Walvoord and Kurylyk, 2016). This phenomenon was modeled by Wang et al. (2018) in the upper Heihe River catchment, for which they reported that the thickening of the active layer increased the ground storage capacity and led to a decrease in runoff and an increase in evapotranspiration. G. Wang et al. (2020) also reported that permafrost thawing accelerated evapotranspiration (1961–2014).

Conversely, Zhang et al. (2003) and Carey and Woo (1999) reported that shallow frozen-ground conditions (such as a shallow active layer) maintain higher water contents close to the surface, promoting higher evaporation. Sjöberg et al. (2021) modeled this phenomenon with a fully coupled cryo-hydrological model including surface energy balance calculation. They modeled a slope with a simplified geometry in 2D for different permafrost coverages. They found that hillslopes with continuous permafrost have rates of evapotranspiration that are twice as high compared to hillslopes with no permafrost.

As such, the interplay between the runoff–evaporation distribution and the ground thermal regime in areas where permafrost coverage shows a spatiotemporal variability is difficult to apprehend (Fig. 10). This complexity is most likely to be due to a strong sensitivity to the drainage conditions (fast flows of steep mountain environments vs. slow flows of lowland catchments) and to the climate setting, both at the annual scale (arid regions vs. wet regions) and at the seasonal timescale (relative timing of temperature variations, rainfall, snowfall, snow melt and ground freeze–thaw).

Because it can both promote evaporation or runoff depending on the setting, the ground thermal regime of the catchment seems to have the possibility to create a positive feedback towards both lake level decreases or increases. Further studies should therefore focus on comparing the thermo-hydrological regime of different Tibetan catchments with contrasting lake level changes and permafrost coverage to test to which extent these differences can contribute to explaining the spatial patterns of lake level changes across the QTP.

5.4 Implications for lake level changes

At the scale of the Paiku catchment and with regard to lake level variations, the results we present highlight the following:

• The sums of the direct precipitation in the lake, the land runoff and the glacier runoff are not enough to compensate for the lake evaporation over the study period, hence driving the observed lake level decrease.

• Long-term hydrological trends in the catchment are led by trends in climate, and precipitation increase, jointly with glacier melt, provides enough water to drive a concomitant increase in runoff and evaporation.

• Ground thermal changes increase the distribution of liquid vs. frozen water in the ground and the duration of seasonal thaw; correlations suggest that these modifications increase evaporation. The warming of the ground is also related to the increase in subsurface runoff towards the lake.

• Ground warming and permafrost thawing promote subsurface runoff over time, contributing to an increase in the runoff / evaporation ratio of the catchment.

• Over the 40 years we studied, the presence of permafrost seems to promote evaporation at the expense of runoff. Yet, this trend appears to be climate dependent, and the cryological state of the ground might shift the runoff–evaporation distribution in the other direction under colder and wetter climates.

At the scale of the QTP, these results have several implications. First, a better understanding of the recent and future lake level variations will come with better knowledge of spatial patterns and temporal trends in precipitation. Second, climate changes are modifying the ground thermal regime of Tibetan catchments. Ground warming may lead to active-layer deepening, permafrost disappearance and/or changes in the seasonal freeze–thaw cycles, affecting evaporation, runoff volumes and pathways and, overall, changing the hydrological functioning of Tibetan catchments (and the water flow provided to the lakes). Finally, the effect of permafrost on the distribution between evaporation and runoff seems to be dependent on the climate settings and the permafrost coverage of the catchment. Further studies should investigate this phenomenon and how it might contribute to explaining the contrasting lake level evolutions across the QTP.

6 Conclusion

Our study quantifies the different terms of the Lake Paiku budget over the past 40 years. Direct precipitation to the lake represents 40 % of the inputs, followed by glacial runoff (35 %) and land runoff (25 %). Glaciers are therefore a particularly important contributor to the runoff towards the lake.

We also confirm that the ground of the Paiku catchment presents different types of cryological states, from seasonally frozen ground to permafrost. Permafrost coverage shrinks from 27 % to 22 % of the land area of the catchment from the 1980s to the 2010s (19 % loss compared to the 1980 permafrost area). The whole catchment warms up at a rate of 0.17 ^∘C per decade (2 m deep), with a substantial elevation-dependent variability. This warming is concomitant with an increase in the duration of the seasonal thaw, mainly supported by a progressive delay of the end date of the thaw period. Where permafrost is present, active-layer deepening is only observed where evaporation is relatively low (<150 mm yr⁻¹).

Over the simulation period, we also report an increase in evaporation (+10.1 mm per decade) and in surface and subsurface runoff (+1.3 and +3.5 mm per decade, respectively). Together, this leads towards an increase in the runoff$/$(runoff + evaporation) ratio of +1.2 % per decade. Our results also highlight the strong interdependence between the ground thermal and hydrological regimes and the necessity to jointly represent them to accurately quantify evaporation and runoff in this type of environment.

Over the last 40 years, the presence of permafrost seems to have promoted evaporation at the expense of runoff. Yet, this trend appears to be climate dependent, and the cryological state of the ground might shift the runoff–evaporation distribution in the other directions under colder and wetter climates. Further studies should investigate this phenomenon and how it might contribute to explaining the contrasting lake level evolutions across the QTP.

Appendix A: Model parameters

Table A1 Parameters of the model.

Appendix B: Geological map of the catchment

Figure B1 Geology of the catchment. Modified from Aoya et al. (2005). The red contour indicates the limits of the Paiku catchment.

Appendix C: TopoSUB subsampling of the catchment

Figure C1 Application of the TopoSUB clustering method (Fiddes and Gruber, 2012) in the Paiku catchment. (a) Number of the TopoSUB points. Strong color changes reflect the footprint of the eight ERA5 pixels that the catchment intersects. Small color changes within a given zone show the distribution of the 50 TopoSUB points covering each tile (Sect. 3.2.2). (b) Topographic map reconstructed using the TopoSUB approach.

Appendix D: Evaluation of forcing data

Figure D1 Comparison between the AWS data and the model forcing data downscaled from ERA5 with the TopoSCALE and TopoSUB approaches. Based on the AWS data, a monthly correction factor is applied to the downscaled data so that monthly data match the observed period for each variable (methodological details in Sect. 3.2.2).