Stochastic Parsimonious Hydrologic Partitioning Model under East Asia Monsoon Climate and Its Application to Climate Change

: A conceptual hydrologic partitioning model suitable for the East Asia monsoon climate region is constructed parsimoniously, and the variability of Horton index, which is the ratio of water vaporization and wetting in the watershed, is investigated. Numerical simulations in the study area show that the inter-annual variability of Horton index is reduced to around 60% of the inter-annual variability of annual precipitation, and there is a strong inverse correlation between Horton index and annual precipitation. Using cumulant expansion theory, the probability distribution function of soil water with various hydro-meteorological variables and watershed characteristics is derived. Using the steady-state soil water probability distribution function, the sensitivity of Horton index to hydro-meteorological variables such as precipitation occurrence probability, average rainfall depth at rainy days, and evapotranspiration rate and hydro-geophysical characteristics such as surface runo ﬀ coe ﬃ cients, threshold soil water value to control vaporization, and exponent value to control groundwater recharge is analyzed. Looking at the future Horton index of the study area using a variety of future climate information ensemble, it is projected that the water stress of vegetation in the watershed is likely to increase due to ﬂuctuations in precipitation patterns and increase in potential evapotranspiration even if annual precipitation increases.


Introduction
Climate change affects the hydrological response of the surface through various paths. Since changes in the hydrological cycle, particularly in the hydrologic partitioning process, result in changes in the ecological environment, climate change also affects the structure and function of ecosystems [1]. Therefore, understanding and predicting the impact of climate change on the watershed's water balance has emerged as an important issue in managing and conserving the natural environment [2]. In particular, soil water near the surface plays an important role in hydrologic partitioning, including precipitation, surface flow, wetting, vaporization, and groundwater recharge processes [3,4]. It is known that the ability of the watershed to store soil water supplied to the soil layer by precipitation and to vaporize it back into the atmosphere is mainly controlled by vegetation [5,6]. Therefore, in order to explore the water balance in the watershed, a variety of studies have been conducted from a simple model to a complex model simulating soil water dynamics [7][8][9][10][11][12][13]. As part of the watershed response to land-use change resulting from human activity, changes in the hydrologic partitioning process have been widely discussed [14]. However, based on these discussions, the impact of climate change on the hydrologic partitioning and even vegetation of the watershed has not been fully investigated yet [15]. Since changes in vegetation provide a very important driver in the climate-ecosystem-hydrologic feedback system through the exchange of carbon, water, and energy [16], studies in the field of eco-hydrology to explore the linkage between hydrological variables and ecosystem responses [17] or the impact of ecosystem structures on hydro-meteorological characteristics of regions [18,19] is actively proceeding. However, many studies have been conducted mainly in arid or semi-arid areas, and relatively few studies in the field of eco-hydrology have been conducted mainly in the East Asia monsoon climate (EAMC), including the Korean peninsula [19,20].
One of the most common approaches to effectively understanding a system's response to external forcings is to model the system and then see what happens when various external forcings are given to the model. In the field of hydrology, deterministic models are suitable for this purpose since they can account for the dynamic interactions of various hydrologic components in response to external forcings [21]. However, the deterministic approach can only achieve its intended purpose if the inputs (both model parameters and external forcings) are provided with sufficient accuracy and resolution to investigate the system's response. However, future climate forcings (especially precipitation) that control hydrologic partitioning in watersheds have not yet been reliably provided at the temporal resolution we desire [22].
In particular, the climate characteristics of the EAMC region can be defined as a marine monsoon climate in terms of precipitation and a continental climate in terms of temperature [23]. EAMC, which dominates the summer rainfall on the Korean peninsula, is highly influenced by the location and activity of polar fronts formed at the boundary between the air masses of the East Asian continent and oceanic air masses [24]. In particular, since more than two-thirds of annual precipitation occurs in summer, securing available water resources in summer is a very important social issue [25]. However, the projected future precipitation in the EAMC region is very difficult to attain credibility because it shows very large variability depending on what climate model is applied [26]. In fact, even when a model to reproduce typhoon-induced rainfall in the past is constructed, very different rainfall simulations are often obtained depending on how the physical options of the regional climate model are combined under the same initial conditions and boundary conditions [27,28]. Uncertainty of the climate model results weakens the reliability of future climate data (especially rainfall) composed of sub-daily resolution, and it is difficult to obtain the reliability of various future hydrological components derived by using a complex deterministic water balance model using such future uncertain climate data as input data [29].
From the viewpoint of modeling, due to the interaction of various geophysical processes that make up the soil water dynamics, the high non-linearity of each process, and the temporal and spatial variability of the hydrological and ecological environment itself, simplification assumptions at various accuracy levels are required [30]. It is based on the basic philosophy that describing an observed process need not be more complicated than necessary [31]. In other words, simple things are sometimes enough. The key motivation for this study is the observation that the relationship between yearly precipitation and yearly stream flow in the Andong dam watershed located in the southern part of the Korean peninsula in East Asia is very simple [32]. Despite the fact that the hydrological cycle of the watershed is composed of numerous interactions between climate, topography, and vegetation, the system's response to input can be explained by a relatively simple linear relationship. There are, of course, many limitations to minimal descriptions and models, since minimal models are lacking in a variety of processes to describe a number of situations outside of average conditions. For example, to illustrate the yearly stream flow of the Andong dam watershed, a model consisting of two parameters of the linear regression formula may not be the best. This implies that adequate conceptualization is required along with the necessary simplification to construct a minimal model [33,34]. From this point of view, the explanatory possibility of the parameters making up the model is another reason for Water 2020, 12, 25 3 of 22 constructing deterministic models. There are contradictory aspects between simplicity and explicability. For example, various types of Budyko curves developed to look at the annual water balance may not have a clear physical meaning of the applied parameters, depending on the structure of the models [35].
Given the difficulty of obtaining a reliable future precipitation at an appropriate resolution, a simple stochastic hydrologic partitioning model that includes realistic data requirements, broad applicability, accuracy to fit the purpose, and physically descriptive model parameters can be applied appropriately to understand how watershed systems respond to climate change. This approach is based on the expectation that the shortcomings resulting from the simple model structure can be complemented by stochastic considerations [36]. The purpose of this study is to develop a conceptual soil water balance model that simulates the hydrologic partitioning process (wetting, surface flow, vaporization, and groundwater recharge) at monthly or annual level in watershed units. Based on the partitioning processes of precipitation, that is, precipitation falling on the surface of a natural watershed first wets the soil layer, and if the surface is saturated, it is transported to the river in the form of surface runoff, and the water infiltrated into the soil layer is vaporized back into the atmosphere, and at the same time, it percolates into deeper soil layers and eventually flows out into the rivers in the form of subsurface flows, the proposed model was developed with a focus on making the precipitation partitioning process more intuitive. At this time, interception, surface depression, and subsurface flow were not included in the model structure. In other words, the focus was on vertical one-dimensional water balance, which determines hydrologic partitioning such as precipitation, surface flow, wetting, vaporization, and groundwater recharge. The main characteristics of the proposed East Asia Monsoon Climate-Stochastic Parsimonious Hydrologic Partitioning (EAMC-SPHP) model are as follows: (1) To clearly identify the hydrologic partitioning process in the watershed using two widely available climate variables, precipitation and potential evapotranspiration (PET) statistics, and watershed characteristics parameters. (2) To reflect seasonality of EAMC by explicitly considering seasonality of precipitation and PET.
(3) Parameterizing the watershed with five physically explicable land elements: -Surface runoff coefficient φ that governs primary partitioning of precipitation, that is, the precipitation is divided into overland flow and soil wetting. -Effective soil water storage depth nZ r . -Critical soil water value s * to control vaporization. -Saturated hydraulic conductivity K s and percolation exponent β controlling groundwater recharge.
(4) Minimum observational data for application of unmeasured areas. A deterministic ordinary differential equation conceptualizing the essential hydrologic partitioning process will be transformed into a stochastic partial differential equation with the state variable of the probability distribution function (PDF) of soil water combined with the statistical properties of precipitation and PET. Steady-state soil water PDF is used to obtain the Horton index representing the watershed. The sensitivity of the Horton index to the parameters of the constructed model is analyzed and the changes in the Horton index with seasonal characteristics of the EAMC area will be explored. Finally, we will look at future Horton index changes using various future climate projections.

