**Preliminary Characterization of Underground Hydrological Processes under Multiple Rainfall Conditions and Rocky Desertification Degrees in Karst Regions of Southwest China**

**Guijing Li 1,2,**†**, Matteo Rubinato 3, Long Wan 1,2,**†**, Bin Wu 2, Jiufu Luo 1,2, Jianmei Fang <sup>2</sup> and Jinxing Zhou 1,2,\***


Received: 22 January 2020; Accepted: 18 February 2020; Published: 21 February 2020

**Abstract:** Karst regions are widely distributed in Southwest China and due to the complexity of their geologic structure, it is very challenging to collect data useful to provide a better understanding of surface, underground and fissure flows, needed to calibrate and validate numerical models. Without characterizing these features, it is very problematic to fully establish rainfall–runoff processes associated with soil loss in karst landscapes. Water infiltrated rapidly to the underground in rocky desertification areas. To fill this gap, this experimental work was completed to preliminarily determine the output characteristics of subsurface and underground fissure flows and their relationships with rainfall intensities (30 mm h<sup>−</sup>1, 60 mm h−<sup>1</sup> and 90 mm h<sup>−</sup>1) and bedrock degrees (30%, 40% and 50%), as well as the role of underground fissure flow in the near-surface rainfall–runoff process. Results indicated that under light rainfall conditions (30 mm h<sup>−</sup>1), the hydrological processes observed were typical of Dunne overland flows; however, under moderate (60 mm h<sup>−</sup>1) and high rainfall conditions (90 mm h−1), hydrological processes were typical of Horton overland flows. Furthermore, results confirmed that the generation of underground runoff for moderate rocky desertification (MRD) and severe rocky desertification (SRD) happened 18.18% and 45.45% later than the timing recorded for the light rocky desertification (LRD) scenario. Additionally, results established that the maximum rate of underground runoff increased with the increase of bedrock degrees and the amount of cumulative underground runoff measured under different rocky desertification was SRD > MRD > LRD. In terms of flow characterization, for the LRD configuration under light rainfall intensity the underground runoff was mainly associated with soil water, which was accounting for about 85%–95%. However, under moderate and high rainfall intensities, the underground flow was mainly generated from fissure flow.

**Keywords:** rock–soil contact area; fissure flow; karst rocky desertification; runoff; rainfall simulation

#### **1. Introduction**

Karst landscapes represent a crucial feature of the earth's geodiversity [1], and they account for 12% of the total land area in the world [2]. Southwest China is highly characterized by karst areas. In fact, the degraded vegetation and widely exposed limestone in Southwest China induced

severe rocky desertification, generating a fragile ecosystem that represent a severe environmental and social issue [3–5]. Furthermore, human activities and forest clearance in western Ireland caused severe soil loss [6]. To date, several studies [7–17] have been conducted to provide a solution to the management of these areas, and water has been identified as the major drive of the hydrological processes investigated [18–21]. In actual fact, due to the nature of these environmental layers typical of karst areas, rainfall is largely lost through underground fissures, karst caves and pipelines [22,23], causing major concerns to local authorities on how to prevent droughts and better manage water distribution within the region.

Hydrological processes generated by surface runoff in karst areas have been identified to be linked to land degradation [24–26] or vegetation restoration [27], and in recent years, thanks to the improvement of existing technology and a wider understanding of local phenomena, most studies have focused on the characterization of groundwater in epikarst zones [28–30].

Multiple experimental tests and field investigations of rock and soil structures have been conducted as well as numerical simulations of underground water flow induced by rainfall events [31–33]. However, these studies did not carry out an accurate experimental campaign incorporating the soil–rock interaction, typical of the natural situation; hence, all results obtained are limited to artificial systems that do not completely reflect the real case scenario. Fu et al. [20] made interesting field observations on the dynamic change of the water balance and the underground water flow component by digging trenches at the subsurface boundary of karst developed over dolomite, characterized by a flat depression surrounded by overlapping hills and ridges [20]. Furthermore, karsts systems without soil are different and have a dissimilar hydrogeological behavior because the soil is the main source for CO2 production, and the surface of the carbonate rocks directly contacting the overlying soil suffers the most intense chemical energies of karstification [34]. These surficial karst processes in turn enhance solutional enlargement of fissures in carbonate rocks and produce an irregular pitted and etched epikarst subsurface [34]. Hence these karsts areas were characterized by a specific rock–soil ratio and permeability of rock–soil fractures, which are important factors influencing underground runoff [35]. If compared with dolomite, carbonate rock and barren soil provide a better material basis for rock desertification, enabling a richer fissure development and causing more serious issues related to soil erosion [36]. Therefore, the karst landforms developed by limestone more typical in karst areas in Southwest China have been the rock chosen to be investigated in this paper in order to provide a better understanding of underground hydrological processes in relation with rainfall events.

Rainfall is in fact another very active dynamic factor in the hydrological cycle [37]. The rainfall amount and the rainfall intensity have a great impact on hydrological dynamics, especially in karst areas [38,39]. Regional climate and hydrological conditions have been demonstrated to affect the soil infiltration and water storage in epikarst zones [40], and the response of infiltration and runoff to rainfall intensity is different for each type of rocky desertification. Hence, it is crucial to further study these underground water behaviors associated with each rainfall intensity to determine how water could be regulated or controlled in karst areas.

To fulfil this gap, this work included the simulation of a unique "dual structure" karst area focusing on the soil–rock integration and the experimental representation of underground runoff with the purposes of (1) providing a better understanding of underground runoff under different rainfall intensities and (2) discussing the influence of rocky desertification on underground runoff. Results obtained by this study have vital theoretical and practical significance to reveal hidden hydrological mechanisms, and can help to reduce water loss supporting healthy and sustainable growth in karst regions.

The paper is organized as follows: Section 2 introduces the description of the experimental setup describing the methods applied and the analysis conducted. Section 3 presents the effects of rainfall intensity and rocky desertification on the dynamics of underground runoff. Section 4 provides a discussion of the results obtained and Section 5 produces a brief summary and concluding remarks of the whole study.

#### **2. Materials and Methods**

#### *2.1. Experimental Setup*

The experimental apparatus used to conduct this research is composed of a pond, previously designed by [40], and a rainfall generator. The pond is made by a reinforced concrete structure with a length of 100 cm, a width of 100 cm and a height of 120 cm. The bottom of the pond is a 20-cm-thick permeable layer, and the permeable layer is filled with gravels (diameter ranging between 5 cm and 10 cm). The remaining 100 cm depth of the pond is used to fill soil and limestone blocks.

The rainfall generator is composed of a water tank which is 110 cm long, 110 cm wide and 25 cm high. A total of 441 holes with a diameter of 1 cm are located at the bottom of the water tank.

The test soil was collected from a 0–100 cm depression located in the typical karst area—Jianshui County, Yunnan province, China (Figure 1). This area is affected by a subtropical monsoon climate, hence characterized by wet and dry seasons. The annual average temperature of this area is 19.7 ◦C, and the annual average precipitation is 828.3 mm. Records states that the rainfall was unevenly distributed throughout the year, happening mainly during the rainy season (from May to October).

**Figure 1.** Geographical location of the study site and bed rock configurations tested (1, 2 and 3), whose characteristics are displayed in Table 1.

The distance between each hole is 5 cm, and the holes are covered with rubber plugs to facilitate the insertion of the needles (inner diameter 0.57 mm) to release the water. The rainfall intensity can be adjusted according to the pressure generated by the different water levels in the water tank. To prevent blockage of the needles, water used to simulate the rainfall is filtered in advance. Furthermore, 9 measuring glasses are located under the sprinkler with a distance between each other of 30 cm, and the volume of the water in the glass can be measured every 10 min to calibrate the water level in the tank and the corresponding rainfall intensity.

Due to the complexity of this phenomenon, we may have simplified the real phenomenon in terms of dimensions or configurations, but we have made reasonable assumptions to clearly distinguish the phenomena observed and how each parameter can affect the flows generated. If considering the main site, there may be issues in terms of operability of tests and too many variables could interfere with the results, making it very challenging to be accurately quantified. In order to accurately simulate the complex karst fissure environment, reducing artificial aspects, standard cuts were made on the limestone rock with the size 10 cm × 10 cm × 10 cm (length × width × height). Soil and rocks were filled in the pond according to the measured soil bulk density in the field from bottom to top. The soil was not sieved before filling and only large soil clumps were dispersed.

In order to reduce and limit the abnormal effects on the edges, the boundaries were artificially compacted manually. Regular sweeping was conducted with a wooden board designed on purpose for this study after each layer was filled. The soil volumetric water content (VWC) at 5 cm, 15 cm, 30 cm, 50 cm and 70 cm of the vertical profile was monitored using an automatic monitoring system because it can be used to calculate the time required for water to penetrate to a certain depth. The probe was buried horizontally perpendicular to the corresponding layer of the soil profile. After completing each layer of soil, probes were buried and calibrated, and the process was repeated for the next layer. Campbell CS616 soil moisture probes (Campbell Scientific, Inc.–Logan, UT, USA) (with an error of ±2.5%) were used for data acquisition, and the acquisition frequency was once every 5 min. The calibration method for the probes can be summarized as follows:


$$M\_{\rm w} = \frac{M\_{\rm a} - M\_{\rm b}}{M\_{\rm b} - M\_{\rm c}} \times 100\% \tag{1}$$

where *Ma* is the total weight of the aluminum box and the wet soil; *Mb* is the total weight of the dried aluminum box and the soil; and *Mc* is the weight of the aluminum box.

• A correlation analysis was then performed between the average value of VWC and the value measured by the sensor at each time. The results obtained showed that there was a linear relationship between these two parameters, and the correlation coefficient was above 0.90.

#### *2.2. Experimental Testing Conditions*

Based on the field investigation of 15 different rocky desertification profiles, and considering existing literature published to date, three configurations were set up with varying parameters: bedrock exposure, vegetation coverage, average soil thickness and rock–soil contact area (Table 1). Each configuration was repeated three times to confirm the repeatability of the tests and reduce experimental errors. According to the grading evaluation criteria for the degree of rocky desertification in the "Technical Regulations for Monitoring Rocky Desertification in the Southwest Karst Area" issued by the National Forestry and Grass Bureau in January 2005, the three configurations tested within this study represented light rocky desertification (LRD), moderate rocky desertification (MRD) and severe rocky desertification (SRD).

The bedrock exposure represented (Figure 1) the vertical projected area of exposed rocks per unit area. The parameter values were 30%, 40% and 50%, which were equivalent to 30, 40 and 50 limestone rocks exposed on the soil surface.

Ryegrass (*Lolium perenne* L.), a common herbaceous plant in the study area, was selected for this study and the vegetation coverage reproduced within the experimental model was referred to the proportion of the vertical projection area of the vegetation to the land surface area. The average thickness of the soil layer represented the ratio of the total volume of the filled soil to the base area of the pond. The rock–soil contact area represented the surface area per unit volume of rock in contact with soil. The larger the rock-soil contact area is, the greater is the degree of rock fragmentation and the karst fissures are more fully developed.

The rainfall intensity was designed as three gradients typical of light rain intensity (30 mm h<sup>−</sup>1), moderate rain intensity (60 mm h−1) and heavy rain intensity (90 mm h−1). The simulated rainfall lasted 1h for each condition. Each rainfall event was repeated 3 times.

The underground runoff was collected every 15 min with the use of a plastic bucket. This procedure was conducted every 1 h after the underground runoff rate was stabilized. To calculate the underground runoff rate Ur <sup>=</sup> (L·h<sup>−</sup>1), the following equation was used:

$$\text{Ur} = \frac{\text{V}}{\text{AT}} \tag{2}$$

where V is the volume of underground runoff (L); A is the base area of the pond (m2); and T is the duration of underground runoff (h).


**Table 1.** Experimental configurations.

#### **3. Results**

#### *3.1. E*ff*ects of Rainfall Intensity and Rocky Desertification on the Dynamics of Underground Runo*ff

Figure 2 shows the dynamics of underground runoff for each of the three configurations tested (LRD (Figure 2a,d,g), MRD (Figure 2b,e,h) and SRD (Figure 2c,f,i)) under different rainfall intensities. It is possible to notice that despite under dissimilar rainfall intensity, underground runoff can be divided into four stages: (i) runoff rising stage, (ii) rapid declining stage, (iii) slow declining stage and (iv) recession stage.

**Figure 2.** Dynamics of underground runoff for light rocky desertification (LRD) (**a**,**d**,**g**), moderate rocky desertification (MRD) (**b**,**e**,**h**) and severe rocky desertification (SRD) (**c**,**f**,**i**) under 30 mm/h (**a**–**c**), 60 mm/h (**d**–**f**) and 90 mm/h (**g**–**i**).

Additionally, different degrees of rocky desertification have an influence on the dynamic process of the groundwater flow. Considering, for example, a moderate rainfall intensity (60 mm/h), the response time to rainfall of underground runoff can be simplified as LRD < MRD < SRD (shown in Figure 3).

**Figure 3.** Time to underground runoff under different rainfall intensities. Note: In the figure, the capital letter represented the significant difference between different rainfall intensity under a certain degree of rocky desertification (*p* < 0.05), and the lowercase letter represented the significant difference between different degrees of rocky desertification under a certain degree of rainfall intensity (*p* < 0.05).

The initial time for the generation of underground runoff of LRD was recorded to be the shortest, 51 min. For MRD and SRD, this time was 18.18% and 45.45% later than the one recorded for LRD. The time to reach the runoff "flood" peak was 105 min, 75 min and 45 min, respectively, indicating that the time to reach the flood peak gradually shortened with the increase of the severity of rocky desertification. Furthermore, as shown in Figure 4, the maximum rate of underground runoff was achieved with the highest rocky desertification. The maximum runoff rate (8.24 L/h) of SRD was 1.20 times and 1.11 times higher than the ones recorded for the LRD and MRD scenarios, respectively.

**Figure 4.** Maximun underground runoff rate under different rainfall intensities.

As initially stated in Section 1, rainfall was considered the main dynamic factor to drive the underground runoff [20], and the rainfall intensity confirmed to have an important impact on the characteristics of the underground runoff.

It can be seen from Figure 2 that with the increase in rainfall intensity, the duration of the rising stage to reach the maximum peak of flow decreased. In parallel, the underground runoff velocity reached the peak value in a shorter time, and the maximum rate of underground runoff significantly increased with the increase in rainfall intensity. Furthermore, the maximum rate of underground runoff in severe rocky desertification was the most affected by rainfall intensity (Figure 4).

For SRD, the maximum rate of underground runoff recorded under the condition of light rain increased 2.71 times under the condition of moderate rain and 24.26 times under the condition of heavy rain. The initial runoff generation duration of the underground runoff decreased significantly with the increase in rainfall intensity (Figure 3).

The underground runoff of LRD was mostly affected by rainfall intensity. Compared with the initial runoff duration under the condition of light rain, the initial runoff duration under the condition of medium rain and heavy rain was reduced by 50.00% and 81.82%, respectively.

#### *3.2. Volume and Percentage of Underground Runo*ff *Based on Di*ff*erent Degrees of Rocky Desertification*

The amount of cumulative underground runoff for different rocky desertification was SRD > MRD > LRD (Figure 5 below). Under light rain conditions, the cumulative underground runoff of SRD was 2.11% and 0.77% larger than that of LRD and MRD. Under moderate rain condition, the cumulative underground runoff of SRD was 17.56% and 9.53% larger than that of LRD and MRD. Under heavy rain condition, the cumulative underground runoff of SRD was 5.47% and 2.67% larger than that of LRD and MRD. The cumulative underground runoff significantly increased with the increase in rainfall intensity, and the cumulative underground runoff of LRD was the most affected by the rainfall intensity (Figure 5). The cumulative underground runoff of LRD increased 1.75 times under heavy rain conditions than that of under light rain conditions.

**Figure 5.** Cumulative underground runoff under different rainfall intensities.

The relationships identified between underground runoff duration and cumulative underground runoff indicated that there was a logarithmic function between the parameters considered (rainfall and rock desertification, as shown in Table 2) and the correlation coefficients were all above 0.85 (Table 2), confirming the trustworthiness of these results.


**Table 2.** Regression analysis of cumulative underground runoff and runoff duration.

*3.3. Characteristics of VWC and Runo*ff *Sources for Di*ff*erent Degrees of Karst Rocky Desertification*

The saturated VWC of each soil layer was between 0.45 and 0.52, as shown in Figure 6. Under the conditions of light, medium and heavy rain, the timings for soil on top (5 cm) to reach saturation were 20–30 min, 15–20 min and 10–15 min, respectively. Meanwhile, the timings for the bottom soil (70 cm) to reach saturation were about 100–150 min, 80–100 min and 45–55 min, respectively. Under light rain conditions, the time to generate underground runoff was 60 min after the soil reached the saturation, typical of Dunne overland flows [41,42]. However, under moderate and heavy rain conditions, the initiation of underground runoff happened before the soil reached the saturation, which is typical of Horton overland flows. The VWC began to decrease after the rainfall stopped for 5–10 min. For example, the VWC at 5 cm and 15 cm of LRD significantly decreased after reaching saturation, decreased by 0.05 in 20 min and finally continued to steadily decline. However, the VWC at 30 cm, 50 cm and 70 cm decreased slowly after saturation.

**Figure 6.** Dynamics of volumetric water content (VWC) for LRD (**a**,**d**,**g**), MRD (**b**,**e**,**h**) and SRD (**c**,**f**,**i**) under 30 mm/h (**a**–**c**), 60 mm/h (**d**–**f**), 90 mm/h (**g**–**i**).

The underground runoff (V) was mainly produced by the content of water within the soil and the infiltration rate through the fissures (F). The volume of fissure flow VF (in liters) can be obtained by applying a water balance, neglecting the effects of evaporation as follows:

$$\text{VF} = \text{P} - \Delta \text{Vsoil} - \text{V} \tag{3}$$

where P is the precipitation (mm); ΔVsoil is the change of soil water content (l); and V is the underground runoff (l).

The proportion of fissure flow increased significantly with the increase in rainfall intensity (Figure 7). Under the conditions of light and moderate rain intensities, the rate of water within the soil water was higher than the fissure flow, which only accounted for 10%–13% and 25%–35% of the total precipitation. However, under the condition of heavy rain, the proportion of fissure flow was relatively high, accounting for 47%–50% of the total precipitation. Moreover, the proportion of fissure flow increased with the increase of rocky desertification degree. The proportion of fissure flow of SRD was 3%–10% larger than that of LRD.

**Figure 7.** Fissure flow ratio under different rainfall intensities.

From the start of the underground runoff to its peak, the water balance analysis was continuously performed to obtain the ratio of soil water and fissure flow vs. the underground runoff. Under light rain conditions, underground runoff was mainly generated by soil water, (85%–95%), while fissure flow only accounted for 5%–15%. Under moderate rain conditions, runoff was mainly produced by fissure flow (72%–87%), while soil water accounted only for 13%–27%. Under heavy rain conditions, the soil was not saturated due to the short time (20 min) in which the underground runoff was formed, hence fissure flow was the main source of the underground runoff.

#### **4. Discussion**

Experimental studies conducted confirmed that underground fissures, which are unique to the karst areas around the world, create migration channels for groundwater and soil losses [32]. Results have established that soil and water losses were caused by factors such as the exposed rate of the bedrock, the degree of karst fissures and the rainfall intensity [43,44]. However, in order of impact, rainfall is the most important driving factor, followed by the degree of karst fissures [45,46].

There were significant differences in underground runoff characteristics between different degrees of rocky desertification; results confirmed that the response time of the underground runoff to the rainfall increased with the increase in degree of rocky desertification.

As shown in Table 1, with the increase in the degree of rocky desertification, the rock–soil contact area and the rock–soil ratio gradually increased, which indicated that the main channel of runoff migration, the soft–hard interface between rock and soil [47], gradually increased; meanwhile, the width of the cracks between the rocks gradually decreased. The size of the cracks between the rocks has demonstrated to have significant influence on the velocity of water moving through the rocks [48], and this velocity in the fractures was much slower than the pipelines [49].

As shown in Figure 6, the water infiltration rate of LRD was faster than the one measured for MRD and SRD. For example, under light rain conditions, at 50 cm, soil VWC for LRD, MRD and SRD started to increase after 65 min, 70 min and 85 min, respectively, while at 70 cm, soil VWC for LRD, MRD and SRD started to increase after 85 min, 90 min and 115 min, respectively. Furthermore, the rock–soil contact area of SRD was 45.69% higher than LRD, but the time needed for the water to move from soil at 50 cm to soil at 70 cm was 10 min longer. This behavior also indicated that although the rock–soil contact area and migration channels for groundwater were the largest, these limited the velocity of the water infiltrating. Therefore, the authors believe that this was a critical point for quantifying the effect of rock–soil contact area on the water infiltration rate. When the rock–soil contact area was greater than the critical value, the water infiltration rate decreased with the increase in water transport channels and the decrease in crack width, which consequently affected the initiation of underground runoff. However, once the underground runoff was produced, a preferred flow path had been established connecting the soil and rock interface. As the results showed, SRD firstly reached the maximum underground runoff rate due to the largest rock–soil contact area. The maximum underground runoff rate and the total amount of underground runoff were also higher than those recorded for LRD and MRD, confirming that the rock–soil contact area had an important influence on the characteristics of underground runoff. The increase in rock–soil contact area will worsen the soil and water losses in the karst area, reducing the time to generate underground runoff and increasing the water leakages.

#### **5. Conclusions**

Due to heavy rainfall events becoming more frequent due to the climate change across the world, it is necessary to carefully consider strategies to guarantee the safety of water resources and ecological systems, especially in karst landform areas, typical in Southwest China, which have not been addressed completely to date. This study included the simulation of a unique "dual structure" karst area focusing on the soil–rock integration, and the experimental representation of underground runoff with the purposes of (1) providing a better understanding of underground runoff under different rainfall intensities and (2) discussing the influence of rocky desertification on underground runoff. Results can be summarized as follows:


The severe loss of groundwater in Southwest China has severely restricted local ecological restoration and economic development, but which can, however, be controlled through sustainable engineering measures. Resources should be directed to reduce the amount of runoff leakage and preserving underground runoff, facilitating water storage. The authors suggest that techniques such as filling cracks with gravels [50] or filling cracks with mud [34] should be implemented to expedite this task. Although engineering measures can be effective, costs are high, and operations are difficult. Planting vegetation could definitely be another solution to improve soil characteristics and enable

rainfall interception. At the same time, the distribution of plant roots in underground fissures could reduce the channels for the underground runoff to a certain extent [51].

