Research ArticleCLIMATOLOGY

Quantifying the evaporation amounts of 75 high-elevation large dimictic lakes on the Tibetan Plateau

See allHide authors and affiliations

Science Advances  26 Jun 2020:
Vol. 6, no. 26, eaay8558
DOI: 10.1126/sciadv.aay8558


Lake evaporation can influence basin-wide hydrological cycles and is an important factor in loss of water resources in endorheic lakes of the Tibetan Plateau. Because of the scarcity of data, published lake evaporation values are inconsistent, and their spatial distribution has never been reported. Presenting a plausible hypothesis of energy balance during the ice-free seasons, we explored the multiyear (2003–2016) average ice phenology and evaporation amounts of 75 large dimictic lakes by using a combination of meteorological and satellite data. Evaporation amounts show large variability in spatial distribution, with a pattern of higher values in the south. Lakes with higher elevation, smaller area, and higher latitude are generally associated with a shorter ice-free season and lower evaporation. The total evaporated water amounts have values of approximately 29.4 ± 1.2 km3 year−1 for the 75 studied lakes and 51.7 ± 2.1 km3 year−1 for all plateau lakes included.


Featuring the largest high-elevation inland lake zone and the highest volumes of glaciers outside the Antarctic and Arctic regions and serving as the birthplace of many major rivers, the Tibetan Plateau (TP) is considered the “Water Tower of Asia”; however, it is now in danger due to rapid climate warming (1). Influenced by the elevation-dependent intensive warming in this area, water resources over the TP are showing marked changes, including an increase in precipitation; spatial inconsistency in snowmelt and snow cover (SC) change; thickening of the active soil layer depth; degradation of permafrost; notable reduction of glacier lengths, areas, and masses, resulting in glacier retreat; and many other hydrological and meteorological variations (2, 3). These intensive hydrological changes can enhance meltwater from the cryosphere, thereby affecting runoff and influencing the level, area, and water storage of lakes (4, 5). TP lakes, which number to more than 32,000 and account for approximately 50% of the lake area in China, have shown marked increases in number and area in recent decades (5). By a combination of in situ measurements, satellite data, model simulations, and hydrological empirical relationships, many previous studies have focused on water balance analyses of these high-elevation lakes (6, 7), for which lake evaporation amount is the most important loss factor. Thus, quantifying the evaporation amounts of these high-elevation lakes is important for water resource management and has high priority in the basin-wide hydrological cycle and water budget analyses.

Lake evaporation amounts of high-elevation lakes have been estimated by various methods; however, the reported values have large discrepancy in their seasonal variations and annual amounts (8, 9). Lake evaporation amounts estimated by different methods can lead to contradictory water budget conclusions for a lake [i.e., water seepage (6) and much higher lake evaporation and water volume change values relative to the water supply (7) in Lake Nam Co]. For ice-covered seasons, ice sublimation was simply assumed to be negligible in previous studies (4, 10); however, in situ eddy covariance (EC) measurements (9) and model simulations (11) contradicted this assumption. For ice-free seasons, remarkable differences in the lake-air interaction processes have been reported through in situ EC measurements over a small lake and a large lake, which have quite different ice phenology, meteorological and environmental conditions, and seasonal variation of evaporation amounts (9). Thus, lakes with different areas, elevations, and latitudes may show quite different patterns in their ice phenology and evaporation amounts. However, most previous publications on evaporation amount estimations for the TP have focused on several specific lakes, and the evaporation amounts over numerous other lakes and their spatial distributions have not been reported.

Recently, EC measurements have been used to quantify the evaporation amounts of several high-elevation lakes, including Lake Nam Co (9), Lake Qinghai (12), Lake Serling Co (13), and Lake Ngoring (14). However, energy budget imbalances have been obtained by EC-based in situ measurements over waters and lands worldwide. For example, energy budget closure (EBC) ratios show values of 0.859 from July to November in Lake Nam Co (9) and 0.836 from 10 April to 30 September in a small lake adjacent to Lake Nam Co (15), whereas the EBC value of a small boreal lake in Finland with all energy budget components estimated by in situ measurements only has a value of approximately 0.77 for two natural years (16). Moreover, other studies that have performed in situ EC measurements on energy budget analyses over lakes also show obvious imbalances. The energy budget imbalance results over lake surfaces are similar to those estimated at land stations (17), which are likely due to the undetected large eddies and secondary circulation by EC measurements (18). Compared with land stations, the EBC values over lake surfaces can be significantly influenced by the large amount of heat storage (release) during spring and summer (autumn and winter) and show large seasonal variations. The EBC values are lower than the unit value during heat-storage seasons, larger than the unit value during heat-release seasons (14, 15), and much closer to the unit value during the complete ice-free seasons (15). Thus, lake evaporation amounts directly from EC measurements can be more convincing with a cross-validation via energy budget methods, where the Bowen ratio (Bo) derived from EC measurements has been found to be close to that estimated from meteorological variables.

The high-elevation lakes of the TP experience high solar radiation and low air temperature; thus, most of these lakes are dimictic and have an obvious ice-covered season in winter and two instances of fully water-column mixing with a constant temperature profile of approximately 4°C, with one instance occurring shortly after ice melt and the other occurring shortly before ice formation. Thus, a plausible hypothesis is that all the energy stored after the ice melt stage can be completely released before the ice-forming stage, and the evaporation amounts during ice-free seasons of these dimictic lakes can be simply estimated by the net radiation and Bo. In this study, Bo is used to allocate the ratio of heat loss by conduction to that by evaporation and further validated with existing in situ EC measurements over the TP (9, 1214). The objectives of this study are to: (i) quantify evaporation amounts of 75 high-elevation large dimictic lakes on the TP and (ii) explore the spatial distribution of the evaporation amounts and the influencing factors. The 75 large dimictic lakes show specific characteristics in their environmental and meteorological conditions and can be classified into three groups (see details in table S1): (i) the “Southern Lakes Group” (SLG; 34 in total) is characterized by the longest duration of ice-free season (DIF) and located in the south (mainly alongside 31°N); (ii) the “Northern Lakes Group” (NLG; 32 in total) is characterized by high-elevation and low lake surface temperature (Ts) and air temperature (Ta) and located in the north (mainly alongside 35°N); and (iii) the “Other Lakes Group” (OLG; 9 in total) is characterized by low elevation and high Ts and Ta. Our study provides an estimation of the evaporation amounts of dimictic lakes on the TP and may provide deep insights in revealing climate change impacts and reducing uncertainties for research on catchment-scale water budgets, hydrological cycles, and water resources management.


