Development and Application of an Interactive Coupling Rainfall-Runoff Model According to Soil Texture

: In this study, we developed a model that simulates surface and subsurface ﬂows for rainfall-runoff analyses using an interactive coupling method. To characterize the interaction between the surface and subsurface ﬂows, we studied the coupling analysis model. The proposed model was designed following comparisons with existing models. For the analysis of the surface and subsurface ﬂows, a governance equation was constructed. The goodness of ﬁt of the model was also tested. To examine the sensitivity of the input parameters, simulations were performed while changing the major parameters of the model according to the soil texture. The developed model showed high applicability to actual watersheds after adjustment of the parameters. This model can be applied to the extension module for channel analysis; therefore, it can be used efﬁciently in urbanized watersheds wherein the upstream and downstream parts are pervious and impervious watersheds, respectively.


Introduction
Rainfall-runoff models are mathematical models that physically represent runoff phenomena in watersheds. Such models are classified into lumped or distributed models, depending on whether the spatial variations in parameters are considered. A lumped model directly transforms rainfall to runoff without expressing the flow distribution in the watershed, such as that produced using the unit hydrograph method [1]. Semi-distributed models are variants of lumped models that have some characteristics of a distributed model. They consist of a set of clustered parameters applied in a quasi-spatially distributed manner. Semi-distributed models are classified according to their input. A model is considered semi-distributed if the input contains set and variance input parameters. Most models are semi-distributed because of data availability differences between the lumped and distributed models. The model process divides the watershed into smaller basins, each using different parameters. Sub-basins represent important features of the catchment and combine the advantages of batch and distributed models [2]. A distributed model can more effectively consider the spatiotemporal variations than a lumped model and can be used to analyze flood runoff relatively accurately, although, such models have some disadvantages including difficulties in composing input data and excessive calculation times. Additionally, one should analyze the applicability to virtual or actual watersheds when using the existing models because the effects of the input data composition on the sensitivity of parameters cannot be excluded [1].
There are several models for simulating surface and subsurface flows based on grids, but most models generally focus only on the flow of one layer, and the process for other layers constitutes the boundary conditions. For example, in the case of irrigation systems, only the surface flows and infiltration are analyzed because detailed analyses of groundwater are not important. Groundwater analysis models such as MODFLOW do not analyze surface flows. Some studies have attempted to conduct a simultaneous analysis of the parameters along with the finite element method in space and the finite difference method in time. The model uses an annular water level-flow relation curve equation at the subdistrict element to take into account the flow of negative currents in a steep slope [10,11].
In this study we aimed to develop a model that can interpret the interaction of surface and subsurface flows for rainfall-runoff analyses by types of coupling. First, rainfallrunoff models were searched for and compared, and each model to which the coupling type can be applied was compared to determine the direction of model development. Thereafter, governing equations for surface and subsurface flow analysis were constructed, and potential models were developed. During model development, rainfall and runoff data from a measured watershed were used to compare the estimated runoff of the models and the measured runoff. The comparison of the calculated and measured runoff was evaluated based on four goodness of fit measures: the root mean squared error (RMSE), percent bias (PBIAS), Nash-Sutcliff efficiency (NSE), and persistence model efficiency (PME). To verify the improvement in the applicability of the chosen model, a more refined fitness criterion was proposed, and the improvement in the applicability of the model to input parameter changes was tested.

Coupling System
The coupling system of surface and subsurface flows must include a numerical process for the surface and subsurface flows and a numerical process for the internal and external conditions. The surface flow system consists of a partial differential equation in the form of a hyperbola or parabola, whereas the subsurface flow system consists of a partial differential equation in the form of a parabola. These two types of coupling processes for the surface and subsurface flow systems can be classified into the following four types, as shown in Figure 1: uncoupled (UC), degenerated uncoupled (DU), iterative coupling (IC), and fully coupled (FC) [12].
Water 2021, 13, x FOR PEER REVIEW 3 of 20 The Vflo™ model is a physically based distribution model developed in 2002 by Vieux of the University of Oklahoma, a world-renowned authority on flood forecasting using rainfall radar. It is capable of both heavy-duty and continuous outflow simulation. The Vflo™ model developed by Vieux is a short-term runoff model that uses distributed parameters along with the finite element method in space and the finite difference method in time. The model uses an annular water level-flow relation curve equation at the subdistrict element to take into account the flow of negative currents in a steep slope [10,11].
In this study we aimed to develop a model that can interpret the interaction of surface and subsurface flows for rainfall-runoff analyses by types of coupling. First, rainfall-runoff models were searched for and compared, and each model to which the coupling type can be applied was compared to determine the direction of model development. Thereafter, governing equations for surface and subsurface flow analysis were constructed, and potential models were developed. During model development, rainfall and runoff data from a measured watershed were used to compare the estimated runoff of the models and the measured runoff. The comparison of the calculated and measured runoff was evaluated based on four goodness of fit measures: the root mean squared error (RMSE), percent bias (PBIAS), Nash-Sutcliff efficiency (NSE), and persistence model efficiency (PME). To verify the improvement in the applicability of the chosen model, a more refined fitness criterion was proposed, and the improvement in the applicability of the model to input parameter changes was tested.

