Sediment transport due to flooding from levee breach: a numerical application of the 2004 Niigata flood

Abstract

A fatal flood disaster due to levee breach of Kariyata River occurred in the Chuetsu area of Niigata Prefecture on 13 July 2004. Many houses were swept away and substantial areas of Sanjo City, Mitsuke City and Nakanoshima Town were flooded by muddy water. Removal of mud and other sediment deposits in residential areas in the aftermath of the flood was difficult. To better understand and prevent flood disasters, a two-dimensional numerical model for unsteady flow and non-uniform sediment transport was developed using the finite difference method on a staggered grid. The model consists of two main components: (1) a simplified hydraulic scheme that predicts floodwater flow depth and velocities; and (2) a sediment transport scheme that simulates multi-class sediment accumulation. The sediment model addressed the dominant process of sediment transport, settling, deposition and re-suspension, and bed changes in the flooded regions around the break point on the Kariyata River. The model was calibrated using a combination of field data and calibration parameters derived from the existing literature. The simulated results show good agreement with field data.

Keywords:

• Multi-class sediment model
• overland flow
• heavy rainfall
• sediment transport
• numerical study

1 Introduction

Changing meteorological conditions have made it difficult to predict the occurrence time and scale of floods. The construction of river improvements cannot safeguard against larger-than-expected rainfall. River conditions are continuously changing due to sediment deposition and the growth of plants and trees. Occasionally, floods can reach a higher level than expected for a particular amount of rainfall. Additionally, adverse weather conditions caused by global warming may result in unusual rainfall patterns in certain regions. Inundation due to dam and levee breaches has caused loss of life and property (Ying et al. 2003). Consequently, the continuous study of flooding events is essential to prevent the occurrence of fatal disasters. Numerical simulation is a common approach for studying flood disasters. The present study proposes a numerical simulation method that considers sediment transport during floods.

The present paper is arranged into the following sections. Section 2 describes the 2004 flood that occurred in central Niigata Prefecture, Japan. Section 3 discusses previous studies of the 2004 Niigata flood and provides a literature review on sediment transport. Section 4 presents the mathematical foundation for the numerical model, and the basic assumptions and parameters used in the model development process. Section 5 evaluates the simulated results and presents a comparison of these results with actual on-site measurements.

2 2004 Niigata flood disaster

Severe flooding occurred in the Ikarashi and Kariyata River basins near Niigata, Japan on 13 July 2004 (Okuma 2004). Both rivers are tributaries of the downstream segment of the Shinano River, which flow across the Chuetsu region of Niigata Prefecture (Figure 1).

Figure 1 Heavy rain disaster in Niigata Prefecture (X: Break points of the Ikarashi and Kariyata Rivers)

The flooding was caused by a heavy flow discharge (1700 m³/s) in Kariyata River, which exceeded the design discharge (1550 m³/s) and broke the embankments in downstream urban areas, significantly damaging Nakanoshima Town. Large-scale flooding devastated Sanjo City, Mitsuke City and Nakanoshima Town where 15 people died, 2500 ha of land were inundated, 30 homes completely destroyed, 129 homes half destroyed, and 95 homes partially destroyed by the flow forces and the current of the flood. In addition, although the building structures remained intact, the basements of 6332 housing units were severely damaged by floodwater and mud. An additional 7260 buildings were similarly damaged, although to a lesser degree (Sato et al. 2006). In most housing units, mud was deposited as high as 150 cm. Figure 2 delineates the inundated region in Nakanoshima Town. The Kariyata River is shown in Figure 2, flowing from the bottom to top with a meandering channel along the lower half of the figure and a straight channel in the upper half of the figure.

Figure 2 Break point and inundated region at Nakanoshima Town

Even when the flooded water receded, the removal of deposited mud inside buildings was a major challenge. The massive deposition of mud also caused the floodwater to be trapped in some areas because manholes and water outlets were blocked. Removal of mud and other sediment deposits in residential areas was common and a high priority problem for the cleanup work by volunteers and civil defence forces (Figure 3).

Figure 3 Mud accumulation in Nakanoshima Town in the aftermath of the 2004 flood

3 Previous studies

