A Quantification of Heat Storage Change-Based Evaporation Behavior in Middle–Large-Sized Lakes in the Inland of the Tibetan Plateau and Their Temporal and Spatial Variations

Abstract

A large number of different-sized lakes exist in the inland area of the Tibetan Plateau (TP), which are examples of the important connection between the atmosphere and hydrosphere through the analysis of lake surface convergence and evaporation processes. The evaporation level changes that occur in middle–large-sized lakes (surface area > 50 km²) in the area directly influence the regional mass and energy balance values, atmospheric boundary layer heat and humidity structures, and weather processes occurring in the lower-reach areas. The studies conducted in the literature at present, concerning lake evaporation processes, generally overlook the differences in lake heat storage behavior due to the reduced amount of data in the literature concerning lake bathymetry. According to the in situ bathymetric data obtained for 68 middle–large-sized lakes in the inner basin of the TP, in this study, we calculated their heat storage (G) change values by using the different vertical-depth water-temperature-change integral method, and we established a regression equation for the heat storage and lake surface net radiation values for 68 lakes. The evaporation rates of 134 middle–large-sized lakes larger than 50 km² in the inland are of the TP were calculated by obtaining the G regression result and adopting it into the Penman model, as well as estimating the evaporation losses of theses 134 lakes from 2002 to 2018. The result shows that the annual average evaporation rate for these lakes is 927.39 mm/year, with an insignificant upward trend (0.10 mm/year). This method achieved good accuracy compared with the Bowen ratio method, which estimates the evaporation rate during the ice-free season, with a high correlation coefficient (R) value of 0.95 and least root mean square error (RMSE) value of 61 mm. The annual mean evaporation rate can be divided into the southern and northern lake groups along a 34°N line with a difference of 314.41 mm/year. The annual average evaporation volume of these lakes was 25.02 km³ and showed an upward trend of 0.35 km³/year. Among them, the annual average evaporation volume contribution ratio of level-1 lakes (50 km² ≤ lake’s area < 100 km², 61 lakes) was 14.04%, showing an upward trend, and the contribution of level-3 lakes (lake’s area ≥ 500 km², 10 lakes) was 41.50%, showing a downward trend. There were no obvious changes in the level-2 lakes (100 km² ≤ lake’s area < 500 km², 63 lakes), which maintained at the same level in approximately 44.46%. Air temperature is the most important factor affecting the evaporation rate of lakes, while the lake surface area is the main factor affecting lake evaporation volume. Our study, considering the actual lake heat storage value, provides a useful reference for further improving lake water budget balance values and watershed hydrologic features in the inland closed lakes located in the TP.

1. Introduction

Water resources are an important fundamental and strategic natural resource for the survival and development of human beings that, through precipitation, evaporation, runoff, soil water, and other natural water cycle processes, allow us to realize the evolution and transformation of lakes [1]. Lakes are an important part of ground–surface water bodies and can provide sensitive and clear feedback regarding the occurrence of local hydrological cycle and regional climate change, especially in the Tibetan Plateau (TP), an environment with a high altitude, high radiation rates, and dry climate [2,3]. The TP, known as the “Third Pole” and “Asia Water Tower” of the Earth [4,5], possesses a total area of lakes accounting for more than 50% of the total area of lakes in China [3]. At the same time, the TP contains more than 1200 lakes with an area larger than 1 km², and they have shown an obvious expansion during the past 30 years [6]. Many studies in the field have qualitatively and quantitatively analyzed the impacts on the lake expansion rates in the TP from glacial melt water and precipitation levels [7,8,9,10,11,12,13]; however, the quantitative analyses conducted on the contribution of lake surface evaporation behavior to water balance levels in the TP are still not detailed enough. Therefore, it is difficult to accurately assess the response and feedback factors of lake water storage changes caused by climate change.

Previous studies have shown that the annual average evapotranspiration rate occurring over the TP is approximately 328 mm/year, and there is an obvious upward trend [14,15,16]. Compared with other land-cover types, lake water has higher transparency, lower albedo, lower surface roughness, and higher heat capacity and thermal conductivity values [17]; therefore, a higher evaporation rate compared to other land cover types exists in the TP [18]. At present, there are two main methods used in the research to study the evaporation rates of lakes located in the TP: in situ observations and remote sensing estimations.

The observation method mainly relies on evaporation pans and the eddy covariance (EC) observation systems to monitor water surface evaporation levels [19]. Li et al. (2001) calculated the annual average evaporation rate of Lake Zigetang by inputting the long-term meteorological observation data into the Penman formula, and it was approximately 950 mm/year exhibiting a downward trend [20]. Li et al. (2007) estimated that the average annual evaporation rate of Lake Qinghai from 1959 to 2000 was about 924 mm/year using the pan evaporation method [21]. By adding the meteorological observation data into Penman–Monteith equation, Zhu et al. (2010) and Zhang et al. (2011) estimated that the annual evaporation rates of Lake Nam Co were 1430 and 1184 mm/year [22,23]. Due to the difference in heat storage levels between lakes and land, considerable errors between evaporation pan observation, land meteorological observation, and real lake evaporation results were evident; the EC observation system combined with the energy balance method was gradually used to monitor lake evaporation values [24]. Li et al. (2015), through a flux analysis of Lake Ngoring, using the EC observation system combined with the Bowen ratio energy balance method, observed that the intrusion of dry and cold air with strong wind could result in significant increases in sensible heat flux (H) and latent heat flux (LE) values that were approximately 2.0–4.5 times as much as those during normal days, and the heat storage change in the lake water levels dramatically decreased and could even provide 70% of the energy for lake surface H and LE [25]. The same method was used to monitor the evaporation rate for Lake Qinghai, and the results showed that the annual evaporation level for Lake Qinghai was approximately 830 mm/year [26].

