**1. Introduction**

Precipitation comes from the antecedent atmospheric moisture, lateral advection, and local evapotranspiration [1,2]. Research has shown that lateral advection is the major moisture source of precipitation globally, with *fTr* (proportion of plant transpiration water vapor in precipitation) and *fEv* (proportion of surface evaporation water vapor in precipitation) being the second and third. However, there is considerable spatiotemporal variation among the three types of moisture in the world [3,4], which profoundly affects the global and local water cycle [5]. Stable hydrogen and oxygen isotopes can play an important role in the quantitative research of the hydrologic cycle [6,7]. The contribution of different sources of moisture to precipitation has been a hot topic in isotope hydrology. In general, it is feasible to observe and calculate the antecedent atmospheric water vapor and lateral advection directly. The difficulty is to determine the contribution of recycled moisture. Recycled moisture mainly comes from *fTr* (proportion of plant transpiration water vapor in precipitation) and *fEv* (proportion of surface evaporation water vapor in

**Citation:** Zhang, Z.; Zhu, G.; Pan, H.; Sun, Z.; Sang, L.; Liu, Y. Quantifying Recycled Moisture in Precipitation in Qilian Mountains. *Sustainability* **2021**, *13*, 12943. https://doi.org/10.3390/ su132312943

Academic Editors: Alban Kuriqi and Luis Garrote

Received: 15 October 2021 Accepted: 20 November 2021 Published: 23 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

precipitation) [8–11]. The linear mixing model for isotopes is an effective method to study the contribution of recycling moisture in different regions. Linear mixing model has been applied around the world [12], such as the Great Lakes region in North America [13], the Slave River Delta of Canada [14], the Amazon Basin of South America [15], the Nam Co Basin in the Qinghai-Tibetan Plateau [16], the Lake Shorty in Madagascar [17] and Lake Kasumiguara [18].

Most research on the moisture cycle focuses on the source of advected water vapor and its transport in different regions [8,19]. Recently, more studies have considered the transport and conversion of basin recycled moisture [4,6]. Research showed that the contribution of recycled moisture to precipitation varies greatly in time and space. In arid inland river basins, the contribution of recycled moisture is less than 30% [6,20], In some small oases, the recycled moisture contribution is less than 5% [21,22]. However, the contribution of recycled moisture can reach 62% in the Tibetan Plateau, but is much lower in winter than that in summer [6,23]. Some studies have focused on the causes of the spatial and temporal differences in recycled moisture, such as relative air humidity [19], soil water content [24], land use types and land cover changes [25].

Many studies have estimated the contribution rate of recycled moisture and its influencing factors in the watershed or regional scale [4,6,21,22,26]. Due to the obvious difference between climate and environment in the mountainous-oasis-desert regions in arid areas, studies based on basins or regional scales tend to overlook the uniqueness of recycling moisture in mountain areas. The upper Shiyang River in the northeast of Qilian Mountain is the transition zone between the Qinghai-Tibet Plateau and the arid zone. Precipitation affects the development of the oasis and desert in the middle and lower reaches. Clarifying the characteristics and influencing factors of upstream water resources changes will contribute to a reasonable solution to the demands of the middle and lower reaches of the water resources. Therefore, this study used the stable isotope data of precipitation, soil water, plant water, surface water, and groundwater from 2016 to 2018 in the upper Shiyang River in the eastern part of the Qilian Mountains to calculate the proportion of plant transpiration water vapor (fTr), surface evaporation water vapor (fEv) and advection water vapor (fAdv) in precipitation. We try to explore the source of moisture in precipitation in the mountain areas of the arid inland river basin, and to reveal the characteristics of the mountain water cycle, reasonable assessment of regional water resources.

#### **2. Study Area and Observation Network**

### *2.1. Study Area*

The study area is located in the Xiying River (XYR) basin in the Qilian Mountains (Figure 1). XYR is the main tributary of the Shiyang River. It originates below the glaciers in the northern slopes of the Qilian Mountains and eventually disappears into the desert. The elevation of the XYR basin ranges from 1510 to 4874 m above sea level. The average annual temperature is 6.3 ◦C and the average annual precipitation ranges from 200 to 700 mm [8,27]. The upper reaches of the Shiyang River basin are located in the East Asian monsoon crisscross zone, controlled by the East Asian monsoon and plateau monsoon. Cold high pressure appears on the plateau in winter, and the air flows from the plateau to the surroundings. In summer, hot low air pressure appears on the plateau, and the air flows from all directions to the plateau [6], which is a typical continental alpine climate. Vegetation is mainly distributed in areas between 2000 m and 3600 m above sea level; the basin is affected by multiple sources of moisture [6,27].