Correction and sensitivity analysis of the forcing data

The meteorological data from the China Meteorological Forcing Dataset (CMFD) and satellite data are vital variables for estimating Bo and evaporation amounts. Thus, the correction equations for these variables are established using in situ measurements from 14 Tibetan Observation and Research Platform [TORP (19)] stations. Air pressure shows an underestimation in 8 of 11 stations, with an average bias value of −5.6 hPa. Ta shows high consistency, with all correlation coefficients higher than 0.97 and an average mean bias (MB) value of 0.71°C. The slopes of the fitting lines higher than the unit value in 8 of 12 stations indicate an underestimation in the CMFD. Thus, a linear equation (y = 1.028 × x − 0.68; fig. S1A) can improve Ta. Similarly, specific humidity (Qa) shows good correspondence and presents correlation coefficient, mean absolute error (MAE), and root mean square error (RMSE) values of 0.976, 0.00059, and 0.00082, respectively. The higher than unit values of the slopes of fitting lines indicate an underestimation of the CMFD, and a linear fitting equation (y = 1.105 × x − 0.00001; fig. S1B) can improve the performance of Qa. Downward longwave radiation (Rl↓) has a high correlation coefficient of 0.927, and a linear fitting equation of y = 1.22 × x − 338.3 (fig. S1C) can decrease the MB, RMSE, and MAE values for the majority of the stations. For downward shortwave radiation (Rs↓), the slope of the linear fitting line passing through the origin has a value of 0.986 (fig. S1D), and it is quite close to the 1:1 line. Furthermore, 9 of 13 stations have an average relative error of less than 5%. Thus, the original Rs↓ is used without correction.

Sensitivity analysis of these meteorological variables to the evaporation amounts indicates that biases of ±20 hPa in air pressure and ±10% in Qa can only lead to approximately ±10 mm variations (approximately 1 to 2%) in open-water evaporation amounts (Eo); biases of ±1.5 ° C in Ta or Ts correspond to approximately ±35 mm (approximately 4 to 5%) of Eo; and variations in Rs↓ and Rl↓ of ±5% have the most remarkable impacts on Eo and can result in variations of approximately ±65 mm (8 to 9%) in Eo. Thus, the CMFD after correction by the aforementioned equations has an acceptable residual bias to the estimated Eo, and more information about the evaluation metrics can be found in tables S2 to S6.

Analysis of ice phenology and related meteorological variables

Table S7 lists the multiyear (2003–2016) average ice phenology of freeze onset (FO), freeze up (FU), break up (BU), and water clean of ice (WCI) events in the 75 studied lakes, and examples of the identified ice phenology events in four lakes with different areas are shown in the top panel of Fig. 1. The seasonal variations of ice phenology are quite clear. For example, in Lake Qinghai (Fig. 1B1), the percentage value of ice (or water) coverage can clearly distinguish the dates of FO [day of year (DOY) = 337], FU (DOY = 361), BU (DOY = 75), and WCI (DOY = 101), and similar results were observed in Lake Ngangla Ring Co (Fig. 1B2), Lake Aksayqin (Fig. 1B3), and Lake Xiangyang (Fig. 1B4). The percentage values of ice and water are complementary, and the results from Terra MODIS (Moderate Resolution Imaging Spectroradiometer) and Aqua MODIS are both similar. The derived multiyear average ice phenology events are similar to those derived by Kropáček et al. (20), where various satellite observations have been used for validation and acceptable accuracy has been obtained. Lake Yamzhog Yum Co, with an area of 638 km2, is an exception due to its complex shoreline, which can result in a limited footprint and lead to large amounts of land-contaminated measurements, not only in ice phenology but also in lake surface temperature.

Fig. 1 The ice phenology events identification and Bo evaluation of eight studied lakes.

(A) The locations of the 75 studied lakes and the 14 TORP stations on the TP. The small and red numbers represent lakes’ serial numbers, by which the attributes of each lake can be found in Table S7. The “black dots” indicate positions of the 14 TORP stations, with names listed close to each dot. The studied lakes are classified into three groups: the “blue box” indexed with blue number “1” indicates the NLG, the red box indexed with large red number “2” indicates the SLG, and the other lakes outside indexed with black number “3” indicates the OLG. The identified multiyear average ice phenology events are marked for Lake Qinghai (B1), Lake Ngangla Ring Co (B2), Lake Aksayqin (B3), and Lake Xiangyang (B4), where the seasonal variations of percentage values of ice and water are plotted with uncertainties marked. In addition, the dates of BU, WCI, FO, and FU are marked as symbols of “is-greater-than,” “circle,” “star,” and “is-smaller-than”, respectively. The seasonal variations of in situ measured and satellite-derived Bowen ratios are shown for Lake Nam Co (C1), Lake Serling Co (C2), Lake Qinghai (C3), and Lake Ngoring (C4). The “black square line” indicates the in situ EC measurements during their observational periods, the “red star line” indicates satellite-derived values during the year of in situ measurements, and the “green circle line” indicates the satellite-derived multiyear (2003–2016) average values.

