Evaluation of lacustrine groundwater discharge, hydrologic partitioning, and nutrient budgets in a proglacial lake in the Qinghai–Tibet Plateau: using ²²²Rn and stable isotopes

Proglacial lakes are good natural laboratories to investigate groundwater and glacier dynamics under current climate conditions and to explore biogeochemical cycling under pristine lake status. This study conducted a series of investigations of ²²²Rn, stable isotopes, nutrients, and other hydrogeochemical parameters in Ximen Co Lake, a remote proglacial lake in the east of the Qinghai–Tibet Plateau (QTP). A radon mass balance model was used to quantify the lacustrine groundwater discharge (LGD) of the lake, leading to an LGD estimate of 10.3±8.2 mm d⁻¹. Based on the three-endmember models of stable ¹⁸O and Cl⁻, the hydrologic partitioning of the lake is obtained, which shows that groundwater discharge only accounts for 7.0 % of the total water input. The groundwater-derived DIN and DIP loadings constitute 42.9 % and 5.5 % of the total nutrient loading to the lakes, indicating the significance of LGD in delivering disproportionate DIN into the lake. This study presents the first attempt to evaluate the LGD and hydrologic partitioning in the glacial lake by coupling radioactive and stable isotopic approaches and the findings advance the understanding of nutrient budgets in the proglacial lakes of the QTP. The study is also instructional in revealing the hydrogeochemical processes in proglacial lakes elsewhere.

1 Introduction

High-altitude and high-latitude areas are intensively influenced by the melting of glaciers due to climatic warming. Of particular importance are the proglacial areas, such as proglacial lakes and moraines, because they are particularly affected by climatic-change-induced glacier retreating and thawing of permafrost (Heckmann et al., 2016; Barry, 2006; Slaymaker, 2011). The proglacial lakes are usually located close to the ice front of a glacier, ice cap, or ice sheet, with the vicinity to the ice front sometimes defined as the areas with subrecent moraines and formed by the last significant glacier advances at the end of the Little Ice Age (Heckmann et al., 2016; Barry, 2006; Slaymaker, 2011; Harris et al., 2009). These lakes are located in the transition zones from glacial to non-glacial conditions, and can serve as natural laboratories to explore hydrological processes, biogeochemical cycles, and geomorphic dynamics under current climatic conditions (Dimova et al., 2015; Heckmann et al., 2016).

The Qinghai–Tibet Plateau (QTP), the third pole of the world, serves as the water tower of most of the major rivers in Asia (Qiu, 2008). Unique landscapes such as endorheric lakes, permafrost, glaciers, and headwater fluvial networks are developed due to the intensive interaction between the atmosphere, hydrosphere, biosphere, and cryosphere (Lei et al., 2017; Zhang et al., 2017a, b; Yao et al., 2012, 2013). Distributed mountainous glaciers and lakes are the most representative landscapes and are highly sensitive to the climate changes. In the past decade, the lakes in the interior of the QTP show overall expansion with respect to an overall increase in precipitation, accelerated glacier melting, and permafrost degradation (Zhang et al., 2013, 2017a, b; Heckmann et al., 2016; Yang et al., 2014). Some recent studies have made efforts to depict the hydrologic partitioning of the majority of the lakes in the QTP based on long-term observation of climatological parameters and remote sensing approaches. However, so far a quantitative evaluation of the water balance and hydrologic partitioning, especially the groundwater component of the lakes in the QTP, has been limited due to the scarcity of observational data. Therefore, there is a great need to conduct refined and systematical field observation to provide a groundtruth dataset and tenable models to depict the water balance and hydrologic partitioning of the lakes, especially proglacial lakes in the QTP (Yang et al., 2014; Zhang et al., 2017b).

Mountainous proglacial lakes, formed by glacial erosion and filled by melting glaciers, are widely distributed in the Qinghai–Tibet Plateau (QTP), especially along the substantial glacier retreating areas of the Himalaya Mountains (MT.), Qilian MT., Tienshan MT., etc. Characterized by higher elevations, small surface areas but relatively large depths, mountainous proglacial lakes in QTP lack systematic field-based hydrological studies due to their remote locations and difficulty in conducting fieldwork (Yao et al., 2012; Farinotti et al., 2015; Bolch et al., 2012).

There has been extensive recognition of the importance of groundwater discharge to various aquatic systems for decades (Dimova and Burnett, 2011; Valiela et al., 1978; Johannes, 1980). Very recently, the topic of lacustrine groundwater discharge (LGD), which is comprehensively defined as groundwater exfiltration from lakeshore aquifers to lakes (Lewandowski et al., 2013, 2015; Rosenberry et al., 2015; Blume et al., 2013), has been introduced. LGD is analogous to submarine groundwater discharge (SGD) in coastal environments. LGD also plays a vital role in lake hydrologic partitioning, which is defined as the separation of groundwater discharge/exfiltration, riverine inflow, riverine outflow infiltration, surface evaporation, and precipitation for the hydrological cycle of the lake (Luo et al., 2017; Good et al., 2015). LGD also serves as an important component in delivering solutes to lakes since groundwater is usually concentrated in nutrients, CH₄, dissolved inorganic/organic carbon (DIC/DOC), and other geochemical components (Paytan et al., 2015; Lecher et al., 2015; Belanger et al., 1985; Dimova et al., 2015). Nutrients and carbon loading from groundwater greatly influence ratios of dissolved inorganic nitrogen (DIN) to dissolved inorganic phosphate (DIP) (referred to as N:P ratios thereafter), ecosystem structure, and the primary productivity of the lake aquatic system (Nakayama and Watanabe, 2008; Belanger et al., 1985; Hagerthey and Kerfoot, 1998).

The approaches to investigate LGD include (1) direct seepage meters (Shaw and Prepas, 1990; Lee, 1977), (2) geo-tracers such as radionuclides and stable ²H and ¹⁸O isotopes (Gat, 1995; Kluge et al., 2007; Kraemer, 2005; Lazar et al., 2008), (3) heat and temperature signatures (Liu et al., 2015; Sebok et al., 2013), (4) numerical modelling (Winter, 1999; Smerdon et al., 2007; Zlotnik et al., 2009, 2010), and (5) remote sensing (Lewandowski et al., 2013; Wilson and Rocha, 2016; Anderson et al., 2013). Recently, some researchers started to investigate groundwater dynamics in periglacial and proglacial areas, mostly based on the approaches of numerical modelling (Lemieux et al., 2008a, b, c; Andermann et al., 2012; Scheidegger and Bense, 2014). However, the quantification of groundwater and surface water exchange in proglacial lakes is still challenging due to limited hydrogeological data and extreme seasonal variability of aquifer permeability (Dimova et al., 2015; Callegary et al., 2013; Xin et al., 2013).

