High-resolution distributed model to simulate erosion and sedimentation in a steep basin: a case study of the Akatani River Basin, Kyushu, Japan

: In this study, we developed a distributed rainfall-runoff and sedimentation model based on one-dimensional kine‐ matic wave equations. Physically-based rainfall-runoff and erosion-sediment processes were coupled and solved for each spatial grid, whilst the spatially distributed grids were connected to each other to allow for space-and-time move‐ ments of water and sediment. The model was applied to the Akatani River basin of the Chikugo River in Kyushu, Japan using a 10 m high-resolution digital elevation model and eXtended RAdar Information Network (XRAIN) data as a time-and-space distributed rainfall input of the northern Kyushu heavy rainfall event in July 2017. Our results indi‐ cate that the rainfall-runoff hydrograph and sediment flow results are in agreement with the collected field data, and elevation of the river bed after the disaster was successfully reproduced by applying a sediment theory to estimate river bed variation. In addition, we found that sediment transport results are sensitive to model spatial resolution. Our simu‐ lation model is intended for use with basins that feature steep slopes and are prone to erosion and shear strength reduction after heavy rainfall events. Hence, this model can be applied to give early warnings by identifying critical erosional areas during forecasted heavy rainfall events.


INTRODUCTION
Recent extreme rainfall events during the rainy season in Japan have resulted in flash floods and sediment-related disasters. The extreme rainfall event in Kyushu in July 2017 has caught the attention of the local authorities and scientific community after the floods and sediment-related effects were reported. After this extreme event, a sedimentation analysis based on laser profiler and aerial photographs before and after the event was conducted by the River and Sabo Restoration Technology Review Committee on the right bank of the Chikugo River (RSRTC, 2018). The report by the committee includes a thorough geophysical assessment of the considerable volume of sand which caused an increase in the elevation of the river bed.
Numerical simulations in most hydrological models rely on accurate rainfall data, topography, and soil mapping information. Regarding accurate flood simulations, the volume of transported and deposited sand in a river bed must also be considered as a critical factor (Kawaike et al., 2018). Although currently available hydrological models, such as the Soil & Water Assessment Tool (SWAT) (Arnold et al., 2012), the TOPography-based hydrological model (TOPMODEL) (Beven and Freer, 2001), and the Water Erosion Prediction Project (WEPP) model (Laflen et al., 1991) that were developed in the last 10 years have brought satisfactory results, there have been discussions regarding the physical assumptions and parameter calibration processes (Hajigholizadeh et al., 2018). The conceptual models listed above simulate sediment routing applying physical process, nevertheless the transport mechanism for high erosive and steep rivers requires a specific model approach. Egashira and Matsuki (2000) proposed a method to resolve sediment transport for a complex channel network into a single channel unit with two inflow points and one outflow point, which modeling is considered to be relevant given satisfactory results in hydrological and sediment predictions. Thus, we designed a fully-distributed model composed of water flow and sediment transport with a very high-resolution digital elevation model. Furthermore, we applied an extended algorithm to allow the connection of every differential element according to the digital drainage model (DIR) as it will be discussed in the next section. After the described background research, our study explores the applicability of the aforementioned distributed hydrological model to study the 2017 Northern Kyushu extreme rainfall event at the Akatani River basin (21 km 2 ) with reference to the RSRTC assessment results (RSRTC, 2018). This paper describes a new approach in the development of a spatial distributed model that considers hillslope rainfall-runoff erosional processes as well as routing and transport processes for a river channel system. This approach takes considerations from the works by Ferro and Porto (2000), who emphasize the relevance of modeling hillslope and channel erosion as two processes with different governing equations. For this purpose, the estimation of soil erosion rate, sediment transport, and deposition at allocated critical spots is developed.
Firstly, we described the theory, methodology, and equations for the distributed rainfall and sediment runoff model. The interaction of the erosion process with the rainfallrunoff model is a key component of the modeling, and is numerically solved by using a finite element method (FEM). Secondly, we calibrated model parameters for rainfall-runoff and sediment production-transport simulations. Finally, we evaluate the sediment yield simulations by the change in elevation of the river bed before and after the 2017 Northern Kyushu extreme rainfall event. Furthermore, we analyzed the sediment simulation results using 10 m and 100 m resolution DEM (digital elevation model). Based on the simulation results, we discuss the sensitivity of model resolutions on hydrological and sediment predictions.

