Simulation of Hydrologic Change of Linggo Co during 1979–2012 Using Hydrologic and Isotopic Mass Balance Model

Abstract

The Tibetan Plateau (TP) and surrounding areas contain the largest number of glaciers outside the polar regions. The region affects downstream water supply and food security, thereby directly influencing one-third of the world’s population. The lakes in the central TP expanded rapidly in recent decades, which has attracted growing attention. Glacier meltwater was considered as a major component in the water balance of TP lakes, although few studies quantified its contribution. Stable isotope analysis is a powerful tool to trace hydrologic circulation, while its interpretation in paleoclimate records has been controversial. To bridge the gap between hydrologic and paleoclimatic studies, we performed a hydrologic and isotopic mass balance model to simulate the lake level change of Linggo Co in the central TP. The model was forced by the meteorological data, calibrated through observed lake level changes, and validated by oxygen isotope compositions (δ¹⁸O) of lake water. Our results indicated that glacier meltwater contributed 73.94% of the inflow water to Linggo Co before 1993 but decreased thereafter. Increasing glacier meltwater together with positive water balance (precipitation/evaporation) in the catchment contributed to the rapid expansion of Linggo Co after the mid-1990s. Lake water δ¹⁸O in Linggo Co was more sensitive to changes in the precipitation amount and precipitation δ¹⁸O than temperature. Our findings could shed light on the usage of δ¹⁸O proxy in future paleoclimate research on the TP.

1. Introduction

The Tibetan Plateau (TP) and surrounding areas contain the largest number of glaciers outside the polar regions [1]. The water resources on the TP affect downstream water supply and food security, thereby directly influencing one-third of the world’s population [2,3,4]. In recent decades, the TP has undergone significant warming, along with remarkable changes in the hydrosphere [2,5]. Major lakes in the central TP have expanded since the 1990s, whereas many lakes in the southern TP have shrunk [5,6,7,8,9]. The reasons for lake variations in the TP have been widely discussed [1,2,7]. Among various factors (e.g., glacier melt, precipitation, and potential evaporation), glacier melt was the most frequently invoked contributor [1,6,9,10]. However, only a few studies quantified its contribution to the lake water budget [7,11,12], which is critical for understanding the interactions between cryosphere and hydrosphere and the response of lake dynamics to climate change.

Stable isotope analysis has long been considered a powerful tool to reconstruct the paleoclimate [13,14,15,16,17] and hydrological change [18,19,20,21,22,23]. For a lake system, oxygen isotope compositions (δ¹⁸O) of lacustrine carbonates were controlled by isotopic compositions of lake water, which in turn responded to hydrological changes driven by climate change [18,24]. Increasing temperature would enhance evaporation on the lake surface and catchment, leading to ¹⁸O enrichment of lake water [18]. Increasing precipitation amount would alter the isotopic composition of inflowing river, thus the δ¹⁸O of lake water [18]. Moreover, δ¹⁸O of precipitation was associated with temperature and precipitation amount due to the so-called “temperature effect” and “amount effect” [18]. The interpretation of δ¹⁸O record on the TP is even more complex since many lakes on the TP are located in a close basin with strong evaporation, and received multi-source water (glacier meltwater, summer monsoon, or mid-latitude westerlies induced moisture) with distinguished isotopic compositions [25,26]. Quantitatively assessing the influence of these factors, such as temperature, precipitation amount, and isotopic compositions of precipitation on the isotopic composition of lake water, can shed light on the interpretation of stable isotope proxies in paleoclimate reconstruction on the TP.

Linggo Co (Figure 1) is a close basin lake in the central TP [27]. It is fed primarily by meltwater from the Purogangri ice sheet [28]. Remote sensing-based investigations have suggested large lake level variations in recent decades [10]. δ¹⁸O records of ostracod shells have revealed a sharp 4‰ shift at the beginning of the Holocene [29]. For a better understanding of this isotopic variation, we performed a hydrologic and isotope mass balance model of Linggo Co. The model was forced by meteorological data from 1979 to 2012, calibrated by the lake area change, and validated by lake volume and lake water δ¹⁸O measurements. Sensitivity analyses were conducted on temperature, precipitation amount, and precipitation δ¹⁸O. The main objectives of this study are to (1) quantify the contribution of glacier meltwater to the expansion of Linggo Co and (2) access the effects of temperature, precipitation amount, and precipitation δ¹⁸O on lake water δ¹⁸O.

