Structural improvement of a kinematic wave-based distributed hydrologic model to estimate long-term river discharge in a tropical climate basin

: A distributed hydrologic model based on a kinematic wave approximation with surface and subsurface flow com‐ ponents is applicable to basins that have temperate climatic conditions similar to those in Japan. However, it is difficult to present long-term river discharge using the existing model structure in basins with different climatic conditions. This study aims to improve the model structure for better estimates of long-term discharge in the Nam Ngum River, the main tributary of the Mekong River, by incorporating bedrock aquifers as part of the slope flow component of the original model structure. Three bedrock groundwater struc‐ tures are configured to incorporate the original model structure. The results show that a combination of the origi‐ nal model component and one unconfined aquifer structure are the best representations of the river flow regime from the original model structure, in which the rate of infiltration from the layer into the bedrock aquifer was calculated using vertical hydraulic conductivity. The Nash–Sutcliffe efficiency coefficient of the original and improved models increased from 0.80 to 0.86 during the calibration period and from 0.56 to 0.62 during the validation period. The results of this study show that the improved model structure is applicable for long-term hydrologic predictions in South‐ east Asian catchments with distinct dry and rainy seasons.


INTRODUCTION
Rapid regional growth and energy demands from neighboring countries have prompted plans to build numerous dams along the mainstream and tributaries of the Mekong River (Kummu and Varis, 2007). A dam reservoir, which controls rivers for both water use and flood control, can play a significant role in effectively managing water resources (Nohara et al., 2016). Thus, the efficient operation of large-scale water infrastructure such as dam reservoirs requires a good estimation of long-term river discharge for reservoir operation plans, water resource management, and flood control.
The distributed hydrologic model based on a kinematic wave approximation with surface and subsurface flow com-Correspondence to: Thatkiat Meema, Department of Civil and Earth Resources Engineering, Kyoto University, C1, Nishikyo-ku, Kyoto 615-8540, Japan. E-mail: meema.thatkiat.42r@st.kyoto-u.ac.jp ponents (DHM-KWSS) is a well-known model used to describe rainfall-runoff processes for many river basins in Japan where most of the basins have steep slopes (Takasao and Shiiba, 1988;Sayama et al., 2006;Kim et al., 2011). The DHM-KWSS is applicable to basins with climatic conditions similar to those of such Japanese basins. However, describing hydrological behavior under different conditions, such as arid basins, is difficult (Hunukumbura et al., 2012). Tanaka and Tachikawa (2015) applied 1K-DHM with a DHM-KWSS structure (Tachikawa and Tanaka, 2013) to two river basins in Europe and Australia that have different climatic conditions than those in Japan for a longterm river flow simulation. The results of this study showed that the model structure of the DHMs-KWSS is applicable to river basins in temperate climatic conditions but has difficulty describing rainfall-runoff processes in semiarid basins. Tanaka (2016) has suggested that to improve the long-term river flow estimation, especially in the dry season, future studies should incorporate a groundwater component into the model structure. Katsura et al. (2008) found that the infiltration of water from soil into the weathered bedrock was the dominant factor and that the underlying bedrock makes important contributions to water flow. Ferket et al. (2010), Samuel et al. (2012), and Luo et al. (2012) adopted a reservoir or multireservoir with linear, non-linear, or combined relationships to generate baseflow in hydrologic models.
In this study, the 1K-DHM was applied to estimate longterm river discharge in a tropical climate basin -the Ngum River basin in Lao People's Democratic Republic. The results show that explaining the hydrological behavior of a tropical climate basin using the original 1K-DHM model structure that includes only surface conditions (i.e. topography and land cover) and surface soil layer properties is difficult. Therefore, the main objective of this paper is to improve the DHM-KWSS structure for estimating longterm river discharge in a tropical climate basin by incorporating three bedrock groundwater structures, including unconfined and confined aquifers into the 1K-DHM. The parameter values for each model structure are identified using the SCE-UA algorithm, and the performance of models with three aquifer structures was evaluated for five years, not including identification periods.