As a result of the limited coverage ability and high maintenance cost of the observation stations, the effect of lake microclimate and heat storage change activity in middle–large-sized lake on lake surface evaporation levels could not be adequately reproduced in this research [27]. The remote sensing estimation method can make up for the shortcomings of the observation method through the use of modularized data, a lake model, and a parameter scheme [28]. Lazhu et al. (2016) used the China meteorological forcing dataset (CMFD)-driven Flake model and Bowen ratio energy balance method to estimate the annual evaporation rates for Lake Nam Co, which was about 832 mm/year and presented an obvious upward trend [29]. Ma et al. (2016) estimated that the annual average evaporation rate for Lake Nam Co was approximately 632 mm/year by using the CMFD data-driven CRLE model, and a downward trend was exhibited [30]. Wang et al. (2020) used CMFD data and the Bowen ratio energy balance model to quantitatively analyze the evaporation rated for 75 large lakes in the TP, and the results showed that the average annual evaporation rate was about 29.4 ± 1.2 km³/year [31].

The abovementioned studies show that, regardless of whether the Penman or Bowen ratio methods are used by the researchers, lake heat storage change has a considerable influence on the accuracy of lake evaporation outcomes. Since middle–large-sized lakes have a greater impact on evaporation levels compared to small lakes due to their higher heat storage capacity, it is key to estimate the heat storage levels of these lakes for accurately quantifying lake evaporation levels in this research [32]. However, due to the limited availability of in situ bathymetric data for most lakes located in the TP, it is difficult to accurately estimate lake heat storage capacity according to the lake water-temperature vertical integral method. In this study, we establish an empirical procedure that can rapidly estimate the lake heat storage change value for a group of lakes with different characteristics based on the relationship between lake heat storage calculated using the lake water-temperature vertical integral method and available remote sensing meteorological variables in the lakes with bathymetric data. Using this method, we analyze the evaporation levels of an increased number of middle–large-sized lakes and elucidate the spatio-temporal-variation characteristics of evaporation occurrence in different types of lakes located in the inland area of the TP.

2. Study Area and Data

2.1. Study Area

The inland area of the TP that we studied covered a total area of approximately 7.08 × 10⁵ km² with an average elevation level greater than 4960 m, and it accounts for 1/3 of the TP’s total area [33,34]. The lakes in the area account for the highest lake density level, and most lakes are closed without outflow runoff activity. Water loss level is mainly composed of evaporation, which has a significant impact on climate change [35,36]. Human activity is reduced in this area; therefore, lake evaporation levels can reflect the impact and response of natural processes [37]. By 2018, 1032 lakes had an area greater than 1 km², and the total area of these lakes was 35,243.27 km² in the inland area of the TP. Among them, the total area of 143 middle–large-sized lakes (lake surface area greater than 50 km²) was 28,738.75 km²; this accounted for approximately 81.54% of the total lake area in the inland of the TP and 56.92% of the total lake area in the entire TP.

Previous studies have discovered that the evaporation process of lakes is affected by many factors, such as the latitude and longitude, the climate condition in the lake basin, lake water physicochemical property, and human activities [38]. However, the similar geographical and climatic environments of lake groups composed of neighboring lakes leads to results of similar evaporation characteristics occurring in the same region [39]. Therefore, according to the spatial positions, the main land-cover types in the lake basin, and the average elevations of middle–large-sized lakes in the inland area of the TP, we divided all the middle–large-sized lakes into seven relatively concentrated lake groups to estimate and analyze lake evaporation levels. The summary information and spatial distribution patterns of lake group zones are shown in

Table 1. Summary of lake group zones.

Zone	Proportion of Land-Cover Types					Elevation (m)
	Water	Open Shrublands	Grasslands	Snow and Ice	Barren Vegetated
S01	1.38%	0.00%	4.90%	0.27%	93.45%	4020
S02	1.06%	0.24%	4.47%	2.50%	91.74%	5119
S03	2.47%	0.00%	15.84%	1.09%	80.59%	4979
S04	2.14%	0.27%	62.90%	0.89%	33.80%	5096
S05	2.45%	4.19%	50.03%	0.18%	43.16%	5056
S06	5.19%	0.23%	86.01%	0.02%	8.55%	4984
S07	7.34%	0.00%	87.97%	0.06%	4.64%	4847

and Figure 1. In S01, S02, and S03 regions, the main land-cover type is barren vegetated; however, the areas have large average elevation differences between them. S04 and S05 regions have similar elevations; however, they have different latitudes, and the main land-cover type is a mixture of grasslands and barren vegetation with different ratios. S06 and S07 regions possess the same main land-cover type of grasslands and the same latitudes; however, they present a 137 m difference in elevation.

In this study, the middle–large-sized lakes in the inland area of the TP are defined as lakes with water surface areas greater than 50 km². The lake area data of the TP were obtained from Zhang et al. (2019) [40] and are available from the Tibetan Plateau Scientific Data Center. In order to further clarify the target lakes and match the remote sensing data, 134 lakes with areas greater than 50 km² for five consecutive years from 2002 to 2018 were selected as the final research target lakes. Detailed information of all 134 lakes can be found in the Supplementary Materials (Table S1). The spatial distribution values are shown in Figure 1. Additionally, variations in the middle–large-sized lakes in the inland area of the TP from 1979 to 2018 are presented in the Supplementary Materials (Figure S1).

2.2. Data

2.2.1. In Situ Meteorological Station Data

The in situ meteorological data were used to calibrate the downscaled data obtained from the China Meteorological Forcing Dataset (CMFD) [41] and MODIS land surface temperature datasets (MODIS LST) [42], which were used to drive the Penman formula and lake models. The spatial distribution values of the meteorological observation stations presented in this paper are shown in Figure 2, and the detailed information of each station is shown in

Table 2. List of the geographic characteristics of the 11 Tibetan Observation and Research Platform (TORP) sites and observation variables.

