Flood Inundation Assessment in the Low-Lying River Basin Considering Extreme Rainfall Impacts and Topographic Vulnerability

Abstract

This study aims to evaluate the change in flood inundation in the Chitose River basin (CRB), a tributary of the Ishikari River, considering the extreme rainfall impacts and topographic vulnerability. The changing impacts were assessed using a large-ensemble rainfall dataset with a high resolution of 5 km (d4PDF) as input data for the rainfall–runoff–inundation (RRI) model. Additionally, the prediction of time differences between the peak discharge in the Chitose River and peak water levels at the confluence point intersecting the Ishikari River were improved compared to the previous study. Results indicate that due to climatic changes, extreme river floods are expected to increase by 21–24% in the Ishikari River basin (IRB), while flood inundation is expected to be severe and higher in the CRB, with increases of 24.5, 46.5, and 13.8% for the inundation area, inundation volume, and peak inundation depth, respectively. Flood inundation is likely to occur in the CRB downstream area with a frequency of 90–100%. Additionally, the inundation duration is expected to increase by 5–10 h here. Moreover, the short time difference (0–10 h) is predicted to increase significantly in the CRB. This study provides useful information for policymakers to mitigate flood damage in vulnerable areas.

1. Introduction

Extreme weather events associated with climate changes have rapidly become one of the global concerns threatening natural environments and human life in recent decades. According to the IPCC Fifth Assessment Report [1], the global average surface temperature increased by about 0.73 °C over the period of 1891 to 2017, and this trend is likely to continue into the new millennium. In addition, precipitation is one of the most critical and useful indicators reflecting changes in climate conditions [2]. The IPCC Fifth Assessment Report (2013) [1] also indicated that rainfall is expected to increase in the future, especially extreme rainfall patterns.

Among many physical factors influencing river flooding, precipitation is the main factor affecting runoff formation and flow regimes [2,3]. Additionally, other factors caused by human interactions such as land-use changes and rapid development of socio-economic activities also promote floods [4]. Due to the complex nature of these factors on river flooding, all of them are considered a challenging task for hydrologists. However, evaluating flood risk by considering the main factors affecting river floods is possible for each river basin’s characteristics.

The risk of severe natural events, such as flood inundation, landslides, and sediment-related disasters caused by extreme rainfall, is likely to increase in the future [5,6]. During summer, severe flood risk from extreme rainfall events is a common phenomenon, particularly throughout Asia [7]. Moreover, large-scale damage caused by flood inundation has been reported around the world in recent years [8,9]. Hence, many attempts have been made to evaluate and predict flood inundation risk along many important river basins in the world to minimize damage due to loss of life and property [4,10,11,12,13,14,15,16]. For example, evaluating the river flood risk in small and ungauged river basins is currently one of the most common challenges for hydrologists. Hence, several studies have been conducted to investigate and map the flood hazards for these river basins to minimize flood damage [13,14,15].

In Japan, a new record has been reached for temperature and precipitation almost every year in recent times [17]. According to the Japan Meteorological Agency (2017) [18], the annual mean surface temperature over Japan has increased by 1.19 from 1898 to 2017, which is higher than the global average surface temperature increase (0.73 °C). Additionally, the annual number of days when precipitation is equal to or greater than 200 mm has increased (statistically significant at a confidence level of 99%) during the period from 1901 to 2017 in Japan [18]. Flood inundation disasters due to extreme rainfall have occurred frequently in Japan and caused considerable damage in recent years. For example, an extreme rainfall event on 5–6 July 2017 caused severe floods in many parts of Kyushu, Japan, and resulted in extensive damage [19]. One year later, during 5–8 July 2018, widespread flood disaster caused by torrential rainfall occurred in western Japan. The Oda River and its three tributaries exhibited the “backwater phenomenon”, where the tributary river floods intersected with the main river floods causing levees to break at eight points [20]. This backwater phenomenon is projected to occur anywhere across Japan because of the presence of several river systems with similar conditions. In the northernmost Japanese island, Hokkaido, extreme rainfall caused severe flood inundation in August 2016; this disaster reportedly caused approximately USD 260 million in damage, in addition to agricultural losses over 40,258 ha of land throughout Hokkaido [21]. Most recently, on 1–13 October 2019, severe flooding and landslides occurred in many river basins in the central–northern parts of Japan, causing great losses to life and property [22]. In this context, to mitigate water-related disasters due to climate changes, especially flood inundation disasters, many researchers have attempted to predict extreme rainfall and severe river flooding events along many important river basins in Japan [23,24,25,26]. These attempts could provide useful information to policymakers for climate change adaptations as well as flood mitigation strategies.

Prediction of changes in the flow regime is a challenge [10,27] as it demands future climate datasets under various climate change scenarios. Popular methods of future projection are the general circulation models (GCMs), whose rainfall and evapotranspiration data could be used as inputs to hydrological models to assess changes in river floods. However, GCMs generally have a coarse resolution for simulation on a regional scale; such coarse resolution could lead to large uncertainties in the simulations [10,28]. Therefore, high-resolution climate data obtained through downscaling approaches, such as statistical and dynamical downscaling, have been investigated in recent years. For example, in a study by Sato (2012) [29], the impact of climate change on river discharge was investigated for several major river basins in Japan using the latest version of an atmospheric GCM (AGCM) with a 20 km horizontal resolution. In a study by Yamada et al. (2018) [30], the super-high-resolution 5 km d4PDF dataset was downscaled from the 20 km resolution data via a nonhydrostatic regional climate model [31]. After downscaling, the rainfall amount and spatiotemporal distribution of rainfall were similar to the observed rainfall events. Furthermore, heavy rainfall patterns could be better observed after downscaling, and the 5 km resolution results represented the topography of the study area more accurately [23,30]. This suggests that the 5 km downscaled results may be appropriate for estimating changes in heavy rainfall as well as flood inundation on the basin scales.

