Seasonal thermal regime and climatic trends in lakes of the Tibetan highlands

The hydrology of the lake-rich Tibetan Plateau is important for the global climate, yet little is known about the thermal regime of Tibetan lakes due to scant data. We (i) investigated the characteristic seasonal temperature patterns and recent trends in the thermal and stratification regimes of lakes on the Tibetan Plateau and (ii) tested the performance of the one-dimensional lake parameterization scheme FLake for the Tibetan lake system. For this purpose, we combined 3 years of in situ lake temperature measurements, several decades of satellite observations, and the global reanalysis data. We chose the two largest freshwater Tibetan lakes, Ngoring and Gyaring, as study sites. The lake model FLake faithfully reproduced the specific features of the high-altitude lakes and was subsequently applied to reconstruct the vertically resolved heat transport in both lakes during the last 4 decades. The model suggested that Ngoring and Gyaring were ice-covered for about 6 months and stratified in summer for about 4 months per year with a short spring overturn and a longer autumn overturn. In summer the surface mixed boundary layer extended to 6–8 m of depth and was about 20 % shallower in the more turbid Gyaring. The thermal regime of the transparent Ngoring responded more strongly to atmospheric forcing than Gyaring, where the higher turbidity damped the response. According to the reanalysis data, air temperatures and humidity have increased, whereas solar radiation has decreased, since the 1970s. Surprisingly, the modeled mean lake temperatures did not change, nor did the phenology of the ice cover or stratification. Lake surface temperatures in summer increased only marginally. The reason is that the increase in air temperature was offset by the decrease in radiation, probably due to increasing humidity. This study demonstrates that air temperature trends are not directly coupled to lake temperatures and underscores the importance of shortwave radiation for the thermal regime of high-altitude lakes.


Introduction
The hydrological regime of the Tibetan highlands is extremely complex and highly sensitive to climate changes (Gu et al., 2005;Ma et al., 2012;Yang et al., 2014).The Tibetan Plateau plays a crucial role in the global water cycle because it directly affects the monsoon system and is the origin of major Asian rivers (Su et al., 2015) including the Yellow, Yangtze, Mekong, Salween, Brahmaputra (Yarlung Tsangpo), and Indus (Fig. 1).The plateau contains thousands of lakes, the hydrological regime of which is closely connected to glaciers, permafrost, and influent and effluent (Zheng, 2011;Zhang et al., 2013).The particular climatic environment of the Tibetan Plateau, with low air pressure and intense solar radiation, creates a unique land-atmosphere interaction (Ma et al., 2009) where the lake system is a crucial component of the regional heat and mass balance (Liu et al., 2009;Gerken et al., 2013;Li et al., 2015;Wen et al., 2015;Dai et al., 2016).Therefore, changes in the regional water budget of the Tibetan Plateau driven by global warming are of key importance for climatic changes on continental scales.Recent evidence indicates that the number of lakes on the Tibetan Plateau is increasing, highlighting global warming effects and local anthropogenic activity, such as infrastructure construction on permafrost (Cheng and Wu, 2007;Liu et al., 2009;Yang et al., 2010).Hence, to estimate how lakes may affect the regional climate, it is necessary to quantify the seasonal thermal dynamics of lakes and ponds with particular reference to the vertical thermal stratification, the ice regime, and the lake interactions with the lower atmosphere.
Due to the importance of lakes for understanding the global-scale changes on the Tibetan Plateau, several studies have been performed in recent years on the lake dynamics and lake-atmosphere interactions.The amount of evaporation and heat flux and their change was revealed in lakes in the Tibetan Plateau during the pre-and post-monsoon periods (Haginoya et al., 2009;Xu et al., 2009).Biermann et al. (2014) measured the turbulent fluxes over wet grassland and a shallow lake with the eddy covariance method at the shoreline of the Nam Co Lake basin.Li et al. (2015) recently studied the characteristics of the energy flux and the radiation balance over Ngoring Lake.With the development of lake schemes and atmospheric models that include the lake surface (Sun et al., 2007), other studies have investigated the effects of the Tibetan Plateau lakes on the local climate (Li et al., 2009;Wen et al., 2015;Yang et al., 2015).Nevertheless, compared to other lake-rich regions of the world, like Fennoscandia, northern Canada, or the eastern Siberian tundra, the lakes of the Tibetan Plateau are poorly studied, despite their extraordinary thermal regimes due to strong seasonal and synoptic variations.
A major problem in resolving the energy and moisture budget of high-altitude lakes is the lack of long-term observations.Observational data on the physics of Tibetan Plateau lakes are scarce and mostly confined to remote-sensing data on lake surface characteristics, such as lake surface temperature (Lin et al., 2011;Zhang et al., 2014).Observations of air-lake fluxes using eddy covariance methods over Tibetan Plateau lakes are too fragmentary to estimate the seasonal variations (Biermann et al., 2014;Li et al., 2015).Furthermore, lack of knowledge of the heat transport within the lake water column, including the seasonal formation of thermal stratification and the thermal dynamics under the ice cover, have limited our understanding of lake-atmosphere interactions on climatic scales.
While the first brief reports were recently presented (Wang et al., 2014;Wen et al., 2016) on the mixing conditions and vertical heat transport within the lake water bodies, there are no observational data covering the annual variability or the long-term trends in thermal conditions and stratification in Tibetan lakes.No regular lake monitoring was performed in this area in the past.In this situation, the modeling of decadal variability in lake dynamics can provide insight into the ongoing changes on seasonal to climatic time scales, provided that the model performance is tested against available observational data for shorter periods.On the other hand, the validation of models for parameterizing lakes in weather prediction and climate studies is of significant value, particularly because it is crucial to properly account for the heat budget of one of the largest lake systems worldwide to adequately model heat and mass transport on a planetary level.
In view of these challenges, the present study aimed to characterize the seasonal stratification regime and ice cover formation in lakes of the Tibetan Plateau using the two largest freshwater lakes of the region as examples; test the ability of the lake model FLake (Mironov, 2008;Kirillin et al., 2011) to simulate the main features of the lake thermodynamics in high-altitude conditions as a potential tool to parameterize the Tibetan lake system in regional climate modeling; and quantify the mean trends in the thermal and ice regime of Tibetan lakes during the last half-century to estimate the climate-driven changes in the seasonal stratification pattern.
For this purpose, we combine 3-year observations on the surface heat budget from Ngoring Lake (Wen et al., 2016) with the surface temperatures from the AVHRR sea surface temperature (SST) product (Casey et al., 2010), the one-dimensional lake temperature and mixing model FLake, and the long-term forcing from the ERA-Interim and NCEP/NCAR global atmospheric reanalyses (Kalnay et al., 1996;Dee et al., 2011).In situ and satellite data on the surface temperatures are primarily used to validate the model Therefore, a comparison of the model results allows us to not only quantify the seasonal and climatic characteristics of the mixing regime, but also to estimate the response of the lake thermal structure to variations in the subsurface radiation regime and total depth of the water column.