Site	Latitude, Longitude	Elevation (m)	Period	Land Cover	Variables (Units)	References
MAWORS	38.42°N, 75.03°E	3668	2005–2016	Alpine desert	2 m-Air temperature (°C) 0 m-Soil temperature (°C) 10 m-Wind speed (m/s) 2 m-Humidity (%)Radiations (w/m2)	Ma et al., 2020 [43]
NADORS	33.39°N, 79.7°E	4270	2005–2016	Alpine desert
QOMS	28.36°N, 86.95°E	4298	2005–2016	Alpine desert
NAMORS	30.77°N, 90.96°E	4730	2005–2016	Alpine steppe
NAQU	31.37°N, 91.9°E	4509	2005–2016	Alpine meadow
SETORS	29.77°N, 94.74°E	3327	2005–2016	Alpine meadow
SHENZHA	30.95°N, 88.7°E	4750	2016–2018	Wetland		Wei et al., 2021 [44]
BATANG	32.85°N, 96.95°E	4003	2017–2018	Meadow
DASHALONG	38.84°N, 98.94°E	3739	2015	Wetland
AROU	38.03°N, 100.45°E	3033	2015	Meadow
YAKOU	38.01°N, 100.24°E	4148	2015	Meadow

. The applied meteorological data include temperature at a 2 m altitude, relative humidity at a 2 m altitude, wind speed at a 10 m altitude, and surface radiation data (including downward short-wave radiation, downward long-wave radiation, and net radiation). The in situ daily meteorological data were used to calibrate the CMFD, and a monthly mean value of 0 m soil temperature was used to calibrate MODIS LST (MOD11A2 and MYD11A2).

2.2.2. Lake Spatial Grid and Meteorological Datasets

A global surface water dataset (GSWD) [45] was used to create lake basin grids for middle–large-sized lakes with bathymetric data to perform water bathymetric data interpolations in the inland area of the TP. In order to obtain lake heat storage change results, resampled CMFD data were used to create a lake model for each lake to obtain the water temperature level results for different water depths. Resampled MODIS land surface temperature (LST) datasets (i.e., MOD11A2 and MYD11A2) were combined with CMFD data to estimate lake surface net radiation levels. Detailed information about remote sensing datasets is summarized in

Table 3. Data used for lake surface evaporation levels.

Data Name	Spatial Resolution	Temporal Resolution	Purpose	Web Link
Lakes larger than 1 km2 in TP dataset [40]	Shapefile	1–10 year	Lake mask	http://data.tpdc.ac.cn (accessed on 1 May 2023)
Global surface water dataset occurrence (GSWD) [45]	30 m (resample to 0.01°)	Time invariant	Water area grid extraction and interpolation boundary	https://global-surface-water.appspot.com (accessed on 1 May 2023)
China Meteorological Forcing Dataset (CMFD) [41]	0.1°	Daily	Driving lake and Penman models	http://data.tpdc.ac.cn (accessed on 1 May 2023)
MODIS Terra LST (MOD11A2) [42]	1 km (resample to 0.01°)	8 days	Water surface temperature and net radiation	https://search.earthdata.nasa.gov (accessed on 1 May 2023)
MODIS Aqua LST (MYD11A2) [42]	1 km (resample to 0.01°)	8 days	Water surface temperature and net radiation	https://search.earthdata.nasa.gov (accessed on 1 May 2023)

.

The spatial meteorological data collected from the CMFD, including 2 m temperature (CMFD TEMP), 2 m relative humidity (CMFD RHUM), downward long-wave radiation (CMFD LRAD), and downward short-wave radiation (CMFD SRAD), are key variables of the Penman formula and lake models that are used to estimate lake evaporation rate and simulate lake water temperature results. Figure 3 shows the correction equations for these variables using in situ measurements obtained from 11 TORP stations. The accuracy test indexes mainly include correlation coefficients (R), root mean square error (RMSE), and percent BIAS (PBIAS). The specific calculation methods used can be observed in the Supplementary Materials.

Because the TP is characterized by high altitude, low temperature, and long winter period, lake ice can block the energy exchange occurring between the water and atmosphere and thus influence lake surface evaporation levels. The ice-cover degree can be divided into ice-free, ice and water mixing, and freeze-up periods with a decrease in the temperature [31]. Therefore, MODIS LST datasets (MOD11A2 and MYD11A2) were resampled as 0.01° × 0.01° of spatial resolution, and the monthly time resolution was adjusted and used to determine whether the lake ice blocked the energy exchange occurring between the water column and atmosphere. Similar to the CMFD meteorological data, the accuracy values of MODIS LST data were also tested and corrected by using in situ measurements of land surface temperature data. The results show a correlation coefficient of 0.94, RMSE value of 3.04 K, and PBIAS of 0.05%, and a linear fitting equation of y = 0.9334x + 18.6346 (Figure 4) can decrease the bias of the monthly mean MODIS LST dataset.

2.2.3. Bathymetric Data of Lake

Numerous studies used water-depth values to calculate lake heat storage change occurrence using the Penman formula due to a lack of spatial bathymetric data for the lake; however, these depth values can lead to uncertain errors in the real heat storage and evaporation values [24,29,31,40,42]. In this paper, we used the bathymetric data of 68 lakes with areas greater than 50 km², which were surveyed by our research group to estimate lake water heat storage values. The bathymetric device used was a Lowrance HDS5, which had a vertical accuracy of 0.01 m, and bathymetric data were recorded once per second. The measured route examples for 14 out of the 68 lakes are presented in Supplementary Materials Figure S2. The underwater lake topography was established with grid units of 0.01° × 0.01° using spatial analyst tools (Topo to Raster) in ArcGIS 10.2 for the 68 lakes. Topo to Raster is an iterative finite difference interpolation technique designed for the creation of a hydrological correct digital elevation model (DEM), which was developed by Hutchinson (1988, 1989) [46]. Previous studies have shown that the accuracy of this method can achieve average error values of less than 1.6 m [39,47], which was suitable for use in our research to estimate lake heat storage changes as the core parameter of algorithms concerning open water surface evaporation levels.

3. Methodology

3.1. Algorithm for Evaporation Rate

The Penman formula is a classical method with good precision and physical significance outcomes used in this research to estimate evaporation rate levels and is widely used in open water all over the world, including the TP [42,48]. Therefore, the Penman formula was selected as the best method to estimate the evaporation rate levels in middle–large-sized lakes located in the inland area of the TP. The evaporation volume of lakes was obtained through multiplying the evaporation rate by the lake area. The core equation, calculation of important parameters, and value of constant term of the Penman formula are shown in the following equation:

$$E_{rate}=\underset{E_R}{\underbrace{\frac{\Delta \left(R_n-G\right)}{{\lambda }_{\nu }\left(\Delta +\gamma \right)}}}+\underset{E_A}{\underbrace{\frac{\gamma f\left(\mu \right)\left(e_s-e_a\right)}{{\lambda }_{\nu }\left(\Delta +\gamma \right)}}}$$

(1)

where $E_{rate}$ is the lake or reservoir evaporation rate (mm); $\Delta$ is the slope of the saturation vapor pressure curve (kPa/°C); $R_n$ is the net radiation (W/m²); $G$ is the heat storage change in the lake water (W/m²); $\gamma$ is the psychrometric constant (kPa/°C); $f\left(\mu \right)$ is the wind function; $e_s$ is the saturated vapor pressure at air temperature (kPa); $e_a$ is the air vapor pressure at air temperature (kPa); and ${\lambda }_{\nu }$ is the latent heat of vaporization.

$E_{rate}$ consists of two components: the aerodynamic component $E_A$ and the radiative component $E_R$. For the calculation of $E_A$, the wind function $f\left(\mu \right)$ is the horizontal aerodynamics function. Here, we used a TP generally applicable wind function Equation (2), which was presented by Lin et al. [18]:

$$f\left(u\right)=0.26\left(1.0+0.536u_2\right)$$

(2)

where $f\left(u\right)$ is the wind function; $u_2$ is the wind speed at the height of 2 m (m/s).

The meteorological dataset used in this paper was obtained from the CMFD and the wind speed data height of the CMFD was 10 m. Equation (3) was used as the different height wind speed conversion formula in this paper:

$$u_2=u_z·\frac{4.87}{ln\left(67.8z-5.42\right)}$$

(3)

where $u_z$ is the wind speed at the height of z m (m/s). $R_n$ and $G$ are two core parameters in the $E_R$ section. In this paper, we used the radiation balance method to obtain the net radiation value, which is expressed as Equation (4) [49]:

$$R_n=\left(1-a\right)R_{s\downarrow }+R_{l\downarrow }-R_{l\uparrow }$$

(4)

where $R_{s\downarrow }$ is the total downward short-wave radiation (W/m²), a = 0.055 is the albedo of water, $R_{l\downarrow }$ is the downward long-wave radiation amount (W/m²), and $R_{l\uparrow }$ is the upward long-wave radiation amount (W/m²), $R_{l\uparrow }=\epsilon \sigma T_S^4$; ε = 0.98 is the emissivity coefficient of water, $\sigma =5.67\times {10}^{-8}$ (W/m²K²) is the Stephen–Boltzmann constant, and $T_S$ is the lake water surface temperature value.

3.2. Algorithm for Lake Heat Storage Change

In the Penman formula, meteorological driving parameters (such as air temperature, air pressure, wind speed, humidity, and radiation) can be obtained in a wide variety of ways; however, it is very difficult to accurately obtain the change in lake heat storage (G) values, which is one of the core parameters required. The energy balance method (G = Rn − H − LE) replacing lake heat storage change occurrence can reduce the error of evaporation estimation results to a certain extent, when there is a lack of water depth data [50,51]. However, this method is still unprecise due to the energy including not only lake water, but also lake sediments and groundwater [49,50]. At present, the most widely used method in the world still attempts to obtain lake heat storage change values based on the vertical integral of lake water temperature change at different depths. However, on the one hand, a lack of detailed bathymetric data of lakes in the inland of the TP causes the water temperature levels at different lake depths to be hardly stimulated by combining remote sensing data and the lake model. On the other hand, it is difficult to conduct in situ surveys of vertical water temperature in more middle–large-sized lakes. Considering the abovementioned factors, we used two schemes to calculate the lake heat storage values of all 134 middle–large-sized lakes (68 lakes with and 66 lakes without bathymetric data) in the study area. The time-series data for detailed lake heat storage changes were implemented in the following 4 steps, and the created flowchart is presented in Figure 5.

First, the time-series lookup table of vertical water temperature profiles of 68 middle–large-sized lakes with bathymetric data was obtained by using the CMFD data to drive the LAKE model (V2.3). Then, the freezing periods and ice-coverage percentages of these lakes were obtained by using the CMFD data to drive the Hostetler model.

Second, by using each lake water temperature profile lookup table result combined with the in situ bathymetric interpolation data, each water column heat storage value was estimated by the vertical integrated lake water temperature using a one-by-one pixel. This is represented by Equation (5) [52]. In winter, lake ice will block the energy exchange between water and atmosphere because the sublimation process will replace the evaporation process. The process of determining lake ice is divided into three steps. Firstly, the Hostetler model was forced by CMFD data in order to obtain the monthly average whole lake ice cover percent (ICP). An ICP of less than 10% of the lake state is determined to be open water; an ICP greater than 90% of the lake state is determined to freeze up; an ICP between 10% and 90% of the lake state is determined to be an ice–water mixing state. When the lake state is determined as open water and freeze up, evaporation and sublimation rates were introduced in Section 3.1 and Section 3.4, respectively. Then, the corrected MODIS monthly average LST data product was used to determine the ice-cover state of the water column surface by using the threshold LST value 273 K. A water column surface LST less than 273 K is determined to be ice and greater than 273 K is determined to be water. Finally, the change in heat storage capacity was defined as the difference between the sums of each pixel lake water column heat storage capacity in the first and final months. The LAKE and Hostetler models used in this step are introduced in Section 3.3.

$$G={\rho }_w{c}_{pw}\frac{d}{dt}{{\displaystyle \int }}_{i=0}^{i=n}{T}_{wi}dz$$

(5)

Here, $t$ is the current time step (d); ${\rho }_w$ is the density of the water (kg·m⁻³); $c_{pw}$ is the specific heat of the water (J·kg⁻¹·°C⁻¹); and $T_{wi}$ is the water temperature in the water layer i (°C).

Third, for the seven lake groups (S01–S07), the time series of the average lake heat storage changes in each lake group obtained in the second step were used as the training samples to obtain the regression equation related to the meteorological parameters. The general regression equation of each lake group was obtained to calculate the heat storage changes in the lakes without bathymetric data. The regression equation is resented in Equation (6):

$$\mathrm{G}=\mathrm{A}·X+B$$