Coupling System
The coupling system of surface and subsurface flows must include a numerical process for the surface and subsurface flows and a numerical process for the internal and external conditions. The surface flow system consists of a partial differential equation in the form of a hyperbola or parabola, whereas the subsurface flow system consists of a partial differential equation in the form of a parabola. These two types of coupling processes for the surface and subsurface flow systems can be classified into the following four types, as shown in Figure 1: uncoupled (UC), degenerated uncoupled (DU), iterative coupling (IC), and fully coupled (FC) [12]. In general, the coupling system of the surface and subsurface flows includes mathematical descriptions of the surface, subsurface, external boundary conditions, and internal boundary conditions. The surface flow systems are often described as hyperbolic or parabolic differential equations, and the subsurface flow systems are described as parabolic partial differential equations (PDEs). Theoretically, the higher the level of coupling, the higher the accuracy. However, this is not necessarily correct because the results can differ numerically depending on the degree of analysis. The following presents descriptions of the four levels of coupling. In general, the coupling system of the surface and subsurface flows includes mathematical descriptions of the surface, subsurface, external boundary conditions, and internal boundary conditions. The surface flow systems are often described as hyperbolic or parabolic differential equations, and the subsurface flow systems are described as parabolic partial differential equations (PDEs). Theoretically, the higher the level of coupling, the higher the accuracy. However, this is not necessarily correct because the results can differ numerically depending on the degree of analysis. The following presents descriptions of the four levels of coupling.
The first is the uncoupled condition. This means that each system operates independently at each stage, and the surface flow is calculated according to the mechanical relationship. When the internal boundary condition is determined based on the solution, the rest of the system is calculated. There is no mutual feedback because the surface flow and subsurface flow are independent. However, because the boundary conditions are applied to the two flow systems, an approximation is required for mathematical analysis. Therefore, it is convenient to estimate the boundary condition from the surface flow system and use it as the condition for the previous time step. The first coupling level can be divided into two subgroups. The first subgroup is described above, and the second subgroup involves mimicking either the surface or subsurface flow system in an algebraic formulation. This is the approach most often used in indicator irrigation models and is referred to as DU.
The third level of coupling involves feedback between the two flow systems. In the first step, similar to the uncoupled level, one flow system is calculated, and when the interfacial boundary conditions are formed, the remaining flow systems are solved using these conditions. The solution of the second flow system is used to update the internal boundary conditions within the same time step. The first flow system is resolved through the updated internal boundary conditions, and the entire process is repeated until the convergence criterion is achieved. This coupling level is referred to as alternating iterative or IC.
The fourth coupling level, the most complete form, solves both flow systems and internal boundary conditions together. That is, two PDEs and an internal equation (generally a differential equation) are solved simultaneously. This coupling level is called FC. Notably, it is not clear whether the FC approach is superior to the IC approach. FC involves the separation and solving of parabolic and hyperbolic equations in a single set. This is because the characteristics of different equations can lead to difficulties in the numerical analysis. In addition, the FC approach means that a larger system has to be solved. In general, numerical solutions involve iterative interpretations of spatially and temporally discretized equations, from a computational point of view, and the most basic difference arises between external (for IC) and internal (for FC) iterations.