Study sites
Ngoring (Big Blue Lake) and Gyaring (Big Gray Lake), sometimes referred as the "twin lakes", are two large freshwater lakes in the source region of the Yellow River (Huang He) on the eastern Tibetan Plateau (Fig. 1).The lakes are located at 34.5-35.5 Ngoring Lake has a surface area of 610 km 2 .The mean and maximum depths are 17 and 32 m, respectively.The lake is oligotrophic; fish in the lake are rare and aquatic plants grow only in the riparian areas.The reported lake water transparency (Secchi depth) does not exceed 3 m (Kar, 2014).Gyaring Lake has a similar surface area (526 km 2 ) but is shallower than Ngoring Lake with an average depth of 8.9 m and a maximum depth of 13.1 m.The lake water is more turbid than in Ngoring, with a Secchi depth of only 1 m (Kar, 2014).Due to the close geographical location and similar surface areas, both lakes are subject to virtually the same atmospheric forcing, but differ in their depths and water transparency, which are two major factors determining radiative heat transport within the water column and the formation of thermal stratification.Also, both lakes are fresh with a salt content of < 1 g L −1 (Kar, 2014), so the effects of salt on density stratification can be neglected in the modeling.

Lake temperature data
The lake-atmosphere heat exchange and the characteristics of the air-lake boundary layer were systematically observed over Ngoring Lake during the summer open water seasons of 2011-2013 (Li et al., 2015(Li et al., , 2016;;Wen et al., 2015Wen et al., , 2016)).The bulk surface temperatures obtained in these observations comprise the primary dataset for verifying the lake model in this study.The bulk temperatures of the lake surface water were recorded by 109 L probes (Campbell Scientific, Inc., Logan, UT, USA) at water depths of 0.05, 0.20, 0.40, and 0.60 m during the observation period.All sensors in the 0.6 m thick surface layer recorded essentially the same temperatures within the probes' accuracy, so the four time series were vertically averaged for the subsequent analysis.The incoming and outgoing shortwave and longwave radiation was measured with a net radiometer (CNR 1 and CNR 4; Kipp and Zonen, Delft, The Netherlands) 1.5 m above the lake.