Hydrological modeling
A distributed hydrological model consists of rectangular slopes connected to each other as shown in Figure 1b. The equations governing flow for each slope are based on the kinematic wave equations as follows (Tanaka and Tachikawa, 2015): where h is water depth (m); q is discharge per unit slope width (m s -1 ); r is rainfall intensity (mm h -1 ); α = √i/n c B -3/2 ; and m = 3/5. The parameters i, n c and B represent the slope gradient (m m -1 ), Manning's roughness coefficient, and the width (m) of river channel, respectively. This dischargedepth relationship for hillslope flow was extended to integrate overland flow with subsurface flow simulation according to saturation level and layer thickness. The mathematical relations for subsurface flow are as follows: where h is the water depth (m) equivalent to the water content in the effective porosity, d a is water depth (m) equivalent to the maximum water content in the effective porosity; d c is the water depth (m) equivalent to the maximum water content in capillary pores; k c relates to the hydraulic conductivity (m s -1 ) when the capillary pore is saturated; β is an exponent parameter that describes the relationship between hydraulic conductivity and saturation; k a is saturated hydraulic conductivity (m s -1 ); and k a = βk c according to the continuity of the relationship shown in Equation (3). Furthermore, n s is the Manning's roughness coefficient for overland flow in the slope runoff cells; and i represents slope gradient (m m -1 ). As the subsurface soil layer of hillslopes conductivity is considerable, we assumed that the surface soil layer is unsaturated before a rainfall event.
Overland flow is computed after the surface soil layer is saturated. The above relationships were discretized and applied as differential FEM equations to be computed under a four-point implicit scheme. For the distributed model, a digital elevation model (DEM) takes a relevant role for topographical representation of the study basin. The DEM used as the base topography before simulations is referred to throughout this paper as "Before (MLIT)". After the extreme event in July 2017 in the northern Kyushu region, the RSRTC carried an aerial survey to collect raw data by a Laser Profiler (LP) equipment at Akatani River basin. After this survey a 10 m spatial resolution DEM was produced, this was used as the topographical reference to analyze our results; it is referred to as "After (MLIT)" for Figures 5 and 6. The DEM and DIR are shown in Figure 1a and Figure S1 at 10 m resolution respectively. These two digital models generate flow accumulation information according to the steepest downstream direction. The digital accumulation grid "flowacc" counts the number of connected upstream elements. The element with the highest "flowacc" number was assigned as the basin outlet (P1 in Figure 1(a)). The flow accumulation map for this simulation is presented in Figure S2. Rainfall data used for the simulation was obtained from XRAIN with a 10 m spatial resolution, as provided by MLIT (Ministry of Land, Infrastructure, Transport and Tourism) and Transport and Tourism Kyushu Regional Development Radar rainfall data (XRAIN), as collected by the MLIT was used for the rainfall-runoff simulations of the model. The rainfall data used ranges from July 5 th at 7:00 to July 6 th at 23:00 for a total interval of 40 hours. The maximum rainfall recorded in the catchment was 107 mm h -1 on July 5 th , which is considered to fall within the extreme rainfall category.