Comparison among Models
In this study, well-known lumped models and semi-distributed models that only analyze infiltration by separating the effective rainfall without conducting direct simulations for the subsurface flow, were excluded from our comparison of analysis methods. In total, we compared eight models used for distributed rainfall-runoff analyses based on grids by conducting research on the analysis of the surface flow and the infiltration and subsurface flow.
As shown in Table 1, there are various analysis methods for surface flow and for infiltration and subsurface flow according to the model research. Most models use empirical or theoretical formulas as originally presented or improved forms achieved through modifications.
The AGNPS model is useful for simulating the spatial distribution of watershed properties, providing immediate responses compared to experiments, and understanding erosion processes. However. it requires extensive preparatory data and is labor and time intensive. It requires high levels of programming ability, and cannot evaluate in-stream processes. In addition, a number of SCS runoff curves and other empirical relationships are used to calculate the runoff and maximum flow, and flood waves cannot be represented. Most models using physically based flow management equations use approximate numerical solutions of equations that are affected by numerical instability problems and limited space and time, growth, and watershed size [13]. The ANSWERS model uses a storage-based flow equation and cannot represent flood waves, although the program code is easy to modify, and the suitability for the analysis results is high. It is difficult to verify as it requires detailed input data preparation, is computationally intensive, and the results are sensitive to minute changes in input variables [13]. The DWSM model is a physically based flow management equation that can represent flood waves arising from a single rainfall event. In addition, DWSM is advantageous when applied to large watersheds because it uses approximate analysis solutions and is not limited by spatial and temporal increments [14].
The GSSHA model has the ability to compute the remaining water on the surface, which cannot be estimated when using the conceptual models; the latter leading to poor estimation of the observed hydrograph. Also, the flood velocity and flood depths at any location inside the catchment can be easily determined. GSSHA model results provided detailed information about residential areas and roads that were inundated during flood events [15].
The Hydro-BEAM model should be constructed using four categories of parameters: scale, topographic, shape, and drain network, using separate GIS tools. When the rainfall pattern on the drainage basin is uniform, the hydrograph shape has a disadvantage that it can only reflect the influence of specific watershed characteristics [16].
MIKE SHE can be dynamically coupled to MIKE 11, a one-dimensional surface water model that simulates fully dynamic channel flow, to represent river processes and river management. It allows process-based approaches to apply different model structures within the same modeling framework [17]. Vflo TM can run calibrated real-time model engines and generate real-time flow predictions almost anywhere in the basin. However, depending on the size of the study area and the grid cell size, extensive computer power may be required [17].
As the model under development here covers surface and subsurface flows and the effects of the external conditions of rainfall and inflow rate, it corresponds to the IC process in Figure 1. Thus, it must include a process to analyze the numerical solutions of the surface and subsurface flow systems, as well as a numerical analysis for the external boundary conditions as feedback.
For the subsurface runoff analysis, we used the Boussinesq equation based on Darcy's Law and the Dupuit-Forchheimer assumption. Regarding the dynamic flow of surface runoff, we configured a 2D diffusion wave model that simplified the Saint-Venant equation into an approximate solution for the diffusion wave. Furthermore, we analyzed the interaction between the surface and subsurface flows by the conservation of mass, and we conducted numerical simulations using the finite volume method for space and the Crank-Nicolson method for time.
The advantages and disadvantages of these models have been demonstrated. The difference between the eight models is due to various functions, and the important points for model development are as follows: (1) ease of constructing input data (rainfall, runoff, and soil and land uses); (2) formulation of governing equations (including dimensions); (3) interfacial boundary conditions that enforce at least pressure and mass flux continuity at surface/underground interfaces; and (4) numerical approaches for spatial and temporal discretization and coupling.
Our development model is highly responsive to soil types and properties, and by utilizing the water retention properties of soil texture provided by USDA, porosity, residual water content, and supply distribution can be constructed relatively easily. The model can address feedback and interactions between surface and subsurface flows. The models tested here are mass-conserving as they ensure a perfect balance of water between the surface and subsurface systems. We provide a solution that repeatedly calculates the interaction of surface and subsurface flow through the mass conservation framework until the convergence criterion is achieved, while updating the internal boundary condition according to the subsurface saturation condition. In addition, the structured grid method is adopted for efficient parallelization and provides advantages in computational simplicity.
Maxwell et al. [18] presented the results of a cross-comparison study of seven coupled surface-subsurface models. Their simulation results showed good agreement for the simulator test cases, whereas more complicated test cases exhibited some differences in physical process representation and numerical solution approaches between the models. Models with more traditional runoff generation mechanisms such as excess infiltration and saturation demonstrated more consensus among models, while models with heterogeneity and complex water table dynamics emphasized differences in formulation. In general, all the models exhibited the same qualitative behavior, building confidence in hydrological uses.
Rather than finding better solutions through comparison between models, like Maxwell et al. [18] did, it was deemed necessary to determine the infiltration amount inherent in the effective rainfall separation or outflow analysis model, which is the scope of analysis for domestic and foreign models. The purpose and orientation of the model development were established by applying the IC coupling system, and the two outflow processes were combined according to soil saturation conditions. This method considers surface and subsurface runoff separately through continuous analysis of infiltration from the time of a rainfall event based on the constructed terrain.

Developed Model
The developed model has three basic assumptions described below. The conceptual definitions are in Figure 2.

1.
When rainfall occurs on a pervious watershed, surface or subsurface runoff occurs depending on the degree of saturation of the soil layer, and interactions between the surface and subsurface runoff occur ( Figure 2A).

2.
When rainfall occurs on an impervious watershed, only surface runoff occurs ( Figure 2B).

3.
The base flow is not considered.