Modeling
The modeling was performed with the lake temperature and mixing model FLake (Mironov, 2008;Kirillin et al., 2011), a highly parameterized model of the lake thermal regime specifically designed to parameterize inland waters in regional climate models and numerical weather predictions.The integrated approach implemented in FLake allows for the combination of high computational efficiency with a realistic representation of the major physics behind turbulent and diffusive heat exchange in the stratified water column.As a result, FLake is widely used for lake representation in the land schemes of regional climate models, being implemented, among others, into the surface schemes of the Weather Research and Forecasting Model (WRF; Mallard et al., 2014), HTESSEL (ECMWF; Dutra et al., 2010), SUR-FEX (Meteo France; Salgado and Le Moigne, 2010), and JULES (UK Met Office; Rooney and Jones, 2010).Thanks to its robustness and computational efficiency, FLake has become a standard choice in climate studies involving the feedback between inland waters and the atmosphere and is used operationally in the NWP models of the German Weather Service (DWD), the European Centre for Medium-Range Weather Forecasts (ECMWF), the UK Met Office, the Swedish Meteorological and Hydrological Institute (SMHI), the Finnish Meteorological Institute (FMI), and others.
For completeness, we provide a short description of the model concept following Kirillin (2010).For further details, the reader is referred to the model manual (Mironov, 2008).FLake is a one-dimensional model calculating the evolution of the horizontally averaged temperature in lakes.The basic principle underlying the model consists of splitting the water column into horizontal layers of variable depth, using physically motivated assumptions about their properties, and subsequently integrating the governing equations over these layers.In this sense, FLake belongs to the family of bulk models, of which the mixed layer models of the oceanic surface boundary layer (Niiler and Kraus, 1979) are the most famous examples.In contrast to the traditional mixed layer models, FLake additionally involves the "thermocline self-similarity" hypothesis to describe the lower stratified layer.The idea consists of adopting the thermocline thickness, h(t), and the temperature jump across the thermocline, T (t), as universal scales for the temperature profile, T (t, z), (Fig. 2).Using these scales, the dimensionless temperature, ϑ, and the dimensionless vertical coordinate, ζ , can be introduced as where T s is the temperature in the surface mixed layer and h is the mixed layer depth.The self-similarity of the temperature profile implies the universality of the function ϑ(ζ ), or in other words, a universal shape of the temperature profile within the stratified layer of all lakes.Several analytical expressions for ϑ(ζ ) have been proposed by various authors that all represent a similar thermocline shape supported by numerous observations (see the comprehensive review by Kirillin, 2002).Assuming a two-layer lake temperature structure with a homogeneous upper layer and a self-similar stratified layer beneath, the vertical temperature profile across the entire water column can be described as (Fig. 2) The temperature profile represented by Eq. ( 2) does not depend explicitly on the vertical coordinate, which allows the heat transport equation to be integrated across the layers 0 < z < h and h < z < (h + h), turning the partial differential equation for T (t, z) into a set of two ordinary differential equations (ODEs).The same self-similarity concept is used to describe the temperature structure of the thermally active upper layer of the bottom sediments (Golosov and Kirillin, 2010) and of the ice cover (Mironov et al., 2012) so that the model, strictly speaking, is four-layered with a unique self-similarity function ϑ(ζ ) for every layer: ice (if it exists), the mixed layer, the stratified layer, and the thermally active sediment.The integration of the turbulent kinetic energy (TKE) equation across the upper mixed layer (Niiler and Kraus, 1979) yields an additional ODE for the evolution of the mixed layer depth, h(t).The derivation procedure and the exact ODE expressions are given by Mironov (2008).The surface boundary conditions enter the reformulated problem as the time-dependent functions for the surface heat and momentum fluxes; the latter is used in the TKE equation to estimate the impact of the wind-generated mixing on the mixed layer h(t) evolution.The shortwave solar radiation is not included in the boundary condition for the surface heat flux, but is treated as a separate heat source absorbed internally within the upper water column by the one-band exponential law for radiation decay.A separate model block calculates the surface heat and momentum fluxes from the standard meteorological variables and the water surface temperature at the previous time step (Mironov, 2008).
. Schematic of a vertical temperature profile in a thermally stratified lake and the corresponding self-similarity coordinates.See Eqs. ( 1)-( 2) for the meaning of the symbols.

Lake temperatures from Pathfinder
As an additional source of the lake surface temperature estimations, we used the AVHRR Pathfinder Version 5.2 (PFV52) remote-sensing data obtained from the US National Oceanographic Data Center and GHRSST (http://pathfinder.nodc.noaa.gov).The PFV52 data are an updated version of the Pathfinder Version 5.0 and 5.1 collections described in Casey et al. (2010).Among the available lake surface temperature products, the AVHRR Pathfinder data have the longest time coverage  and a 4 km spatial resolution, distinguishing both Gyaring and Ngoring as water surfaces.
The lake mask of the PFV52 includes 25 grid points for Gyaring, 5 of which are "open water" points not adjoining the lake shoreline, and 31 points for Ngoring, 8 of which are "open water" points.In the further analysis, only spatial averages of the water temperatures from the open water points were used (Fig. 1) to reduce the potential influence of the surrounding land on the temperature readings.

Long-term model forcing from reanalysis data
In the absence of direct climatic observations for the study area, we adopted for two widely used global reanalyses of atmospheric fields, the NCEP/NCAR Reanalysis 1 (previously the NCEP/NCAR Reanalysis 40; Kalnay et al., 1996) and the ERA-Interim reanalysis produced by the European Centre for Medium-Range Weather Forecasts (ECMWF; Dee et al., 2011).Both datasets rely on different atmospheric models and data assimilation schemes.Accordingly, we first compared the surface layer meteorology of the Tibetan Plateau from the two reanalysis datasets, and then we investigated how the surface meteorology affected the modeled lake

Model setup and validation
We avoided tuning the model-specific parameters and algorithms to the observed lake data, with the exception of the lake heat content at the initial stage of the open water season.The water temperature immediately after the ice breakup was set to 4 • C independently of the modeled water temperatures during the preceding ice-covered period.This amendment aimed to account for the heating of the water column by solar radiation during the ice-covered period, which is not considered in the current version of FLake.Neglecting the heating by solar radiation penetrating the ice cover did not introduce large errors in FLake applications to low-altitude, seasonally ice-covered lakes (Kirillin, 2010;Bernhardt et al., 2012) where the incoming solar radiation in winter is lower and the relatively thick snow cover damps the radiative heating.However, in the high-altitude, low-latitude lakes of the Tibetan Plateau, the radiation penetrating the ice appears to be an important component of the lake heat budget, so the water temperatures directly after the ice breakup are close to the maximum density value (Wen et al., 2016).
The water transparency characteristics used in the model configuration were adopted from Kar (2014).The light extinction coefficient γ was calculated as the reciprocal of the Secchi depth h S , γ = C γ h −1 S , assuming C γ = 1.8, an approximate mid-value of the range reported in previous studies (see reviews in Kirillin and Shatwell, 2016;Shatwell et al., 2016).
The model performance was assessed by comparing the daily averages of the modeled and observed (in situ and satellite) surface temperature (T s ) between May and October.The satellite data contained mostly only one data point on a given day, which may increase error.Some outliers existed, usually under the ice cover or at subzero skin temperatures.All T s ≤ 0 were excluded from the modeled and observed data, i.e., also the whole ice-covered period.Different aspects of the model performance were examined including model bias (Eq.3), centered root mean square error (RMSE c ; Eq. 4), normalized standard deviation (σ norm ; Eq. 5), and normalized root mean square error (I 2 and I 3 ; Eqs. 6 and 7): 5) Here, m i and o i are the modeled and the observed values, respectively, and m and o are their means.

