Inundation analysis using coupling storage function model with a distributed hydrological model in Kushiro marsh, Japan

: In Kushiro River basin, inundation is likely to occur due to heavy rain, since the river-bed slope is very gentle in the downstream portion, which has natural levees. Conse‐ quently, in past floods, river discharge was observed to increase slowly and a delay of a few days in peak river dis‐ charge was observed compared to the peak precipitation. Therefore, we applied a distributed hydrological model, Geophysical fluid CIRCulation model (GeoCIRC), to reproduce such flood discharge in Kushiro River. GeoCIRC is based on object-oriented programming and various hydrological processes, such as infiltration flow, under‐ ground water flow, surface flow, and river flow can be implemented easily. We proposed a model that can incorpo‐ rate the effect of return flow using storage function by introducing new parameters, such as storage time and time lag. This was done to consider not only the flood inunda‐ tion in Kushiro River but also the return flow from flood inundation to the river flow. As a result, we obtained high Nash-Sutcliffe coefficients for river discharge for two large flood events.


INTRODUCTION
Kushiro River is one of the longest rivers in Hokkaido Island, located in eastern Hokkaido and originates from Lake Kussharo (Figure 1(a)). The main river length is 154 km and the watershed size is 2510 km 2 . In the middle and downstream areas, there exists Kushiro Marsh, which is the largest marsh and wetland in Japan and is the first Ramsar site to be registered under the Ramsar treaty. The main stream of Kushiro River has a very gentle slope in the upstream region of approximately 1/1000, which is smaller than other rivers of Japan. Also, in the region where Kushiro Marsh exists, the gradient of river-bed slope is approximately 1/8000 to 1/3000 and contains meandering streams. Therefore, flood inundation is likely to occur in the area where river-bed slope changes rapidly, such as the Correspondence to: Shino Sakaguchi, Civil Engineering, Kobe University, 1-1 Rokkodai-cho, Nada-ku, Kobe, Hyogo 657-0013, Japan. E-mail: s.vier.no.tn0@gmail.com *Present address: Kyoto City Office, Japan starting point of Kushiro Marsh in the mountainous region. In addition, since there are no artificial dams and levees in Kushiro Marsh, flood inundation is likely to occur in heavy rain (Kushiro Wetland Restoration Committee, 2005).
The largest flood inundation in Kushiro River occurred in March, 1960 and the observed peak discharge reached 778 m 3 /s at Shibecha field observation station due to heavy rains as well as the increase in snow melt discharge. Since no levees exist in Kushiro Marsh and flood inundation occurred easily, 2204 houses were flooded and fields were inundated in 1960 (Kushiro Wetland Restoration Committee, 2005). One of the most significant characteristics of flooding in Kushiro Marsh is a few days of delay in peak river discharge from the peak of precipitation. This phenomenon is not rare in the world (Kobayashi and Takara, 2012), and was found to occur in 2011 in the Monsoon season in Chao Praya River in Thailand (Sayama et al., 2017). Chao Praya River also exhibits similar characteristics as Kushiro River, which has a very gentle slope of approximately 1/50000, without having any artificial levees in the downstream region. Since the delay in peak river discharge extends the duration of high water level, the duration of the flood inundation becomes longer due to the delay in peak river discharge when flood occurs in Kushiro River. The peak river discharge delay due to flood inundation is a serious problem not only in Kushiro River but also in other rivers, such as Chao Praya River. It is necessary to clarify the characteristics of river discharge pattern for Kushiro River during heavy rain to minimize the damage caused due to flood inundation. However, the pattern of the river discharge due to heavy rain cannot be reproduced successfully by a simple hydrological model, such as a kinematic wave model, and computed river discharge does not correspond with field observations. It is apparent that the flood inundation has not been accurately included in the general hydrological model, thus, the river discharge pattern could not be reproduced successfully.
One of the solutions of this problem is the application of a distributed hydrological model, which considers the effect of topography by using rectangular meshes and computes river flow, surface flow, infiltration flow, and so on, in each grid cell (Dutta and Nakayama, 2009;Kobayashi et al., 2012;Kobayashi et al., 2016). Bates et al. (2010) success-fully reduced computational cost of the inertial terms in the long wave equations. The other method to reduce run-time cost is the application of diffusive wave simplification, which neglects the inertial terms of the momentum equation (Yu, 2010;Leandro et al., 2014). On the other hand, Kobayashi et al. (2015) suggested that a super computer should be used to compute flood inundation with fine mesh grid size, such as 5 m, by solving long wave equations directly. Additionally, Bouragba et al. (2019) proposed the Geophysical fluid CIRCulation model (GeoCIRC) based on the long wave equations which was similar to that of Kobayashi et al. (2016). GeoCIRC is a distributed hydrological model, and the Object-Oriented Programming (OOP) is applied to include new variables, ecological models, and different flow layers, such as infiltration layer, river network layer, inundation layer, and so on. OOP was also successfully applied in a three-dimensional hydrodynamic model, Fantom, and its high applicability was shown by Nakayama et al. (2012Nakayama et al. ( , 2014Nakayama et al. ( , 2016Nakayama et al. ( , 2019. In Kushiro River, both disaster prevention and environmental protection are needed, thus, the application of GeoCIRC is one of the solutions proposed to address various kinds of problems in Kushiro River. This study aims to apply GeoCIRC into flood inundation in Kushiro River. In Kushiro River, at least 2 m mesh grid size is needed for resolving the flood inundation from branch rivers and for analyzing the delay of peak river discharge. However, since the target area size considered in this study is 2500 km 2 , the 2 m mesh grid computation becomes too expensive to carry out in practice (broken lines in Figure 1(a)). Therefore, by using the significant characteristics of OOP, we attempted to reproduce discharge by coupling infiltration flow, surface flow, and river network flow with the flood inundation modelled by the storage function method, which is a nonlinear runoff model based on a conceptual model using two parameters (Kimura, 1961). Furthermore, we aim to develop a new runoff analysis model for Kushiro River.

Study site
Our target area is the region from the middle to the downstream reach of Kushiro River that includes Kushiro Marsh (broken lines in Figure 1(a)). A mesh grid size of 500 m was used. The widths of river ranges from 20 m to 170 m from the upstream to downstream of the main stream, and the mean width of branch rivers is 10 m with the maximum width being 40 m. Six field observation stations exist in the target domain for measuring the hourly river discharge (white circles in Figure 1(a)). In this study, we aim to analyze two flood inundation events, summer flood of October, 1979 caused due to a typhoon and spring flood of April, 2013 caused by heavy rains (Table I). Figure  1(b) and (c) shows the time series of river discharge at Hirosato Station for each event. The summer flood was caused by typhoon 197920, which resulted in heavy rains from 18 to 20 October, 1979. Total rainfall measured was 166.4 mm and the maximum discharge measured at  Shibecha Station was 428 m 3 /s, which is the fifth largest discharge since field observations started. The spring flood was caused due to low pressure and snow melt from 7 to 8 April, 2013, and the total rainfall measured was 115.4 mm.

Model structure
GeoCIRC is based on the column object for onedimensional mass transport, which can include arbitrary layers vertically, such as meteorological layer, surface layer, infiltration layer, river network layer, ground water layer, and so on (Bouragba et al., 2019). Since OOP is applied in GeoCIRC, it is easy to combine necessary layers. Each layer included in this study is introduced below. Surface flow layer In this study, surface flow layer was put second from the top of the column. We considered one layer as the surface flow layer. The governing equation used for the surface flow layer was continuity equation, and the velocity was given as a function of the square root of the slope of the ground surface elevation. When the surface flow layer was wet and fulfilled the criteria of dry and wet = 1 mm water depth, excessive water in the surface flow layer was allotted to the closest cell (column) of the river network layer. When the excessive water was transported to the river network cell, the arrival time was calculated based on the distance between the surface cell and the river network cell. The size of horizontal mesh grid was 500 m each in x and y directions and the time step was 600 s. Infiltration flow layer Infiltration flow layer was put under the surface flow layer, and the governing equations used for the infiltration flow layer were the continuity equation (Equation (1)) and Darcy's law (Equation (2)), in which the finite volume method is applied. Five layers were considered in the infiltration layer of 1.0 m each (1.0 m × 5 layers). Permeability coefficient, k I , was applied by Brutsaert (1968) (Equation (4)). Saturated and residual volumetric moisture contents in Kushiro River were found to be 15% and 0%, respectively (Kushiro Wetland Restoration Committee, 2005). The size of mesh grid was 500 m each in x and y directions and the time step was 600 s, which was the same as that of the surface flow layer.
where, t is the time (s), ρ w is the density of water (kg m -3 ), h is the water depth (m), m I is the volumetric moisture content, (u I , v I , w I ) are the velocities in the (x, y, z) coordinates (m s -1 ), ϕ is the velocity potential (m), φ is the suction (m), k 0 is the saturated permeability coefficient (m s -1 ), m min is the residual volumetric moisture content, m max is the saturated volumetric moisture content, and β is the parameter for soil. River network layer A river network layer was put at the top of the column.
One-dimensional long wave equations were applied (Equations (5) to (6)) to solve the river flow. A staggered grid was applied using the finite difference method, and the predictor corrector method was applied to solve the long wave equations stably, such as flow including hydraulic jump (Bouragba et al., 2019). Equations (7) to (11) are the numerical schemes used, and the length of mesh grid was 500 m each in x and y directions, and the time step was 10 s.
where u is the velocity (m/s), ɡ is the gravity acceleration (m/s 2 ), n m is the Manning roughness (= 0.03 s/m 1/3 ), u i n is the velocity at the node i (m/s), n corresponds to time, Δt is the time step (s), û and ũ are the predictor and corrector (m/s), and z 0 is the height from the base level.

Mass conservation in a column
First, computation of each layer was conducted horizontally, and then, vertical transport of mass was computed in the infiltration flow layer. Second, vertical mass transport was computed from the top to the bottom layer by considering the vertical permeability. Last, mass conservation was computed from the bottom to the top layer. For example, when the infiltration flow layer contained excessive water, the excessive water was transported to the surface flow layer. Consequently, if the surface flow layer contained excessive water, the excessive water was transported to the river network layer. This mass conservation process was carried out during computation for each time step.

Determination of initial conditions
Since Kushiro Marsh falls in the downstream region of Kushiro River, volumetric moisture content in the infiltration flow layer was mostly filled with water, and thus, the main source of base flow discharge was given by the infiltration flow layer in Kushiro River. Therefore, we filled all the layers with water and carried out 365 d computation without any precipitation to obtain initial conditions for base flow discharge. Afterwards, we selected the volumetric moisture content for arbitrary base flow discharge as initial conditions for analysis. We decided the saturated permeability coefficient by trial and error by comparing three reduction curves of observed discharge that were not affected by precipitation from May to June in 2013.

Computational conditions
Case 0: no inundation To understand the effect of flood inundation on river discharge, we assumed that no inundation exists, implying that the height of the levee was sufficient to prevent flood inundation. Case 1: inundation without return flow To estimate the volume of flood inundation, we assumed that any excessive overflow of water over the natural river banks never returns to the river. The overflow was considered only in the main stream, and the height of the natural river bank was considered to be 3.0 m. Case 2: inundation using storage function model Since volumetric moisture content in Kushiro Marsh was mostly saturated, excessive water that overflowed the natural river bank returned to the river flow without penetrating into the ground. Therefore, in case 2, first, the excessive water was accumulated in the virtual cell just above the river network cell where excessive water appeared. Afterwards, some of the accumulated water was returned to the river network by using a storage function model (Equations (12) and (13)). In the storage function model, we assumed that discharge is a linear function of storage, and discharge was rewritten as Equation (14) by introducing the new parameter, storage time (= T S ). The solution of Equation (14) is given in Equation (15) when the initial discharge is S 0 , which suggests that T S is the representative time scale for excessive water.
where, S is the accumulated excessive water from the river network (m), Q is the discharge that returns to the river network (m s -1 ), K and P are the parameters for the storage function model, T S is the storage time (s), and S 0 is the integral constant (m). Case 3: inundation using storage function model with time lag In case 2, we considered that the accumulated water returns to the river network without any delay. However, Figure 1(b) and (c) shows the time lag when accumulated excessive water starts to return to the river flow. Therefore, we assumed in case 3 that the accumulated excessive water started returning with the time lag Tl.

RESULTS AND DISCUSSION
In case 0, where we assumed that flood inundation does not occur, river discharge agreed with four field observation stations in the branch rivers, Shimoosobetsu, Shimokuchoro, Setsuri, and Hororo Stations, for both summer and spring events (Figures 2 and 3). Therefore, it was found that the branch rivers did not inundate, since field observation stations in the branch rivers were located in a steeper slope region compared to Kushiro Marsh. On the other hand, the computed peak river discharge was found to be significantly higher than the field observation at Hirosato Station for both summer and spring events, i.e. approximately 1100 m 3 /s and 500 m 3 /s, respectively. In addition, the computed peak river discharge occurred approximately three days earlier than field observations. In contrast to case 0, river discharge at Hirosato Station was very small for both summer and spring events in case 1 (Figure 4(a) and  (b)). Therefore, we suggested that both flood inundation and the return flow from flood inundation to river flow should be taken into account.
In case 2, where flood inundated water returned to the river flow by using the storage function model with storage time as a function, peak river discharge, arrival time of the peak river discharge, and the Nash-Sutcliffe coefficient were computed (Figure 4(c), (d) and Table II) (Nash and Sutcliffe, 1970). The best values of storage time for summer and spring events were decided by trial and error. The storage times were obtained as 10 d and 20 d for summer and spring events, respectively, with the highest Nash-Sutcliffe coefficient. The computed peak river discharge was found to be almost equal to field observations. However, the arrival time of peak river discharge was not similar to field observations and it was found that the arrival time of peak river discharge did not change with the change in the storage time for summer and spring events. Therefore, it was demonstrated that the storage function model, without the time lag Tl, was insufficient to reproduce the arrival time of peak river discharge.
In case 3, where time lag for the storage function model is introduced, Nash-Sutcliffe coefficients were found to be 0.90 and 0.87 for summer and spring events, respectively, when Tl = 1.5 d (Figure 4(e), (f) and Table Ⅱ). Similar to case 2, the storage time was decided by trial and error to obtain the best fit to the time series of observed river discharge. For a summer event, the storage time was calculated as 10 d with the best fit to the peak river discharge and the arrival time of the peak river discharge. For a spring event, the best value of the storage time was obtained as 5 d, however, the entire discharge was found to be slightly smaller than field observations. This is because we did not include additional river discharge as a result of snow melt. Dutta and Nakayama (2009) showed that river flow can respond to rainfall quickly when water surface elevation of the floodplain compartment is closer to the (e) Hirosato station river compartment. As Kushiro River basin holds more water in spring than summer due to snow, the storage time may be given as shorter in spring than in summer. It should be noted that we need to investigate the grid size effect on river discharge in the future study because Dutta and Nakayama (2009) revealed that the spatial grid resolution greatly affects the simulated peak river discharge because of the changes of average slope and river length.
In conclusion, the Kushiro River basin model was successfully developed by applying OOP. The larger the storage time is given, the smaller the peak river discharge is obtained. It can be concluded that the storage time does not have significant influence on the arrival time of the peak  river discharge. Therefore, time lag should be included when applying the storage function model. Time lag may be associated with the flood inundation area and the other topographical characteristics of a river basin. Furthermore, since a summer event consists of a higher total rainfall compared to a spring event, the best fit value of storage time for a summer event is longer than that for a spring event.