**Figure 1.** The distribution of all sampling sites. A: Xiyingwugou (Foothill), B: Huajian (Arbor belt), C: Ningchan (Shrub belt), D: Lenglongling (River source).

#### *2.2. Observation Network, Sampling, and Analysis*

From April 2016 to October 2018, four sampling sites were established in the XYR Basin (Figure 1). A total of 867 samples were collected from the upstream mountain areas (Table 1). Two-hundred forty-five precipitation samples at sampling sites were collected immediately after the end of each precipitation process using a rain gauge. The precipitation samples were put into a 50 mL polyethylene sampling bottle. The bottle cap was tightened, and the bottle mouth was sealed with a sealing film and stored in cold storage until analysis. Surface water samples were sealed and stored in cold storage after each collection. Meanwhile, automatic meteorological observation instruments recorded meteorological elements such as temperature, precipitation, relative humidity and atmospheric pressure [6]. All samples were analyzed for δ2H and δ18O using liquid water. For runoff, the bottle was placed under the water surface until the container was filled. A total of 82 runoff samples were collected at different elevations (Table 1). For soil water, soil samples were collected from 10 cm to 100 cm, using soil spirals every 10 cm. All samples were placed in 100 mL polyethylene vials. A total of 450 soil samples were collected (Table 1). For plant water, typical vegetation at each sampling site was selected, and scissors were used to cut the xylem or branches of the same vegetation. A total of 90 plant samples were collected (Table 1). All samples were stored in a mobile freezer (−5 ◦C) and transferred to a freezer laboratory (−15 ◦C) within a week after collection.

The samples were melted at room temperature (20~25 ◦C) before further analyses. The plant and soil water were extracted using a cryogenic vacuum distillation apparatus (LI-2100, LICA United Technology Limited, China). All samples were analyzed using a liquid water isotope analyzer DLT-100 (Los Gatos Research, Inc.) in the Stable Isotope Laboratory, College of Geography and Environmental Science, Northwest Normal University. The isotope ratios of samples are expressed as parts per mil (‰) relative to the Vienna Standard Mean Ocean Water (V-SMOW) using δ notation: δ(‰) = (RS/RV-SMOW−1) × 1000, where δs is the isotope ratio of the samples relative to V-SMOW and Rs is the ratio of D/H or 18O/16O in the samples. The precision of the measurements was ±0.6% for <sup>δ</sup>2H and ±0.2% for δ18O, respectively.


**Table 1.** Basic information for each sampling site.

Meteorological data were obtained from four weather stations of the XYR basin. These weather stations record meteorological data such as temperature, precipitation, wind speed, evaporation, relative humidity, and soil moisture every 15 min.

#### **3. Method**

#### *3.1. Three-Component Mixing Model*

Moisture from precipitation is derived from local evapotranspiration and advection [28,29]. The linear mixing model can be used to calculate the contribution rate of each water source. The three-component mixing model based on δ18O and δ2H values can be expressed as follows.

$$
\delta^{18}O\_{\mathrm{Pr}} = \delta^{18}O\_{\mathrm{Tr}} \times f\_{\mathrm{Tr}} + \delta^{18}O\_{\mathrm{E}\upsilon} \times f\_{\mathrm{E}\upsilon} + \delta^{18}O\_{\mathrm{Ad}\upsilon} \times f\_{\mathrm{Ad}\upsilon} \tag{1}
$$

$$
\delta D\_{\rm Pv} = \delta D\_{\rm Tr} \times f\_{\rm Tr} + \delta D\_{\rm Ev} \times f\_{\rm Ev} + \delta D\_{\rm Adv} \times f\_{\rm Adv} \tag{2}
$$

$$1 = f\_{Tr} + f\_{Ev} + f\_{Adv} \tag{3}$$

In formula (1)–(3), the contribution rate of each moisture source is expressed in terms of *f* while the subscript indicates the source of moisture. Subscript *Pv* indicates the local atmospheric moisture. Subscripts *Tr, Ev, Adv*, indicate moisture from transpiration, evaporation, and advection, respectively.

The contributions of *fTr*, *fEv*, and *fAdv* can be calculated as follows.