The three lakes with the longest DIFs are Lake Xuru Co, Lake Paiku Co, and Lake Tangra Yum Co, and they have DIF values of 285, 265, and 264 days, respectively. All three lakes are located in the SLG and do not freeze up in certain years (20). In contrast, the three lakes with the shortest DIFs are Lake Hoh Xil (111 days), Lake Yinma (112 days), and Lake Xiangyang (117 days), which are close to each other and situated in the NLG. The dates of freeze and melt events correlate most strongly with air temperature (21); thus, two temperature-related variables (elevation and latitude) are correlated with the ice phenology events of FO and WCI. Elevation and latitude have negative correlation values (−0.319** and −0.721**, “**” indicates correlation coefficients passing the significance level of 0.01) with FO event and positive correlation values (0.669** and 0.468**) with WCI event; thus, it indicates that a lake with a lower elevation and lower latitude will probably correspond to a late date of ice formation event and an early date of fully ice-free event.

Furthermore, to test the reliability of these relevant variables, the relationships of LST3 (lake surface temperature at the warmest 3 months), meteorological variables (Rs↓, Rl↓, and Ta), and lakes’ attributes (DIF, elevation, area, and latitude) during ice-free seasons are explored, and the correlation coefficients of these variables are shown in the lower left of Table 1. LST3 has a significant negative relationship with elevation (−0.787**). For example, three of five lakes with highest LST3 are Lake Tuosu, Lake Sugan, and Lake Qinghai, which have low-elevation values of 2793, 2805, and 3194 m, respectively. The other two lakes (Lake Zabuye and Lake Laguo) with high LST3 mostly result from a high salinity, a small area, and a shallow depth. Significant positive correlations (0.891** and 0.931**) of Rl↓ and Ta with LST3 result from the intensive lake-air interaction, where a higher LST3 also corresponds to a longer DIF (0.695**) occurring in the SLG. Solar radiation is generally higher in lakes at a lower latitude; thus, a negative correlation (−0.28*) exists between latitude and LST3. To remove the influence of elevation and latitude and explore the relationship between LST3 and area, lakes with similar elevations (4500 to 4800 m) and latitudes (30°N to 32°N) are analyzed, and a significant negative correlation coefficient value of −0.574** is obtained. Moreover, the DIF has significant positive correlations with Rl↓, Ta, LST3, and area, and significant negative correlations with elevation and latitude. Thus, a lake with a large area, a high surrounding atmospheric heating, a low elevation, and a small latitude over the TP tends to have a much longer DIF. All these relationships are reasonable and demonstrate the reliability of our results.

Table 1 The correlation coefficients of the variables.

The lower left is for elevation (Ele), area, latitude (Lat), Rs↓, Rl↓, Ta, DIF, and LST3; the upper right is for Bo, Rn, Ele, area, Lat, DIF, and Eo. “*” and “**” indicate the correlation coefficients that pass the significance level of 5 and 1%, respectively.

View this table:

Validation of the Bo and evaporation amounts

The estimated Bo and evaporation amounts are validated by the published EC observations of the four high-elevation lakes (Lake Nam Co, Lake Qinghai, Lake Serling Co, and Lake Ngoring). The increasing trends of Bo during the ice-free season of a large lake can be generally reproduced by our methods, with large errors existing mainly in Lake Serling Co (Fig. 1C2). The average Bo during the observational periods of EC measurements is 0.189 in Lake Nam Co, 0.229 in Lake Qinghai, and 0.360 in Lake Ngoring, while the corresponding satellite-based values are 0.179, 0.211, and 0.384, respectively. Thus, the close values between in situ measurements and the estimation indicate a good performance of the method used in this study. The average seasonal variation of Bo from the 75 studied lakes shows a weakly decreasing trend shortly after ice melt (DOY = 90 to 190) and a continuous increasing trend afterward, and these variations correspond to the seasonal variation of water-air temperature gradients, with higher values in cold months and lower values in warm months. This pattern can also be evidenced by the EC observations in Lake Nam Co and Lake Ngoring.

To obtain the annual evaporation amounts (Ea), winter sublimation during the ice-covered season cannot be ignored. The EC-based LE (unit, W m−2; L = 2.45 MJ kg−1 is the latent heat of vaporization; E is the evaporation amount; unit, mm; 1 W m−2 is equal to approximately 1.059 mm/month) has monthly values of 59.8, 43.1, and 48.6 W m−2 for January, February, and March, respectively, in a small lake adjacent to Lake Nam Co in 2012 and values of 62.1, 34.2, and 34.3 and 55.7, 27.8, and 34.3 W m−2 for February, March, and April in Lake Nam Co in 2016 and 2017, respectively. Furthermore, winter ice sublimation over a small high-elevation lake by a high-resolution thermodynamic ice model reported an average value of 45 W m−2 from November to April, with the smallest monthly value of 31 W m−2 in December and largest monthly value of 54 W m−2 in February (11). Thus, relative to evaporation amounts during ice-free months (9), winter sublimation during the ice-covered season is relatively small. The values from in situ measurements and model simulations are consistent, with a monthly value of approximately 45 W m−2 (≈47.7 mm) by averaging all the aforementioned values. Thus, Ea is summed up by Eo and winter sublimation, with the latter determined by the duration of ice-covered period.