Developed Model
The developed model has three basic assumptions described below. The conceptual definitions are in Figure 2. 1. When rainfall occurs on a pervious watershed, surface or subsurface runoff occurs depending on the degree of saturation of the soil layer, and interactions between the surface and subsurface runoff occur ( Figure 2 (A)). 2. When rainfall occurs on an impervious watershed, only surface runoff occurs ( Figure 2 (B)). 3. The base flow is not considered. According to the three basic assumptions and conceptual definitions above, the numerical modeling concept for runoff analysis can be subdivided as follows ( Figure 3): 1. Uniform rainfall R(t) occurs in an area with a certain area size. 2. The vertical distance from the reference plane to where the pervious soil layer begins is Z (x, y). 3. The saturation depth for the thickness of the pervious soil layer is expressed as hp (x, y). 4. The surface runoff depth when the thickness of the pervious soil layer and the saturation depth hp (x, y) become equal or the rainfall exceeds the infiltration capacity of the pervious soil layer is hs (x, y). 5. The runoff rate per unit width is q (x, y), surface slope is S (x, y), average flow rate is v (x, y), and total head from the reference plane is H (x, y). According to the three basic assumptions and conceptual definitions above, the numerical modeling concept for runoff analysis can be subdivided as follows ( Figure 3):

1.
Uniform rainfall R(t) occurs in an area with a certain area size.

2.
The vertical distance from the reference plane to where the pervious soil layer begins is Z (x, y).

3.
The saturation depth for the thickness of the pervious soil layer is expressed as hp (x, y).

4.
The surface runoff depth when the thickness of the pervious soil layer and the saturation depth hp (x, y) become equal or the rainfall exceeds the infiltration capacity of the pervious soil layer is hs (x, y). 5.
The runoff rate per unit width is q (x, y), surface slope is S (x, y), average flow rate is v (x, y), and total head from the reference plane is H (x, y). Water 2021, 13, x FOR PEER REVIEW 8 of 20 The dynamic flow of the surface runoff can be represented by a 2D diffusion wave model that has simplified the Saint-Venant equation, consisting of a continuous equation and momentum equation, into a diffusion wave approximate solution [19], that consists of the following expressions: Friction slope term: where, H is the total head, S is the slope in the x-direction, and S is the slope in the y-direction. Surface flow partial differential equation: where, h is the surface runoff depth. When surface flow occurs, and h become zero, and the 2D diffusion wave model can be expressed with the above equations as follows: where, R(t) is uniform rainfall, The subsurface runoff was constructed using the Boussinesq equation based on Darcy's law and the Dupuit-Forchheimer assumption, with the continuous equation and momentum equation as basic equations. The process is described below.
Assuming that there is no resistance to the vertical flow, the effects of the movement of the free surface and rainfall can be separated into q and R, respectively. The movement of the free surface in a pervious soil layer over time is given as q = n , and the rainfall can be expressed as R = R(t). Under the assumption of horizontal isotropy for inertial flow and through Darcy's law, the movement can be expressed with the Boussinesq equation, that represents a 2D nonuniform flow in an unconfined pervious layer for a constant flow rate [20]: where, n is the effective porosity, h is the subsurface runoff depth, K is the hydraulic conductivity and H is the total head. In the developed model, eight main parameters in ASCII file format are used as inputs, and these parameters consist of a digital elevation model (DEM), permeability coefficient, porosity, effective soil depth, roughness coefficient, pervious/impervious condition, residual water content, and pore distribution. Each parameter is differentiated by the soil type, land cover type, and drainage rating. The dynamic flow of the surface runoff can be represented by a 2D diffusion wave model that has simplified the Saint-Venant equation, consisting of a continuous equation and momentum equation, into a diffusion wave approximate solution [19], that consists of the following expressions: Friction slope term: where, H is the total head, S f x is the slope in the x-direction, and S f y is the slope in the y-direction. Surface flow partial differential equation: where, h s is the surface runoff depth. When surface flow occurs, dz dt and ∂ ∂t h p become zero, and the 2D diffusion wave model can be expressed with the above equations as follows:

∂H ∂t
where, R(t) is uniform rainfall, D(h sur ) = 1 n h 5/3 s √ S f . The subsurface runoff was constructed using the Boussinesq equation based on Darcy's law and the Dupuit-Forchheimer assumption, with the continuous equation and momentum equation as basic equations. The process is described below.
Assuming that there is no resistance to the vertical flow, the effects of the movement of the free surface and rainfall can be separated into q h por and R, respectively. The movement of the free surface in a pervious soil layer over time is given as q h dx = n e ∂h por ∂t , and the rainfall can be expressed as R = R(t). Under the assumption of horizontal isotropy for inertial flow and through Darcy's law, the movement can be expressed with the Boussinesq equation, that represents a 2D nonuniform flow in an unconfined pervious layer for a constant flow rate [20]: where, n e is the effective porosity, h p is the subsurface runoff depth, K is the hydraulic conductivity and H is the total head. In the developed model, eight main parameters in ASCII file format are used as inputs, and these parameters consist of a digital elevation model (DEM), permeability coefficient, porosity, effective soil depth, roughness coefficient, pervious/impervious condition, resid-ual water content, and pore distribution. Each parameter is differentiated by the soil type, land cover type, and drainage rating.
As mentioned above, the developed model can be used to analyze under the basic assumption that in the event of rainfall in a pervious watershed, surface or subsurface runoff will occur depending on the degree of saturation of the soil layer, and interactions between a pervious watershed and rainfall and between the surface runoff and subsurface runoff occur. Thus, the model reacts to the soil type and characteristics. Among the eight main input parameters, we configured the porosity, residual water content, and pore distribution using the water retention properties of soil textures provided by the USDA (Table A1). For the residual water content among the input data, the minimum value was determined by subtracting the standard deviation from the mean, which resulted in a negative value. Hence, we used zero as the minimum value so that it could be recognized in the model.

