Predicting River Embankment Failure Caused by Toe Scour Considering 1D and 2D Hydraulic Models: A Case Study of Da-An River, Taiwan

: Physically based numerical models can predict scour depth at embankments located in bend reaches. However, methodologies for utilizing these numerical models to assess the risk of reinforced concrete embankment failure are rarely investigated. Therefore, a new assessment methodology is proposed to predict the riverbank failure caused by bend scour. The methodology is primarily based on a bend scour simulation model that integrates a one-dimensional (1D) hydraulic model, a two-dimensional (2D) hydrodynamic ﬁnite-volume model, and an empirical equation of bend scour prediction. The model was ﬁrst applied to the Shuiwei Embankment located in a river bend reach of Da-An River in Taiwan and veriﬁed against results from the 1D hydraulic model and ﬁeld data. The model was then used to simulate 2D ﬂow ﬁeld and the temporal evolution of bend scour depth under di ﬀ erent return period ﬂood events to examine the relationships between river discharge, water level, shear stress, and bend scour depth. The inﬂuence of shear stress on the stability of toe protections was also investigated. The ﬁeld data (from two events) and numerical solutions (four scenarios) were assessed to conceive two empirical equations for predicting shear stress and bend scour depth. A new assessment methodology was proposed using these two equations to predict the risk of river embankment failure during ﬂood periods. The proposed methodology can be easily applied in other disaster-prone regions to mitigate the e ﬀ ects of disasters caused by bend scouring.


Introduction
The principal purpose of river embankments is to reduce flood risk and protect properties. Taiwan is located in the Pacific Ocean and is struck by three to four typhoons per year on average. Typhoon-induced floods frequently cause disasters. For example, Typhoon Morakot hit Taiwan in 2009, causing severe river embankment failures and numerous deaths. Several mechanisms may lead to river embankment failure during flooding periods, including overtopping, toe scouring of the foundation, seepage, or piping and slope sliding of river embankments [1][2][3][4]. To reduce the risk of riverbank failure disasters, studies on river embankment failures that involve considering different failure mechanisms are required.
Riverbank failure is a complex spatiotemporal process involving fluvial vertical and lateral erosion, seepage erosion, and geotechnical failure. Several empirical models and process-based models have been proposed to assess river embankment failure. Langendoen and Simon [5] applied the channel evolution model termed CONCEPTS (conservational channel evolution and pollutant transport system) to simulate the processes of river embankment failure in a bend reach on Goodwin Creek. Hossain et al. [6] employed a limit equilibrium stability analysis and steady-state seepage analysis to investigate river embankment failure in the Manu River. To account for seepage, Polemio and Lollino [7] adopted a finite element method to analyze the failure process of embankment slope stability in Apulia, Italy. Midgley et al. [8] proposed a bank stability and toe erosion model (BSTEM) to predict riverbank stability in the Barren Fork Creek streambank. The BSTEM can be used to estimate the timing of riverbank collapse and riverbank retreat from noncohesive and cohesive fluvial erosion. Huang et al. [9,10] employed a limit equilibrium method to assess river embankment stability in southern Taiwan for various flood return periods. In their adopted method, three failure mechanisms can be studied: overtopping, levee foundation scouring, and slope sliding. Hossain and Haque [11] employed numerical finite element software to evaluate the influence of short-term rainfall on embankment stability in the Jamuna River.
These studies have indicated that river embankment failure models are capable of simulating hydraulic or geotechnical processes involved in embankment failure. Furthermore, combinations of physically based numerical models and riverbank failure models have been used to predict river embankment failure. For example, Faeh [12] proposed a 2D depth-averaged model to simulate the overtopping breaching of embankments in the Elbe River. Jia et al. [13] presented a three-dimensional (3D) model of riverbank erosion in the Shishou bend of the middle Yangtze River. This 3D bank failure model simulates riverbank erosion of noncohesive material. Mizutani et al. [14] proposed a numerical model to simulate river embankment failure caused by overtopping flow that accounted for infiltration effects. Their proposed numerical model comprises three modules: 2D shallow-water flow with nonequilibrium, sediment transport, seepage flow, and 2D slope stability. Lai [15] presented a 2D morphodynamic model coupled with a riverbank erosion model, to model both vertical and lateral riverbank erosion. The model is robust and predicted bank retreat distance in studies in both the Middle Rio Grande and Cho-Shui River. Rousseau et al. [16] coupled the Telemac-Mascaret model with a geotechnical module to simulate riverbank erosion in the meandering river reach of the Medway Creek, London, Ontario. Stecca et al. [17] proposed a simplified algorithm for noncohesive riverbank erosion in the Selwyn River of New Zealand to reduce the 2D morphodynamic model to a cross-sectional model.
These studies have primarily focused on assessing earthen types of riverbank failure. Few studies have assessed embankment failure in reinforced concrete. In Taiwan, the Water Resource Agency (WRA) [18] reported that embankment toe scouring is the most common failure mechanism in reinforced concrete river embankments and frequently occurs in a bend reach. When modeling the bend scouring caused by floods, a challenge is to develop a numerical model that accurately consider physical processes of scour. Therefore, Guo et al. [19] proposed a simplified simulation method for simulating flood-induced embankment toe scouring in bend reaches. The method involves the use of two models: one predicts a 2D flow field and the other estimates the temporal evolution of bend scour depth near the embankment. The agreement between the simulated and field-measured bend scour depth against time at the Shuideliaw Embankment on the Cho-Shui River is good. However, further investigation of other rivers is required to assess the model applicability. Furthermore, field scouring data have become available and could be incorporated in a new bend scour computation equation. Moreover, methods for embankment toe protection could be considered for enhancing the practical applications of scour models in the riverbank failure assessment. In the present study, the simulation method by Guo et al. [19] will be extended to link with 1D model for predicting river embankment failure caused by toe scour.
The purpose of this paper is to present a novel bend scour simulation model based on the simulation method reported by Guo et al. [19]. The present scour model integrates the 1D module of the WASH123D (Watershed Systems of 1-D Stream-River Network, 2-D Overland Regime, and 3-D Subsurface Media) model [20], hereafter referred to as the WASH1D model, the SFM2D (shallow-water finite-volume model) [19], and a predictive equation of bend scour depth [21]. Through the boundary connection of 1D-2D models, it is possible to achieve more accurate flow and scouring at embankment toe and apply the 1D model elsewhere in the river reach. In this way, the computational time can be saved because the number of cross-sections in 1D model are relatively lower than those of computational grids in 2D model. The SFM2D model linked with WASH1D applied to the river flow is thus suitable for modelling the unsteady short-term toe scouring. Two bend scouring events that occurred in 2017 and 2018 in a braided reach of Da-An River, Taiwan, were selected as case study sites for evaluating the applicability of this integrated model. To adequately represent river bed topography in the scour model, an unmanned aerial system (UAS) was created and implemented to collect high-resolution topographic data as the basis for bend scour simulations. Furthermore, the effects of shear stress and bend scour depth on the risk of riverbank failure were studied. Correspondingly, a predictive methodology that considers methods of toe protection is proposed to assess river embankment failure.
The contributions of this study include (a) the application of the bend scour simulation model to the practical river flow simulation in the natural-irregular riverbed, (b) two proposed empirical equations for practical use in predicting shear stress and bend scour depth, and (c) the new proposed assessment methodology for rapidly predicting the risk of river embankment failure prior to field site investigation.