Parsimonious Hydrologic Paritioning Model
The hydrologic partitioning model applied in this study was based on the following dynamic governing equations. where n is the soil porosity (dimensionless) and Z r is the thickness of the soil layer where the vegetation roots are located near the surface (L). The wetting W is the amount of water supplied to the soil from the precipitation P, and the vaporization V is the amount of soil water lost to the atmosphere depending on the soil and vegetation conditions. The groundwater recharge G is the amount of soil water lost from the near-surface soil layer to the deeper soil layer. The normalized soil water s is dimensionless and therefore nZ r s is the soil water depth stored in the soil layer near the surface (its dimension (L)). The dimensions of the hydrologic partitioning components P, V, W, and G are all (L/T). In Equation (1), t is time, which is based on day (that is, calculation time interval ∆t = 1-day) in this study. The daily precipitation P is assumed to follow the mixed exponential PDF as follows [29]: where P m is the average precipitation on the day of precipitation, and λ is the probability of precipitation on some day. Soil water wetting W by precipitation is simulated in one of three forms, depending on the soil water just before the precipitation and the current precipitation as follows [20]: where φ is the surface runoff coefficient and Z s is nZ r (1 − s)/∆t, the maximum free space in which the top soil layer can receive water just before precipitation. Therefore, the conditional PDF f W|S (w s)or of W is as follows: where w m is (1 − φ)P m and M s is defined as exp − Z s w m . This is shown in Figure 1. where is the soil porosity (dimensionless) and is the thickness of the soil layer where the vegetation roots are located near the surface (L). The wetting is the amount of water supplied to the soil from the precipitation , and the vaporization is the amount of soil water lost to the atmosphere depending on the soil and vegetation conditions. The groundwater recharge is the amount of soil water lost from the near-surface soil layer to the deeper soil layer. The normalized soil water s is dimensionless and therefore is the soil water depth stored in the soil layer near the surface (its dimension (L)). The dimensions of the hydrologic partitioning components , , , and are all (L/T). In Equation (1), is time, which is based on day (that is, calculation time interval Δ = 1-day) in this study.
The daily precipitation is assumed to follow the mixed exponential PDF as follows [29]: where is the average precipitation on the day of precipitation, and is the probability of precipitation on some day.
Soil water wetting by precipitation is simulated in one of three forms, depending on the soil water just before the precipitation and the current precipitation as follows [20]: where is the surface runoff coefficient and is (1 − )/Δ , the maximum free space in which the top soil layer can receive water just before precipitation. Therefore, the conditional PDF | ( | ) of is as follows: where is (1 − ) and is defined as exp [− ]. This is shown in Figure 1. Groundwater recharge is simulated as follows [29]: where is the saturated hydraulic conductivity and is the exponent of the groundwater recharge process.
Vaporization is conceptualized depending on the state of soil water as follows [20]: Groundwater recharge G is simulated as follows [29]: where K s is the saturated hydraulic conductivity and β is the exponent of the groundwater recharge process.
Vaporization V is conceptualized depending on the state of soil water as follows [20]: Water 2020, 12, 25

of 22
where PET is potential evapotranspiration, and s * is the critical soil water value to control vaporization. Vaporization V and groundwater recharge G are shown in Figure 2.
Water 2020, 12, x FOR PEER REVIEW 5 of 22 where is potential evapotranspiration, and * is the critical soil water value to control vaporization. Vaporization and groundwater recharge are shown in Figure 2.

Horton Index
According to Troch et al. [37], the Horton index can be defined as: where < • > is an expectation operator and is computed over a period of time (1-year in this study).
In areas that receive very little solar energy, most of the precipitation will be lost to surface runoff because of the lack of energy to allow vaporization, and conversely, in dry or semi-arid climate regions, river runoff will be much smaller than precipitation. For the EAMC region including the Korean Peninsula, < 1. It can be regarded as a measure of the water use efficiency of vegetation in general (Huxman et al., 2004) and can be treated as an eco-hydrologic index that indirectly shows the effect of changes in precipitation on vegetation growth [38][39][40]. Recently, a variety of studies have been conducted on factors affecting the Horton index [41][42][43][44].
Using Equations (2) and (3), one can write Equation (7) to express Horton exponents as follows: where ( ) is the steady-state PDF of soil water following the governing Equation (1). Therefore, the Horton index is summarized as follows: where = < >/ is Budyko index, defined as the ratio of potential evapotranspiration to precipitation [45]. As can be seen from Equation (10), it is possible to quantify the Horton index if the steady-state PDF of soil water conforming to the governing Equation (1) can be obtained.

Horton Index
According to Troch et al. [37], the Horton index H can be defined as: where < · > is an expectation operator and is computed over a period of time (1-year in this study).
In areas that receive very little solar energy, most of the precipitation will be lost to surface runoff because of the lack of energy to allow vaporization, and conversely, in dry or semi-arid climate regions, river runoff will be much smaller than precipitation. For the EAMC region including the Korean Peninsula, H < 1. It can be regarded as a measure of the water use efficiency of vegetation in general (Huxman et al., 2004) and can be treated as an eco-hydrologic index that indirectly shows the effect of changes in precipitation on vegetation growth [38][39][40]. Recently, a variety of studies have been conducted on factors affecting the Horton index [41][42][43][44].
Using Equations (2) and (3), one can write Equation (7) to express Horton exponents as follows: where p(s) is the steady-state PDF of soil water following the governing Equation (1). Therefore, the Horton index is summarized as follows: where B = < PET >/λP m is Budyko index, defined as the ratio of potential evapotranspiration to precipitation [45]. As can be seen from Equation (10), it is possible to quantify the Horton index if the steady-state PDF of soil water conforming to the governing Equation (1) can be obtained.

