Analysis of Detailed Lake Variations and Associated Hydrologic Driving Factors in a Semi-Arid Ungauged Closed Watershed

: Lakes are key factors in maintaining ecosystems in semi-arid regions. However, due to data shortage, most studies used remote-sensing data and water-balance models to analyze lake variations in semi-arid ungauged closed watersheds, resulting in the oversimpliﬁed assessment of lake variations and their associated hydrologic processes. This study aimed to enhance the understanding of the mechanisms behind the water supplement and consumption of lakes and reveal the inﬂuences of hydrological processes on lake variations in such watersheds. Physically based and lake-oriented hydrologic modeling, remote-sensing technology, and results from previous studies were comprehensively integrated to achieve the research objective. The Hongjiannao (HJN) watershed in Northwest China was selected as the study area of this research. The calibration and validation results demonstrated that remote-sensing data and results from previous studies indeed guaranteed the accuracy of the lake-oriented model. Further hydrologic and statistical analyses revealed the linkage between lake variations and their associated hydrologic processes, and the mechanisms behind the linkage. Speciﬁcally, rainfall and snowmelt were found to be the most stable sources of HJN Lake, particularly in dry years. Due to the differences in recession rates, groundwater inﬂow was more stable than upstream inﬂow and inﬂow from the contributing area of HJN Lake. The correlations between hydrologic processes and the storage variation of HJN Lake varied signiﬁcantly at daily and monthly time scales, which can be explained by the generation mechanisms of these processes. This study provided valuable guidance for water resources management and ecosystem protection in the HJN watershed and can be further applied for hydrologic simulations in other similar watersheds.


Introduction
Lakes are essential elements for conserving water sources, purifying water quality, regulating ecological environments, and protecting biodiversity in semi-arid regions [1][2][3][4]. Due to low precipitation and intense evaporation, lakes in such regions are sensitive to environmental changes, and recovery from damage is difficult [5][6][7]. Over the past few decades, climate warming and human activities have significantly affected the water supplement and consumption of lakes in semi-arid regions, resulting in the shrinkage or disappearance of lakes, increased salinization, and changes in local ecosystems [8][9][10]. Therefore, it is crucial to understand the variations of lakes and the impacts of hydrologic processes on the variations to maintain the states of lakes and protect the fragile ecological environment in semi-arid regions [11,12].
Remote-sensing datasets, such as Landsat and MODIS data, have been widely used to extract lake areas and monitor the spatio-temporal variations of lakes [13][14][15]. Furthermore, the extracted lake areas can be combined with meteorological data (e.g., precipitation, temperature, and evaporation) to reveal the climatic driving factors of lake variations by using statistical methods [16][17][18]. Though the development of remote-sensing technology has improved data availability for natural scientific research in large-scale regions [19][20][21][22], the quality and quantity of the data are highly susceptible to climate conditions and record frequency [23,24]. Therefore, using only remote-sensing data may not be sufficient to obtain continuous daily variations of lakes and reveal the hydrological processes associated with the variations.
Integrating remote-sensing data with hydrologic modeling is an effective method to simulate and predict continuous hydrologic processes under various climate conditions [25][26][27]. Compared with remote-sensing data alone, the combined method can depict the internal hydrologic processes between climate and lakes (e.g., surface runoff, groundwater flow, and snowmelt), and quantify the impacts of these processes on lake variations [28][29][30]. Many studies have been conducted to simulate the daily variations of lakes and their associated hydrologic processes with long time periods in gauged watersheds [31][32][33]. However, due to the lack of observed streamflow, most studies in semi-arid and ungauged watersheds utilized simple water-balance models to simulate the water budgets and lake areas and analyze the hydrologic driving factors of lake variations at monthly or annual time scales [34][35][36][37]. These models cannot reveal detailed water-exchange procedures due to their simplified structures and large simulation time scales [38,39]. Therefore, studies using physically based hydrologic models with fine time scales are essential for promoting the simulation qualities in semi-arid and ungauged lake watersheds.
The selection of model evaluation data is critical for applying physically based models in semi-arid and ungauged watersheds. In some gauged watersheds, observed streamflow data and remotely sensed data have been utilized for model calibration and validation [40][41][42], and their combination indeed improved the performance of hydrologic models. For ungauged watersheds, several studies have verified the possibilities of using remote-sensing data (e.g., evapotranspiration, leaf area index, and soil moisture data) for model evaluation [43,44], and their results demonstrated that, without observed streamflow, the model calibrated by remotely sensed data can still simulate hydrologic processes appropriately. Overall, the above studies provide a feasible option for the evaluation of hydrologic models in data-poor watersheds and make the applications of the physically based models to semi-arid and ungauged closed watersheds possible.
The objective of this study was to establish a comprehensive approach to promote the understanding of lake supplement and consumption mechanisms and reveal the influences of hydrological processes on lake variations in semi-arid ungauged closed watersheds. Unlike previous studies that used simple water-balance models, the new approach adopted a modified lake-oriented hydrologic model to simulate hydrologic processes related to lake variations in closed watersheds, and utilized multiple remote-sensed data and results from previous studies to set up and evaluate the model. The new approach was further applied to the Hongjiannao (HJN) watershed in Northwest China, where there are no meteorological and stream gauge stations, to: (1) simulate detailed lake variations and their associated hydrologic processes with limited data; (2) reveal the linkage between lake variations and hydrologic processes at different time scales; and (3) interpret the linkage based on the generation mechanisms of lake-related hydrologic processes.