The EC-based evaporation amount from May to January (Eo) in Lake Nam Co sums up to a value of approximately 915 mm (fig. S2F), which is 6.9% lower than the estimated value of 981 mm (9) using in situ Rn and Bo. The two values are both within the uncertainty of satellite-derived Eo value of 916.4 ± 89.1 mm. Furthermore, the EC-based Ea (1025 mm) in Lake Nam Co is also close to the satellite-derived annual value of 1068.5 ± 63.9 mm. For Lake Serling Co, the multiyear average Ea by model simulation is estimated to be approximately 1066.9 mm, which is 17.6% lower than our estimation (1294.7 ± 59.7 mm). The differences in Lake Serling Co may be related to two factors: (i) winter sublimation is considered in our estimation, while it is ignored by the model simulation; and (ii) obvious positive values of heat storage in the water by the model simulation exist (8), and this energy seems inappropriate and can be allocated to the lake-air turbulent heat flux exchange. The EC-based Ea in Lake Qinghai is approximately 830 mm, which is 20.4% lower than the estimated evaporation value of 1042.1 ±56.4 mm in this study. Furthermore, the EC-based evaporation amount in Lake Ngoring from June to November is estimated to be 436 mm, and it is smaller than the value of 535 mm estimated by a combination of Rn and Bo. Moreover, the average Eo (ranging from 350 to 1200 mm) and Ea (ranging from 650 to 1300 mm) are approximately 865 and 1000 mm, respectively, while evaporation amounts from 1200 to 1450 mm are obtained by the Penman method (4), which does not consider winter sublimation and shows a remarkable overestimation (22). Thus, considering the widely observed energy imbalance issues of EC measurements over all types of different lakes, the nonideal EC observation environments and observation uncertainties in these high-elevation lakes, the estimated Bo, and evaporation amounts by our proposed method show reasonable and good accuracy.

Spatial distributions of the evaporation amounts and related variables

Figure 2 displays the spatial distribution of the evaporation amounts, lake attributes, and meteorological variables. Lakes with an area of larger than 450 km2 are mostly distributed in the SLG and OLG compared with the NLG (Fig. 2A), and the average area (207.9 km2) in the NLG is only approximately half of that (362.3 km2) in the SLG. Compared with the NLG, lakes in the SLG and OLG also show higher values of LST3 (Fig. 2B), which are caused by the characteristics of high thermal forcing and low elevation, respectively. Furthermore, the DIFs in these high-elevation lakes show similar spatial distributions as LST3 (Fig. 2C), with average values of 232.6, 199.8, and 152.4 days for the SLG, OLG, and NLG respectively. Loss of lake ice can contribute to the decreased availability of freshwater because of the increased evaporation rates (23). Thus, Eo has a similar spatial variation as the DIF (Fig. 2F). In addition, the spatial distribution of Ea shows a similar spatial pattern as Eo. The average amounts of Eo in the SLG, NLG, and OLG are 1028.3, 873.9, and 735.2 mm, respectively, while the average amounts of Ea are 1171.9, 1055.7, and 961.5 mm, respectively. Sublimation during the winter season approximately accounts for 12.3, 17.2, and 23.5% of the Ea in the SLG, NLG, and OLG, with these percentage values determined mainly by the lengths of ice-covered seasons. The total water resources evaporated from the lake surface of the 75 lakes are summed to be approximately 29.4 ± 1.2 km3 year−1, where the annual water resource evaporated over a lake is equal to its annual evaporation amount multiply by its area; and considering the area ratio of the 75 lakes to all the lakes, the total water resource is approximately 51.7 ± 2.1 km3 year−1. The spatial distribution of Eo is also positively correlated with the spatial distribution of Rn (Fig. 2E), which is correlated closely with the spatial distribution of DIF; it is negatively correlated with the spatial distribution of Bo (Fig. 2D), which is mostly determined by the thermal regime of these high-elevation lakes.

Fig. 2 The spatial distributions of lake evaporation amounts and related variables.

The spatial distributions of lake area (A), LST3 (B), DIF (C), Bo (D), net radiation (E), and evaporation during the ice-free season (F) of the 75 studied lakes.

To further explore the relationship between the evaporation amounts and related variables, the correlation coefficients of these variables are shown in the upper right of Table 1 and fig. S3. The evaporation amounts have a positive correlation coefficient (0.985**) with Rn and a negative correlation coefficient (−0.782**) with Bo. The negative correlation coefficients of latitude (−0.68**) and elevation (−0.53**) with DIF indicate that a lake located at a higher latitude and a higher elevation will generally have a shorter ice-free season, which will further lead to a lower evaporation amount. For meteorological variables, a larger Rl↓, which has a high correlation with Ta (0.911**), can lead to a higher Ts and then a larger Es. Thus, lakes with higher atmospheric forcing present a positive correlation with the water-air humidity gradient (δE) and negative correlation with the water-air temperature gradients (δT), with the former presenting a positive correlation (0.404**) with evaporation amounts and the latter corresponding to lower Bo and then higher evaporation amounts.


Uncertainties in this study