Bend Scour Simulation Model
According to Melville and Coleman [22], bend scour is defined as the short-term general scour. The special feature of bend flow leads to the presence of secondary currents and causes increased scour near the embankment toe. Field observations at a specific location could be helpful in assessing the magnitude of bend scour. In addition, numerical analysis has become popular in investigating scouring. Therefore, it is desirable to propose a bend scour simulation model that practically predict the flow field and temporal evolution of bend scouring so that simulation time efficiency is achieved. This may be performed by following the framework reported by Melville and Coleman [22].
Scour analysis can be undertaken at appropriate scales within the scale connection. The five appropriate scales are catchment scale, stream section scale, river structure far-field scale, river structure near-field scale, and local scour scale. Different analysis tools are used for different spatial scales. The bend scour simulation model proposed in this study is based on a boundary connection of the WASH1D model, SFM2D model, and empirical-based bend scour equation.
The linking procedure is conducted in the following steps: (a) The WASH1D model for unsteady river flow is employed at the stream section scale to simulate the flow parameters (river discharge, river velocity, and water level) as input data for 2D bend scour modeling at the river structure far-field scale, (b) the SFM2D model is employed to predict the flow velocity and water depth near the river structure, and (c) the collected results values are input into the empirical equation to predict the temporal evolution of the localized scour depth. Summaries of the models and the boundary connection of WASH1D and SFM2D are presented in this section.