Statistics and trend analysis
The lakes were defined as stratified when the difference between the surface and bottom temperature was T s −T b > 0.5 K for summer stratification and T s − T b < −0.5 K or when the ice cover thickness was h ice > 0 for winter stratification.The stratification duration as well as the stratification start and end dates, given in day of the year (doy), refer to the longest uninterrupted stratification event in each year or season.The lakes were defined as mixed when they were not stratified.Spring and autumn overturn were defined as the longest uninterrupted periods of mixing following winter and summer stratification, respectively.In simulations when the lakes froze in winter, they always remained frozen until thawing in spring.Therefore, the ice cover duration as well as the freeze and thaw dates refer to both the longest uninterrupted period of ice cover and the total number of days of ice cover in each winter season.The seasons are defined by the typical hydrological and ice regime; winter: January-March, spring: April-June, summer: July-September, and autumn: October-December.The trends in the data were assessed using linear regressions.The homogeneity of the variance and the normality of residuals were assessed by inspecting the plots of residuals vs. fitted values and normal quantile-quantile plots.Leverage was assessed using Cook's distances.All statistical analyses were performed with R version 3.3.0(The R Project for Statistical Computing, Core Team, Vienna, Austria, 2016).The satellite data on lake surface temperatures, despite the coverage of several decades, were found inappropriate for reliable trend estimation due to numerous data gaps of several months, exceeding the halfperiod of the seasonal variations (Fig. 3)

Model performance
The model simulated the lake surface temperatures during the ice-free period well (Fig. 4; Table 1).The temperature   bias of the model was negligible (< 1 % of the range of temperature variability) when forced by both the NCEP and ERA reanalyses (Fig. 5).The correlation between the modeled and observed values was notably high for the ERA input: r = 0.96, or 92 % of the variance was captured by the model.The NCEP forcing provided slightly worse results: r = 0.93, or 86 % of the correctly simulated variance.The normalized standard deviation (σ norm ) was also better for the ERA forcing, differing from 1 by less than 0.3 %, whereas σ norm was around 1.25 for the NCEP input (Fig. 5).
A comparison of the model results to the satellite data demonstrated slightly worse but still qualitatively acceptable agreement (Table 1; Fig. 3).Among the obvious reasons for the disagreement is the low and irregular temporal resolution of the remote-sensing data as well as the differences between the surface temperatures in situ and the values derived empirically from the remote infrared radiation measurements.The validation of the satellite data against the bulk in situ temperatures yielded similar divergences to the model-satellite comparison, indicating that the uncertainties in the satellite data were the primary reason for their deviation from the modeling results.The satellite data from the Pathfinder SST produced generally higher surface temperatures than both the in situ data and the modeling results, with the bias varying between 0.4 K for the 3 years of observations and 0.7-1.5 K for the 30 years of the modeling results (Table 1).lines are the model predictions; the thin solid (blue) lines are the in situ bulk temperatures at 0.5 m depth; the half-tone (gray) lines are the skin surface temperatures determined from the upward longwave radiation (Li et al., 2016); and the symbols are the available T s estimations from satellite data.
The strongest deviations of the model results from the in situ temperatures and remote-sensing data were observed in the late spring and early summer (Figs. 2, 3).These deviations apparently resulted from the model uncertainties in the prediction of the ice cover breakup.No long-term time series on the ice breakup dates are available for the lakes of the region, making it impossible to quantify the model performance on the ice regime.An indirect estimation of the ice breakup dates from the surface temperatures from 2011-2013 suggests that the model error amounts to several weeks.
Since the model driven by the ERA data performed generally better than with the NCEP input, the results on seasonal characteristics and climatic trends are discussed below based on the ERA-FLake configuration.The outcomes of the NCEP-driven model were overall similar to those of the ERA-driven model with slightly higher water temperatures and weaker thermal stratification.The statistics are presented in Table 2, and an extended discussion of the NCEP results is omitted for the sake of brevity.