(6)

where A and B are the regression coefficients, and X is the meteorological parameter that has the best correlation with G.

Fourth, based on the specific estimation equation of each lake group obtained in the previous step, combined with the lake surface meteorological parameter (the main variable of the equation) and the lake surface temperature (used to determine the state of lake water: frozen up, ice water mixed, or open water), the parameters of lake heat storage variations in the water body of the lakes without bathymetric were obtained.

3.3. The One-Dimensional Numerical Lake Model

Two widely used classical one-dimensional lake models, including the Lake and Hostetler models [53,54,55], were used in the study to simulate lake water temperature profiles and lake the ice-coverage ratio of the middle–large-sized lakes sampled in this study.

The Lake model is a comprehensive, turbulent, closed, one-dimensional lake model developed by V.M. Stepanenko et al. [56]. The model uses the finite difference method to solve lake stratification parameters. The Lake model not only simulates the internal thermodynamic processes of lakes, but also includes biogeochemical modules to simulate the exchange counts of O₂, CO₂, and H₂O between lakes and the atmosphere. The kinetic process of the model is expressed by the horizontal mean Reynolds convection-diffusion equation, as shown in Equation (7):

$$c\frac{\partial \overline{f}}{\partial t}=-\frac{c}{A}{{{\displaystyle \int }}^{ }}_{\Gamma A\left(z\right)}f\left(u_h·n\right)dl+\frac{1}{A}\frac{\partial }{\partial z}\left(A{k}_f\frac{\partial \overline{f}}{\partial z}\right)-\frac{1}{A}\frac{\partial A\overline{F_{nz}}}{\partial z}+\frac{1}{A}\frac{dA}{dz}\left(F_{nz,b}\left(z\right)+F_{tz,b}\left(z\right)\right)+R_f\left(\overline{f},\dots \right)$$

(7)

where $f$ is the target parameter velocity component, such as temperature, gas concentration, turbulence kinetic energy, and dissipation rate. Additionally, the first item on the right of Equation (7) is the convection activity caused by inflow and groundwater; the second item is the turbulent diffusion; the third item is a turbulent flux divergence; the fourth item is the total vertical flux at bottom of the lake; and the final item is a source and sink for other values. In this study, because we only needed to consider the temperature diffusion process of lake water, Equation (7) can be rewritten as Equation (8). When only considering the thermodynamic process, $f$ is replaced by water temperature $T$. $z$ is the depth. The lake surface depth is defined as $z=0$, and the vertical downward value is positive. $h=h\left(t\right)$ is the lake depth at time t. $\xi =z/h$is defined as the coordinate for convenience.

$$c_w\rho \frac{\partial T}{\partial t}=\frac{c}{h^2}\frac{\partial }{\partial \xi }\left(\lambda \frac{\partial T}{\partial \xi }\right)+c_w\rho \frac{dh}{dt}\frac{\xi }{h}\frac{\partial T}{\partial \xi }-c_w\rho \frac{1}{h}\frac{d{h}_0}{dt}\frac{\partial T}{\partial \xi }-\frac{1}{h}\frac{\partial S}{\partial \xi }+M$$

(8)

Here, $c_w$ is the specific heat capacity of the water; $\rho$ is the density of the water; $\lambda$ is the eddy diffusivity; $d{h}_0/dt$ is the water balance on the lake surface; $S$ is solar radiation flux; penetration to depth is $z$; and $M$ is buoyancy mixing behavior (convection).

The Hostetler model is a typical vorticity diffusion model that divides the entire lake into several water layers to obtain the lake water temperature profile by solving the vertical thermal diffusion equation [57,58]. The basic thermal mixing processes occurring in lake water include wind-driven eddy turbulence, parameterized thermal conductivity based on Henderson-sellers, buoyant convection, and molecular thermal diffusion processes. Important factors, such as heat, water, wind force, and radiation flux, were also considered in the model. The governing equation of the Hostetler model is presented in Equation (9). The Hostetler model simulates the thermal diffusion process of lake water by inputting driving and lake water depth data, and it can directly output the water ice phenology change process. Based on the simulated lake ice-coverage ratio (ICR) and the ice thickness, we classified the lake water status as ice-free (ICR ≤ 10%), ice and water mixing (10% ≤ ICR ≤ 90%), and freeze-up (ICR ≥ 90%) periods.

$$\frac{\partial T_{z,t}}{\partial t}=\frac{1}{A_{z,t}}\frac{\partial }{\partial z}\left[A_{z,t}\left({\kappa }_m+{\kappa }_{E,z,t}\right)\frac{\partial T_{z,t}}{\partial z}\right]+\frac{1}{A_{z,t}}\frac{1}{c_{w,z,t}}\frac{\partial \left(K_{z,t}{A}_{z,t}\right)}{\partial z}$$

(9)

Here, $T_{z,t}$ is the water temperature at depth z and time t; $A_{z,t}$ is the water area at depth z and time t; ${\kappa }_m$ is the molecular diffusivity (1.39 × 10⁻⁷ m²s⁻¹); ${\kappa }_{E,z,t}$ is the eddy diffusivities (m²s⁻¹); $c_{w,z,t}$ is the specific heat capacity of lake water at depth z and time t; and $K_{z,t}$ is the transmitted short-wave energy at depth z.

3.4. Empirical Formula Method for Lake Ice Sublimation

Due to the high altitude and low air temperature and pressure prevalent in the inland area of the TP, the lake water bodies in this region present longer ice phenology of freeze onset (approximately 3 months) and ice–water mixing periods (include ice freeze-up and break-up, approximately 2 months) than lakes in low-altitude regions [31]. Therefore, ignoring the sublimation contribution of lake ice may lead to some errors in analyzing the spatial and temporal variations in lake evaporation levels in the inland areas of the TP. Sublimation can occur as long as the relative humidity level is lower than 100% and the temperatures of the ice–vapor interface is less than 273 K [59]. In this paper, we elected an empirical method to estimate the sublimation rate; the detailed formula and calculation process is follows:

$$S_{rate}=3.96\frac{u_2\left(e_s-e_a\right)}{\left(273.16+T\right){\rho }_{ice}}$$