2. Materials and Methods

2.1. Geological Settings

Linggo Co (88°31′–88°39′ E, 33°48′–33°54′ N, 5059 m a.s.l.) is located within a 1830 km² closed basin in the central TP (Figure 1). The Puruogangri ice sheet to the east of the watershed has an extent of ~422.6 km², of which ~135 km² lies in the catchment of Linggo Co [27]. Small seasonal rivers, which are supplied by the meltwater from Puruogangri ice sheet in the summer, flow over the Gobi zone in the piedmont and drain into Linggo Co [28]. The lake area is characterized by a sub-frigid semi-arid climate. Based on the observations from October 2011 to September 2013 [30], the mean annual temperature is about −4 °C, with a high monthly temperature of about 7 °C in July and a low of about −15 °C in January. The annual precipitation, which mainly falls from June to September, is ~333 mm. Actual evapotranspiration is close to the annual precipitation, while the mean annual potential evaporation is ~1000 mm [30,31].

2.2. Data Source and Structure of the Model

The meteorological data of the model mainly came from the China Meteorological Forcing Dataset (CMFD, http://westdc.westgis.ac.cn/ (accessed on 16 May 2015)), which was produced by merging a variety of data sources (including China Meteorological Administration (CMA) station data, TRMM satellite precipitation data (3B42), and the GLDAS precipitation data) with the spatial and temporal resolutions of 3 h and 0.1°*0.1°, respectively [32,33]. In order to assess the data quality, the CMFD data were compared to the High Asia Refined Reanalysis Dataset (HAR, http:// www.klima-ds.tu-berlin.de/har/ (accessed on 14 June 2015)) [34] and CMA station observed data (temperature, wind speed, relative humidity, and precipitation) from five stations around our study area (Shenza, Gaize, Naqu, Bange, and Amduo; Supplementary Figures S1–S4). Considering the overestimation of precipitation in the monsoon season (from May to October) of the CMFD data after 2000 (Supplementary Figure S1), CMFD data were combined with HAR data from May to October during 2000–2012 (Supplementary Figure S5).

Glacier albedo was extracted from HAR data [34]. Glacier meltwater δ¹⁸O was based on the average value of Puruogangri ice core (−12.5‰, [35]). Sunshine hours data came from CMA data [36]. The lake area was extracted using remote sensing and GIS techniques [10]. The bathymetry map was measured in the field in 2011 and digitized using ArcGIS software to extract the area at the surface and bottom. Monthly precipitation δ¹⁸O was extracted using the Online Isotopes in Precipitation Calculator (Version 2.2, [37]). One lake water sample was collected near the northern shore of the lake (red dot in Figure 1) on May 25th, 2015, by Dr. Yanbin Lei from Institute of Tibetan Plateau Research, Chinese Academy of Sciences (ITPCAS), and oxygen isotope was analyzed in the Key Laboratory of Tibetan Environment Changes and Land Surface Processes, ITPCAS (personal communication).

2.3. Model Structure

The hydrologic and isotope mass balance model of Linggo Co is simplified after Steinman [22,23]. The model is described using the following equations:

$$\frac{d{V}_L}{dt}={\sum }^{ }\mathrm{I}-{\sum }^{ }\mathrm{O}$$

(1)

$$\frac{d\left(V_L\ast \delta L\right)}{dt} ={\sum }^{ }\mathrm{I}\times \mathsf{\delta }\mathrm{I}-{\sum }^{ }\mathrm{O}\times \mathsf{\delta }\mathrm{O}$$

(2)

where V_L is the reservoir volume, ∑I and ∑O are the total inflows and outflows from the reservoir, and δ is the isotopic composition of the inflows and outflows. The hydrologic model is defined by three separate equations, each corresponding to a different water reservoir: the glacier, catchment, and lake (Figure 2).

2.3.1. Glacier Reservoir

There are three processes that occur at the glacier surface: melting, accumulation, and sublimation. The glacier mass balance is based on energy balance methods [38]. The energy at the surface–air interface is expressed as:

Q_N + Q_H + Q_L + Q_G + Q_R + Q_M = 0

(3)

where Q_N is the net radiation, Q_H is the sensible heat flux, Q_L is the latent heat flux (Q_H and Q_L are referred to as turbulent heat fluxes), Q_G is the ground heat flux, Q_R is the sensible heat flux supplied by rain, and Q_M is the energy consumed by melt. The Q_N represents the majority of the energy based on the observation in the TP (Zhadang glacier: 82% [39]; Qiyi glacier: 90% [40]). Therefore, we use Q_N as the melt energy in this work.

$${\mathrm{Q}}_{\mathrm{N}} =\mathrm{G} \times (1-\mathsf{\alpha })+\mathrm{L}\downarrow +\mathrm{L}\uparrow$$

(4)

where G is the global radiation, α is the surface albedo, L$\downarrow$ is the downward longwave radiation, and L$\uparrow$ is the upward longwave radiation which can be calculated as:

$$\mathrm{L}\uparrow = \epsilon \sigma T^4$$

(5)

where ε is the emissivity; σ = 5.67 × 10⁻⁸ Wm⁻² K⁻⁴ which is known as Stefan–Boltzmann constant [38]. For melting surface, the L$\uparrow$ = 316W*m⁻². If the melt energy (Q_N) is less than zero, all the precipitation/snow is accumulated, while if the energy is greater than zero, the melt rate, M, is then computed from the available energy by:

$$\mathrm{M}=\frac{Q_M}{\rho W\times Lf}$$

(6)

where ρ_W denotes the density of water; L_f = 334,000 J*kg⁻¹ stands for the latent heat of fusion. The sublimation is also calculated by the above equation. The difference is the energy for sublimation, which is 2,848,000 J*kg⁻¹.

2.3.2. Catchment Reservoir

The catchment is described as a two-layer soil model. Precipitation and evapotranspiration occur within the top soil. The deep soil receives glacier meltwater, which experiences a brief exposure time before evaporating and forming runoff. In the rainy season, if top soil receives more water than its maximum available water capacity, water in top soil may infiltrate into deep soil.

The actual evapotranspiration is calculated following Morton’s CRAE method [41]. This method is based on the hypothesis that, for large homogeneous areas where there is little advective heat and moisture, there is a complementary relationship (CR) between actual evaporation and potential evaporation:

ET_act = 2ET_wet − ET_pot

(7)

where T_act is the actual areal or regional evapotranspiration from an area large enough that the heat and vapor fluxes are controlled by the evaporation power of lower atmosphere and unaffected by upwind transitions. ET_wet is the potential evapotranspiration that would occur under steady-state meteorological conditions in which the soil/plant surface is saturated with abundant water. ET_pot is the (point) potential evapotranspiration for an area so small that the heat and water vapor fluxes have no effect on the overpassing air. The CR relationship has been confirmed in the field observation at Shanghu area that is 60 km south of Linggo Co [42]. The details of the calculation process can be found in the study by McMahon et al. (2013) [43].

2.3.3. Lake Reservoir

Evaporation, precipitation, runoff, and outseepage occurred in the lake reservoir. The evaporation model applied in this paper is the simplified version of the modified Penman equation with elevation calibration proposed by Valiantzas (2006) [44]:

$$\begin{gathered} \mathrm{E}=\left[0.051\times \left(1-\mathrm{ALB}\right)\times R_S\times {\left(T_a+9.5\right)}^{\frac{1}{2}}-2.4\times \left(R_S/R_a\right)2+0.052\times \left(T_a+20\right)\times (1 \\ -\mathrm{RH}/100)\times \left(a_u-0.38+0.54\times \mathrm{WS}\right)+0.00012\times \mathrm{ELE}\right]\times 30 \end{gathered}$$

(8)

where E is the evaporation in mm per month, ALB is the albedo of the lake (0.08), T_a is the average daily temperature in °C, R_S is the solar radiation in MJ m⁻²d⁻¹, and R_a is the extraterrestrial radiation in MJ m⁻²d⁻¹. RH is the average daily relative humidity expressed as a percentage, a_u is the Penman wind function constant, and WS is the average daily wind speed in ms⁻¹. ELE is the elevation of the lake. The evaporation is set to zero if T_a < 0. R_S and R_a are calculated in the following equations:

R_S = R_a × (0.5 + 0.25 × n/N)

(9)

N ≈ 4∅ sin(0.53i − 1.65) + 12

(10)

R_a ≈ 3N sin (0.131N − 0.95∅) for |∅| > 23.5π/180

Or R_a ≈ 118N^0.2 sin (0.131N − 0.2∅) for |∅| < 23.5π/180

(11)

where n is the measured bright sunshine hours per day (h) and N is the maximum possible duration of daylight (h). ∅ is the latitude of the study site (radians), positive for the Northern Hemisphere and negative for the Southern Hemisphere. i is the rank of the month (the first month is January).

The outseepage function is based on the exponential relationship between lake radius and hydraulic conductivity values (K) described by Genereux et al. (2001) [45]:

K = A × e^bR

(12)

where R is the radius of an ideal round lake. A and b are parameters that vary with different settings of lakes. We simplified the Linggo Co hypsographic profile to an upside-down frustum of a right circular cone (Supplementary Figure S6). Thus, the total hydraulic conductivity (K_total) can be expressed as:

$${\mathrm{K}}_{\mathrm{total}}={\int }_{\mathrm{R}0}^{\mathrm{R}1}A\times {\mathrm{e}}^{\mathrm{bx}}\times dx$$

(13)

where R₀ and R₁ are the radii at the bottom and surface of Linggo Co (Supplementary Figure S6). A and b are coupled parameters that needed to be adjusted in the model calibration.

2.3.4. Isotope Mass-Balance Equations

Oxygen isotope value, expressed in standard delta (δ) notation as the per mil (‰) relative to the standard Vienna Standard Mean Ocean Water (VSMOW), is calculated for each reservoir at each time step and is multiplied by the corresponding hydrologic flux to determine the isotope mass balance of each theoretical water mass.

The model of Craig et al. (1965) [46] is used to calculate the isotopic composition of water evaporating from the lake surface:

$${\mathsf{\delta }}_{\mathrm{E}}=\frac{{\mathsf{\alpha }}^{*}{\mathsf{\delta }}_{\mathrm{L}}-{\mathrm{h}}_{\mathrm{n}}{\mathsf{\delta }}_{\mathrm{A}}-{\mathsf{\epsilon }}_{\mathrm{tot}}}{1-{\mathrm{h}}_{\mathrm{n}}+0.001{\mathsf{\epsilon }}_{\mathrm{k}} }$$

(14)

where α* is the reciprocal of the equilibrium isotopic fractionation factor, ${\delta }_L$ is the isotopic composition of the lake water, $h_n$ is the ambient humidity normalized to lake water temperature, ${\delta }_A$ is the isotopic composition of atmospheric moisture, and ${\epsilon }_{tot}$ is the total isotopic separation, which equals to (${\epsilon }_{eq}$+${\epsilon }_k$). ${\epsilon }_{eq}$ is the equilibrium isotopic separation and ${\epsilon }_k$ is the kinetic isotopic separation.

The saturation vapor pressure of the overlying air ($e_{s-a}$) and the saturation vapor pressure at the surface water temperature ($e_{s-w}$) in millibars are used to determine the normalized relative humidity ($h_n$).

$${\mathrm{h}}_{\mathrm{n}} =\mathrm{RH} \times e_{\mathrm{s}-\mathrm{a}}/e_{\mathrm{s}-\mathrm{w}}$$

(15)

$${\mathrm{e}}_{\mathrm{s}-\mathrm{a}}\mathrm{or} {\mathrm{e}}_{\mathrm{s}-\mathrm{w}}=6.108\times exp\left(\frac{17.27\times T}{T+237.7}\right)$$

(16)

Atmospheric moisture (${\delta }_A$) is assumed to be in isotopic equilibrium with precipitation [47]:

$${\mathsf{\delta }}_{\mathrm{A}} ={\mathsf{\delta }}_{\mathrm{A}} -{\mathsf{\epsilon }}_{\mathrm{eq}}$$

(17)

The equilibrium isotopic fractionation factor (α) and the reciprocal of the equilibrium isotopic fractionation factor (α*) for oxygen are calculated as [48]:

$$\mathrm{In}\mathsf{\alpha }=0.35041\times \left(\frac{{10}^6}{T_{\mathrm{w}}^3}\right)-1.6664\times \left(\frac{{10}^3}{T_{\mathrm{w}}^2}\right)+6.7123\times \left(\frac{1}{T_{\mathrm{w}}}\right)-7.685\times {10}^{-3}$$

(18)

where T_w is the temperature (degrees K) of the lake surface water. The lake surface temperature is assumed to be zero if the atmospheric T_a < 0; the lake surface temperature is 0.95× T_a before June, while it is equal to T_a after July. The per mil equilibrium separation (${\epsilon }_{eq}$) of oxygen is calculated as:

$${\mathsf{\epsilon }}_{\mathrm{eq}}=1000\times \left(1-{\mathsf{\alpha }}^{*}\right)$$

(19)

Kinetic fractionation (${\epsilon }_k$) is controlled by molecular diffusion and moisture deficit (1 − $h_n$) over the lake surface:

$${\mathsf{\epsilon }}_{\mathrm{k}}=\mathrm{C}\times \left(1-{\mathrm{h}}_{\mathrm{n}}\right)$$

(20)

where C is the experimentally derived isotopic separation value of 14.3‰ for oxygen that is applied.

2.3.5. Model Calibration and Sensitivity Analyses

All model simulations were conducted with Stella software using fourth-order Runge–Kutta numerical integration method. According to a previous study [10], lake area of Linggo Co changed from 94.75 km² to 97.59 km² from 1976 to 1992. Compared with the rapid lake expand rate thereafter, this period was relatively stable which was treated as hydrological steady status. Therefore, in the first step of calibration, the model was forced by monthly average data of 1979–1992 (Supplementary Table S1) and the adjusted parameters A and b in Equation (13) to make the simulated lake area close to the observed lake surface in September 1992 (97.59 km²). There are several combinations of A and b that corresponded to the single K_total (

Table 1. Parameters in model calibration.

Simulation		Observsed Lake Area *	1st	2rd	3th
A			39.35	1.17	0.013
B			0.0001	0.001	0.002
Lake Area (km2)	Sep-1992	97.59	100.87	100.09	100.09
	Aug-1999	101.35	103.44	102.07	101.57
	Sep-2006	109.33	110.1	107.01	105.03
	Aug-2010	115.21	112.99	108.58	105.7
Mean MSD			0.02	0.03	0.04

Note: * stands for the observation lake area after [10].

). In the second step of calibration, we forced the hydrological steady lake in the model using monthly data from 1979 to 2012. The combination of A and b was selected that made the simulated lake area in 1999 and 2006 as close to the observation as possible (Table 1). The validity of the model was testified by the difference in isotopic values of lake water between observation and simulation.

In order to determine the isotopic sensitivity of Linggo Co to temperature, precipitation amount, and isotopic composition of precipitation, two sets of sensitivity analyses at steady and transient states were conducted.

Simulation at steady state was conducted on monthly step over 500 model years using modern catchment parameters and the average climate data from 1979 to 1992. Between the 1601st and 2008th simulation month, the tested parameter was altered (increased or decreased) by a constant percent, according to its coefficient of variation over 1979–1992 (Supplementary Table S1).

Simulation at transient state was conducted on a monthly step over 234 model years. For the first 200 years, the model was forced by the average data from 1979 to 1992. From the 2401st to 2808th month (January of the 201st year), the model was forced by the monthly data from 1979 to 2012, except for the tested parameters which increase or decrease by the same amount as the previous sensitivity test. As for the isotopic composition of precipitation, it increased by 10% in the summer season (from June to August). The same amount of change and time allow the results of two sensitivity analyses to be compared.

3. Results

3.1. Model Calibration

The adjusted hydraulic conductivity after the first step of calibration is 0.055 m/month, which lies within the range of dense clay (<0.06 m/month [49]). After the second step of calibration, the values of A and b are chosen as (39.35, 0.0001), since the combination led to the closest simulated lake area in August 2010 with observation (Table 1). Moreover, the lake water isotope value at the end of 2012 is −8.22‰, which fits with the observed data (−8.32‰) in May 2015 after stretching (Figure 3), implying the validity of the model.

After two steps of calibration, we can reconstruct the continuous records of lake area changes during 1979–2012 and quantify the amount of each flux input/output to Linggo Co. The hydrological change of Linggo Co can be separated into three stages: insignificant expansion from 1979 to 1992, stable status from 1993 to 1996, and rapid expansion thereafter (Figure 3). Annual precipitation on Linggo Co ranged from 1.49 × 10⁷~4.95 × 10⁷ m³, runoff to Linggo Co had a wider variation from 9.34 × 10⁷~16.14 × 10⁷ m³ per year, evaporation over lake surface ranged from 5.48 × 10⁷~9.04 × 10⁷ m³ per year, and outseepage was relatively stable from 6.37 × 10⁷~7.55 × 10⁷ m³ per year. The net water budget (precipitation minus evaporation) in the catchment varied from 0.00~19.19 × 10⁷ m³ (

Table 2. Simulated annual amount of each flux of Linggo Co.

Year	Precipitation on Lake (107m3)	Runoff into Lake (107 m3)	Evaporation (107 m3)	Outseepage (107 m3)	Glacier Melt (107 m3)	P-E in Catchment (107 m3)	Glacier Melt/Inflow * (%)
1979	2.25	9.34	6.63	6.39	7.52	3.75	64.84
1980	2.86	10.02	5.63	6.37	7.70	3.20	59.82
1981	2.52	10.53	5.92	6.39	10.75	2.43	82.44
1982	2.61	10.87	5.76	6.40	8.80	2.33	65.29
1983	2.47	10.86	5.83	6.42	9.22	0.95	69.12
1984	1.63	10.70	5.88	6.43	8.09	1.46	65.57
1985	3.02	10.60	5.48	6.44	10.01	1.44	73.53
1986	2.36	10.71	6.24	6.46	10.43	0.19	79.82
1987	1.68	10.73	5.94	6.47	10.05	0.52	80.94
1988	1.92	10.63	6.01	6.47	9.84	0.18	78.41
1989	2.01	10.54	5.61	6.47	10.01	0.44	79.76
1990	2.85	10.83	5.59	6.49	11.09	3.72	81.03
1991	2.53	13.93	7.54	6.53	12.17	3.04	73.96
1992	2.05	11.40	5.56	6.55	9.76	0.00	72.56
1993	2.33	11.32	5.50	6.58	11.20	0.36	82.04
1994	1.49	13.46	8.77	6.59	10.15	0.00	67.86
1995	1.85	12.81	9.04	6.58	9.87	2.84	67.38
1996	4.74	12.99	7.68	6.60	9.64	19.19	54.42
1997	2.40	12.72	6.25	6.64	8.68	0.29	57.40
1998	2.71	15.73	8.61	6.69	12.18	5.73	66.02
1999	3.25	15.51	8.62	6.74	7.61	10.94	40.56
2000	3.84	13.88	6.03	6.80	10.26	6.46	57.94
2001	2.96	14.10	6.28	6.88	8.87	0.68	52.01
2002	4.27	13.77	5.86	6.95	9.88	5.78	54.74
2003	3.52	13.73	6.40	7.03	10.43	0.17	60.46
2004	2.97	16.14	8.54	7.09	10.96	3.70	57.37
2005	3.11	12.96	6.63	7.13	10.77	0.25	66.98
2006	3.54	13.02	6.81	7.17	11.81	5.73	71.30
2007	3.98	16.11	8.58	7.23	11.38	4.81	56.64
2008	4.09	15.97	8.09	7.30	9.89	0.69	49.29
2009	4.48	13.30	7.60	7.35	12.10	7.22	68.06
2010	3.91	13.88	7.39	7.40	12.90	3.48	72.51
2011	4.95	14.27	6.53	7.46	9.28	5.74	48.27
2012	4.79	15.01	6.58	7.55	10.56	12.77	53.35

Note: * inflows of water to Linggo Co were calculated as the sum of precipitation on lake and runoff into lake.

).

3.2. Sensitivity Analysis of Model

In response to the 10% increase in precipitation amount, the lake area at steady status expanded by ~5.3 km², while by ~7.2 km² at transient simulation (red curve in Figure 4a,b). Decreasing precipitation by 10% leads to lake area shrinking ~2.4 km² at steady status and ~6.2 km² at transient status (blue curve in Figure 4a,b). As for δ¹⁸O, increase by 10% in precipitation amount causes δ¹⁸O value to decrease by 0.25‰ at steady status and 0.33‰ at transient simulation (red curve in Figure 4c,d). A 10% decrease in precipitation causes δ¹⁸O value to increase 0.14‰ at steady status and 0.33‰ at transient simulation (blue curve in Figure 4c,d). δ¹⁸O value of precipitation in the summer (June–August) increasing by 10% (~0.15‰ in value) leads to a 0.65‰ shift at steady status and 0.67‰ at transient status (pink curve in Figure 4c,d).

A 32.8% increase in temperature led to the lake area shrinking in the first 5 months; after that, the lake expanded ~0.8 km² at steady status. The δ¹⁸O experienced the same trend: increase in the first five months and then decrease. A 32.8% decrease in temperature, that decreased the actual evaporation as much as 0.30× 10⁷ m³ per year, gave rise to the lake expanding 2 km² and δ¹⁸O decreasing 0.10‰. Little variation occurred in lake area and isotopic compositions under temperature changes in the transient simulation (Figure 5).

4. Discussion

4.1. The Role of Glacier Melt to Lake Expand in Recent Decades

Previous studies have suggested that many factors, such as glacier meltwater [10,50], precipitation [51], and potential evaporation [52,53,54], influence the water budget in the central TP. The difference between simulated results at steady and transient states can provide some guidance in calculating the water budget. In response to precipitation amount change, the variation of lake area at transient status was larger than that at steady status (Figure 4). It implies that the increase in precipitation alone cannot fully account for the lake expansion during the past decade. As the temperature changes, the lake area changes at both statuses are limited. It indicated that the increasing temperature can accelerate glacier melt, while its contribution to the water budget may be offset by enhanced evaporation at the same time.

As for the hydrologic change in Linggo Co, the trend that accelerated expansion in the mid-1990s is consistent with other closed lakes in the central TP [55]. The melt of Puruogangri ice sheet also increased from 0.75 × 10⁸ m³ in 1979 to 1.06 × 10⁸ m³ in 2012, which is also evidenced in others glacier in the TP [50]. From 1979 to 1993, glacier meltwater contributed 73.94% of the average inflow of water (runoff into lake plus precipitation on lake) to Linggo Co. From 1993 to 2012, the contribution of glacier meltwater to runoff decreased to 59.67% on average (Table 2, Figure 6). In other words, increasing glacier meltwater, together with positive water balance (precipitation/evaporation) in the catchment, contributed to the rapid expansion of Linggo Co after 1993.

4.2. The Isotopic Response of Linggo Co to Temperature and Precipitation Amount

The δ¹⁸O of lacustrine carbonates (e.g., ostracod shells), representing the δ¹⁸O composition of lake water, was widely used in paleoclimate studies in the TP [24]. The interpretation of this proxy varies, such as temperature variation [56], lake water balance that responds to the amount of monsoonal precipitation [57], and the interaction between Asian summer monsoon and mid-latitude westerlies [13,58]. The fundamental question is how the lake water δ¹⁸O responded to changes in temperature, precipitation amount, and precipitation δ¹⁸O.

Based on our sensitivity analysis results, isotopic composition of lake water in Linggo Co was least sensitive to temperature variation, relatively sensitive to precipitation amount, and sensitive to precipitation isotopic composition. These results could shed light on the interpretation of the ~4‰ shift in δ¹⁸O record of ostracod shells in Linggo Co in the early Holocene [29]. Increasing temperature could not explain this range of δ¹⁸O fluctuation. The precipitation amount would increase in response to enhanced Asian summer monsoon during the early Holocene. However, the increased precipitation amount may not be large enough to support the isotopic variation, considering that a 10% variation in precipitation amount led to 0.33‰ isotopic variation of lake water during 1979–2012. On the other hand, moisture induced by Asian summer monsoon and westerlies has large isotopic differences. For example, the δ¹⁸O difference of summer precipitation between Nagqu and Gaize is 6.4‰ (Nagqu: −17.0‰, Gaize: −10.6‰ [59]). Meanwhile, δ¹⁸O value of summer precipitation varies by ~0.15‰ and can lead to ~0.67‰ change in lake water δ¹⁸O value based on our sensitivity analyses. Thus, it is highly possible that the shift in δ¹⁸O records of Linggo Co was caused by a change in atmospheric circulation systems.

5. Conclusions

In this study, we first carried out the hydrologic and isotopic mass balance model in lakes on the TP. We simulated the hydrologic change of Linggo Co during 1979–2012 to quantify the contribution of glacier meltwater to the expansion of Linggo Co and assess the influences of temperature, precipitation amount, and precipitation isotopic composition on lake water isotopic composition. Our findings are as follows:

1. Hydrologic change of Linggo Co during 1979–2012 can be separated into two stages: before 1993, meltwater contributed 90% of the runoff, while after that, the contribution of glacier meltwater to runoff decreased to at least 49% (in 1999). Increasing glacier melting together with positive water balance (precipitation/evaporation) in the catchment contributed to the rapid expansion of Linggo Co.

2. Isotopic composition of lake water in Linggo Co was not sensitive to temperature variation, relatively sensitive to precipitation amount, and sensitive to precipitation isotopic composition. A large shift in δ¹⁸O record of Linggo Co during the early Holocene was caused by the interplay between Asian summer monsoon and mid-latitude westerlies.