Calibration of the rainfall-runoff model
Although the flow at the outlet of the Akatani River basin is not directly gauged, it was estimated according to the rational method discharge estimations published by Chikugo River Committee reports (RSRTC, 2018). Thus, at first, the model parameters for the 100 m rainfall-runoff model were calibrated with the observed discharge from the Kagetsu River basin (Kagetsu parameters, Table SI), located adjacent to the Akatani River basin. The initial parameter values were set based on prior knowledge of the model applications. Then, the adjusted model parameters of 100 m resolution (identified at the Kagetsu basin) were calibrated by trial and error for the 10 m and 100 m distributed models of Akatani River basin to fit the estimated hydrograph by RSRTC in terms of the maximum flow discharge and time to peak. For these calibrations, the roughness coefficients (slope and channel) were changed to adjust the time to peak, then the depth soil of capillary and finally hydraulic conductivity were tuned in the described order to match the peak discharge. Parameter values related to rainfall-runoff processes such as hydraulic conductivity, soil depth, and roughness coefficient parameters are summarized in Table SI.

Sediment erosion-deposition modeling
The hillslope erosion process was modeled by considering raindrop impact and surface flow as the driving mechanisms for sediment production. Due to the surface flow, the concept of stream power (Govers, 1990) was adopted for the erosion process, in which soil detachment occurs when the stream power exceeds the critical stream velocity.
For river routing processes, the generated sediment from upstream cells was incorporated to the river system. The distributed model has grid cells in which either hillslope erosion − or river processes are dominant. The grid cells are classified to the hillslope-related or river-related grid cells according to the flow accumulation number with a threshold value. Both flow direction and accumulation information derived from the DEM were used as the governing parameters that control the modeling of sediment generation, transport, and deposition (refer to Figure 2 for a graphical description of these two processes).

Governing differential equations for erosion and deposition
The continuity of sediment mass is the basic governing equation for modeling soil dynamics (Morgan et al., 1998). This equation is applied to every grid cell in the developed model as follows: Here, C S is the sediment volumetric concentration (m 3 m -3 ); A is the cross-sectional channel area (m 2 ); and Q is the flow discharge (m 3 s -1 ). The expression q S (x, t) represents the rate of external input or extraction at each grid cell explained in the next sub-section. The net erosion rate e(x,t) in the left-hand side of the equation represents the sum of detachment by surface flow DF and rain drop impact DR according to the following equation:

Detachment by flow
Soil detachment by flow DF is calculated by using the sediment erosion-deposition theory. This takes into account that DF depends on the difference between the current sediment concentration C S (m 3 m -3 ) and the transport capacity T C (m 3 m -3 ), which is the maximum concentration that the surface flow can hold. According to Govers (1990), the expression for DF is presented as follows: where β S is the flow detachment efficiency coefficient with a normal range from 0.4 to 1.0 according to soil cohesion; w is the surface flow velocity (m s -1 ); and v S is particle settling velocity (m s -1 ). According to Equation (6), detachment by surface flow will decrease in proportion to an increase in sediment concentration. When the concentration reaches transport capacity, detachment by flow does not take place.

Detachment by rainfall impact
Soil detachment due to the impact of rainfall DR is a sediment process related to hillslope simulations. Most laboratory experiments give evidence as to how the momentum of a raindrop falling on the hillslope surface generates erosion. One of the effects of rainfall is an initial compacting of the soil surface. Another effect, and the main source of the sediment production, is related to raindrops carrying the kinetic energy necessary for soil detachment (Poesen and Govers, 1985). Figure 2. Schematic representation of the physical sediment transport process. The hillslope process is dominant at upstream cells (green shade zones), before reaching accumulation threshold. Afterwards, the conditional statement assigns the following elements to river process where water flow also produces erosion (blue shaded zone) The following relationship was proposed by Torri et al. (1987) based on previous research: Here, k sd is the soil detachment index ranging from 0.01 to 10 g J -1 according to soil texture and shear strength; ρ s is the soil density (kg m -3 ); KE = 8.95 + log(r); e -Zh is a correction factor that reduces the capacity of raindrop impact for sediment production, r is the rainfall intensity (mm h -1 ), h is the surface runoff depth (m), and Z is the experimental exponent ranging from 0.9 and 3.1. Raindrop detachment is reduced when surface flow is present, and a significant surface flow depth prevents soil detachment due to a reduction in raindrop impact. River bed elevation change was calculated for every cell identified as river element applying the Exner equation, whilst the differential grid element elevation was calculated using the sediment flux (Julien, 2018) as follows: Here, z b is the bed elevation; w is the unit width (m) of the cell element; δx and δt are the differential time and space elements; Q S is the volumetric sediment discharge (m 3 s -1 ); and p o is the soil porosity set at an experimental value of 0.43.