The hydrological model is a fundamental approach to further understand the process of river flow. In recent years, numerous hydrological models have been developed and applied to simulate and predict potential river floods along many river basins worldwide, such as the HEC-HMS [32], IFAS [33], TOPMODEL [34], and SWAT models [35]. However, these models cannot simulate flood inundation along the river basin and must be combined with hydraulic models, such as HEC-RAS [36], FLO-2D [37], and MIKE 21 [38], which are popular. The 1D models cannot simulate the inundation risk in floodplains, while 2D flood models require a long simulation time. Therefore, both these types of models were integrated for comprehensive flood inundation simulations and reduced computation time. Some popular integrated 1D–2D models are SOBEK [39] and the Rainfall–Runoff–Inundation (RRI) model [40,41,42]. In a study by San et al. (2020) [43], the performances of RRI and SOBEK models were compared in terms of user-friendliness, cost, type of output, and correlation between simulated and observed data in the Bago River basin. The results indicated that although the SOBEK model could provide more accurate results, the RRI model is preferable considering the calculation time, cost, and user interface friendliness. In addition, the RRI model is more applicable to the near-real-time inundation mapping and the flood early warning system. Therefore, due to these advantages, flood inundation risk was assessed using the RRI model in this study. It can simulate rainfall–runoff and inundation processes simultaneously in both the floodplain and mountainous zones [40,41,42]. Several studies have used this model to estimate flood inundation for many river basins globally and have demonstrated its good performance [10,11,44,45,46,47,48].

In a study by Nguyen [23], changes in extreme rainfall and changes in river discharges were evaluated in the Ishikari River basin (IRB), Japan, as well as its sub-basins using the IFAS model coupled with a large-ensemble rainfall dataset (d4PDF) with 5 km resolution. The results indicated that extreme rainfall was projected to increase by 21.2% in the IRB and 34.8% in the Chitose River basin (CRB), whereas the river discharge could increase drastically by 33.3% in the IRB and 37% in the CRB. The effect of climate change is significant in the CRB, the downstream lowland area of the IRB. This river basin is considered to be the most flood-prone in the IRB. In addition, the time difference between the peak discharge at the reference station in each tributary and the peak water level at the confluence point for the seven Ishikari River main sub-basins using the IFAS model was also predicted. However, since the IFAS model uses a kinematic wave for river channel flow simulation [33], the prediction might not accurately reflect the predicted results.

Considering the aforementioned information, this study aims to assess the change in flood inundation risk comprehensively (including inundation area, inundation volume, peak inundation depth, inundation frequency, and inundation duration) considering the extreme rainfall impacts and topographic vulnerability in the CRB using the RRI model and rainfall data extracted from d4PDF with a high-resolution of 5 km. Owing to the topographical characteristics of its low-lying area, the CRB is frequently affected by backwater from the Ishikari River and experiences severe flood inundation. Moreover, this study is conducted to improve the time difference prediction result considering the dangerous impacts of backwater phenomenon on flood inundation risk in the CRB using the RRI model, which uses a 2D distributed hydrodynamic model with diffusive wave models. The time difference is evaluated between the peak discharge at the reference station in the CRB and the peak water level at the confluence point intersecting the main river. With this effort, we believe that the result might provide more reliability compared to previous studies of the CRB. Flood inundation is more severe when the time difference is shorter. The study results are expected to provide useful information for developing climate change adaptation measurements and mitigating future flood inundation damage in the floodplain areas.

2. Materials and Methods

2.1. Description of the Study Area

The study area considered herein is the IRB, which is located in Hokkaido, Japan (Figure 1). At a length of 268 km, this is the longest river in Hokkaido and the third longest river in Japan. In addition, this river has the second largest basin area in Japan, with a drainage area of 14,330 km². The river originates from Mt. Ishikari-dake (1967 m above sea level), passes through the west of Hokkaido, and empties into the Sea of Japan. It flows through 48 municipalities (including Sapporo, the prefectural capital, and Asahikawa, the second-largest city), accounting for 52% of the population of Hokkaido. The Ishikari Plain occupies most of the basin’s area and is located around the central and downstream basin area, which is the most productive agricultural area in Hokkaido. The IRB plays a vital role in the socioeconomic development of Hokkaido. The mean annual precipitation in this basin is 1300 mm; the hydrologic peaks generally occur from March to May during the snow-melt period and from August to September during the rainy season [49]. The CRB is one of the main tributaries of Ishikari River (Figure 1). Owing to the topographical characteristics of its low-lying area, the CRB is frequently affected by backwater from the main river, causing severe flood inundation. Therefore, the CRB is considered a vulnerable region in the IRB. The Chitose River has a length of 108 km and a total basin area of 1244 km². The main parts of the New Chitose Airport belong to this basin area [50].

In the past, the IRB has experienced large-scale historical flood events and severe damage. For example, the largest flood event in this basin occurred in August 1981, and severe damage was recorded in the target basin. Figure 2 shows the flood inundation from August 1981 in the CRB. The backwater phenomenon was observed when the water levels at the confluence point surged higher and water from the main river flowed back into the tributary, which combined with inland water and caused severe flood inundation. This flood event caused damages amounting to approximately JPY 115.2 billion, two fatalities, and 22,500 inundated houses throughout the IRB. For the CRB, 192 km² of land was inundated and 2700 inundated houses were recorded for this serious flood event [49,50].

2.2. Overview of the Rainfall–Runoff–Inundation (RRI) Model

In this study, the RRI model was employed to simulate flood inundation and identify flood-prone areas. The RRI model was developed by the International Center for Water Hazard and Risk Management (ICHARM). The RRI model is a two-dimensional model that is capable of simulating rainfall–runoff and flood inundation simultaneously [40,41,42]. At a grid cell in which a river channel is located, the model assumes that both slope and river are positioned within the same grid cell. The channel is discretized as a single line along its centerline of the overlying slope grid cell. The flow on the slope grid cells was calculated with the 2D diffusive wave model, whereas the channel flow was calculated with the 1D diffusive wave model. The RRI model simulates the lateral subsurface, vertical infiltration, and surface flows to better represent flood characteristics. The lateral subsurface flow, which is more important in mountainous regions, was treated in terms of the discharge–hydraulic gradient relationship, including both saturated subsurface and surface flows, whereas the vertical infiltration flow was evaluated using the Green–Ampt model. The details of the model are shown below [41].

