2.1 Distributed Runoff Model

In this study, the distributed runoff model consists of four types of flows: surface flow, subsurface flow, infiltration flow, and river channel flow. The kinematic wave method was applied to the surface flow and the subsurface flow. The surface flow equations and the subsurface flow are described by the following equations:

where h$_{s}$ is the water depth on the local surface, Hs is the height of water from the datum, qsx and qsy are the unit discharges in the x and y directions, respectively, r is rainfall, n is Manning’s roughness parameter, and f is the infiltration.

where h$_{g}$ is the water depth in the ground, H$_{g}$ is the height of water from the datum, q$_{gx}$ and q$_{gy}$ are the unit discharges in the x and y directions, respectively, ${\phi}$ is the porosity, and k$_{a}$ is the permeability coefficient.

For the infiltration f, we applied the Green-Ampt infiltration model:

where kv is the effective saturated hydraulic conductivity, $\theta $ is the initial water volume content, $S_{f}$ is the suction at the vertical wetting front, and F is the cumulative infiltration. To solve Eqs. (1)-(7), the finite element method and the fifth-order adaptive-time Runge-Kutta method were applied to approximate the spatial direction and time direction, respectively. The relationship between the surface water depth, groundwater depth, rainfall, and infiltration in each grid cell is shown in Fig. 1.

For the channel flow, a one-dimensional version of the equation for the surface flow was used. However, it is necessary to consider the cross-sectional area of the channel flow because the amount of change in water depth depends on the width of the channel at each point, even when the same amount of water flows into the channel. In this study, the cross section of the channel was assumed to be a rectangle, and its width and depth were determined by using a of the catchment area ^[4]. The width W and the depth D of the river channel at a certain point are represented by Eqs. (8) and (9):

where Cw, Sw, Cd, and Sd are parameters for each river, and A is the catchment area for each point.

The water supply to the river channel was calculated using the inflow from the surface flow. The outflow from the subsurface flow was also added to the surface water depth as a return flow when it exceeds the depth of the seepage layer so that the final amount of water flowing into the river would be determined by only the surface water depth. We adopted a step-down equation ^[4] in which inflow occurs when the river water level is lower than the ground surface. The lateral inflow into the river channel at that time, q, is calculated by the Eq. (10):

where g is gravitational acceleration.

Data assimilation is also required to update the state of the model each time a value is observed. In this study, optimal interpolation was used as the data assimilation method. An optimal interpolation method is widely used as a simple data assimilation method and in real-time forecasting because it is considered to have a lower computational cost at each time step than other methods ^[9]. In the field of runoff analysis, Miyake et al. ^[6] applied an optimal interpolation method to the runoff calculation of a distributed model and demonstrated its usefulness by calculating and comparing the water level one hour later from the pre-assimilated and assimilated states.

In the optimal interpolation method, the assimilated value (analytical value) x$^{{a}}$ is expressed by Eq. (11).

where x$^{{b}}$ and y are first estimates (pre-assimilation values computed by the model) and observed values, respectively. B is the covariance matrix of background errors, R is the covariance matrix of observed errors, and H is an observation matrix representing interpolation from model space to observation space.

Both R and B need to be determined in advance. We set the background error variance to 1.0 and the observed error variance to 0.5 based on a previous study ^[6]. For the covariance matrix of the observation errors, we assumed that each observation error is independent and uncorrelated, and we set the matrix as a diagonal matrix. The spatial correlation coefficients of the background error covariance matrix were calculated from the water levels in each grid cell at the peak time of the multiple simulations.

The prediction method by the physics-based model in this study has the following steps.

[Step 1.]Generate grid cells as a basis for calculation using elevation, flow direction, and catchment area data and then determine the river channel parameters using the cross-sectional profile of the target river.

[Step 2.]Prepare flood data for calculating the weights of the optimal interpolation and calculate the runoff using the rainfall data at the time. Calculate the weight matrix W of the optimal interpolation using the difference between the calculated water level and the observed water level at the peak time.