NUMERICAL SOLUTION
Firstly, we carried out the rainfall-runoff simulations using Equations (1), (2) and (3). The hydrological model was run using the parameters presented in Table SI. Equations (5), (6) and (7) were used for the numerical solution regarding the sediment process, and the initial sed-iment concentration was set at a very small value near zero. The soil parameters and equation factors are presented in Table SII. For hillslope processes, DR and DF were computed for every grid cell, whereas for river processes, only DF was computed as raindrops do not make physical contact with the riverbed (Figure 2 and Equation 7).
The flowchart in Figure 3 schematically displays the algorithm for parallel processing of both the hydrological and the sediment model. For the entire simulation processes, the flow rate of the next time step for each cell Q n [i] was calculated using the 4-point box scheme and stored in a temporal memory, after which the memory was accessed for the parallel sediment concentration calculations C S[i] , where i represents a spatial discretization index. These routine calculations were performed for a set time step at the given [i] spatial divisions.
The developed model is based on one-dimensional cell distribution, hence the direction for flow discharge Q n and sediment discharge (Q S = Q · C S ) is controlled by the flow direction DIR under the eight-direction algorithm (Lehner et al., 2006). The mathematical relations and the program code have been modeled to reproduce the sediment transport processes to the steepest down-slope neighbor grid cell. A graphical representation of this drainage process can be found in Figure 1b. The model was coded in C++ language, and compiled in Linux environment. Figure 4 shows the predicted discharge hydrographs for the extreme rainfall event using the 10 m and 100 m resolution rainfall-runoff models, and the hydrograph produced by the synthetic rational method by the RSRTC (2018). The maximum rainfall intensity during the flood event is more Figure 3. Model processes flowchart. The output from the kinematic wave equations being current discharge Q C[i] and next time step discharge Q n [i] are stored in a temporal memory. The subscript [i] denotes the current spatial division, whilst the initial simulation loop i = 0 considers the initial surface water depth by using initial runoff height (mm h -1 ). Alongside other inputs at the upper box, the sediment model calculates the concentration CS c and sediment flow Q C ending the first iteration than 100 mm h -1 and the catchment size is 21 km 2 . In these conditions the estimated hydrograph by the RSRCT represented by the dotted line using the rational method is considered to be reasonable.

Runoff simulation
The simulations for 10 m and 100 m give overall similar runoff results. The peak river discharge hydrograph at the outlet (P1) of the Akatani River is 480 m 3 s -1 using the developed model, while the RSRTC discharge is 520 m 3 s -1 . Hydrologic features were similar in time to peak and peak discharge. However, this model is adding subsurface flow and the consequence is that the recession limb will be gentler. This is observed to hydrographs at both resolutions.

