Evolution of Main Water Cycle Fluxes in the Karst Mountain Region of Southwest China

: Distributed hydrological simulation in karst regions has always been a challenging task because of their unique hydrogeological characteristics. The karst mountain region of southwest China (KMRSC), one of the largest continuous karst areas in the world, contributes to about 54 percent of water supply in the basins. In spite of its importance, we have a poor understanding of the evolution laws of hydrological cycle and water resources in KMRSC. We developed a physically-based, distributed hydrological model, called Water and Energy transfer Processes (WEP)-karst model, for KMRSC by introducing the equivalent porous medium approach to the WEP-L model, and dividing the modelling domain into 2021 sub-watersheds. The area of sub-watersheds ranges from 55 to 920 km 2 , with an average value of 170 km 2 . The model showed a good performance in simulating the monthly discharge at 18 representative hydrological stations, with the Nash–Sutcli ﬀ e e ﬃ ciency (NSE) values ranging from 0.71 to 0.94, and the relative error (RE) values from − 9.8% to 8.3% during the validation period (1980–2000). Then, we employed an in-depth analysis of the temporal and spatial variation of main water cycle ﬂuxes, including precipitation, inﬁltration, evapotranspiration, blue water (i.e., river runo ﬀ ), and green water (i.e., vegetation transpiration) over 1956–2015. In addition, the impact of climate change on these ﬂuxes was evaluated under the median emission scenario (RCP4.5). The results showed that: (1) annual average precipitation of KMRSC reached 1506 mm, which is 2.4 times of the national average level, and about 47% (701 mm) of it contributed to river runo ﬀ . The inﬁltration and evapotranspiration were 862 and 870 mm, respectively. The transpiration from plants and trees accounted for 51% of the evapotranspiration. (2) Except for the green water, other ﬂuxes experienced a signiﬁcant decrease over the past 60 years. Blue water showed the largest interannual ﬂuctuation and the strongest sensitivity to climate change. (3) Both precipitation and inﬁltration concentrated from May to August, and blue water increased notably from May to June and peaked in June. Blue water and precipitation were more likely to decrease in the future over 2021–2050 due to the climate change.


Introduction
Water is critical for industrial and agricultural production and ecological sustainability. Due to the increasingly intense human activities and climate change, water scarcity has become a bottleneck hampering sustainable development around the world [1]. Mountain regions are referred to as the world's "water tower" because they provide the lowlands with essential freshwater for human and evapotranspiration were selected for analysis in this paper. Moreover, considering the coordination of water use in production, living, and ecological protection, the blue water and green water were also taken into account in our model. According to the concept proposed by Falkenmark and Rockström [20], freshwater resources can be classified as blue water or green water. The former is traditionally quantified as the runoff that can be directly used for human consumption, while the latter is the summation of the actual evaporation (the nonproductive part) and the actual transpiration (the productive part) [21][22][23]. In some references, only the transpiration is regarded as the green water component [24,25]. In this study, blue water was defined as the river runoff that provides more than 90% of the water used for production and living in KMRSC. Green water was defined as vegetation evaporation (i.e., rainwater used efficiently by the vegetation), which is an important indicator of the ecology of mountain vegetation.

Study Area
The boundary of KMRSC was determined based on the south China karst topographic map as well as the Asian karst topographic distribution map ( Figure 1). The land area is 21.3 × 10 4 km 2 , covering the middle and upper basins of Pearl River and Yangtze River. In terms of administrative divisions, this area involves the Guizhou and Guangxi provinces, as listed in Table 1. The region is characterized by the subtropical wet monsoon climate, with abundant rainfall and high temperature. The annual average values of precipitation and temperature are about 1500 mm and 17.1 • C, respectively. The topography is high in the northwest and low in the southeast. The land-use types in low-, middle-, and high-altitude areas are mainly paddy field and dry land, forest land, and grassland, respectively. The total population was about 48 million by the end of 2015, and its GDP reached 1765 billion Yuan. KMRSC involves multiple closed and unclosed basins, while the hydrological model is generally applied to closed basins. Based on the river system involved in KMRSC, the simulated model domain was appropriately expanded to cover the study area. Considering the integrity of the river system and the demands of integrated river basin water management, China is divided into 10 Class I Water Resources Regions (WRRs), 80 Class II WRRs, and 210 Class III WRRs. The modelling domain in this study covered seven of Class III WRRs ( Figure 2).   Note: P and ET represent the annual average precipitation and evapotranspiration, respectively.

Data Preparation
To quantitatively analyze the evolution of regional water cycle fluxes in this research, three types of data needed to be collected and prepared, i.e., model input data, model validation data and the simulated data of regional climate model.
Firstly, model input data were employed to run the distributed karst simulation model. The input data can be classified into seven categories: (i) hydro-meteorology, (ii) topography, (iii) soil, (iv) river system, (v) hydrogeology, (vi) land use, and (vii) vegetation.
(i) Hydro-meteorology data consisted of five items: rain/snow, air temperature, sunshine hours, vapor pressure/relative humidity, and wind speed with daily temporal resolution during the period 1956-2015. Sixty-nine key national meteorological stations inside and around the modelling domain, which are parts of the National Meteorological Information Center of China, were selected for this purpose (Figure 3a).