WASH1D Model
The WASH123D model consists of three basic modules: a 1D river routing module, a 2D overland inundation module, and a 3D subsurface computation module. The model has been successfully employed for various applications, including channel flood routing [23,24], river discharge estimations [25,26], and 2D inundation simulations [27,28]. For 1D river routing, the WASH1D model solves the following continuity and momentum equations: where t is time; x is the axis along the river direction; A is the cross-sectional area of the river; Q is the flow rate of the river; S S is the human-induced source; S R is the source caused by rainfall; S E is the sink caused by evapotranspiration; S F is the source caused by exfiltration from the subsurface media; S 1 and S 2 are the source terms contributing to overland flow; V is the river velocity; g is gravity; Z 0 is the bottom elevation; h is the water depth; ρ is the water density; ∆ρ = ρ − ρ 0 is the density deviation from the reference density (ρ 0 ), which is a function of temperature and salinity as well as other chemical concentrations; c is the shape factor of the cross-sectional area; F x is the momentum flux caused by eddy viscosity; M S is the external momentum-impulse from artificial sources or sinks; M R is the momentum-impulse from rainfall; M E is the momentum-impulse lost to evapotranspiration; M F is the momentum-impulse from the subsurface by exfiltration; M 1 and M 2 are the momenta-impulses from the overland flow; B is the top width of the cross-section; τ S is the surface shear stress; P is the wet perimeter; τ b is the bottom shear stress, which can be assumed proportional to the flow rate as τ b /ρ = κV 2 where κ = g(n M ) 2 /(R H ) 1/3 , where R H is the hydraulic radius and n M is the Manning roughness.
Equations (1) and (2) are solved using numerical methods, including a hybrid Lagrangian-Eulerian finite element method for dynamical wave model, conventional finite element method for diffusion wave model, and semi-Lagrangian method for kinematic wave model. In this paper, the kinematic wave model in WASH1D is adopted for numerical modelling. In addition, this study focuses on the river routing and bend scour simulation. Compared to upstream watershed, the effect of the rainfall and the exfiltration on river channel flow process is relatively small. For this reason, the source terms relating to rainfall, exfiltration, overland flow, and evapotranspiration are neglected herein. Only the sources of M S and S S are considered as lateral flows.
Three boundary conditions are implemented in the WASH1D model: (1) water depth, (2) river discharge, and (3) rating curves (water stage-discharge relationship). A related study provided a more detailed description of the WASH1D model [20].

SFM2D Model
To simulate typhoon-induced flood flow involving irregular riverbeds, a hydrodynamic SFM2D model was employed in the framework of the cell-centered finite volume method. In this model, the 2D depth-averaged shallow-water governing equations in conservation form are as follows [29,30]: where Q is the vector of conserved variables; F I and G I are the inviscid flux vectors in the x and y directions, respectively; F V and G V are the viscous flux vectors in the x and y directions, respectively; S is the source term; u and v are the depth-averaged velocity components in the x and y directions, respectively; T xx , T xy , and T yy are the depth-averaged turbulent stresses; s 0x and s 0y are the bed slopes in the x and y directions, respectively; and s fx and s fy are the friction slopes in the x and y directions, respectively.
For the divergence theorem and the rotational invariance of the governing equations [30], Equation (3) can be represented by where A cell is the area of the cell; dW cell is the area element; is the control volume; m is the index representing the side of the cell; M is the total number of sides for the cell; T(θ) −1 is the inverse of the rotation matrix corresponding to the m side; θ is the angle between the outward unit vector n and the x axis; n is the outward unit vector normal to the boundary of the control volume; L m is the length of the m side of the cell; T is the vector of conserved variables transformed from the original vector of conserved variables Q, in which u n = u cos θ + v sin θ and v t = v cos θ − u sin θ are respectively the flow velocity components in normal and tangential directions. The SFM2D model involves using appropriate algorithms to solve Equation (8). For the time integration and treatment of the source terms, the fractional splitting technique is used to obtain a well-balanced finite-volume scheme: where H (∆t) and I (∆t) are operators associated with solutions to the homogeneous and inhomogeneous parts, respectively; n is the time index; ∆t is the time step; Q n is the vector of conserved variables for a cell at time index n. Equation (9) is then solved with the Euler method for a homogeneous part and the hydrostatic reconstruction method for source terms. The upstream flux-splitting finite-volume (UFF) scheme proposed by Lai et al. [29] has been verified and is especially suitable for resolving shallow-water free-surface flow involving discontinuities as well as hybrid flow regimes. Therefore, the UFF scheme was utilized to compute numerical inviscid fluxes of mass and momentum. Further detailed descriptions of the SFM2D model have been reported [29,30].