(10)

where $S_{rate}$ is the sublimation rate (cm), T is the air temperature (°C), and ${\rho }_{ice}$ is the density of ice (0.92 g cm⁻³).

3.5. Attribution Analysis

The “detrending” method [60,61] was used to analyze the contributions of change in related variables (such as meteorological variables and surface reservoir area) to the trend of evaporation volume/rate. First, the base case of the evaporation volume/rate ($E{T}_{base}$) was calculated using all the detrended variables. Then, the evaporation volume/rate affected by variable X ($E{T}_X$) was calculated using the original variable X, and the rest of the detrended variables, for example, the evaporation volume/rate affected by Ta, were calculated using the original Ta and the rest of the detrended variables. The contribution of variable X to the evaporation volume/rate was calculated as follows:

$$Co{n}_X=\frac{Tren{d}_{E{T}_X}-Tren{d}_{E{T}_{base}}}{Tren{d}_{E{T}_{original}}}\times 100\%$$

(11)

where $Co{n}_X$ is the contribution of variable X to the evaporation volume/rate, and $E{T}_{original}$ was calculated using all the original variables. All trends and significance indexes used to achieve process parameters and evaporation results in this study were calculated using the nonparametric Mann–Kendall test.

4. Result

4.1. The Lake Heat Storage Change

For the lakes without available bathymetric data, it was difficult to estimate lake heat storage changes using the Penman formula. Here, according to the method presented in Section 3.2, we used the monthly heat storage change value results of the lakes with bathymetric data obtained from 2002 to 2018 as the model training sample data to calculate the correlation values between the heat storage changes in lakes and meteorological parameters, and we created a grid diagram of the correlation coefficients among all the parameters (Figure 6).

As can be observed from Figure 6, since the heat exchange between the lake surface and atmosphere was almost negligible during the lake surface freeze-up period, the lake heat storage change value presented a low correlation for almost all the meteorological parameters. However, there was a good correlation between the lake heat storage change and lake surface net radiation in the whole, ice and water mixing, and ice-free periods, and the correlation coefficient R values were 0.68, 0.89, and 0.79, respectively. Therefore, it was feasible to establish the one-dimensional linear regression equation of lake heat storage changes based on the net radiation of lake surface. The one-dimensional linear regression equation of each lake group is shown in Equation (6) X in Section 3.2. The parameter list obtained through the regression analysis is shown in

Table 4. Summary of the locally calibrated coefficients of the linear regression equation for the heat storage change and the statistical agreement between linear regression and vertical integrate methods.

Lake Groups	No. Lakes with Measured Depth/No. Lakes without Measured Depth	No. of Data Pairs	G = a × Rn + b		Statistical Agreement
			a	b (w/m2)	R2	RMSE (w/m2)	NSE
S01	2/1	196	0.97	−77.57	0.66	28.78	0.45
S02	9/14	673	1.00	−80.66	0.59	26.88	0.30
S03	8/22	543	1.03	−89.78	0.65	27.13	0.45
S04	14/11	1275	0.85	−84.75	0.51	29.13	0.30
S05	8/8	968	1.02	−107.84	0.62	31.25	0.40
S06	19/8	2429	1.15	−117.80	0.71	28.89	0.55
S07	9/1	1178	1.09	−107.28	0.77	23.43	0.52

. The lake heat storage change value that was estimated by the water temperature vertical integrate method was used as verification data. Compared with the linear regression results, it can be observed that linear regression method can rapidly estimate the heat storage changes in lakes without in situ bathymetric data and presents good accuracy, with R² values ranging from 0.51 to 0.77 and root mean square error of values 23.43 to 31.25 w/m².

The change in monthly heat storage values between different lake groups from 2002 to 2018 showed that all lake groups had the same seasonal characteristics (Figure 7). In spring, because the air temperature began to rise and lake surface ice begins to break up, the cold lake water begins to absorb heat, which causes the lake heat storage change to positively increase and reach the maximum heat absorption level in summer. In the early autumn, due to the decrease in the air temperature, the heat absorption rate of lake water decreases; the heat absorption process of lake water remains positive for a certain period of time. In late autumn, the air temperature and sun radiation levels are below a certain threshold value, and the temperature of the lake water is higher than that of the atmosphere. The lake starts to release heat from lake water surface into atmosphere until the lake water and atmosphere temperature become relatively balanced. In winter, as the lake surface freezes, the heat exchange between the lake surface and atmosphere reduces until the lake water surface is completely frozen and the heat exchange almost ceases. The heat storage change value is close to 0 w/m², and a complete cycle is formed. Although the general pattern is similar, there are still some differences among different lake groups. Lake groups S05, S06, and S07 have a shorter glaciation period due to their lower latitude, and the absorption and release peak interval of ±80 w/m² is significantly higher than that of lake groups S02, S03, and S04, which is ±50 w/m².

4.2. Spatial–Temporal Variations in the Lake Evaporation Rate

For the 134 middle–large-sized lakes in the inland are of the TP as a whole, the lake evaporation rate showed an insignificant upward trend of about 0.10 mm/year from 2002 to 2018. The maximum evaporation rate was 976.41 mm/year in 2010, and the minimum was 859.12 mm/year in 2002. The average annual evaporation rate in 2018 was 933.83 mm/year. For the seven lake groups, S05 had the largest annual mean evaporation rate value of 1148.05 mm/year in 2018. The minimum average annual evaporation rate in the S03 region was 718.73 mm/year. Spatially, the evaporation rate increased from north to south ranging from 594.17 to 1309.27 mm/year, with a median value of 918.24 mm/year.

It can be observed from the lake evaporation rate spatial distribution diagram that there is an obvious north–south difference in lake evaporation rate values in the inland area of the TP (Figure 8). When using the 34° N latitude as the dividing line, the 134 lakes in the inland area of the TP were divided into north and south lake groups for which the time-series difference changes in lake evaporation rate are shown in Figure 9. There were 59 middle–large-sized lakes located in the north and 75 lakes in the south. From 2002 to 2018, the annual evaporation rate in the southern lakes showed an upward trend of about 0.62 mm/year, while that in the northern lakes showed a downward trend of about −0.80 mm/year. From 2002 to 2018, the annual evaporation rate for the northern lakes was 750.09 mm/year; the annual evaporation rate for the southern lakes was 1064.50 mm/year; and the difference between the annual evaporation rates of the northern and southern lakes was 314.41 mm/year. The maximum, minimum, and median differences between the evaporation rates for south and north middle–large-sized lakes in the inland area of the TP were 350.49, 268.65, and 308.52 mm/year, which appeared, respectively, in 2009, 2013, and 2012.