Soil Water PDF
Using the cumulant expansion theory [20,46,47], the Fokker-Planck Equation representing the temporal behavior of PDF of soil water p(s, t) from Equation (1) can be derived as follows: where p(s, t) is the PDF of soil water and η is ( is the covariance operator. Assuming that soil water s is given at time t, the terms necessary for the calculation of < η t > in Equation (11) are derived as follows: < V|s > = <PET> s * , for 0 < s ≤ s * .
An approximation to the advection and dispersion coefficients of the Fokker-Planck Equation (10) can be made by assuming that the relationship between η and ∂η ∂s is delta-correlation, and by using the fact that s is given at time t. Note that G and ∂G ∂s need not be considered in the calculation of the covariance terms in Equation (10) since they are deterministic term when s is given at time t. It is also assumed that W and V are mutually independent. Under this framework, the second advection coefficient reduces to where where θ is the scale of fluctuation of the daily precipitation time series [48]. In order to calculate C, σ 1 and σ 2 in Equations (16) and (17), the relationship between W and ∂W ∂s , and the relationship between V and ∂V/∂s need to be summarized as follows: 1 When W = 0,
In such a condition, the probability of ∂W ∂s = nZ r /∆t is 1 − λ, and the probability of ∂W ∂s = nZ r (P 1 − P o )/P 1 ∆t is λ, where P 1 is the day's precipitation and P o is the precipitation that occurred before ∆t. 3 When W = Z s ,

•
The probability of W = Z s at some day is λM s .

•
In such a condition, the probability of ∂W ∂s = nZ r /∆t is 1 − λ, and the probability of  (16) and (17) are derived as follows: In more concise form, Equation (11) can be written as follows: where A(s) and D(s) are called respectively the advection and dispersion coefficients, as the form of Equation (26) closely resembles the advection-dispersion equation. Equation (26) is basically a continuity equation and the state variable of the equation is the probability density. To obtain a steady-state PDF of soil water from Equation (26), we applied the Chang and Cooper scheme [49].
The key property of the Chang and Cooper scheme is that the numerical probability currents (that is, J = A(s)p(s, t) − D(s)∂p(s, t)/∂s) used to discretize the right-hand-side of Equation (26) are equal to zero when the distribution p(s, t) is equal to the equilibrium distribution p(s). Considering the equilibrium solution of the Fokker-Planck Equation (26), and assuming the advection and dispersion coefficient are independent of p(s), the solution of J = 0, which is the steady-state PDF of soil water under stochastic precipitation and potential evapotranspiration forcings, can be written as follows: where p o is a constant of integration. Therefore, the steady-state PDF of soil water can be obtained numerically using Equations (12)- (25). As can be seen from Equation (27), the steady-state PDF of soil water is expressed as a function of hydrologic meteorological statistics such as frequency of occurrence of precipitation, mean precipitation of precipitation days, scale of fluctuation of daily precipitation time series, and mean and variance of potential evapotranspiration and parameters related to the soil and vegetation of the watershed such as surface runoff coefficient, effective soil depth, saturated hydraulic coefficient and loss exponent, and critical soil water value controlled by vaporization. The derivation process of Horton index, which has been examined so far, is as follows: (1) The Horton index, which can represent the hydrologic partitioning of the watershed, can be estimated using Equation (10). However, to calculate Equation (10), a PDF of soil water is needed. (2) Construct the soil water balance Equation (1) in the watershed clearly expressing hydrologic partitioning. (3) Using cumulant expansion theory, Equation (27) can be derived from Equation (1) to calculate steady-state PDF of soil water. (4) Equations (10) and (27) can be linked to calculate the Horton index of the watershed.
In fact, using observed climate data or projected climate data under future climatic conditions, the corresponding soil water PDF can be obtained by running the deterministic model (i.e., Equation (1)) directly, and then Horton index can be calculated using the numerically obtained soil water PDF and Equation (10). Since the deterministic model also simulates wetting and vaporization directly, we could calculate Horton index directly from wetting and vaporization simulation data without calculating soil water PDF. However, in this study, the reliability of the future daily precipitation data itself was not relatively high, but the basic statistics such as the mean and the variance of the data were judged to be the reliable future prospects. Therefore, we tried to examine the change of future Horton index by constructing stochastic model that could obtain soil water PDF using basic statistics of climate data.

Study Area and Model Parameters Estimation
For the application of the proposed model, Andong dam watershed (exactly the upstream catchment of the dam), located in south side of Korean peninsula of 128 • 43 -129 • 12 in longitude and 36 • 32 -37 • 12 in latitude, was selected, and the geophysical and meteorological data of this site were used. The area of Andong dam watershed is 1584 km 2 . The mean annual temperature is 9.7 • C, and the mean annual precipitation is 1166 mm. The selected area has the typical inland climate of long dry cold winter and long hot humid summer.
In order to estimate parameters of the proposed deterministic model, the data from Korea Meteorological Administration (www.kma.go.kr) were used. For the calculation of the PET, one of the model input data, the data of surface air temperature, relative humidity, sun light, and wind speed from 1981 to 2018 were obtained, and the Penman-Monteith method was applied. The daily precipitation data from 1981 to 2018 was also used. Daily inflow into Andong dam from 1981 to 2018 was obtained from Korean water resources information system (www.wamis.go.kr).
The number of watershed parameters for the proposed model is five, including the surface runoff coefficient φ that divides the precipitation into overland flow and soil wetting, the soil water storage depth nZ r of the watershed, the critical soil water value s * to control vaporization, and saturated hydraulic conductivity K s and exponent β. Of these, nZ r was determined to be 25400/CN-I − 254 = 207.8 mm due to the fact that the CN-I value of the Andong dam watershed is 55. Equation 25400/CN-1 − 254 is a formula for the potential maximum retention depth (mm) of the NRCN-CN method found in almost all hydrologic textbooks (i.e., [50]). Since this parameterization is applied to many hydrologic models including SWAT and SWMM, it is considered to be applicable to Andong dam catchment. For reference, there are limitations, but NRCS-CN empirical equation developed for small catchments and hill slopes in the US has been the most widely used in modeling watersheds in research areas. Therefore, it was judged that such parameterization is realistic and reasonable for conceptualizing the one-dimensional vertical hydrologic partitioning in the study area. In order to estimate the remaining parameters, a genetic algorithm, which is one of the optimization techniques, was used to find the parameter combination that minimizes the objective function as follows: where N m is the number of observational monthly data (N m = 456 in this study) applied to the parameter estimation, and Q s m and Q o m are the monthly dam inflow calculated by the proposed model and observed monthly dam inflow, respectively. The proposed hydrologic partitioning dynamics have a threshold behavior as shown in Equation (3). In other words, a discontinuous governing equation is used to mathematically explain the wetting process. This discontinuity, however, adds unnecessary non-smoothness to the objective function for parameter estimation. Several studies recommend applying smoothing to threshold based model Equations (i.e., [51]). However, the model proposed in this study consists of relatively simple governing equations, and its purpose is to implement average hydrologic partitioning behavior on a yearly basis (or monthly and seasonal basis) using probabilistic techniques. The smoothing of the discontinuity was considered to be beyond the scope of this study. Thus, no model smoothing strategies were applied in this study. During parameter estimation procedure, the range of parameters was set as follows: φ = [0, 1]: all physically possible ranges, s * = [0, 1]: all physically possible ranges, K s = [30,270]: the range of NRCS hydrological soil groups Type A, B, and C (mm/day), β = [1, 10]: subjectively determined by authors as realistically possible range. We thought that the estimated parameters were all considered to be in the appropriate range. Looking at the estimated parameter values, φ = 0.26134, s * = 0.48844, K s = 87.5183 mm/day, and β = 4.7581, respectively. Using these parameters, the observed monthly dam inflow data from 1981 to 2018 were compared to the simulated data (see Figure 3). As can be seen from the figure, the modeled monthly dam inflows properly re-generated the observed inflows with R 2 of 0.8868 and NSE (Nash-Sutcliffe model efficiency coefficient [52]) of 0.8861.
to Andong dam catchment. For reference, there are limitations, but NRCS-CN empirical equation developed for small catchments and hill slopes in the US has been the most widely used in modeling watersheds in research areas. Therefore, it was judged that such parameterization is realistic and reasonable for conceptualizing the one-dimensional vertical hydrologic partitioning in the study area. In order to estimate the remaining parameters, a genetic algorithm, which is one of the optimization techniques, was used to find the parameter combination that minimizes the objective function as follows: where is the number of observational monthly data ( = 456 in this study) applied to the parameter estimation, and and are the monthly dam inflow calculated by the proposed model and observed monthly dam inflow, respectively. The proposed hydrologic partitioning dynamics have a threshold behavior as shown in Equation (3). In other words, a discontinuous governing equation is used to mathematically explain the wetting process. This discontinuity, however, adds unnecessary non-smoothness to the objective function for parameter estimation. Several studies recommend applying smoothing to threshold based model Equations (i.e., [51]). However, the model proposed in this study consists of relatively simple governing equations, and its purpose is to implement average hydrologic partitioning behavior on a yearly basis (or monthly and seasonal basis) using probabilistic techniques. The smoothing of the discontinuity was considered to be beyond the scope of this study. Thus, no model smoothing strategies were applied in this study. During parameter estimation procedure, the range of parameters was set as follows: = [0, 1]: all physically possible ranges, * = [0, 1]: all physically possible ranges, = [30,270]: the range of NRCS hydrological soil groups Type A, B, and C (mm/day), = [1,10]: subjectively determined by authors as realistically possible range. We thought that the estimated parameters were all considered to be in the appropriate range. Looking at the estimated parameter values, = 0.26134, * = 0.48844, = 87.5183 mm/day, and = 4.7581, respectively. Using these parameters, the observed monthly dam inflow data from 1981 to 2018 were compared to the simulated data (see Figure 3). As can be seen from the figure, the modeled monthly dam inflows properly re-generated the observed inflows with of 0.8868 and NSE (Nash-Sutcliffe model efficiency coefficient [52]) of 0.8861. For reference, the yearly mean wetting estimated using these parameters was 859.2 mm/year, and the yearly mean vaporization was 562.2 mm/year. Thus, the Horton index of the Andong dam watershed was 0.654. Figure 3 might not be sufficient for the evaluation of the model. For example, the simulation results in Figure 3 show significant scatter along the line. That is, the proposed model sometimes shows a tendency to underestimate (or overestimate) the monthly stream flow. There can be many causes. The proposed model does not have an explicit process for lateral subsurface flows, and a simple description of the study catchment as a single soil water tank might be the reason. However, For reference, the yearly mean wetting estimated using these parameters was 859.2 mm/year, and the yearly mean vaporization was 562.2 mm/year. Thus, the Horton index of the Andong dam watershed was 0.654. Figure 3 might not be sufficient for the evaluation of the model. For example, the simulation results in Figure 3 show significant scatter along the line. That is, the proposed model sometimes shows a tendency to underestimate (or overestimate) the monthly stream flow. There can be many causes. The proposed model does not have an explicit process for lateral subsurface flows, and a simple description of the study catchment as a single soil water tank might be the reason. However, as shown in Figure 3, the average monthly stream flow is relatively well reproduced using the minimum number of parameters. As mentioned earlier, the model representing the hydrologic partitioning process requires simplification assumptions at various levels of accuracy depending on the application. This is because it does not have to be more complex than describing the observed process. Constructing complex models can increase the reproducibility of unusual observations that deviate from average behavior. In return, however, one should bear the increased number of parameters to be estimated, including those that have no physical meaning or are difficult to measure routinely in practice. For the purpose of investigating hydrologic partitioning at the annual (or at best seasonal or monthly) level, the results in Figure 3 may not be unreasonable parameter estimates. However, for more accurate stream flow prediction, it is necessary to add more complex topographical components and components that are not considered in the current model structure.

Inter-Annual Variability of the Horton Index
The annual Horton index obtained by averaging wetting W and vaporization V simulated daily using the proposed deterministic model is shown in Figure 4 along with annual precipitation and Budyko index. For reference, the horizontal line in Figure 4 is the average value during the data period. A notable result is that the coefficient of variation of precipitation is about 23%, while the coefficient of variation of the Horton index is about 14%, which is about 60% of the coefficient of variation of precipitation. Looking at the range of minimum and maximum values relative to the mean, the precipitation was in the range of 0.59-1.6 for 38 years, while the Horton index was in the range of 0.75-1.2. These results make it possible to infer that the Horton index could be usefully used as an index representing the complex characteristics of meteorological-ecological-hydrological conditions in an area [53]. In the case of the Budyko index, which is a similar index, and which is widely applied in the meteorological field, the inter-annual variability is about 25% for 38 years and the inter-annual variability of the Budyko index is larger than that of the precipitation. Therefore, it would be difficult to utilize Budyko index as an index representing the hydro-meteorological characteristics of the Korean Peninsula with the characteristics of EAMC.
Water 2020, 12, x FOR PEER REVIEW 10 of 22 as shown in Figure 3, the average monthly stream flow is relatively well reproduced using the minimum number of parameters. As mentioned earlier, the model representing the hydrologic partitioning process requires simplification assumptions at various levels of accuracy depending on the application. This is because it does not have to be more complex than describing the observed process. Constructing complex models can increase the reproducibility of unusual observations that deviate from average behavior. In return, however, one should bear the increased number of parameters to be estimated, including those that have no physical meaning or are difficult to measure routinely in practice. For the purpose of investigating hydrologic partitioning at the annual (or at best seasonal or monthly) level, the results in Figure 3 may not be unreasonable parameter estimates. However, for more accurate stream flow prediction, it is necessary to add more complex topographical components and components that are not considered in the current model structure.

Inter-Annual Variability of the Horton Index
The annual Horton index obtained by averaging wetting W and vaporization V simulated daily using the proposed deterministic model is shown in Figure 4 along with annual precipitation and Budyko index. For reference, the horizontal line in Figure 4 is the average value during the data period. A notable result is that the coefficient of variation of precipitation is about 23%, while the coefficient of variation of the Horton index is about 14%, which is about 60% of the coefficient of variation of precipitation. Looking at the range of minimum and maximum values relative to the mean, the precipitation was in the range of 0.59-1.6 for 38 years, while the Horton index was in the range of 0.75-1.2. These results make it possible to infer that the Horton index could be usefully used as an index representing the complex characteristics of meteorological-ecological-hydrological conditions in an area [53]. In the case of the Budyko index, which is a similar index, and which is widely applied in the meteorological field, the inter-annual variability is about 25% for 38 years and the inter-annual variability of the Budyko index is larger than that of the precipitation. Therefore, it would be difficult to utilize Budyko index as an index representing the hydro-meteorological characteristics of the Korean Peninsula with the characteristics of EAMC. Another thing that can be seen in Figure 4 is the inverse relationship between the precipitation and the Horton index. The relationship between annual precipitation and annual Horton index can be seen in more detail from Figure 5. The cross-correlation coefficient between the annual precipitation and the Horton index was about −0.92, indicating strong correlation. From the regression analysis, the relationship between precipitation and Horton index was derived and the empirical relation was obtained as follows: Another thing that can be seen in Figure 4 is the inverse relationship between the precipitation and the Horton index. The relationship between annual precipitation and annual Horton index can be seen in more detail from Figure 5. The cross-correlation coefficient between the annual precipitation and the Horton index was about −0.92, indicating strong correlation. From the regression analysis, the relationship between precipitation and Horton index was derived and the empirical relation was obtained as follows: where P is the annual precipitation (mm) and HI is the Horton index calculated over the same period. When the precipitation increased, the Horton index decreased. When the precipitation decreased extensively, the Horton index became 1, and when the precipitation increased excessively, the Horton index converged to 0. The special relationship between the Horton index and PET was difficult to find (see Figure 5b). The relationship between the Horton index and the Budyko index shows a distinct pattern between the case where the Budyko index was less than 1 and the case where the Budyko index was larger than 1 (see Figure 5c). As can be inferred from Equation (7), an increase in wetting with increasing precipitation under PET (or vaporization), which was relatively constant (or independent of Horton index, see Figure 5b), would lead to a decrease in Horton index. In addition, from Figure 5c, Horton index could be treated as a re-scaled version of the Budyko index. where is the annual precipitation (mm) and HI is the Horton index calculated over the same period. When the precipitation increased, the Horton index decreased. When the precipitation decreased extensively, the Horton index became 1, and when the precipitation increased excessively, the Horton index converged to 0. The special relationship between the Horton index and PET was difficult to find (see Figure 5b). The relationship between the Horton index and the Budyko index shows a distinct pattern between the case where the Budyko index was less than 1 and the case where the Budyko index was larger than 1 (see Figure 5c). As can be inferred from Equation (7), an increase in wetting with increasing precipitation under PET (or vaporization), which was relatively constant (or independent of Horton index, see Figure 5b), would lead to a decrease in Horton index. In addition, from Figure 5c, Horton index could be treated as a re-scaled version of the Budyko index. The Budyko index has been used as an index representing the climatic conditions of the watershed, and has been discussed in the field of hydrology [35]. However, it has been pointed out in many studies that the high inter-annual variability of the Budyko index is the limit of the application of the Budyko index [42,53]. On the other hand, as the results of this study, the Horton index shows a relatively low inter-annual variability, suggesting that it could be a good alternative to the Budyko index. In addition, the Horton index is known to be an index that can be considered ecological aspects, because it can reflect the response of vegetation through the vaporization process [37]. However, the Horton index has the disadvantage that it cannot be measured by observation data alone. Examining Horton index through (eco-) hydrologic models may be an important research topic in the field of hydrology. In this respect, this study is also a small part of the attempt to solve such problems. The Budyko index has been used as an index representing the climatic conditions of the watershed, and has been discussed in the field of hydrology [35]. However, it has been pointed out in many studies that the high inter-annual variability of the Budyko index is the limit of the application of the Budyko index [42,53]. On the other hand, as the results of this study, the Horton index shows a relatively low inter-annual variability, suggesting that it could be a good alternative to the Budyko index. In addition, the Horton index is known to be an index that can be considered ecological aspects, because it can reflect the response of vegetation through the vaporization process [37]. However, the Horton index has the disadvantage that it cannot be measured by observation data alone. Examining Horton index through (eco-) hydrologic models may be an important research topic in the field of hydrology. In this respect, this study is also a small part of the attempt to solve such problems.

Steady-State Soil Water PDF
The performance of the proposed stochastic parsimonious hydrologic partitioning model was evaluated by comparison with the results from corresponding numerical simulations. For this purpose, parameters λ, P m , θ, < PET > and < PET 2 > representing the hydro-meteorological characteristics of the watershed were estimated from the data in Section 3.1. λ = 0.2721, P m = 11.74 mm/day, and θ = 1.728 day were calculated from the daily rainfall time series, respectively. < PET > = 2.692 mm/day and < PET 2 > = 9.259 (mm/day) 2 also were estimated from the daily potential evapotranspiration time series calculated from the various types of weather data using the Penman-Monteith method.
The steady-state PDF of soil water can be calculated from the stochastic parsimonious hydrologic partitioning model proposed in Section 2.3 using the estimated hydrologic meteorological statistics. This is compared with the histogram calculated from the simulated soil water time series by the deterministic hydrologic partitioning model described in Section 2.1 to verify the validity of the proposed method. Figure 6 shows the time series of the simulated soil water using a deterministic model. The steady-state PDF of the soil water calculated from the stochastic model and the histogram of the soil water calculated from the deterministic model is shown in Figure 7. The performance of the proposed stochastic parsimonious hydrologic partitioning model was evaluated by comparison with the results from corresponding numerical simulations. For this purpose, parameters λ , ,· , < > and < > representing the hydro-meteorological characteristics of the watershed were estimated from the data in Section 3.1. λ = 0.2721, = 11.74 mm/day, and = 1.728 day were calculated from the daily rainfall time series, respectively. < PET > = 2.692 mm/day and < > = 9.259 (mm/day) 2 also were estimated from the daily potential evapotranspiration time series calculated from the various types of weather data using the Penman-Monteith method. The steady-state PDF of soil water can be calculated from the stochastic parsimonious hydrologic partitioning model proposed in Section 2.3 using the estimated hydrologic meteorological statistics. This is compared with the histogram calculated from the simulated soil water time series by the deterministic hydrologic partitioning model described in Section 2.1 to verify the validity of the proposed method. Figure 6 shows the time series of the simulated soil water using a deterministic model. The steady-state PDF of the soil water calculated from the stochastic model and the histogram of the soil water calculated from the deterministic model is shown in Figure 7.    Figure 6b shows the lowest Horton index (0.504) in the year 2003. The year 1994 was the year of the most extreme drought in Korea (annual precipitation 878 mm). Even in the EAMC area, vegetation is in a state of stress on water in the same climate conditions as in 1994, and water use efficiency (i.e., Horton index) similar to that in the semi-arid region is shown [54]. For reference, the year 2003 was the year in which the The performance of the proposed stochastic parsimonious hydrologic partitioning model was evaluated by comparison with the results from corresponding numerical simulations. For this purpose, parameters λ , ,· , < > and < > representing the hydro-meteorological characteristics of the watershed were estimated from the data in Section 3.1. λ = 0.2721, = 11.74 mm/day, and = 1.728 day were calculated from the daily rainfall time series, respectively. < PET > = 2.692 mm/day and < > = 9.259 (mm/day) 2 also were estimated from the daily potential evapotranspiration time series calculated from the various types of weather data using the Penman-Monteith method. The steady-state PDF of soil water can be calculated from the stochastic parsimonious hydrologic partitioning model proposed in Section 2.3 using the estimated hydrologic meteorological statistics. This is compared with the histogram calculated from the simulated soil water time series by the deterministic hydrologic partitioning model described in Section 2.1 to verify the validity of the proposed method. Figure 6 shows the time series of the simulated soil water using a deterministic model. The steady-state PDF of the soil water calculated from the stochastic model and the histogram of the soil water calculated from the deterministic model is shown in Figure 7.    Figure 6b shows the lowest Horton index (0.504) in the year 2003. The year 1994 was the year of the most extreme drought in Korea (annual precipitation 878 mm). Even in the EAMC area, vegetation is in a state of stress on water in the same climate conditions as in 1994, and water use efficiency (i.e., Horton index) similar to that in the semi-arid region is shown [54]. For reference, the year 2003 was the year in which the   Figure 6b shows the lowest Horton index (0.504) in the year 2003. The year 1994 was the year of the most extreme drought in Korea (annual precipitation 878 mm). Even in the EAMC area, vegetation is in a state of stress on water in the same climate conditions as in 1994, and water use efficiency (i.e., Horton index) similar to that in the semi-arid region is shown [54]. For reference, the year 2003 was the year in which the most serious typhoon damage occurred in Korea, and the annual precipitation was the largest (annual precipitation 1864 mm). It can be seen that the soil water was almost saturated at the time of typhoon occurrence, and the lowest Horton index was recorded in the simulated period . On the other hand, the bar in Figure 7 represents the relative frequency of simulated daily soil water for 38 years, and the solid line with circles is PDF derived from Equation (26). It can be seen that the PDF derived using five meteorological variable statistics relatively well reproduced the statistical characteristics of numerically simulated soil water. For the numerically simulated soil water, the mean was 0.283 and the standard deviation was 0.124. For soil water derived from the proposed PDF, the mean was 0.281 and the standard deviation was 0.114. This agreement means that the stochastic model applied in this study reproduced the numerically simulated soil water behavior. Therefore, it would be possible to examine the behavior of the Horton index using the PDF of soil water obtained using the proposed stochastic model. As mentioned previously, since it is possible to directly simulate daily soil water using the deterministic model, it is also possible to obtain a soil water PDF using numerically simulated soil water time series. However, running the deterministic model requires daily climate date. In other words, if there is observed (or simulated) daily climate data, a soil water PDF can be obtained and thus Horton index can be calculated without the stochastic model. However, considering the variability of climate data (e.g., 10% increase or decrease in mean or 10% increase or decrease in variance), soil water PDFs using the proposed stochastic model can be useful for calculating the Horton index. The proposed stochastic model does not require the input of daily climate data. Since only representative statistics of the climate date are needed to drive the proposed stochastic model, the methodology of calculating Horton index from the soil water PDF using the proposed stochastic model can be useful for examining the effects of future climate change scenarios. Figure 8 shows the results of examining whether the proposed stochastic model reproduces the seasonality of the EAMC region. The precipitation of the EAMC region is concentrated in two-thirds of the annual precipitation in the summer, and the dry season lasts in the winter. The potential evapotranspiration is highest in spring. Thus, the EAMC region has a climatic character, with distinct characteristics of each of the four seasons (see Table 1). Therefore, whether this seasonality of climate can be clearly reflected in the model is also one of the important considerations. As can be seen in Figure 8, the proposed stochastic model can be used to derive the steady-state PDF of soil water that reflects the seasonal characteristics of the EAMC area. For reference, it can be found that the Horton index (H (probabilistic) in Table 1) calculated by the stochastic model as shown in Table 1 also reproduces the Horton index (H (numerical) in Table 1) calculated by the numerical model relatively well.

Sensitivity Analysis
As shown in the above analysis, the Horton index of a certain watershed is particularly sensitive to the change in precipitation, and it is considered that it is desirable to grasp the sensitivity level by two characteristics because of the intermittent characteristic possessed by precipitation. In other words, we tried to understand the sensitivity of the Horton index by dividing the two major processes that constitute the precipitation phenomenon into the frequency of precipitation and the amount of precipitation when it occurred. Figure 8 shows the sensitivity of the Horton index to meteorological variables. In addition, from the results of the annual fluctuation of precipitation and potential evapotranspiration (see Figure 4), the precipitation has an annual variability of 25%, and the potential evapotranspiration has about 10% annual variability. Therefore, the sensitivity of the Horton index to the meteorological variables was examined considering the range of variability. Figure 9a shows the changes in the Horton index for cases where precipitation was 25% less and 25% more than the average 38 years. The Horton index was in the range of 0.589-0.729, while the Budyko index was in the range of 0.674-1.124. This was consistent with the results in Figure 4. The effect of intermittent precipitation on the Horton index is shown in Figure 9b. While the total precipitation remained

Sensitivity Analysis
As shown in the above analysis, the Horton index of a certain watershed is particularly sensitive to the change in precipitation, and it is considered that it is desirable to grasp the sensitivity level by two characteristics because of the intermittent characteristic possessed by precipitation. In other words, we tried to understand the sensitivity of the Horton index by dividing the two major processes that constitute the precipitation phenomenon into the frequency of precipitation and the amount of precipitation when it occurred. Figure 8 shows the sensitivity of the Horton index to meteorological variables. In addition, from the results of the annual fluctuation of precipitation and potential evapotranspiration (see Figure 4), the precipitation has an annual variability of 25%, and the potential evapotranspiration has about 10% annual variability. Therefore, the sensitivity of the Horton index to the meteorological variables was examined considering the range of variability. Figure 9a shows the changes in the Horton index for cases where precipitation was 25% less and 25% more than the average 38 years. The Horton index was in the range of 0.589-0.729, while the Budyko index was in the range of 0.674-1.124. This was consistent with the results in Figure 4. The effect of intermittent precipitation on the Horton index is shown in Figure 9b. While the total precipitation remained constant, we examined the response of the Horton index in the situation that the probability of wet day, λ, increased (or decreased) by x% and the mean precipitation depth of wet day, P m , decreased (or increased) by x%. The Horton index was in the range of 0.626-0.684. The change in the Horton index to a 25% change in total precipitation was 0.14, and the change in the Horton index to a 25% change in the probability of a wet day (where total precipitation is constant) was 0.058. Therefore, in a rough sense, it can be seen that the relative contribution of the intermittency of precipitation to the sensitivity of the Horton index to precipitation was about 40%. Obviously, the Budyko index for the intermittency of precipitation was not found to be sensitive.
(or increased) by x%. The Horton index was in the range of 0.626-0.684. The change in the Horton index to a 25% change in total precipitation was 0.14, and the change in the Horton index to a 25% change in the probability of a wet day (where total precipitation is constant) was 0.058. Therefore, in a rough sense, it can be seen that the relative contribution of the intermittency of precipitation to the sensitivity of the Horton index to precipitation was about 40%. Obviously, the Budyko index for the intermittency of precipitation was not found to be sensitive. The sensitivity of the Horton index to PET was about 45% of the sensitivity to total precipitation as shown in Figure 9c, and the sensitivity of the Horton index to the second moment of daily PET was almost insignificant. Although not shown in the figure, the sensitivity of the Horton index to the scale of fluctuation of daily precipitation was not large. Figure 10 shows the sensitivity of the stead-state soil water PDF and its corresponding Horton index to land surface characteristics such as surface runoff coefficient and effective soil depth n . It can be recognized that the mode of stead-state soil water PDF moved in a decreasing direction as more surface runoff occurred during rainfall, and vegetation thus tended to use water more efficiently (i.e., increased the Horton index). The effect of effective soil depth, which is the range of influence of vegetation roots, on the stead-state soil water PDF mode was negligible and had a significant influence on the spread of stead-state soil water PDF. It was also found that the sensitivity of the Horton index to effective soil depth was limited. Obviously, the sensitivity of the Budyko index to land surface characteristics could not be considered in the proposed model. The sensitivity of the Horton index to PET was about 45% of the sensitivity to total precipitation as shown in Figure 9c, and the sensitivity of the Horton index to the second moment of daily PET was almost insignificant. Although not shown in the figure, the sensitivity of the Horton index to the scale of fluctuation of daily precipitation was not large. Figure 10 shows the sensitivity of the stead-state soil water PDF and its corresponding Horton index to land surface characteristics such as surface runoff coefficient φ and effective soil depth nZ r . It can be recognized that the mode of stead-state soil water PDF moved in a decreasing direction as more surface runoff occurred during rainfall, and vegetation thus tended to use water more efficiently (i.e., increased the Horton index). The effect of effective soil depth, which is the range of influence of vegetation roots, on the stead-state soil water PDF mode was negligible and had a significant influence on the spread of stead-state soil water PDF. It was also found that the sensitivity of the Horton index to effective soil depth was limited. Obviously, the sensitivity of the Budyko index to land surface characteristics could not be considered in the proposed model.  Figure 11 shows the sensitivity of the Horton index to vegetation and soil properties such as * , , and . The smaller the * and , and the larger the , the higher the Horton index. This sensitivity can be explained using Equation (1), which is the governing equation of soil water dynamics. That is, as the loss of soil water due to vaporization increases and the loss of soil water due to groundwater recharge decreases, the Horton index moves in an increasing direction. Given that the Horton index is an eco-hydrological proxy representing the water use efficiency of vegetation, the activation of vegetation by vaporization will cause the Horton index to increase. In addition, the reduction of groundwater recharge will increase the Horton index because vegetation makes the environment more favorable for using soil water [55]. In terms of the magnitude of the sensitivity, it can be seen that the sensitivity of the Horton index to * and is much greater than . In addition, the proposed model does not take into account the interaction between vegetation and soil properties and local meteorology, indicating that the sensitivity of the Budyko index to vegetation and soil properties was not significant.

Future Projection of the Horton Index
Global climate models (GCMs) have low spatial resolution. Therefore, in the case of a peninsula region with a small land area and a large marine influence such as Korea, it is difficult to use the GCMs outputs to project the future climate. In these regions, it is reasonable to look at future climate using outputs from regional climate models (RCMs), which dynamically down-scale the outputs of GCMs. In this study, we used RCM outputs (KOR-11) that dynamically down-scaled the East Asia region including the Korean Peninsula at a horizontal resolution of 12.5-km. In KOR-11, representative concentration pathways (RCP) 4.5 and 8.5 was applied to future climate change scenarios, and used two GCMs including MPI-ESM-LR (Max Plank Institute Earth System Model-Low Resolution, ML) and HadGEM2-AO (Hadly Center Global Environmental Model Version 2  Figure 11 shows the sensitivity of the Horton index to vegetation and soil properties such as s * , K s , and β. The smaller the s * and K s , and the larger the β, the higher the Horton index. This sensitivity can be explained using Equation (1), which is the governing equation of soil water dynamics. That is, as the loss of soil water due to vaporization increases and the loss of soil water due to groundwater recharge decreases, the Horton index moves in an increasing direction. Given that the Horton index is an eco-hydrological proxy representing the water use efficiency of vegetation, the activation of vegetation by vaporization will cause the Horton index to increase. In addition, the reduction of groundwater recharge will increase the Horton index because vegetation makes the environment more favorable for using soil water [55]. In terms of the magnitude of the sensitivity, it can be seen that the sensitivity of the Horton index to s * and β is much greater than K s . In addition, the proposed model does not take into account the interaction between vegetation and soil properties and local meteorology, indicating that the sensitivity of the Budyko index to vegetation and soil properties was not significant.  Figure 11 shows the sensitivity of the Horton index to vegetation and soil properties such as * , , and . The smaller the * and , and the larger the , the higher the Horton index. This sensitivity can be explained using Equation (1), which is the governing equation of soil water dynamics. That is, as the loss of soil water due to vaporization increases and the loss of soil water due to groundwater recharge decreases, the Horton index moves in an increasing direction. Given that the Horton index is an eco-hydrological proxy representing the water use efficiency of vegetation, the activation of vegetation by vaporization will cause the Horton index to increase. In addition, the reduction of groundwater recharge will increase the Horton index because vegetation makes the environment more favorable for using soil water [55]. In terms of the magnitude of the sensitivity, it can be seen that the sensitivity of the Horton index to * and is much greater than . In addition, the proposed model does not take into account the interaction between vegetation and soil properties and local meteorology, indicating that the sensitivity of the Budyko index to vegetation and soil properties was not significant.

Future Projection of the Horton Index
Global climate models (GCMs) have low spatial resolution. Therefore, in the case of a peninsula region with a small land area and a large marine influence such as Korea, it is difficult to use the GCMs outputs to project the future climate. In these regions, it is reasonable to look at future climate using outputs from regional climate models (RCMs), which dynamically down-scale the outputs of GCMs. In this study, we used RCM outputs (KOR-11) that dynamically down-scaled the East Asia region including the Korean Peninsula at a horizontal resolution of 12.5-km. In KOR-11, representative concentration pathways (RCP) 4.5 and 8.5 was applied to future climate change scenarios, and used two GCMs including MPI-ESM-LR (Max Plank Institute Earth System Model-Low Resolution, ML) and HadGEM2-AO (Hadly Center Global Environmental Model Version 2

Future Projection of the Horton Index
Global climate models (GCMs) have low spatial resolution. Therefore, in the case of a peninsula region with a small land area and a large marine influence such as Korea, it is difficult to use the GCMs outputs to project the future climate. In these regions, it is reasonable to look at future climate using outputs from regional climate models (RCMs), which dynamically down-scale the outputs of GCMs. In this study, we used RCM outputs (KOR-11) that dynamically down-scaled the East Asia region including the Korean Peninsula at a horizontal resolution of 12.5-km. In KOR-11, representative concentration pathways (RCP) 4.5 and 8.5 was applied to future climate change scenarios, and used two GCMs including MPI-ESM-LR (Max Plank Institute Earth System Model-Low Resolution, ML) and HadGEM2-AO (Hadly Center Global Environmental Model Version 2 coupled with the Atmosphere-Ocean, H2) and four RCMs (MM5, RegCM4, RSM, and WRF). Therefore, a total of 16 future climate information ensembles were produced.
Once the input meteorological variables λ, P m , θ, < PET >, and < PET 2 >, which are needed in the proposed model, were created for the control (1981-2010) and scenario (2021-2050) periods, future down-scaled variables were constructed by applying climate change factors as follows: where V s , V c , V o , and V F are the input meteorological variable (i.e., λ, P m , θ, and < PET >) for the scenario and control runs, the observations, and the perturbed observations representing the future climate conditions, respectively. In the case of < PET 2 >, the following climate change factor was applied in consideration of the dimension of the variable.
For reference, only five meteorological variables are needed in the proposed model. Figure 12 shows the future projections of five meteorological variables in the form of a box-plot.
Water 2020, 12, x FOR PEER REVIEW 17 of 22 coupled with the Atmosphere-Ocean, H2) and four RCMs (MM5, RegCM4, RSM, and WRF). Therefore, a total of 16 future climate information ensembles were produced. Once the input meteorological variables λ, , θ, < >, and < >, which are needed in the proposed model, were created for the control (1981-2010) and scenario (2021-2050) periods, future down-scaled variables were constructed by applying climate change factors as follows: where , , , and are the input meteorological variable (i.e., , , , and < >) for the scenario and control runs, the observations, and the perturbed observations representing the future climate conditions, respectively. In the case of < >, the following climate change factor was applied in consideration of the dimension of the variable.
For reference, only five meteorological variables are needed in the proposed model. Figure 12 shows the future projections of five meteorological variables in the form of a box-plot. In all climate information ensembles, more precipitation occurred during the rainy days, but the probability of precipitation was projected to decrease. This means that the pattern of future precipitation was likely to proceed toward more intermittent and concentrated precipitation. Since the scale of fluctuation of the future daily precipitation time series appeared to decrease, the temporal variability of precipitation was expected to increase further. In other words, it can be seen that future precipitation would change in the direction of increasing variability in terms of precipitation occurrence, intensity and persistence. Annual mean precipitation was projected to increase by 0%-8.1%, and annual potential evapotranspiration was expected to increase by 2.5%-7.1% to reflect rising surface air temperatures. Figure 13 shows the steady-state soil water PDFs simulated using future meteorological variables. For reference, the thick line is the soil water PDF simulated using observed meteorological variables, and the thin line is the future soil water PDF simulated using future meteorological variables derived from various climate model combinations. At this time, it is assumed that the land surface characteristics do not change. One can find that the mode of PDF of future soil water is not very different from the mode of present PDF, but depending on which GCM-RCM-RCP combination is applied, the dispersion of future soil water is increased or decreased more than the present dispersion without consistent direction. In all climate information ensembles, more precipitation occurred during the rainy days, but the probability of precipitation was projected to decrease. This means that the pattern of future precipitation was likely to proceed toward more intermittent and concentrated precipitation. Since the scale of fluctuation of the future daily precipitation time series appeared to decrease, the temporal variability of precipitation was expected to increase further. In other words, it can be seen that future precipitation would change in the direction of increasing variability in terms of precipitation occurrence, intensity and persistence. Annual mean precipitation was projected to increase by 0%-8.1%, and annual potential evapotranspiration was expected to increase by 2.5%-7.1% to reflect rising surface air temperatures. Figure 13 shows the steady-state soil water PDFs simulated using future meteorological variables. For reference, the thick line is the soil water PDF simulated using observed meteorological variables, and the thin line is the future soil water PDF simulated using future meteorological variables derived from various climate model combinations. At this time, it is assumed that the land surface characteristics do not change. One can find that the mode of PDF of future soil water is not very different from the mode of present PDF, but depending on which GCM-RCM-RCP combination is applied, the dispersion of future soil water is increased or decreased more than the present dispersion without consistent direction. Water 2020, 12, x FOR PEER REVIEW 18 of 22 Figure 13. Future soil water PDF. Figure 14 shows the degree of change in the future Budyko and Horton indexes. The present Budyko index was 0.843, and the rage of change in the future Budyko index was projected to range from −2.13% to +7.78%. The ensemble average Budyko index using 16 future climate ensembles was 0.862, with a potential increase of 2.27%. In other words, the future climate of the study area was likely to be slightly drier than the present climate. On the other hand, the present Horton index was 0.655, and the rate of change in the future Horton index was projected to change in the range of −1.96% to +4.09%. The average of the future ensemble of Horton index, simulated using 16 future climate ensembles, was 0.663, which means that it was likely to have a slightly larger value than the present. Future precipitation in the applied area increased, but Horton index was likely to increase further. Increasing Horton index means that vegetation was required to use water in the soil more efficiently than it was now. Thus, an increase in Horton index could be interpreted to indicate that the stress the vegetation was subjected to by water was likely to increase slightly more than now. This result revealed that it was a wrong approach to interpret future increases in precipitation as a reduction in vegetation water stress in the watershed from the relationship in Figure 5. Due to the variation in precipitation patterns and the increase in potential evapotranspiration, the water stress of vegetation in the watershed was more likely to increase. However, there was a clear limitation that the results of this study were based on steady-state soil water PDF. In other words, the current and future soil water PDFs were relatively similar, but when analyzed seasonally or monthly (see [29]), it was expected that other interesting facts would be revealed. For this analysis, the dynamic behavior of soil water PDF should be analyzed, so further research will need to be conducted based on this study.
In addition, this study focused on vertical one-dimensional water balance, which determines hydrologic partitioning such as precipitation, wetting, vaporization, and groundwater recharge,  Figure 14 shows the degree of change in the future Budyko and Horton indexes. The present Budyko index was 0.843, and the rage of change in the future Budyko index was projected to range from −2.13% to +7.78%. The ensemble average Budyko index using 16 future climate ensembles was 0.862, with a potential increase of 2.27%. In other words, the future climate of the study area was likely to be slightly drier than the present climate. On the other hand, the present Horton index was 0.655, and the rate of change in the future Horton index was projected to change in the range of −1.96% to +4.09%. The average of the future ensemble of Horton index, simulated using 16 future climate ensembles, was 0.663, which means that it was likely to have a slightly larger value than the present. Future precipitation in the applied area increased, but Horton index was likely to increase further. Increasing Horton index means that vegetation was required to use water in the soil more efficiently than it was now. Thus, an increase in Horton index could be interpreted to indicate that the stress the vegetation was subjected to by water was likely to increase slightly more than now. This result revealed that it was a wrong approach to interpret future increases in precipitation as a reduction in vegetation water stress in the watershed from the relationship in Figure 5. Due to the variation in precipitation patterns and the increase in potential evapotranspiration, the water stress of vegetation in the watershed was more likely to increase.  Figure 14 shows the degree of change in the future Budyko and Horton indexes. The present Budyko index was 0.843, and the rage of change in the future Budyko index was projected to range from −2.13% to +7.78%. The ensemble average Budyko index using 16 future climate ensembles was 0.862, with a potential increase of 2.27%. In other words, the future climate of the study area was likely to be slightly drier than the present climate. On the other hand, the present Horton index was 0.655, and the rate of change in the future Horton index was projected to change in the range of −1.96% to +4.09%. The average of the future ensemble of Horton index, simulated using 16 future climate ensembles, was 0.663, which means that it was likely to have a slightly larger value than the present. Future precipitation in the applied area increased, but Horton index was likely to increase further. Increasing Horton index means that vegetation was required to use water in the soil more efficiently than it was now. Thus, an increase in Horton index could be interpreted to indicate that the stress the vegetation was subjected to by water was likely to increase slightly more than now. This result revealed that it was a wrong approach to interpret future increases in precipitation as a reduction in vegetation water stress in the watershed from the relationship in Figure 5. Due to the variation in precipitation patterns and the increase in potential evapotranspiration, the water stress of vegetation in the watershed was more likely to increase. However, there was a clear limitation that the results of this study were based on steady-state soil water PDF. In other words, the current and future soil water PDFs were relatively similar, but when analyzed seasonally or monthly (see [29]), it was expected that other interesting facts would be revealed. For this analysis, the dynamic behavior of soil water PDF should be analyzed, so further research will need to be conducted based on this study.
In addition, this study focused on vertical one-dimensional water balance, which determines hydrologic partitioning such as precipitation, wetting, vaporization, and groundwater recharge, However, there was a clear limitation that the results of this study were based on steady-state soil water PDF. In other words, the current and future soil water PDFs were relatively similar, but when analyzed seasonally or monthly (see [29]), it was expected that other interesting facts would be revealed. For this analysis, the dynamic behavior of soil water PDF should be analyzed, so further research will need to be conducted based on this study.
In addition, this study focused on vertical one-dimensional water balance, which determines hydrologic partitioning such as precipitation, wetting, vaporization, and groundwater recharge, rather than horizontal water movement such as subsurface flow. Therefore, there was a limit to the investigation of future water resources availability using the method proposed in this study. However, it was expected to contribute to improving understanding of hydrologic partitioning in one-dimensional vertical water movement and to investigate the degree of change in rainwater use efficiency of the catchment's ecosystems under future climate change scenario conditions.

Conclusions
In this study, a conceptual soil water balance model was proposed to simulate the hydrologic partitioning process (wetting, surface flow, vaporization, and groundwater recharge) of precipitation at the monthly or annual time scale in natural watersheds. Using the model, the variability of the Horton index was investigated. After estimating the appropriate model parameters for vegetation and soil in the study watershed, the observed climate data were entered into the model to numerically simulate the major hydrological components of the watershed (surface runoff, wetting, vaporization, groundwater recharge, and soil water). The simulations yielded Horton index with inter-annual variability that is much less than the inter-annual variability of annual precipitation and Budyko index. This means that Horton index can be usefully applied as a surrogate variable representing hydro-meteorological characteristics in the East Asia monsoon region. The numerical results also show a strong inverse correlation between annual precipitation and Horton index. Recalling that the Horton index is related to the efficiency with which vegetation uses precipitation in watersheds, this inverse correlation indicates that vegetation in natural watersheds tended to survive by increasing utilization efficiency when available precipitation decreased. This was consistent with several studies.
Meanwhile, the deterministic ordinary differential equation, which parsimoniously conceptualizes the hydrologic partitioning process, was converted to the Fokker-Planck equation using the cumulant expansion theory with the probability distribution function of soil water as the state variable. The converted Fokker-Planck equation takes the form of a typical advection-dispersion partial differential equation, with advection and dispersion coefficients accounted for the statistics of precipitation and latent evaporation and hydro-geophysical characteristics. Using the steady-state soil water PDF, the Horton index representing the watershed was obtained. In addition, it was confirmed that the proposed stochastic model satisfactorily yields the steady-state PDF of soil water that reflects the distinct seasonality of the East Asia monsoon climate region. From the sensitivity analysis of the proposed stochastic model, it was found that the relative contribution of the intermittent nature of precipitation to the Horton index's overall sensitivity to precipitation in the study watershed was about 40%. Besides, the effect of potential evapotranspiration on the Horton index was less than half the effect of precipitation. Sensitivity analysis of Horton index on land surface characteristics revealed that surface runoff coefficient, critical soil water value controlling vaporization, and exponent value controlling groundwater recharge had the almost same sensitivity as that of the Horton index to precipitation. This means that unlike Budyko index, Horton index can be used as an eco-hydrologic surrogate variable that reflects the hydro-geophysical characteristics of the watershed.
The 16 future climate change ensembles available for the study watershed showed that future precipitation is likely to increase. However, the results of applying the future climate ensemble to the proposed model revealed that even if precipitation increases, the water stress of vegetation in the watershed would increase due to the increased variability of precipitation patterns and the increase in potential evapotranspiration.
Finally, it is important to note that this study examined the watershed response to climate forcings using the Horton index, but the Horton index is highly dependent on the results of the model used.
In addition, although this study was based on the results of several studies that the Horton index can be used as an eco-hydrological surrogate in a watershed, it is still necessary to discuss whether the Horton index can really be such a surrogate. In view of the uncertainty of the model, insufficient support of the results from this study could be pointed out as a matter of supplementation. Uncertainties from various future climate change ensembles, uncertainties from model structures (or parameters), and uncertainties from short observations also need to be investigated. Nevertheless, in this study, a novel methodology that examines the watershed response to climate change by deriving the Horton index from the stochastic analysis of a model that explicitly implements hydrologic partitioning in the East Asian monsoon climate region with minimal parameters was proposed. In this regard, the academic contribution of this study was expected to be given.