Applied Watershed and Rainfall
To examine the goodness of fit of the developed model, we selected and analyzed a real watershed, Guryang watershed (Donghyang-myeon, Jinan-gun, Jeolabuk-do, Korea) located at a latitude of 35 • 46 03 -35 • 53 35 N and a longitude of 127 • 30 24 -127 • 44 47 E; K-water and the hydrological data, such as rainfall, runoff, and evapotranspiration, have been recorded at the watershed. All measurements were subjected to quality control by verifying their range according to climatological values, as well as spatiotemporal consistency. The Guryang watershed is a typical steep mountain watershed with an area of 172.37 km 2 and a large pervious area that accounts for 97% of the total watershed area ( Figure 4).
Water 2021, 13, x FOR PEER REVIEW 9 of 20 As mentioned above, the developed model can be used to analyze under the basic assumption that in the event of rainfall in a pervious watershed, surface or subsurface runoff will occur depending on the degree of saturation of the soil layer, and interactions between a pervious watershed and rainfall and between the surface runoff and subsurface runoff occur. Thus, the model reacts to the soil type and characteristics. Among the eight main input parameters, we configured the porosity, residual water content, and pore distribution using the water retention properties of soil textures provided by the USDA (Table A1). For the residual water content among the input data, the minimum value was determined by subtracting the standard deviation from the mean, which resulted in a negative value. Hence, we used zero as the minimum value so that it could be recognized in the model.

Applied Watershed and Rainfall
To examine the goodness of fit of the developed model, we selected and analyzed a real watershed, Guryang watershed (Donghyang-myeon, Jinan-gun, Jeolabuk-do, Korea) located at a latitude of 35°46′03″-35°53′35″ N and a longitude of 127°30′24″-127°44′47″ E; K-water and the hydrological data, such as rainfall, runoff, and evapotranspiration, have been recorded at the watershed. All measurements were subjected to quality control by verifying their range according to climatological values, as well as spatiotemporal consistency. The Guryang watershed is a typical steep mountain watershed with an area of 172.37 km 2 and a large pervious area that accounts for 97% of the total watershed area ( Figure 4). Among the input data, a DEM ( Figure 5) was constructed from numerical data from the National Geographic Information Institute. The effective soil depth (Figure 6), roughness coefficient (Figure 7), and impervious areas (Figure 8) are obtained from land use and detailed soil maps provided by the Ministry of Environment and the National Institute of Agricultural Sciences. In addition, permeability (Figure 9), porosity ( Figure 10), residual moisture content (Figure 11), and pore size distributions ( Figure 12) were set to indicate moisture-preserving properties for various soil textures provided by the USDA.                                  (Table 2).

Goodness of Fit
The goodness of fit determination method for each method is as follows. (1) In the case of the RMSE, the goodness of fit was determined to be high if the standard deviation of the predicted error of the model was close to zero. (2) The PBIAS is an estimation of the average trend of the difference between the model's calculated runoff hydrograph and the measured runoff hydrograph. This metric is used to evaluate model performance with respect to the tendency to exceed or underestimate the simulation flow. The optimal value of the goodness of fit for this measure is 0.0. A positive value means that the calculated runoff hydrograph was underestimated compared with the measured runoff hydrograph; a negative value means that it was overestimated, and the goodness of fit was low. If the PBIAS value is greater than 20%, it is considered an unacceptable value. (3) The NSE is frequently used to determine the relative size of the residual variance with respect to the measured data variance. The NSE value ranges from −1 to 1. The closer to 1, the more accurate the model is. If the NSE value is 0, it means that the model prediction is as accurate as the average value of the observed data. A value <0 indicates that the average value of the observed data is a better predictor than the model results. (4) In the case of the PME, a value of 1.0 means that the calculated runoff hydrograph and the measured runoff hydrograph were the same. If the value is higher than 0.0, the model is considered to display a goodness of fit [21][22][23].
where q sim t is the flow of the calculated runoff hydrograph at time t; q obs t is the flow of the measured runoff hydrograph at time t; and q mean is the average flow of the measured runoff hydrograph.