²²²Rn, a naturally occurring inert gas nuclide highly concentrated in groundwater, can be more applicable in fresh aquatic systems and has been widely used as a tracer to quantify groundwater discharge in freshwater lakes (Luo et al., 2016; Corbett et al., 1997; Dimova et al., 2013, 2015; Dimova and Burnett, 2011; Kluge et al., 2007, 2012; Schmidt et al., 2010) and terrestrial rivers and streams (Burnett et al., 2010; Cook et al., 2003, 2006; Batlle-Aguilar et al., 2014). Of particular interest are investigations of temporal ²²²Rn distribution in lakes, since it can be used to quantify groundwater discharge and reflect the locally climatological dynamics (Dimova and Burnett, 2011; Luo et al., 2016). Temporal radon variations give high-resolution estimates of groundwater discharge to lakes over diel cycles, allowing evaluation of LGD and the associated chemical loadings. However, there has been no study of radon-based groundwater discharge in mountainous proglacial lakes, especially for those lakes in the QTP.

This study aims to investigate the groundwater–surface water interactions for the proglacial lake of Ximen Co, by estimating the LGD and evaluating the hydrologic partitioning of the lake. LGD is estimated with the ²²²Rn mass balance model, and the hydrologic partitioning of the lake is obtained with the three-endmember model coupling the mass balance of water, stable isotopes, and Cl⁻. Moreover, LGD-derived nutrients are estimated and the nutrient budgets of the lake are depicted. This study, to our knowledge, makes the first attempt to quantify the LGD, hydrologic partition, and groundwater-borne nutrients of the proglacial lake in the QTP and elsewhere via the approach of integrating multiple tracers. This study provides insights into hydrologic partitioning in a typical mountainous proglacial lake under current climate conditions and reveals groundwater-borne chemical loadings in this proglacial lake in the QTP and elsewhere.

2 Methodology

2.1 Site descriptions

The Nianbaoyeze MT., located at the eastern margin of the QTP and being the easternmost part of NW–SW trending Bayan Har Shan, is situated at the main water divide of the upper reaches of the Yellow River and Yangtze River (Fig. 1). With a peak elevation of 5369 m, the mountain rises about 500–800 m above the surrounding peneplain and displays typical Pleistocene glacial landscapes such as moraines, U-shaped valleys, and cirques (Lehmkuhl, 1998; Schlutz and Lehmkuhl, 2009; Wischnewski et al., 2014). The present snow line is estimated to be at an elevation of 5100 m (Schlutz and Lehmkuhl, 2009). Controlled by the South Asia and East Asia monsoons, the mountain has an annual precipitation of 975 mm in the southern part and 582 mm in the northwestern part, with 80 % occurring during May and October (Yuan et al., 2014; Zhang and Mischke, 2009). The average temperature gradient is about 0.55 ^∘C per 100 m, and the closest weather station, located in Jiuzhi (N: 33.424614^∘, E: 101.485998^∘) at the lower plains of the mountain, recorded a mean annual temperature of 0.1 ^∘C. Snowfalls occur in nearly 10 months of the entire year and there is no free frost all year around (Böhner, 1996, 2006; Schlutz and Lehmkuhl, 2009). The precipitation, daily bin-averaged wind speed, and temperature in August 2015 were recorded to be 90 mm, 0.7 m s⁻¹, and 9.5 ^∘C from Jiuzhi weather station (Fig. 2). The water surface evaporation was recorded to be 1429.8 mm in 2015 from Jiuzhi weather station.

Figure 1 The geological and topographic map of the Yellow River source region, Nianbaoyeze glacial mountains (a), and sampling settings of Ximen Co Lake (b), with the bathymetry map of the lake (d). (c) Photograph of Ximen Co Lake and the surrounding geomorphic settings looking in a northeasterly direction on 22 August 2015, showing the late-laying snowpack in the U-shaped valleys of the northern part of the Nianbaoyeze mountains.

Figure 2 The climatological parameters (wind speed, air temperature, and precipitation) in August 2015 recorded from Jiuzhi weather station.

Among the numerous proglacial lakes developed in the U-shaped valleys of the Nianbaoyeze MT., Ximen Co Lake is located at the northern margin of the mountain with an elevation of 4030 m a.s.l., and is well studied and easily accessible (Lehmkuhl, 1998; Zhang and Mischke, 2009; Schlutz and Lehmkuhl, 2009; Yuan et al., 2014). The lake was formed in a deep, glacially eroded basin with a catchment area of 50 km², and has a mean and a maximum depth of 40 and 63.2 m, and a surface area of 3.6 km². The vegetation around the lake is dominated by pine meadows with dwarf shrubs, rosette plants, and alpine cushion (Schlutz and Lehmkuhl, 2009; Zhang and Mischke, 2009; Yuan et al., 2014). Mostly recharged by the glacial and snowpack melting water and regional precipitation, the lake is stratified with an epilimnion depth of about 4.4 m in the summertime. The lake is usually covered by ice in the wintertime (Zhang and Mischke, 2009). The superficial layer within the U-shaped valley is characterized by peat, clay, and fluvial gravels with a depth of about 1–3.5 m. Discontinuous and isolated permafrost is present at the slope of the valley above an elevation of about 4150 m. The maximum frozen depth is about 1.5 m for the seasonal frozen ground around the lake. The seasonal frozen ground serves as an unconfined aquifer during the unfrozen months from July to October, and groundwater discharges into the epilimnion of the lake (Wang, 1997; Schlutz and Lehmkuhl, 2009; Zhang and Mischke, 2009).

2.2 Sampling and field analysis

The field campaign to Ximen Co Lake was conducted in August 2015, when it is warm enough to take the water samples of different origins as the studied site is seasonally frozen. A ²²²Rn continuous monitoring station was set up in the southeast part of the lake, where is fairly flat for setting up our tent and monitoring system. Surface water samples were collected around the lake and rivers at the upstream and downstream. Porewater samples were collected at one side of the lake as the other side is steep and rocky. The basic water quality parameters of conductivity (EC), dissolved oxygen (DO), TDS, ORP, and pH in the water were recorded with the multi-parameter meter (HANNA, Co.). Relative humidity was recorded with a portable thermo-hydrometer (KTH-2, Co.). Lake water samples were taken with a peristaltic pump into 2.5 L glass bottles for ²²²Rn measurement with the Big Bottle system (Durridge, Co.). Surface water samples were filtered with 0.45 µm filters (Advantec, Co.) in situ and taken into 5, 15, 15, and 50 ml Nalgene centrifugation tubes for stable isotope, major anion, cation, and nutrient analysis. Porewater samples were taken from the lakeshore aquifers with a push point sampler (M.H.E, Co.) connected to a peristaltic pump (Solinst, Co.) (Luo et al., 2014, 2016); 100 ml raw surface water or porewater was titrated with a 0.1 µM H₂SO₄ cartridge (Hach, Co.) in situ to measure total alkalinity (Hasler et al., 2016; White et al., 2016; Warner et al., 2013). Porewater was filtered with 0.45 µm syringe filters in situ and taken into 5, 15, 15, and 50 ml Nalgene centrifugation tubes for stable isotope, major anion, cation, and nutrient analysis; 250 ml porewater was taken for ²²²Rn measurement with RAD7 H₂O (Durridge, Co.). Samples for major cation analysis were acidified with distilled HNO₃ immediately after the sampling.