*fTr* <sup>=</sup> *<sup>δ</sup>*<sup>18</sup>*OPv* <sup>×</sup> *<sup>δ</sup>DEv* <sup>−</sup> *<sup>δ</sup>*<sup>18</sup>*OPv* <sup>×</sup> *<sup>δ</sup>DAdv* <sup>+</sup> *<sup>δ</sup>*<sup>18</sup>*OEv* <sup>×</sup> *<sup>δ</sup>DAdv* <sup>−</sup> *<sup>δ</sup>*<sup>18</sup>*OEv* <sup>×</sup> *<sup>δ</sup>DPv* <sup>+</sup> *<sup>δ</sup>*<sup>18</sup>*OAdv* <sup>×</sup> *<sup>δ</sup>DPv* <sup>−</sup> *<sup>δ</sup>*<sup>18</sup>*OAdv* <sup>×</sup> *<sup>δ</sup>DEv δ*<sup>18</sup>*OTr* × *δDEv* − *δ*<sup>18</sup>*OTr* × *δDAdv* + *δ*<sup>18</sup>*OEv* × *δDAdv* − *δ*<sup>18</sup>*OEv* × *δDTr* + *δ*<sup>18</sup>*OAdv* × *δDTr* − *δ*<sup>18</sup>*OAdv* × *δDEv* (4)

$$f\_{\rm Ez} = \frac{\delta^{18}O\_{\rm Pv} \times \delta \rm D\_{\rm Adv} - \delta^{18}O\_{\rm Pv} \times \delta \rm D\_{\rm T} + \delta^{18}O\_{\rm Tv} \times \delta \rm D\_{\rm PV} - \delta^{18}O\_{\rm Tv} \times \delta \rm D\_{\rm Adv} + \delta^{18}O\_{\rm Adv} \times \delta \rm D\_{\rm T} - \delta^{18}O\_{\rm Adv} \times \delta \rm D\_{\rm Ev}}{\delta^{18}O\_{\rm v} + \delta^{18}O\_{\rm Tv} \times \delta \rm D\_{\rm dev} - \delta^{18}O\_{\rm Fv} \times \delta \rm D\_{\rm Fv} + \delta^{18}O\_{\rm Adv} \times \delta \rm D\_{\rm Ev} + \delta^{18}O\_{\rm Adv} \times \delta \rm D\_{\rm Fv}}} \tag{5}$$

$$f\_{\rm Adv} = \frac{\delta^{18} \mathcal{O}\_{\rm Pv} \times \delta \mathcal{D}\_{\rm T} - \delta^{18} \mathcal{O}\_{\rm Pv} \times \delta \mathcal{D}\_{\rm Ev} + \delta^{18} \mathcal{O}\_{\rm Tv} \times \delta \mathcal{D}\_{\rm Ev} - \delta^{18} \mathcal{O}\_{\rm Tv} \times \delta \mathcal{D}\_{\rm Pv} + \delta^{18} \mathcal{O}\_{\rm Ev} \times \delta \mathcal{D}\_{\rm Pv} - \delta^{18} \mathcal{O}\_{\rm Ev} \times \delta \mathcal{D}\_{\rm Tv}}{\delta^{18} \mathcal{O}\_{\rm Ev} + \delta^{18} \mathcal{O}\_{\rm Ev} \times \delta \mathcal{D}\_{\rm Adv} - \delta^{18} \mathcal{O}\_{\rm Ev} \times \delta \mathcal{D}\_{\rm Tv} + \delta^{18} \mathcal{O}\_{\rm Adv} \times \delta \mathcal{D}\_{\rm Ev} - \delta^{18} \mathcal{O}\_{\rm Adv} \times \delta \mathcal{D}\_{\rm Ev}}} \tag{6}$$

The *δ*<sup>18</sup>*Opv* and *δDpv* value can also be derived from the value of the stable isotope from local precipitation and the equilibrium fractionation factor:

$$
\delta \, ^{18}O\_{\mathcal{P}v} \stackrel{\simeq}{=} \delta \, ^{18}O\_{\mathcal{P}} - 10^3 \times \left( \mathfrak{a}\_{\mathcal{W}-V}^{18} - 1 \right) \tag{7}
$$

$$
\delta D\_{\mathcal{P}v} \cong \delta D\_P - 10^3 \times (\alpha\_{\mathcal{W}-V}^2 - 1) \tag{8}
$$

The specific formula is as follows [30,31]:

$$10^3 Ina\_{W-V}^{18} = 1.137 \times \left(\frac{10^6}{T^2}\right) - 0.4156 \times \left(\frac{10^3}{T}\right) - 2.0667\tag{9}$$