Comparison between the Measured Flows and Runoff Analysis Results from the Model
When the measured runoff hydrographs of the four rainfall events in the Guryang watershed were compared with those calculated using the developed model, the shapes of the hydrographs for the four rainfall events were similar. However, the calculated peak flows of the simulated runoff hydrographs were smaller (range: 1.12-9.84 m 3 /s) than those of the measured runoff hydrographs by 11% on average (Table 3, Figure 13). Table 3. Comparison of peak flows (m 3 /s) of the runoff hydrographs for the initial input data.

Comparison between the Measured Flows and Runoff Analysis Results from the Model
When the measured runoff hydrographs of the four rainfall events in the Guryang watershed were compared with those calculated using the developed model, the shapes of the hydrographs for the four rainfall events were similar. However, the calculated peak flows of the simulated runoff hydrographs were smaller (range: 1.12-9.84 m 3 /s) than those of the measured runoff hydrographs by 11% on average (Table 3, Figure 13).

Goodness of Fit Results
Dongfang and Defu [24] proposed a coupling model that fully integrates surface and subsurface flows using 2D kinematics and the Richard equation, and this model was verified through a sandbox experiment and soil column tests. Krzhizhanovskaya et al. [25] recommended the use of a coupling-based distributed model of surface and subsurface flows as a decision support system for issuing early flood warnings and disaster management. In addition, they provided modeling results for a flood in St. Petersburg, Russia, that were used for decision making. Sinha et al. [26] developed the HYPE model, which is a semi-distributed hydrological model, for verifying the spatiotemporal changes in the runoff coefficient and applied it to the Selke region in the central part of Germany. The results showed that the precipitation upstream of Silberhutte and the runoff coefficient continuously increased in response to topographical effects, whereas evaporation decreased.
Additionally, He et al. [27] applied the MIKE-SHE hydrologic model together with the ensemble transform Kalman filter data assimilation technique to improve the performance of surface and subsurface water simulation models within a small urban watershed in Denmark. They revealed that although the quality of rainfall data is important in Figure 13. Comparison of runoff hydrographs.

Goodness of Fit Results
Dongfang and Defu [24] proposed a coupling model that fully integrates surface and subsurface flows using 2D kinematics and the Richard equation, and this model was verified through a sandbox experiment and soil column tests. Krzhizhanovskaya et al. [25] recommended the use of a coupling-based distributed model of surface and subsurface flows as a decision support system for issuing early flood warnings and disaster management. In addition, they provided modeling results for a flood in St. Petersburg, Russia, that were used for decision making. Sinha et al. [26] developed the HYPE model, which is a semi-distributed hydrological model, for verifying the spatiotemporal changes in the runoff coefficient and applied it to the Selke region in the central part of Germany. The results showed that the precipitation upstream of Silberhutte and the runoff coefficient continuously increased in response to topographical effects, whereas evaporation decreased.
Additionally, He et al. [27] applied the MIKE-SHE hydrologic model together with the ensemble transform Kalman filter data assimilation technique to improve the performance of surface and subsurface water simulation models within a small urban watershed in Denmark. They revealed that although the quality of rainfall data is important in rainfall predictions, the initial conditions are more important when predicting the groundwater head. Visser et al. [28] proposed a framework for evaluating the effects of climate change on rivers using coupled hydrological and hydroecological models. They applied the SRES A1F1 climate change scenario to the River Nar, Norfolk, UK and found that a reduction in diversity of in-stream organisms is occurring.
Lastly, Guzha and Hardy [29] developed a coupling model that combined MODFLOW and TOPNET, and applied it to the Big Darby Waters region in Ohio, U.S. For calibration and verification purposes, they analyzed the average goodness of fit, PBIAS, and NSE between the measured flows and simulation results. The results obtained were 0.83, 11.15%, and 0.83, respectively, which indicated good agreement. By calibrating the model, they were able to adjust the average goodness of fit to 0.84, the PBIAS to 8.75%, and the NSE to 0.85. Because the coupling analysis method is based on numerical analyses, the goodness of fit for the calculated results may vary depending on the application of the model. Therefore, it may be necessary to assess the performance of the developed model for specific applications. Based on the results of the goodness of fit for the model developed in this study, we investigated the possibility of achieving further improvements in the applicability of the model through parameter adjustment.
To verify the applicability of a developed model, the goodness of fit must be determined between the measured rainfall-runoff data for a watershed and the calculation results derived from the developed model by using characteristic data for the watershed. For the goodness of fit determination in this study, we used four measures that are appropriate for quantitative comparisons of the measured runoff hydrographs and the calculated runoff hydrographs of the model, namely, the RMSE, PBIAS, NSE, and PME.
The results showed that the RMSE values were higher than 5.0 for rainfall events 1 and 4, indicative of a lack of fit. The PBIAS values were higher than 0.1 for all the rainfall events, again indicating a lack of fit. The NSE values were close to 1.0 for all the rainfall events, and thus, an unfit judgment was not made based on these results. The PME values were lower than zero for rainfall events 2 and 4, again indicative of a lack of fit. In general, the goodness of fit results indicated that improvements would be required for all rainfall events (Table 4).