[Step 3.]Optimize the parameters using the training data. The initial water level of the channel mesh is calculated using the observed water level W at each station in Step 2 and Eq. (11). A quasi-Newtonian method is used for optimization and is updated based on the errors between the calculated water level and the observed water level.

[Step 4.]Calculate the runoff using the optimized parameters. Every 10 minutes, the calculation step proceeds using the observed rainfall data. When the water level is observed every hour, the data are assimilated using the optimal interpolation method. However, the assimilated values are only used in Step 5 and are not reflected in the next calculation steps.

[Step 5.]Make predictions using the assimilated values. Apart from Step 4, simulations are run using the assimilated values as the river water depths. In this simulation, the surface and subsurface water depths are taken from the data in Step 4, and the observed rainfall is used as the predicted rainfall data. The calculated water level is regarded as the predicted value.

2.2 Neural Network Model

For the machine learning model, we constructed a model based on a DNN ^[7]. A neural network is a model that consists of many layers of artificial neurons that imitate neurons in the human brain. They are generally called an artificial neural network when there is one intermediate layer and a DNN when there are two or more layers. An artificial neuron is a simple model in which the output value is determined by a function called the activation function and the weighted value of each input signal. The basic structure of a DNN is shown in Fig. 2, and an artificial neuron is represented by Eq. (12):

where x is an input variable, w is a weight parameter, $\theta $ is a threshold value, and u is an activation function.

In this study, we used the ReLU function as the activation function of the intermediate layer and the constant function as the activation function of the output layer. The ReLU function outputs the input signal as it is when the input signal is positive and 0 when the input signal is negative. It is widely used as an activation function for the intermediate layer because of its good computational efficiency and stable learning. The water level prediction method by machine learning in this study has the following steps.

[Step 1.]Determine the network structure according to the forecast condition. The input layer is determined based on the data to be used for the forecast, such as the number of stations in the watershed. The output layer is determined by the forecast target, such as the number of time points to forecast.

[Step 2.]Train the model using multiple flood cases.

[Step 3.]The observed water level, observed rainfall, and predicted rainfall are given to the trained model at each time, and the amount of change from the water level at the current time is calculated. The predicted values are obtained by adding them to the current water level. Since it is necessary to consider the model performance and the predicted rainfall separately, the rainfall prediction is assumed to be a perfect prediction, and the observed rainfall is used for training and forecasting.

4.1 Target Watersheds and Used Data

The target watershed is Hiwatari water level station in the Oyodo River basin, Miyazaki Prefecture, Japan ^[7]. Hiwatari water level station has a watershed area of 816 km$^{2}$, as shown in Fig. 3, and is surrounded by mountains and hills. There are 5 water level stations and 14 rainfall stations in the upstream area.

In the machine learning model, rainfall amounts and water levels obtained at the stations were obtained from the sluicegate water-quality database. Elevation data and radar rainfall data were used to construct a distributed model. C-band radar was used for the radar rainfall data, and the functions of the RRI model were used to obtain the elevation data and to partition the mesh.

4.2 Experimental Conditions

Manning's roughness coefficient n and porosity $\phi $ for the slope and river have a great influence on the velocity of the downstream flow, flood arrival time, and peak water level, so they were considered for the optimization. Since the approximate ranges of these values has been clarified by previous studies ^[11,12], we used them for optimization. For other parameters, common values in studies of distributed runoff models were used. Let n$_{s}$ be the roughness coefficient in the channel and n$_{\mathrm{r}}$ be the roughness coefficient on the slope. The parameters are summarized in Table 1.

For the calculation of the background error variance for data assimilation, the peak water levels of four cases in September 2011, June 2012, July 2012, and June 2016 were used. The case of July 2007 was used for parameter optimization. The predicted water level was calculated by correcting the calculated water level with the optimal interpolation method and doing the runoff calculation separately. Therefore, the observed water levels do not affect the original runoff calculations and are reflected in only the predictions. When constructing the deep learning model, we used a DNN with two intermediate layers based on another model ^[7]. The detailed settings of each layer and the learning conditions were determined by trial and error, as shown in Table 2.