Study Area
HJN Lake is an irregular triangle-shaped lake that is located at the junction of Shaanxi and Inner Mongolia provinces in Northwest China. It is the largest freshwater desert lake in China and the largest habitat and breeding ground for relict gulls globally [45,46]. The lake was formed in 1929 with an area of 1.3 km 2 , which expanded to 60-67 km 2 during the 1960-1970s and gradually decreased to 35 km 2 in 2018 [45,46]. Currently, the east-west width, north-south length, and depth of HJN Lake are approximately 8 km, 11.2 km, and 8.2 m, respectively [45]. Before 2000, HJN Lake was supplied by seven seasonal rivers, which were the Zhasake River, Manggaitu River, Qibosu River, Mudushili River, Erlintu River, Songdaogou River, and Miaohao River [46] (Figure 1a). However, due to climate change and human activities, the Miaohao River, Erlintu River, and Songdaogou River have stopped supplying water to HJN Lake. Furthermore, in 2005, the construction of the Zhasake reservoir blocked the upstream section of the Zhasake River. Our field survey indicates that, in recent years, the Manggaitu River and Qibosu River have streamflow during dry periods, while the Zhasake River and Mudushili River have streamflow after rainfall events.

Study Area
HJN Lake is an irregular triangle-shaped lake that is located at the junction of Shaanxi and Inner Mongolia provinces in Northwest China. It is the largest freshwater desert lake in China and the largest habitat and breeding ground for relict gulls globally [45,46]. The lake was formed in 1929 with an area of 1.3 km 2 , which expanded to 60-67 km 2 during the 1960-1970s and gradually decreased to 35 km 2 in 2018 [45,46]. Currently, the east-west width, north-south length, and depth of HJN Lake are approximately 8 km, 11.2 km, and 8.2 m, respectively [45]. Before 2000, HJN Lake was supplied by seven seasonal rivers, which were the Zhasake River, Manggaitu River, Qibosu River, Mudushili River, Erlintu River, Songdaogou River, and Miaohao River [46] (Figure 1a). However, due to climate change and human activities, the Miaohao River, Erlintu River, and Songdaogou River have stopped supplying water to HJN Lake. Furthermore, in 2005, the construction of the Zhasake reservoir blocked the upstream section of the Zhasake River. Our field survey indicates that, in recent years, the Manggaitu River and Qibosu River have streamflow during dry periods, while the Zhasake River and Mudushili River have streamflow after rainfall events. The HJN watershed is a closed lake watershed with a total area of 1343.79 km 2 , and an altitude of 1210 to 1518 m. There are no meteorological and hydrological stations in the watershed. Based on the record data of nearby meteorological stations, the HJN watershed receives most of the precipitation in summer, with annual precipitation ranging from 108.6 mm to 819.1 mm and annual temperature ranging from −28.1 ℃ to 38.9 ℃. In this study, ArcGIS Pro 2.8 and a surface delineation algorithm, HUD-DC [47], were used to identify rivers, HJN Lake, the Zhasake reservoir, and the boundaries of each sub-basin ( Figure  1b) based on the ALOS Global Digital Surface Model (https://www.eorc.jaxa.jp/ALOS/en/dataset/aw3d30/aw3d30_e.htm (accessed on 16 May 2022)). The HJN watershed is a closed lake watershed with a total area of 1343.79 km 2 , and an altitude of 1210 to 1518 m. There are no meteorological and hydrological stations in the watershed. Based on the record data of nearby meteorological stations, the HJN watershed receives most of the precipitation in summer, with annual precipitation ranging from 108.6 mm to 819.1 mm and annual temperature ranging from −28.1°C to 38.9°C. In this study, ArcGIS Pro 2.8 and a surface delineation algorithm, HUD-DC [47], were used to identify rivers, HJN Lake, the Zhasake reservoir, and the boundaries of each sub-basin ( Figure 1b) based on the ALOS Global Digital Surface Model (https://www.eorc.jaxa.jp/ ALOS/en/dataset/aw3d30/aw3d30_e.htm (accessed on 16 May 2022)).

Research Approach
To accurately simulate the hydrologic processes and variations of HJN Lake, it is necessary to select a physically based model suitable for simulating the drainage from multiple rivers to a lake. The selected model should be able to simulate traditional hydrologic processes (e.g., surface runoff and runoff routing), and the variations of lakes and the hydrologic processes directly related to the lakes (e.g., lake evaporation and inflows from rivers and groundwater). Additionally, the model must be set up and evaluated based on limited available data. Input data of the model include meteorological data, land surface data, soil data, and groundwater data, which can be obtained from remote-sensing datasets and previous studies in the absence of meteorological and hydrological stations. To ensure model performance, key hydrologic processes must be evaluated using available data. According to the hydrologic characteristics of the HJN watershed, this study selected four valuation variables for model evaluation, including daily evapotranspiration (ET), average annual surface runoff of all sub-basins, average annual inflow from groundwater to HJN Lake, and the area of HJN Lake. Figure 2 outlines the detailed procedures of the proposed approach. Remote-sensing data were processed to generate meteorological data, land surface data, and soil data required for the selected model, and groundwater data were obtained from previous studies. Then, the selected model was used to calculate hydrologic processes and the variations of HJN Lake. The simulated values of the evaluation variables were compared with their observed or reference values, and the input parameters of the model were adjusted until all evaluation variables were well-fitted.
To accurately simulate the hydrologic processes and variations of HJN Lake, it is necessary to select a physically based model suitable for simulating the drainage from multiple rivers to a lake. The selected model should be able to simulate traditional hydrologic processes (e.g., surface runoff and runoff routing), and the variations of lakes and the hydrologic processes directly related to the lakes (e.g., lake evaporation and inflows from rivers and groundwater). Additionally, the model must be set up and evaluated based on limited available data. Input data of the model include meteorological data, land surface data, soil data, and groundwater data, which can be obtained from remote-sensing datasets and previous studies in the absence of meteorological and hydrological stations. To ensure model performance, key hydrologic processes must be evaluated using available data. According to the hydrologic characteristics of the HJN watershed, this study selected four valuation variables for model evaluation, including daily evapotranspiration (ET), average annual surface runoff of all sub-basins, average annual inflow from groundwater to HJN Lake, and the area of HJN Lake. Figure 2 outlines the detailed procedures of the proposed approach. Remote-sensing data were processed to generate meteorological data, land surface data, and soil data required for the selected model, and groundwater data were obtained from previous studies. Then, the selected model was used to calculate hydrologic processes and the variations of HJN Lake. The simulated values of the evaluation variables were compared with their observed or reference values, and the input parameters of the model were adjusted until all evaluation variables were well-fitted.

Hydrologic Modeling
HYDROL-DC is a daily and semi-distributed hydrologic model that utilizes a physically based approach to simulate hydrologic processes related to lakes and rivers at watershed scales [48]. As illustrated in Figure 3a, HYDROL-DC divides a watershed into several sub-basins, and each sub-basin can be further divided into one main channel and several puddle-based units (PBUs) and channel-based units (CBUs). CBUs are further classified into on-stream CBUs, which interact directly with the main channel, and off-stream CBUs, which have no direct interaction. Each PBU is composed of a ponding area and its Sustainability 2023, 15, 6535 5 of 20 associated contribution area, and each CBU is composed of a channel segment (a main channel or tributary channel) and its associated contribution area. Surface runoff in a PBU transfers from the contribution area to the ponding area and is eventually released to its downstream unit when the storage capacity of the ponding area is exceeded. Surface runoff in a CBU runs directly into its tributary channel or the main channel.
HYDROL-DC is a daily and semi-distributed hydrologic model that utilizes a physically based approach to simulate hydrologic processes related to lakes and rivers at watershed scales [48]. As illustrated in Figure 3a, HYDROL-DC divides a watershed into several sub-basins, and each sub-basin can be further divided into one main channel and several puddle-based units (PBUs) and channel-based units (CBUs). CBUs are further classified into on-stream CBUs, which interact directly with the main channel, and offstream CBUs, which have no direct interaction. Each PBU is composed of a ponding area and its associated contribution area, and each CBU is composed of a channel segment (a main channel or tributary channel) and its associated contribution area. Surface runoff in a PBU transfers from the contribution area to the ponding area and is eventually released to its downstream unit when the storage capacity of the ponding area is exceeded. Surface runoff in a CBU runs directly into its tributary channel or the main channel. Vertically, each unit is divided into five zones including the canopy zone, snow zone, surface zone, upper soil zone, and lower soil zone. Additionally, all units in a sub-basin share the same shallow groundwater zone. The upper soil zone can be further divided into multiple cells to simulate the soil water content at different depths, whereas the lower soil zone is lumped together with a reservoir function controlling percolation to the shallow groundwater zone. The major processes in each zone are listed in Table 1. To account for the closed characteristic of the study area, HYDROL-DC was modified by allowing a lake to serve as the outlet of a sub-basin ( Figure 3b). In this study, the HJN watershed was divided into 8 sub-basins and each sub-basin contains only one unit.  Vertically, each unit is divided into five zones including the canopy zone, snow zone, surface zone, upper soil zone, and lower soil zone. Additionally, all units in a sub-basin share the same shallow groundwater zone. The upper soil zone can be further divided into multiple cells to simulate the soil water content at different depths, whereas the lower soil zone is lumped together with a reservoir function controlling percolation to the shallow groundwater zone. The major processes in each zone are listed in Table 1. To account for the closed characteristic of the study area, HYDROL-DC was modified by allowing a lake to serve as the outlet of a sub-basin ( Figure 3b). In this study, the HJN watershed was divided into 8 sub-basins and each sub-basin contains only one unit.  Table 2 lists the major input data utilized in HYDROL-DC, which include topographic data (e.g., sub-basin area and average elevation, the maximum ponding area and storage of each lake, and the length and slope of each river), precipitation, air temperature, leaf area index, potential ET, SCS curve number, and soil properties (e.g., saturated hydraulic conductivity, wilting point, and field capacity). The topographic data were calculated using ArcGIS based on the ALOS Global Digital Surface Model. The daily precipitation data were obtained from the China Meteorological Forcing Dataset [49] to calculate the average daily precipitation of each sub-basin. The leaf area index and air temperature were downloaded from the Global Land Surface Satellite dataset [50]. The potential ET of each sub-basin was estimated using the Hargreaves method, which considered the minimum and maximum temperatures [51]. The curve numbers of all sub-basins were calculated based on the Global Curve Number data [52]. Soil properties of all sub-basins were calculated using the Saxton Equation [53] based on the percentages of clay and sand recorded in the Harmonized World Soil Database v1.2 [54]. Groundwater routing was handled using a linear reservoir function, and the parameters of the groundwater zone were determined based on Yang et al. (2019) [55].

Evaluation Data
The daily observed ET of each sub-basin was calculated using a Harmonized Global Land Evaporation dataset from model-based products covering 1980-2017 [56]. To determine lake areas, Landsat-8 data was downloaded from the Geospatial Data Cloud (http://www.gscloud.cn/search (accessed on 23 August 2022)) and the Normalized Difference Water Index (NDWI) method [57] was applied to calculate lake areas based on the band 3 and band 6 (i.e., green light and near-infrared bands) of the Landsat-8 data (Equation (1)). Table 3 lists all the lake areas used in this study and their imaging dates.
where Green and N IR are the reflectance of green light and near-infrared bands, respectively. Yang (2014) [58] calculated the natural average annual surface runoff and the consumptions caused by land cover change, climate change, reservoir construction, and domestic usage in the HJN watershed. In this study, the actual average annual surface runoff of each sub-basin was calculated based on the above values (Table 4). In addition, Yang (2014) [58] estimated the average annual groundwater discharge, which was 649 × 10 4 m 3 /a.

Model Setup and Performance Evaluation
The construction of the Zhasake reservoir has a significant impact on the inflow from the Zhasake River to HJN Lake; therefore, we selected the data from  [59,60] was utilized for sensitivity analysis. Since this study focused on lake variations, the simulation output in the sensitivity analysis was specified as the mean of lake areas, and the percentage adjustment of each calibrated input parameter was set to 10%.
where ∆P is the percentage adjustment of an input parameter; S 0 , S 1 and S 2 are the simulation outputs when the input parameter equals P 0 , P 0 − ∆P and P 0 + ∆P, respectively. After the sensitivity analysis, three statistical parameters, including the Nash-Sutcliffe efficiency (NSE) coefficient, percent bias (PBIAS), and Pearson correlation coefficient (r), were used to evaluate the performance of HYDROL-DC in simulating ET of all sub-basins and the area of HJN Lake based on the performance ratings summarized by Moriasi et al. (2007) [61]. The NSE, PBIAS, and r can, respectively, be expressed as: where O i is the value of the i th observed data, S i is the value of the i th simulated data, O is the mean value of the observed data, and N is the dataset size.

Analyses of Lake Variation and Impact Factors
The hydrologic processes that directly control the variations of HJN Lake include: (1) rainfall and snowmelt over the lake surface; (2) inflow from the contribution area of HJN Lake (i.e., CA inflow); (3) inflow from upstream rivers (i.e., upstream inflow); (4) the groundwater discharge to the lake (i.e., groundwater inflow); and (5) evaporation from the lake surface. In this study, the daily and annual variation characteristics of these processes were analyzed to reveal the mechanisms behind the dynamic changes in HJN Lake. The ratio of the annual inflow from each river to the annual upstream inflow of HJN Lake was calculated to quantify the contribution of each river to the lake variations. For each river, the relative standard deviation of the ratios in different years was calculated to analyze the stability of the inflow from each river to HJN Lake (Equation (6)). In addition, to access the impact of each hydrologic process on lake variations, the entire simulation period was divided into an increase stage (when the lake storage increased) and a decrease stage (when the lake storage decreased), and the linear regression method and Pearson correlation coefficient (Equation (5)) were used to analyze the correlations between each hydrologic process and the variations of lake storage at daily and monthly time scales.
where RSD is the relative standard deviation, x i is the value of the i th sample point, x is the mean value of a sample, N is the sample size. Table 5 lists the calibrated parameters and their sensitivity index values (i.e., I). The curve number and lake evaporation coefficient were extremely sensitive to the area of HJN Lake, followed by the shallow groundwater percolation coefficient, initial abstraction, and an empirical evapotranspiration parameter (i.e., C2 et ). According to the sensitivity analysis, surface runoff and lake evaporation had significant influences on the variations of HJN Lake, and the shallow groundwater percolation coefficient impacted the groundwater inflow of the lake. Thus, these parameters were mainly adjusted in the calibration process. In addition, the other less sensitive parameters in Table 5 were also calibrated to improve the simulation results.  Figure 4 compares the simulated and observed ET for all sub-basins. Throughout the simulation period, the simulated ET shows a similar trend to the observed data. The coefficient of correlation (r) ranged from 0.82 to 0.87 in the calibration period and from 0.79 to 0.83 in the validation period, indicating a high degree of agreement between the simulated and observed trends in the variations of daily ET. Due to the underestimation of potential ET by using the Hargraves equation, the simulated ET from July to October was lower than the observed values. The ET of sub-basin eight was overestimated from April to June, which can be attributed to the constant lake evaporation coefficient used in HYDROL-DC. In general, HYDROL-DC slightly overestimated the ET of sub-basin eight and underestimated the ET of other sub-basins. The PBIAS values ranged from −8.52% to +18.21% in the calibration period and from −2.08% to 23.61% in the validation period, indicating the deviation between the simulated and observed ET. The ranges of NSE values were 0.63-0.74 and 0.56-0.67 for the calibration and validation periods, respectively. According to Moriasi et al. (2007) [61], the simulation results were rated between satisfactory and good.

Evaluation of Model Performance
The average annual surface runoff of each sub-basin was calculated and compared with the reference value from Yang (2014) [59] (Table 6). The simulation error ranged from −9.54% to +9.64%, with the largest error in sub-basin eight and the smallest in sub-basin six. In general, HYDROL-DC slightly overestimated the surface runoff in most sub-basins, except for sub-basins four and eight. For groundwater evaluation, the simulated average annual groundwater discharge was 605.13 × 10 4 m 3 , with a simulation error of −6.67% compared to the reference value (i.e., 649 × 10 4 m 3 ). Overall, HYDROL-DC successfully simulated the annual surface runoff for all sub-basins and the groundwater discharge to HJN Lake. underestimated the ET of other sub-basins. The PBIAS values ranged from −8.52% to +18.21% in the calibration period and from −2.08% to 23.61% in the validation period, indicating the deviation between the simulated and observed ET. The ranges of NSE values were 0.63-0.74 and 0.56-0.67 for the calibration and validation periods, respectively. According to Moriasi et al. (2007) [61], the simulation results were rated between satisfactory and good. The average annual surface runoff of each sub-basin was calculated and compared with the reference value from Yang (2014) [59] (Table 6). The simulation error ranged from −9.54% to +9.64%, with the largest error in sub-basin eight and the smallest in sub-basin six. In general, HYDROL-DC slightly overestimated the surface runoff in most sub-basins, except for sub-basins four and eight. For groundwater evaluation, the simulated average annual groundwater discharge was 605.13 × 10 4 m 3 , with a simulation error of −6.67% compared to the reference value (i.e., 649 × 10 4 m 3 ). Overall, HYDROL-DC successfully simulated the annual surface runoff for all sub-basins and the groundwater discharge to HJN Lake.   Figure 5 shows the simulated and observed areas of HJN Lake in the calibration and validation periods. The values of r, NSE, and PBIAS in the calibration period were 0.93, 0.71, and +1.89%, respectively, and those values in the validation period were 0.96, 0.65, and +1.11%, respectively. As shown in Figure 5, the variations of the simulated lake areas followed a similar trend to the observed areas with some underestimation in dry years (e.g., 2012 and 2015) due to the relatively simple groundwater module of the HYDROL-DC model. The high r values indicated a close similarity between the simulated and observed trends in the variations of the lake area. Furthermore, according to Moriasi et al. (2007) [61], the values of NSE and PBIAS were rated between good and very good for both the calibration and validation periods.
In summary, the r values were all greater than 0.79, indicating that the simulated ET and lake area had significant linear correlations with their observed values, respectively. The values of NSE mostly ranged from 0.65 to 0.74, most of the PBIAS values were within ±20%, and all the errors were smaller than ±10% during the entire simulation period. Therefore, according to Moriasi et al. (2007) [61], the overall performance of the HYDROL-DC in simulating hydrologic processes and lake areas can be rated as good. Considering the data shortage in the HJN watershed, the proposed approach yielded reasonable simulation results.

Hydrologic Process Analysis
The daily rainfall and snowmelt, evaporation, and lake storage are shown in Figure 6a, and Figure 6b depicts the daily CA inflow, upstream inflow, groundwater inflow, and lake area. Both figures were combined to analyze the dynamic characteristics of all hydrologic processes and the influence of these processes on lake variations. According to Figure 6, rainfall was concentrated from July to September, and the CA inflow, upstream inflow, and groundwater inflow were mainly generated after rainfall events. Additionally, inflows also could be generated in early spring (e.g., March 2016) and the end of fall (e.g., November 2015) due to snowmelt. In terms of the recession rate, the groundwater inflow had the lowest rate, followed by upstream inflow and CA inflow (Figure 6b). The storage and area of HJN Lake continued to increase after rainfall and snowmelt, and this increase continued into winter if there was sufficient supplementation (e.g., in the winter of 2013). The reduction in the storage and area typically occurred in spring and early summer due to increased evaporation (Figure 6a) and inadequate water supply (Figure 6b).  Figure 5 shows the simulated and observed areas of HJN Lake in the calibration and validation periods. The values of r, NSE, and PBIAS in the calibration period were 0.93, 0.71, and +1.89%, respectively, and those values in the validation period were 0.96, 0.65, and +1.11%, respectively. As shown in Figure 5, the variations of the simulated lake areas followed a similar trend to the observed areas with some underestimation in dry years (e.g., 2012 and 2015) due to the relatively simple groundwater module of the HYDROL-DC model. The high r values indicated a close similarity between the simulated and observed trends in the variations of the lake area. Furthermore, according to Moriasi et al. (2007) [61], the values of NSE and PBIAS were rated between good and very good for both the calibration and validation periods. In summary, the r values were all greater than 0.79, indicating that the simulated ET and lake area had significant linear correlations with their observed values, respectively. The values of NSE mostly ranged from 0.65 to 0.74, most of the PBIAS values were within ±20%, and all the errors were smaller than ±10% during the entire simulation period. Therefore, according to Moriasi et al. (2007) [61], the overall performance of the HYDROL-DC in simulating hydrologic processes and lake areas can be rated as good. Considering the data shortage in the HJN watershed, the proposed approach yielded reasonable simulation results.   Table 7 show the annual water balance of HJN Lake. The source terms of HJN Lake included rainfall and snowmelt, CA inflow, upstream inflow, and groundwater inflow. Evaporation was the sink term of HJN Lake, with an average value of 125 mm. Based on the storage variations of HJN Lake, the entire simulation period can be divided into a decrease period (2008-2011), a fluctuation period (2012)(2013)(2014)(2015), and an increase period (2016-2017). During the decrease period, the evaporation was greater than its average value, and lake storage decreased by 346.84 mm due to the significant evaporation and insufficient water supply. Rainfall and snowmelt were the main factors maintaining the lake storage during this period. As the drought continued, the proportion of rainfall and snowmelt in the total water supply of HJN Lake gradually increased from 52.54% to 95.25%. In the fluctuation period, due to the changes in the total water supply and evaporation, lake storage increased by 69.88 mm from 2012 to 2013 and then decreased by 124.32 mm from 2014 to 2015. The proportion of each source term in the total supply changed alternately, with the variation in groundwater inflow lagging behind other source terms. In the increase period, lake storage increased by 256.49 mm. Compared with rainfall and snowmelt, surface and groundwater inflows contributed the most water to HJN Lake. In 2016 and 2017, the proportion of the inflows in the total water supply was 82.80% and 78.83%, respectively.
The relative standard deviations of the annual rainfall and snowmelt, CA inflow, upstream inflow, and groundwater inflow were 0.33, 1.41, 1.12, and 1.10, respectively. Therefore, rainfall and snowmelt provided the most stable supplement to HJN Lake, followed by groundwater inflow, upstream inflow, and CA inflow. lake area. Both figures were combined to analyze the dynamic characteristics of all hydrologic processes and the influence of these processes on lake variations. According to Figure 6, rainfall was concentrated from July to September, and the CA inflow, upstream inflow, and groundwater inflow were mainly generated after rainfall events. Additionally, inflows also could be generated in early spring (e.g., March 2016) and the end of fall (e.g., November 2015) due to snowmelt. In terms of the recession rate, the groundwater inflow had the lowest rate, followed by upstream inflow and CA inflow (Figure 6b). The storage and area of HJN Lake continued to increase after rainfall and snowmelt, and this increase continued into winter if there was sufficient supplementation (e.g., in the winter of 2013). The reduction in the storage and area typically occurred in spring and early summer due to increased evaporation ( Figure 6a) and inadequate water supply (Figure 6b).   Table 7 show the annual water balance of HJN Lake. The source terms of HJN Lake included rainfall and snowmelt, CA inflow, upstream inflow, and groundwater inflow. Evaporation was the sink term of HJN Lake, with an average value of 125 mm. Based on the storage variations of HJN Lake, the entire simulation period can be divided into a decrease period (2008-2011), a fluctuation period (2012-2015), and an increase period (2016-2017). During the decrease period, the evaporation was greater than its average value, and lake storage decreased by 346.84 mm due to the significant evaporation and insufficient water supply. Rainfall and snowmelt were the main factors maintaining the lake storage during this period. As the drought continued, the proportion of rainfall and snowmelt in the total water supply of HJN Lake gradually increased from 52.54% to 95.25%. In the fluctuation period, due to the changes in the total water supply The proportion of each source term in the supply changed alternately, with the variation in groundwater inflow lagging be other source terms. In the increase period, lake storage increased by 256.49 mm. pared with rainfall and snowmelt, surface and groundwater inflows contributed the water to HJN Lake. In 2016 and 2017, the proportion of the inflows in the total water ply was 82.80% and 78.83%, respectively. The relative standard deviations of the an rainfall and snowmelt, CA inflow, upstream inflow, and groundwater inflow were 1.41, 1.12, and 1.10, respectively. Therefore, rainfall and snowmelt provided the mos ble supplement to HJN Lake, followed by groundwater inflow, upstream inflow, an inflow. Figure 7. Annual water balance of HJN Lake (RSIN = rainfall and snowmelt over the lake su EVA = evaporation from the lake surface; CAIN = inflow from the contributing area of HJN GWIN = inflow from the groundwater zone; USIN = inflow from the upstream rivers). Table 7. Annual water balance table of Hongjiannao Lake (RSIN = rainfall and snowmelt ov lake surface; EVA = evaporation from the lake surface; CAIN = inflow from the contributing a Hongjiannao Lake; GWIN = inflow from the groundwater zone; USIN = inflow from the ups rivers; DS = lake storage variation).

Figure 7.
Annual water balance of HJN Lake (RSIN = rainfall and snowmelt over the lake surface; EVA = evaporation from the lake surface; CAIN = inflow from the contributing area of HJN Lake; GWIN = inflow from the groundwater zone; USIN = inflow from the upstream rivers). Table 7. Annual water balance table of Hongjiannao Lake (RSIN = rainfall and snowmelt over the lake surface; EVA = evaporation from the lake surface; CAIN = inflow from the contributing area of Hongjiannao Lake; GWIN = inflow from the groundwater zone; USIN = inflow from the upstream rivers; DS = lake storage variation).  Figure 8 shows the ratio of the annual inflow from each river to the annual upstream inflow of HJN Lake. Based on the surface delineation results shown in Figure 1, the contribution areas of Mudushili River and Qibosu River are sub-basins three and six, respectively. The contributing area of the Zhasake River is composed of sub-basins one, two, and seven, and the contribution area of the Manggaitu River is composed of subbasins four and five. The ratio of the inflow from the Manggaitu River to the upstream inflow ranged from 0.41 to 0.50, while for the Qibosu River, Zhasake River, and Mudushili River, the ranges were 0.30-0.54, 0.02-0.16, and 0.01-0.07, respectively. The average ratio of the inflow from the Manggaitu River to the upstream inflow was 0.45, followed by the Qibosu River (0.44), Zhasake River (0.07), and Mudushili River (0.04). Thus, the Manggaitu River and Qibosu River contributed the most surface runoff to HJN Lake throughout the simulation period. Furthermore, the supply stability of each river was quantified by calculating the relative standard deviation of the ratios in different years. Accordingly, the Manggaitu River was the most stable (RSD = 0.06), followed by the Qibosu River (RSD = 0.19), Mudushili River (RSD = 0.58), and Zhasake River (RSD = 0.73).

Sustainability 2023, 15, x FOR PEER REVIEW
14 of 20 Figure 8 shows the ratio of the annual inflow from each river to the annual upstream inflow of HJN Lake. Based on the surface delineation results shown in Figure 1, the contribution areas of Mudushili River and Qibosu River are sub-basins three and six, respectively. The contributing area of the Zhasake River is composed of sub-basins one, two, and seven, and the contribution area of the Manggaitu River is composed of sub-basins four and five. The ratio of the inflow from the Manggaitu River to the upstream inflow ranged from 0.41 to 0.50, while for the Qibosu River, Zhasake River, and Mudushili River, the ranges were 0.30-0.54, 0.02-0.16, and 0.01-0.07, respectively. The average ratio of the inflow from the Manggaitu River to the upstream inflow was 0.45, followed by the Qibosu River (0.44), Zhasake River (0.07), and Mudushili River (0.04). Thus, the Manggaitu River and Qibosu River contributed the most surface runoff to HJN Lake throughout the simulation period. Furthermore, the supply stability of each river was quantified by calculating the relative standard deviation of the ratios in different years. Accordingly, the Manggaitu River was the most stable (RSD = 0.06), followed by the Qibosu River (RSD = 0.19), Mudushili River (RSD = 0.58), and Zhasake River (RSD = 0.73).

Correlation Analysis
The correlations between the storage variation of HJN Lake and the water balance terms at monthly and daily time scales are shown in Figure 9. The storage variation was categorized into two conditions: an increase condition (variation greater than zero) and a decrease condition (variation less than zero). At the monthly time scale, the two major factors that increased the lake storage were upstream inflow (r = 0.97) and CA inflow (r = 0.95) (Figure 9a). In contrast, rainfall and snowmelt had a weaker correlation with the

Correlation Analysis
The correlations between the storage variation of HJN Lake and the water balance terms at monthly and daily time scales are shown in Figure 9. The storage variation was categorized into two conditions: an increase condition (variation greater than zero) and a decrease condition (variation less than zero). At the monthly time scale, the two major factors that increased the lake storage were upstream inflow (r = 0.97) and CA inflow (r = 0.95) (Figure 9a). In contrast, rainfall and snowmelt had a weaker correlation with the storage variation (r = 0.56). Evaporation was found to be the factor causing the decrease in lake storage (Figure 9b) and had a strong negative linear relationship with the variation in lake storage (r = −0.87). At the daily time scale, rainfall and snowmelt were the most significant factors that increased the lake storage (r = 0.80), followed by the CA inflow (r = 0.64) and upstream inflow (r = 0.58) (Figure 9c). Similar to the results at the monthly time scale, evaporation had a strong negative linear relationship with the variation in the lake storage at the daily time scale, with a coefficient of correlation of −0.89 (Figure 9d). The groundwater inflow had little correlation with the storage variation, which might be due to the small quantity and the low recession rate of the groundwater inflow.

Criteria and Results of Model Evaluation
This study aimed to solve the contradiction between data shortage and the requirement for detailed hydrological information in semi-arid and ungauged closed watersheds. Previous studies have confirmed the effectiveness of utilizing remote-sensed ET data for model evaluation in such areas [42][43][44]. Based on the hydrological characteristics of lakes in semi-arid regions [13,30,58], in addition to the remote-sensed ET data, we also used remotely sensed lake areas and the results from previous studies for model evaluation. The calibration and validation results indicated that the simulated ET was between satisfactory and good in general and underestimation mainly occurred from July to October. This phenomenon has been observed in previous studies and can be attributed to the underestimation of potential ET in semi-aired regions by using the Hargraves equation [62][63][64]. Throughout the simulation period, the simulated lake area corresponded well with the observed values. However, due to the simplification of groundwater, the simulated lake area decreased faster than the observed values in some dry years (e.g., 2011 and 2015). HYDROL-DC employs a linear reservoir function to simulate groundwater discharge based on soil percolation, resulting in possible underestimation of discharge when there is no sufficient replenishment. Previous studies using similar functions had the same issue in simulating groundwater discharge under dry climate conditions [65][66][67].

Comparison with Previous Studies
Most previous studies in semi-arid and ungauged closed watersheds used remotesensing data and water-balance models to analyze monthly and annual lake variations [25][26][27][28][29][30]. This study provides a valuable supplement to these previous studies and can address the gaps in simulating detailed hydrologic processes related to lakes. In the HJN watershed, previous studies identified precipitation and evaporation (or temperature) as the meteorological driving factors influencing lake variations [68][69][70], and surface runoff and evaporation were identified as the hydrologic driving factors [30]. These results were consistent with the annual mass balance in our study. In addition, our study utilized the daily simulation results to reveal the linkage between hydrologic processes and lake variations: (1) rainfall and snowmelt were the most stable supplements to HJN Lake due to their instantaneous supply mechanism; (2) compared with rainfall and snowmelt, the generation of groundwater inflow, upstream inflow, and CA inflow were relatively more susceptible to climate and soil conditions, and therefore had greater effects on the water supplement of HJN Lake in wet years rather than dry years (Figure 7 and Table 7); (3) the stability of groundwater inflow was followed by upstream inflow and CA inflow due to the differences in their recession rates ( Figure 6).
Liang and Li (2019) [69] reported that precipitation and ET had a positive and negative correlation, respectively, with the area of HJN Lake. The coefficient of correlation varied at monthly, seasonal, and annual time scales. Our study yielded similar results, which can be attributed to the generation mechanisms of hydrologic processes. The correlation analysis showed that evaporation dominated the decrease of lake storage (Figure 9b,d), and the coefficient of correlation was slightly greater at the daily time scale than at the monthly time scale due to the instantaneous generation mechanism of lake evaporation. At the daily time scale, rainfall and snowmelt showed a stronger correlation with the increase of lake storage compared to the CA inflow and upstream inflow, while the correlation was weaker at the monthly time scale. The upstream inflow showed a stronger correlation with the increase of lake storage compared to the CA inflow, whereas the correlation was weaker at the daily time scale. These findings can be explained by the following reasons: (1) rainfall and snowmelt over the lake have an immediate impact on increasing lake storage; (2) the CA inflow originates in the contribution area of HJN Lake and supplies the lake after surface runoff routing; and (3) in comparison to the CA inflow, the upstream inflow requires channel routing before supplying HJN Lake, in addition to the runoff generation and surface runoff routing processes. Additionally, Liang and Li (2019) [69] found the correlations between meteorological driving factors, especially precipitation, and lake area had a lag period of approximately 1-3 months at the monthly time scale. This can be observed in Figure 6 and explained by the generation mechanisms of hydrologic processes related to HJN Lake.

Guidance for Local Water Management
The new approach can enhance our comprehension of the water cycle in semi-aired and ungauged closed watersheds, and guide local water resource management. Though HJN Lake expanded in 2016 and 2017, previous studies (e.g., [30,70]) and our research show that the expansion was due to increased precipitation. Therefore, it is still essential to implement scientific management policies to protect HJN Lake and its associated ecological system. Precipitation and surface runoff are the main sources of water supply for HJN Lake, and the Manggaitu River and Qibosu River are the most significant and stable rivers that contribute water to HJN Lake. Therefore, conserving surface runoff, particularly the streamflow of the Manggaitu River and Qibosu River, may be an efficient way to prevent the shrinkage of HJN Lake.

Limitations and Future Work
The results of this study are subject to the influences of input data and model functions. Since most of the input data were processed from remote-sensing data with resolutions of 250 m and 1 km, using high-resolution data can enhance the accuracy of this study. It has been observed in many studies that the Hargreaves equation tends to underestimate the potential ET in semi-arid regions, indicating that the equation should be adapted for specific locations [62][63][64]. However, due to the lack of observed data, this study utilized the original equation to calculate the potential ET of the HJN watershed. Therefore, obtaining more observed meteorological data can improve the calculated potential ET. Moreover, the groundwater module of the HYDROL-DC model is relatively basic, which may cause the underestimation of groundwater flow during long-term drought conditions. To address this issue, future studies can consider enhancing the HYDROL-DC model by integrating it with groundwater models such as MODFLOW.

Conclusions
In this study, a comprehensive approach was proposed to simulate lake variations in semi-arid and ungauged closed watersheds and analyze potential hydrologic driving factors of the variations. Unlike previous studies that used water-balance models, the new approach integrated lake-oriented and physically based hydrologic modeling, remotesensing data, and results from previous studies to reveal detailed hydrologic processes related to lake variations.
The application in the HJN watershed demonstrated that using remote-sensed ET and lake area data and the results from previous studies effectively solved the data shortage issue in applying physically based hydrologic models to the HJN watershed. During the entire simulation period, rainfall and snowmelt were the most stable supplements to HJN Lake, followed by groundwater inflow, upstream inflow, and CA inflow, due to the differences in their generation mechanisms. In the decrease period (i.e., 2008-2011), rainfall and snowmelt were the major supplements of HJN Lake, while the CA inflow, upstream inflow, and groundwater inflow dominated lake variations in the increase period (2016-2017). Different rivers displayed distinct contributions to HJN Lake, with the Manggaitu River and Qibosu River being the primary and most stable rivers. The generation mechanisms also influenced the correlations between each hydrologic process and the storage variation of HJN Lake at different time scales. The correlation analysis revealed that evaporation dominated the shrinkage of HJN Lake, with the coefficient of correlation at the daily time scale slightly greater than that at the monthly time scale. Rainfall and snowmelt were the key factors for lake storage increase at the daily time scale, while the upstream inflow and CA inflow were the major factors at the monthly time scale.
Overall, the approach proposed in this study effectively simulated the detailed variations and associated hydrologic processes of HJN Lake. The results of this study can not only enhance our comprehension of the hydrological processes related to HJN Lake but also provide valuable guidance for local water resources management and ecological protection. Moreover, this approach can serve as an alternative option for hydrological simulation in other semi-arid and ungauged closed watersheds.