Evaluation of Land-Use Changes Impact on Watershed Health Using Probabilistic Approaches

: This study evaluated watershed health (WH) change using reference values for environmental changes at various times. Land use in 1985 was deﬁned as the reference value under the most natural conditions, and the WH for the years 1995 to 2019 was calculated in comparison to 1985. The proposed method was used to assess the WH of 78 standard subbasins in South Korea’s Geum River Basin (GRB), where complex land-use change has occurred since 1995. For evaluating hydrology and water quality (WQ) health index, Soil and Water Assessment Tool (SWAT) and four land-use maps (1985, 1995, 2008, and 2019) were used to simulate the hydrology and WQ. A multivariate normal distribution (MND) from poor (0) to good (1) was used to assess WH based on SWAT modeling results. Based on the reference value, the WQ health from 1995 to 2019 changed to within 0.1, while the range of changes in the hydrology index was analyzed over 0.18. As a result of WH changes from 1985 to 2019, hydrological health deteriorated in high-density urbanized subbasins, while WQ health deteriorated in upland-cultivation-increased subbasins. This study provides useful information for recognizing potential WH issues related to long-term environmental changes. GRB deteriorates over time as a result of human activities (i.e., urbanization and upland reclamation). The ﬁndings of this study indicate that the derived approach has the advantage of providing an absolute comparison of WH change. Based on such an analysis, it can be extended to the analysis of WH for various indicators and data, such as climate change, socioeconomic change, vulnerabilities, and recovery priorities of the basin. To conclude, this study can provide the decision-making data for basin management.


Introduction
Watersheds have been the bases of many local communities and an important source of livelihood for people [1]. Although many countries have achieved rapid economic growth, water quality (WQ) and quantity and the ecological environmental conditions in watersheds are facing serious problems [2,3]. Recently, stream watersheds have undergone increases in impermeable areas and water intake due to climate change and urbanization, and these watersheds are also affected by river improvement projects. In this situation, the exiting normal water cycle system of the watershed collapses, lowering the groundwater level and reducing the stream flow. Consequently, streams lose their normal functions. In addition, the shortage of maintenance flow caused by the drying of streams results in an increase in water pollutants, especially in streams [4][5][6][7].
Such a change in the watershed environment and the inefficient utilization of water resources seriously distort the hydrologic health of watershed. Increased surface flow, decreased groundwater flow, increased flooding, soil erosion, water pollution, and urban desertification are all side effects of this crisis [8]. Accordingly, to maintain a sustainable and healthy water cycle, it is desirable to accurately identify the overall water cycle systems of each watershed and develop an analysis technology that enables the water availability