Sensitivity Analysis of the Input Parameters and Parameter Calibration
In the developed model, the initial conditions are classified by land use, land type, and drainage class. Because the setting of the soil moisture property value for the 11 soil textures was difficult, the water-retention properties based on soil texture were reset for the six drainage classes based on the sensitivity analysis to improve the results and applicability of the model.
To adjust for the tendency to underestimate the peak flow of the calculated runoff hydrograph compared with the measured runoff hydrograph and to improve the goodness of fit of the model, the sensitivity of the developed model was examined by changing the initial conditions. Among the major input parameters for the developed model, the porosity, residual water content, and pore distribution were selected as the three parameters for the sensitivity analysis. For each parameter, the averages of the USDA water-retention properties by soil texture were set as the initial input data, and the input data were divided into six blocks of 25% in the range of ±75%. When the residual water content was negative in the variation range, we arbitrarily adjusted it to zero to prevent numerical errors in the model.
The sensitivity analysis of the porosity, residual water content, and pore distribution for the four rainfall events showed that the porosity and pore distribution were inversely related to the peak flow of the runoff hydrograph. The residual water content showed a proportional relationship with the peak flow of the runoff hydrograph. Specifically, the residual water content showed a relatively small sensitivity when it was negative based on the initial condition, but the average sensitivity was higher than 15% when it was positive, and the maximum sensitivity was 54.21%. Thus, larger residual water contents corresponded to a higher sensitivity in the peak flow for the runoff hydrograph ( Figure 14). These findings suggest the water-retention properties based on soil texture for the developed model should be derived from normal rainfall-runoff analyses when possible.
Water 2021, 13, x FOR PEER REVIEW 15 of 20 was negative in the variation range, we arbitrarily adjusted it to zero to prevent numerical errors in the model. The sensitivity analysis of the porosity, residual water content, and pore distribution for the four rainfall events showed that the porosity and pore distribution were inversely related to the peak flow of the runoff hydrograph. The residual water content showed a proportional relationship with the peak flow of the runoff hydrograph. Specifically, the residual water content showed a relatively small sensitivity when it was negative based on the initial condition, but the average sensitivity was higher than 15% when it was positive, and the maximum sensitivity was 54.21%. Thus, larger residual water contents corresponded to a higher sensitivity in the peak flow for the runoff hydrograph ( Figure 14). These findings suggest the water-retention properties based on soil texture for the developed model should be derived from normal rainfall-runoff analyses when possible. In the case of the Guryang watershed, the soil drainage classes mostly consist of the slightly good and good classes. Therefore, the appropriate porosity and pore distribution is +25% in the slightly good class and +50% in the good class, and the appropriate residual water content is −50% in the good class and −25% in the slightly good class.
Based on the sensitivity analysis results that considered the measured runoff hydrographs for the four rainfall events of the Guryang watershed and the water-retention properties by soil texture, smaller porosity and pore distribution values will result in larger amounts of runoff. Furthermore, a larger residual water content will result in a larger value for the calculated runoff. Therefore, the proportional and inversely proportional relationships of the residual water content, porosity, and pore distribution with the peak flow can be applied for making adjustments. In the case of porosity and pore distribution, the input value and peak flow showed an inversely proportional relationship, with the average value assumed as the average drainage class. The parameters were calibrated by applying −75% as very poor, −50% as poor, −25% as slightly poor, +25% as slightly good, +50% as good, and +75% as very good. In the case of the residual water content, the input value and peak flow had a proportional relationship. Thus, with the average value as the average drainage class, the parameters were calibrated by applying In the case of the Guryang watershed, the soil drainage classes mostly consist of the slightly good and good classes. Therefore, the appropriate porosity and pore distribution is +25% in the slightly good class and +50% in the good class, and the appropriate residual water content is −50% in the good class and −25% in the slightly good class.
Based on the sensitivity analysis results that considered the measured runoff hydrographs for the four rainfall events of the Guryang watershed and the water-retention properties by soil texture, smaller porosity and pore distribution values will result in larger amounts of runoff. Furthermore, a larger residual water content will result in a larger value for the calculated runoff. Therefore, the proportional and inversely proportional relationships of the residual water content, porosity, and pore distribution with the peak flow can be applied for making adjustments. In the case of porosity and pore distribution, the input value and peak flow showed an inversely proportional relationship, with the average value assumed as the average drainage class. The parameters were calibrated by applying −75% as very poor, −50% as poor, −25% as slightly poor, +25% as slightly good, +50% as good, and +75% as very good. In the case of the residual water content, the input value and peak flow had a proportional relationship. Thus, with the average value as the average drainage class, the parameters were calibrated by applying −75% as very good, −50% as good, −25% as slightly good, +25% as slightly poor, +50% as poor, and +75% as very poor (Table A2).