Due to the nature of the stones used for this study, which are not irregular, the authors suggest that this limitation can be solved by considering this scenario within a new bespoke experimental setup, incorporating also larger field areas. Furthermore, future research should focus on (i) the effects and mechanisms of plant measures on groundwater loss control; (ii) the specific impact of the rock–soil interface; (iii) the definition of fissure characteristics, such as the density values or widths; and (iv) further strengthening the model research of the karst dual hydrological unit to replicate field observations and experimental simulations, crucial to provide a more comprehensive understanding of the complex underground runoff processes in the karst areas.

**Author Contributions:** Conceptualization, G.L., L.W., J.F. and J.Z.; data curation, G.L., L.W., J.Z. and M.R.; formal analysis, G.L., L.W., J.L., B.W., J.F. and J.Z.; funding acquisition, L.W. and J.Z.; investigation, G.L., L.W., J.L., J.F. and M.R.; supervision, G.L., L.W., J.L., B.W. and J.Z.; validation, G.L., L.W., J.L., B.W., M.R. and J.Z.; writing—original draft, G.L. and L.W.; writing—review and editing, G.L., L.W and M.R. All authors have read and agreed to the published version of the manuscript.

**Funding:** The research was supported by the Fundamental Research Funds for the Central Universities of China (BLX201709), the National Natural Science Foundation of China (31870707; 31700640), and the National Basic Research Program of China (2016YFC0502502; 2016YFC0502504).