Boundary Connection of 1D and 2D Models
The 2D model is suitable for providing accurate hydraulic resolutions at the river reach containing river structures. When the length of river reach is quite long, the 2D model requires more computational grids and thus is time-consuming. Moreover, if the gauging stations failed, the measured data may not be obtained. Therefore, the measured data could be achieved from 1D numerical model, which is quite suitable for longer river reach with quite less computational time consumed.
The integration of 1D and 2D numerical models is becoming an increasingly efficient approach for creating river flood simulations [31][32][33][34][35]. The flux-based connection method in flow direction proposed by Bladé et al. [33] is employed in this study for its practical advantage. When the 1D river reaches flow into a 2D domain, the flux treatment at the connection boundary is necessary to ensure the mass conservation. As displayed in Figure 1, the variable with subscript L is the known state and the value with subscript R is the unknown solution at the boundary. The variables at initial time step condition represent the known variables (h 2D L , u n  Four conditions are implemented as follows: 1. When the inflow at the connection boundary is subcritical, the unknown unit width discharge (hu n ) 2D R has to be given. The three unknown variables (h 2D R , u n 2D R , and v t 2D R ) can be solved simultaneously from the following three equations [36,37]: in which the unit width discharge (hu n ) 1D R could be given from the simulated results of the WASH1D model. Thus, the value of (hu n ) 1D R at the flow distribution in 2D computation cells can be achieved from the following equation calculations: Q c,ol,or = Q K c,ol,or K c + K ol + K or , K c,ol,or = A c,ol,or where K is the conveyance, and the variables with subscript c, ol, and or denote the values at the main channel, left overbank, and right overbank, respectively.

2.
When the outflow at the connection boundary is subcritical, the unknown water depth h 2D R should be given. Thus, other two unknown variables (u n 2D R and v t 2D R ) can be obtained from the following two relationships: in which water depth h 2D R is obtained from the 1D simulated results (e.g., h 2D When the inflow at the connection boundary is supercritical, all three unknown variables (i.e., h 2D R , u n 2D R and v t 2D R ) should be given from measured data. However, as the field data are not available, one can use the relation (h 2D R = h 1D R ) and employ Equations (11) and (12) to achieve variables u n 2D R and v t 2D R .

4.
When outflow at the connection boundary is supercritical, the transmissive boundary conditions are specified as follows: Employing the integration of 1D and 2D models, the main advantage of the above flux-based connection method is that the numerical fluxes and source terms are preserved to maintain stable solutions. For the numerical modelling, the Courant-Friedrichs-Lewy (CFL) condition is required to ensure numerical stability. Thus, the time step must be limited by [36,37]: where d L,R denotes the distances between the centroids of the two adjacent cells (e.g., L and R).

Empirical Equation for Bend Scour Prediction
Several empirical equations were developed and can be used to estimate bend scour depth [38][39][40]. However, most equations were obtained using laboratory data. Therefore, these laboratory-based equations may not be suitable for predicting bend scouring in the natural river.
The development of more reliable equations using field data are required. Lin and Lin [41] and Lu [42] have presented highly efficient numbered-brick method for directly measuring bend scour depths at embankment sites of Cho-Shui River and Da-Chia River. Furthermore, Hong [21] employed sand pots and wireless tracers to measure the short-term field evolution of bend scour at specific embankments at Da-Chia River and Da-An River. The field data obtained from these measurements enable numerical model verifications. The field data can also be used to develop more reliable empirical equations than those based on laboratory scour data.
To develop a reliable equation for the maximum scour depth in a river bend, Guo et al. [19] collected field data for short-term bend scour in Taiwan and proposed a practical computation equation. More bend scouring field data have become available; thus, the computation equation from Guo et al. [19] can be reconstructed. Hong [21] introduced a concept of effective flood peak duration (tp) to determine the following physical parameters: where ϕ is an unknown function; q is the unit width discharge; d bs is the bend scour depth; S 0 is the channel slope; σ g is the geometric standard deviation of the particle size distribution; is the centerline radius of the bend, where R o is the outer radius of the bend and W is the water surface width. Su and Lu [43] defined that the effective peak flow duration is the time spent from critical flow discharge (Q c ) to peak flow discharge (Q p ) in the rising limb of a flood event. For practical application, Q c can be expressed as the width-full flow discharge. From the field study by Su and Lu [43], the higher duration value tends to cause embankment damage during flood-induced scouring. Thus, the effective peak-flow duration is one of the significant factors in affecting embankment toe failure.
On the basis of field data collected from Cho-Shui River, Da-Chia River, and Da-An River from 11 flood events, the reconstructed scour estimation equation can be expressed as follows [21]: Equation (17) should enable more accurate predictions than laboratory-based equations because it is based on specific field data in Taiwan. In addition, it is easy to apply in 1D and 2D numerical algorithms because of the concept of approach flow [19]. The data ranges for Equation (17) are as follows: S 0 = 0.486%-1.11%, R c = 500-985 m, σ g = 3.28-19.96, tp = 1.5-44 h, q = 8.90-23.88 m 2 /s, and W = 60-450 m, respectively. For modelling bend scour, physical variables expressed in Equation (17) are used to predict bend scour depth. The values of W and q can be obtained from WASH1D and SFM2D models. Others (e.g., σ g , S 0 , R c , tp) are achieved from the field data. Figure 2 shows the methodological flow chart of bend scour simulations. Section 3.1 describes the study site, field data collection, and model setup for WASH1D and SFM2D models. Sections 3.2 and 3.3 present the numerical modelling, including the performance verifications of WASH1D model and SFM2D model. In Section 3.3, comparisons of WASH1D and SFM2D concerning boundary connection are presented. In Section 4, the application of SFM2D model for river embankment failure assessment is given and the proposed assessment methodology is also discussed.

Study Site, Field Data, and Model Setup
This study was conducted in the Da-An River basin, which has a drainage area of 758 km 2 and is located in central Taiwan. As displayed in Figure 3, the riverbed elevation ranges from 0 m to 3877 m. The average bed slope in the main channel reach of Da-An River is approximately 0.011. During typhoons, bend scours may form rapidly because of river flooding, causing severe damage to the river embankment toe. The survey data from the WRA [21] indicated that over three embankment failure incidents have been reported at the Shuiwei Embankment. Therefore, the Shuiwei Embankment was chosen as the study site for bend scour modeling. The 1D hydraulic model is suitable for simulating the main channel processes in a longer river reach. As shown in Figure 3, the selected 1D simulation domain is a 65-km reach of Da-An River between Xiangbi and Taizhong gauging stations. The WASH1D model was employed to simulate the 1D hydraulic flow conditions. The geometric data of 74 cross-sections were collected for model setup.
Because the study area contains the Shuiwei Embankment located in a river bend reach, the SFM2D model was employed to resolve the further detailed 2D flow field and the bend scouring in a 10-km river reach from Zhuolan gauging station to 500 m downstream of the Yili gauging station. To accurately represent riverbeds, a UAS was employed in March 2019 to produce a high-resolution digital elevation model (DEM) with a spatial resolution of 1.0 m for bend scour modeling. As shown in Figure 4, the bed-elevation contour produced from the UAS DEM represents actual river bed topography.  The study area has different types of gauging stations. The measured data from these gauging stations were used for model inputs and verifications. Therefore, different field data were collected, including data on hydraulic structures, river discharge, river water level, tidal stage, bend scour depth, and cross-sectional bed elevation.
Two extreme typhoons, Saola (July 2012) and Soulik (July 2013), were selected for the verification of the WASH1D model. In the 1D simulation, the dynamic river discharge measured at Xiangbi gauging station was set as the upstream boundary and the tidal hydrograph measured at the Taizhong gauging station was set as the downstream boundary. Figure 5a,b display the measured hydrographs in the test events of Saola and Soulik, respectively, as the boundary conditions of simulations. Furthermore, three lateral flows converge onto the main river reach. However, measured data are not available for these lateral locations. Thus, the rainfall runoff model [44,45] was employed to produce the discharge hydrographs of lateral flows. For bend scour modelling, two flood-induced bend scours that occurred in recent years (Typhoon Nesat, July 2017 and Typhoon Maria, July 2018) were selected and simulated. Two types of field data were collected [21,46]: water levels, measured at the Shuiwei water stage gauging station and bend scour depth, measured at the Shuiwei scour gauging station. The measured water levels and bend scour depths were adopted for model verifications. For upstream boundary conditions of 2D simulation, the measured river discharge hydrographs for these two events were available for Zhuolan gauging station. Therefore, using Equation (13) only the water level hydrographs created by the WASH1D model were imposed at the downstream boundary (500 m downstream of the Yili gauging station) of 2D simulation. Figure 6 presents the boundary conditions for the test events of Typhoons Nesat and Maria. Once Equation (17) is used to predict bend scour depth, several physical variables should be given. Table 1 lists the field data collection for short-term bend scour induced by Typhoons Nesat and Maria. From survey field data, it is obviously shown that the bed slope and geometric standard deviation of the particle size distribution did not change significantly during these two floods. In addition, the effective flood peak duration and the centerline radius of the bend did change significantly before and after floods. All of the 1D and 2D simulations were performed on an ASUS computer with Intel ® Core™ i5-8500 CPU and 8.0 GB RAM from Taipei, Taiwan.

Performance Verification of WASH1D Model
Before bend scour modeling, the accuracy of the WASH1D model was verified to ensure the acquired flow solutions at the boundaries of the SFM2D simulation domain were of high quality.
The simulated water levels obtained from the WASH1D model were compared with the measured data, as displayed in Figures 7 and 8, in which the abscissa represents the simulation time and the ordinate denotes the water level. The results reveal that the simulated water levels had acceptable agreement overall with the measurements.  To further evaluate the numerical performances of the WASH1D model for two flood events, the simulated water levels are transferred as simulated water depths using the lowest bed elevation in the cross-section and three criteria are employed [19]: in which Eh p is the peak water depth error; ET p is the error of time to peak water level; R 2 denotes the coefficient of determination; h sim p and h mea p denote the simulated and measured peak water depths, respectively; T sim p and T mea p are the simulated and measured time to the peak water depth, respectively; From the quantitative results listed in Table 2, the simulation performance at Xiangbi station for the Typhoon Saola event had the smallest Eh p (1.64%). The overall performance for Eh p was less than 18.24%. Furthermore, the overall performance regarding the R 2 coefficient for two flood events was larger than 0.75. The quantitative results for the three criteria indicated that the WASH1D model achieved satisfactory overall numerical accuracy in modeling 1D river flow.

Performance Verification of SFM2D Model
Since the bed-elevation data was measured in March 2019, the measured water levels in Typhoons Saola (2012) and Soulik (2013) may not be suitable for calibration due to the river bed changes in 2012 to 2019. Therefore, the measured data in Typhoon Nesat (July 2017) was selected for the verification of the SFM2D model. The dimensions of the SFM2D simulation domain is 10 km in along flow direction and 0.75 km in cross-section direction. The SMS (Surface-water Modelling System) software was used to produce the nonrectangular computational mesh with 7465 elements and 7910 nodes. The minimal and maximum grid size in cross-section direction is 18.5 m and 42.6 m, respectively. The influence of CFL on simulated result was tested and the result indicated that CFL may not affect the simulated water levels significantly. Thus, the CFL is set to be 0.8 herein (i.e., computational time step = 0.1 s).
The main parameter of the SFM2D model which needs calibration is the Manning roughness. The measured water levels at Yili station and Shuiwei scour gauging station in Typhoon Nesat are employed to calibrate the value of the Manning roughness coefficient. Figure 9 shows the influence of the Manning roughness on simulated water levels. The final calibrated value of the Manning roughness is 0.025, resulting in an acceptable match between the simulated and measured water-level hydrographs. Table 3 shows the simulated results using the Manning roughness value of 0.025. The overall performance of the R 2 coefficient was larger than 0.8. The overall performance of the Eh p may smaller than 19%.  Furthermore, the resolution of the adopted computational grid could affect the accuracy of bend scour modeling. Hence, the effect of computational mesh on the simulated result is performed. Two different computational meshes are tested, including coarse (1230 elements and 1364 nodes) and fine meshes (7465 elements and 7910 nodes). The time step is taken as 0.1 s. Figure 10 shows the contour plots of simulated water depths based on two computational meshes for Typhoon Nesat. Employing the fine computational mesh, the simulated water depth approaches to more reasonable resolution. The simulated results and the CPU time are summarized in Table 4. The CPU time consumed by coarse mesh is relatively smaller but the simulated water depth by coarse mesh does not produce accurate resolution, leading to larger errors with Eh p of 57.65% and 34.13% at Yili and Shuiwei scour stations, respectively.

Comparisons of WASH1D and SFM2D Models for Scour Prediction
After the accuracies of the WASH1D and SFM2D models were verified, bend scour modeling was performed. The prediction bend scour depth at the embankment toe generally depends on the approach flow conditions considered with different numerical models. Therefore, the following three approaches are presented and compared: 1.
WASH1D-HR approach, formed using high-resolution bed-elevation data. The high-resolution DEM was used to recreate the cross-sectional topography for the WASH1D model. 3.
SFM2D model based on high-resolution bed-elevation data.
Equation (17) is employed in all three approaches. The boundary connection of Equation (13) is used in the SFM2D model herein for bend scour simulations. In addition, employing Equation (17) for bend scour prediction requires the approach flow conditions (q and W) as inputs. According to the report by Hong [21], the design purpose of Shuiwei water stage gauging station is to provide the approach flow condition. Therefore, the location at Shuiwei water stage gauging station could be approach flow cross section. The numerical performances of WASH1D-CS, WASH1D-HR, and SFM2D in resolving approach flow cross section have to be verified.
For SFM2D model, the computational time step is 0.1 s, the Manning roughness is 0.025, and the fine mesh mentioned in Section 3.3 is employed herein. As shown in Figure 11, the SFM2D model has the best numerical performance of the three presented approaches. Furthermore, the WASH1D-CS approach has the most overprediction. These differences may be partly due to the uncertainty of model resolution. Figure 11. Comparisons of the measured and simulated water levels obtained using three approaches at the Shuiwei gauging station for (a) Typhoon Nesat and (b) Typhoon Maria. Table 5 compares the simulated results with the measured data based on three criteria. For these two flood events, the SFM2D model had the most favorable numerical performance for Eh p , ET p , and R 2 , whereas the WASH1D-CS approach obtained the most erroneous solution. Figure 12 displays the simulated velocity contours under the peak flood condition of Typhoon Nesat obtained using the SFM2D model. The SFM2D model achieves a reasonable velocity component in the main channel, which is faster than those in the floodplain. These results indicate that incorporating high-resolution data into simulations could improve numerical performance. Table 5. Simulated results at the Shuiwei water gauging station obtained using the three presented approaches.

Events
Approaches Three Criteria   Figure 13 shows simulated unit width discharges based on three approaches. For SFM2D model, the output solution of simulated unit width discharge is located at Shuiwei water stage gauging station. However, the numerical framework of WASH1D-CS and WASH1D-HR are not based on the 2D computational grids. Hence, the approach flow cross section shown in Figure 12 is adopted as the output sections of WASH1D-CS and WASH1D-HR. After obtaining unit width discharges, one substitutes the approach flow conditions (q and W) with parameters (σ g , S 0 , R c , tp) into Equation (17) to obtain the bend scour depth at each simulation time step. Comparisons of the measured and simulated bend scour depths over time at the Shuiwei scour observation station for Typhoons Nesat and Maria are presented in Figure 14a,b, respectively. For Typhoon Maria, all three approaches overpredicted results. The differences between the 1D and 2D models may depend on the typhoon event. To further investigate the differences between the 1D and 2D models, an error norm was employed, defined as follows: where (d bs ) sim i and (d bs ) mea i are the simulated and measured bend scour depths at each simulation time, respectively. Among the three approaches tested, the results displayed in Table 6 demonstrated that the SFM2D was the most preferable because of its optimal numerical performance. In addition, the CPU time consumed by WASH1D and SFM2D for Typhoon Maria are 0.042 h and 1.27 h, respectively.

Model Application for River Embankment Failure Assessment
In disaster-prone regions, accurate and realizable predictions of the riverbank failure risk posed by scouring are vital to mitigating the effects of disasters. However, riverbank failures caused by scouring often occur suddenly and without warning. Directly monitoring the safety of embankment during floods is very difficult. Several numerical models can be employed to predict the risk of failure; however, these models require the inputting of certain parameters as well as the preprocessing of the model setup. Therefore, an alternative and efficient methodology would be valuable for rapidly assessing the risk of embankment failure prior to field site investigation.
The simulated results presented in Section 3 demonstrate that the proposed SFM2D model with fine mesh achieves acceptable performance for bend scour simulations. The model was applied to investigate the influence of river discharge or water level on the risk of riverbank failure. On the basis of the analysis results, an assessment methodology without calibrated parameters was proposed. The procedures of this new assessment methodology consist of two approaches, as displayed in Figure 15.

New Assessment Methodology: Approach 1
Approaches to river embankment toe protection are crucial for preventing embankment toe scouring. Once typhoon-induced floods cause the failure of embankment toe protections, the embankment toe becomes directly damaged by scouring. If the embankment toe is eroded until the scoured bed level is smaller than the riverbed level of the embankment toe foundation, the embankment toe could fail. Therefore, embankment toe protections must be considered in riverbank failure assessments.
When the shear stress acting on riverbank toe protections exceeds critical shear stress, the riverbank toe protections erode and may fail. The shear stress considers the actual force of the river flow on the embankment toe foundation. Hence, it is a more reasonable predictor of scour potential than river velocity. According to non-uniform flow theory, the shear stress can be simplified as follows [15]: (20) in which the water depth and friction slopes in the x and y directions in every computational grid are simulated using SFM2D model. Figure 16 illustrates the river shear stress contour plot under the peak flow condition for Typhoon Maria. The results reveal that the simulated value of shear stress at the Shuiwei Embankment under peak-flood condition is 345.91 N/m 2 , which is smaller than the critical value of 588.6 N/m 2 [47]. Thus, the embankment toe protections were stable. To further investigate the relationships between shear stress and river discharge, four discharge values for 1.11-, 2-, 5-, and 10-year return period floods with peak discharges of 950, 2630, 4780, and 6490 m 3 /s were employed to conduct scenario simulations. As displayed in Figure 17, the simulated results of shear stress and the river discharge for four scenarios were combined with the field data under two events, revealing the relationships between river discharge, shear stress, and bend scour depth. In Figure 17, two field data of bend scour depth for Typhoons Nesat and Maria are considered. It is noted that the field data for shear stress is not available, thus the shear stress for Typhoons Nesat and Maria and four scenario cases are obtained from the simulated results under peak-flood condition. On the basis of this relationship curve, the critical flow discharge at the embankment toe protection works was estimated to be 2800 m 3 /s for a critical shear stress of 588.6 N/m 2 . Furthermore, the warning value of river discharge for riverbank failure was selected based on the failure event during Typhoon Soulik, which caused serious damage to the Shuiwei Embankment under peak flow with a discharge of 3120 m 3 /s. Approach 1 in the proposed assessment methodology is principally based on river discharge. The two threshold values of river discharge detailed in Figure 15 were used as the indicators to aid with the risk identification of riverbank failure. If river discharge is larger than Q c1 of 2800 m 3 /s, the embankment toe protection works may fail. Therefore, without protection, the river embankment is at risk. Furthermore, sustained river discharge higher than a Q c2 value of 3150 m 3 /s may lead to failure of the riverbank toe foundation.

New Assessment Methodology: Approach 2
Water depth data are more easily obtained than river discharge. Therefore, the second approach proposed employs water depth as the input. The aim of this study was to develop two deterministic equations that could be considered useful tools for estimating shear stress and bend scour depth, respectively. Because the collected field data are limited, the simulated results from four numerical experiments in Section 4.1 were integrated into the model as well as field data of two events. As described in Section 4.1, the shear stress for Typhoons Nesat and Maria is not available, thus is obtained from the proposed model at peak flow simulation. In addition, the water depth h u in the upstream straight reach is a suitable information as input that can be obtained from approach flow cross section, as shown in Figure 12.
Employing the four simulated and two field values comprised of h u , R c , and shear stress, a predictive regression equation to determine shear stress based on water depth is proposed as follows.
Equation (21) provides acceptable results regarding measured versus computed data (R 2 value of 0.95), as illustrated in Figure 18. Applying the four simulated and two field data comprised of h u , R c , W, and bend scour depth, the local scour depth, created by river flood flow in a bend reach, can be proposed as follows: The results detailed in Figure 19 demonstrate that the predictions agree well with the integrated data, producing an R 2 value of 0.87. Following Equations (21) and (22), the procedure for approach 2 to assess riverbank failure is as follows: 1.
Combine Equations (21) and (22) with water depth values to estimate the shear stress and bend scour depth.

2.
The estimated shear stress can then be compared with the critical value.

3.
If the estimated shear stress is larger than the critical value of 588.6 N/m 2 , the riverbank toe protection works may fail.

4.
The estimated bend scour depth can be compared with the designed value. If the estimated bend scour depth is larger than the designed value, the riverbank may fail.
In summary, the proposed methodology incorporates the integration of the available field data and numerical results from the SFM2D model. The methodology avoids the use of calibrated parameters by only inputting flow conditions (e.g., river discharge or water depth). Therefore, it can be used to provide a rapid assessment of the risk of embankment failure. The methodology is expected to provide an effective strategy to assist river managers in decision-making.

Conclusions
A new methodology is proposed that can be used to assess the risk of river embankment failure caused by bend scouring prior to field site investigations. The proposed bend scour simulation model integrates the WASH1D and SFM2D models with high-resolution bed elevation data. The proposed scour model was applied to simulate bend scour at the Shuiwei Embankment in Da-An River, Taiwan, and was verified against the field data from Typhoons Nesat and Maria. Furthermore, two empirical equations were proposed for predicting shear stress and bend scour depth. A practical assessment methodology was therefore developed on the basis of these two equations. The methodology involves considering the effects of shear stress and bend scour depth on the risk of riverbank failure.
The main conclusions and contributions of this paper are as follows. 1.
The integrations of WASH1D, SFM2D, and bend scouring computation equations create a bend scour simulation model. The applications of the integrated model indicated that the model was able to successfully resolve the temporal process of bend scour caused by typhoon-induce floods.

2.
Based on high-resolution bed-elevation data, three approaches (namely WASH1D-CS, WASH1D-HR, and SFM2D) were presented and compared. The results demonstrate that the SFM2D model produces reasonable scour predictions with acceptable accuracy for water level and bend scour depth.

3.
A new methodology for assessing embankment failure is proposed. The novel aspects of the proposed methodology include the consideration of toe protections and two new predictive equations. The former is implemented to enhance the reliability of riverbank failure assessments. The latter provides more specific deterministic results.
Due to the limitation of computing resources, the resolution of the used computational mesh is much lower than that of bed-elevation data. However, from the simulated results presented in Figures 10b and 12, the proposed model achieves the reasonable water depth and velocity solutions in the main channel and floodplain. As shown in Tables 5 and 6, the analyses of three criteria from simulated results indicate that the numerical performance of the model is acceptable. The more detailed test on the sensitivity analysis of computational mesh requires further study.
The advantage regarding the proposed methodology in terms of applicability is that it can easily be implemented in other river reaches. The method is particularly helpful for river managers in Taiwan and could be integrated into monitoring systems.