$$10^3 \text{Im}\alpha\_{W-V}^2 = 24.844 \times \left(\frac{10^6}{T^2}\right) - 76.248 \times \left(\frac{10^3}{T}\right) - 52.612\tag{10}$$

The *δ*<sup>18</sup>*OEv* and *δDev* can be expressed by the value of the stable isotope from advection (*δ*<sup>18</sup>*OAdv* and *δDAdv*) and local the surface water (*δ*18*Os* or *δDs*), mean relative humidity (*h*), the sum of equilibrium (*εeq*) and kinetic (Δ*ε*):

$$
\delta^{18}O\_{Ev} = \frac{\delta^{18}O\_S - h \times \delta^{18}O\_{Adv} - \varepsilon^{18}}{1 - h} \tag{11}
$$

$$
\delta D\_{Ev} = \frac{\delta D\_S - h \times \delta D\_{Adv} - \varepsilon^2}{1 - h} \tag{12}
$$

$$
\varepsilon^{18} = \varepsilon\_{eq}^{18} + \Delta \varepsilon^{18} \tag{13}
$$

$$
\varepsilon^2 = \varepsilon\_{eq}^2 + \Delta \varepsilon^2 \tag{14}
$$

$$
\varepsilon\_{c\eta}^{18} = 1000 \times (1 - \frac{1}{a\_{W-V}^{18}}) \tag{15}
$$

$$
\varepsilon\_{c\eta}^2 = 1000 \times (1 - \frac{1}{\mathfrak{a}\_{W-V}^2}) \tag{16}
$$

$$
\Delta \varepsilon^{18} = 14.2 \times (1 - h) \tag{17}
$$

$$
\Delta \varepsilon^2 = 12.5 \times (1 - h) \tag{18}
$$

Based on previous studies on *δ*<sup>18</sup>*OAdv* and *δDAdv*, we used the following formula to investigate the characteristics in the XYR Basin.

$$
\delta \, ^{18}O\_{Adv} \cong \delta \, ^{18}O\_{pv} + (a\_{w-v}^{18} - 1) \times \ln F \tag{19}
$$

$$
\delta D\_{Adv} \cong \delta D\_{pv} + (a\_{w-v}^2 - 1) \times \ln \text{F} \tag{20}
$$

In the above formulas, the *F* indicates the ratio between the initial and the final vapor, which is estimated by the precipitable water amount in the two sites. Previous research in the Tianshan Mountains concluded that precipitable water correlates with moisture pressure (*c* = 1.657e, where *c* indicates atmospheric moisture content in mm and e indicates surface moisture pressure in hPa, r<sup>2</sup> = 0.94) [32]. Hence, the surface moisture pressure ratio between the two sampling sites is equal to the value of *F*. Since the isotope ratios in precipitating vapor at the sampling sites are much more depleted than the upwind station, the Rayleigh distillation equation was also applied. In this research, formula 6 was used to calculate the stable isotope ratios in advection vapor. If there is no significant depletion of isotopes between the sampling site and the upwind station, then the stable isotope ratios in advection vapor at the sampling site are considered to be the same as that in the precipitating vapor at the upwind station.

Because there is no fractionation during the transport of moisture from *fTr* to the atmosphere [33], the value of *δ*<sup>18</sup>*OTr* and *δDTr* is the same as that of local water used by plants (*δ*18*Ow* and *δDw*). This research calculated the average value of *δ*18*O* (*δ*2*H*) from the soil surface to a depth of 40cm below and xylem, to obtain *δ*<sup>18</sup>*OTr* and *δDTr*. The results for each sampling site in this study area are shown in Table 2.


**Table 2.** The data needed to calculate the recycled water vapor contribution rate and the calculation results.

— It cannot be calculated because the data are missing.

#### *3.2. Hysplit Model and the Upper Wind Direction*

We applied the HYSPLIT model to simulate the moisture sources in the Qilan Mountains [34–37]. We found that westerly winds, southeast monsoons, and plateau monsoons all affect the Qilian Mountains in summer. In winter, westerly winds mainly affect the Qilian Mountains.

According to the clustering of air mass in different seasons, the air mass gathered at the northern foot of Qilian Mountain and then moved from a low elevation to a high elevation along the valley. Therefore, sampling site A was used as an upwind station for spring, summer, and autumn. In winter, the study area was dominated by a westerly wind, and Urumqi and Hotan (GNIP) in Central Asia were regarded as the upwind direction stations.