A ²²²Rn continuous monitoring station was set up in the northwest of the lake, close to the downstream of the lake (Fig. 1b). Lake water (about 0.5 m in depth) was pumped with a DC pump (12 V) driven by lithium batteries (100 Ah) and sprinkled into the chamber of RAD7 AQUA with a flow rate >2 L min⁻¹, where ²²²Rn in water vapour was equilibrated with the air ²²²Rn. The vapour in the chamber was delivered into two large dry units (Drierite, Co) to remove the moisture and circulated into a RAD7 monitor, where ²²²Rn activities were recorded every 5 min. A temperature probe (HOBO^@) was inserted into the chamber to record the temperature of the water vapour. The monitoring was performed from 11:31, 22 August to 06:30, 24 August 2015. During the period of 13:50–16:30 on 22 August, a sudden blizzard occurred, leading to an hourly precipitation of about 0.6 mm onto the lake area. Daily and hourly climatological data such as wind speed, air temperature, and precipitation were retrieved from the nearest weather station in Jiuzhi (N: 33.424614^∘, E: 101.485998^∘). Moreover, another RAD7 was placed at the lake shore to measure ²²²Rn in the ambient air around the lake. Due to extremely low activities, the monitoring period was conducted only for 4 h, and the mean activity was adopted as the background radon-222 activity to be used in the mass balance model. Water level and temperature fluctuations were recorded with a conductivity–temperature–depth diver (Schlumberger, Co.) fixed at about 20 cm below the lake surface and calibrated with local atmospheric pressure recorded by a baro-diver (Schlumberger, Co.) above the lake. To correct for dissolved ²²⁶Ra supported ²²²Rn, one radium sample was extracted from 100 L lake water with MnO₂ fibre as described elsewhere (Luo et al., 2014; Moore, 1976).

2.3 Chemical analysis

Major ions were measured with ICS-1100 (Dionex. Co.) in the Department of Earth Sciences, University of Hong Kong. The uncertainties of the measurements are less than 5 %. Nutrients, DIN, and DIP were analysed with flow injection analysis equipped with an auto-sampler (Lachat. Co.) in the School of Biological Sciences, University of Hong Kong. Stable ¹⁸O and ²H isotopes were measured with a MOA-ICOS laser absorption spectrometer (Los Gatos Research (LGR) Triple Isotope Water Analyzer (TIWA-45EP)) at the State Key Laboratory of Marine Geology, Tongji University, Shanghai. The stable isotopic standards and the recovery test have been fully described elsewhere (Luo et al., 2017). The measurement uncertainty is better than 0.1 % for ¹⁸O and 0.5 % for ²H. ²²⁶Ra was detected with RAD7 with the method described elsewhere (Kim et al., 2001; Lee et al., 2012; Luo et al., 2018).

2.4 Radon transient model

Previous studies employed a steady-state radon-222 mass balance model to quantify LGD to lentic systems such as lakes and wetlands (Dimova and Burnett, 2011; Luo et al., 2016). This model assumes that radon inputs derived from groundwater inflow, diffusion, and river inflow are balanced by the radon losses of atmospheric evasion, decay, and river outflow. However, recent studies revealed that the steady state is mainly reached after 2–15 days of constant metrological conditions, and most lentic systems cannot be treated as steady state due to rapid radon-222 degassing to the atmosphere driven by wind-induced turbulence (Gilfedder et al., 2015; Dimova and Burnett, 2011).

Ximen Co Lake has been demonstrated to be highly stratified, with an epilimnion of 4.4 m (Zhang and Mischke, 2009). The lake was formed by glacier erosion and the lake bed is characterized by granite bedrock with a thin sedimentary clay layer. Previous studies have indicated that sediment consisting of clay, soils, and gravels has developed on the bedrock and forms the lakeshore aquifer with a thickness of 0.7–3.3 m (Schlutz and Lehmkuhl 2009). Porewater sampled in the aquifer immediately behind the lake shore can well represent groundwater discharging into the lake, as suggested previously (Lewandowski et al., 2015; Rosenberry et al., 2015; Schafran and Driscoll, 1993). LGD has been widely considered to occur within the first few metres of the lake shore (Schafran and Driscoll, 1993; Rosenberry et al., 2015; Lee et al., 1980) and groundwater is considered to predominately discharge into the epilimnion since deep groundwater flow is highly limited by the Precambrian bedrock (Einarsdottir et al., 2017). Due to a negligible hydrological connection between the epilimnion and hypolimnion, the ²²²Rn mass balance model is established to quantify LGD to the epilimnion from the lake shore.

The governing equation of the radon-222 transient mass balance model within $\mathrm{1}\times \mathrm{1}\times z$ cm (where z is the depth in centimetres) can be expressed as (Gilfedder et al., 2015)

$$\begin{array}{}\text{(1)}& z{\displaystyle \frac{\partial I_{\mathrm{w}}}{\partial t}}=F_{\mathrm{gw}}+(I_{{}^{\mathrm{226}}\mathrm{Ra}}-I_{\mathrm{w}})\times z\times {\mathit{\lambda }}_{\mathrm{222}}+F_{\mathrm{diff}}-F_{\mathrm{atm}},\end{array}$$

where F_gw, F_diff, and F_atm (Bq m⁻² d⁻¹) are ²²²Rn loadings from LGD, water–sediment diffusion, and water–air evasion, respectively; z (m) is the lake water level depth recorded by the diver. λ₂₂₂ is the decay constant of ²²²Rn with a value of 0.186 d⁻¹. ${\mathit{\lambda }}_{\mathrm{222}}\times I_{{}^{\mathrm{226}}\mathrm{Ra}}$ and λ₂₂₂×I_w account for the production and decay of ²²²Rn (Bq m⁻² d⁻¹) in the water column, respectively. I_w and $I_{{}^{\mathrm{226}}\mathrm{Ra}}$ (Bq m⁻²) represent ²²²Rn and ²²⁶Ra inventories in the epilimnion, and are expressed as $I_{\mathrm{w}}=H\times C_{\mathrm{w}}$ and $I_{{}^{\mathrm{226}}\mathrm{Ra}}=H\times C_{{}^{\mathrm{226}}\mathrm{Ra}}$, respectively, where H (m) is the depth of the epilimnion; C_w and $C_{{}^{\mathrm{226}}\mathrm{Ra}}$ are the ²²²Rn and ²²⁶Ra activity (Bq m⁻³), respectively.

The model is valid under the following assumptions. (1) The epilimnion is well mixed, which is the actual condition for most natural boreal and high-altitude glacial lakes (Zhang and Mischke, 2009; Åberg et al., 2010). (2) ²²²Rn input from riverine water inflow, and loss from the lake water outflow and infiltration into the lakeshore aquifer, are negligible compared to the groundwater-borne ²²²Rn, because ²²²Rn concentration of groundwater is 2–3 orders of magnitude larger than that of lake water (Dimova and Burnett, 2011; Dimova et al., 2013). Generally, ²²²Rn in the epilimnion is sourced from LGD and decay input from the parent isotope of ²²⁶Ra under secular equilibrium, and is mainly lost via atmospheric evasion and radioactive decay.