First, the biases of the monthly Bo are relatively large in January in Lake Nam Co (Fig. 1C1) and in December in Lake Qinghai (Fig. 1C3) when the lakes’ surfaces are in the transition season of ice formation. Ice forms much earlier in the lake shore area, where the EC measurements are located. Thus, observation uncertainties of EC measurements over ice rather than over water and the water-ice pixels by satellite observations will both contribute to the large uncertainties during the transition seasons. Second, these high-elevation lakes have quite high variability in the ice phenology events, which may influence the representative of the spatial distribution of Bo (Fig. 2D). Thus, to verify the spatial distribution of Bo during the whole open-water season, we obtained the related meteorological variables and Bo for the warmest 3 months of these lakes, including LST3, Ta3, Es3, and Ea3, and Bo3. The spatial distributions of these variables during the warmest 3 months are similar to those during the whole ice-free seasons, with correlation coefficients of the corresponding variables being 0.671**, 0.830**, 0.895**, 0.772**, and 0.889**, respectively. Thus, the spatial distribution of Bo during the whole open-water seasons of these high-elevation lakes is reliable. Third, because of the existence of lakes in the high-elevation and cold region, air temperature in a region with lakes is relatively higher, and this effect is called as warm lake effect. The warm lake effect of these high-elevation large lakes is corrected using in situ measurements at Lake Nam Co (with an area of 2021 km2); thus, the derived correction relationship may present a certain degree of bias for lakes with a smaller area. As warm lake effect is weaker in smaller lakes relative to Lake Nam Co, the lake-air temperature gradients in smaller lakes should be higher. Thus, the spatial distribution of higher Bo and lower evaporation amounts in the NLG than the SLG can be even more obvious. Fourth, the uncertainties in Eo are shown higher than those in Ea (table S2), and it is caused by the uncertainty definitions of the two variables in ice phenology (details of ice phenology definitions are shown in Materials and Methods). For Eo, the uncertainty of ice phenology is determined by the difference between two DIFs (one from WCI to FO and the other from BU to FU); and for Ea, the difference of Eo between two DIFs can be partially offset by the inverse lengths of ice-covered periods. Further, although winter sublimation only accounts for approximately 10 to 25% of the annual total water evaporated to the atmosphere, winter sublimation value of approximately 63.5 mm in February of a small lake adjacent to Lake Nam Co is approximately twice the value in March and May in Lake Nam Co. Sublimation over snow or over ice is an important influencing factor for such variability, which needs to be addressed by detailed in situ experiments, satellite measurements, and lake-atmosphere modeling in the future. Last, the total water evaporated from the lake surfaces is estimated at approximately 51.7 ± 2.1 km3 year−1. The value is roughly estimated by a combination of the evaporation amounts of the 75 large lakes and the relative proportion of their area to the total lake area of the whole TP. Because the DIFs of small lakes are usually shorter than those of large lakes, the value of 51.7 ± 2.1 km3 year−1 likely represents a slight overestimation if the impact of lake size is considered. In addition, lake depth, salinity, lake shoreline complexity, and other environmental and meteorological conditions will all influence the total evaporation estimated by our method. However, our estimation can provide roughly validated amounts of water evaporated from the lake surface, and the total evaporation value from 1990 to 2013 over all these high-elevation lakes is approximately 9.3 times of the total increase in lakes’ water storage change (4). The total amount of water resources over the TP and variations of the water resources have not been clarified, and these issues remain important topics that need to be addressed.

Factors influencing ice phenology

Ice phenology describes the processes of ice formation, growth, and decay and is the most important factor for variation in the evaporation amounts. It is determined by the energy surplus or deficit of the water surface. Thus, the factors affecting the ice phenology are mainly from the climatic and nonclimatic controls on the surface energy budget, including meteorological forcing (Ta, wind, radiation, etc.), lakes’ attributes (elevation, latitude, area, depth, salinity, shoreline complexity, lake bathymetry, etc.), and environmental conditions (snow, cloud cover, surface runoff, etc.) (20, 21, 23). Rs↓ of these high-elevation lakes are high; thus, the mean annual air temperature value of 8.4°C, which has been considered the most important factor for classifying lakes into annual or intermittent winter ice cover for low-elevation lakes (23), is not appropriate for lakes over the TP. Among all the influencing factors in Table 1, air temperature has the highest correlation coefficient with ice phenology events. However, the influence of air temperature on ice phenology is found to be highly nonlinear, and other factors will also influence the ice phenology via the surface energy budget. Significant correlations of elevation and latitude with the DIF reflect their influences on the spatial distribution of the 0°C isotherm, which has significant correlations with ice phenology events (24).

Lake attributes (area, depth, shoreline complexity, etc.) and wind speed are essential influencing factors on vertical heat transfer of lakes to cool the surface water temperatures to 0°C (25). Lakes with a larger area and depth have a much higher thermal capacity and thus require a longer period for ice formation (23). A significant positive correlation coefficient occurs between area and DIF in our results. A shorter duration of the ice cover season is expected in lakes with deeper depths (23); however, the depth information of these high-elevation lakes is still limited. Further, lakes with complex shorelines most probably have a longer ice cover season compared with those with simpler shorelines due to the influence of such complexity on fetch (23). Radiation and cloud cover can influence the heat exchange with the atmosphere and thus affect the ice phenology. A specific example in lakes at the polar region during the darkness period shows that longwave radiation is the dominant factor for ice formation and decay under such conditions (21). Furthermore, clouds are found to produce a warming effect by trapping longwave radiation and a cooling effect by reflecting solar energy away from the ice, with the dominant process being their overall effects (26). Rs↓ dominates the energy forcing of these high-elevation lakes on the TP, and they can lead to large diurnal and seasonal variations of the lake surface temperature relative to the high-latitude lakes.

SC can influence the freeze and thawing rates of ice and control the breakup dates through its insulation effect (21). For example, the winter snow accumulation over the ice surface can reflect the Rs↓, delay the ice melt process, and thus extend the duration of ice-frozen period. Another important factor for BU events in lakes is spring runoff (21). Because runoff recharges to these high-elevation lakes from glacier melt and precipitation mostly appear in the monsoon season, the dynamic effect of spring runoff can generally be ignored relative to their thermal effect. In addition to these aforementioned factors, salinity can influence the freezing temperature of a lake and is another important factor affecting the ice phenology of these endorheic lakes on the TP. For example, certain lakes do not freeze up completely in some years, and high salinity is considered an important factor (20). Currently, salinity information of these high-elevation lakes is still limited due to difficulties associated with the harsh environment and remote location. Although high variability exists in ice phenology of these lakes, the spatial distribution of climatic average ice phenology is identified clearly, and it can justify the spatial distribution of the evaporation amounts of these high-elevation dimictic lakes on the TP.

Future change and significance of lakes’ evaporation amounts