4.3. Spatial–Temporal Variations in the Lake Evaporation Volume

For the 134 middle–large-sized lakes located in the inland area of the TP, the total average annual evaporation volume from 2002 to 2018 was 25.02 km³. The maximum evaporation volume of 27.53 km³ occurred in 2016, and the minimum evaporation volume of 20.27 km³ occurred in 2002. In 2018, the total evaporation volume for middle–large-sized lakes was 27.38 km³. Among the lake groups, S06 presented the maximum evaporation volume with 9.20 km³ and accounted for 34.85% of the total lake evaporation volume in the same year; S01 presented the minimum evaporation volume with 1.13 km³ accounting for 4.18%. For the 134 middle–large-sized lakes in the inland area of the TP, their annual average evaporation volumes had large spatial heterogeneity ranging from 0.03 to 2.61 km³, with a median value of 0.09 km³ (Figure 10).

In order to explore the contribution of the lakes’ evaporation volumes from different-sized lakes, we divided 134 middle–large-sized lakes into three levels based on the water area, including level 1: 50 km² ≤ lake’s area < 100 km² (61 lakes), level 2: 100 km² ≤ lake’s area < 500 km² (63 lakes), and level 3: 500 km² ≤ lake’s area (10 lakes), and performed a statistical numerical analysis of the evaporation volume. Among them, from 2002 to 2018, the annual average value of total evaporation volume of all level-1 lakes was 3.53 km³; the maximum value was 4.14 km³; and the minimum value was 2.47 km³. The mean value of total evaporation of level-2 lakes was 11.13 km³; the maximum value was 12.26 km³; and the minimum value was 9.04 km³. The annual mean, maximum, and minimum values of total evaporation volume of level-3 lakes were 10.37, 11.13, and 8.76 km³, respectively. The three levels’ maximum and minimum values of annual total evaporation volume all occurred in 2016 and 2002, respectively (Figure 11). We also performed a statistical analysis of the different levels’ lake evaporation volume contribution ratios to the total middle–large-sized lakes’ evaporation volumes in the inland area of the TP. The results show that, from 2002 to 2018, level-1, –2, and -3 lakes contributed 14.04%, 44.46%, and 41.50% of the total annual average evaporation volume levels, respectively. At the same time, the three levels of lakes showed different trends. The contribution of level-1 lakes showed an upward trend from 12.18% in 2002 to 15.07% in 2018, and the contribution of level-3 lakes showed a downward trend from 43.21% to 40.22% in 2002 to 2018. There was no obvious change in the level-2 lakes, which was maintained at 44.50% (Figure 11).

4.4. Magnitudes and Trends of Lake Evaporation

The spatio-temporal variations in lake evaporation rates in the seven lake groups are shown in Figure 12. S02 and S07 showed an obvious upward trend, with increased values of 5.81 and 5.95 mm/year, respectively. S01, S03, and S04 showed a decreasing evaporation rate trend, and S01 showed the fastest decreasing rate of −8.10 mm/year, while S04 showed the slowest decreasing rate of −3.16 mm/year. Meanwhile, although the annual evaporation rates of S05 and S06 showed a decreasing trend, the trends change were not obvious, which were −0.36 and −0.67 mm/year, respectively. Among the 134 middle–large-sized lakes, 71 showed a decreasing trend in the evaporation rate, while 63 lakes showed an increasing trend. Although the number of lakes with decreasing evaporation rates was higher than that with increasing evaporation rates, the increasing trend was more significant than the decreasing trend in the intensity of change. The lake with the greatest decreasing trend in evaporation rate was Lake Ayakkum Co (−18.28 mm/year), and the lake with the greatest increasing trend was Lake Dong Co (10.71 mm/year).

Due to the increasing trends in lake area and lake surface evaporation rate from 2002 to 2018, the total evaporation volume of middle–large-sized lakes in the inland area of the TP also showed an obvious upward trend, with a value of about 0.35 km³/year. For the seven lake groups, the annual total volume of lake evaporation rates showed varying upward-trend levels, among which the maximum upward trend was 0.10 km³/year in S03 and the minimum upward trend was 0.002 km³/year in S05. A total of 108 lakes showed an upward trend in lake evaporation volume, and 26 lakes showed a downward trend (Figure 13). Lake Zhuonai showed the greatest decreasing trend in evaporation volume of −0.0054 km³/year, and Lake Serling Co showed the greatest increasing trend of 0.0362 km³/year.

4.5. Attribution Analysis

For each middle–large-sized lake in the inland area of the TP, the “Attribution analysis” method was used to analyze the contributions of the changes occurring in all parameters, which participated in the evaporation rate estimating method, such as lake heat storage (G), net radiation of lake surface (Rn), air pressure (Pres), air specific humidity (Shum), air temperature (Temp), vapor pressure deficit (VPD), and air wind speed (Wind), to the trend of evaporation rate. The result showed that the median values of contributions of the seven variables (G, Rn, Pres, Shum, Temp, VPD, and Wind) to the evaporation rate trend were −21.05%, 60.44%, −0.82%, 0%, 90.18%, 0%, and 0%, respectively (Figure 14a). Among the 134 lakes studied, the numbers of lakes with G, Rn, Shum, and Temp as main contribution variables were 9, 34, 1, and 90, respectively (Figure 14b). Therefore, air temperature is the main dominant factor affecting evaporation rate, followed by the net radiation of the lake’s surface and lake heat storage value.