In the aftermath of the Niigata flood disaster, a number of studies were completed. Hosoyamada (2005) performed a numerical analysis of this flood disaster. Fukudome et al. (2006) performed a field survey of evacuation activities in response to the flood. Hiramitsu et al. (2004) of Niigata University studied this flood event in their annual report. The previous studies had not considered sedimentation factors. By evaluating sedimentation behaviour, the present study is able to more accurately analyse the effects of the flood.

The present model includes the effects of sediment transport with re-suspension from the ground and deposition on to the ground surface in the flooded area. The accurate behaviour of flood flow is studied interactively with changing ground features and elevation. Our model identifies potential areas of massive sediment deposition, assisting with disaster prevention and post-disaster management.

Sediment transport modelling has been developed and widely applied to river management since the 1970s. In recent years, many two-dimensional (2D) and three-dimensional (3D) sediment transport models (Shimizu et al. 1990, Wu et al. 2000a) have been created to simulate complex channel evolution in more detail. Liu et al. (2002), Lumborg (2005), Lumborg and Windelin (2003), and Lumborg and Pejrup (2005) studied numerous cohesive sediment transport models. Rasmeemasmuang and Sasaki (2006a,b) performed accumulation of cohesive and non-cohesive sediment transport in Tokyo Bay. Wu (2004a) calculated the depth-average 2D flow and sediment transport model in Lower Yellow River. Mehta (1986) studied cohesive sediment transport process in estuaries. One and two dimensional non-uniform sediment transport models in open channels and channel network were studied by Wu et al. (2004) and Wu (2004b).

In the study of Kawaike et al. (2002), they developed an inundation flow model in hillside cities to estimate sediment yield using unstructured grid. The major strength of using unstructured grid is that each cell can closely represent the ground feature of the study area. However, unstructured grid approach introduces complexity in using second-order derivative term (diffusion term). In this study, we used rectangular and structured grid and incorporation of second-order derivative term is taken into consideration. The incorporation of second-order derivative term reduces numerical error and makes the model more stable. In order to closely represent the ground feature of the study area, we use 1000 × 1100 grids making every cell represent 1 m². By using detail rectangular grid, highly accurate results can be obtained.

Most sediment models have very similar structures, differentiated by different components. One of the major components evaluates bed shear stress. In most previous studies, only the contribution of water motion due to current was included in the computation of bed shear stress and several studies in coastal engineering considered the effects of current and wave motions. In the present study, current-induced bed shear stress is used because sediment is mainly transported by flood current.

The approach of the models to sediment settling velocity has another distinguishing concept. Many studies adopt constant settling velocity calculations based on Stokes' Law, while other models apply a concentration-dependent formula in which settling velocity is a function of sediment concentration. The current study calculates settling velocity of cohesive sediments using the sediment concentration.

Models are commonly used to simulate the long-term alteration of bays, estuaries, rivers and channel networks. However, these models generally assume that sediment particles are uniform with one representative size. In contrast, our numerical simulation of floodwater and sediment transport calculates impacts to residential areas due to a levee breach on the Kariyata River during the 2004 flood event. The present model establishes multi-class sediment transport models with short time-scales whereas former models used longer time-scales.

In order to effectively classify the sediment, field data from the River Management Division, Niigata Prefecture and data from the Nagaoka University of Technology ‘Niigata Heavy Rain Disaster Survey Report’ (Niigata Disaster Research Group 2005) were used in the current study.

4 Numerical method

In the present paper, hydrostatic approximation is described using three major equations, including (1) 2D non-linear shallow water equations for flood flow propagation; (2) sediment transport equations with the inclusion of dominant processes such as sediment deposition and re-suspension; and (3) bathymetry change equations for bed changes in flooded regions. The tendency for levee breaks is not evaluated in this study. This numerical simulation starts at the time of levee breach and its inflow boundary condition is based on the study conducted by Hosoyamada (2005) and Nozaki (2009). Inflow discharge hydrograph and sediment concentration hydrograph are shown in Figures 4 and 5, respectively. Inflow discharge sharply increased at the beginning because levee breach occurred in short period of time.

Figure 4 Inflow discharge hydrograph

Figure 5 Sediment concentration hydrograph

4.1 Governing equations for non-linear shallow water waves

For the development of the mathematical model, some simplifications were assumed: the flow was considered to be 2D, incompressible and hydrostatic. The density was considered to be temporally and spatially constant. The governing equations can be defined by Eqs 1–3.