Seasonal temperature and mixing regime
The modeled seasonal course of the thermal stratification and ice regime was typical for a dimictic regime (Hutchinson, 1975;Kirillin and Shatwell, 2016) in both Tibetan lakes (Table 2).The ice-covered period lasted for 26-29 weeks and the ice thickness reached its maximum of 0.8-1.0m in early March.The ice melted on 16 May ± 1 week, followed by a 2-3 week period of full mixing (spring overturn).These phenological characteristics of the seasonal stratification were similar in Ngoring and Gyaring since they were not significantly affected by the differences in the morphology and water transparency.Both lakes stratified for 16-19 weeks in summer (Fig. 6).The autumn overturn (period of complete mixing) started around 14 October in Ngoring and about 3 weeks earlier in the shallower Gyaring.The autumn mixing also lasted longer in Gyaring, so summer stratification was shorter.The full mixing period lasted 4-6 weeks, which is 2 times longer than the spring overturn; this is a remarkable feature of the seasonal mixing regime with important effects on the lake biogeochemistry as discussed below.The summer surface temperatures in Ngoring (12.01 • C) and Gyaring (12.16 • C) were similar.The temperatures at the lake bottom were, in turn, appreciably higher in the shallower Gyaring ), which also implies generally higher mean temperatures in this lake.The differences in the water transparency between Ngoring and Gyaring were evident primarily in the thickness of the surface mixed layer h mix (epilimnion) in summer: in Ngoring h mix was 1.5-2.0m deeper than in Gyaring (Fig. 7).

Trends in atmospheric forcing
The mean trends in the principal atmospheric variables differed in the NCEP and ERA reanalyses (Table 3).The ERA data showed significant trends only in air temperature (warming at ∼ 0.28 K decade −1 ) and in air humidity.
A positive humidity trend was present in the NCEP data, too, whereas no trend existed in the air temperature.In turn, the NCEP data suggested a weak but statistically significant increase in cloud cover and a barely distinguishable decrease in wind speeds significant at a ∼ 97 % confidence level.The picture was clearer when the trends were split into seasonal bins: both reanalyses showed an increase in air temperatures in summer of 0.2-0.5 K decade −1 (with higher trend values in the ERA data).T a , however, remained unchanged in spring and sank in autumn (statistically significant in NCEP only at −0.5 K decade −1 and insignificant in ERA at −0.15 K decade −1 ).The most significant increase in T a occurred in winter (slightly stronger in ERA at 0.7 K decade −1 ), resulting in a positive trend in annual T a in the ERA data, but an insignificant trend in annual T a in the NCEP data.The decrease in the solar radiation in summer, as the major driver of the regional heat budget, was remarkable in both datasets (−3.8 W m −2 decade −1 in ERA vs. −1.2W m −2 decade −1 in NCEP).

Long-term lake response to atmospheric trends
Several characteristics of the seasonal thermal regime did not reveal any significant long-term trends in either of the two lakes based on model simulations with both reanalysis datasets.These characteristics included the dates of onset and breakup of the ice cover, ice thickness, dates of full mixing in spring, dates of the maximum ice thickness and the maximum surface temperature, bottom temperature, and mean lake temperatures in summer.The fact that the mean water temperature remained essentially unaffected by longterm atmospheric trends in both lakes is rather unexpected in view of the numerous reports on lake warming throughout the world.In the ERA-driven results, the surface temperature, T s , did not increase significantly in either lake, whereas according to the NCEP reanalysis, T s increased at 0.24 and 0.18 K decade −1 in Ngoring and Gyaring, respectively (Table 4).