APPLYING 1K-DHM TO THE NAM NGUM RIVER BASIN
The watershed boundary of the Nam Ngum River at the Naluang station, the depth to bedrock (DTB) map and the map of shallow groundwater storage (Viossanges et al., 2018) are shown in Figure 1. Text S1 describes how the 1K-DHM was applied to the Nam Ngum River basin. The simulation optimized soil depth parameters -such as effective soil porosity depth (d a ) -to be much larger than the average DTB of the basin as shown in Figure 1 so as to maintain river discharge in the dry season. The range of porosity for an unconsolidated deposit is approximately 0.25-0.7 (Freeze and Cherry, 1979); thus, the range of d a should be from 0.8-2.2 m. However, the calibrated soil depth of the original 1K-DHM structure is much larger (the value is shown in a later section). The thick soil layer also resulted in an underestimation of the river discharge during the transition from the dry to the wet season. Thus, the existing model structure that considers only surface conditions (i.e. topography and land cover) and the properties of the soil surface layer does not explain the hydrological behavior of a tropical climate basin that has distinct dry and rainy seasons. Therefore, this study aims to improve the model structure of the 1K-DHM for a better estimation of long-term river discharge in the basin.

IMPROVEMENT OF THE MODEL STRUCTURE
Shallow and deep bedrock aquifers are analyzed in this study based on the literature on conceptual aquifer models described in the introduction section. Three aquifer structures are configured to incorporate the 1K-DHM by adding a single or multi aquifer of bedrock into the surface-subsurface component (soil surface) of the original model structure. The schematic of each slope flow model and discharge-storage relationship for each component are summarized in Figure 2, and the description of each structure is provided below.