Data Preparation
To quantitatively analyze the evolution of regional water cycle fluxes in this research, three types of data needed to be collected and prepared, i.e., model input data, model validation data and the simulated data of regional climate model.
Firstly, model input data were employed to run the distributed karst simulation model. The input data can be classified into seven categories: (i) hydro-meteorology, (ii) topography, (iii) soil, (iv) river system, (v) hydrogeology, (vi) land use, and (vii) vegetation.
(i) Hydro-meteorology data consisted of five items: rain/snow, air temperature, sunshine hours, vapor pressure/relative humidity, and wind speed with daily temporal resolution during the period 1956-2015. Sixty-nine key national meteorological stations inside and around the modelling domain, which are parts of the National Meteorological Information Center of China, were selected for this purpose (Figure 3a).
(ii) Topography data were obtained from the Shuttle Radar Topography Mission (SRTM) with the spatial resolution of 90 m (available at the Geospatial Data Cloud website of Computer Network Information Center of the Chinese Academy of Sciences: http://www.gscloud.cn.) (iii) Soil data including map of the 1km grid soil types ( Figure 3b) and the corresponding characteristic parameters were obtained from the National Second Soil Survey Data.
(iv) The river system included a river network map and river cross-sections. A surveyed river network with a scale of 1:250,000 was obtained from the National Geography Database. The cross-section information was obtained from the hydrological yearbook.
(v) The hydrogeology data of aquifers including hydrogeology unit division, permeability, and storage coefficient were obtained from the Distribution Map of Hydrogeology in China.
(vi) The land-use information was based on the national 1km grid land-use type data provided in seven periods of 1980, 1990, 1995, 2000, 2005, 2010, and 2015. This data set was obtained from the Chinese Academy of Sciences ( Figure 3c).
(vii) The monthly vegetation information (leaf area index and area fraction of vegetation) from 2001 to 2015 was obtained using the Moderate Resolution Imaging Spectroradiometer (MODIS) satellite data ( Figure 3d).
Secondly, for validation of the model, we used the monthly discharge data at the main river cross-sections (the so-called model validation data). Since the model simulates natural hydrological processes without considering the water use, the simulated monthly discharges should be compared with the statistical ones, which are the traditional naturalized river discharges obtained by reverting statistical water consumption to the observed ones. To ensure the reliability of the validation data, the results of the second and most recent national water resources survey completed in 2004 were used. We selected 18 representative hydrological stations that are located throughout the main rivers in KMRSC, and computed their naturalized monthly discharges from 1956 to 2000, as shown in Figure 3a.
Thirdly, the simulated data of the regional climate model RegCM4.0 under the Representative Concentration Pathway (RCP) 4.5 scenarios (i.e., the median emission scenarios) was used to assess the impact of climate change on water cycle fluxes in KMRSC. This data set was provided by the National Climate Center of China Meteorological Administration, including the historical simulation from 1950 to 2005 and the climate change simulation from 2006 to 2099, with a grid spacing of 50 km [26].
(iv) The river system included a river network map and river cross-sections. A surveyed river network with a scale of 1:250,000 was obtained from the National Geography Database. The crosssection information was obtained from the hydrological yearbook.
(v) The hydrogeology data of aquifers including hydrogeology unit division, permeability, and storage coefficient were obtained from the Distribution Map of Hydrogeology in China.
(vi) The land-use information was based on the national 1km grid land-use type data provided in seven periods of 1980, 1990, 1995, 2000, 2005, 2010, and 2015. This data set was obtained from the Chinese Academy of Sciences ( Figure 3c).
(vii) The monthly vegetation information (leaf area index and area fraction of vegetation) from 2001 to 2015 was obtained using the Moderate Resolution Imaging Spectroradiometer (MODIS) satellite data (Figure 3d).
Secondly, for validation of the model, we used the monthly discharge data at the main river cross-sections (the so-called model validation data). Since the model simulates natural hydrological processes without considering the water use, the simulated monthly discharges should be compared with the statistical ones, which are the traditional naturalized river discharges obtained by reverting statistical water consumption to the observed ones. To ensure the reliability of the validation data, the results of the second and most recent national water resources survey completed in 2004 were used. We selected 18 representative hydrological stations that are located throughout the main rivers in KMRSC, and computed their naturalized monthly discharges from 1956 to 2000, as shown in Figure 3a.
Thirdly, the simulated data of the regional climate model RegCM4.0 under the Representative Concentration Pathway (RCP) 4.5 scenarios (i.e., the median emission scenarios) was used to assess the impact of climate change on water cycle fluxes in KMRSC. This data set was provided by the National Climate Center of China Meteorological Administration, including the historical simulation from 1950 to 2005 and the climate change simulation from 2006 to 2099, with a grid spacing of 50 km [26].

Description of WEP-L Model
The Water and Energy transfer Processes model for Large river basins (WEP-L model) is a physically-based distributed hydrological model developed by Jia et al. for dynamic assessment of (c) (d)

Description of WEP-L Model
The Water and Energy transfer Processes model for Large river basins (WEP-L model) is a physically-based distributed hydrological model developed by Jia et al. for dynamic assessment of water resources [27]. The WEP-L model is a physically based distributed hydrological model, and the schematic illustration of its distribution structure is shown in Figure 4. Instead of grid cells, contour belts inside small sub-watersheds are used as computation units. In mountain region, each sub-watershed is further divided into several contour belts based on the difference of maximum and minimum elevation as well as considering the equalization of contour belt area. Because of little topographic change, no further division is performed for sub-watersheds in plain areas. The spatial discretization method has some advantages in reflecting the vertical zonality heterogeneity of meteorology, hydrology and vegetation variables, including their impact on mountain water movement. There are nine vertical layers within a contour belt, namely an interception layer, a depression layer, three upper soil layers, a transition layer, an unconfined aquifer, an aquitard and a confined aquifer. In order to reflect the heterogeneity of the land cover within contour belt, the land use is further divided into five groups by adopting the mosaic method [28]. It consists of water body, soil-vegetation, impervious area, irrigated farmland, and non-irrigated farmland. The soil-vegetation group is further classified into bare soil, tall vegetation (forest or urban trees) and short vegetation (grass or crops). The impervious area group consists of impervious urban cover and urban canopy.

Transition Layer
Groundwater Lift

Division and Codification of Computation Units
Using the Digital Elevation Model (DEM), resampled into 1km spatial resolution, the virtual river networks of the catchment were extracted, and the sub-watersheds were divided. The entire modelling catchment was divided into 2021 sub-watersheds, 1205 of which were in KMRSC, as shown in Figure 5. The area of sub-watersheds ranges from 55 to 920 km 2 , with an average value of   The coupling simulations of natural hydrological processes and energy transfer processes are taken into account. Evapotranspiration is calculated using the Penman-Monteith equation. During periods of heavy rainfall (>10 mm/h), the multi-layered Green-Ampt model [29] is used to calculate the infiltration excess runoff generation, and the saturation excess runoff is simulated through a water balance analysis of the unsaturated soil layers. The simulated energy transfer processes include short-wave radiation, long-wave radiation, latent heat flux, sensible heat flux and soil heat flux. We classify the land cover within the computation unit using the land use type data during the modelling period, and then estimate the water and heat fluxes of each land cover. The areal average of water and heat fluxes from all land covers in a contour belt produces the averaged fluxes in the computation unit. This feature improves the determination of the actual evapotranspiration. A one-dimensional kinematic wave approach is adopted for overland and riverine runoff routings. Additionally, the interactions between surface water and groundwater are considered using the Boussinesq equations to simulate groundwater flow. Snowmelt is simulated by employing a temperature index method.
For more details of this model, refer to Jia et al. [27,30].

Division and Codification of Computation Units
Using the Digital Elevation Model (DEM), resampled into 1km spatial resolution, the virtual river networks of the catchment were extracted, and the sub-watersheds were divided. The entire modelling catchment was divided into 2021 sub-watersheds, 1205 of which were in KMRSC, as shown in Figure 5. The area of sub-watersheds ranges from 55 to 920 km 2 , with an average value of 170 km 2 . The coincidence degree between simulated and actual boundary of the domain reached 99%. Furthermore, each sub-watershed was divided into 1-10 contour belts. Thus, 12,792 contour belts inside small sub-watersheds were obtained as computation units of the model with the average area of 26.9 km 2 . The codification of computation units was carried out based on stem-branch topological relationship of river networks [31].

Improvement of Water Movement Simulation in Karst Vadose Zone
In traditional distributed hydrological models such as WEP-L, the vadose zone is simply characterized as a soil-bedrock structure, where soil water movement parameters are determined solely based on the texture information of different soil types, such as sand soil, loam soil, and clay soil. However, these parameters fail to properly represent the soil water dynamics in the vadose zone in the karst region. Specifically, the karst vadose zone has a particular structure that includes soil, epikarst zone and unsaturated zone [8,11]. The uppermost soil layer is usually thin and even missing. Epikarst zone plays a key role in water storage and movement, and transpiration of vegetation root systems, where the fissure network is strongly developed. Moreover, as the extent and frequency of widening diminishes gradually with depth, the epikarst permeability diminishes with depth [32]. The unsaturated zone has a poor fissure development and its water capacity and permeability are significantly reduced compared to that in the epikarst zone. Therefore, the integrated soil-epikarst system mainly controls the infiltration and redistribution of precipitation [33][34][35].

Improvement of Water Movement Simulation in Karst Vadose Zone
In traditional distributed hydrological models such as WEP-L, the vadose zone is simply characterized as a soil-bedrock structure, where soil water movement parameters are determined solely based on the texture information of different soil types, such as sand soil, loam soil, and clay soil. However, these parameters fail to properly represent the soil water dynamics in the vadose zone in the karst region. Specifically, the karst vadose zone has a particular structure that includes soil, epikarst zone and unsaturated zone [8,11]. The uppermost soil layer is usually thin and even missing. Epikarst zone plays a key role in water storage and movement, and transpiration of vegetation root systems, where the fissure network is strongly developed. Moreover, as the extent and frequency of widening diminishes gradually with depth, the epikarst permeability diminishes with depth [32]. The unsaturated zone has a poor fissure development and its water capacity and permeability are significantly reduced compared to that in the epikarst zone. Therefore, the integrated soil-epikarst system mainly controls the infiltration and redistribution of precipitation [33][34][35].
As shown in previous studies, with the increase of basin area, the local influence caused by karst fissure, pipeline and sinkhole diminishes, and an equivalent porous medium approach (EPM) has proved adequate for simulating the karst system in a large basin [36][37][38]. Therefore, this paper attempted to introduce the equivalent porous medium approach into the existing WEP-L model. In this way, the basic laws obtained from small-scale experimental observations in KMRSC were effectively generalized, and the WEP-karst model was developed.
Firstly, the vadose zone in KMRSC was divided into four layers: soil, upper epikarst, lower epikarst and transition layer ( Figure 6). The water capacity and permeability of upper epikarst was lower than that of lower epikarst, and the transition was the lowest. Secondly, based on the field investigations of karst experimental basin, the equivalent soil moisture movement parameters of the four layers were approximately estimated. Investigations by Chen et al. [17] revealed that in KMRSC the medium-permeability fractures were widely distributed, and the hydraulic conductivities in the near-surface epikarst layer were of the order 10 −3 m/s. This is much larger than 10 −5 m/s of the low-permeability zones. Thickness of the epikarst zone varied greatly in space, roughly between 2 and 10 m. Finally, the evaporation and transpiration processes were further modified to depend on vegetation species, and soil and rock fractures. According to Xiang et al. [39], most vegetation roots concentrated in shallow soil and underlying rock fractures in a depth less than 1 m. Soil water and epikarst water are main water sources for plant uptake in the karst forest area. Consequently, evaporation of bare soil was assumed to be caused by the soil layer. Considering root depth, water absorption and transpiration of grass and crop only occurred in the soil and upper epikarst layer, while that of trees covered the soil and the whole epikarst layer. There was no evapotranspiration in the transition layer because it is difficult for tree roots to extend into this layer due to the hard rock. Table 2 listed the initial values of main soil water movement parameters, which were roughly estimated based on the field investigations of karst experimental basin and literature research. Moreover, vertical distribution of saturated hydraulic conductivity in the epikarst and transition layers was represented by an exponentially decreasing function, namely: where, k 0 is the saturated hydraulic conductivity at the surface of epikarst layer; l is the (positive) distance taken in the downward direction normal to the surface of epikarst layer; and a is the parameter.
The default values of k 0 and a are set to 0.5 cm/s and 0.25, respectively.
Considering the distribution of hydrological stations and rivers, KMRSC was divided into nine parameter tuning zones, as shown in Figure 7. The upper and lower limits of the parameters tuned were two times and 0.5 times of their initial values respectively, and 0.01 times of the initial values Water 2020, 12, 2262 9 of 23 was selected as the tuning step. Finally, these parameters were tuned according to runoff process simulation of representative hydrological stations.
where, k0 is the saturated hydraulic conductivity at the surface of epikarst layer; l is the (positive) distance taken in the downward direction normal to the surface of epikarst layer; and a is the parameter. The default values of k0 and a are set to 0.5 cm/s and 0.25, respectively. Considering the distribution of hydrological stations and rivers, KMRSC was divided into nine parameter tuning zones, as shown in Figure 7. The upper and lower limits of the parameters tuned were two times and 0.5 times of their initial values respectively, and 0.01 times of the initial values was selected as the tuning step. Finally, these parameters were tuned according to runoff process simulation of representative hydrological stations.

Criteria for Model Calibration and Validation
Using the WEP-karst model, continuous simulations over 60 years from 1956 to 2015 were conducted in 2021 sub-watersheds and 12792 contour belts around the modelling catchment for natural hydrological processes. All of the parameters were initially estimated according to land cover information, observation data and remote sensing data. Some other parameters were selected for model calibration by sensitivity analysis. The calibration was performed based on "trial and error". The calibration parameters were chosen according to the sensitivity analysis by Jia et al. [27] and included maximum depression storage depth of land surface, soil saturated hydraulic conductivity,

Criteria for Model Calibration and Validation
Using the WEP-karst model, continuous simulations over 60 years from 1956 to 2015 were conducted in 2021 sub-watersheds and 12792 contour belts around the modelling catchment for natural hydrological processes. All of the parameters were initially estimated according to land cover information, observation data and remote sensing data. Some other parameters were selected for model calibration by sensitivity analysis. The calibration was performed based on "trial and error". The calibration parameters were chosen according to the sensitivity analysis by Jia et al. [27] and included maximum depression storage depth of land surface, soil saturated hydraulic conductivity, hydraulic conductivity of unconfined aquifer, permeability of riverbed material, Manning roughness, snow melting coefficient, and critical air temperature for snow melting. The calibration and validation periods were chosen as 1965-1980 and 1981-2000, respectively. Moreover, two frequently-used criteria for model calibration and validation were employed: (1) minimizing the relative error (RE) of annually averaged river runoff, and (2) maximizing the Nash-Sutcliffe efficiency (NSE) of monthly discharge.
where Q sim,i and Q obs,i are the simulated and observed values of monthly runoff, m 3 /s; N is the number of months of simulation period; and Q obs,i is the annual average of the measured monthly runoff, m 3 /s.

Improved Performance of the WEP-Karst Model
The overall improvement effect of soil moisture movement was obvious in KMRSC. The runoff processes simulated by the WEP-karst model were more consistent with statistical results with the higher NSE and lower RE. Through the model improvement, the NSE values of the validation period in KMRSC were increased by 13% and the average of the absolute value of RE was reduced to less than 6%. Furthermore, the downstream control station of Longjiang River basin, the Sancha station, was selected as a representative station to analyze the improvement effect of water movement simulation in the karst vadose zone. Longjiang River basin (LJ), located in the southeast of the study area (Figure 3a), is a typical karst basin in Southwest China [40]. By comparing the simulation results of the WEP-L model and WEP-karst model, it was found that whether the water capacity and permeability of epikarst zone were considered or not has a great influence on the simulation of hydrological processes in the karst area, as shown in Table 3 and Figure 8. When the traditional WEP-L model that does not consider epikarst zone was used, the water storage capacity of the watershed was underestimated, and runoff yield was overestimated. Moreover, the peak discharge was overestimated and river base flow was underestimated. This resulted in a steeper monthly discharge hydrograph. By accounting for epikarst zone, the monthly discharge hydrograph simulated by the WEP-karst model was in a better agreement with the observations from the Sancha station, as demonstrated by the higher NSE and lower RE values.

Model Applicability in the Karst Area
Calibration and validation results of the simulated monthly discharge obtained by the WEP-karst model at 18 representative hydrological stations were reported in Table 4 and shown in Figure 9. For the calibration period, the NSE values were in the range 0.73-0.94, and the RE values ranged from −6.4% to 8.5%. There were 14 stations (i.e., 80% of the total number of stations) with an NSE value larger than 0.8 and an absolute value of RE less than 5%. Only the Jiangbianjie and Zhexiang stations have NSE values smaller than 0.8, and the Sinan, Jiangbianjie and Baise stations have the absolute value of RE larger than 5%. Similarly, for most stations, the NSE values were mostly above 0.8, and the RE values were controlled between −5% and 5% during the validation period. In general, both the NSE and RE values were quite encouraging, which showed that WEP-karst model has a good applicability across KMRSC.  Table 4 and shown in Figure  9. For the calibration period, the NSE values were in the range 0.73-0.94, and the RE values ranged from −6.4% to 8.5%. There were 14 stations (i.e., 80% of the total number of stations) with an NSE value larger than 0.8 and an absolute value of RE less than 5%. Only the Jiangbianjie and Zhexiang stations have NSE values smaller than 0.8, and the Sinan, Jiangbianjie and Baise stations have the absolute value of RE larger than 5%. Similarly, for most stations, the NSE values were mostly above 0.8, and the RE values were controlled between −5% and 5% during the validation period. In general, both the NSE and RE values were quite encouraging, which showed that WEP-karst model has a good applicability across KMRSC.

Vertical Water Cycle Fluxes
The annual average precipitation, infiltration, and evapotranspiration from 1956 to 2015 were determined in all 2021 sub-watersheds using the WEP-karst model. Furthermore, the amounts were aggregated to the seven Class III WRRs and KMRSC, as shown in Figure 10 and Table 5. The annual average precipitation in the modelling domain ranged from 900 to 2240 mm. It gradually increased from northwest to southeast, affected by the southwest and southeast monsoon. Moreover, the annual average precipitation of less than 1300 mm formed a belt distribution in the northwest. There existed three high precipitation centers larger than 1800 mm, concentrated in the east and south. For the seven Class III WRRs, LJR had the largest precipitation, followed by the ZYM and HSR, with the annual average of nearly 1600 mm and higher. In USS and SPR, the values were the lowest (~1200 mm). In the case of KMRSC, the annual average precipitation over 1956-2015 reached 1506 mm, which is about 2.4 times the national average. Affected by meteorological, hydrological, vegetation, and soil conditions, the annual average evapotranspiration generally showed a spatial distribution pattern with higher values in the south and low values in the north. The evapotranspiration exhibited large spatial variability in the modelling domain, varying from 550 to 1300 mm. The lowest value was in USS located in the Yangtze River basin, only 732 mm. However, in HSR and YJR, it reached about 1000 mm. There were some areas with high evaporation in SPR, LJR, HSR, and YJR. The annual average evapotranspiration in KMRSC was 870 mm, which is roughly equivalent to the infiltration.

Blue Water and Green Water
The annual average blue water over 1956-2015 showed a similar spatial pattern as precipitation, which had a range of 100-1500 mm (Figure 11a). The areas with abundant blue water resources were concentrated in the east of LJR, the middle of HSR and the northwest of ZYM. The river runoff yield in the upper reaches of SPR was relatively small (less than 400 mm). At the Class III WRRs scale, the amount of blue water was the largest in LJR (918 mm), which was significantly higher than that in other regions. In contrast, USS, YJR, and SPR had low blue water resources, which varied from 417 to 560 mm. In KMRSC, the annual average blue water reached 701 mm, which means that the local per capita river runoff was about 3100 m 3 .
Green water resource is an important part of regional evapotranspiration, especially in areas (c) In the whole modelling domain, spatial distribution of the annual average infiltration was mainly reflected in the difference between the upstream and downstream, with a range of 293-1934 mm. In most basins, the amount of infiltration in the upper reaches was generally smaller than that in the lower reaches. The reason may be that the thickness of the vadose zone in the upper reaches is relatively small, resulting in a smaller soil water storage capacity and infiltration. Moreover, the closer the area to the main channel, the greater the infiltration values. There was no significant difference in the annual average infiltration among these Class III WRRs. The infiltration in USS was the smallest, 749 mm, while that in HSR was the largest, 954 mm. KMRSC had an annual average infiltration of 862 mm, which accounted for 57% of the precipitation.
Affected by meteorological, hydrological, vegetation, and soil conditions, the annual average evapotranspiration generally showed a spatial distribution pattern with higher values in the south and low values in the north. The evapotranspiration exhibited large spatial variability in the modelling domain, varying from 550 to 1300 mm. The lowest value was in USS located in the Yangtze River basin, only 732 mm. However, in HSR and YJR, it reached about 1000 mm. There were some areas with high evaporation in SPR, LJR, HSR, and YJR. The annual average evapotranspiration in KMRSC was 870 mm, which is roughly equivalent to the infiltration.

Blue Water and Green Water
The annual average blue water over 1956-2015 showed a similar spatial pattern as precipitation, which had a range of 100-1500 mm (Figure 11a). The areas with abundant blue water resources were concentrated in the east of LJR, the middle of HSR and the northwest of ZYM. The river runoff yield in the upper reaches of SPR was relatively small (less than 400 mm). At the Class III WRRs scale, the amount of blue water was the largest in LJR (918 mm), which was significantly higher than that in other regions. In contrast, USS, YJR, and SPR had low blue water resources, which varied from 417 to 560 mm. In KMRSC, the annual average blue water reached 701 mm, which means that the local per capita river runoff was about 3100 m 3 . not a single closed basin, but consists of several unclosed basins. Consequently, there existed a certain amount of lateral groundwater recharge from outside KMRSC during the simulation period , which is about 1.38 million m 3 . Note: Amount and Proportion represent the annual average green water and its proportion in evapotranspiration, respectively.

Annual Variability
Focused on KMRSC, the annual change trends of water cycle fluxes from 1956 to 2015 were analyzed through the nonparametric Mann-Kendall test [41]. Besides, the coefficient of variation (CV) Green water resource is an important part of regional evapotranspiration, especially in areas with high precipitation and dense vegetation cover. In the modelling domain, the annual average of green water over 1956-2015 was from 36 to 886 mm. Figure 11b showed that spatial variation of the green water seems to be consistent with that of the evapotranspiration and Normalized Difference Vegetation Index (NDVI). Affected by the precipitation and vegetation conditions, there existed some differences in the annual average green water and its proportion in the evapotranspiration among the seven Class III WRRs, as listed in Table 6. The annual average in YJR was the largest with 562 mm, while in NPR it was only 346 mm. Moreover, the green water resources in USS and NPR were similar, but their proportions were significantly different. As seen in Figure 10a, it may be mainly due to the precipitation difference in the two regions.
By comparing the spatial distribution of blue and green water, some interesting facts were revealed. Firstly, in YJR and SPR, the green water resource was relatively high, while the amount of blue water was low. This may be related to the vegetation capacity to hold soil moisture and reduce runoff. Especially in the upper reaches of the YJR, most of the precipitation was consumed by vegetation transpiration, so the reasonable vegetation coverage threshold needs to be further explored in this regard. Secondly, in the southwest of KMRSC, we observed low green water resource and high river runoff. Due to the poor vegetation conditions, this area faces relatively high flood risk. Thirdly, in some areas such as the east of LJR and the lower reaches of HSR, the precipitation was large and the vegetation conditions were good. Thus, both the blue and green water were relatively high. In this case, one should focus on the balance between water uses in industrial and agricultural production and ecology.
It can be noted that the sum of the annual average values of evapotranspiration and runoff depth in KMRSC is 1571 mm, which is larger than the annual average precipitation (1506 mm). KMRSC is not a single closed basin, but consists of several unclosed basins. Consequently, there existed a certain amount of lateral groundwater recharge from outside KMRSC during the simulation period , which is about 1.38 million m 3 . Note: Amount and Proportion represent the annual average green water and its proportion in evapotranspiration, respectively.

Annual Variability
Focused on KMRSC, the annual change trends of water cycle fluxes from 1956 to 2015 were analyzed through the nonparametric Mann-Kendall test [41]. Besides, the coefficient of variation (CV) was adopted to analyze the annual fluctuation of these fluxes. It can be seen from Table 7 and Figure 12 that all fluxes showed a decreasing trend, of which the annual precipitation, infiltration, and evapotranspiration decreased at the p = 0.01 significant level with a change rate of −2.59, −1.48, and −1.65 mm/year, respectively. The amount of blue water resource decreased at the p = 0.05 significant level with a change rate of −2.07 mm/year, while green water indicated an insignificant trend. It should be noted that due to extreme drought events in southwest China around 2010, precipitation and river runoff decreased to historically low levels [42]. Furthermore, there existed significant differences in the annual fluctuations of the five fluxes. In the case of annual blue water, the CV value reached 0.28 and the fluctuation range was large. This brings more uncertainty to regional water supply. The annual fluctuation of the precipitation was close to the average level in southern China with the CV value of 0.13. In contrast, the annual infiltration showed a similar fluctuation range with the precipitation, while the evapotranspiration and green water had a weak fluctuation. This may be due to the rich and stable groundwater in the karst region, resulting in the limited impact of annual fluctuations of precipitation and runoff on the evapotranspiration and green water. reached 0.28 and the fluctuation range was large. This brings more uncertainty to regional water supply. The annual fluctuation of the precipitation was close to the average level in southern China with the CV value of 0.13. In contrast, the annual infiltration showed a similar fluctuation range with the precipitation, while the evapotranspiration and green water had a weak fluctuation. This may be due to the rich and stable groundwater in the karst region, resulting in the limited impact of annual fluctuations of precipitation and runoff on the evapotranspiration and green water.    Table 8 lists the monthly average of water cycle fluxes in KMRSC from 1956 to 2015. Moreover, the box diagram was chosen presentation of the intra-annual variability of water cycle fluxes in KMRSC, as shown in Figure 13. It reflected not only the monthly distribution of these fluxes, but also the change of monthly data over the past 60 years. A sudden increase of precipitation was recorded in April and May. The highest and lowest monthly precipitation was 273 mm in June and 36 mm in January, respectively. The precipitation is concentrated in May to August, accounting for 62.2% of the annual total precipitation. During the rainy period from June to August, the precipitation showed a relatively strong fluctuation. In dry months, the precipitation fluctuation happened to the largest in October. Similar to the case of the precipitation, 58.6% of annual total infiltration was distributed from May to August. However, in contrast, the fluctuation range of infiltration was relatively large in each month. Compared to the precipitation and infiltration, monthly evapotranspiration data seemed to be more evenly distributed and showed smaller fluctuation. The maximum monthly evaporation appeared in July with 129 mm, which happens to be one month later than that of precipitation and infiltration. Blue water increased strongly from May to June and peaked in June with 136 mm. In the first three months of the year, the blue water was the smallest, accounting for only 6.6% of the annual blue water resource. Interestingly, the larger monthly blue water, the stronger fluctuation it seems. Because of the close relationship between green water and vegetation growth, the distribution difference of green water in each month was obviously expected. The green water in July and August accounted for 38.9% of the total annual value, which was significantly higher than that in other months. The fluctuation range of green water between May and September was larger than that of other months.

Assessing the Impact of Climate Change
Considering the rough prediction data, the impact of climate change on water cycle fluxes was assessed and discussed at the Class III WRRs scale using the WEP-karst model, as shown in Table 9. In this analysis, the period 1956-2015 was chosen as the base period, and the period 2021-2050 was chosen as the prediction period. As shown in Figure 14, a clear trend with the decrease in precipitation and blue water resource under the median emission scenario (RCP4.5) was observed. The infiltration, evapotranspiration and green water changed slightly compared to the base period results. Among the Class III WRRs, the change rate of annual average precipitation varied between −4.7% and 0.5%. Regarding the spatial variability, the river runoff yield in USS, SPR, and NPR, which are located in the upper reaches of KMRSC, seemed to be more sensitive to climate change and had a more dramatic decrease with a change rate of about −16%. As a result, increasing water scarcity risk is expected in KMRSC. The change in annual average infiltration, evapotranspiration and green water in HSR was larger than that in other regions. It may be that HSR is in the downstream area, with good vegetation conditions and thick soil layer. Therefore, it may be more influenced by precipitation reduction and temperature increase.

Assessing the Impact of Climate Change
Considering the rough prediction data, the impact of climate change on water cycle fluxes was assessed and discussed at the Class III WRRs scale using the WEP-karst model, as shown in Table 9. In this analysis, the period 1956-2015 was chosen as the base period, and the period 2021-2050 was chosen as the prediction period. As shown in Figure 14, a clear trend with the decrease in precipitation and blue water resource under the median emission scenario (RCP4.5) was observed. The infiltration, evapotranspiration and green water changed slightly compared to the base period results. Among the Class III WRRs, the change rate of annual average precipitation varied between −4.7% and 0.5%. Regarding the spatial variability, the river runoff yield in USS, SPR, and NPR, which are located in the upper reaches of KMRSC, seemed to be more sensitive to climate change and had a more dramatic decrease with a change rate of about −16%. As a result, increasing water scarcity risk is expected in KMRSC. The change in annual average infiltration, evapotranspiration and green water in HSR was larger than that in other regions. It may be that HSR is in the downstream area, with good vegetation conditions and thick soil layer. Therefore, it may be more influenced by precipitation reduction and temperature increase.

Water Scarcity Challenges in KMRSC
The results on the evolution characteristics of water cycle fluxes contribute to better utilization and management of water resources in KMRSC. The Falkenmark indicator has been widely used for measuring water scarcity. This indicator is based on the annual river runoff available for human use (i.e., the blue water resource) and makes use of the 1700 and 1000 m 3 per capita thresholds to distinguish between the water stressed and scarce areas [43]. In KMRSC, the blue water provides almost all the water resource used for the local production of industry, agriculture, and people's daily life. However, due to the significant spatial variability in water resource conditions and human activity intensity, several regions in KMRSC are faced with different levels of water stress. We used the statistical population data at the level of county, and interpolated them into 1km grid cell using the GIS technique based on land-use distribution. Then the local per capita blue water in 2021 sub-watersheds was calculated and used to identify water stressed areas ( Figure 15).

Water Scarcity Challenges in KMRSC
The results on the evolution characteristics of water cycle fluxes contribute to better utilization and management of water resources in KMRSC. The Falkenmark indicator has been widely used for measuring water scarcity. This indicator is based on the annual river runoff available for human use (i.e., the blue water resource) and makes use of the 1700 and 1000 m 3 per capita thresholds to distinguish between the water stressed and scarce areas [43]. In KMRSC, the blue water provides almost all the water resource used for the local production of industry, agriculture, and people's daily life. However, due to the significant spatial variability in water resource conditions and human activity intensity, several regions in KMRSC are faced with different levels of water stress. We used the statistical population data at the level of county, and interpolated them into 1km grid cell using the GIS technique based on land-use distribution. Then the local per capita blue water in 2021 subwatersheds was calculated and used to identify water stressed areas ( Figure 15).
There existed many regions with small local per capita blue water resource of less than 1700 m³. These regions were widely distributed in the west, north, and southeast of the modelling domain. In some areas, water scarcity was very prominent, because the per capita blue water resource was less than 1000 or even less than 500 m 3 . USS, NPR, and SPR may be subjected to a more severe water stress, because these areas are the water source areas and lack of upstream water supply. Meanwhile, the ecological environment in the upstream areas is generally fragile. Once the environment has been deteriorated, it is difficult to recover the weakened ecosystem. As for KMRSC, shaped like a dumbbell, the water stressed areas were located mainly at the two ends of the dumbbell. KMRSC had 19,437 km 2 of water stressed area and 44,625 km 2 of water scarce area. These values accounted for 9.7% and 22.2% of the total area, respectively. Considering the strong possibility of precipitation reduction and decrease in water resources due to the impact of climate change, the utilization of water resources may face greater challenges. It is suggested that the water stressed areas in the northwest KMRSC need to give priority to the protection of water source and reasonable control of local population growth. The economic and social development is mainly concentrated in the southeast, where the regional water resource conditions are poor. Thus, it is necessary to pay attention to the improvement of water use efficiency in industrial and agricultural sectors, and invest on the construction of water conservancy infrastructure to improve the water supply capacity. Besides, efficient use of green water resources could also be an alternative solution facing the increasingly serious water resources scarcity.  There existed many regions with small local per capita blue water resource of less than 1700 m 3 . These regions were widely distributed in the west, north, and southeast of the modelling domain. In some areas, water scarcity was very prominent, because the per capita blue water resource was less than 1000 or even less than 500 m 3 . USS, NPR, and SPR may be subjected to a more severe water stress, because these areas are the water source areas and lack of upstream water supply. Meanwhile, the ecological environment in the upstream areas is generally fragile. Once the environment has been Water 2020, 12, 2262 20 of 23 deteriorated, it is difficult to recover the weakened ecosystem. As for KMRSC, shaped like a dumbbell, the water stressed areas were located mainly at the two ends of the dumbbell. KMRSC had 19,437 km 2 of water stressed area and 44,625 km 2 of water scarce area. These values accounted for 9.7% and 22.2% of the total area, respectively. Considering the strong possibility of precipitation reduction and decrease in water resources due to the impact of climate change, the utilization of water resources may face greater challenges. It is suggested that the water stressed areas in the northwest KMRSC need to give priority to the protection of water source and reasonable control of local population growth. The economic and social development is mainly concentrated in the southeast, where the regional water resource conditions are poor. Thus, it is necessary to pay attention to the improvement of water use efficiency in industrial and agricultural sectors, and invest on the construction of water conservancy infrastructure to improve the water supply capacity. Besides, efficient use of green water resources could also be an alternative solution facing the increasingly serious water resources scarcity.

Research Progresses and Prospects
Several previous studies have provided related conclusions on the evolution characteristics of water resources in the karst regions of southwest China. For instance, Wu et al. estimated that Guizhou and Yunnan have been in a warmer and drier conditions in the period from 1956 to 2000 [44]. Wang et al. also observed increasing drought occurred in Southwest China during 1962-2012 by SPEI_pm [45]. Liu et al. concluded that the precipitation and potential evapotranspiration in this region exhibit decreasing trends [9]. Overall, the results of this study agree well with those of previous studies. In this paper, the boundary of KMRSC was further clarified, and the spatiotemporal variation characteristics of main water cycle fluxes were systematically revealed and analyzed by the WEP-karst model.
The WEP-karst model owns a distribution structure composed of sub-watershed, contour belt and land cover mosaic, which has several advantages in reflecting the heterogeneity of meteorology, hydrology and underlying surface conditions. Moreover, the spatial discretization method of the WEP-karst model is similar but different from that of the Soil and Water Assessment Tool (SWAT) Model. SWAT is a public domain model developed by a group of scientists from the USDA-Agricultural Research Service, which has been widely used in hydrological research. In the model, a sub-watershed is further divided into several hydrologic response units (HRUs) that consist of homogeneous land use, management, and soil characteristics. However, the model ignores the hydraulic connection among the HRUs, and the flows are summed directly from all HRUs to the sub-watershed level. Just as Gassman et al. [46] pointed out, "the nonspatial aspect of the HRUs is a key weakness of the model". While in the WEP-karst model, there exists hydraulic connection among the contour belts inside small sub-watersheds. Overland flow is routed from the upstream contour belt to downstream ones in each sub-watershed using the kinematic wave method.
Because the hydrological phenomenon in the karst region is complicated, revealing it precisely is a process that requires continuous improvement. By introducing the equivalent porous medium approach, the hydrogeological complexity and hydraulic heterogeneity caused by micropores, small fissures, large fractures and conduits have been effectively diminished in the WEP-karst model. However, the inconsistency between the boundaries of surface watershed and aquifer failed to be considered in the current study, which may affect the simulation effect of several local areas, especially small-scale watersheds. It is necessary to strengthen the description of the inconsistency between the boundaries of surface watershed and aquifer, which is an important direction to improve the simulation accuracy of the WEP-karst model. Furthermore, the current study focuses on the simulation of the natural hydrological processes for various land cover conditions, but fails to consider the effects of reservoir, water withdrawal, irrigation, and wastewater reuse. The next research plans to couple social water cycle system into the WEP-karst model to improve model applicability in the areas affected by strong human activities.

Conclusions
In this study, the WEP-karst model, a distributed hydrological model for the karst mountain region of Southwest China, was developed. Dynamic assessment of regional water cycle fluxes was conducted. The main conclusions were summarized in the following.
By introducing the equivalent porous medium approach, the regulation and storage effects of epikarst zone on precipitation infiltration was considered, and then the simulation accuracy of hydrological processes in KMRSC was significantly improved. At the 18 representative hydrological stations located throughout the modelling domain, the WEP-karst model showed a good performance in simulating the monthly discharge. A quantitative assessment during the validation period of 1981-2000 showed that the range of NSE and RE values was 0.71-0.94 and −9.8%-8.3%, respectively.
The annual average precipitation over 1956-2015 in KMRSC was 1506 mm, which reached about 2.4 times the national average. It gradually increased from northwest to southeast, with a range of 900 to 2240 mm. Moreover, the low precipitation area less than 1300 mm formed a belt distribution in the northwest of KMRSC. The annual average infiltration and evapotranspiration were 862 and 870 mm, both accounting for about 57% of the precipitation. The annual average blue water reached 701 mm. This means that the local per capita river runoff was about 3100 m 3 . The green water was only 52% of blue water, and showed similar spatial variability with the evapotranspiration and NDVI.
Over the past 60 years, except for green water which experienced a slight decrease, other four fluxes showed a statistically significant decrease trend at 95% confidence level. Blue water was active, which showed notable sensitivity to climate change. With respect to the intra-annual variability, both the precipitation and infiltration concentrated from May to August, while evapotranspiration seemed to be more evenly distributed over the year. Blue water increased strongly from May to June and peaked in June, while green water was mainly distributed in July and August. More importantly, the hydrological predictions over 2021-2050 showed that climate change under the median emission scenario of RCP4.5 is very likely to cause the reduction of precipitation and water resources.
The research can provide a reference for water security in KMRSC and is valuable for water resources regulation and control in other similar mountainous regions around the world. Future studies will concentrate on reflecting the influences of water utilization and flow regulation on hydrological processes into WEP-karst to improve its applicability.