The climatic average evaporation amounts have been determined in this study; however, future changes in the evaporation amount of these high-elevation lakes remain to be clarified. Under global climate warming, later freezing at a rate of 5.8 days per 100 years and earlier breakup at 6.5 days per 100 years have been observed by long-term (1846–1995) in situ ice phenology measurements over 39 lakes and rivers in the Northern Hemisphere (27). In Canada, trends toward earlier breakup dates for most lakes from 1951 to 2000 have also been reported (24). Moreover, 3.7% of the world’s lakes in the southern reaches of the winter ice cover region may change from ice-covered systems to open-water systems, and lakes with annual ice cover may shift to intermittent ice cover; thus, warming-induced northward expansion will occur (23), and the ice cover loss may occur at an accelerated rate. However, the trend of the DIF in lakes of the TP appears to be highly variable from MODIS data (20). In our results, a high positive correlation coefficient (0.708**) between air temperature and DIF has been reported; thus, the durations of the ice-free season of these high-elevation lakes will likely present an increasing trend due to the accelerated warming. The SC days and snow depth decrease in the north and northwest regions of the TP (28), and it is an additional evidence corresponding to the earlier breakup dates of these high-elevation lakes. Thus, an increasing rate of lake ice loss together with a growing DIF and increasing evaporation amounts can be deduced for these high-elevation lakes.

Further, a perspective from physical explanation is that more heat can be absorbed in lakes with a longer ice-free season because of the much lower albedo of water relative to snow and ice, thus leading to a higher net radiation for heat allocation. The enhanced solar forcing can lead to lake warming, i.e., in Lake Nam Co (29), which may further prolong the nonfreezing season, thus forming a positive cyclic feedback. These conditions have been proved to promote the shift of perennially frozen lakes to temporarily frozen lakes (23). A model simulation has also suggested that accelerated lake evaporation amounts can occur by the combined effects of shorter periods of ice cover in a warmer climate and decreased Bowen ratios in water bodies (30). In addition to the increased evaporation amounts of lakes, the increased heat gains can accelerate regional temperature rise through the warm lake effect, provide more water vapor for regional weather variations, and result in more intensive hydrological cycles. The effects of lakes on surrounding climate, including increased evaporation, lake-effect snow, and thermal moderation in the surrounding areas, are expected to be more prominent not only in high-latitude areas (21) but also on the high-elevation lake zones.

Given the importance of TP water resources to billions of people and numerous ecosystems downstream, the variations and total amounts of hydrological elements (including permafrost, glacier melt, precipitation and snowfall, river discharge, atmospheric water vapor, evapotranspiration over land and evaporation over lakes, etc.) are being explored by scientists all over the world with the support of Chinese Academy of Sciences Strategic Priority Research Program entitled “The Pan-Third Pole Environment Study for a Green Silk Road (Pan-TPE)” and the Second Tibetan Plateau Scientific Expedition and Research program. Recently, the ground ice volume (liquid water equivalent) in the permafrost is estimated to be up to 1.27 × 104 km3 (31), while it is approximately 4 × 103 km3 for the total liquid water equivalent of the glacier (32). Further, approximately 80 km3 of ground ice in the permafrost will melt every 10 years and then participate in the regional water cycle, and a total annual mass budget of −15.6 ± 8.9 km3 year−1 over 80% of the glacier area was obtained by using ICESat data during 2003 to 2009 (33). Our research is the first to provide the total evaporation amount of these high-elevation lakes, with a high accuracy value of 51.7 ± 2.1 km3 year−1. Urgent efforts are being given to the estimation of total amounts of other hydrological components, such as rivers’ discharge, snow melt, precipitation, atmosphere water vapor transportation, and evapotranspiration over the heterogeneous land surface, to achieve catchment-scale water balance and even over the whole TP. Moreover, our method has a good prospective to be applied to other large ice-covered lakes worldwide. Precise estimation of lake evaporation over lakes on the TP and other regions will definitely be important for regional and global water balance analysis and have important implication for future climate change.


Based on a plausible hypothesis that the change in heat storage during the ice-free season will be close to zero and the winter sublimation during the ice-covered season will have a constant average value, we estimated the evaporation amounts of the 75 high-elevation large dimictic lakes by combining meteorological and satellite data. The ice phenology events in these lakes can be identified through satellite data, which can clearly distinguish icing and non-icing periods. A lake with a smaller area at higher elevation and higher latitude generally has a later fully ice-free date. After evaluation and correction of the meteorological data and satellite products, we used in situ EC measurements of four high-elevation lakes to validate our data with the estimated Bo and evaporation amounts. The estimated Bo shows similar seasonal variations and close average values against in situ EC measurements, while the evaporation amounts also show acceptable accuracy. The spatial distribution of ice phenology, meteorological variables, and evaporation amounts shows large variability, and higher evaporation amounts are spatially distributed in the SLG and OLG compared with the NLG. The total water evaporated from the lake surfaces to the atmosphere is estimated at 29.4 ± 1.2 km3 year−1 for the 75 studied lakes and 51.7 ± 2.1 km3 year−1 for all lakes of the TP. Our results provide a reference lake evaporation dataset for these high-elevation large dimictic lakes and can reduce the evaporation-induced bias in lakes’ water budget analysis. The proposed method also shows quite good application prospects for other large dimictic lakes worldwide.


Satellite products