Model performance
FLake demonstrated good abilities to capture the thermal characteristics of high-altitude lakes.In particular, the good simulation of the surface temperature indicates that the model is suitable for lake representation in coupled landatmosphere modeling of the Tibetan Plateau and other alpine regions.The performance of the model for the two Tibetan lakes is comparable to that reported previously in temperate (Kirillin, 2010;Stepanenko et al., 2014;Shatwell et al., 2016) and tropical (Thiery et al., 2014) lakes.The results underscore the importance of the atmospheric input for the cor-   rect simulation of lake thermodynamics.The ERA-Interim, a more recent reanalysis dataset with updated data assimilation, provides an appreciably better FLake output, while the older generation of the NCEP/NCAR dataset demonstrates an appreciable bias in lake temperatures as well as a slightly higher root mean square error.A better performance of the ERA reanalyses over the Tibetan highlands, particularly in reproducing the observed long-term trends in air temperature, has been also reported in earlier studies (Bao and Zhang, 2013).The major drawback of the model, which should be taken into account when simulating the thermal regime of Tibetan lakes as well as other high-altitude lake systems with an extended ice-covered season and low snow precipitation in winter, is an oversimplified representation of the winter conditions during the ice-covered period.The performance of FLake at simulating ice cover in temperate freshwater lakes was previously analyzed by Bernhardt et al. (2012).FLake, as well as all other lake parameterization schemes currently used in land-atmosphere coupling of climate models (Hostetler et al., 1993;Subin et al., 2012), neglects the heating of the water column by solar radiation penetrating the ice cover.One of the negative effects of this simplification is the error in simulating the ice thickness and the ice breakup date.It was not possible to quantitatively estimate this effect in the present study due to the lack of long-term observations on ice-on and ice-off dates.Nevertheless, the modeled duration of the ice cover from October-November to May-June and the ∼ 1 m maximal ice thickness in Tibetan lakes agrees with the existing occasional field observations (Wen et al., 2016).The modeling results on the ice thickness may be assumed as more reliable than the ice breakup dates since ice growth is mainly governed by the ice-atmosphere heat budget (Aslamov et al., 2014;Leppäranta, 2015), which is well captured by the model.Another apparent error introduced by neglecting the under-ice radiative heating is that the mean heat content of the lake after ice-off is underestimated.The result is a strong underestimation of the lake surface temperatures in early summer.Taking into account the strong constraints on the computational costs of the lake parameterization schemes, a remarkable improvement in regional climate models can be achieved by extending the lake modules to include one of the bulk algorithms developed earlier for convective mixing based on integration across the mixed layer of the equation for mixing energy production with solar radiation as a volumetric energy source (Mironov et al., 2002;Oveisy and Boegman, 2014).
The surface temperatures provided by the AVHRR Pathfinder were generally higher than both the in situ and the modeled temperatures.The reason for the positive bias is not clear; we may tentatively relate the error to an overestimation of the cool skin effect by the current AVHRR algorithm and hypothesize that corrections for the atmospheric radiation, used in the empirical Pathfinder SST algorithms, do not perform well over high-mountain regions with low atmospheric density and radiation.The cool skin effect, while not considered directly in our study, is expected to be even stronger in high-altitude lakes than in the ocean due to strong radiative heating and permanent cooling at the lake (Li et al., 2015).

Seasonal thermal regime and ice conditions
Both studied lakes are clearly dimictic and share the typical features of this type of mixing regime found in other regions: two periods of weak vertical mixing lasting several months separated by shorter periods of full vertical mixing (overturns) in spring and in autumn.As follows from the modeling results, a particular difference of the dimictic regime in Tibetan lakes from those in lowland temperate regions is the comparably short period of spring overturn (SO), so that summer stratification forms only about 1 week after ice-off.This is apparently due to the alpine climatic conditions under which the low air temperatures produce relatively long icecovered periods that last until late May or early June.The solar radiation flux in the early summer is, however, much stronger than at lower altitudes and quickly heats the upper water column directly after ice-off.
Since spring overturn was short, the major biogeochemical processes that take place when oxygenated surface water comes into contact with the highly mineralized sediment should occur in autumn, when overturn is longer.The characteristics of the biogeochemical processes in Tibetan lakes are largely unknown, but some key features, particularly the oxygen regime, may be deduced from the seasonal stratification pattern (Golosov et al., 2007).Due to the low eutrophication of both lakes, the oxygen consumption rates are probably low.Strong vertical oxygen gradients and an oxygen deficit in the deep parts of the lakes probably develop during icecovered periods, when the lake surface is isolated from the atmosphere for several months.This oxygen deficit, in turn, may trigger biogeochemical processes, such as methane production, at the interface between the water and the highly mineralized lake sediment.The spring overturn may be too short to completely replenish the depleted oxygen after winter, so the anoxic conditions may persist in the deep layers of lakes throughout the summer season, isolated by the stratification.In this case, the most intensive exchange between the lake hypolimnion and the atmosphere, including potential emissions of methane and carbon dioxide, should most probably appear during the autumn overturn.