For the features, we used the water level at 1 hour before and at the current time of the target station, the water level fluctuation up to 4 hours before from each station, and the rainfall from 4 hours before to 9 hours after. The output layer corresponds to the water level change from 1 hour to 9 hours after. As training data, we used 5 days of data from the top 20 cases for the period of 2007-2016.

4.3 Results

Figs. 3 and 4 show the predicted water level from 1 to 9 hours after the observation time. Fig. 3 shows the results of predictions by the physics-based model, and Fig. 4 shows the results of predictions by the machine learning model. Figs. 3(a-1) and (b-1) were obtained by using MSE in parameter optimization, while Figs. 3(a-2) and (b-2) were obtained by using WHL in parameter optimization. Figs. 5(a-2) and (b-2) show the results of the model predictions with parameters optimized by WHL, superimposed on the observed water levels. Figs. 5(c-1) and (c-2) show the results of the physics-based model and the machine learning model, respectively.

The prediction by the model optimized by MSE had a small error in the rising part of the water level, but the underestimation near the peak was significant, which indicates that the parameters were optimized so that the overall error is small. For the model optimized by WHL, the peak water level and the rising part of the water level were overestimated, and there was some increase in the overall error. However, the peak flood was predicted earlier and with a slight overestimation. As mentioned, the prediction for the extended lead time converged to a certain value, but this value was also overestimated compared to that of MSE.

In Fig. 5, we can see that the prediction by the physics-based model converged to a certain value even when the lead time was extended. This was due to correcting only the depth of the river channel sequentially using the observed water level. When the lead time was short, the effect of the correction was strong because the upstream water corrected by the optimal interpolation method flowed downstream, but the effect of the correction became smaller as the time step was increased. As time passed after the correction, the effect of the runoff component due to soil moisture increased, yet the values were almost constant because they were not corrected by the optimal interpolation method. These constant values were almost consistent with those obtained in the overall simulation, and the correction was effective only for prediction up to two or three hours ahead.

The prediction results by the machine learning model show that the prediction values were generally scattered and that the accuracy varied with the lead time. According to the prediction results for each lead time, the prediction accuracy worsened as the lead time increased, although it was almost error free in the short-term prediction. It is also considered that the contribution of rainfall increases as the lead time increases, while the change of water level in the upper stream contributes significantly to the prediction accuracy in the short-term prediction.

Table 3 shows the values of various evaluation indices for the prediction results. In terms of the overall error magnitude, the prediction by the physics-based model had a larger error than that by the machine learning model, and the error was about twice as large as that by the machine learning model, regardless of the lead time. This is considered to result from the influence of the channel geometry in the physics-based model. Since the physics-based model calculates the amount of water movement between meshes using the difference in elevation and channel depth, it often maintains a constant value when the fluctuation of water depth is moderate, such as before the rise of water level or after the peak. This constant value was determined by the difference in elevation from the surrounding mesh and the depth of the river channel, but in the distributed model used in this study, the channel cross-section was assumed to be rectangular, and there was a discrepancy with the actual cross section. As a result, the model stabilized at a higher position than the actual water level at many points in time, resulting in a large overall error.

In terms of the peak time, the predictions by the physics-based model were stable, while those by the machine learning model were unstable, and the largest delay was three hours. It is considered that the peak time of the predictions by the physics-based model did not change even if the water depth in the river channel was corrected by the optimal interpolation method because the water depth at the ground surface and underground determines the discharge into the river channel. The peak time of the prediction by the machine learning model was unstable because the model does not know the storage volume of the entire watershed, although it uses rainfall data before and after the starting time of the prediction.

For the difference by the objective function, MSE had a noticeable underestimation error in the peak water level, although the overall error was smaller. In contrast, WHL had little underestimation of the peak water level, and both models showed improvement. However, WHL had a larger overall error than MSE and increased RMSE by up to 25%.