In addition to the variables related to the evaporation rate, lake area was added to the “Attribution analysis” method to analyze the contributions of change in each variable to the trend of evaporation volume for each lake. The median values of contributions of the eight variables (G, Rn, Pres, Shum, Temp, VPD, Wind, and lake surface area (Area)) to the trend of evaporation volume in the 134 lakes were 2.83%, −4.05%, −0.13%, 0%, 12.20%, 0%, 0%, and 72.62%, respectively (Figure 15a). Among the 134 lakes, the numbers of lakes with G, Rn, Temp, and Area as the main contribution variables were 3, 10, 39, and 81, respectively (Figure 15b). Thus, the change in lake surface area dominated the trend of evaporation volume.

5. Discussion

5.1. Comparison between Evaporation Rate in This Paper and Previous Studies

To date, few studies focus on the evaporation rates of the middle–large-sized lake groups in the TP. The most detailed analysis of lake evaporation in the TP were published by Wang et al. (2021), whose study covers 75 different-sized lakes’ evaporation rates over the whole TP based on the Bowen ratio method combined with CMFD datasets [31], and among them, there were 57 middle–large-sized lakes that overlapped with our studied 134 lakes. The result determined by Wang et al. (2021) an be used as a reference: the average annual evaporation rates compared to our results and the reference for the overlapped 57 middle–large-sized lakes show that the correlation coefficient (R), root mean square error (RMSE), percent bias (PBIAS), and Nash–Sutcliffe efficiency (NSE) are 0.95, 61.00 mm, 1.29%, and 0.90, respectively (Figure 16). The good agreement indicates that our method is reliable to calculate middle–large-sized lakes’ evaporation rates in the inland of the TP. On the one hand, the good accuracy between the two methods indicates that it is feasible to estimate the middle–large-sized lakes’ heat storage change values by using the regression equation. On the other hand, the Penman method can better estimate lake evaporation levels without in situ bathymetric data and thus obtain evaporation spatio-temporal variation characteristics.

However, compared with the reference results, 2.71% and 0.38% underestimations are evident our results for the lake groups S01–S04 and S05–S07. These differences were mainly due to the following reasons: (1) different locations of the lakes produce different evaporation rates due to water depth differences. In order to quantitatively analyze the evaporation rate influenced by water heat storage value in Lake Serling Co, a previous study found that the deeper the lake depth, the lower the evaporation rate of the lake’s surface [62]. This produces some errors when calculating the evaporation rate by only using the center point water depth or average lake depth values. (2) The duration of the ice-free season regarding the water column surface influences the evaporation rate directly. Reference data and our study showed that the ice-free period was about 152.4 and 124.3 days in average in the lake groups S01–S04. This was probably the main reason for the evaporation underestimation that occurred in lake groups S01–S04. (3) The difference in correction methods of the basic data can also cause some differences between the reference results and our study. By the reference data, we determined the average land surface temperature and air temperature values as 283.5 K and 279.1 K in lake groups S01–S04, and 285.6 K and 282.5 K in lake groups S05–S07; however, they were lower in lake groups S01–S07 in our study. This can explain the difference between reference result and evaporation rate values in our study.

5.2. Uncertainty

Some uncertainties exist in this study due to the data availability and model parameters design. First, because of the extreme natural environment of the TP, the maintenance of 11 TORP observation station instruments and the in situ bathymetric survey of 68 lakes were difficult to conduct, resulting in some errors in the measured data. It was necessary to conduct the timely meteorological instrument calibration and maintenance of the observation stations and supplementation of lake surveys. At the same time, the remote sensing data errors were inevitable, although CMFD and MODIS data which were used in this paper were corrected. So, in the future we will find more suitable data to replace them. Secondly, the application prerequisite of the one-dimensional lake model was that the lake water body was fully mixed, and the difference of water bodies in the vertical direction was more significant than that it was in the horizontal direction. However, the horizonal differences of the water bodies still existed due to strong wind in the inland area of the TP. It was therefore necessary to develop 3D lake models considering the in situ lake bathymetrical data to reduce this potential uncertainty. Finally, compared with the previous studies, the sublimation model of lakes during the freezing-up and ice-water mixing periods was considered in this paper. However, the sublimation model used in this paper is still a typical semi-empirical model used in the research; with its lack of detailed physical significance, it causes some uncertainties regarding lake evaporation rate and volume data.

6. Conclusions

Combined with the detailed in situ bathymetric data of the lakes, the difference in lake heat storage changes was considered in the Penman formula to estimate lake evaporation occurrence in a large number of middle–large-sized lakes located in the inland area of the TP. Our results show the following outcomes:

(1)

The lake heat storage changes in lakes without available bathymetric data may be estimated by the regional regression model, which is established between the lake heat storage changes in lakes with available water depth and lake surface net radiation values.

(2)

There is a high degree of accordance of evaporation rates estimated by both the Bowen ratio and Penman methods in 57 middle–large-sized lakes. The correlation coefficient (R) is 0.95; the root mean square error (RMSE) is 61 mm; the percentage bias (PBIAS) is 1.29%; and the Nash–Sutcliffe efficiency (NSE) is 0.90. The good accuracy of the evaporation results make it abundantly clear that the heat storage change regression equations are credible.

(3)

From 2002 to 2018, the average annual evaporation rates of the 134 middle–large-sized lakes in the inland of the TP show an insignificant upward trend at about 0.10 mm/year but presented obvious spatial differences. The annual evaporation rate in the southern lakes was 1064.50 mm/year with an upward trend of about 0.62 mm/year, while that in the northern lakes was 750.09 mm/year with a downward trend of about −0.80 mm/year. Temperature is the main factor affecting the evaporation rate of 134 middle–large-sized lakes.

(4)

For the 134 middle–large-sized lakes we studied, the average annual evaporation volume from 2002 to 2018 was 25.02 km³ with an obvious upward trend of about 0.35 km³/year. The lake evaporation results did not show obvious spatial distribution differences; however, the total evaporation volume of different levels lakes showed obvious results. The results show that, from 2002 to 2018, the middle–large-sized lakes annual average evaporation volume contribution ratios were 14.04%, 44.46%, and 41.50% from 61 lakes with an area of 50–100 km², 63 lakes with an area of 100–500 km², and 10 lakes with an area greater than 500 km², respectively. Lake area is the key influence factor of the total lake evaporation volume occurring in the inland area of the TP.