**Acknowledgments:** We would like to thank Lechuan Qiu, Jianghua Liao, Jiatong Zhang and Weixin Zhang for their invaluable help with measurements taken in the field and in the laboratory.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 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 (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Analysis of Long-Term Water Level Variations in Qinghai Lake in China**

#### **Jianmei Fang 1, Guijing Li 1, Matteo Rubinato 2, Guoqing Ma 3, Jinxing Zhou 1, Guodong Jia 1, Xinxiao Yu 1,\* and Henian Wang <sup>4</sup>**


Received: 17 September 2019; Accepted: 9 October 2019; Published: 14 October 2019

**Abstract:** Qinghai Lake is the largest inland saline lake on the Tibetan Plateau. Climate change and catchment modifications induced by human activities are the main drivers playing a significant role in the dramatic variation of water levels in the lake (Δ*h*); hence, it is crucial to provide a better understanding of the impacts caused by these phenomena. However, their respective contribution to and influence on water level variations in Qinghai Lake are still unclear and without characterizing them, targeted measures for a more efficient conservation and management of the lake cannot be implemented. In this paper, data monitored during the period 1960–2016 (e.g., meteorological and land use data) have been analyzed by applying multiple techniques to fill this gap and estimate the contribution of each parameter recorded to water level variations (Δ*h*). Results obtained have demonstrated that the water level of Qinghai Lake declined between 1960 and 2004, and since then has risen continuously and gradually, due to the changes in evaporation rates, precipitation and consequently surface runoff associated with climate change effects and catchment modifications. The authors have also pinpointed that climate change is the main leading cause impacting the water level in Qinghai Lake because results demonstrated that 93.13% of water level variations can be attributable to it, while the catchment modifications are responsible for 6.87%. This is a very important outcome in the view of the fact that global warming clearly had a profound impact in this sensitive and responsive region, affecting hydrological processes in the largest inland lake of the Tibetan Plateau.

**Keywords:** climate change; water levels; causes and implications; Qinghai Lake, Tibetan Plateau

#### **1. Introduction**

The observed climate changes [1] had a significant impact on physical and natural processes on Earth during the past decade. The IPCC's (Intergovernmental Panel on Climate Change) latest report pointed out that if global warming will continue at its current rate, it could reach an increase in temperature up to 1.5 ◦C between 2030 and 2052 [2], causing the rising of sea levels as well as warming of water surfaces in oceans and lakes. Furthermore, human activities also had a strong impact on hydrological processes considering increased water consumption and situations of water shortage recorded around the world during the last decade. Hence, recent research focused on the impact of climate change and human activities on hydrological processes because this topic was identified by many researchers as a priority [3,4] when planning for the future and making new developments

more sustainable. Lake ecosystems usually provide indicators (i.e., water temperature, water levels, dissolved organic carbon (DOC)) of climate change, either directly or indirectly [5]. Recently, many studies have also been completed to investigate what has caused water level variations of lakes, such as the ones conducted on (i) the North American Great Lakes [6–8]; (ii) Lake Chad [9]; (iii) the Salton Sea [10] in the United States; (iv) Lake Lisan, Dead Sea rift [11]; and (v) Poyang Lake in China [12–14]. Summarizing the results achieved to date, hydrological conditions of each lake are affected by the lake's location, upstream boundaries, geographical climate and specific human activities undertaken on it such as residential developments, industry and irrigation tasks; hence, it is necessary to figure out the key factors that affect water levels to develop methods and procedures that can regulate those that alter or alleviate hydrological extreme processes in lakes.

The Tibetan Plateau (TP), known as "the roof of the world", "the third pole" and "the water tower of Asia", is the largest plateau in China and the highest plateau in the world [5,15], and it is considered the perfect location to identify the effects of global climate change [16–18]. Qinghai Lake is the largest inland saline lake on TP, and it has attracted extensive attention due to its special geographical location and its wide area characterized by fragile ecosystems. Over the last years, researchers have reported a drastic change in water levels in Qinghai Lake, indicating that the ecological environment around it is undergoing a rapid evolution [19–22]. Typically, inland closed lakes with no outlet streams are ideal to distinguish hydrologic processes and phenomena affecting the water balance because changes in water levels can result from limited factors such as precipitation, evaporation, groundwater infiltration [23–25] and the presence of specific vegetation [19].

The most effective way to estimate water levels in lakes is by applying the water-balance equation model, where the gain or loss of water directly reflects the changes in water levels [26,27]. Multiple inland lakes, due to the lack of funding and therefore resilient and accurate equipment to monitor basic data, are considered not appropriate to identify and quantify objectively the factors affecting the water balance. Nevertheless, Qinghai Lake, being a closed one with inlet river streams and without outlet river streams, is an ideal place to study, especially having the availability of meteorological and hydrological monitored data. To date, the study conducted by Li et al. [28] calculated the main water balance estimation of Qinghai Lake, while Cui et al. [20] preliminarily analyzed the climatic factors that affect the water level variations of Qinghai Lake. Despite this, there is still a need for new studies to fully distinguish and assess the relative contribution of anthropogenic activities and climate variability to water level variations in Qinghai Lake and their impacts on the corresponding water balance.

This paper present the analysis of water level variations recorded in Qinghai Lake during the last 57 years, examining the evolution and interpreting the impacts of driving factors to better understand the hydrological process of this inland lake basin in the northeast of TP, enhancing the present understanding of climatic variations on surface changes to provide a reference for local and regional water management.

The paper is organized as follows: Section 2 presents the description of the study area, introducing the data monitored and describing the statistical analysis used. Section 3 includes the estimation of long-term variations in Qinghai Lake as well as the impacts of climate change and catchment modifications for inflow runoff to Qinghai Lake. Section 4 provides a discussion of the results obtained, and Section 5 produces a brief summary and concluding remarks of the whole study.

#### **2. Materials and Methods**

#### *2.1. Study Area and Data Availability*

Qinghai Lake basin (97◦50 ~101◦20 E, 36◦15 ~38◦20 N) is located in the northeast of TP, covering an area of 29,664 km2. The average annual air temperature ranges from <sup>−</sup>0.8 ◦C to 1.1 ◦C, and the average annual precipitation ranges from 327 to 423 mm. However, the annual precipitation is unevenly distributed, decreasing from the east and south to the west and north. The total surface runoff of local main rivers including Buha River, Shaliu River, Haergai River, Heima River and Daotang River

accounts for 83% of the total surface runoff into Qinghai Lake [28,29], with the first two (Buha & Shaliu) constituting the 64% of the surface runoff for the entire basin [30]. The main vegetation types in the basin are alpine meadows and alpine grasslands.

Qinghai Lake is the largest inland saline lake in this basin with an area of 4400 km<sup>2</sup> (in 2016), and it is located at an altitude of 3193 m above sea level. As previously mentioned in Section 1, this lake is a closed one with no surface water outflow. It is about 106 km in length from east to west, and 63 km in width from north to south, and 360 km in circumference [30]. The average annual air temperature above the lake is about 1.2 ◦C, and the average annual precipitation near its proximity is about 357 mm.

The datasets available used for this study can be summarized as follows:

	- i. Daily meteorological data of 14 national meteorological stations from 1960 to 2016, obtained by the China Meteorological Information Center.
	- ii. Monthly meteorological data from 1960 to 2010 at three meteorological stations, obtained by Qinghai Meteorological Bureau in China.
	- iii. Daily precipitation data of Buha River rain station from 1962 to 2016 obtained by ICQHB.
	- iv. Daily evaporation data from 1984 to 2016 at Xiashe station obtained from ICQHB.
	- v. Yearly evaporation data from 1960 to 1988 obtained from the literature [29]. (Figure 1 and Table 1).

**Table 1.** Detailed information of the meteorological stations in and around Qinghai Lake basin used to collect the datasets previously described.


**Figure 1.** The location of Qinghai Lake basin in China (top right corner), physical characteristics of Qinghai Lake in Qinghai Lake basin (centre) and the spatial distribution of the hydrological stations (red triangles) and meteorological stations (blue dots).

#### *2.2. Governing Equations*

#### 2.2.1. Lake Water Balance Model

As Qinghai Lake is a closed-catchment with no surface water outflow, the annual hydrological water balance equation can be expressed as follows:

$$
\Delta \mathbf{h} = P\_I - E\_I + R\_{l\mathbf{s}} + R\_{l\mathbf{g}} \pm \varepsilon \tag{1}
$$

where Δ*h* is the yearly water level variation, (mm); *Pl* is the yearly precipitation on the lake surface, (mm); *El* is the yearly evaporation from the lake surface, (mm); *Rls* is the yearly surface runoff into the lake, (mm); *Rlg* is yearly underground runoff on the lake bottom, (mm); ε is the error, (mm). For this watershed scenario, the surface runoff is almost equivalent to the river runoff and the slope surface runoff can be considered negligible. Δ*h* can be quantified as well as follows:

$$
\Delta h = h\_i - h\_{i-1} \tag{2}
$$

where *hi* and *hi*−<sup>1</sup> are the lake level at *i* year and at *i* − 1 year.

The yearly average water level of the lake was obtained from the daily water level data at Xiasha station. *Pl* was calculated by applying the Thiessen Polygon method focusing in the area between Buha River station, Gangcha station, Haiyan station and Gonghe station, which are the nearest stations to the lake. *El* was obtained from evaporation pan (type of E20) data at Xiashe station [29]. The yearly total surface runoff (*Rls*, mm) into the lake was obtained from the surface runoff (*Qls*, m3) of Buha River and Shaliu River by using the proportion amplification method [31].

In this paper, Δ*h* is subject to a combination of climate and human activities effects (such as farmland reclamation, grazing, afforestation which indirectly influence runoff and catchment characteristics). *Pl* and *El* represents the climate variability, *Rls* and *Rlg* are the results of the combination of climate and catchment modifications. Hence, to correctly quantify the contribution rate of climate and catchment characteristics (human induced) to Δ*h*, it is necessary to calculate accurately the contribution rate

of climate and catchment characteristics (human induced) to *Rls* and *Rlg*. Therefore, the calculation method can be applied as follows:

$$\left(\boldsymbol{R}(\Delta\boldsymbol{h})\_{\boldsymbol{\varepsilon}} = \boldsymbol{R}(\Delta\boldsymbol{h}\_{\text{Pl}})\_{\boldsymbol{\varepsilon}} + \boldsymbol{R}(\Delta\boldsymbol{h}\_{\text{El}})\_{\boldsymbol{\varepsilon}} + \boldsymbol{R}(\Delta\boldsymbol{h}\_{\text{Rls}})\_{\boldsymbol{\varepsilon}} + \boldsymbol{R}(\Delta\boldsymbol{h}\_{\text{Rl}\oplus\boldsymbol{\varepsilon}})\_{\boldsymbol{\varepsilon}}\right) \tag{3}$$

$$R(\Delta h)\_{\mathfrak{u}} = R(\Delta h\_{\mathbb{R}ls})\_{\mathfrak{u}} + R(\Delta h\_{\mathbb{R} \lg \pm \mathfrak{e}})\_{\mathfrak{u}} \tag{4}$$

where *R*(Δ*h*)*c*, *R*(Δ*h*)*<sup>u</sup>* respectively represent the contribution rate of climate and catchment modifications to Δ*h*. *R*(Δ*hPl*)*c*, *R*(Δ*hEl*)*<sup>c</sup>* respectively represent the climate change contribution of *Pl* and *El* to Δ*h*. *R*(Δ*hRls*)*c*, *R*(Δ*hRls*)*<sup>u</sup>* respectively represent the *Rls* contribution rate of climate change and catchment change to Δ*h*. *R*(Δ*hR*lg±ε) *c* , *R*(Δ*hR*lg±ε) *<sup>u</sup>* respectively represent the *Rlg* ± ε contribution rate of climate change and catchment change to Δ*h*.

#### 2.2.2. Land Use Dynamic Index

Land use change can reflect the effect as well as the intensity of human activities. The Land Use Dynamic Index was proposed by Chen et al. [32] and was adopted in this study to describe the change of land use types in the research area for a certain period (1980–2015). The calculation method was completed as follows:

$$LC = \frac{\sum\_{i=1}^{n} \Delta l I\_{\text{lin}}}{2 \sum\_{i=1}^{n} l I\_{i0}} \times \frac{1}{T} \times 100\% \tag{5}$$

where *LC* is the Land Use Dynamic Index in a certain period of time in the research area (%), Δ*Ui*in refers to the area of type *i* land use converted into the non- *i* type land use within a certain period of time in the research area (km2), *Ui*<sup>0</sup> is the area of type *i* land use at the beginning of the study period (km2), *T* is the research period (years).

#### 2.2.3. Statistical Analysis

#### (1) The Non-parametric Mann-Kendall Test

The non-parametric Mann-Kendall test [33,34] (M-K test) and the cumulative anomaly method were used to detect any point of abrupt changes in the variables considered. The M-K test has been widely applied to identify the point at which the hydrological processes change significantly due to the climate [35,36]. The details about this statistical method can be obtained in the relevant literature [37].

First, the partial M-K test statistics are calculated as:

$$S\_k = \sum\_{i=1}^k \sum\_{j=1}^{i-1} a\_{ij} \quad (k = 2, 3, 4, \dots, n) \tag{6}$$

$$\alpha\_{i\bar{j}} = \begin{cases} 1 & \mathbf{x}\_{\bar{i}} > \mathbf{x}\_{\bar{j}} \\ 0 & \mathbf{x}\_{\bar{i}} \le \mathbf{x}\_{\bar{j}} \end{cases} \quad 1 \le j \le i \tag{7}$$

Statistical variable *UF* is adopted and defined as:

$$
\Delta IF = \begin{array}{c c} \frac{\text{S} - E(\text{S}\_k)}{\sqrt{Var(\text{S}\_k)}} & (k = 1, 2, 3, \dots, n \text{ }) \end{array} \tag{8}
$$

$$E(S\_k) = \frac{k(k-1)}{4} \tag{9}$$

$$Var(S\_k) = \frac{k(k-1)(2k+5)}{72} \tag{10}$$

Proceed to Equation (11) putting the data sequence *x* in reverse order:

$$\begin{cases} \text{LIB}\_k = -\text{LIF}\_{k'}\\ k' = n + 1 - k \end{cases} \quad (k = 1, 2, 3, \dots, n \ )\tag{11}$$

In the M-K curve, if the value of the intersection of the curve forward (*UF*) or the curve backward (*UB*) is greater than 0, this suggests that the record sequence shows an upward trend; less than 0 suggests a downward trend. When the record exceeds the critical line (Given the significance level α = 0.05, the critical lines *U*0.05 = ±1.96), this suggests that an increase or decrease in the trend may be significant. The range of exceeding the critical line is the time zone in which the abrupt change occurs. If there is an intersection between the curves of *UF* and *UB* in the range of the critical lines, the time of the intersection is the time of the abrupt change started [38].

#### (2) The Cumulative Anomaly Method

The cumulative anomaly method is widely used to indicate the runoff [39], precipitation and other factors that deviate from the normal situations, focusing on the difference between a certain value and the average value of a series [38].

#### (3) The Principal Component Regression Analysis

The principal component regression (PCR) analysis [40] is a combination of principal component analysis and regression analysis. Typically, this method considers regressing the outcome on a set of covariates based on a standard linear regression model, using PCA (principal component analysis) for estimating the unknown regression coefficients. Generally, only a subset of all the principal components for regression is used; hence, PCR tends to act as a regularized procedure.

#### (4) The Grey Relational Analysis

The grey relational analysis [41] is adopted in this study to solve uncertain problems such as limited data and incomplete information by calculating the grey correlation degree γ*i*, quantifying the correlation degree among the influential factors of underground runoff.

#### (5) The Least Square Method

The least square method [42] is applied to procure unknown data and minimize the sum of squared errors between the obtained data and the actual data. The least square method can also be used for curve fitting.

#### (6) The Partial Least Squares Regression Method

The partial least squares regression method [43] is a combination of multiple linear regression analysis, canonical correlation analysis and principal component analysis, reflecting the influence of the sample population on the predicted values and fully considering the influence of the comprehensive effect between individual factors on the predicted ones.

#### 2.2.4. Sensitivity Analysis Based on the Budyko Framework

Climate change and human activities are the most important drivers to determine the river hydrological process of the catchment [44]. In this study, the sensitivity coefficient method [45] based on Budyko Theory [46] was used to quantitatively separate the impacts of climate change and human activities on the variations of streamflow into Qinghai Lake. The theoretical equation of Budyko curve [47] can be applied as follows:

$$\frac{ET}{P} = 1 + \frac{ET\_0}{P} - \left[1 + \left(\frac{ET\_0}{P}\right)^{\omega}\right]^{1/\omega} \tag{12}$$

where *ET* is the evapotranspiration of the upper catchment area, (mm); *P* is the precipitation of the catchment area, (mm); *ET0* is the potential evapotranspiration of the catchment area, (mm); the empirical parameter ω represent catchment characteristics, such as human activities, land use, vegetation, topography, and properties of soil [48,49], [ω ∈ (1, ∞)].

The change of surface runoff in a given basin can be characterized by climate and human activities changes as follows:

$$
\Delta Q = \Delta Qc + \Delta Qu \tag{13}
$$

where Δ*Qc* and Δ*Qu* represent the surface runoff variation caused by climate change and human activities changes, respectively. The surface runoff variation caused by climate change can be expressed by the following formula [45]:

$$
\Delta \mathbf{Q} \mathbf{c} = \frac{\partial \mathbf{Q}}{\partial P} \times \Delta P + \frac{\partial \mathbf{Q}}{\partial ET\_0} \times \Delta ET\_0 \tag{14}
$$

The surface runoff variation caused by human activities can be expressed by the following formula:

$$
\Delta Q \mu = \frac{\partial Q}{\partial \omega} \times \Delta \omega \tag{15}
$$

where Δ*P* is the variation of precipitation, Δω is the variation of the empirical parameter ω of a given catchment; <sup>Δ</sup>*ET*<sup>0</sup> is the potential evapotranspiration variation; <sup>∂</sup>*<sup>Q</sup>* <sup>∂</sup>*<sup>P</sup>* , <sup>∂</sup>*<sup>Q</sup>* ∂*ET*<sup>0</sup> , <sup>∂</sup>*<sup>Q</sup>* ∂ω respectively represent the sensitivity coefficient of runoff to precipitation, runoff to potential evapotranspiration, runoff to precipitation, runoff to the empirical parameter represent catchment characteristics. All of the sensitivity coefficients can be calculated as follows:

$$\frac{\partial Q}{\partial P} = \left[1 + \left(\frac{ET\_0}{P}\right)^{\omega}\right]^{(1/\omega - 1)}\tag{16}$$

$$\frac{\partial \mathcal{Q}}{\partial ET\_0} = \left[1 + \left(\frac{P}{ET\_0}\right)^{\omega}\right]^{(1/\omega - 1)} - 1\tag{17}$$

$$\frac{\partial Q}{\partial \omega} = \left[P^{\omega} + ET\_0^{\omega}\right]^{1/\omega} \cdot \left[\left(-\frac{1}{\alpha^2}\right) \cdot \ln\left(P^{\omega} + ET\_0^{\omega}\right) + \frac{1}{\alpha} \cdot \frac{1}{P^{\omega} + ET\_0^{\omega}} \cdot \left(\ln P \cdot P^{\omega} + \ln ET\_0 \cdot ET\_0^{\omega}\right)\right] \tag{18}$$

In this paper, the potential evapotranspiration at the meteorological stations was calculated by applying the FAO56 method, Penman-Monteith model [50], because previous literature [51,52] has demonstrated how these methods are reliable to estimate potential effects of climate change on the calculation of the evaporation as well as the influence of climate change on water cycles. *P* and *ET0* of the entire basin were obtained by applying the area-weight method of Tyson Polygon.

#### **3. Results and Analysis**

#### *3.1. Long-Term Variations in Water Levels and the Hydro-Climatic Factors*

#### 3.1.1. Long-Term Variations in Water Levels

Figure 2 shows the annual water levels of Qinghai Lake recorded during the period 1960–2016. Overall, comparing datasets within other locations in the semi-arid areas of Western China, the water level varied significantly (≈3.5 m) and it is possible to notice a clear inflection point recorded in 2004. Herein, the analysis of the graph was divided into two periods to simplify the procedure: period I was selected to be between 1960 and 2004, while period II was selected to be between 2005 and 2016. The annual water level of the lake declined at the rate of 7.84 cm/year (*P* < 0.001), with a total decrease of 3.46 m in period I, while the annual water level of the lake has risen at the rate of 13.80 cm/year (*P* < 0.001), with a total increase of 1.49 m in period II.

The variation of the water levels of the lake (Δ*h*) reflected the acquisition and loss of water volume over the years due to multiple factors, and the trend is shown in Figure 3. According to the results obtained, Δ*h* tended to increase during period I (*R*<sup>2</sup> = 0.0105, *P* < 0.01) as well as during period II (*R*<sup>2</sup> = 0.0291, *P* < 0.01). The increasing rate of the water level in period II was notably faster than the one in period I, and the water level of the lake in 1960 could be reached again by 2030 if the present increasing rate continues constantly.

**Figure 2.** The annual water levels of Qinghai Lake recorded during the period 1960–2016.

**Figure 3.** The trend of water level variation (Δ*h*) during the period of 1960–2016.

#### 3.1.2. Analysis of Hydro-Climatic Factors Influencing Water Levels

Δ*h* was dependent on the lake hydro-climatic conditions *Pl*, *El*, *Rls* and *Rlg* ± ε. All together, these variables affected the rise or fall of the water level in the lake, and the relationship between them and the corresponding variations in the water levels of the lake are shown in Figure 4 and Table 2.

**Figure 4.** Variations of the hydro-climatic factors of Qinghai Lake during the period of 1960–2016. *Pl* is the precipitation over the lake (mm), *El* is the evaporation from the lake surface (mm), *Rls* is the surface runoff into the lake (mm), *Rlg* ± ε is the underground runoff and study error (mm), Δ*h* is lake level variation (mm).



\* The bracketed values refer to the percentage of total input or output represented by average yearly volumetric flux (mm) changes at different periods. Sign + represents water in, while sign - represents water out.

During the period under investigation, surface runoff into the lake (*Rls*) was mainly due to the Buha River and Shaliu River, which contributed about 48.70% of the total water in the lake. The underground runoff on the lake bed and the associated error (*Rlg* ± ε) accounted for 5.56% of the total water intake of the lake, and this value fluctuated significantly. The precipitation on the lake surface (*Pl*) contributed about 45.74% of the total water in the lake while the evaporation from the lake surface (*El*) contributed about 100% of the total water removed from the lake. As possible to notice from Figure 4, *Pl*, and *Rls* had similar trends, while this could not be confirmed for *El* and *Rlg*. The peak of Δ*h* often corresponded to the peak of *Pl*, *Rls*.

The annual average values for each hydro-climatic variable *Pl*, *Rls*, *Rlg* ± ε, and *El* were 381.59 (mm), 406.32 (mm), 46.42 (mm), and 869.77 (mm), respectively. In period I, approximately 45.67% of the total water input into the lake came via *Pl*, with 45.18% of water input coming from *Rls*, and a small

fraction of water input was due to *Rlg* ± ε. The whole outflow was estimated to be associated with *El*. In period II, approximately 43.38% and 56.62% of the total water input was associated with *Pl* and *Rls*, respectively, while *El* contributed 93.5% to the outflow with a small fraction of water escaping the lake attributed to *Rlg* ± ε (6.5%). This indicated that the water balance of Qinghai Lake was mainly determined by *Pl*, *Rls* and *El*. Therefore, the authors can confirm that *Rlg* ± ε*,* always being a small percentage, accounted for a small proportion of the water balance of Qinghai Lake.

#### *3.2. Causes of Changes in Water Levels of the Lake*

#### 3.2.1. Impact of Climate Change on Water Levels

Figure 4 shows the relationship between *Pl*, *El* and Δ*h*. The correlation coefficient between *Pl* and Δ*h* calculated was 0.356 (*P* < 0.01). The fluctuation range for *Pl* was estimated between 277.2 and 546.4 (mm), where the mean value was obtained equal to 381.6 (mm), showing an upward trend at the rate of 1.4347 mm/year (*R*<sup>2</sup> = 0.164, *P* < 0.01). The correlation coefficient between *El* and Δ*h* was −0.705 (*P* < 0.01) and the fluctuation range for *El* was between 702.6 and 1070.5 (mm), where the mean value was 869.8 (mm), showing a downward trend at the rate of 2.2823 mm/year (*R*<sup>2</sup> = 0.290, *P* < 0.01). Changes in behavior for *Pl* were consistent with the fluctuations of Δ*h*, and the peaks of Δ*h* were noticed to be delayed by 1 year when comparing them with those associated with *Pl*. However, *El* was generally contrary to the fluctuations of Δ*h*.

Based on the results achieved, *El* and *Pl* were identified as the main important climate factors affecting the water level changes, and the correlation relationships between *El* and Δ*h* estimated were of higher quality than those obtained between *Pl* and Δ*h*.

It was clear that at peaks in the rising periods for Δ*h* corresponded to higher values of *Pl* but lower values of *El* (e.g., 1968, 1989, and 2012). Furthermore, decreasing peaks of Δ*h* corresponded to lower values of *Pl* but higher values of *El* (e.g., 1980 and 1991). Precipitation rates were quantified to affect both the runoff of the inflow rivers and underground runoff acting on the water level changes. Finally, evaporation was selected as the only factor of climate influencing water "exiting" the lake, playing a significant role in the fluctuation of the water level.

3.2.2. Impact of Human Activities on Catchment Modifications and Consequently on Water Levels

Figure 5 illustrates 7 different years spanning 1980 and 2015, to highlight the land use changes from 1980 to 2015 obtained by the superposition function fitted within ArcGIS10.2. Despite being present and noticeable, land use changes observed were not particularly significant as possible to notice in Figure 5 (Tables 3 and 4).


**Table 3.** The land use dynamic attitude (LC) from 1980 to 2015.



**Figure 5.** The changes in land use in Qinghai Lake basin from 1980 to 2015.

As the overall topography of Qinghai Lake basin is relatively gentle, the production and confluence of the runoff typically have longer durations, highly dependent on the catchment. The empirical parameter representing land surface characteristics of the basin (ω) was calculated by applying the least square method [45] according to Equation (6) and results are displayed in Figure 6. The correlation coefficient between ω and Δ*h* was −0.262 (*P* < 0.05). Figure 6 and Table 4 confirmed as previously stated that there were little changes in the land use, and the main reasons causing ω changes were not associated with the land change use, but were probably due to changes in local vegetation and soil conditions.

**Figure 6.** Variation the empirical parameter representing catchment characteristics (ω) and Δ*h* in Qinghai Lake basin during the period of 1960–2016.

By summarizing all the contributions, results have also confirmed that affecting factors Δ*h*, *Pl* and *El* were entirely attributable to climate change, while *Rls* and *Rlg* ± ε were the results of the joint action between climate and catchment modifications. Therefore, it became necessary to focus on *Rls* and *Rlg* ± ε and to quantify their impact due to climate change and catchment change.

#### 3.2.3. Impact of Climate and Catchment Modifications on the Surface Runoff (*Rls*)

Surface runoff is a crucial variable to consider when completing any lake water balance [20,53] and in this study it accounted for 45.18~56.62% of the total lake inflow during the study period, which demonstrates how this parameter was a key factor affecting the water level variations in the lake. The correlation coefficient between *Rls*, and Δ*h* was calculated to be 0.590 (*P* < 0.01).

Surface runoff was generated from the surrounding catchment area of about 25,000 km<sup>2</sup> (obtained by subtracting the lake area from the total basin area). Surface runoff (mm) was obtained by dividing the annual total runoff of the basin (m3) by the annual catchment area. The catchment area of annual maximum, annual minimum and mean annual values was 25,439.71 km2 (in 2004), 25,136.70 km2 (in 1960), and 25,308.56 km2, respectively. In this paper, the Mann-Kendall test method (M-K test) and the cumulative anomaly method were used to identify remarkable changes in the variables' behaviors (i.e., surface runoff, precipitation and potential evapotranspiration). As shown in Figure 7a, the precipitation in the basin is always rising (*UF* > 0), and the precipitation trend shows a significant change from 2005 (*UF* > 1.96). The intersection of *UF* and *UB* curves indicates an abrupt change point in 2003. Furthermore, when focusing on the evapotranspiration trends in Figure 7b, and intersection point was noticed between the two curves in 1998. The M-K test failed to identify any abrupt change point in the trend of the surface runoff recorded, while the cumulative anomaly analysis method correctly estimated it as observed in 2004 (Figure 7c). The results indicated that surface runoff, precipitation and lake

water levels were closely related. Additionally, the year corresponding to the anomalies noted in the variation of the lake water levels and the surface runoff slightly lagged behind the year corresponding to the abrupt change point related to the precipitation. Based on the cumulative anomaly analysis method, the runoff series were divided into two periods like the variation of water levels: period I (1960–2004) and period II (2005–2016), which enabled the authors to calculate the basin characteristic parameters and sensitivity coefficients for period I and period II that are presented in Table 5.

**Figure 7.** *Cont.*

**Figure 7.** Mann-Kendall test of annual precipitation from 1960 to 2016 (**a**), Mann-Kendall test of annual potential evapotranspiration (**b**), and Cumulative anomaly of the annual surface runoff (**c**) in Qinghai Lake basin.

**Table 5.** Basin characteristic parameters and sensitivity coefficients in each study period.


According to the results, the total surface runoff variation was measured as +33.9 mm. Contribution rates of climate change and catchment modifications to this variation were calculated using Equations (14) and (15) and the results show that the main cause of runoff change from 1960 to 2016 was climate change (producing an increased surface runoff by 25.54 mm for a corresponding contribution rate of 80.19%). On the other hand, the effect caused by catchment modifications was not to be considered a negligible factor, considering that it generated an increase in surface runoff of 6.31 mm and hence its contribution rate was 19.81%. Error was estimated to be 2 mm, corresponding to 5.91% of the total variation (Table 6).

**Table 6.** The results of attribution analysis of runoff change.


3.2.4. Impact of Climate and Catchment Modifications on the Underground Runoff (*Rlg* ± ε)

*Rlg* ± ε was included in the water balance equation model but due to limitations in data availability, it was more challenging to accurately estimate the effects on its variations due to climate and catchment

modifications. However, by applying the equations described in Section 2, the correlation between *Rlg* ± ε and Δ*h* was calculated to be 0.143.

Underground runoff could be affected directly and indirectly by multiple factors; therefore, the authors believed it was not appropriate to complete a comprehensive and systematic analysis of its dynamic changes by using a single factor analysis. Hence, the principal component regression method [40] was used. The factors *xi* identified to influence the underground runoff (*y*) can be summarized as follows:


The correlation analysis is shown in Table 7; the multi-collinearity among the influential factors, and the grey relational degree analysis (Table 8) show that each factor had a closed relationship with *y*, and the grey relational degree of each factor *xi* is displayed as γ*i*.


**Table 7.** Correlation coefficient matrix of influential factors.

**Table 8.** Grey relational degree matrix of influential factors on subsurface runoff.


It can be seen from Table 9 that the cumulative contribution rates of the first, second and third principal components (*F1*, *F2*, *F3*) were more than 80%, indicating that they basically contained all the information of the original impact factors. Subsequently, the three principal components were used to analyse *y*. According to Table 10, the linear equations were obtained as follows:

$$F\_1 = 0.3589\mathbf{x}\_1 - 0.5745\mathbf{x}\_2 + 0.3797\mathbf{x}\_3 + 0.2533\mathbf{x}\_4 + 0.5738\mathbf{x}\_5 + 0.0598\mathbf{x}\_6\tag{19}$$

$$F\_2 = 0.4464\mathbf{x}\_1 + 0.5176\mathbf{x}\_2 + 0.1273\mathbf{x}\_3 - 0.0852\mathbf{x}\_4 + 0.119\mathbf{x}\_5 + 0.7037\mathbf{x}\_6\tag{20}$$

$$F\_3 = 0.4448\mathbf{x}\_1 - 0.0788\mathbf{x}\_2 - 0.2944\mathbf{x}\_3 + 0.6988\mathbf{x}\_4 - 0.4700\mathbf{x}\_5 - 0.0069\mathbf{x}\_6\tag{21}$$

Principal component *F*<sup>1</sup> could be almost interpreted as precipitation on the lake surface area (*x*2) and precipitation across the entire basin area (*x*5), principal component *F*<sup>2</sup> as the empirical parameter representing catchment characteristics (*x*6), and principal component *F*<sup>3</sup> as lake surface evaporation (*x*4). The *catchment factor* represents the catchment change; the precipitation and lake evaporation factor represents the climate change. Taking *F*<sup>1</sup> (precipitation factor), *F*<sup>2</sup> (catchment factor), and *F*<sup>3</sup> (evaporation factor) as independent variables and *y* as a dependent variable, multiple linear regression was carried out, and the partial least squares [43] was adopted to obtain the following equation:

$$y = 384.0088 - 0.3556F\_1 - 0.3061F\_2 + 1.9858F\_3 \tag{22}$$

<sup>\*</sup> *P* < 0.05, \*\* *P* < 0.01.


**Table 9.** Eigenvalue of the correlation coefficient matrix and variance contribution rate.

**Table 10.** Eigenvectors of the correlation coefficient matrix.


By converting the influence of principal component factors on underground runoff into percentages, it could be known that evaporation had the largest influence on underground runoff (75%), while precipitation and catchment modifications had a significant influence on the underground runoff (−13.43% and −11.57%, respectively) with an inverse relationship.

Therefore, the contribution rate of climate change to underground runoff can be estimated to be 83.43%, while catchment modifications correspond to −11.57%.

#### 3.2.5. Summary

The partial least squares method was used to analyze the contribution rate of *Pl*, *Rls*, *El* and *Rlg* ± ε to Δ*h* according to Equation (1), and the relative contribution rate obtained was 20.86%, 29.58%, −41.28% and 9.36%, respectively.

The contribution rates of climate change and catchment variability to (Δ*h*) were obtained by the Equations (3) and (4). It was concluded that the contribution rate to the lake water level variations caused by climate and catchment factors was 93.13%, and 6.87%, respectively (Table 11).

$$R(\Delta h)\_c = 20.64\% + 40.84\% + 29.26\% \times 80.19\% + 9.26\% \times 88.44\% = 93.13\% \tag{23}$$

$$R(\Delta \text{h})\_u = 29.26\% \times 19.81\% + 9.26\% \times 11.56\% = 6.87\% \tag{24}$$


**Table 11.** Contribution rate of the hydro-climatic factors.

#### **4. Discussion**

#### *4.1. Relationship between the Hydro-Climatic Factors and Lake Water Level Variations*

In line with global warming consequences, high temperatures enhanced water vapor transport and redistribution across the entire catchment area, increasing the precipitation rates on the TP (the precipitation recorded in the basin under investigation increased by 1.4347 mm/year). Furthermore, datasets demonstrated that when temperatures increased, the potential evapotranspiration showed a decreasing trend, which confirmed the theory of the "Evaporation Paradox" [54]. Furthermore, according to datasets recorded, the annual maximum temperature and annual minimum temperatures

in the basin had risen during the period of study, while the solar duration and wind speed had significantly decreased [55,56]; therefore, these factors may have had a crucial impact on the reduction of potential evapotranspiration [57]. The evaporation rate in the lake decreased at the beginning (1960–1967), then increased (1968–1979), then declined again towards the end of the period under investigation (1980–2016). During the period from 1960 to 2004, the main reasons for the decline of water levels were the overall strong evaporation, the lack of rain and runoff, and the decrease of evaporation since 1980 could not reverse this negative trend. On the other hand, since 2005, the water level of the lake increased and this could have been due to increased precipitation recorded, combined with more runoff and lower evaporation rates.

Song et al. pointed out that runoff in most parts of the world has been decreasing significantly [58], such as in southern Australia, southern Europe, the southern region of South America and the western region of North America [59], as well as in most areas in the North of China (such as the Huaihe River [60]). While previous studies completed by Zhao et al. [61] showed that the annual runoff reduction at four main hydrologic stations in the Yellow River basin (a neighborhood area adjacent to the one investigated by this study) ranged from 17.93% to 40.79%, the results of this study showed that runoff in Qinghai Lake basin presented an upward trend, which was similar to the research of Wang et al. [62]. The results of runoff evolution attribution analysis showed that the increase in precipitation and the decrease of evaporation are the main factors leading to the increase in runoff. The trends of surface runoff and water level variations of lake were strongly consistent, and water level variations were largely affected by the effects of the climate factors. The change in precipitation had a more obvious influence on the runoff in the basins of TP, which are relatively arid, than in the humid area.

It was found that fluctuation of the annual underground runoff was not only affected by precipitation, evaporation and infiltration of surface runoff in the lake area and surrounding areas, but was also related to the fluctuation of water levels of the lake [29], and there was a noticeable connection between surface water and groundwater [63,64], showing that the runoff into the lake had positive and negative values. Since 2005, the decrease of evaporation and the increase in precipitation changed the conversion process of surface runoff and underground runoff, and the negative values increased significantly, indicating that more and more water in the lake was replenishing groundwater.

#### *4.2. Relationship between the Catchment Modifications and Water Level Variations*

In general, water level variations in the lake were the result of combined effects due to climate change and human activities. Among them, direct water intake (e.g., agricultural irrigation and drinking water for livestock) only affected the inflow rates into the lake for 4.8% of the total river discharge [28]; hence, it can be considered a negligible factor. By also developing farming areas and reducing forests (especially with local projects started in 2000), direct water intake dropped even more. Therefore, this paper did not consider the influence of direct water intake on water level changes but mainly focused on the influence of climate change and catchment change on water levels.

In those areas potentially affected by major human activities, the changes in ω were mainly manifested in land use changes and vegetation changes. It was found that 72% of the total grassland showed significant improvement in Qinghai Lake basin [65]; however, Qinghai lake basin is located at a high altitude and it is affected by a cold climate, and has low population density (4.08 people/km2), so there was little impact due to the land use changes. With the implementation of the returning pasture (farmland) to grass project since 1999 and the comprehensive management project in Qinghai Lake basin since 2008, the vegetation condition had been improved, and changes are reflected in ω trends.

According to the research conducted by Yuan [66], the annual average ground temperature in Qinghai Lake area increased by a rate of 0.74 ◦C/10 years. The depth of the annual average maximum permafrost region was then reduced by the rate of 11.7 cm/10 years and the change of permafrost layers [25,67] could definitely change the hydrologic processes under investigation. By becoming smaller, the permafrost area could not contribute consistently as previously to regulate the runoff of the catchment.

#### *4.3. Uncertainty*

The range of hydrological processes typical of great lakes is inherently uncertain, plus data scarcity adds uncertainty and methodological limitations. Firstly, this study assumed that climate and catchment were independent of each other, but the two factors were interacting in nature [68], and the effects could have also cancelled each other out. Secondly, the mathematical statistics method was used to obtain the contribution rate of climate and catchment modifications to underground runoff, but these methods had some limitations within the assumptions. Finally, the presence of permafrost complicated the investigation of hydrological processes and the characterization of their anomaly behaviors associated with climate warming.

Despite these limitations, the main purpose of this study was to use existing monitoring data to analyze the evolution law of Qinghai Lake level, separating the contribution rate of climate and catchment change to the water level variation, and better guide the future water resource management and rational utilization. From 1960 to 2016, the maximum lake area of Qinghai Lake was 4527.3 km<sup>2</sup> (in 1960), while the minimum was 4224.3 km2 (in 2004), with a difference of 303 km2, which is equivalent to the size of Co Nag Lake in China, the highest fresh water lake in the world. Therefore, this study can be very useful as a pilot case to associate with other behaviors recorded in lakes with similar and different conditions.

#### **5. Conclusions**

This study analyzed the trend of water level variation and hydro-climatic factors in Qinghai Lake Basin from 1960 to 2016 and revealed the main causes affecting the lake water levels. The paper provided a reference base for the development and management of water management in this region and provided important insights that could be applied to other basins.

Conclusions can be summarized as follows:


**Author Contributions:** J.F. and G.L. and J.Z. designed the study. J.F. performed the analysis and wrote the first draft of the manuscript. M.R., G.M., X.Y. and H.W. contributed to reviewing and editing the final version of the manuscript, software: G.J.

**Funding:** This research was supported by the National Key Research and Development Program of China (2016YFC0500802), and the Beijing Municipal Education Commission (CEFF-PXM2019\_014207\_000099), and Special Funds for Scientific Research of Forestry Public Welfare Industry (201404308), and Driving Analysis of Extreme Climate Events on Variation of Runoff and Sediment Discharge in Jinsha River Basin, the National Natural Science Foundation of China (Grant No. 41601279).

**Acknowledgments:** The authors would like to give special thanks to the Institute of Information Center of Qinghai Hydrographic Bureau, China (ICQHB) for providing data on water levels, evaporation and surface runoff.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 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 (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Modelling E**ff**ects of Rainfall Patterns on Runo**ff **Generation and Soil Erosion Processes on Slopes**

#### **Qihua Ran 1, Feng Wang <sup>1</sup> and Jihui Gao 1,2,3,\***


Received: 19 September 2019; Accepted: 22 October 2019; Published: 25 October 2019

**Abstract:** Rainfall patterns and landform characteristics are controlling factors in runoff and soil erosion processes. At a hillslope scale, there is still a lack of understanding of how rainfall temporal patterns affect these processes, especially on slopes with a wide range of gradients and length scales. Using a physically-based distributed hydrological model (InHM), these processes under different rainfall temporal patterns were simulated to illustrate this issue. Five rainfall patterns (constant, increasing, decreasing, rising-falling and falling-rising) were applied to slopes, whose gradients range from 5◦ to 40◦ and projective slope lengths range from 25 m to 200 m. The rising-falling rainfall generally had the largest total runoff and soil erosion amount; while the constant rainfall had the lowest ones when the projective slope length was less than 100 m. The critical slope of total runoff was 15◦, which was independent of rainfall pattern and slope length. However, the critical slope of soil erosion amount decreased from 35◦ to 25◦ with increasing projective slope length. The increasing rainfall had the highest peak discharge and erosion rate just at the end of the peak rainfall intensity. The peak value discharges and erosion rates of decreasing and rising-falling rainfalls were several minutes later than the peak rainfall intensity.

**Keywords:** rainfall patterns; rainfall-runoff; soil erosion; slope length; slope gradient; InHM

#### **1. Introduction**

Rainfall patterns and landform characteristics are controlling factors of the runoff and soil erosion processes in natural catchments [1,2]. Due to climatic change and climatic variability, rainfall events commonly show great temporal variation in intensity, especially in hilly areas [3,4] and the peak rainfall rates within an event may dozens of times higher than the mean event rate [1,5]. Although the temporal distribution of an individual rainfall event is diverse, some patterns of such distribution in a region can be derived based on historical data (e.g., [6,7]).

Previous studies have recognized that the rainfall patterns greatly affect the runoff generation and soil erosion processes (e.g., [8,9]). Parsons and Stone [10] adopted five rainfalls with different patterns but the same total kinetic energy to the soil surface. They found that the soil erosion amount under a constant-intensity storm are reduced by about 25% compared to varied-intensity storms, and that the eroded sediments are coarser under the constant-intensity pattern. An et al. [8] used the similar rainfall patterns and indicated that, although the total runoff was nearly not affected by the rainfall pattern, the varied intensity patterns yield 1–5 times more soil losses than even-intensity patterns and the rising pattern resulted in a consistently higher soil loss relative to the other four rainfall patterns. Conversely, Dunkerley [3] performed rainfall simulations of varying intensity profile in a dryland intergrove (runoff source area) and discovered that the late peak events showed runoff ratios that were more than double those of the early peak events and the constant rainfall yielded the lowest total

runoff, the lowest peak runoff rate. The reason was inferred to be the reductions in soil infiltration capacity during late rainfall. Zhai et al. [11] applied a distributed hydrological model at the basin scale, and found that the rainfall patterns have significant impact on the rainfall threshold of flood warning, which the flood rainfall threshold of advanced rainfall is the highest.

However, in most studies on rainfall pattern at plot scale, spatially distributed results of infiltration and soil erosion processes were not carefully considered. The temporal variation of precipitation can lead to corresponding spatial and temporal variations of infiltration, overland flow generation, and further soil erosion. Only considering runoff and soil erosion data at plot outlet, like many previous study did, will miss some important information (e.g., distributed cumulative infiltration or erosion depth) within the study area for comprehensive interpreting the influence of rainfall pattern on runoff and soil erosion processes.

In recent years, many studies on slopes have reported that observed runoff coefficient in Hortonian runoff processes decreases with increasing slope length, and variance of runoff reduces as slope scale increases (e.g., [12,13]). A reason was that the runoff generated upslope can infiltrate in downslope areas, which was called the run-on infiltration [14] or the re-infiltration [15]. Although rainfall characteristics such as duration were one of the major factors affecting runoff generation at different slope scales (e.g., [16]), it is still unknown how slope length influences the effect of temporal rainfall pattern on rainfall-runoff and soil erosion processes.

Slope steepness was an important topographic factor of hillslope rainfall-runoff and soil erosion processes. At plot scale, contradictory results were derived regarding slope effects on infiltration (e.g., [17,18]) and soil erosion (e.g., [19,20]). Besides, some researchers observed that runoff volume and soil loss on slopes increases with increasing slope angle till a critical slope angle of 20◦–30◦ (e.g., [21]), while others reported that soil erosion is not correlated with slope gradient in tilled fields (e.g., [22]). However, a majority of the studies focusing on slope steepness neglected the influence of rainfall temporal variation. There is a lack of systematically studies on the effect of slope gradient under different rainfall patterns.

Numerical modelling is an effective approach to reveal spatial and temporal impacts of rainfall patterns on infiltration, overland flow and soil erosion processes at slopes with wide ranges of steepness and length, which can broaden the limitation of the artificial rainfall experiment (e.g., a plot with a few meters long [8]). Further, strictly controlling factors such as initial condition and soil property, the effect of rainfall patterns can be specifically focused. As a mature hydraulic model, Integrated Hydrology Model (InHM) can quantitatively simulate surface (2D) and subsurface (3D) hydrologic responses to rainfall in a fully coupled approach [23,24]. Previously, InHM has been successfully applied in the simulations of hillslope hydrology and slope failure (e.g., [23,25]). As this physics-based hydrological model employs fundamental physics laws to describe natural processes [26], its output results have clear physical meanings and can be used to generalize our understanding of rainfall pattern effects on runoff and soil erosion processes.

The main purpose of this study is to investigate the impact of rainfall temporal patterns on infiltration, runoff generation and soil erosion on slopes with a range of slope lengths and gradients, using a physically based modelling approach. These modelling results are expected to improve the theoretical basis for hillslope runoff and soil erosion prediction, which will be further helpful in soil conservation planning and land management.

#### **2. Methodology**

#### *2.1. InHM Model*

The Integrated Hydrology Model (InHM) was originally developed by VanderKwaak [27], which exceeds the specifications of the hydrologic-response model proposed by Freeze and Harlan (1969) that the model based upon the numerical solution to an almost-complete set of coupled partial differential equations which describe water movement processes at surface and in unsaturated and

saturated subsurface [26]. With the advantage of the model that it doesn't need a priori assumption of a dominant runoff-generation mechanism [28,29], InHM is capable of accurately simulating dynamic infiltration, runoff and sediment processes under temporal varying rainfall. Previous studies have shown that the calibrated model reproduced accurately measured runoff and soil erosion results on semiarid hillslopes during constant-intensity rainfall-simulation events [28,30]. The equations and a detailed description of InHM can be found in VanderKwaak [27], VanderKwaak and Loague [26] and Heppner et al. [30].

#### *2.2. Model Setup*

Runoff and soil erosion processes were simulated and analysed on slopes with four horizontal projective slope lengths (25 m, 50 m, 100 m and 200 m), which were all 40 m wide and 3 m deep. For each horizontal projective length, nine slope gradients from 5◦ to 40◦ in 5◦ increments were considered, and identical rainfall amount revieved on slope surfaces was ensured for different slope gradients due to the constant horizontal projective length. In total, for each rainfall scenario the runoff and erosion processes were simulated for 36 slopes. The schematic representation of the 200 m slope used in the simulation was shown in Figure 1 as an example. To avoid the influence of the downstream outlet boundary, the overall projective slope length of the 3D finite element meshes were 220 m. The vertical nodal spacing (D*z*) in the mesh varies from 0.01 to 0.1 m; the horizontal nodal spacing (D*x*, D*y*) is 0.5 m. The 3-m mesh depth was sufficient for short-time simulation of rainfall-runoff events and deep groundwater movement was not considered in this study. The boundary conditions contain impermeable boundaries (A-E-H-D, B-F-G-C, E-F-G-H), flux boundary (A-E-F-B) and permeable boundary (A-B-C-D, C-G-H-D). The total numbers of the nodes and the elements of the 200 m slope 3D meshes are 35,721 and 70,400, respectively.

**Figure 1.** The schematic representation of the 200 m slope used the InHM modeling, and A-H represents each of the boundary nodes.

In this modelling study, the parameters of the plots (Table 1) were obtained from Ran et al. [28] who calibrated and validated the InHM parameters via the plot-scale experiments of Horton overland flow and surface erosion on the silty clay loam slopes within the Los Alamos National Laboratory [31]. The slope gradient of their experimental plot was 25.8◦ and its vegetation coverage was 61%, which was

an ideal condition for hillslope runoff and soil erosion study. A 1-h 40 mm h−<sup>1</sup> rainfall event is approximately equivalent to a 5-year return period event in that area.



#### *2.3. Rainfall Scenarios*

In this study, five temporal patterns of rainfall intensity were designed: constant rainfall intensity, increasing rainfall intensity, decreasing rainfall intensity, rising-falling rainfall intensity and falling-rising rainfall intensity (Figure 2).

**Figure 2.** Five different rainfall patterns. (**a**) constant, (**b**) decreasing, (**c**) increasing, (**d**) rising-falling and (**e**) falling-rising.

Similar rainfall scenarios were also adopted by other researchers foucusing on temporal patterns of rainfall (e.g., [3,33]). All the five rainfall patterns had a 1-h duration and 40 mm rainfall depth. Previously, many studies focusing on rainfall patterns adopted extremely high rainfall intensities in their rainfall simulations [1]. For instance, Flanagan et al. [34] used rainfall rates that peaked at 250 mm h−<sup>1</sup> and Parsons and Stone [10] adopted rainfall rates in the range 46.4–170.8 mm h−1. The rainfall intensities of this study were within the range of 1–80 mm h<sup>−</sup>1, which reprensents a more general condition.

#### **3. Results**

#### *3.1. Hydrological Responses*

The modelling results of runoff in the rainfall scenarios at the different projective slope lengths were summarized in Table 2. Generally, the increasing and rising-falling rainfalls had the largest total runoff, and the constant and falling-rising rainfalls had the least total runoff. The constant rainfall had the lowest total runoff when the slope length was shorter than 100 m. At a same slope gradient, the relative difference between the total runoff of different rainfall patterns was up to 111% (Table 2). The runoff coefficient at different projective slope lengths were shown in Figure 3. The runoff coefficient increased with increasing slope gradient until 15◦ and then decreased to the lowest value at 40◦ slope. The runoff coefficients of the rising-falling rainfalls were close to those of the increasing rainfalls at the 25 m slope, then they gradually became higher as slope length increased. As the projective slope length increased, the runoff coefficient of all the five rainfall patterns decreased, and that of the falling-rising rainfall decreased most greatly.

The peak discharge of the constant rainfall was the lowest (e.g., Figure 4). For increasing rainfall, the discharge rate kept increasing with time and reached the highest peak value among the five rainfall patterns at the end of the event. For the decreasing and rising-falling rainfalls, the peak discharges (Figure 4) were several minutes after the rainfall intensity decreased (Figure 2). The peak discharge of the falling-rising rainfall dropped a lot compared with other rainfalls as the projective slope length increased from 25 m to 200 m.


**Table 2.** Simulated results of total runoff (m3) at different projective slope lengths.

**Figure 3.** The runoff coefficient at different slopes with projective slope lengths of (**a**) 25 m, (**b**) 50 m, (**c**) 100 m and (**d**) 200 m.

**Figure 4.** The hydrographs of five different rainfall patterns with 25 m and 200 projective slope lengths at 15◦ slope (**a**,**c**) and 40◦ slope (**b**,**d**), respectively.

#### *3.2. Soil Erosion*

The soil erosion results of the rainfall scenarios at different projective slope lengths were summarized in Table 3. Similar to the hydrologic-response results, the increasing and rising-falling rainfalls had the largest total soil erosion, and the constant and falling-rising rainfalls had the least total soil erosion. This is due to the fact that rainfall-runoff was the controlling factor for soil erosion at the plot scale. Besides, at a same slope gradient, the relative difference between the soil erosion amounts of different rainfall patterns was up to 381% (Table 3). In general, relative differences in soil erosion among the five rainfall patterns were higher than relative differences in runoff under the same condition. The soil erosion amount at the different projective slope lengths were shown in Figure 5. Different from the runoff coefficient, the slope gradients that had the peak values of soil erosion amount were around 25◦–40◦, and this critical slope gradient decreased as projective slope length increased. Soil erosion amount of the rising-falling rainfalls became higher than that of the increasing rainfalls when slope length over 50 m. When the slope lengths were 25–100 m, soil erosion amount of the constant rainfalls were the lowest among the five rainfall patterns; while it became higher than that of the falling-rising rainfall when the slope length was 200 m. The sedigraphs were similar with the corresponding hydrographs. However, the erosion rate at the 40◦ slope with 25 m slope length was much higher than that at the 15◦ slope (Figure 6a,b), even the hydrographs at the two slopes only had small differences (Figure 4a,b). An example under rising-falling rainfall with 25 m slope length was shown in Figure 7.


**Table 3.** Simulated results of total soil erosion (kg) at different projective slope lengths.

**Figure 5.** The soil erosion amount at different slopes with projective slope lengths of (**a**) 25 m, (**b**) 50 m, (**c**) 100 m and (**d**) 200 m.

**Figure 6.** The sedigraphs of five different rainfall patterns with 25 m and 200 projective slope lengths at 15◦ slope (**a**,**c**) and 40◦ slope (**b**,**d**), respectively.

**Figure 7.** The comparison between the hydrographs and the sedigraphs under rising-falling rainfall at 25 m slope.

#### **4. Discussion**

#### *4.1. E*ff*ect of Rainfall Patterns on Total Runo*ff *and Soil Erosion at Di*ff*erent Slope Lengths*

For constant rainfalls, their total runoff and soil erosion were lower than those of increasing rainfalls, decreasing rainfalls and rising-falling rainfalls. It was consistent with previous studies. Dunkerley [3] observed that the runoff ratios of varying intensity rainfalls were 85–570% larger than that of constant rainfall and Wang et al. [35] found that the constant rainfall produced the lowest sediment yield at around 61.8% of the average soil loss for the increasing rainfall. Comparing the cumulative infiltration of increasing rainfalls, decreasing rainfalls and rising-falling rainfall with that of constant rainfalls (Figure 8), their gaps almost reached the highest value when slope was around 12 m, and then gradually stabilized as projective slope length increased. Thus, the differences in total runoff and soil erosion between these inconstant and constant rainfalls increased with increasing slope length. However, the total runoff and soil erosion amount of the constant rainfall become larger than those of falling-rising rainfall when the projective slope length was over 100 m (Figures 3 and 5). Previous experimental studies did not find this as their plots were much shorter than 100 m. For instance, Wang et al. [35] adopted 2 m-long flume and Parsons and Stone [10] used 2.45 m-long flume in their rainfall exerpiments.

The cumulative infiltration of the constant rainfall along the slope axis was the highest, until the projective slope length was around 50 m as the cumulative infiltration of falling-rising rainfall became larger (Figure 8). When the slope was short, compared with the constant rainfall, there was not much water infiltrated downstream the slope during the low rainfall intensity period (i.e., rainfall intensity = 1 mm h−1) under the falling-rising rainfall, because the runoff generated during the first rainfall peak quickly flowed out of the slope. Thus, the cumulative infiltration of falling-rising rainfall was lower than the constant rainfall when projective slope length was shorter than 50 m (Figure 8), leading to larger total runoff, peak discharge (Figures 3a and 4a) and erosion depth (Figure 9) than those of the constant rainfall. As slope length increased, for the falling-rising rainfall, the runoff generated during the first half of the event lasted longer before flowed out so that more water infiltrated downstream the slope during the low rainfall intensity period. Meanwhile, the recession period of the falling-rising rainfall became much longer than that of the constant rainfall (e.g., Figure 4c,d), which dramatically increased the cumulative infiltration (Figure 8). Due to these reasons, when slope length was over 100 m, the runoff coefficient of the constant rainfall became higher than that of the falling-rising rainfall (Figure 3c,d). Because of less infiltration downstream, less sediment deposited downstream the slope under the constant rainfall than the falling-rising rainfall (Figure 9), resulting in higher soil erosion amount when the slope length was over 100 m (Figure 5c,d). The results of the

falling-rising rainfall indicated that, the runoff and soil erosion amount of such multi-peak rainfall may even lower than those of the uniform rainfall, especially when the slope length was long.

**Figure 8.** The cumulative infiltration distribution along the slope axis at the 15◦ slope.

**Figure 9.** The erosion depth distribution along the slope axis at the 15◦ slope.

The total runoff and soil erosion of decreasing rainfall were much lower than those of increasing rainfall and rising-falling rainfall. It was in agreement with the experimental research on small plots by Dunkerley [3], which found the runoff ratio of late peak rainfall was double that of early peak rainfall, and numerical simulation research by Zhai et al. [11], which reported that the delayed rainfall pattern yield higher flood volume and peak than the early peak pattern. The reason was that the soil infiltrability remained high in the early part of the event under decreasing rainfall [36]. Higher cumulative infiltration of decreasing rainfall compared with those of increasing rainfall and rising-falling rainfall was obvious, especially when the projective slope length was shorter than 50 m (e.g., Figure 8). Moreover, the smaller total runoff and peak discharge under decreasing rainfall (e.g., Figure 4a,c) led to shallower erosion depth (Figure 8), due to much lower stream power and sediment transport capacity.

The simulation results also indicated that the rising-falling rainfall had the highest runoff and soil erosion amount than other rainfall patterns when projective slope length was over 50 m (Tables 2 and 3), which was not consistent with previous studies. Dunkerley [3] indicated that the late peak rainfall had the highest peak runoff rate and runoff ratio. An et al. [8] reported that in their rainfall experiments the soil loss under increasing rainfall were the highest. The main reason was that, compared with rising-falling rainfall, the increasing rainfall had much longer recession period when slope was long (e.g., Figure 4c,d), leading to larger amount of infiltration and sediment deposition.

#### *4.2. The Impact of Slope Gradient on Total Runo*ff *and Erosion under Five Rainfall Patterns*

From Figures 3 and 5 it can be seen that, for all the five different rainfall patterns, total runoff or total soil erosion showed a trend that it increased with increasing slope gradient and then gradually decreased after a critical slope. The critical slope of the total runoff was 15◦, which was independent of rainfall pattern and projective slope length. Taking the rising-falling rainfall as an example, Figure 10 shows the cumulative infiltration distribution along the slope axis at 10◦–20◦ slopes. For slopes lower than 15◦ (e.g., 10◦), overland flow velocity was slower than that on the 15◦ slope, thus leading to more infiltration, especially when the projective slope length was over 30 m (Figure 10). For slope steeper than 15◦ (e.g., 20◦), because the slope length was longer than the 15◦ slope, overland flow had to travelled longer path to reach the outlet and caused more infiltration. It can be seen in Figure 10 that cumulative infiltration difference between 15◦ and 20◦ slopes was mainly lay in area 10–50 m and 150–200 m from the slope top. Wu et al. [37] also found the critical slope for runoff rate was around 11◦ regardless of rainfall duration and slope length through a modified Green-Ampt model. The critical slope of total runoff may be affected by the surface condition (e.g., vegetation coverage, surface roughness) and the soil property (e.g., permeability, soil surface sealing), which worth further investigation.

**Figure 10.** The cumulative infiltration distribution along the slope axis under rising-falling rainfall at 10◦–20◦ slopes.

*Water* **2019**, *11*, 2221

The critical slope of soil erosion amount decreased from 35◦ to 25◦ when projective slope length increased under five different rainfall patterns (Figure 5), except for constant rainfall at 25 m slope length. Such simulation result was close to the range of critical slope of soil loss often observed in the field, which was 20◦–30◦ (e.g., [21,38]). The smaller critical slope maybe because their rainfall experiments adopted slopes with equal length. Generally, the critical slope for the constant rainfall was 5◦ larger than those of other rainfalls. Taking the rising-falling rainfall as an example, Figure 11 shows the erosion depth distribution along the slope axis at 25◦–35◦ slopes. At 35◦ slope, the erosion depth curve rose more quickly than other slopes as projective slope length increased to 25 m, due to higher flow velocity and shear stress. Thus, the critical slope of soil erosion amount was 35◦ for the slope shorter than 25 m. As mentioned above, for slope steeper than 15◦, the increase of slope gradient resulted in longer slope length and more infiltration. As the projective slope length increased, for each rainfall pattern, the reduction of runoff from 25◦ slope to 35◦ slope became larger (Figure 3) so that more sediment deposited on the slope. In consequence, the erosion depth curve at 25◦ slope finally reached the highest peak and decreased much slower than the other slopes when the projective slope length over 50 m (Figure 11). The critical slope of soil erosion amount was 25◦ for the slope longer than 100 m.

**Figure 11.** The erosion depth distribution along the slope axis under rising-falling rainfall at 25◦–35◦ slopes.

#### *4.3. E*ff*ect of Rainfall Patterns on Runo*ff *and Soil Erosion Peaks*

The time and value of the peak discharge as well as peak erosion rate were greatly influenced by the rainfall pattern. For the non-constant rainfall patterns, the increasing rainfall had the highest peak discharges and peak erosion rates, which was also mentioned in previous studies (e.g., [33]). Under increasing rainfall, as the rainfall intensity gradually increased and the surface gradually became saturated, the discharge rate and soil erosion rate kept increasing and reached the highest peak discharge and erosion rate (e.g., Figures 4 and 6) [39].

Because the infiltrability of the surface soil was high in the early part of the event, the decreasing and rising-falling rainfalls generally had lower peak discharges and peak erosion rates than the increasing rainfall (Figures 4 and 6). The peak discharge and erosion rate of increasing rainfall were reached just at the end of the peak rainfall intensity, while those of decreasing rainfall were several minutes later than the end of the peak rainfall intensity. High infiltrability of the surface soil in the early part of the event may be also the reason for the delay of the peak discharge and erosion rate, which dramatically slowed down the runoff generation process. Under rising-falling rainfall, the time

of the peak discharge and erosion rate was also later than the end of the peak rainfall intensity, but the time was shorter than that under decreasing rainfall as its peak time was during the middle of the event.

For the falling-rising rainfall, as the two high rainfall intensity periods were separated by the low-rainfall-intensity period (Figure 2), the rainfall amount for peak discharge was much less than other non-constant rainfalls. Thus, the peak discharges and erosion rates were lower than those the increasing rainfall (Figures 4 and 6). As projective slope length increased, the effect of rainfall amount was more important so that the peak discharge and erosion rate under falling-rising rainfall was even lower than those of the decreasing rainfall on the 200 m slope.

#### *4.4. Benefits and Future Work*

This research work provided comprehensive theoretical studies on effects of rainfall patterns at slope scale. Even though it lacked field measurements as validation, the parameters of the slope that used in this study were well validated previously so the simulation results were rational and realistic for runoff and sediment research. The lumped and distributed simulation results showed how rainfall patterns affected runoff generation and soil erosion processes on the wide ranges of slope gradient (5◦ to 40◦) and length (25–200 m), which can improve the accuracy of hillslope runoff and soil erosion prediction and be helpful for catchment flood management.

Table 4 illustrates the comparison between this study and previous studies, aiming at identifying the differences and emphasizing the findings of this study. In the future, the effect of rainfall patterns on hydrological responses at catchment scale will be explored. This study indicated that slope length and steepness may have great influence on the impact of rainfall patterns, and different features of hillslope length and steepness in a natural catchment will be carefully considered. Besides, rainfall patterns with multi-peak (e.g., falling-rising rainfall) showed great variety in runoff and soil erosion amount at different slope length compared with other rainfall patterns so more research effort can be put into it.


**Table 4.** Comparison of the findings between the current and previous studies.


**Table 4.** *Cont.*

Additionally, in the modelling study a stable slope surface was assumed, which meant that evolution of rill was not considered on the surface. On steep and long hillslopes, rill may generate under heavy rainfall (e.g., [40]). As the surface flow and related soil erosion characteristics in rills are different from those in an interrill area [41], it may influence the runoff generation and erosion dynamics at various scales. Although addressing the influence of rill was beyond the scope of this study, it is worth further investigations.

#### **5. Conclusions**

In this study, the effect of rainfall pattern on runoff generation and soil erosion processes on slopes were analysed through numerical modelling. The modelling work provides infiltration, runoff and soil erosion differences among five rainfall patterns on wide ranges of slope gradient (5◦ to 40◦) and slope length (25–200 m). The simulation result indicated that the rising-falling rainfall generally had the largest total runoff and soil erosion amount. The constant rainfall did not have the lowest total runoff and soil erosion amount when the projective slope length was over 100 m, which was higher than the falling-rising rainfall. The critical slope of the total runoff was 15◦, which was independent of rainfall pattern and slope length. However, the critical slope of the soil erosion amount varied, which decreased with increasing projective slope length from 35◦ to 25◦. And the critical slope for the soil erosion of the constant rainfall was generally 5◦ larger than that of other rainfalls. The increasing rainfall had the highest peak discharge and erosion rate just at the end of the peak rainfall intensity, while those of the decreasing and rising-falling rainfalls were lower and were several minutes later than the end of peak rainfall intensity.

These findings are helpful to improve the knowledge of the characteristics in runoff generation and soil erosion processes under various rainfall patterns at slopes, and they may be also beneficial for further understanding of hillslope morphology and ecology. Further work will be required for adequate meteorological and hydrological data to gain a more comprehensive understanding of rainfall pattern effects on hydrological processes at larger scale.

**Author Contributions:** Conceptualization, Q.R. and F.W.; Investigation, F.W. and J.G.; Writing—Original Draft Preparation, F.W. and J.G.; Writing—Review & Editing, Q.R. and J.G.

**Funding:** This research was funded by National Key Research and Development Program of China (2016YFC0402404) and National Natural Scientific Foundation of China (51679209).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 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 (http://creativecommons.org/licenses/by/4.0/).