Lake response to climatic trends in external forcing
The surface temperatures of lakes are recognized as important indicators of regional climate change driven by global warming (Adrian, 2009).In temperate climate regions, warming trends in lakes are reported to be close to the air temperature trends (Arhonditsis et al., 2004;Danis et al., 2004;Kirillin, 2010;Schneider and Hook, 2010).Some studies even suggest that temperate lakes warm faster than the atmosphere due to the shortening of the ice-covered period and the resulting higher storage of incoming solar radiation (Austin and Colman, 2007;Kintisch, 2015).On the other hand, tropical and alpine lakes tend to warm more slowly than the air above (Livingstone, 2003;Vollmer et al., 2005;Coats et al., 2006).A recent synthesis of available lake temperature data (O'Reilly et al., 2015) has revealed high geographic variability in lake warming trends over the globe.In particular, no significant warming trends were determined for lakes on the Tibetan Plateau.Our results indeed suggest that Tibetan lakes have not been warming during the last several decades and provide insight into the mechanisms behind this.
Atmospheric warming trends were present in both reanalysis datasets, although these were unevenly distributed over the seasons.Air temperature increased in summer (more strongly in the ERA reanalysis at 0.5 K decade −1 than in NCEP at 0.2 K decade −1 ), remained unchanged in spring, decreased in autumn (significantly only in NCEP at −0.5 K decade −1 and insignificantly in ERA at −0.15 K decade −1 ), and increased in winter (slightly more in ERA at +0.7 K decade −1 vs. +0.4K decade −1 in NCEP).The winter warming trend was stronger than the summer trend, but its effect on the lakes is weak due to isolation from the atmosphere by the ice cover.The statistically significant but rather small tendency for ice thickness to decrease at ∼ 1 cm decade −1 simulated in both Ngoring and Gyaring may be related to the overall effect of a warmer atmosphere in winter.
Despite the overall warming of the atmosphere, which was especially apparent in the ERA dataset, there was no significant trend in the lake surface temperatures in the ERAforced results and only a slight increasing trend in the NCEP results.This fact is closely connected to the negative trend in the solar radiation, which decreased in the ERA dataset at 3.8 Wm −2 decade −1 , more than 3 times faster than in the NCEP dataset.A characteristic feature of lakes as compared with lowland lakes is their stronger heating by solar radiation and, as a result, a lake surface that is persistently warmer than the atmosphere (Wen et al., 2016).Hence, the effects of the air temperature increase on the lake heat budget were counteracted by the reduced radiative heating, resulting in no significant long-term change in the mean lake temperatures (Fig. 8).This core result of the present study indicates that the trends in the lake temperatures are not coupled directly to the increase in the air temperature, which underscores the crucial role of shortwave radiation for the thermal conditions in alpine lakes.
In general, the response of the stratification regime to changes in external atmospheric forcing in the deep and clear Ngoring was stronger than in the shallow and more turbid Gyaring, which was most apparent in the NCEP-forced results.A plausible explanation for this effect is that the high turbidity of the Gyaring water damped the effects of the higher air temperatures and weaker solar radiation.In the transparent waters of Ngoring, the effect of the radiation decrease is distributed over a water column several meters thick, whereas the increase in the air temperature affects the heat transport at the upper air-water interface.This decoupling of the two counteracting trends produced a temperature increase at the upper boundary with a simultaneous decrease in the heat content in the water column beneath, resulting in a shallower mixed layer and stronger stratification, as found in the modeling results for Ngoring, but not for Gyaring (see the Results section).In Gyaring, on the other hand, both the decrease in radiation and the increase in air temperature are concentrated at the lake surface; the daytime heat is stored in a shallower surface layer and is more quickly released to the atmosphere by nighttime cooling.Thus, the opposing trends in air temperature and radiation effectively cancel each other at the lake surface with little effect on the thermal characteristics of the deeper waters.
To quantify the counteracting effects of the increasing air temperature and the decreasing shortwave radiation, we performed model runs using artificial shortwave radiation input constructed from ERA values (which had a significant negative trend) with the summer radiation trend removed.The trend was removed based on the daily averaged values, and the diurnal cycle of solar radiation was subsequently restored to the detrended values based on daytime duration and the daily maxima of solar radiation.Additionally, a "control" model run was performed with the radiation input reconstructed from the daily averages using the same procedure, but without trend removal.The model results based on the modified "control" radiation input did not differ significantly from those based on the original ERA input (not shown).
In turn, the model driven by the detrended radiation input produced a positive trend in the mean summer surface temperature of 0.26 K decade −1 and a maximum summer surface temperature of 0.34 K decade −1 .The trend in surface water temperatures was significant at a > 99 % level of confidence (p = 0.003), whereas the original ERA input produced no significant long-term increase in surface temperatures (0.15 K decade −1 at a < 93 % level of confidence; see Table 4).Hence, the shortwave radiation decrease did indeed cancel the effect of the air temperature increase on lake temperatures in ERA-driven scenarios.The hypothetical surface temperature trend of 0.26 K decade −1 is still lower than the summer air temperature increase of 0.50 K decade −1 and is accompanied by significant changes in the summer stratification: a later breakdown of the stratification due to autumn cooling (1.8 days decade −1 ; p = 0.034) and a shallower surface mixed layer depth (0.14 m decade −1 ; p = 0.006).Several previous studies devoted to alpine or tropical lakes also reported lake temperature trends lower than the regional atmosphere warming (Livingstone, 2003;Vollmer et al., 2005;Coats et al., 2006), while lowland temperate lakes are often found to be in thermal equilibrium with the atmosphere (Arhonditsis et al., 2004;Danis et al., 2004;Kirillin, 2010).
The radiation decrease can be hypothetically related to the direct effect of anthropogenic aerosols at high altitudes.An-other possible reason for the lower solar radiation flux in the Tibetan highlands might be increased evaporation and the resulting higher air humidity, which is also supported by positive trends in the humidity and cloud amounts present in the reanalysis data.In particular, Shen et al. (2015) report local cooling over the Tibetan Plateau due to increased evaporation resulting from increased vegetation growth.The exact mechanisms of the decrease in the solar radiation of the Tibetan Plateau require further investigation.Our model experiments are well supported by the observational data from the largest freshwater lakes of the Tibetan Plateau.They compare the dynamics of lakes with different depths and water transparencies under the same atmospheric forcing that is characteristic for high-altitude conditions.The presented results provide the first detailed insight into the specific mechanisms of internal lake mixing and lake-atmosphere coupling; these are specific to the Tibetan lakes and extendable to other highaltitude hydrological systems.It should be noted that freshwater lakes comprise only about 10 % of the Tibetan lake system (Jiang and Huang, 2004), while ∼ 75 % of the lakes have brackish or saline water (the salt content of the rest of the lakes is unknown).The salt stratification can significantly contribute to the seasonal mixing regime as well as to the lake response to long-term changes in the atmospheric input.The effects of salinity as well as the lake physics during the icecovered period, as discussed above, are the two major facets of lake dynamics for further clarification with regard to the specific features of the Tibetan lake system.
Data availability.The results of modeling and of statistical analysis can be downloaded from http://www.flake.igb-berlin.de/tibet/tibet_lakes_model.zip(Kirillin et al., 2017) and http://www.flake.igb-berlin.de/tibet/tibet_lakes_statistics.zip(Shatwell et al., 2017) respectively.The temperature measurements from Ngoring Lake are part of of the ongoing monitoring program and are available from Lijuan Wen (wlj@lzb.ac.cn) by request.NCEP reanalysis data were provided by the NOAA/OAR/ESRL PSD in Boulder, Colorado, USA from the website at https://www.esrl.noaa.gov/psd/.ERA-Interim data were provided courtesy of the ECMWF.AVHRR data were provided by GHRSST and the US National Oceanographic Data Center.
Author contributions.G. Kirillin conceived and designed the study and performed the modeling and analysis of the remote-sensing data.L. Wen performed the collection and analysis of the field measurements.T. Shatwell evaluated the model performance and performed the statistical trend analysis.All authors contributed to the final analysis of the results.G. Kirillin and T. Shatwell wrote the paper with contributions from L. Wen.
Competing interests.The authors declare that they have no conflict of interest.

Figure 1 .
Figure 1.The study site map and the points of the AVHRR Pathfiner grid used for estimation of the surface temperature from satellite data.The green shaded region is the Tibetan Plateau, shown as land above 2000 m a.s.l.
physics.The parameters used to force the lake model were 6 h surface layer characteristics: 2 m air temperature T a [recalculated to • C], 2 m air humidity e a [mb], incoming solar radiation I R [W m −2 ], 10 m wind speed U 10 [m s −1 ], and total cloud cover N [0-1].The period 1975-2014 was applied for the modeling forced by NCEP/NCAR, and a shorter period of 1979-2014 was used for the ERA-driven modeling because the ERA-Interim data have been available since 1979.

Figure 3 .
Figure 3.The modeled lake surface temperatures (lines) and the Pathfinder SST data (symbols) for 1982-2012 in (a) Ngoring and (b) Gyaring.The thick lines are from the FLake model forced by the ERA-Interim reanalysis, and the thin gray lines are from the NCEP/NCAR forcing.

Figure 4 .
Figure 4.A comparison of the results from the FLake model forced by (a) NCEP/NCAR and (b) ERA-Interim reanalyses against the in situ lake surface temperatures in Ngoring.The thick solid (red)lines are the model predictions; the thin solid (blue) lines are the in situ bulk temperatures at 0.5 m depth; the half-tone (gray) lines are the skin surface temperatures determined from the upward longwave radiation(Li et al., 2016); and the symbols are the available T s estimations from satellite data.

Figure 5 .
Figure 5.A comparison of the observed and modeled lake surface temperatures in Ngoring (a, c) and Gyaring (b, d) forced by the ERA-Interim reanalysis (a, b) and the NCEP/NCAR reanalysis (c, d).The blue circles are derived from satellite data and the green squares from in situ measurements.Only the values from May to October are shown.

Figure 6 .
Figure 6.The modeled seasonal thermal stratification pattern and ice cover development in (a) Ngoring and (b) Gyaring.The model was forced by the ERA-Interim reanalysis.Only the values from May to December 1979 are shown.The black line is the ice cover thickness.The white line is the bottom of the surface mixed layer defined as the depth of the maximum vertical temperature gradient.

Figure 7 .
Figure 7.The typical modeled vertical profiles of temperature in Ngoring (stars) and Gyaring (circles) in the middle of the summer stratification period (12 August 1979).

Figure 8 .
Figure 8.The trends in the annual mean air temperature (a) and solar radiation (b) from ERA-Interim, as well as the trends in the modeled mean summer lake temperature (c) in Ngoring (blue) and Gyaring (green).
performance at short and long time scales, respectively.The model is forced by reanalysis data in two different configurations corresponding to the two neighboring large freshwater lakes located in the northeastern part of the plateau, Ngoring and Gyaring.Both lakes undergo virtually the same atmospheric forcing, but differ in mean depth and water transparency.

Table 1 .
The goodness of fit statistics for the simulations in Ngoring and Gyaring based on surface temperatures.The simulations were forced by either the ERA-Interim or the NCEP/NCAR reanalysis.The last row shows a comparison of the satellite-derived surface temperatures with in situ measurements in Ngoring.The statistics are described in the Methods section (Eqs.3-7).

Table 2 .
The mean characteristics of the seasonal thermal regime in Ngoring and Gyaring forced by the NCEP/NCAR and ERA-Interim reanalysis datasets.T s , T b , and T m are temperatures at the surface, bottom, and averaged over the whole water column, respectively; h mix and h ice are the thickness of the mixed layer and the ice cover respectively; and doy is day of the year.

Table 4 .
The trends in the characteristics of the lake thermal regime; abbreviations as for Table2.Marginally significant trends (0.05 < p < 0.1) are marked with italics.