Sediment transport simulation
Regarding sediment transport and sedimentation results, the high-resolution DEM and DIR allowed for higher convergence in results as we found that the model is sensitive to spatial resolution. The peak sediment concentration is 0.061, 0.695 and 0.080 for control points P1, P2 and P3 respectively. It is expected that a river cell in the lower or middle basin will take a higher volumetric concentration but also the highest sediment transport rate in m 3 .
The highest intensity in mm h -1 for the distributed rainfall is located in the upper catchment; thus, C S will surge due to DR and reach its peak value at the upper or middle basin before or while the peak rainfall was recorded. This is according to a physical process where sub-basin outflows enter the river system. Indeed, point (P2) in Figure 1(a) takes the peak value in this simulation at 18:00 h on July 5 th . According to Equation 7 the sediment concentration will follow a similar pattern to the rainfall intensity if the water depth is close to zero. Afterwards, C S is exponentially reduced by the right-hand side expression e -Zh , where h is the surface water depth. Following the initial stage, the total erosion process e(x, t) from Equation 5 is controlled by detachment by flow and the sediment continuity Equation 4. Finally, after a rainfall-runoff cycle has finished and the water depth h at the slopes returns to a minimum value, a new rainfall event can trigger sediment production as shown in the final hours of Figure 5a. We verified that the DR and DF run accordingly at separate threads, afterwards we fully integrated to the complete expression of Equation 4.
The simulated main river profile was generated taking the reference from the original bed elevation "Before (MLIT)" DEM. The bed elevation change for every differential element is according to the Exner relations. For instance, Figure 5b shows the elevation relationship with time for three selected points (P1, P2, and P3) where the level change is stable after 20:00 h on the first day. To generate the profile for the full river extension, the same algorithm was used for every element to determine the final elevation profile shown in Figure 6. In this figure the simulated river elevation profile is compared to the laser profiler topography "After (MLIT)" produced by RSRTC. The bed profile at the downstream river supports the expected aggradation trend.
For a deeper understanding of the results, a comparative sediment transport balance among the simulations with the 10 m and 100 m model and RSRTC LP processed data is summarized in Figures S4 and Table SIII. In this figure we can understand that the simulations are sensitive to the DEM resolution. The 10 m results show sediment transport rates are closer to the upstream Otsushi and middle Akatani basins. The estimated sediment transport using the 100 m resolution model is lower in comparison.
Consequently, we found that the developed model is able to reproduce the sediment transport and change of topography using a spatially distributed rainfall at the upper catchment, a steep basin topography and soil parameters set to unconsolidated mountain region. The critical sediment deposit was located at the Akatani River middle basin according to Figure S4 and the profile estimation according to the results in Figure 6. The sediment transport balance for both resolutions after the event and comparison to RSRTC Report can be found on Table SIII.

CONCLUSION
A flow discharge and sediment transport parallel model is an effective approach to simulate erosion/deposition for hillslopes and rivers. It was able to reproduce according to a physical process, and field data collected before and after the event (Figure 6) supports these results. The developed model shows an improvement in sediment simulations using a 10 m resolution DEM. Compared to lower resolutions, the slope represented by the developed 10 m model allowed for higher spatial divisions and therefore, higher convergence. Furthermore, the model allows for the estima- tion of a sediment budget, which describes the volume of sand from the river bed that was eroded or deposited. Figures 5, 6, S3 and S4 provide an understanding of how the main river section downstream suffered accretion and the upstream basin the most erosion.
From Figure 5a, we can conclude that the theoretical assumption adopted for the sediment production for slopes corresponds to the real physical process, as the sediment production rate started to decline after reaching the peak discharge.
A limitation of the developed model is that it cannot account for large diameter sediment particles regarding river processes, such as boulders on the river bed. Another consideration is that the flow direction (DIR) is processed before the simulations making the drainage path constant throughout the simulation, thereby restricting the model to a one-dimensional change after sediment yield. This is the reason we recommend using this kinematic wave model for steep basins where the drainage path is constant in time.
In summary, it is necessary to collect more data before and after sediment related disasters to improve the simulation of the physical processes, especially for the location of critical sediment deposits or erosive areas.

ACKNOWLEDGMENTS
The 10 m resolution DEM and DIR was provided by the River Department, Kyushu Regional Development Bureau, Ministry of Land, Infrastructure, Transport and Tourism, Japan.  Figure S2. Flow accumulation number plot Figure S3. Akatani River sub basin list and location Figure S4. Sediment flow balance per sub basin and basin outlet Table SI. Hydrological model parameters  Table SII. Sediment model parameters. Particle size, settling velocity, and conductivity are based on average values for mountain regions in Fukuoka. Detachability and stream power are based on empirical values according to Morgan et al. (1998)