MODIS is a multispectral sensor onboard the sun-synchronous near-polar orbital satellites Terra and Aqua, and it can monitor Earth’s surface twice per day and thus generate four observations per day at local overpass times of 10:30 and 22:30 for Terra and 13:30 and 1:30 for Aqua. Terra and Aqua products have been available since March 2000 and July 2002, respectively. A number of MODIS products have been developed, including land surface temperature (LST) products (MOD11A2 and MYD11A2, level 3, 8-day average) and SC products (MOD10A2 and MYD10A2, level 3, 8-day average). MODIS LST products are produced using the generalized split-window algorithm or the day/night algorithm, which have good accuracy for LST estimates over lake surfaces. Thus, the spatial variations of LST over different lakes have been studied worldwide. MODIS SC products are produced by a Normalized Difference Snow Index and several other criteria tests that minimize the impact of cloud cover (34). The classification scheme takes advantage of visible and near-infrared wavelengths with high reflectivity for ice and low reflectivity for water and further divides lake surface into five classes: snow, no snow, water, ice, and cloud. The ice phenology of lakes on the TP has been identified with an acceptable accuracy (20). Moreover, the digital elevation model of the ASTER GDEM (Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model, version 3) with a spatial resolution of 30 m is used for lake boundary shape extraction.

CMFD dataset

A high spatial resolution (0.1°) forcing dataset (namely, CMFD) covering China was developed by the Institute of Tibetan Plateau Research, Chinese Academy of Sciences (35). The CMFD can provide meteorological variables, including air temperature, relative humidity, air pressure, Rs↓, and Rl↓, at a temporal resolution of 3 hours from 1979 to 2016, and it merges a variety of data sources, including (i) 740 meteorological stations organized by the China Meteorological Administration, (ii) Global Land Data Assimilation Systems, (iii) GEWEX Surface Radiation Budget products, and (iv) Princeton forcing data. Thus, after a proper correction of CMFD data in the studied area, the CMFD has been widely used for forcing in climate change studies over the TP (10, 29).

In situ measurements

EC measurements can directly observe lake-air sensible and latent heat flux (H and LE, W m−2), and the published values in Lake Nam Co, Lake Qinghai, Lake Serling Co, and Lake Ngoring are collected to evaluate the seasonal variations of Bo and the evaporation amounts during their observational periods. To further estimate winter sublimation, EC measurements during ice-covered seasons from 2016 to 2017 in Lake Nam Co (9) and from 2012 to 2013 for a small lake adjacent to Lake Nam Co (36, 15) are presented, where LE from the wind sectors of ice surface are averaged with a consideration of the influences of the footprint and data quality. In situ measurements of water surface temperature at a depth of 0.5 m by a temperature probe (HOBO water temperature Pro v2 021001, Onset, USA) from July to November in 2015 and 2016 are collected for evaluation and correction of satellite LST data. Furthermore, air temperature at a height of approximately 10 m by an Automatic Weather Station on the island of Lake Nam Co is collected and can be used for correcting the warm lake effect. Configurations of these in situ measurements can be found in (9). Moreover, in situ measurements of air temperature, humidity, pressure, Rs↓, and Rl↓ are collected from 14 stations (Mustugh, Ali, Qomos, Linzhi, Shuanghu, Nam Co, MS3608, BJ, NPAM, Amdo, D66, D110, D105, and NewD66) from TORP (locations are shown in Fig. 1A) for evaluation and correction of the CMFD. The 14 TORP stations are independent stations, and the instrument configurations of the TORP sites can be found in (19).

The energy budget of a lake can be expressed as Rn = H + LE + Gs + Gb + Ga, where Gb (the heat transfer between water and bottom sediments, W m−2) and Ga (the net energy gained or lost by a lake due to the exchange of water masses resulting from inflow-outflow balance, W m−2) are minor for these endorheic large lakes and negligible relative to other items. The heat storage in the water (Gs, W m−2) can also be ignored for the complete ice-free season, when the heat storage after ice melt is mostly released before ice forming (15). Thus, lake evaporation amounts during the ice-free season can be roughly estimated by using net radiation (Rn, W m−2) and Bo: LE=Rn(1+Bo), where Rn=(1α)Rs+RlεσTs4, and Bo=HLE=γTsTaEsEa. Rs↓ (W m−2) and Rl↓ (W m−2) are the downward shortwave radiation and downward longwave radiation, respectively; α = 0.06 is the albedo of water; ε = 0.98 is the emissivity of water; σ is the Stephen-Boltzmann constant; Ts (K) is the water surface temperature; γ is a psychrometric constant; Ta (K) is the air temperature; Es (hPa) is the saturated vapor pressure at Ts; and Ea (hPa) is the vapor pressure of air. Thus, meteorological variables (Rs↓, Rl↓, Ta, and Ea), the Ts, and the lake ice phenology are key parameters for estimating lake evaporation amounts during the ice-free season, and these variables can be obtained from the CMFD and satellite data.

The data processing procedures

In this study, 75 lakes on the TP are selected following (37). The lakes, all with areas larger than 88 km2, have a total area of approximately 26,450 km2, which accounts for approximately 56.9% of the lake surfaces over the whole TP. The multiyear (2003–2016) average ice phenology and Ts are processed by using MODIS SC and LST products, respectively. The following seven procedures are included.

The first procedure is the “lake boundary shape extraction.” The boundary shape of a lake is constructed by using digital elevation model of ASTER GDEM v3, which is mosaicked over the whole TP by ArcGIS software and then resampled into a spatial resolution of 1 km. The boundary shape of each lake, characterized by same elevation value and clustered into a specific shape, is extracted manually one by one. The boundary shapes of the studied lakes are shown in blue color in Fig. 1A and can be used to precisely derive the meteorological variables, MODIS SC, and LST data of each lake.