Study Area
The Geum River Basin (GRB) is located in the most western area of South Korea (35 • 08 00" N-37 • 09 00 N and 126 • 07 00 E-128 • 00 00 E). The basin serves 12% of South Korea's population and has the largest proportion of land use dedicated for agriculture (29.5% of the basin area). The metropolitan city of Daejeon occupies midstream of the Geum River and urbanization has expanded around this city. The downstream is an important agricultural production area and in recent years, the rice paddy area has been changed to upland crop area. The water management issues in this basin include water scarcity and high rates of water pollution incidents in comparison with other basins. In the GRB, two multipurpose dams (Daecheong (DCD) and Yongdam (YDD)) and three multifunctional weirs (Backje (BJW), Gongju (GJW), and Sejong (SJW)) have been installed and are currently in operation. GRB is under temperature monsoon climates; annually, the temperature is 12.2 • C. Annual rainfall is 1261 mm and 50-60% of rainfall has occurred during summer (June-August) for 35 years (1985-2019). Figure 1a shows GRB and 78 standard subbasins that were designed for modeling of SWAT and WH evaluation, and Figure 1b illustrates two mid-subbasins for comparison of the WH changes by land use. The hydrology and WQ health of 78 standard subbasins were evaluated based on the SWAT simulation results. In addition, two subbasins with rapid changes in land use among the evaluation results were analyzed in detail.

Study Area
The Geum River Basin (GRB) is located in the most western area of South Korea (35°08′00″ N-37°09′00″ N and 126°07′00″ E-128°00′00″ E). The basin serves 12% of South Korea's population and has the largest proportion of land use dedicated for agriculture (29.5% of the basin area). The metropolitan city of Daejeon occupies midstream of the Geum River and urbanization has expanded around this city. The downstream is an important agricultural production area and in recent years, the rice paddy area has been changed to upland crop area. The water management issues in this basin include water scarcity and high rates of water pollution incidents in comparison with other basins. In the GRB, two multipurpose dams (Daecheong (DCD) and Yongdam (YDD)) and three multifunctional weirs (Backje (BJW), Gongju (GJW), and Sejong (SJW)) have been installed and are currently in operation. GRB is under temperature monsoon climates; annually, the temperature is 12.2 °C . Annual rainfall is 1261 mm and 50-60% of rainfall has occurred during summer (June-August) for 35 years (1985-2019). Figure 1a shows GRB and 78 standard subbasins that were designed for modeling of SWAT and WH evaluation, and Figure 1b illustrates two mid-subbasins for comparison of the WH changes by land use. The hydrology and WQ health of 78 standard subbasins were evaluated based on the SWAT simulation results. In addition, two subbasins with rapid changes in land use among the evaluation results were analyzed in detail.

WH Evaluation Method
Hydrologic factors (surface runoff (SQ), ET, SW, and groundwater flow (GWQ)) and WQ factors (suspended solids (SS), total phosphorus (T-P), and total nitrogen (T-N)) can be derived by the following equations for WH (Equations (1)-(3)), which were proposed by Ahn and Kim [18,19]. For each of the four hydrologic factors and three WQ factors, Equation (1) was utilized to derive the normalized component values. Based on these values, two subindices for hydrology and WQ were obtained (Equation (2)). A comprehensive health evaluation for hydrology and WQ was conducted by deriving one index from the two subindices and Equation (3). The normalized values were calculated based on percentile ranks; each index ranged from 0 to 1. An index closer to 1 represents a healthier

WH Evaluation Method
Hydrologic factors (surface runoff (SQ), ET, SW, and groundwater flow (GWQ)) and WQ factors (suspended solids (SS), total phosphorus (T-P), and total nitrogen (T-N)) can be derived by the following equations for WH (Equations (1)-(3)), which were proposed by Ahn and Kim [18,19]. For each of the four hydrologic factors and three WQ factors, Equation (1) was utilized to derive the normalized component values. Based on these values, two subindices for hydrology and WQ were obtained (Equation (2)). A comprehensive health evaluation for hydrology and WQ was conducted by deriving one index from the two subindices and Equation (3). The normalized values were calculated based on percentile ranks; each index ranged from 0 to 1. An index closer to 1 represents a healthier The MND is a generalization of the univariate normal distribution. For correlated variable random vectors, it is a distribution with a univariate normal distribution for each vector element. There is no connection between variables in the simplest case, and the vector components are separate univariate normal random variables. The MND has the mean vector µ and the covariance matrix Σ as parameters. These are similar to the mean parameter µ and the variance parameter σ 2 of the univariate normal distribution. The diagonal element of Σ includes the variance for each parameter, while the nondiagonal element of Σ includes the covariance of variables. As a megavariate distribution that can be easily dealt with, the MND is usually used as a model for multivariate data [23]. If the bivariate normal distribution is defined as two discrete random variables X and Y, which are different from each other, the bivariate probability function expressed as P (X = x, Y = y) shows a joint probability where the random variables X and Y have the particular values of x and y, respectively.
where, x = x 1 x 2 is the variable following the bivariate normal distribution, µ = µ 1 µ 2 is the mean value of each variable, σ x 1 is the variance of the variable x 1 , σ x 2 is the variance of , and E denotes the mean value [24,25]. To improve the existing evaluation index for WH, which is applicable only to relative evaluations, this study calculated the reference value by applying the MND probability function. Based on the SWAT modelling results using the 1985 land use, the hydrology and WQ health were evaluated using Equations (1)-(3). The evaluation results were defined as the reference values under the most natural conditions. The SWAT simulation results, which were necessary to derive WH indices from the reference values, were set as variables. A python program was implemented for the variables. Thus, a probability density function was derived from the MND [26].

SWAT Description
The SWAT, which was released by the USDA-ARS [27], was conducted to simulate the long-term runoff from the GRB. The SWAT is a physically dependent continuous model that can simulate hydrology and WQ for various types of land use. The model simulates SQ, SW, and GWQ using the water balance equation for HRU. The SCS-CN method is used in SWAT to calculate daily streamflow. The MUSLE was used to simulate phosphorus and nitrogen transport and concentration estimation. The water bodies considered included streams, dams and weirs, and streamflow; SS and nutrients were considered [28].  (Figure 2). WAMIS classifies land cover into the following eight categories: urban, crop, paddy, grass, forest, wetland, bare field, and water [29]. The land-use map was composed of Landsat TM, Landsat ETM, and Korea Multi-Purpose Satellite-2 (KOMPSAT-2) images. Before analyzing the hydrologic conditions and WQ according to the change in land use, this study gathered seven categories of land-use maps from the Water Resources Management Information System (WAMIS) for the years 1985 and 1995, and land-use maps of 2008 and 2019 were provided by the Environmental Geographic Information Service (EGIS) (Figure 2). WAMIS classifies land cover into the following eight categories: urban, crop, paddy, grass, forest, wetland, bare field, and water [29]. The land-use map was composed of Landsat TM, Landsat ETM, and Korea Multi-Purpose Satellite-2 (KOMPSAT-2) images. (b) Figure 2. Land-use data of (a) Geum River Basin and (b) two test subbasins.

GIS and Measured Data for SWAT Evaluation
The 30 m DEM and soil map from WAMIS were conducted as the GIS spatial data for the SWAT. Sandy loam occupied high proportions of 58% of the soil types in GRB. The 35-year period (1985-2019) SWAT weather inputs were daily data from nine weather stations [30].
The two dams (DCD and YDD) and three weirs (BJW, GJW, and SJW) operation data were collected from Korea Water Resources Corporation (K-water). The dam storage and daily inflow data for DCD collected over 35 years were used. For YDD, the dam daily inflow and storage data were collected from January 2002 to December 2019. For BJW, GJW, and SJW, the daily weir storage and inflow data were collected from August 2012 to December 2019. The groundwater level data observed by the National Groundwater Information Center (GIMS) groundwater level observation stations were collected. Five sites (BYBY, JSJS, BEMR, CASS, and OCCS) in the GRB were selected (Figure 1a) [31]. The flux tower DY, operated by K-water, is in the upper GRB within the forested area [31]. The DY daily ET (mm) and soil moisture (%) data were collected from April 2011 to December 2019. For 11 years (2005-2019), dam weir, groundwater, ET, and soil moisture were used to calibrate and validate the watershed hydrology.
Point source data for SWAT and daily sewage treatment plant monitoring data were aggregated from the National Institute of Environmental Research (NIER).

GIS and Measured Data for SWAT Evaluation
The 30 m DEM and soil map from WAMIS were conducted as the GIS spatial data for the SWAT. Sandy loam occupied high proportions of 58% of the soil types in GRB. The 35-year period (1985-2019) SWAT weather inputs were daily data from nine weather stations [30].
The two dams (DCD and YDD) and three weirs (BJW, GJW, and SJW) operation data were collected from Korea Water Resources Corporation (K-water). The dam storage and daily inflow data for DCD collected over 35 years were used. For YDD, the dam daily inflow and storage data were collected from January 2002 to December 2019. For BJW, GJW, and SJW, the daily weir storage and inflow data were collected from August 2012 to December 2019. The groundwater level data observed by the National Groundwater Information Center (GIMS) groundwater level observation stations were collected. Five sites (BYBY, JSJS, BEMR, CASS, and OCCS) in the GRB were selected (Figure 1a) [31]. The flux tower DY, operated by K-water, is in the upper GRB within the forested area [31]. The DY daily ET (mm) and soil moisture (%) data were collected from April 2011 to December 2019. For 11 years (2005-2019), dam weir, groundwater, ET, and soil moisture were used to calibrate and validate the watershed hydrology.
Point source data for SWAT and daily sewage treatment plant monitoring data were aggregated from the National Institute of Environmental Research (NIER). From 2005 to 2019, SS, T-P, and T-N weekly data were collected from three WQ monitoring stations (Jewon (JW), Gomnaru (GNR), and Jeongdong (JD)) for evaluation of stream WQ. The fertilizers applied were the measured fertilizers of upland crops and paddy rice in South Korea [32].

Data Reconstruction for Evaluating Hydrology and Water Quality Health
The EPA [17] uses the ratio of the water storage of a dam to the annual average stream flow before the development of the watershed to evaluate the WH at the side of the flood gate. However, the water resource characteristics of watersheds in South Korea require that stream flow, ET, SW, and GWQ be considered very important factors. Accordingly, this study evaluated WH by considering all hydrologic and water cycle factors for watershed management focusing on water quantity. This study utilized the SWAT verification and calibration results of Lee et al. [31] to evaluate hydrologic WH. The evaluation results were analyzed with respect to the following four factors: ET, SQ, SW, and GWQ.
For the WH evaluation with respect to WQ, the distributions of the SS, T-P, and T-N were calculated as evaluation factors for each standard watershed. Among the SWAT simulation result files, Output.rch was used to arrange the WQ simulation results. As Output.rch converted the simulation results of WQ to loads (ton/day or kg/day), the concentrations (mg/L) were calculated based on the simulation result of stream flow.

SWAT Evaluation for Hydrology and Water Quality
This study validated and calibrated SWAT before simulating the hydrology and WQ according to the land-use trend from 1985 to 2019. For the hydrologic validation and calibration of SWAT for the GRB, the parameters proposed by Lee et al. [31] were applied. As shown in Table 1, validation and calibration were conducted for two dams (YDD and DCD), and three weirs (SJW, GJW, and BJW). The NSE of dam inflow was analyzed and found to range from 0.57 to 0.67. The RMSE was 0.82 mm/day-1.66 mm/day. The PBIAS ranged from −5.65% to 8.02%. Thus, the calibration results were consistent with the SWAT calibration guidelines (NSE ≥ 0.5, PBIAS ≤ 28.0%, and R 2 ≥ 0.6) and were found to be satisfactory [19]. Based on the hydrologic validation and calibration results, the WQ at 3 points (JW, GNR, JD) in GRB was validated and calibrated. For the validation and calibration of WQ, the stream WQ measurement data between 2005 and 2019 were collected from the water environment information system. As nutritive substances are affected by the movement and discharge of sediment load in the SWAT, the SS load was calibrated before nutritive substances were calibrated [33,34]. The applicability of the SWAT applicability was evaluated using the R 2 . It was calculated that the statistical values of SS, T-P, and T-N were 0.72-0.81, 0.64-0.78, and 0.69-0.76, respectively (Table 1).

Land-Use Change Analysis
In GRB, there are dominant land-use types of area forest and paddy, which encompass up to 73% of the total area over the 35 years. In 1985, forests accounted for 63.6% of the total area, and this value decreased to 54.1% by 2019. In contrast, the urban area increased every year and accounted for 9.0% in 2019 ( Table 2). Comparing the land-use data of 1985,1995,2008, and 2019, the obvious change occurred in urban, agricultural, and forest land-use classes. Figure 3 represents the proportions of urban and upland crop area for each of the 14 mid-subbasins in the GRB. The mid-subbasins where urban area expanded by more than 5% compared to 1985 were 3009, 3010, 3011, 3012, 3013, and 3014 mid-subbasins ( Figure 3a). In particular, the urban area of the 3009 mid-subbasins increased by up to 12.3% compared to 1985. From 1985 to 2019, the upland crop area increased by more than 5% at eight mid-subbasins, and in 3013 mid-subbasin, where the change was the greatest, 11.6% of the watershed was converted into upland crop. (Figure 3b). for each of the 14 mid-subbasins in the GRB. The mid-subbasins where urban area ex-panded by more than 5% compared to 1985 were 3009, 3010, 3011, 3012, 3013, and 3014 mid-subbasins (Figure 3a). In particular, the urban area of the 3009 mid-subbasins increased by up to 12.3% compared to 1985. From 1985 to 2019, the upland crop area increased by more than 5% at eight mid-subbasins, and in 3013 mid-subbasin, where the change was the greatest, 11.6% of the watershed was converted into upland crop. (Figure  3b).   The areas of urban and upland crop increased while the area of forest and paddy decreased. As the forest area was generally reduced for the 14 mid-subbasins as shown in Figure 2, in this study, we selected the test areas where urbanized and upland crop cultivation increased rapidly. Table 2 describes occupied land use classification in the most urbanized (3009) and upland-cultivation-increased (3013) subbasins. The urban area increases in 3009 subbasin was at the expense of paddy field. Daejeon city is located in this subbasin, and the city was expanded to high density every year, resulting in the largest increase in urban area within the total watershed. The urban expanded rate in 3009 mid-subbasin was 75% higher than the average expanded rate in the GRB. The 3013 mid-subbasin is a representative agricultural area in the GRB. This subbasin is a representative agricultural area in the GRB, with paddy area accounting for 42.8% of the basin in 1985. The paddy area of 2019 decreased by 12.8% while upland crop area increased by 11.6%. This change is due to policies that encourage the cultivation of economic crops other than rice. Table 3 shows the hydrologic and WQ simulation results of SWAT according to the conservation of land use.  7%, 9.3%, and 11.0%, respectively. However, the SQ increased by 81.4% (Table 3). The change of the hydrology changed more severely in the area where high-density urbanization was concentrated. In 3009 mid-subbasin, the SCS-CN value of 2019 was 9.3 greater than the 1985 CN value and surface runoff of 2019 was 248.3% greater than 1985's. Due to the high-density urban area increase, the ET, SW, and GWQ reduced by 5.8%, 17.7%, and 25.9%, respectively. In 3013 mid-subbasins, where the rice paddy area was changed to upland crop area, the SQ increased; however, that was lower than the increase rate of the total watershed. SW and GWQ also decreased by 8.6% and 9.6%, respectively. The WQ simulation results of total area showed an increasing pattern over time. Between 1985 and 2019, the SS, T-P, and T-N increased by 13.5 mg/L, 0.009 mg/L, and 0.031 mg/L. However, 3009 urbanized area showed that the T-P and T-N concentration decreased by 2.7% and 7.4%, respectively. This phenomenon could be because N and P were not leaching to groundwater in high-density urban areas and were washed off to downstream of 3009 mid-subbasin. The 3013, which is an upland-cultivation-increased area, showed that the SS, T-P, and T-N concentration increased by 24.9 mg/L, 0.013 mg/L, and 0.063 mg/L, respectively. In particular, the rate of increase in WQ in 3013 uplandcrop-increase area was greater than the overall average of the GRB and the increase rate in urbanized areas. Figure 4 shows the monthly average changes in hydrology (SQ and GWQ) and WQ (SS, T-P, and T-N) according to urbanized (Figure 4a) and upland-cultivation-increased area (Figure 4b). In the figure, the grey area is the simulation result of GRB average, the yellow triangle line represents the 1985 results of each test subbasin, and the blue square line represents the 2019 results of test subbasins. In urbanized watershed, the increase in the overall impervious area increased the SQ, and the decreases were in GWQ. The SS concentration was shown to increase compared to the reference year (1985), which caused the increase in the SQ. The T-P and T-N concentration decreased compared to the reference year; in particular, the T-P decreased the most from June to September when the SQ increased (Figure 4a). In 3013 subbasin, which is an upland-cultivation-increased area, the increase is in SQ and WQ concentration. The T-P and T-N increased in spring and winter when the SQ decreased (Figure 4b).

Evaluation of WH Change
This study defined the result obtained from the land use of 1985 as the reference value under the most natural conditions and derived the hydrologic and WQ health for the years from 1995 to 2019 compared to that in 1985. Using the SWAT results, WH analyses of each component were implemented for 78 standard watersheds as the GRB. Figure 5 shows the four hydrologic and three WQ normalized health components for the GRB in the reference year. Hydrology subindices were evaluated by four hydrologic normalized components, and WQ subindex was calculated by SS, T-P, and T-N health components. An index value closer to 1 indicates a healthier watershed.

WH Evaluation Results for the Reference Year
Throughout the GRB, SQ was evaluated as healthy in upstream (Figure 5a), while SW and GWQ were evaluated as healthy in downstream of the GRB (Figure 5c,d). The T-P and T-N health index was evaluated as healthy in upstream (Figure 5g,h), and the WQ subindex was evaluated in the same pattern (Figure 5i). In the reference year, the WH of total watershed was evaluated between 0.17 and 0.80. The scores of ET, SQ, SW, and GWQ health evaluations before the high-density urban expansion in 3009 mid-subbasin were 0.67, 0.67, 0.38, and 0.38, respectively, with ET and SQ healthy above the overall average of GRB. In the 3013 mid-subbasin, where paddy fields account for 42.8%, the health evaluation results of ET and SQ were 0.42 and 0.31, respectively, lower than the overall average of GRB, while SW and GWQ were healthy at 0.65. WQ health for 3009 and 3013 mid-subbasins was assessed at 0.54 and 0.60, respectively. component were implemented for 78 standard watersheds as the GRB.
3.3.1. WH Evaluation Results for the Reference Year Figure 5 shows the four hydrologic and three WQ normalized health components for the GRB in the reference year. Hydrology subindices were evaluated by four hydrologic normalized components, and WQ subindex was calculated by SS, T-P, and T-N health components. An index value closer to 1 indicates a healthier watershed. Throughout the GRB, SQ was evaluated as healthy in upstream (Figure 5a), while SW and GWQ were evaluated as healthy in downstream of the GRB (Figure 5c,d). The T-P and T-N health index was evaluated as healthy in upstream (Figure 5g,h), and the WQ subindex was evaluated in the same pattern (Figure 5i). In the reference year, the WH of total watershed was evaluated between 0.17 and 0.80. The scores of ET, SQ, SW, and GWQ health evaluations before the high-density urban expansion in 3009 mid-subbasin were 0.67, 0.67, 0.38, and 0.38, respectively, with ET and SQ healthy above the overall average of GRB. In the 3013 mid-subbasin, where paddy fields account for 42.8%, the health evaluation results of ET and SQ were 0.42 and 0.31, respectively, lower than the overall average of GRB, while SW and GWQ were healthy at 0.65. WQ health for 3009 and 3013 midsubbasins was assessed at 0.54 and 0.60, respectively.

WH Change Based on Reference Year
Compared to the reference value, changes in hydrological and WQ health according to land-use changes from 1995 to 2019 were absolutely evaluated. Figure 6

WH Change Based on Reference Year
Compared to the reference value, changes in hydrological and WQ health according to land-use changes from 1995 to 2019 were absolutely evaluated. Figure 6 presents the boxplot of the evaluated health index for the 78 standard subbasins of GRB from 1985 to 2019. The total average of the 2019 SQ index evaluated value decreased to 0.17 as a result of the increase in impervious area. The SQ health index for the reference year was distributed between 0 and 1; this index from 1995 to 2019 was determined to range from 0 to 0.71 (Figure 6b), indicating that the health was worse than that in the reference year. The increase in impervious area caused a decrease in SW (Figure 6c) which resulted in the GWQ (Figure 6d) index tending to decrease every year. In 2019, the integrated hydrological health index decreased by an average of 0.18 compared to the reference year (Figure 6e).
Compared to that in the base year, the WQ health from 1995 to 2019 changed to within 0.1, while the range of changes in the hydrology index was analyzed very largely. The T-P health index for the reference year was distributed between 0 and 1, but the SS index from 1990 to 2019 was determined to range from 0 to 0.95 (Figure 6f), indicating that the health was worse than that in the reference year. The T-N and T-P also showed poor health compared to the reference year (Figure 6g,h). The integrated WQ health index for the three subindices (SS, T-P, and T-N) has been shown to deteriorate every year compared to the reference year. As a result of evaluating the health of the watershed by combining four hydrological health indicators and three WQ health indicators, WH has deteriorated since the reference year (Figure 6i). In particular, hydrological health tended to decrease compared to WQ (Figure 6j). Figure 7 presents the WH indices in two study subbasins (3009 and 3013) between 1985 and 2019, which were obtained based on the SWAT results for 1985. In Figure 7, the change in the value of each subindex and the changes compared with the reference year are shown on the basis of each average value of the mid-subbasin. When the WH changes were analyzed according to two test subbasins, the hydrologic WHIs decreased in the analysis periods compared to the reference value of 1985. In the urbanized subbasin (3009), the hydrology health index decreased to 0.26. In particular, the health of the SQ showed a deteriorating trend from 1995 to 2019. In comparison with the upland-cropcultivation-increased subbasin (3013), the SQ health worsened in 2019. Similar to the study, it was shown that the expansion in urban area and the decrease in forest area reduced SW, resulting in a decrease in GWQ [35][36][37].
Water 2021, 13, x FOR PEER REVIEW 13 of 19 of the increase in impervious area. The SQ health index for the reference year was distributed between 0 and 1; this index from 1995 to 2019 was determined to range from 0 to 0.71 (Figure 6b), indicating that the health was worse than that in the reference year. The increase in impervious area caused a decrease in SW (Figure 6c) which resulted in the GWQ (Figure 6d) index tending to decrease every year. In 2019, the integrated hydrological health index decreased by an average of 0.18 compared to the reference year (Figure 6e). Compared to that in the base year, the WQ health from 1995 to 2019 changed to within 0.1, while the range of changes in the hydrology index was analyzed very largely.
The T-P health index for the reference year was distributed between 0 and 1, but the SS index from 1990 to 2019 was determined to range from 0 to 0.95 (Figure 6f), indicating that the health was worse than that in the reference year. The T-N and T-P also showed poor health compared to the reference year (Figure 6g,h). The integrated WQ health index for the three subindices (SS, T-P, and T-N) has been shown to deteriorate every year compared to the reference year. As a result of evaluating the health of the watershed by combining four hydrological health indicators and three WQ health indicators, WH has deteriorated since the reference year (Figure 6i). In particular, hydrological health tended to decrease compared to WQ (Figure 6j). Figure 7 presents the WH indices in two study subbasins (3009 and 3013) between 1985 and 2019, which were obtained based on the SWAT results for 1985. In Figure 7, the change in the value of each subindex and the changes compared with the reference year are shown on the basis of each average value of the mid-subbasin. When the WH changes were analyzed according to two test subbasins, the hydrologic WHIs decreased in the analysis periods compared to the reference value of 1985. In the urbanized subbasin (3009), the hydrology health index decreased to 0.26. In particular, the health of the SQ showed a deteriorating trend from 1995 to 2019. In comparison with the upland-crop-cultivation-increased subbasin (3013), the SQ health worsened in 2019. Similar to the study, it was shown that the expansion in urban area and the decrease in forest area reduced SW, resulting in a decrease in GWQ [35][36][37]. The SS health index deteriorated rapidly in both study watersheds; however, the T-P and T-N health showed a conflicted trend. In urbanized subbasin (3009), T-P and T-N concentration was decreased from 1985 to 2019, and consequently the evaluated T-P and T-N health was healthier than the reference year. Meanwhile, the evaluated T-P and T-N health of upland-cultivation-increased watershed (3013) was slightly poorer than the reference year. Compared to the reference year, the distribution range of health indices The SS health index deteriorated rapidly in both study watersheds; however, the T-P and T-N health showed a conflicted trend. In urbanized subbasin (3009), T-P and T-N concentration was decreased from 1985 to 2019, and consequently the evaluated T-P and T-N health was healthier than the reference year. Meanwhile, the evaluated T-P and T-N health of upland-cultivation-increased watershed (3013) was slightly poorer than the reference year. Compared to the reference year, the distribution range of health indices changed. In the case of the urbanized subbasin, which was found to have undergone the most serious health degradation, the change in land use alone under the same weather conditions caused a significant change in hydrology and WQ health.

Land-Use Change Impact on WH
Land-use changes influence ecosystems by altering water cycles and WQ [38][39][40]. Many studies have confirmed that the land-use change alters regional hydrology [41][42][43], and WQ conditions and cities amplify ecosystem changes in natural watersheds [43][44][45] and urban sprawl impact on the agricultural area [46]. The results of the WH assessment show that its present land use is worse than its natural state in the past ( Figure 6). Although only satellite land-use data were analyzed from 1985 to 2019 under the same meteorological conditions, the hydrology and WQ simulation results for these years are greatly different. The surface runoff increased due to the increase in impervious area, but groundwater decreased, resulting in a change in the hydrology in the watershed (Figure 4). Our results in hydrological cycle and WQ corresponded with a considerable number of studies [47][48][49][50][51].
The newly revealed 12.3 % increase in 3009 mid-subbasin urban areas was shown to increase surface runoff in the 2019 SWAT simulation (Tables 2 and 3). Therefore, it affected the decrease in ET, soil water, and groundwater flow. The decreased ET corresponded to mid-subbasin, where the major land use changes from paddy to urban. This phenomenon can be explained as follows: as urbanization increased impervious areas, surface runoff increased and reduction ET, the ground inflow of rain, also decreased, which affected the groundwater level, and decreased the groundwater flow, while the runoff rate increased during flood seasons [52][53][54][55].
While urban areas increased, the rainfall runoff and pollutant emissions in residential and commercial areas showed various characteristics depending on the population density and economic activities in each area. However, as impervious areas, including parking lots and roads, increased in most urban areas, the rainfall runoff and flow and the WQ displayed rapid changes. In particular, because an urban area shows a remarkable initial runoff, by which pollutants are emitted at the beginning of rainfall events, a maximum outflow of pollutants occurs, which has a great impact on the downstream of urban areas [56]. However, in this study, it was found that the T-P and T-N health in the 3009 mid-subbasin, which was the most high-density urbanized area (Figure 2), improved compared to the reference year (Figure 7). The increase in the impermeable area caused a rapid increase in surface runoff and the first flush of water pollutants, which caused a decrease in infiltration and spread of groundwater depletion, preventing N and P from leaching into the soil. Similar findings are presented in some other studies. Hobbie et al. [57] confirmed that urban watersheds, characterized by high-density urbanized areas, lose phosphorus to surface runoff. Wang et al. [58] also found the effect of urbanization on stream WQ and high correlations according to the groundwater results, and that high-density urban growth was more efficient in reducing NO3-N and T-P concentration than sprawl development.
The 3013 mid-subbasin was a representative area where paddy rice fields were converted to upland cultivation area. In comparison to paddy rice fields, upland crops have different discharge characteristics due to the slope, scale, and shape of the field, as well as the characteristics of the rainfall. Because of the characteristics of the land used for paddy field, most of the upland reclamation is flat land with a slope close to 0%, so the characteristics of runoff are different from those of conventional farming on the sloping land [59]. In 3013 mid-subbasin, the amount of SQ increased while SW and GWQ decreased due to increased crop water demand (Figure 7). Compared to paddy fields, SQ health was poor due to the increase of SQ loss due to frequent water use for cultivating field crops. SS health s rapidly became poor due to increased SQ, and the T-P and T-N health score decreased due to frequent fertilizer use (Figure 7). The results of this phenomenon are analogous to those of previous studies about crop fertilization. Jang et al. [60] reported that the reclamation of upland areas has been expanded with the use of fertilizers to improve the productivity of crops. They have since had a negative impact on WQ issues in streams and groundwater. Zhang et al. [61] point out that the factors influencing T-N in crops and paddy fields differed. This was primarily due to variations in crops and paddy fields runoff, N fertilization rates, and management strategies. Ni et al. [62] published similar findings about SS, T-P, and T-N concentration increases in crop land as a result of improvements in crop management activities, including fertilization. Pollutant content in agricultural areas is a critical environmental concern due to the high risk of pollutants reaching nearby bodies of water by leaching or runoff [63].

Assessment of WH Using MND
Hydrologists and ecologists interested in encouraging sustainable activities are heavily involved in assessing watersheds in terms of human and ecosystem health [64]. The US EPA has moved toward integrated watershed assessment at the state level; they suggested that managers use integrated watershed evaluations to help them determine WH and building priority of watersheds for protecting and restoring [17]. Ahn and Kim [18,19] have proposed an application strategy to evaluate WH in the central area of South Korea. They assessed six indices of WHIs for 237 subbasins and found that the WH decreased further downstream in the watershed. In this study, we evaluated WH according to the consequence in land use and the findings showed that WH decreases in high-density urbanized area and upland-cultivation-increased area. These results are similar to previous studies for WH evaluation. Alilou et al. [65] used fuzzy logic theory to assess WH in the Khoy watershed and suggest a method for sub-watershed prioritization. To assess WH, they focused on climate vulnerability, relative surficial rock erosion rate, topographic indexes, thirteen morphometric characteristic features, and slope-weighted K-factor possible nonpoint source pollution. Human-caused activities and habitat changes had degraded the WH in the Khoy watershed (e.g., such as agriculture and expanding urbanization). Moreover, Mosaffaie et al. [22] suggest a framework for assessment for the WH during 2004-2018 in the Gorganroud watershed. In this framework, the key issues for WH are the lack of soil erosion rate, supplies of groundwater, and flood capacity. According to this study, the WH of the Gorganroud watershed deteriorates gradually as a result of socioeconomic activities such as agriculture and industry. These results have clearly shown that changes in land use such as urban expansion and an increase in crop area affect watershed health. They also imply that a long-term evaluation of health is necessary.
Using the proposed previous method to calculate distributed WHIs between 0 and 1 for every analysis, it was possible to compare only the WH of each sub-watershed for every analysis. This method is effective in comparing WH among sub-watersheds in the present state. However, when the change in the watershed environment takes place over a long period of time from the past to the future, ratings change at different times. In this case, the method could not provide any reference value. This study solved this problem by ensuring that changes in WH could be absolutely analyzed by using reference values for environmental changes at different times. Reference values could be provided for analyzing WH for each period or each evaluation period. In addition, an absolute comparison could now be possible for the same watershed. Thus, the environmental change in a watershed caused by human activity could be identified and compared to past natural conditions, and influential factors on WH could be analyzed.
We found that the changing hydrology and WQ of the basin could be assessed by the health index through the absolute comparison. The results of hydrological modeling, which utilizes land use as input, were used, and the results of the land use for the homogenous year were calculated as the reference year to derive the absolute value for health. Although only land-use data were analyzed under the same meteorological conditions using SWAT modelling, it can be extended to the analysis of WH for various indicators and data, such as climate change, socioeconomic change, vulnerabilities, and recovery priorities of the basin. The derived method from this study has the advantage of being able to intuitively judge how much the current WH has improved and deteriorated compared to the past.

Limitations of this Study
Although this study selected and analyzed only two out of six subindices evaluated by existing studies, it is also possible to analyze the change in WH by comparing reference values with land cover, vegetation, stream, habitat, aquatic ecology, organic matter, and heavy metal. Based on such an analysis, decision-making data for basin management could be produced. The bias and incompleteness of WH indices will cause uncertainty for this study. Moreover, because the modeling results for reference values and the data used for identifying the change in land use contain uncertainties, additional studies are required to produce reliable evaluation results. In particular, the national portal site provides observational hydrological and WQ data for 1990 and after. If a reference value is calculated by considering the period, more reliable results regarding hydrology and WQ health can be obtained. In addition, this study did not sufficiently verify the statistical adequacy of the probability distribution method that was adopted to calculate the reference value. For this reason, various statistical methods also need to be verified. If uncertainties are improved and nationwide hydrology and WQ models are constructed on a national scale, the long-term patterns of WH in terms of hydrology and WQ could be studied.

Conclusions
This study was aimed at establishing a reference condition for WHI based on MND and calculated WH change by land use. Among the previous study variables, two WH evaluation indices (hydrology and WQ) were selected. The urbanized area caused a rapid increase in surface runoff and the first flush of water pollutants, which caused a decrease in infiltration and spread of groundwater depletion, limiting N and P from leaching into the soil. In upland-crop-cultivation-increased area, the amount of SQ increased while SW and GWQ decreased due to increased crop water demand. This study indicated that the health of GRB deteriorates over time as a result of human activities (i.e., urbanization and upland reclamation). The findings of this study indicate that the derived approach has the advantage of providing an absolute comparison of WH change. Based on such an analysis, it can be extended to the analysis of WH for various indicators and data, such as climate change, socioeconomic change, vulnerabilities, and recovery priorities of the basin. To conclude, this study can provide the decision-making data for basin management.

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