A storage cell-based inundation model [51] was used to calculate lateral flows on slope grid cells. The governing equations for the 2-dimensional unsteady flow model include the mass balance equation (Equation (1)) and momentum equations (Equations (2) and (3)):

$$\frac{\partial h}{\partial t} + \frac{\partial q_x}{\partial x} + \frac{\partial q_y}{\partial y} =r-f$$

(1)

$$\frac{\partial q_x}{\partial t} + \frac{\partial u{q}_x}{\partial x} + \frac{\partial v{q}_x}{\partial y} = -gh\frac{\partial H}{\partial x} - \frac{{\tau }_x}{{\rho }_{\omega }}$$

(2)

$$\frac{\partial q_y}{\partial t} + \frac{\partial u{q}_y}{\partial x} + \frac{\partial v{q}_y}{\partial y} = -gh\frac{\partial H}{\partial y} - \frac{{\tau }_y}{{\rho }_{\omega }}$$

(3)

where h is the water level from the local surface; $q_x$ and $q_y$ are the unit width discharges in the x and y directions, respectively; $u$ and $v$ are the flow velocities in the x and y directions, respectively; $r$ is the rainfall intensity; $f$ is the infiltration rate; H is the water level from the datum; $p_{\omega }$ is the density of water; $g$ is the gravitational acceleration; and ${\tau }_x$ and ${\tau }_y$ are the shear stresses in the x and y directions, respectively.

The second terms of the right side of Equations (2) and (3) were calculated with Manning’s equation.

$$\frac{{\tau }_x}{{\rho }_{\omega }} = \frac{g{n}^{2}u\sqrt{u^2+v^2}}{h^{1/3}}$$

(4)

$$\frac{{\tau }_y}{{\rho }_{\omega }} = \frac{g{n}^{2}v\sqrt{u^2+v^2}}{h^{1/3}}$$

(5)

where n is the Manning’s roughness parameter.

The RRI model spatially discretized the mass balance (Equation (1)) as follows:

$$\frac{d{h}^{i,j}}{dt} + \frac{q_x^{i,j-1}- q_x^{i,j}}{\Delta x} + \frac{q_y^{i-1,j}- q_y^{i,j}}{\Delta y} = r^{i,j} - f^{i,j}$$

(6)

where $q_x^{i,j}$, $q_y^{i,j}$ are the x and y direction discharges from a grid cell at ($i, j$), respectively.

Equations (7) and (9) describe the saturated subsurface flow based on the Darcy law, while Equations (8) and (10) describe the combination of the saturated subsurface flow and the surface flow. For the kinematic wave model, the hydraulic gradient was assumed to be equal to the topographic slope, whereas the RRI model assumed the water surface slope as the hydraulic gradient.

$$q_x = -k_{a}h\frac{\partial H}{\partial x}, \left(h\le d_a\right)$$

(7)

$$q_x = -\frac{1}{n}{\left(h-d_a\right)}^{\frac{5}{3}}\sqrt{\left|\frac{\partial H}{\partial x}\right|}sgn\left(\frac{\partial H}{\partial x}\right) - k_{a}h\frac{\partial H}{\partial x}, (d_a<h)$$

(8)

$$q_y = -k_{a}h\frac{\partial H}{\partial y}, \left(h\le d_a\right)$$

(9)

$$q_y = -\frac{1}{n}{\left(h-d_a\right)}^{\frac{5}{3}}\sqrt{\left|\frac{\partial H}{\partial y}\right|}sgn\left(\frac{\partial H}{\partial y}\right)-k_{a}h\frac{\partial H}{\partial y}, (d_a<h)$$

(10)

where $k_a$ is the lateral saturated hydraulic conductivity and $d_a$ is the effective porosity of the soil depth times.

Equations (11) and (12) were used to simulate the effect of unsaturated subsurface flows, saturated subsurface flows, and the surface flow with the single variable of $h$.