F_atm is the key sinking component of the transient model and is finally a function of wind speed and water temperature, both of which are temporal variant variables (information in the Supplement). Lake water level z is also a temporal variant variable which represents the fluctuations of water volume of the epilimnion. This equation is discretized by the forward finite difference method, and the groundwater flux at each time step can be solved as follows:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{[}^{\mathrm{222}}{\mathrm{Rn}}_{t+\mathrm{\Delta }t}]=\\ \text{(2)}& {\displaystyle }& {\displaystyle }{\displaystyle \frac{[z{\times }^{\mathrm{222}}{\mathrm{Rn}}_t+[F_{\mathrm{diff}}+F_{\mathrm{gw}}-F_{\mathrm{atm}}{-}^{\mathrm{222}}{\mathrm{Rn}}_t\times \mathit{\lambda }\times z]\times \mathrm{\Delta }t}{z}},\end{array}$$

where ²²²Rn_t+Δt and ²²²Rn_t+Δt (Bq m⁻³) are the ²²²Rn activity at the current time step and at the previous time steps, respectively, and Δt (min) is the time step which is set to be 5 min, consistent with the ²²²Rn record interval. With the inverse calculation based on Eq. (2), the groundwater inflow at each time step can be obtained. However, large errors of the final LGD calculation will be induced by even a small amount of noise in the measured ²²²Rn data due to the ${}^{\mathrm{222}}{\mathrm{Rn}}_{t+\mathrm{\Delta }t}{-}^{\mathrm{222}}{\mathrm{Rn}}_t$ term being with the measured uncertainty. To reduce the random errors of the measured ²²²Rn concentrations, the time window with a width of 1 h is proposed to smooth the curve (information in the Supplement).

Figure 3 The temporal distributions of ²²²Rn (a), water temperature (b), water level fluctuation recorded by the divers (c), and hourly wind speed and air temperature recorded in Jiuzhi weather station (d).

3 Results

3.1 Time series data

Figure 2 shows the basic climatological parameters of the lake catchment during the campaign month. There are discrete rainfall events occurring throughout the month with an average rainfall of 3.1 mm d⁻¹. The temperature during the month ranges from 5.0 to 12.5 ^∘C, with an average of 9.3 ^∘C. The daily averaged wind speed ranges from 0.7 to 2.5 m s⁻¹, with an average of 1.7 m s⁻¹. ²²²Rn temporal distribution and other time series data are shown in Fig. 3a and listed in Table S1 in the Supplement. Generally, ²²²Rn concentration varies from 32.2 to 273 Bq m⁻³, with an average of 144.2±27.7 Bq m⁻³. ²²²Rn over the monitoring period shows a typical diel cycle, much higher at nighttime and lower in the daytime. Figure 3b–d show the time series data of temperature (5 min interval), nearshore lake water level (1 min interval), and wind speed (1 h interval). Temperature and lake water level also show typical diel cycles, but with antiphase fluctuations with each other. Temperature is higher during the daytime and lower at nighttime. However, a sudden decrease in temperature was recorded due to the sudden blizzard (Fig. 3b). Water level is higher at nighttime and lower during the daytime, with a strong fluctuation due to the turbulence caused by the blizzard (Fig. 3c). The variability might reflect the dynamics of groundwater input and surface water inflow. The air temperature of the lake area is in phase with the water temperature. Wind speed is normally higher during the daytime and lower at nighttime (Fig. 3d).

The variation of ²²²Rn is nearly in antiphase with the fluctuations of lake water temperature and air temperature, indicating that the dominated controlling factors of ²²²Rn fluctuations are water temperature and wind speed (Fig. 3a). This phenomenon is reasonable as lake water ²²²Rn is predominately lost via atmospheric evasion, which is the function of wind speed and water temperature (Dimova et al., 2013, 2015; Dimova and Burnett, 2011). High water temperature and wind speed lead to elevated atmospheric evasion and cause the decline of ²²²Rn concentration in the lake water. However, there is a sudden reduction of radon activity from 14:00 to 16:00 (GMT + 8) on 22 July 2015, when the snow event led to a sudden decrease in water temperature, increase in wind speed, and large surface water turbulence as indicated by water level fluctuations (Fig. 3a–d). ²²²Rn in the porewater is 2–3 orders of magnitude larger than ²²²Rn in the lake water, suggesting that ²²²Rn is an ideal tracer to estimate the LGD (Table S1). ²²²Rn concentrations in surface water range from 22.2 to 209 Bq m⁻³, with an average of 92.5 Bq m⁻³ (n=12), which is in the range of ²²²Rn continuous monitoring results, suggesting reliable ²²²Rn measurements (Table S2).

Figure 4 The cross plots of Cl⁻ versus Na⁺ (a), ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ versus Cl⁻ (b), and Ca²⁺ versus Cl⁻ (c). The relations of ²H versus ¹⁸O (d), Cl⁻ versus ²H (e), and Cl⁻ versus ¹⁸O (f).

Figure 5 Cross plots of ²²²Rn versus DIN (a) and DIP (b).

3.2 Geochemical results

The results of major ions, nutrients, and stable isotopes in different water endmembers are shown in Figs. 4 and 5. Cl⁻ ranges from 0.6 to 2.1 mg L⁻¹ in the surface water (including riverine inflow water, lake water, and downstream water) and 0.4 to 2.7 mg L⁻¹ in porewater, and has a much higher concentration of 5.9 mg L⁻¹ in rainfall water. Na⁺ ranges from 1.6 to 3.4 mg L⁻¹ in the surface water and 1.2 to 4.4 mg L⁻¹ in porewater and has a concentration of 4.4 mg L⁻¹ in rainfall water. ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ ranges from 1.2 to 2.3 mg L⁻¹ in the surface water and 0.4 to 1.7 mg L⁻¹ in porewater, and has a significant low concentration of 0.01 mg L⁻¹ in rainfall water. Ca²⁺ ranges from 3.0 to 12.4 mg L⁻¹ in lake water and 3.4 to 12.5 mg L⁻¹ in porewater, and has a significantly high concentration of 20.5 mg L⁻¹ in rainfall water. Other concentrations of major ions are listed in Table S2. As shown in Fig. 4d and Table S2, δ¹⁸O in the lake water ranges from −13.06 ‰ to −12.11 ‰, with an average of −12.41 ‰ (n=7), and δ²H ranges from −91.83 ‰ to −87.47 ‰, with an average of −89.0 ‰ (n=7). δ¹⁸O in the riverine inflow water ranges from −13.44 ‰ to −13.29 ‰, with an average of −13.37 ‰ (n=2), and δ²H ranges from −93.25 ‰ to −91.92 ‰, with an average of −92.59 ‰ (n=2). δ¹⁸O in the downstream water ranges from −12.51 ‰ to −12.18 ‰, with an average of −12.35 ‰ (n=3), and δ²H ranges from −88.96 ‰ to −87.1‰, with an average of −87.98 ‰ (n=3). δ¹⁸O in the porewater ranges from −12.66 ‰ to −11.52 ‰, with an average of −11.97 ‰ (n=8), and δ²H ranges from −91.3 ‰ to −82.87 ‰, with an average of −85.5 ‰ (n=8). DIN in the surface water (including riverine inflow water, lake water, and downstream water) ranges from 6.6 to 16.9 µM, with an average of 10.3 µM, and DIP from 0.36 to 0.41 µM, with an average of 0.38 µM. DIN for the porewater ranges from 0.7 to 358.8 µM, with an average of 92.8 µM, and DIP from 0.18 to 0.44 µM, with an average of 0.31 µM (Fig. 5).

4 Discussion

4.1 Proglacial hydrologic processes and geochemical implications

Generally, major ion concentrations in the lake water and porewater of Ximen Co Lake are significantly lower than those in major rivers, streams, and other tectonic lakes in the QTP (Yao et al., 2015; Wang et al., 2010, 2016b), and are similar to those of snow and glaciers (Liu et al., 2011), suggesting that the lake water mainly originated from glacier and snow melting. Ion concentrations in the lake and porewater of Ximen Co Lake are much lower than those of rainfall collected in Jiuzhi. This suggests that lake water is less influenced by precipitation (Fig. 4a–c). The concentrations of major ions in the porewater are high compared to the lake water, indicating weathering effects from the aquifer grains. The ratios of Ca²⁺ ∕ Na⁺ in the porewater and groundwater are >1, also suggesting influences of weathering digenesis of major ions from the seasonal frozen ground at the lakeshore aquifer (Weynell et al., 2016; Yao et al., 2015; Wang et al., 2010).

The isotopic compositions of the lake water and porewater are significantly isotopically depleted, with values close to the compositions of glaciers and surface snow in the QTP, suggesting the lake is dominantly recharged from snow and glacier melting (Cui et al., 2014; Wang et al., 2016a; Zongxing et al., 2015). The relation of δ¹⁸O versus δ²H for the lake water is ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}=\mathrm{4.25}\times {\mathit{\delta }}^{\mathrm{18}}\mathrm{O}-\mathrm{35.99}$, with a slope much lower than that of the global meteoric water line (GMWL) (Fig. 4d), suggesting the effects of lake surface evaporation. The relation of δ¹⁸O versus δ²H for the porewater is ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}=\mathrm{6.93}\times {\mathit{\delta }}^{\mathrm{18}}\mathrm{O}-\mathrm{2.67}$, overall on the GWML (Fig. 4d). Deuterium excess is defined as $\mathrm{\Delta }\mathrm{D}=\mathit{\delta }\mathrm{D}-\mathrm{8}\times {\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ (Dansgaard, 1964). The value of ΔD is dependent on airmass origins, altitude effect, and the kinetic effects during evaporation (Hren et al., 2009). Global meteoric water has a ΔD of +10 ‰. In the QTP, glacier/snowpack melting water usually has a large positive ΔD, while the precipitations derived from the warm and humid summer monsoon have a lower ΔD (Ren et al., 2013, 2017). In this study, ΔD of surface water, lake, and porewater ranges from +8.5 to +11.8 ‰, close to the glacier melting water but much smaller than that of the local precipitation of +18.8 ‰. This indicates the stream and lake waters mainly originated from glacial/snowpack melting rather than precipitation (Gat, 1996; Wang et al., 2016a; Lerman et al., 1995). The slopes of δ²H versus δ¹⁸O in lake water and porewater are 4.25 and 6.93, both of which are lower than that of the GMWL due to surface evaporation. Lake water is more intensively influenced by evaporation compared to porewater. The plots of δ¹⁸O versus Cl⁻, and δ²H versus Cl⁻, are well clustered for porewater endmember (orange area), lake water endmember (blue area), riverine inflow water endmember (yellow area), and precipitation water (Fig. 4e and f), suggesting stable δ¹⁸O and δ²H isotopes and Cl⁻ can serve as tracers to quantify the hydrologic partitioning of the lake by setting three-endmember models.

The concentrations of DIN and DIP are all within the ranges of other glacial melting water and proglacial lake water (Hawkings et al., 2016; Hodson, 2005, 2007; Hudson et al., 2000; Tockner et al., 2002). Briefly, rainfall and upstream lake water such as YN-4 have the highest DIN concentrations, indicating the glacier melting and precipitation could be important DIN sources in proglacial areas (Dubnick et al., 2017; Anderson et al., 2017). DIN in porewater is overall higher compared to the lake water, suggesting the porewater to be an effective DIN source, and DIP concentration is higher in the lake water compared to porewater, suggesting the porewater is a DIP sink (Fig. 5). The N:P ratios in the lake water and porewater are averaged to be 27.1 and 320.5, respectively, both much larger than the Redfield ratio ($N:P=\mathrm{16}:\mathrm{1}$) in water and organisms in most aquatic systems, and within the range of other proglacial lakes (Anderson et al., 2017). This also suggests that the lake water and porewater are under phosphate-limited conditions. The N:P ratio in the rainfall water is 30.4, similar to the lake water. The average N:P ratio of porewater is much higher than that of lake water, indicating DIN enrichment in the lakeshore aquifers (Fig. 5). In pristine groundwater, ${\mathrm{NO}}_{\mathrm{3}}^{-}$ is the predominant form of dissolved nitrogen and is highly mobile within the oxic aquifers, leading to much higher DIN concentrations in the porewater; DIP has a high affinity with the aquifer grains, resulting in much lower DIP concentrations in the porewater (Lewandowski et al., 2015; Rosenberry et al., 2015; Slomp and Van Cappellen, 2004). Thus, analogously to surface runoff from glacier/snowpack melting, LGD can also be regarded as an important DIN source for the proglacial lakes. Because of very high DIN and N:P ratios in the porewater, a relatively small portion of LGD delivers considerable nutrients into the glacial lake, shifting the aquatic N:P ratios and affecting the proglacial aquatic ecosystem (Anderson et al., 2017).

Figure 6 The conceptual model of the ²²²Rn transient model (a) and the three-endmember model (b).

4.2 Estimation of LGD

Figure 6a shows all the sinks and sources of radon with the epilimnion of the lake. Within the ²²²Rn transient mass balance model, the dominant ²²²Rn loss is atmospheric degassing/evasion. Generally, the ²²²Rn degassing rate is the function of the radon-222 concentration gradient at the water–air interface and the parameter of gas piston velocity k, which is finally the function of wind speed and water temperature (Dimova and Burnett, 2011; Gilfedder et al., 2015). To evaluate the ²²²Rn evasion rate, this study employs the widely used method proposed by MacIntyre et al. (1995) which is also detailed in the Supplement. Based on the field data of ²²²Rn concentration in the lake water, wind speed, and temperature log, the radon degassing rate is calculated in a range of 0.8 to 265.2 Bq m² d⁻¹, with an average of 42.0 Bq m² d⁻¹.

Table 1 Parameters used to establish the mass balance model of ²²²Rn in Ximen Co Lake.

Figure 7 The results of the final LGD derived from the ²²²Rn transient model.

In addition to the atmospheric loss and sedimentary diffusion inputs, ²²²Rn is also sinked via radioactive decay, and sourced from decay of the parent isotope of ²²⁶Ra. The decay loss of ²²²Rn fluctuates in phase with the distribution of ²²²Rn concentration monitored by RAD 7 AQUA. The equations to estimate benthic fluxes are shown in the Supplement. The decay loss is calculated to be 26.4 to 223.4 Bq m⁻² d⁻¹, with an average of 118.0±22.7 Bq m⁻² d⁻¹. ²²⁶Ra concentration is 0.01 Bq m⁻³ for the lake water. Under secular equilibrium, the ²²⁶Ra decay input can be calculated by multiplying ²²⁶Ra concentration in the lake water by λ₂₂₂ (Corbett et al., 1997; Kluge et al., 2007; Luo et al., 2016). ²²⁶Ra decay input is calculated to be 0.83 Bq m⁻² d⁻¹, which is significantly low compared to other ²²²Rn sources to the epilimnion.

With the obtained sinks and sources of ²²²Rn in the lake, and the constants given in Table 1, LGD rate can be obtained by dividing the groundwater-derived ²²²Rn by its concentration in the groundwater endmember. The obtained LGD rate ranges from −23.7 to 90.0 mm d⁻¹, with an average of 10.3±8.2 mm d⁻¹ (Fig. 7). The LGD rate range is relatively smaller than the daily lake water level variations (≈50 mm), indicating that the lake water level variation could be a combined effect of surface runoff and LGD (Hood et al., 2006). The negative values of LGD rate reflect the return groundwater flow due to infiltration into the porewater. Normally, the dominant values are positive, indicating LGD rate is significant compared to water infiltration into the lakeshore aquifer. The temporal variation of LGD rate could be attributed to the fluctuations of the hydraulic gradient in the proglacial areas (Hood et al., 2006; Levy et al., 2015). As indicated by ΔD (mostly >10) of surface water, the lake and the upstream water are considered to be mainly recharged from glacial/snowpack melting rather than other precipitations.

To assess the magnitude of uncertainty of the ²²²Rn transient model, the sensitivity of estimated LGD to changes in other variables is examined. A sensitivity coefficient f is proposed to evaluate this uncertainty according to Langston et al. (2013):

$$\begin{array}{}\text{(3)}& f=(\mathrm{\Delta }{F}_{\mathrm{LGD}}/F_{\mathrm{LGD}})/(\mathrm{\Delta }{y}_i/y_i),\end{array}$$

where ΔF_LGD is the amount of change in F_LGD from the original value. Δy_i is the amount of change in the other variable of y_i from the original value. Thus, higher f indicates a large uncertainty of the final LGD estimate. The uncertainty mainly stems from ²²²Rn measurements in different water endmembers, the atmospheric loss, and the water level record. The uncertainties of ²²²Rn measurement are about 10 % and 15 %–20 % in groundwater and lake water endmembers, respectively. The uncertainty of atmospheric loss is derived from uncertainty of ²²²Rn in lake water (with an uncertainty of 15 %–20 %), temperature (with an uncertainty of ≈5 %), and wind speed (with an uncertainty of ≈5 %). Thus, the final LGD estimate has an integrated uncertainty of 35 %–40 %.

4.3 Hydrologic partitioning

Compared to the groundwater labelled radionuclide of ²²²Rn, stable ¹⁸O ∕ ²H isotopes are advantageous in the investigation of evaporation processes due to their fractionations from water to vapour and have been widely used to investigate the hydrologic cycle of lakes in various environments (Stets et al., 2010; Gat, 1995; Gonfiantini, 1986; Gibson et al., 1993). With the field data of stable isotopic composition and Cl⁻ concentrations in different water endmembers, groundwater input, surface water input, lake water outflow and infiltration, and evaporation can be partitioned by coupling the stable isotopic mass balance model with the Cl⁻ mass balance model (Fig. 6b).

The model, consisting of the budgets of stable isotopes and Cl⁻, and water masses for the epilimnion, is used to quantify riverine inflow, lake water outflow and infiltration, and evaporation (LaBaugh et al., 1995, 1997; Gibson et al., 2016). The model is valid under the following assumptions: (1) constant density of water; (2) no long-term storage change in the reservoir; (3) well mixed for the epilimnion (Gibson, 2002; Gibson et al., 2016; Gibson and Edwards, 2002; LaBaugh et al., 1997). The above assumptions are reasonably tenable during the short monitoring period. The model can be fully expressed as

$$\begin{array}{ll}\text{(4)}& {\displaystyle }& {\displaystyle }{F}_{\mathrm{in}}+F_{\mathrm{LGD}}+F_{\mathrm{p}}=F_{\mathrm{E}}+F_{\mathrm{out}},{\displaystyle }& {\displaystyle }{F}_{\mathrm{in}}\times {\mathit{\delta }}_{\mathrm{in}}+F_{\mathrm{LGD}}\times {\mathit{\delta }}_{\mathrm{gw}}+F_{\mathrm{p}}\times {\mathit{\delta }}_{\mathrm{p}}\\ \text{(5)}& {\displaystyle }& {\displaystyle }=F_{\mathrm{E}}\times {\mathit{\delta }}_{\mathrm{E}}+F_{\mathrm{out}}\times {\mathit{\delta }}_{\mathrm{L}},{\displaystyle }& {\displaystyle }{F}_{\mathrm{in}}\times [{\mathrm{Cl}}^{-}{]}_{\mathrm{in}}+F_{\mathrm{LGD}}\times [{\mathrm{Cl}}^{-}{]}_{\mathrm{gw}}+F_{\mathrm{p}}\times [{\mathrm{Cl}}^{-}{]}_{\mathrm{p}}\\ \text{(6)}& {\displaystyle }& {\displaystyle }=F_{\mathrm{out}}\times [{\mathrm{Cl}}^{-}{]}_{\mathrm{L}},\end{array}$$

where F_in (mm d⁻¹) is the surface water inflow to the lake; F_gw (mm d⁻¹) is the LGD rate. F_p (mm d⁻¹) is the mean daily rainfall rate during the sampling period. F_E (mm d⁻¹) is the lake evaporation. F_out (mm d⁻¹) is the lake water outflow via runoff and infiltration into the lakeshore aquifer. δ_in, δ_gw, δ_E, and δ_p are the isotopic compositions of surface water inflow, LGD, and evaporative flux, respectively. The values of δ_in, δ_gw, and δ_p are obtained from field data and the composition of δ_E are calculated as shown in the Supplement. [Cl⁻]_in, [Cl⁻]_gw, [Cl⁻]_L, and [Cl⁻]_p are the chloride concentrations in the inflow water, porewater, lake water, and precipitation, respectively.

Table 2 Input parameters for the three-endmember model of Ximen Co Lake.

The components of the mass balance model can be obtained from the field data of isotopic composition and Cl⁻ concentrations in different water endmembers. The average ¹⁸O composition −13.37 ‰ of riverine inflow water is taken as the value of the input parameter δ_in. δ¹⁸O and δ²H in the groundwater endmember and lake water endmember are calculated to be −12.41 ‰ and −87.18 ‰, respectively. δ¹⁸O and δ²H in the rainfall are measured to be −5.47 ‰ and −24.98 ‰, respectively. With the measured values of δ_L, h, and δ_in, and the estimated ε and δ_a, the isotopic composition of δ_E is calculated to be −35.11 ‰, which is in line with the results of alpine and arctic lakes elsewhere (Gibson, 2002; Gibson et al., 2016; Gibson and Edwards, 2002). The values of [Cl⁻]_in, [Cl⁻]_gw, and [Cl⁻]_L are calculated to be 0.91 mg L⁻¹, 1.48 mg L⁻¹, and 1.02 mg L⁻¹, respectively. All the parameters used in the model are shown in Table 2.

According to Eqs. (4)–(6), the uncertainties of calculations of F_in, F_out, and E are mainly derived from the uncertainty of F_LGD and the compositions of Cl⁻, ΔD, and δ¹⁸O in different water endmembers as suggested in previous studies (Genereux, 1998; Klaus and McDonnell, 2013). The compositions of Cl⁻, ΔD, and δ¹⁸O in surface water–groundwater endmembers have an uncertainty of 5 %. The uncertainty of δ_E is reasonably assumed to be ≈20 %. Thus, considering the uncertainty propagation of all the above parameters, the uncertainties of F_in, F_out, and E would be scaled up to 70 %–80 % of the final estimates.

Figure 8 The hydrologic partition of the proglacial lake of Ximen Co (a) and the budgets of DIN (b) and DIP (c).

4.4 The hydrologic partitioning of the glacial lake

Based on the three-endmember model of ¹⁸O and Cl⁻, the riverine inflow rate was calculated to be 135.6±119.0 mm d⁻¹, and the lake outflow rate is estimated to be 141.5±132.4 mm d⁻¹; the evaporation rate is calculated to be 5.2±4.7 mm d⁻¹. The summary of the hydrologic partitioning of the lake is shown in Fig. 8a. Generally, the proglacial lake is mostly recharged by the riverine inflow from the snowpack or the glacier melting. The groundwater discharge contributes about only 7.0 % of the total water input to the lake, indicating groundwater input does not dominate water input to the proglacial lake. The recent review of LGD rate by Rosenberry et al. (2015) suggests that the median of LGD rate in the literatures is 7.4 mm d⁻¹ (0.05 to 133 mm d⁻¹), which is about 2∕3 of the LGD rate in this study. This difference may be due to the hydrogeological setting of the lakeshore aquifer. This aquifer is formed by grey loam, clayey soil, and sand (Lehmkuhl, 1998; Schlutz and Lehmkuhl, 2009), which has a relatively high permeability.

Previous studies have indicated that groundwater forms a key component of proglacial hydrology (Levy et al., 2015). However, there have been limited quantitative studies of groundwater contribution to the hydrologic budget of proglacial areas. This study further summarizes the groundwater discharge studies over the glacial forefield areas. Based on long-term hydrological and climatological parameter monitoring on Nam Co Lake in the QTP, Zhou et al. (2013) estimated the LGD to be 5–8 mm d⁻¹, which is comparable to the surface runoff input and LGD of this study. Brown et al. (2006) investigated the headwater streams in the proglacial areas of Taillon Glacier in France and found that groundwater contributes 6 %–10 % of the stream water immediately downwards of the glacier. Using the water mass balance model, Hood et al. (2006) show that groundwater inflow is substantial in the hydrologic partitioning of the proglacial Lake O'Hara in front of Opabin Glacier in Canada and comprised 30 %–74 % of the total inflow. Roy and Hayashi (2008) studied the proglacial lakes of Hungabee Lake and Opabin Lake at the glacier forefield of Opabin Glacier and found that the groundwater component is a predominant water source of the lakes and consisted of 35 %–39 % of the total water input of the lakes. Langston et al. (2013) further investigated a tarn immediately in front of Opabin Glacier and indicated the tarn is predominantly controlled by groundwater inflow/outflow, which consisted of 50 %–100 % of the total tarn volume. Magnusson et al. (2014) studied the streams in the glacier forefield of the Damma Glacier, Switzerland, and revealed that groundwater contributed only 1 %–8 % of the total surface runoff. Groundwater contribution in this study is similar to those obtained from the mountainous proglacial areas in Europe, but much lower than those obtained in the proglacial areas of polar regions. It is concluded that proglacial lakes/streams in front of mountainous glaciers are mainly recharged by surface runoff from glacier/snowpack melting. This might be due to well-developed stream networks and limited deep groundwater flow (Einarsdottir et al., 2017; Brown et al., 2006; Magnusson et al., 2014). However, proglacial tarns and lakes in the polar areas are predominantly controlled by groundwater discharge, due to less connectivity of surface runoff and high shallow and deep groundwater connectivity (Langston et al., 2013; Hood et al., 2006; Roy and Hayashi, 2008).

The evaporation constitutes a relatively small ratio (≈3.5 %) of total water losses. The annual evaporation rate was recorded to be 1429.8 mm (equivalent to 3.92 mm d⁻¹) in 2015 by the Jiuzhi weather station, lower than the obtained evaporation in this study. This may be due to much higher evaporation in August during the monitoring period. The estimation of evaporation in this study generally represents the upper limit of the lake, as the sampling campaign was conducted during the summer time when the highest evaporation might occur. The lake surface evaporation derived from the pan evaporation in the QTP ranges from ∼700 mm yr⁻¹ in the eastern QTP to over 1400 mm yr⁻¹ in the interior lakes of the QTP (Zhang et al., 2007; Ma et al., 2015; Yang et al., 2014). The evaporation of this study is rather in line with the previous evaporation observation in the eastern QTP, stressing the tenability of evaporation in this study. The runoff input is the predominant recharge component (>90 %) compared to other components, with an area normalized value comparable to previous studies of runoff input in other glacial melting dominant lakes in the QTP (Zhou et al., 2013; Zhang et al., 2011; Biskop et al., 2016). The runoff input and the lake evaporation of the study area, however, are subject to high daily, seasonal, and inter-annual variability as indicated by previous studies in the QTP (Zhou et al., 2013; Lei et al., 2017; Ma et al., 2015; Lazhu et al., 2016). Therefore, further investigations of long-term and high-resolution climatological and isotopic data are required to provide precise constraints on hydrologic partitioning of the lakes in the QTP.

4.5 LGD-derived nutrient loadings, nutrient budget, and ecological implications

Compared to extensive studies of SGD-derived nutrient loadings in the past decade (Luo and Jiao, 2016; Slomp and Van Cappellen, 2004), studies of LGD-derived nutrient loadings have received limited attention, even given the fact that groundwater in lakeshore aquifers is usually concentrated in nutrients (Lewandowski et al., 2015; Rosenberry et al., 2015). Even fewer studies focus on chemical budgets in the proglacial lakes which are often difficult to access for sampling. Groundwater-borne DIN and DIP across the sediment–water interface in this study are determined with an equation coupling the advective, or LGD-derived, and diffusive solute transport (Lerman et al., 1995; Hagerthey and Kerfoot, 1998):

$$\begin{array}{}\text{(7)}& F_j=-n{D}_j^{\mathrm{m}}{\displaystyle \frac{\mathrm{d}{C}_j}{\mathrm{d}x}}+v_{\mathrm{gw}}{C}_j,\end{array}$$

where $-n{D}_j^{\mathrm{m}}\frac{\mathrm{d}{C}_j}{\mathrm{d}x}$ is the diffusion input and v_gwC_j is the LGD-derived fluxes; F_j (µM m⁻² d⁻¹) is the mol flux of nutrient species j (representing DIN or DIP). n is the sediment porosity. $D_j^{\mathrm{m}}$ is the molecular diffusion coefficient of nutrient species j, which is given to be $\mathrm{4.8}\times {\mathrm{10}}^{-\mathrm{5}}$ m² d⁻¹ for DIP (Quigley and Robbins, 1986) and $\mathrm{8.8}\times {\mathrm{10}}^{-\mathrm{5}}$ m² d⁻¹ for DIN (Li and Gregory, 1974), respectively. C_j (µM) is the concentration of nutrient species j. x (m) is the sampling depth. v_gw is LGD rate estimated by the ²²²Rn mass balance model and has a value of 10.3±8.2 mm d⁻¹. $\frac{\mathrm{d}{C}_j}{\mathrm{d}x}$ is the concentration gradient of nutrient species j across the water–sediment interface.

Substituting the constants and the field data of DIN and DIP into Eq. (6), LGD-derived nutrient loadings are calculated to be 954.3 and 3.2 µmol m⁻² d⁻¹ for DIN and DIP, respectively. Riverine inflow brings 1195.0 µmol m⁻² d⁻¹ DIN and 52.9 µmol m⁻² d⁻¹ DIP into the lake. Lake water outflow-derived nutrient loss is estimated to be 1439.9 and 54.7 µmol m⁻² d⁻¹ for DIN and DIP, respectively. Nutrients in the lake can also be sourced from atmospheric deposit (mostly in the form of precipitation). With the nutrient concentrations in the rainwater during the monitoring period, the wet deposit is calculated to be 76 and 2.5 µmol m⁻² d⁻¹ for DIN and DIP, respectively. The loadings of DIN to the lakes are mainly from surface runoff and LGD, which comprised 42.9 % and 53.7 % of the total DIN loadings. Groundwater-derived DIP input, however, constitutes only 6.3 % of the total DIP inputs to the lake, indicating groundwater-borne DIP is less contributive to the nutrient budget of the lake compared to DIN. Very recent studies on polar regions have indicated that the glacier/snowpack water is the main N source to the proglacial lakes (Anderson et al., 2013; Dubnick et al., 2017). However, they do not consider the contribution of groundwater-borne N, in spite of the high groundwater connectivity in the proglacial areas (Roy and Hayashi, 2008). This study stresses that groundwater-borne DIN could be comparable to the surface runoff-derived DIN.

Based on nutrient results, the lake is considered to be an oligotrophic lake, similar to other glacier melting dominant lakes in the QTP (Mitamura et al., 2003; Liu et al., 2011). Phytoplankton are good dissolved organic phosphate (DOP) recyclers and will overcome inorganic P limitation though DOP cycling in most template lakes (Hudson et al., 2000). However, this may be not applicable for the glacial melting water and the periglacial/proglacial lake water. Previous studies show that phosphate nutrients are dominated by DIP and particulate phosphate, and the DOP contributes less than 10 % of the dissolved phosphate (Cole et al., 1998; Hawkings et al., 2016; Hodson, 2007). Thus, DOP recycling is not likely to lower the N:P ratio under these conditions. Thus, the primary production (PP) is considered to be controlled by the DIP loadings. The sum of DIN and DIP inputs minus the calculated DIN and DIP outputs leads to surpluses of 785.4 and 3.9 µmol m⁻² d⁻¹ for DIN and DIP, respectively. The surpluses are expected to be consumed by the phytoplankton and converted into the PP under phosphate-limited conditions. As primary producers in the fresh lacustrine system consume the nutrient under variant N:P ratios (7.1 to 44.2, mean: 22.9) (Downing and McCauley, 1992), the biological uptake of DIN is roughly estimated to be 89.3 µM m² d⁻¹. Therefore, the nutrient budgets for DIN and DIP can be finally conceptualized in Fig. 8b and c.

4.6 Implications, prospects, and limitations

Mountainous proglacial lakes are readily developed in glacier forefields of the QTP and other high mountainous glacial areas such as the European Alps and Pamir in central Asia (Heckmann et al., 2016). The proglacial lakes are always trapping systems of sediment and sinks for water and chemicals originating from glacier/snowpack melting and groundwater. Analogously to cosmogenic isotopes such as ¹⁰Be serving as a tool to quantify the sediment sources, approaches integrating ²²²Rn and stable isotopes provide both qualitative and quantitative evaluations of groundwater contributions and hydrologic partitioning in these remote and untapped lacustrine systems. Thus, it is expected that the multiple aqueous isotopes will be considered to be effective tools to investigate the LGD and hydrologic partitioning in other proglacial lakes. This study is mainly limited by the relatively short sampling and monitoring period. As a special hydrologic regime, the lakeshore aquifers of the proglacial lakes are experiencing frozen–unfrozen transition seasonally, and the dominant recharge of glacial melting could fluctuate significantly due to air temperature variation. Therefore, future groundwater and hydrological studies can be extended to longtime sampling and monitoring of stable isotopes and ²²²Rn in different water endmembers to reveal the seasonal hydrological and hydrogeological dynamics and their impacts on local biogeochemical cycles and ecological systems. Special concerns would be placed on how surface–groundwater interactions and the associated biogeochemical processes respond to the seasonal frozen ground variations and glacier/snowpack melting intensity.

5 Conclusions

A ²²²Rn continuous monitoring is conducted at Ximen Co Lake, a proglacial lake located in the eastern QTP. A dynamic ²²²Rn mass balance model constrained by radium mass balance and water level fluctuation is used to quantify temporal distribution of the LGD of the lake. The obtained LGD over the monitoring time ranges from −23.7 to 80.9 mm d⁻¹, with an average of 10.3±8.2 mm d⁻¹. Thereafter, a three-endmember model consisting of the budgets of water, stable isotopes, and Cl⁻ is used to depict the hydrologic partitioning of the lake. Riverine inflow, lake water outflow via surface runoff, and surface evaporation are estimated to be 135.6, 141.5, and 5.2 mm d⁻¹, respectively. LGD-derived nutrient loading is estimated to be 785.4 and 3.2 µmol m⁻² d⁻¹ for DIN and DIP, respectively. This study also implies that LGD constitutes a relatively small portion of the proglacial hydrologic partitioning; however, it delivers nearly half of the nutrient loadings to the proglacial lake.

This study presents the first attempt to quantify LGD and the associated nutrient loadings to the proglacial lake of the QTP. To our knowledge, there is almost no study on the groundwater–lake water interaction in the high-altitude proglacial lakes in the QTP. This study demonstrates that the ²²²Rn-based approach can be used to investigate the groundwater dynamics in the high-altitude proglacial lakes. The method is instructional to similar studies in other proglacial lakes in the QTP and elsewhere. For a comprehensive understanding of the hydrological and biogeochemical dynamics in the QTP, interdisciplinary and multi-approach integrated studies are greatly needed. Of particular importance are the lake hydrology and groundwater–surface water interaction studies based on multiple approaches such as remote sensing products, long-term and high-resolution observation of climatological parameters, and isotopic data.