Validation of the Applicability to Real Watersheds
When the runoff hydrographs for the measured rainfall events of the Guryang watershed were compared with the parameter-calibrated hydrographs derived by applying the porosity, residual water content, and pore distribution values that were calibrated by soil drainage class, the results showed improvements compared with those before parameter calibration ( Figure 15).

Validation of the Applicability to Real Watersheds
When the runoff hydrographs for the measured rainfall events of the Guryang tershed were compared with the parameter-calibrated hydrographs derived by apply the porosity, residual water content, and pore distribution values that were calibrated soil drainage class, the results showed improvements compared with those before par eter calibration (Figure 15). When the peak flow of the measured runoff hydrograph was compared with the p flow of the hydrograph after parameter calibration, the calculated peak flow was sma by 1% on an average (−0.66 to +3.64 m 3 /s). The calibrated parameters improved the d by 10% on an average in all rainfall events compared to the data before parameter cali tion ( Table 5). As shown in the comparison results for the runoff hydrograph type peak flow, the applicability of the model improved after parameter calibration accord to the six drainage classes.  When the peak flow of the measured runoff hydrograph was compared with the peak flow of the hydrograph after parameter calibration, the calculated peak flow was smaller by 1% on an average (−0.66 to +3.64 m 3 /s). The calibrated parameters improved the data by 10% on an average in all rainfall events compared to the data before parameter calibration ( Table 5). As shown in the comparison results for the runoff hydrograph type and peak flow, the applicability of the model improved after parameter calibration according to the six drainage classes. When the four types of quantitative goodness of fit measures were evaluated for the runoff hydrographs after parameter calibration and the measured runoff hydrographs, the data remained unfit. However, the RMSE (43%), PBIAS (34%), NSE (2%), and PME (46%) improved after parameter calibration based on the goodness of fit measures. Consequently, the applicability of the model improved after parameter calibration (Table 6).

Conclusions
In this study, methods for rainfall-runoff analyses were reviewed, and a model with interactive coupling was developed to analyze both surface and subsurface flow. The developed model was then applied to a real watershed. Four goodness of fit tests were conducted to evaluate the applicability of the model. The main results of this study can be summarized as follows:

1.
With the developed model, which is a grid-based distributed rainfall-runoff model, the measured and simulated runoff hydrographs were compared for four rainfall events using the initial input data for the runoff parameters including the porosity, residual water content, and pore distribution. The results showed that the peak flows of the simulated runoff hydrographs were 11% lower on an average than those of the measured runoff hydrographs.

2.
The sensitivities of the runoff parameters, namely, the porosity, residual water content, and pore distribution, were analyzed in the range of ±75% based on the average USDA water-retention properties categorized by soil texture for each parameter. The results showed that the porosity and pore distribution were inversely proportional to the peak flow, whereas the residual water content is directly proportional to the peak flow. The sensitivity of the residual water content for the peak flow was higher than 15% on an average, and it reached a maximum of 54.21%. Thus, among the runoff parameters, the residual water content had the largest sensitivity.

3.
Based on the sensitivity analysis results, the appropriate criteria for the water-retention properties by soil texture were calibrated in accordance with the watershed characteristics of South Korea for six drainage classes (very poor, poor, slightly poor, slightly good, good, and very good) within the Guryang watershed.

4.
When the runoff analysis was performed again with the calibrated parameters, the peak flow was smaller by 1% on an average than the results for the initial input data, and these data showed a higher goodness of fit with the peak flows of the measured runoff hydrographs.

5.
The analysis results of the quantitative goodness of fit of the model showed that the RMSE, PBIAS, NSE, and PME improved by 43%, 34%, 2%, and 46%, respectively, for all four rainfall events.
The results indicate that rainfall-runoff analyses will be improved when the waterretention properties by soil texture are calibrated and applied to the developed model, and applicability to real watersheds will also be high. The applicability of the proposed model as a distributed model could be improved in the future if a support feature to allow the construction of spatial information in the model and an auto evaluation feature achieved through sensitivity analyses of the simulation results are added.
Author Contributions: C.-J.K.: investigation; conceptualization; data curation; formal analysis; software; validation; visualization, and roles/writing-original draft; J.-S.K.: funding acquisition; methodology; project administration; supervision; and writing-review & editing. All authors have read and agreed to the published version of the manuscript. Data Availability Statement: Data sharing not applicable to this article as no datasets were generated or analyzed during the current study.

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