$$q_x=\left\{\begin{array}{c}-k_m{d}_m{\left(\frac{h}{d_m}\right)}^{\beta }\frac{\partial H}{\partial x}, \left(h\le d_m\right)\\ -k_a\left(h-d_m\right)\frac{\partial H}{\partial x}- k_m{d}_m\frac{\partial H}{\partial x}, (d_m<h\le d_a)\\ -\frac{1}{n}{\left(h-d_a\right)}^{\frac{5}{3}}\sqrt{\left|\frac{\partial H}{\partial x}\right|}sgn\left(\frac{\partial H}{\partial x}\right)-k_a\left(h-d_m\right)\frac{\partial H}{\partial x}-k_m{d}_m\frac{\partial H}{\partial x}, (d_a<h)\end{array}\right.$$

(11)

$$q_y=\left\{\begin{array}{c}-k_m{d}_m{\left(\frac{h}{d_m}\right)}^{\beta }\frac{\partial H}{\partial y}, \left(h\le d_m\right)\\ -k_a\left(h-d_m\right)\frac{\partial H}{\partial y}- k_m{d}_m\frac{\partial H}{\partial y}, (d_m<h\le d_a)\\ -\frac{1}{n}{\left(h-d_a\right)}^{\frac{5}{3}}\sqrt{\left|\frac{\partial H}{\partial y}\right|}sgn\left(\frac{\partial H}{\partial y}\right)-k_a\left(h-d_m\right)\frac{\partial H}{\partial y}-k_m{d}_m\frac{\partial H}{\partial y}, (d_a<h)\end{array}\right.$$

(12)

Infiltration is an important hydrological process. For relatively flat areas, the vertical infiltration process during the first period of rainfall has more impact of large-scale flooding. Therefore, the vertical infiltration can be treated as loss for event-based simulation. In RRI, infiltration loss $f$ is calculated with the Green–Ampt infiltration model. The Green–Ampt infiltration model is a simplified physical model and based on the Richard equation. It relates the rate of infiltration to measurable soil properties such as the porosity, hydraulic conductivity, and soil water content of a particular soil column [41,52].

$$f= k_v⌈1+\frac{\left(\varnothing - {\theta }_i\right)S_f}{F}⌉$$

(13)

where $k_v$ is the vertical saturated hydraulic conductivity; $\varnothing$ is the soil porosity; ${\theta }_i$ is the initial water volume content; $S_f$ is the suction at the vertical wetting front and $F$ is the cumulative infiltration depth.

The flow interactions between the river channel and slope were estimated using different overflow formulas [41].

(1)

Flow from slope to river under the normal condition: when the river water level is lower than the ground level, the discharge from the slope into the river is defined by Equation (14):

$$q_{sr}={\mu }_1{h}_s\sqrt{g{h}_s}$$

(14)

where $q_{sr}$ is discharge from slope to river; ${\mu }_1$ is the constant coefficient (=0.544), and $h_s$ is the water depth on a slope cell.

(2)

No water exchange between slope and river: when the river water level is higher than the ground level and both the river and slope water levels are lower than the levee height.

(3)

Overtopping flow from river to slope: when the river water level is higher than the levee height and the slope water level:

$$q_{rs}={\mu }_2{h}_1\sqrt{2g{h}_1}$$

(15)

where $q_{rs}$ is discharge from river to slope; ${\mu }_2$ is the constant coefficient (=0.35), and $h_1$ is the difference between the river water level and the levee height.

(4)

Overtopping flow from slope to river: when the slope water level is higher than the levee height and river water level, and the river water level is higher than the levee height. The flow exchange was calculated by Equation (15), with $h_1$ presenting the difference between the two water levels and $q_{rs}$ being replaced with $q_{sr}$.

To solve the diffusive equations, the RRI model uses the fifth-order Runge–Kutta method with adaptive time-step control. A detailed explanation is given by Sayama et al. (2012) [41]. The RRI model provides the output of river discharge, river water level, inundation area, and inundation depth simultaneously. A schematic of this model is shown in Figure 3. This model has been applied to several important river basins to evaluate river floods as well as flood inundation, and good simulation results have been obtained [10,11,44,45,46,47,48].

2.3. Data Inputs

2.3.1. Topographic Data

The topography data, including digital elevation model (DEM), flow direction (DIR), and flow accumulation (ACC), were obtained from the HydroSHEDS project provided by the U.S. Geological Survey (USGS) [53]. HydroSHEDS is derived primarily from elevation data of the Shuttle Radar Topography Mission (SRTM) data. HydroSHEDS offers a suite of geo-referenced datasets in raster and vector format, including stream networks, watershed boundaries, drainage directions. The HydroSHEDS supports regional and global watershed analyses, hydrological modelling, and freshwater conservation planning. Available resolutions of HydroSHEDS are 5 min, 30 arc-seconds, 15 arc-seconds, and 3 arc-seconds (approximately 10 km, 1 km, 500 m, and 90 m at the equator, respectively). The topography of the HydroSHED data is available at https://www.hydrosheds.org/ (accessed on 1 February 2021). The higher DEM resolution could provide a finer resolution of inundation mapping. However, this requires greater computational power compared to the lower DEM resolution. The HydroSHEDS data at 3 arc-seconds (approximately 90 m) were used in this study. However, to reduce computational time due to the large study area, a topography data scale-up algorithm was used in the RRI model to convert the original topography data to 24 arc-seconds (approximately 720 m) for the entire IRB, and to 9 arc-seconds (approximately 270 m) for the CRB. The DEM, DIR, and ACC data were converted to ASCII format, which is the required format for the RRI model.

2.3.2. Land Cover Data

Global Land Cover Characterization (GLCC-V2) data with 1 km resolution provided by the USGS were used in this study. Land cover has a crucial impact on the hydrological process, and the original land cover data for the target basin has 14 distinct land cover types. However, these 14 types are too detailed to assign different parameters in the RRI model; therefore, similar land cover types were merged and three land cover types, namely forest, sparse vegetation, and cropland/floodplain, were finalized using the ArcGIS function. The reclassified land cover data were resampled to the same resolution as the topography data to input the RRI model.

2.3.3. River Cross-Section Data

The river cross-section (river width and depth) is an important parameter to input the RRI model to minimize uncertainties in the output of the hydrologic simulation. In the RRI model, a rectangular cross-section is applied to the cross-section of the river channel. The river width (W) and river depth (D) were approximated by the following equations:

$$W=C_W{A}^{S_W}$$

(16)

$$D=C_D{A}^{S_D}$$

(17)

where A is the upstream contributing area (km²), W is the channel width (m), and D is the channel depth (m). C_w, S_w, C_D, and S_D are geometrical parameters. The parameters of the equations were estimated by regression analysis with the cross-section data at 42 locations spanning the IRB. In addition, the parameters of the equations were estimated at 33 locations across the CRB. The obtained parameters for the Ishikari River cross-section are C_w = 1.73, S_w = 0.60, C_D = 0.18, S_D = 0.45 and the obtained parameters for the Chitose River cross-section are C_w = 14.45, S_w = 0.34, C_D = 0.022, and S_D = 0.80.

In addition, to increase the reliability of future flood inundation hazard assessments under the various climate change scenarios in the vulnerable areas, the Ishikari River and Chitose River cross-sections of future plans were provided by the River Management Division of the Sapporo Development and Construction Department, Hokkaido Regional Development Bureau and future flood inundation simulations were applied in this study. The obtained parameters for the Ishikari River cross-section plan are: C_w = 1.68, S_w = 0.60, C_D = 0.17, and S_D = 0.46. The Chitose River cross-section future plan remained unchanged compared to the present river cross-section; therefore, the Chitose River cross-section future plan parameters are C_w = 14.45, S_w = 0.34, C_D = 0.022, and S_D = 0.80.

2.3.4. Gauged Rainfall and Runoff Data

In the past, the IRB experienced severe large-scale historical flood events. In particular, the flood event of August 1981 was the largest in the target basin. Notably, during this flood event, the CRB was severely affected—backwater phenomenon integrated inland flooding was observed and caused tremendous damage. Therefore, in this study, we used this large-scale flood event to calibrate the RRI model. Additionally, three severe historical flood events in September 2001, September 2011, and August 2016 were chosen to validate the model.

To evaluate the performance of the RRI model, a comparison between the observed and simulated river discharges was employed at the Ishikari Ohashi station, approximately 26.6 km upstream from the river mouth (Figure 1). In addition, the inundated area in the CRB was evaluated and compared with the observed inundated area for the flood event from August 1981 to calibrate the RRI model.

The gauged rainfall data and observed river discharge for the August 1981 flood event were provided by the Hokkaido Regional Development Bureau (HRDB). The observed inundated area in the CRB for the August 1981 flood event was also provided by the HRDB. The gauged rainfall for the September 2001, September 2011, and August 2016 flood events and the observed river discharges for the validation process were obtained from the Water Information System managed by the Ministry of Land, Infrastructure, and Transport, Japan (MLIT) [54].

2.3.5. Evaluation of RRI Model Performance

The RRI model calibration and validation process performances were evaluated using four statistical indicators, including the Nash–Sutcliffe coefficient (NSE), coefficient of determination (R²), relative volume error (VE), and peak discharge error (E_p). In addition to the common statistical indicators NSE and R², the VE and E_p statistical indicators were defined as follows:

$$VE= \frac{{{\displaystyle \sum }}_{i=1}^n{Q}_{Si\left(i\right)}-{{\displaystyle \sum }}_{i=1}^n{Q}_{Ob\left(i\right)}}{{{\displaystyle \sum }}_{i=1}^n{Q}_{Ob\left(i\right)}}$$

(18)

$$E_p = \frac{Q_{MOb}-Q_{MSi}}{Q_{MOb}}$$

(19)

where $Q_{Si}$ is the simulated discharge (m³/s); $Q_{Ob}$ is the observed discharge (m³/s); $Q_{MOb}$ is the peak observed discharge (m³/s); $Q_{MSi}$ is the peak simulated discharge (m³/s); and n is the total number of measurements.

2.4. d4PDF Dataset

This study used the large-ensemble rainfall dataset from the database of Policy Decision making for Future climate change (d4PDF) [55], which consists of large-ensemble climate experiments from the global atmospheric model, with a horizontal resolution of 60 km, and regional downscaling simulations covering the Japan area, with a horizontal resolution of 20 km. The d4PDF consists of historical climate simulations (1951–2010) with 50 ensembles (historical simulation: 50 ensembles × 60 years = 3000 events), and the +4 K future climate simulations (2051–2110) with 90 ensembles, including six sea surface temperature (SST) patterns and 15 ensembles for each SST (+4K future simulation: 6 SST × 15 ensembles × 60 years = 5400 events). The +4 K future climate was predicted for a global mean air temperature 4 °C warmer than the preindustrial period. The +4 K future climate simulation corresponded to the representative concentration pathway 8.5 (RCP8.5) experiments under phase 5 of the CMIP5 [56]. This study used the d4PDF rainfall data with a 5 km resolution, which were downscaled from the 20 km resolution data [23]. The target rainfall data were defined from 1 June to 1 December, where the annual maximum rainfall (mm/72 h) was obtained over the target basin for each of the 3000 historical and 5400 future simulation events. The rainfall data with a high resolution of 5 km reflect the topography of the study area more accurately [23,30]. In this study, we used the top 20 and top 36 rainfall events from the previous study of Nguyen et al. (2020) [23] to evaluate changes in flood inundation risk comprehensively in the CRB using the RRI model.

2.5. Analytical Procedure

This study was conducted to assess changes in severe river flooding in the IRB as well as changes in flood inundation in the CRB using the top 20 and top 36 rainfall events from the d4PDF dataset [23] as the input rainfall data to the RRI model. This study was performed according to the following steps. First, the RRI model was calibrated to a large historical flood event in August 1981 using the “trial and error” method. Then, the tuned parameters were used for the validation process of the flood events in September 2001, September 2011, and August 2016. The optimal parameters after calibration and validation were applied to the RRI model to evaluate future changes in severe river discharges in the IRB and flood inundation in the CRB. To evaluate the flood inundation in the CRB, the time series of water levels at the confluence point (Figure 1) was used as a boundary condition in the RRI model. The simulations were conducted for an extra period of two weeks before the main period for model warmup and also for initial moisture conditions in the target basin. Changes in flood inundation in the CRB were investigated comprehensively in terms of the inundation area, inundation volume, peak inundation depth, inundation frequency, and inundation duration. The statistical Kolmogorov–Smirnov (K–S) test, which is a nonparametric test of two samples, was used to examine the variations in flood inundation between the top 20 rainfall events (historical simulation) and top 36 rainfall events (future simulation). The value of the test statistic for the two samples is defined by

$$D=sup\left|F_1\left(x\right)-F_2\left(x\right)\right|$$

(20)

where F₁ and F₂ are the empirical distribution functions based on the two samples, and sup is the supremum function. The hypothesis of the test is as follows:

H₀: 

The two samples have no significant differences in their cumulative distribution functions.

H₁: 

Two samples have different trends in the cumulative distribution functions.

The null hypothesis is rejected if the p-value is less than the significance level α of 0.05.

Finally, the time difference prediction results between peak discharge at the reference station in the CRB and the peak water level at the confluence point intersecting the Ishikari River were improved compared to the previous study using the RRI model.

3. Results

3.1. RRI Model Calibration and Validation

The RRI model was calibrated for the largest flood event in the IRB in August 1981. The optimal parameters of the RRI model after calibration are listed in

Table 1. Optimal parameter settings of the RRI model.

Parameters	Notation	Unit	Mountains	Plains
Manning’s coefficient for river	nriver	(m−1/3s)	0.035	0.035
Manning’s coefficient for slope	n	(m−1/3s)	0.15	0.15
Soil depth	d	(m)	1.0	1.0
Soil porosity	Φa	(-)	0.12	0.472
Vertical hydraulic conductivity	kv	(cm/h)	-	0.06
Wetting front soil suction head	Sf	(-)	-	0.273
Lateral saturated hydraulic conductivity	ka	(m/s)	0.17	-
Unsaturated porosity	Φm	(-)	0.03	-
Parameter of unsaturated hydraulic conductivity	β	(-)	8.0	-

. The simulation results are shown in Figure 4. As shown in Figure 4, the RRI model closely reproduced the wave shape and flood duration in most cases. Figure 4a compares the observed and simulated discharges at the Ishikari Ohashi station for the flood event of August 1981 after calibration. The simulated discharges show close agreement with the observed discharges. This result is determined by reasonable statistical indicators NSE, R², and VE with values of 0.96, 0.98, and −0.006, respectively. The observed peak discharge (11,330 m³/s) has good synchronization with the simulated one (11,310 m³/s), as indicated by an E_p value of 0.002. The optimal parameters after the calibration process were used for the validation process of the flood events in September 2001, September 2011, and August 2016, as shown in Figure 4b–d, respectively. It can be seen that the simulated discharges matched well with the observed discharge in most cases, as indicated by the reasonable NSE, R², and VE values. However, the simulated peak discharges were found to be slightly higher than the observed peak discharges for these flood events. This result is demonstrated by the E_p values of −0.234, −0.1, and −0.17 for the flood events in September 2001, September 2011, and August 2016, respectively. The detailed statistical indicators are listed in

Table 2. Statistical indicators of the RRI model performance.

	Flood Events	NSE	R2	Ep	VE
Calibration	August 1981	0.96	0.98	0.002	−0.006
Validation	September 2001	0.81	0.99	−0.234	0.246
	September 2011	0.88	0.91	−0.10	0.0007
	August 2016	0.78	0.89	−0.17	−0.002

. Overall, these results suggest that the RRI model can perform reasonably well for further analysis in the next sections.

Figure 5 shows the simulated maximum flood inundation depth and inundation areas for the August 1981 flood event in the CRB. As shown in Figure 5, the simulated inundation area defined with a water depth threshold of 0.5 m overlaps with the observed inundation area for the flood event in August 1981 in the downstream area of the CRB. Owing to the limitations of the observed inundation data, inundation depth data were not available. A comparison between the simulated and observed inundation areas was conducted using the basic photo provided by the Hokkaido Regional Development Bureau. The observed flood depth data were not available; therefore, in this study, we selected the threshold of water depth is 0.2 and 0.5 m to classify the inundated and noninundated areas. According to the simulation results, the simulated inundation areas are 252 and 240 km² for the 0.2 and 0.5 m water depth thresholds, respectively. Additionally, the observed inundation area was 192 km² [50]. The simulated inundation area for the threshold of 0.5 m water depth was more similar to the observed ones compared with the simulated inundation area for the threshold of 0.2 m water depth. Additionally, the simulated inundation area and peak inundation depth for the August 1981 flood event are quite close to the simulated results estimated by Oki et al. (2017) [50] of 209 km² for the inundated area. Furthermore, the peak inundation depth is similar to the value simulated by Oki et al. (2017) [50]. As a result, the threshold of 0.5 m water depth was chosen to classify the inundated and noninundated areas in the CRB in the following sections of this study.

3.2. Evaluation of Climate Change Impacts on River Flood in the IRB

The optimal parameters of the RRI model after calibration and validation were used to evaluate the changes in severe river floods in the IRB, as well as the flood inundation in the CRB using the top 20 and top 36 rainfall events extracted from the d4PDF rainfall dataset as input data to the RRI model. The simulations were conducted for an additional period of 72 h before the main period as an initial condition for model warmup.

Changes in severe river flooding were examined at three stations along the main river of the IRB. The locations of the three stations are shown in Figure 1, including the Osamunai station (located in the upstream river), Naieoohashi station (located in the middle river), and Ishikari Ohashi station (located in the downstream river). From the simulation results at three stations along the main river, the percentage differences in the means of the peak discharge between the top 20 and top 36 rainfall events indicated that severe river floods are expected to increase by 21–24% in the target basin. These results indicate that severe river flooding is projected to increase by a similar extent throughout the entire IRB. The greatest increase of 24% is observed at the Naieoohashi station, which is located in the middle of the Ishikari River.

3.3. Evaluation of Climate Change Impacts on Flood Inundation in the CRB

Future changes in severe flood inundation were evaluated for the CRB. To further understand the changes in the flood inundation risk in the CRB, the flood inundation area and flood inundation volume were evaluated for an inundation depth equal to or larger than 0.5 and 3.0 m in this study. Figure 6 shows the inundation area averaged for the top 20 rainfall events (historical simulation) and top 36 rainfall events (future simulation) for flood depths equal to or larger than 0.5 (Figure 6a) and 3.0 m (Figure 6b).

Table 3. Percentage difference of the mean values of the inundation area and inundation volume for inundation depths of ≥0.5 and ≥3.0 m between the top 20 rainfall events (historical simulation) and top 36 rainfall events (future simulation).

	Inundation Area	Inundation Volume
Inundation depth ≥0.5 m	24.5	46.5
Inundation depth ≥3.0 m	124.6	134.2

shows the percentage difference of the mean values of the inundation areas and inundation volumes for an inundation depth ≥0.5 m and inundation depth ≥3.0 m between the top 20 rainfall events (historical simulation) and the top 36 rainfall events (future simulation). The results indicate that flood magnitudes are projected to increase significantly in the future in the CRB. For flood depths ≥0.5 m, the inundation area and inundation volume are expected to increase by 24.5 and 46.5%, respectively. As shown in Figure 6a, the downstream area of the CRB is the area with more severe flood inundation in the future. In addition, flood inundation was observed on the main sub-basins of the CRB. Severe flood inundation events for flood depths ≥3.0 m are expected to increase drastically in the CRB by 124.6 and 134.2% for the inundation area and inundation volume, respectively. In particular, an inundation depth equal to or larger than 3.0 m was observed only in the downstream area of the target basin. Owing to the low-lying characteristics, this area is the most vulnerable area in the CRB. The results of the K–S tests of the inundation area with inundation depths equal to or larger than 0.5 and 3.0 m demonstrate significant differences between the top 20 and top 36 rainfall events, with p-values of 0.0003 and 0.0001 (i.e., much smaller than alpha of 0.05), respectively. A similar K–S test result for the inundation volume with the inundation depths equal to or greater than 0.5 and 3.0 m also revealed remarkable differences between the top 20 and top 36 rainfall events, with p-values of 0.0006 and 0.0001, respectively.

Figure 7 shows the boxplots of the peak inundation depth (Figure 7a), inundation volumes with an inundation depth ≥0.5 m (Figure 7b), and inundation volumes with an inundation depth ≥3.0 m (Figure 7c) for the top 20 past and top 36 future rainfall events in the CRB. In addition to the significant increase in inundation area and inundation volume, the peak inundation depth is predicted to significantly increase by a factor of 13.8% in the future. Notably, the peak inundation depth can reach 8.45 m in the future in the CRB. This inundation depth is a very serious level of classification based on the risk categories for inundation depths in Japan. In addition, there was a huge variation in peak inundation depth in the CRB. The first quartile value of the CRK peak inundation depth for the future simulation is 6.02 m, which was greater than the third quartile value for the historical simulation (6.0 m). Similarly, the inundation volume with the flood depth $\ge$0.5 m for the future simulation is 0.55 km³, which was equal to the third quartile value from the historical simulation. For the flood depth $\ge$3.0 m, the first quartile value for the inundation volume is 0.13 km³, which is greater than the third quartile value from the historical simulation of 0.12 km³. In addition, the median inundation volume with a flood depth $\ge$0.5 m is expected to increase significantly by 52.9% in the CRB, whereas a drastic increment of 330% was projected in the CRB for the flood depth $\ge$3.0 m. The K–S test also shows a significant difference for the peak inundation depth between the top 20 past and top 36 future rainfall events at a significance level of 5% (p-value = 0.0001). A similar K–S test result for the inundation volume with inundation depths equal to or greater than 0.5 and 3.0 m also revealed remarkable differences between the top 20 and top 36 rainfall events.

Figure 8 shows the spatial distribution of the inundation frequency for the top 20 past and top 36 future rainfall events, and the inundation frequency difference (IFD) between these rainfall events. As shown in Figure 8a,b, the downstream area of the CRB is clearly seen to be the most severe flood inundation. This result is suitable for the terrain characteristics of the CRB. It can be seen that flood inundation is likely to occur in this area in the future with a frequency of 90–100%. Moreover, the flood inundation was observed in the tributaries of the Chitose River for both historical and future simulations. This would increase the risk of flood inundation in the downstream area of the CRB, in which the inland flooding integrated to the backwater phenomenon from the main river. In addition, in Figure 8c, the IFD indicates a positive value for almost all areas in the CRB. The proportion for the area with an inundation frequency smaller than 0 is only 4.2%. Therefore, 95.8% of the area had an inundation frequency larger than 0, in which the area proportional to 0 < IFD ≤ 0.3 is 74.7%, and the area proportional to 0.3 ≤ IFD ≤ 0.6 is 21.1%. On the other hand, this result clearly indicates that the flood-prone areas downstream of the CRB will be more widely spread in the future. Hence, this area should be specially considered in the future to mitigate severe flood inundation damage.

In addition to evaluating the inundation area, inundation volume, peak inundation depth, and inundation frequency, the future changes in inundation duration were also investigated for the top 20 past and top 36 future rainfall events in this study. Figure 9 shows the spatial distribution of the inundation duration for the top 20 past and top 36 future rainfall events as well as for the difference between them. This study focused on assessing the changes in extreme flood inundation during the short-term period (72 h). Each event from the top 20 past and top 36 future rainfall events is the annual maximum rainfall (mm/72 h) [23]. It can be seen that the longer durations of inundation were observed in the downstream area of the CRB for both historical and future simulations. The inundation in most tributaries lasts for a short duration of 0.03 to 10 h, whereas the inundation duration in the downstream area of the CRB lasts from 20 to 40 h (approximately 1–2 days) for these rainfall events. Overall, the increment of flood inundation was observed at the tributaries and floodplain areas in the CRB. The proportion of the area with an inundation duration difference of more than 0 h is 99.7% in the entire CRB, and only 0.3% of the areas will remain unchanged in terms of inundation duration in the future. The results indicate that the inundation duration is expected to increase by approximately 5–10 h in vulnerable areas in the CRB, as shown in Figure 9c. Additionally, the results reveal that the area with high inundation frequency corresponds to the long inundation duration and vice versa in the CRB.

3.4. Evaluation of Time Difference Impacts on Flood Inundation in the CRB

Owing to the topographical characteristics of its low-lying area, the CRB is frequently affected by backwater from the main river. The CRB has experienced severe flood inundation events caused by extreme short-term rainfall events. The backwater phenomenon integrated inland flooding, which has caused serious damage in the past, was observed. When the water level increases at the confluence point, the backwater phenomenon can occur in this basin. Furthermore, flood inundation is more severe when the time difference is shorter. Therefore, considering the dangerous impacts of the backwater phenomenon on flood inundation risk in the CRB, the differences between the time of peak discharge at the reference station in the CRB and that of peak water level at the confluence point intersecting the Ishikari River was investigated using the RRI model. Figure 10 shows the frequency of the time differences in the CRB for the top 20 past and top 36 future rainfall events. The results reveal that the time difference has a large range from 0 to 50 h in the future, whereas the time difference from the historical simulation has a range from 0 to 30 h in the CRB. However, it can be observed that the short time difference from 0 to 10 h is projected to significantly increase in the future. This result indicates that the flood inundation risk due to the backwater phenomenon is likely to be more severe in the future in the CRB. In addition, the results indicate that the time difference from 10 to 20 h does not change in the future, whereas the time difference from 20 to 30 h is predicted to increase remarkably in the future. These useful results will help policymakers enact proactive measurements to mitigate flood inundation damages in vulnerable areas in the target basin.

4. Discussion

This study investigated the impact of climate change on potential river flood hazards along the IRB, and in particular flood inundation risk in the CRB, a tributary of the Ishikari River, using a large-ensemble rainfall dataset with a high resolution of 5 km d4PDF as input data to the RRI model. The RRI model performance was examined by a comparison between the simulated and observed river discharges for four large-scale flood events in the IRB. In addition, although the flood inundation depth data were not available in the CRB, an attempt was made to reproduce the greatest historical flood event of August 1981 in this basin. In addition to appropriate river discharge simulation results, the simulated inundation results also indicated that the flood inundation area was similar to the observed inundation area for this flood event. In addition, the simulated results were close to the calculated flood inundation areas estimated by Oki et al. (2017) [50]. Therefore, this study indicated that the RRI model could provide reasonable results not only for river discharge along the IRB but also flood inundation along the CRB.

The prediction results of future river floods from this study clearly indicated that the river discharge in the IRB was significantly affected by climate change impacts. The extreme river flood is expected to increase significantly by 21–24% along the IRB. This result is in close agreement with several studies conducted on the IRB. In a study by Nguyen et al. [23], the river peak discharge was evaluated in the IRB, as well as in its main sub-basins using the IFAS model. They indicated that the river peak discharge is expected to increase by 33.3% in the entire IRB. Duan et al. (2017) [57] evaluated the impacts of climate change on the hydroclimatology of the upper IRB using the downscaled large-scale Hadley Centre Climate Model 3 Global Circulation Model A2 and B2 scenarios data. They found out that the annual mean river discharge is predicted to increase for all three periods (2020–2039), (2050–2069), and (2080–2099), except the 2090s under the A2a scenario, and the largest increase is approximately 7.56%. However, a study by Sato et al. [29] evaluated the impacts of climate change on river discharge in several major river basins in Japan using the super high-resolution atmospheric general circulation model (AGCM) with a resolution of 20 km for the present climate (1980–1999) and future climate (2080–2099). The results showed that although the high temperature is projected to increase in the IRB, the monthly river discharge will not change significantly.

The CRB is the main tributary located in the downstream low-lying area of Ishikari River. Hence, this river basin is considered to be a flood-prone area of the IRB. In the past, the CRB has experienced severe flood inundation events due to short-term heavy rainfall and caused serious damage. Therefore, assessment of flood inundation risk in this basin is a crucial target for flood damage mitigation. To the best of our knowledge, there have been no studies evaluating the impacts of climate changes on flood inundation hazards in this basin area. As a result, this study can be considered a pioneering work to comprehensively assess the impacts of climate changes on flood inundation risk in the CRB. According to the simulation results, the flood inundation magnitude would be more severe in the CRB in the future compared to the historical climate simulations. The inundation magnitude is expected to increase significantly by 24.5, 46.5, and 13.8% for the inundation area, inundation volume, and peak inundation depth, respectively. In addition to the investigation changes in the inundation area, inundation volume, and peak inundation depth in this area, the inundation frequency, and inundation duration were assessed for the historical and future simulations. The results indicate that the area with the high inundation frequency corresponds to a long inundation duration and vice versa in the target basin. The downstream area of the CRB with low-lying characteristic is the most flood-prone area and should be specially considered to reduce flood damages in the future.

In addition, compared to the results of Nguyen et al. [23], this study is conducted to improve the predicted time difference between the time of peak discharge at the reference station in the CRB and that of peak water level at the confluence point intersecting the Ishikari River using the RRI model. The RRI model uses a 1D diffusive wave for the river channel and a 2D diffusive wave for the slope, which provides more reliable time difference results in the CRB. The result indicated that the short time difference (0–10 h) is predicted to increase significantly in the CRB. Hence, the flood inundation risk due to the backwater phenomenon is likely to be more severe in the future.

5. Conclusions

This paper presents the impacts of climate change on severe river discharges along the IRB and in particular the flood inundation in the CRB, a tributary of the Ishikari River, considering the extreme rainfall impacts and topographic vulnerability. We used a large-ensemble rainfall dataset with a high resolution of 5 km (d4PDF) as the input rainfall data for the RRI model. The following results were found:

-

Extreme river flooding in the IRB is significantly affected by climate changes. River discharge is expected to increase by 21–24% in the IRB. The greatest increase was observed in the middle of the Ishikari River.

-

The flood inundation is expected to be severe and higher in the CRB, with increments of 24.5, 46.5, and 13.8% for the inundation area, inundation volume, and peak inundation depth, respectively. The downstream area of the CRB with low-lying characteristics is the most flood-prone area and should be specially considered for minimizing flood damages in the future. This area is also likely to experience flood inundation with a frequency of 90–100% in the future. The inundation duration is expected to increase by approximately 5–10 h. In addition, the results indicate that the area with high inundation frequency corresponds to long inundation durations and vice versa in the CRB.

-

The predicted time difference is improved in the CRB using the RRI model compared to that of the previous study. The short time difference from 0 to 10 h is predicted to increase significantly in the CRB, which indicates that the flood inundation risk due to the backwater phenomenon is likely to be more severe in the future in this basin.

Overall, these results are expected to provide useful information for river basin management, particularly for climate change adaptation and flood damage mitigation in floodplain areas.

However, there still remains considerable uncertainty in this study. We focused on assessing the changes in natural flood inundation risk in the target basin considering extreme rainfall impacts and topographic vulnerability. Therefore, existing flood control facilities such as dam operations, pump stations, and retarding basins were not implemented. Furthermore, land-use activities as well as the socio-economic rapid development were not considered. The impacts of these factors should be considered along the target basin in future studies as the rapid development would cause uncertainty in the flood inundation prediction. In addition, the effectiveness of flood control facilities should also be considered for flood inundation hazards in the IRB as well as in the CRB to obtain more reliable projections. This study was conducted to assess changes in the future flood inundation hazards using the top 20 past and top 36 future rainfall events (corresponding to return periods equal to or larger than 150 years) [23] as input rainfall data for the RRI model. Therefore, assessing the future flood inundation risk based on yearly rainfall data for 3000 historical and 5400 future climate change scenarios extracted from the d4PDF dataset would provide more reliable future projections in the target basin.