where t is time; x and y are the horizontal Cartesian coordinates; h is the flow depth (m); M and N are the vertically integrated velocity components in the x and y directions (m² s⁻¹); z* is the altitude of the ground (m); A _h is the horizontal viscous coefficient (m² s⁻¹); g is gravitational acceleration (m s⁻²); and n is the Manning roughness coefficient. The horizontal viscous coefficient is calculated by , where is the bed shear velocity ; and α is the empirical coefficient 0.13. The Manning roughness coefficient depends upon the type of land use classification, such as farm or road.

For cases in which water is discontinuous, as shown in Figure 6, overflow formulae (4) are used.

Figure 6 Overflow conditions

4.2 Governing equation for sediment transport

River sediment grain sizes vary based on geographical conditions. Figure 7 shows the graph of Kariyata River's grain size distribution from field measurement (Niigata Disaster Research Group 2005). Based on this graph, we selected four sediment grain sizes that cover the range of all possible grain sizes found in Kariyata River embankment as shown in Table 1. Sediment flowing with flood discharge was classified into four categories: fine sand, very fine sand, median silt and clay particles (Yang 2003). Fine sand and very fine sand are assumed to be dominant with non-cohesive characteristics, whereas median silt and clay particles exhibit cohesive characteristics.

Figure 7 Grain size distribution of Kariyata River

Table 1 Sediment classification in the present model.

Parameter	Fine sand	Very fine sand	Median silt	Clay
Diameter (mm)	0.25	0.15	0.313	0.0039
% of sediment in floodwater	10	15	32	43

Transport of the sediment materials is governed by the advection–diffusion equation, applied with a settling term. The transport equation can be written as

where C is the sediment volume concentration of each particle (m³ m⁻³); h is water depth (m); u and v are the velocities in the x and y directions (m s⁻¹); D _h is the diffusivity coefficient of sediment (m² s⁻¹); E is the erosion rate (kg m⁻² s⁻¹); D is the deposition rate (kg m⁻² s⁻¹).

4.3 Settling velocity