The second procedure is the “ice phenology identification.” The MODIS Reprojection Tool is used to mosaic and resample selected tiles of MODIS products to a GeoTIFF format at a pixel size of 1 km and with a geographic projection. The lake ice phenology can be identified by MODIS SC 8-day composite data, which can distinguish the lake surface into classes of ice/snow and water. The ice phenology of a specific lake can be obviously distinguished into four events: FO (the date when a detectable ice appears), FU (the date when the lake is fully ice covered), BU (the date when a detectable water surface appears), and WCI (the date when the lake is fully ice-free). In general, these ice phenology events can be determined by percentage values of approximately 3 and 80% by ice and snow pixels to the highest value of ice and snow pixels in majority of the lakes. However, because of the influences of cloud, water, and filling pixels in MODIS SC products, an additional visual interpretation of the multiyear ice phenology events of a specific lake is needed. The DIF can be determined by two periods, with one from WCI to FO and the other from BU to FU, and the former one is used as the DIF in this study without additional information. The ice phenology events can be cross-validated by the variation of water percentage values and further visually checked over each lake. The multiyear (2003–2016) average values are used to represent the climatic average to enhance the confidence of the ice phenology results over these large dimictic lakes.

The third procedure is the “lake surface temperature extraction.” Ts is an important variable in describing lake-air interaction processes and can be evaluated with an acceptable accuracy to monitor its spatial, diurnal, and seasonal variations in lakes where in situ measurements are rarely available. In this study, the representative Ts of each lake is obtained as follows: selecting all the pixels of Ts with the best data quality using the boundary shape file, sorting these values by numerical size, removing the lowest and highest 10% of the data to exclude the influences of land-contaminated pixels or pixels with abnormal values, and averaging all the remaining values. Because of high thermal capacity of large lakes, the spatial variability and the seasonal variations of Ts should not present marked changes. However, marked decreases of Ts can be found in MODIS night observations but not in MODIS day observations, which is similar to the results (29). Thus, MODIS day observations are averaged to represent Ts, and the daily average in situ Ts, which can minimize the impacts of surface warming during the day and surface cooling at night, is used to correct satellite-derived LST data.

The fourth procedure is the “correction of CMFD.” The meteorological variables Rs↓, Rl↓, Ta, and Ea and air pressure from the CMFD are corrected by using in situ measurements from the 14 TORP sites. The inconsistencies of elevation differences between the CMFD and ASTER GDEM are first corrected by using the adiabatic lapse rate for air temperature and the equation [P = 1000 * (2.406 − 0.0000534 * h)5.26, where h is the elevation value, in meters] for air pressure, respectively. Because Ta in a large lake is generally different from that alongside the lake. Ta at a height of approximately 10 m in Lake Nam Co is compared between one station on the island and another station alongside the lake, and the pattern of differences in seasonal variations is used to correct for the warm lake effect of these large lakes.

The fifth procedure is the “determination of the uncertainty of evaporation amounts.” To estimate the uncertainty of evaporation amounts during the ice-free season (Eo) and annually (Ea) at each lake, we summed the biases caused by uncertainties in ice phenology and the biases from uncertainties in CMFD forcing. The SDs of Eo (or Ea) estimated by two DIF are considered the uncertainty from the ice phenology. The uncertainties in CMFD forcing result from the biases from air pressure (±10 hPa), Ta (±0.75 ° C), Ts (±0.75 ° C), air humidity (±5%), Rs↓ (±2.5%), and Rl↓ (±2.5%), with all the values in parentheses representing half of the largest bias in the CMFD before correction. The uncertainties from all these error sources are summed in quadrature to estimate the overall uncertainty.

The last procedure is the “EC data processing in Lake Nam Co.” EC measurements from the island of Lake Nam Co from July 2016 to July 2018 are processed to construct a long evaluation data series excluding the influences from the footprint, quality, and instrument malfunction. By distinguishing EC measurements over different wind sectors, with one sector for observations over the lake surface (Hlake, LElake, Tlake, Ulake, and Elake) and the other sector for observations influenced by land (Hother, LEother, Tother, Uother, and Eother), the results show that the monthly turbulent heat flux and meteorological variables over the two wind sectors are quite different (fig. S2). Compared with the monthly Hother, the monthly Hlake is much lower during May to August and much higher during the winter seasons (fig. S2A). The monthly LElake is much higher than the monthly LEother, especially from August to January (fig. S2B). The monthly Ulake behaved similarly as the monthly LElake (fig. S2C). The monthly air temperature and humidity from the two wind sectors are quite similar (fig. S2, D and E), except for the air temperature during cold months. H (or LE) is mostly determined by the wind speed multiplied by the temperature (or humidity) gradients, and the seasonal variation of Ts should have no marked changes because of the large thermal capacity of large water bodies. Thus, in situ measurements of the turbulent heat flux from land sector can be simply corrected by the following equations: Hcor=(TsTother)×Uother(TsTlake)×UlakeHlake and LEcor=(EsEother)×Uother(EsElake)×UlakeLElake, where Ts is from satellite measurements and Es is the saturated water vapor pressure at Ts. The corrected monthly H and LE can be lastly determined by summing up Hcor (or LEcor) and Hlake (or LElake), considering their relative proportions. The seasonal variations of monthly H and LE are shown in fig. S2F. The turbulent heat fluxes from July to November are close to those estimated by a combination of the bulk transfer method and in situ measurements (9).


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: Funding: This research was funded by the Second Tibetan Plateau Scientific Expedition and Research (STEP) program (grant no. 2019QZKK0103), the Strategic Priority Research Program of Chinese Academy of Sciences (XDA20060101), the National Natural Science Foundation of China (41705005, 91837208, and 41830650), the China Postdoctoral Science Foundation. The ASTER GDEM v3 products and MODIS products can be downloaded freely from and, respectively. We would like to thank the anonymous reviewers and the Editorial team of Science Advances for their good suggestions and comments. Author contributions: B.W., Y.M., and Z.S. conceived this study and wrote the manuscript. Y.W. and W.M. helped for conducting in situ measurements and processing the data. All authors contributed ideas to the project. Competing interests: All authors declare that they have no competing interests. Data and materials availability: The meteorological data of the Nam Co station are accessible from the Third Pole Environment Database ( All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article