Model structure 1 (M1)
M1 consists of the surface-subsurface flow, shallow aquifer, and deep aquifer components on the slope unit. Surface-subsurface flow component The surface-subsurface flow component has the same structure as the original 1K-DHM (M0) (explained in Text S2). Thus, the discharge-storage relationship applies the same formula as M0, as expressed in Equation (S1). To consider the vertical infiltration from the surfacesubsurface flow component into the shallow aquifer bedrock, the continuity equation of the soil surface layer is modified as follows: where t is time, x is the space coordinate, r is rainfall intensity, e is actual evapotranspiration (AET), h s is the water stage in the soil surface component, q s is the discharge per unit width of the soil surface component, and p u is the vertical infiltration rate from the surface soil layer into the shallow aquifer. Unconfined aquifer component The shallow groundwater aquifers of river basins are predominantly unconfined, suggesting that the relationship between baseflow and storage is nonlinear (Wittenberg, 1999). Banks et al. (2009) described the flow through bedrock as being controlled by the fracture network and connectivity. This study assumed that some part of the groundwater flowed out from the soil layer boundary, infiltrated into the bedrock, flowed through the underlying bedrock and the connected rock fracture, and contributed to Figure 1. General physical condition of the Nam Ngum River Basin (a) watershed boundary with rain gauges location; (b) DTB map of Naluang catchment (Shangguan et al., 2017); (c) aquifer storage map of Naluang catchment (Viossanges et al., 2018) IMPROVEMENT OF KINEMATIC WAVE-BASED HYDROLOGIC MODEL the river discharge. According to Darcy's law, the total discharge per unit width of the unconfined aquifer q u (h u ) is expressed as: where d u is the total effective depth of rock fracture in the unconfined bedrock aquifer as shown in Figure S3; h u is the total water stage in the fracture of the aquifer; k u is the hydraulic conductivity that corresponds to the actual crosssectional area of flow in the rock fracture. In the unconfined aquifer, the hydraulic conductivity is assumed to decrease when h u decreases; therefore, the hydraulic conductivity is assumed to be k u (h u ) = k u (h u /d u ) β u , where β u is a parameter that corresponds to the reduction of hydraulic conductivity in the unconfined aquifer (in this study , β u is assumed to be 1), i is the gradient of the hillslope, and α u = k u i/d u .
The unconfined aquifer in a slope unit is calculated by the following continuity equation: where p c is the vertical infiltration from the unconfined aquifer to the confined aquifer. Confined aquifer component A deep aquifer can be defined as a confined aquifer with a linear discharge-storage relationship. For the confined aquifer, d c is the total effective depth of the rock fracture in the confined bedrock aquifer, and h c is the total water stage in the rock fracture of the aquifer. k c is the hydraulic conductivity that corresponds to the actual cross-sectional area of the flow in the rock fracture. Thus, the total discharge per unit width q c (h c ) can be calculated by: where α c = k c i. The confined aquifer in a slope unit is calculated by the following continuity equation: The total discharge from all slope flow components on both sides of the river channel (2q s , 2q u , and 2q c ) contributes to the river channel component as the lateral discharge per unit length. Consequently, the continuity equation of a river channel flow model was modified as follows: where A is the river cross-section area, Q is the river discharge, and the relation Q = αA m (α = √i o /(nB m-1 ) where n is the Manning's roughness coefficient of the channel; i o is the slope gradient of the channel and B is the channel width, m = 5/3) is assumed based on the kinematic wave flow assumption.

Model structure 2 (M2)
The structure of M2 is similar to that of M1 but consists of only an unconfined aquifer. Thus, M2 consists of two main components on the slope unit -the surface-subsurface flow and the unconfined aquifer components.
The surface-subsurface flow component is calculated with the same discharge-storage relationship as M1 with the continuity equation in Equation (1). For the unconfined aquifer component, the discharge per unit width (q u ) is calculated using Equation (2) and the continuity equation in Equation (3), where the vertical infiltration from the unconfined aquifer to the confined aquifer (p c ) is set to 0. For a river channel flow model, the total lateral discharge from the confined aquifer component (q c ) in Equation (6) is set to 0.

Model structure 3 (M3)
The structure of M3 is similar to that of M2 and consists of two main components on the slope unit: the surfacesubsurface flow and the unconfined aquifer components. The difference between M3 and M2 is the estimation of the vertical infiltration from the soil surface layer to the uncon-

Estimation of vertical infiltration
To estimate the amount of vertical infiltration into the bedrock, the storage routing technique (Arnold et al., 1998) was adopted. The storage routing technique is based on the equation where S (t=0) and S (t=1) are the water contents at the beginning and end of the calculation time steps, respectively, Δt is the time step, and T is the infiltration time of water from the upper layer. Thus, by subtracting S (t=1) from S (t=0) , the amount of vertical infiltration (P) can be calculated with the equation: The estimation of the vertical infiltration rate from the soil layer to the unconfined aquifer (p u ) in this model follows the assumption that the infiltration process occurs when the water content in the soil layer exceeds the unsaturated flow condition (h s > d m ). Thus, by adopting the storage routing technique, p u can be calculated as follows: where d m is the capillary depth, and T is calculated by T = (h s -d m )/k. For M1 and M2, the saturated hydraulic conductivity in the soil layer (k a ) is used for k; the vertical hydraulic conductivity (k v ) is used for M3. Only the M1 structure considers the confined aquifer, and the vertical infiltration from the unconfined aquifer to the confined aquifer (p c ) is calculated as follows: where T = h u /k u . The infiltration from the upper component does not occur when the storage in the aquifer component is filled full of water.

MODEL CALIBRATION PROCESS
The simulation is divided into the calibration and validation stages. Due to data availability and the condition of the basin (without upstream dam reservoirs), 1990 was selected as the calibration period, and 1991-1995 was selected as the validation period. For each model, the SCE-UA algorithm (Duan et al., 1994) was applied to optimize the model parameters by searching for the parameter with the highest Nash-Sutcliffe efficiency coefficient (NSE). The important preparation to identify the parameter set is to design the range of the parameters. The details of the model calibration process are described in Text S3.

RESULTS AND DISCUSSION
According to the optimized parameter sets for the Nam Ngum basin in Table I, the effective soil depth parameters such as d a and d m are the primary difference between M0 and M1-M3. M0 attempted to simulate the base flow in the dry season by reflecting a soil depth value much larger than that of M1-M3. The comparison of optimized parameters between M1-M3 simulated in the Nam Ngum Basin and M0 simulated in Japanese basins shows no large difference in the effective soil depth parameters. Nevertheless, M0 for Japanese basins was simulated based on large flood events. Therefore, for the long-term simulation of M1-M3 in the

IMPROVEMENT OF KINEMATIC WAVE-BASED HYDROLOGIC MODEL
Nam Ngum Basin, β is set much higher to maintain the base flow during the dry season. The observed and simulated discharge hydrographs during the calibration and validation stages obtained from each model structure at the Naluang station are compared in Figure 3. Figure 4 shows the model performance indices, such as the NSE, Root mean square error (RMSE), and Percent bias (PBIAS), of each model during the calibration and validation stages.
All the models presented the characteristics of long-term river discharge in the basin with distinct dry and rainy seasons. However, M0 attempts to maintain the base flow during the dry season by placing an emphasis on both soil and bedrock groundwater processes in only the soil layer component. This resulted in the overly thick d a value shown in Table I.

Improvement by M1 and M2
According to the simulation results, the M1 and M2 hydrograph and performance indexes showed similar tendencies. The combination of discharge generated by the confined and unconfined aquifers in M1 is similar to the discharge generated by the unconfined aquifer in M2. Moreover, the hydrologic processes in the soil layer component of the model structure are the same. On the one hand, the result of the river flow hydrograph at the basin outlet in M1 is very similar to that of M2. On the other hand, the incorporation of the bedrock aquifer into the 1K-DHM structure in M1 and M2 yielded better model performance indices than M0 for NSE and RMSE for all periods (annual, wet, and dry seasons) during the calibration stage. NSE and RMSE also showed a better performance in M1 Figure 3. Hydrograph of observed (black line) and simulated discharge (red line) at Naluang station including discharge generated by the unconfined aquifer (dash red line; in M1 means combination of unconfined and confined aquifers) and confined aquifer (dot red line) and M2 than in M0 during the validation stage.

Disadvantages of M1 and M2
For M1 and M2, the optimized processes control the amount of infiltration by reflecting a high value of k a to maintain unsaturated flow conditions in the soil layer. When high rainfall intensity causes the water stage in the soil layer (h s ) to exceed d m , saturated flow conditions occur, and the vertical infiltration process begins. Most rainwater was infiltrated into the bedrock aquifer because of high rates of k a that caused the simulated river discharge at the beginning of the rainy season to be underestimated during the calibration stage ( Figure 3). Additionally, the recession part of simulated hydrograph during the calibration period was higher than the observed hydrograph because of a slow discharge rate in the aquifer and the time required to release the infiltrated water that is stored in the aquifer components. Therefore, the calculation of the vertical infiltration component is required to improve the simulation.

Improvement by M3
As mentioned above, the contribution of the confined aquifer in M1 is substituted by that of the unconfined aquifer in M2 for this study area. This suggests that developing M3 based on the M2 structure has more advantages than M1 because, to consider the confined aquifer in the model as M1, extra parameters are required which least to a more complicated calibration process.
To appropriately estimate the amount of infiltration, k v was introduced into M3. The reduction in the amount of infiltration not only improved the simulation of base flow during the recession period, but also captured the peak flow at the beginning of the rainy season during the calibration stage ( Figure 3). Therefore, M3 provided the best perfor- mance for all indices such as NSE, RMSE, and PBIAS for all periods (annual, wet, and dry seasons) among the modified models during the calibration period. For the validation periods, the average annual NSE and RMSE in M3 are quite similar to those in M1 and M2; however, PBIAS is smaller, indicating that M3 improved the simulation of river discharge obtained in M1 and M2.

Prediction uncertainties due to limited data accuracy
Accurate data is necessary to evaluate the performance of the model. One of the reasons for the universally poor performance of the models in 1994 may be the uncertainty in the observed data. From Aug 29-Sep 15, the amount of rainfall in the basin was 69 mm; however, the amount of runoff in the basin reached 290 mm (more than 4 times the rainfall amount) (Figure 3). These questionable results suggested underestimation of river discharge and explained the poor performance of the model.
Moreover, the lack of rainfall stations in the study basin may directly affect the performance of the model. By using the nearest neighbor method to distribute the amount of rainfall from 4 stations over the Naluang catchment, we found that only 3 main stations contributed rainfall values over the catchment as input for the model. This may lead to uncertainty in the amount and spatial distribution of rainfall estimation − one of the most important factors in the distributed hydrologic model.
The overestimation in the dry season by all the models may be caused by the AET estimation. In this study, AET was roughly estimated by adopting the difference between annual rainfall and discharge in the basin and using the average of this amount as daily AET input into the model. Consequently, AET might be underestimated during the dry season, which, in turn, may have resulted in overestimation of low flow during the dry season.

CONCLUSIONS
The model that combined the soil layer and the unconfined aquifer component with the estimation of the vertical infiltration based on k v (M3) best reproduced the river in this study. The NSE values for M3 at the calibration and validation stages were 0.86 and 0.62, respectively, and notably improved the results of the original model (M0).
Therefore, by incorporating the bedrock aquifer into the DHM-KWSS model structure, the estimation of long-term river discharge − especially low flow in the dry season − was explained and improved remarkably. The model improved river discharge estimation and produced a reasonable set of parameters that agreed with physical data sets such as DTB and groundwater storage in the study area. Furthermore, the result of k v value seemed to be reasonable. Marechal et al. (2004) determined that a horizontal-to-vertical anisotropy ratio for hydraulic conductivity close to 10 (k h /k v~1 0) in the fractured rock aquifer. In this study, k u /k v is approximately 70 (consider k u through the total aquifer depth) which shows a large difference between k u and k v . This may be caused by the fact that k v responded as the infiltration threshold related to the input rainfall intensity (if rainfall intensity is smaller than k v , all of rain water is infiltrated). Due to data limitations, the input rain-fall intensity is estimated by using 24-hour mean rainfall (in which k v of 1.16 × 10 -7 m/s is approximately 12 mm/24 hours). On the other hand, assuming that the 12 mm of rainfall occurs in 4 hours (assume that this value is the infiltration threshold), the input rainfall intensity might increase 6 times; then k v values could be increased and k u /k v potentially close to 10.
The estimation of rainfall-runoff in the hydrologic model is subject to uncertainties in the accuracy of the data, such as the discharge and rainfall amount measurements, and the number of rain gauge stations in the basin. Furthermore, a better estimation of AET is required to improve long-term river flow prediction in basins with distinct dry and wet conditions. Additionally, further study is required to test the application of the improved model in other basins that have similar or different climatic conditions.

ACKNOWLEDGMENTS
We acknowledge the TEAM consulting engineering and management for providing the hydrological data. This work was supported by MEXT.

SUPPLEMENTS
Text S1. Description of applying 1K-DHM to the Nam Ngum River basin Text S2. General description of 1K-DHM Text S3. Details of the model calibration process Figure S1. Estimation of actual evapotranspiration Figure S2. Schematic of flow simulation in 1K-DHM Figure S3. Conceptual of soil and bedrock aquifer models Figure S4. Finite difference mesh used in kinematic model