The settling velocity of fine-grained cohesive particles is significantly different from that of coarse-grained non-cohesive particles. For cohesive sediments, we adopted the standard expression of concentration-dependent settling velocity shown in Eq. 6. For non-cohesive particles, individual particle settling velocities, which are determined from the balance of gravitational and drag forces (Stokes' Law) as shown in Eq. 7, can be adopted.

For cohesive particles,

For non-cohesive particles,

where d is the particle diameter (m); υ is kinematics viscosity (m² s⁻¹); c is sediment concentration (kg m⁻³); c _p is the standard sediment concentration (3.5 kg m⁻³); and are coefficients: m ranges from 1 to 2; ranges between 3 and 5; and and (Wu 2007).

4.4 Deposition and re-suspension

Cohesive sediments are deposited in the form of flocs, if the flocs are strong enough to settle through the bottom region of high shear (Metha and Partheniades 1975). Where the shear is stronger, some flocs may break up as they approach the bed and return to the flow. Krone (1962) and Metha and Partheniades (1975) investigated the deposition process of fine sediments and proposed approaches to determine deposition rates.

where a is the deposition probability coefficient between 0 and 1, which is related to the bed shear stress .

where is the critical bed shear stress below which all sediment particles will deposit on the bed, and is the critical bed shear stress above which all sediment particles remain in suspension, yielding a zero deposition rate. Mehta and Partheniades (1975) found that might be larger than zero and might vary from 0.18 to 1.1 N/m². This study used critical bed shear parameters of this range.

The depositional flux rate is computed using the following sediment formulation for non-cohesive sediment (Sekine 2004).

Erosion of cohesive sediments is affected by flow conditions, sediment properties and bed configurations. According to Partheniades (1965), the surface erosion rate is a linear function of the dimensionless excess shear stress:

where E ₀ is the empirical erosion rate constant, and τ_e is critical shear stress for erosion.

The erosion flux rate is computed using the following sediment formulation for non-cohesive sediment (Sekine 2004):

where E _s is rolling up speed calculated by with and in which rep is calculated by .

The critical shear stress values for deposition and for re-suspension are dependent on the sediment classifications. Coarser-grained sediments are more difficult to re-suspend and deposit much easier than the finer-grain sediments. The parameters were selected carefully according to suggestions of Li et al. (1994). For sand bottom, critical shear stress for erosion is set to be 2.0 N/m². For the reach covered by cohesive sediments, the values were given ranging from 1.3 to 1.5 N/m². Critical bed shear stress for re-suspension is shown in Table 2.

Table 2 Critical shear stress parameter used in the present study.

Parameter (critical bed shear stress)	Fine sand	Very fine sand	Median silt	Clay
Re-suspension (N/m^2)	2.15	2.1	1.5	1.3

Bed load transport is neglected because suspended load is dominant in flood waves in residential areas. Although the process of sediment re-suspension closely relates to the local bed material composition, the adjustment process of bed materials are not considered in this study because of the following reasons: (1) a large portion of the study region are covered by roads and buildings, (2) the whole simulation is done within a short period (3 h), (3) deposition is more dominant than re-suspension in this study. The temporal change in land height is calculated in Eq. 13:

where λ is porosity of bed soil, and is height of bed.

4.5 Calculation conditions

In the Kariyata River calculation area, there are 1000 rectangular grids in the x-direction and 1100 rectangular grids in the y-direction. Each grid is 1 m × 1 m. The north–south length is 1100 m and the east–west length is 1000 m. Figure 8 shows the altitudes of all numerical grids covering the entire calculation domain for Nakanoshima Town. The fine land height data were obtained from an airborne survey using a laser profiler. Small rectangles in the figure show the roofs of houses and buildings.

Figure 8 Altitude in Nakanoshima Town (X-direction: 1000 m, Y-direction: 1100 m)

To determine the floodwater and sediment deposition depth, this simulation uses the inflow discharge hydrograph based on the maximum flood discharge – 400 (m³ s⁻¹) – shown in Figure 4 and inflow boundary conditions for the sediment concentration hydrograph, based on an 8% maximum value shown in Figure 5, are given at the levee break point. The simulation duration is 3 hours. Each time step is 0.01 s. The time step is adjusted based on the Courant number. A mixture of all four sediment types is used in the ratio shown in Table 1.

5 Simulation results

The present model is capable of simulating the time series for mud accumulation, sediment concentration, water depth, and floodwater flow velocities in Nakanoshima Town. The simulation results for ‘mud accumulation, sediment concentration, water depth and flow velocities’ at 1 h intervals are shown in Figure 9. Figure 9(a) shows the evolution of mud accumulation during the computational period. The sediment accumulation depth is high in the middle part of the study region, particularly near the break point. The total concentration is also high near the break point, as shown in Figure 9(b). After a 30 min simulation time, floodwater spread throughout the entire study region, as is seen in Figure 9(c). In Figure 9(d), the arrowhead shows the flow direction and concentrated vector shows the high flow velocity. According to the vector map, high flow velocity can be observed in the region near the break point and along the roads.

Figure 9 Time series of accumulation results

Computational results for each sediment concentration during flooding are shown in Figure 10. The concentration of coarse-grained materials is highest near the break point, as shown in Figure 10(a). Coarse-grained particles mostly settle near the break point because of their higher settling velocity, while fine-grained materials with lower settling velocities are transported further by currents and settle farther from the break point, as shown in Figure 10(b)–(d).

Figure 10 Sediment concentrations for particle size classification

Mud deposits along the residential area roads and near the campus of houses. In contrast, re-suspension occurs in large open areas such as farms. This simulation does not consider the conditions of floodwaters entering buildings or mud depositing under houses because housing structures differ and details are difficult to determine. In this simulation, re-suspension is lower than deposition as a result of residential land use or roads in the study area.

At the nearest place of the breaking point, less accumulation and higher erosion occurs because of the high flow rate, as shown in Figure 9(a). Figure 11 shows that a depression region occurred in the simulation as well as under actual conditions. According to Survey Report, the depression depth is nearly 1.6 m at the flat region and 7 m at the levee (Niigata Heavy Rain Disaster Survey Committee, Niigata Prefecture, River Management Division 2004). The simulation result and field measurement are shown in Figure 12.

Figure 11 Depression region along the Kariyata River

Figure 12 Depression depth

To analyse the effects of meteorological and hydrodynamic factors on sediment accumulation, the computational results at representative study points were collected at 3.75 min intervals and compared with the field data collected by Niigata Prefecture River Management Division.

The location of study points and actual field data are shown in Figure 13. The comparison of accumulated mud depths (from simulation) and field measurement is shown in Figure 14. The lines show increasing simulated mud depths during 3 h. The symbols show the depths of accumulated mud at respective points after flood disaster.

Figure 13 Actual sediment accumulation data (unit in cm) (Niigata Prefecture River Management Division)

Figure 14 Accumulated sediment depth at seven points

The simulation reveals that sediments discharging from the break points gradually deposit within the flooded urban area. The rate of mud accumulation varied with time and distance. Mud accumulation depth sharply increased during the initial simulation time. However, after exceeding a simulation time of 30 min, the mud accumulation rate became nearly constant.

Points A, B, C and D are located at building campuses. Points E, F and G are located along a road. Among them, Point A is nearest to the break point and lies at the flow direction as shown in Figure 13. Coarse-grained sediment could be transported to Point A, at which coarse-grained sediment was in greater abundance. Point G also lies within the flow direction and the current transferred coarse-grained sediment to these points, as shown in Figure 15. Among the points along the roads, Point E shows the highest deposition because floodwater reached this point faster than the other points that lie along the roads. At all points on the building campuses, fine-grained deposits were more abundant than coarse-grained sediments due to the distance of the flow path from the break point, from which currents could not transport much coarse-grained.

Figure 15 Sediment height at each point

The pattern of sediment accumulation is explained by the fact that a higher rate of accumulation occurs near the break point and the rate decreases with increased distance from the break point. This description is also consistent with the field investigations.

The positions of seven points and actual maximum water level during flooding are shown in Figure 16. Figure 17 shows simulated water depth of seven points around the flooded region. It can be clearly seen that the current simulation results for maximum water depth at the collection points are in good agreement with the interpolation results of field traces data of Niigata Disaster Research Group (Niigata Disaster Research Group 2005).

Figure 16 Maximum water level estimated using actual water level field traces during the 2004 Niigata flood (unit in metre) (X-direction: 515 m; Y-direction: 594 m)

Figure 17 Water depth at seven points

6 Conclusions

This case study demonstrates the applicability of combined modelling tools to the characterization of spatial and temporal variation in sediment transport in an urban-flooded area. The present paper describes the sediment accumulation, water depth, flow velocities and sediment concentrations in the urban area of Nakanoshima Town in Niigata Prefecture, simulated using a multi-class sediment model and integrated with a hydrodynamics model. The proposed model adopts the finite difference method to solve 2D shallow water equations for unsteady flow. The sediment transport equation is also discretized using the finite difference method for a staggered grid. The model can account for the sediment transportation, settling, deposition and re-suspension processes. Numerical results show that the temporal and spatial characteristics of sediment accumulation are related to metrological and hydrodynamic features, particularly the local bed shear stress. Computed sediment concentration (8%) used in this simulation produces results that are in good agreement with the field data. The simulated water elevation and mud accumulation depths are consistent with those of the field investigation.

This model can provide consistent information regarding sediment transport and floodwater behaviour, providing the complete picture of the flood needed to optimize flood management. Moreover, unlike other numerical models, the behaviour of floodwater is simulated by incorporating dynamic erosion and deposition of multi-class sediments, which alter ground features simultaneously. These dynamic parameters provide the most realistic view of flood behaviour. Consequently, the proposed model can be successfully applied to the prediction of flood flows resulting from the overflowing of riverbanks, dam breaks or breached levees caused by heavy rainstorms. The results could also be used to support the development of flood forecasting and warning systems.

Notation

The following symbols are used in this paper.

h	=	flow depth
M, N	=	vertically integrated velocity components in the x and y directions
A h	=	horizontally viscous coefficient
z*	=	bed elevation
u *	=	bed shear velocity
w s	=	settling velocity of sediment
u, v	=	components of flow velocity in x and y directions
D h	=	diffusivity coefficient of sediment
c, c p	=	sediment concentration
τb	=	bed shear stress
υ	=	kinematics viscosity
λ	=	bed material porosity
Δz b	=	total bed change at time step Δt
t	=	time
x,y	=	spatial coordinates
g	=	acceleration due to gravity
n	=	Manning roughness coefficient