Improving the 2D Numerical Simulations on Local Scour Hole around Spur Dikes

: Local scour is a common threat to structures such as bridge piers, abutments, and dikes that are constructed on natural rivers. To reduce the risk of foundation failure, the understanding of local scour phenomenon around hydraulic structures is important. The well-predicted scour depth can be used as a reference for structural foundation design and river management. Numerical simulation is relatively efﬁcient at studying these issues. Currently, two-dimensional (2D) mobile-bed models are widely used for river engineering. However, a common 2D model is inadequate for solving the three-dimensional (3D) ﬂow ﬁeld and local scour phenomenon because of the depth-averaged hypothesis. This causes the predicted scour depth to often be underestimated. In this study, a repose angle formula and bed geometry adjustment mechanism are integrated into a 2D mobile-bed model to improve the numerical simulation of local scour holes around structures. Comparison of the calculated and measured bed variation data reveals that a numerical model involving the improvement technique can predict the geometry of a local scour hole around spur dikes with reasonable accuracy and reliability.


Introduction
Steep slopes and severe bed changes during floods are the common characteristics of the rivers in Taiwan. Floods often cause channel incision, bank erosion, and meandering migration and endanger the safety of structures or river protection works. Spur dikes are protections commonly used for river engineering to lead water flow and reduce bank erosion. However, the local scour around spur dikes influences the stability of the dikes' foundation. To reduce the risk of foundation failure, the understanding of local scour phenomenon around hydraulic structures in rivers is important.
One can understand the formation mechanism of scour hole by the flow field and sediment movement behavior of a single cylinder pier. Figure 1 displays the scouring process around a pier structure to the front of the scour hole. Generally, the sediment entrainment and transport around the scour hole were most related to the turbulence in the approach flow. Through strong vertical flows, the turbulence fluctuations in the water are brought into the sediment layer. The turbulence fluctuation near the sediment layer enhanced the ability of the water to entrain sediment. When scour hole occurs gradually around the hydraulic structure, turbulence fluctuation in the approach flow contributes to the scouring process where there is sediment particles' pick up, trajectory, deposition, and bed evolution processes with a mixed non-equilibrium sediment transport. Moreover, the characteristics of sediment, bed load with the gravity projection, and the bed geometry of a slope have crucial effects on incipient sediment motion. The sliding effect of scour hole occurs when a local bed slope is higher than the angle of repose. The angle of repose is the steepest angle of descent relative to the horizontal bed to which the bed material can be piled without slumping. At this angle, the bed material on the slope face is on the verge of sliding. Figure 1. Scouring process around a pier structure. A laboratory experiment was conducted in the Hydraulics Laboratory of the Department of Civil Engineering, Dokuz Eylul University, Turkey [1]. (a) Initial state before the scour occurred. (b) Local scour hole that occurred around the pier. (c) Occurrence of the sliding effect in the front of the scour hole when the local bed slope was higher than the angle of repose.
Most mobile-bed models proposed in the past few decades have focused on the sediment transport of an alluvial channel. Jia et al. [2,3] developed a two-dimensional (2D) mobile-bed model, known as CCHE2D numerical model, which is a numerical system for modeling 2D non-steady turbulent river flow, sediment transport, and water quality prediction. The model was also designed for use in multiprocess simulation of surface water with complex geometry, such as riverbed geometry; of bank erosion with both uniform and non-uniform sediment; and of meander river migration. The model was validated through many physical experiments and field cases in the development process. The CCHE2D numerical model can be used for river engineering and design of hydraulic structures, such as grade control structures and dikes. In the review of past research on mobile-bed models, Lai et al. [4] proposed a 2D model, known as SRH-2D numerical model that was employed to simulate the bed evolution that occurred downstream of the Chi-Chi weir in the Choushui River, Taiwan. Engineering plans pertaining to scour prevention and treatment were proposed and simulated using the model. Liao et al. [5] developed a 2D mobile-bed model, known as EFA2D numerical model that considers the non-equilibrium suspended sediment concentration profile with a bedrock erosion mechanism for natural rivers.
In the review of the research on the mechanism of modifying sediment transport behavior, Juez et al. [6][7][8] discuss the numerical assessment of bedload formulations for transient flow. The dam break flows over dry/wet initial conditions and erodible channel of experimental cases are simulated by considering the improved model, called Smart CFBS (Smart Combined Friction and Bed Slope model). The study has allowed a detailed analysis of the relative behavior to be performed in 2D situations of different sediment discharge formulas that were derived from 1D laboratory cases, and use of the Smart CFBS formula is suggested, regardless of the hydro-morphodynamic situation. The bedload transport over steep slopes with gravity projection has been studied and embedded in a non-trivial way, not only in the sediment transport rate but also in hydrodynamic equations. It is noted that the gravity effects play an important role on the sediment transport magnitude.
Horvat et al. [9] developed a 2D model that couples the active layer and multiple size-class approaches for modeling sediment transport by using an enhanced advection algorithm. The focus of their study was improvement of advection computation from the point of view of numerical methods by reducing the limitation of simulation time. The study proposed the advection scheme that allows for larger time steps in simulation. This restriction was overcome by introducing modifications to the characteristic method. The model was assessed and validated using field measurements conducted on the Danube River.
The 3D mobile-bed model has the capability to simulate the local scouring process. Jia and Wang [10] simulated 3D flow in a preformed scour hole. A detailed comparison of their simulation results and measurements were presented. Nagata et al. [11] developed a 3D model to simulate flow and bed deformation around river hydraulic structures. The model solved the fully 3D, Reynolds-averaged Navier-Stokes equation that was expressed in a moving boundary-fitted coordinate system to calculate the flow field of water and bed surfaces at various times. A non-linear k-ε turbulence model was employed to predict the flow near the structure, where 3D flow was dominant. Burkow and Griebel [12] conducted a numerical simulation of a fully 3D fluid flow to reproduce the current-driven sediment transport processes. The scouring at a rectangular obstacle was investigated, and the typical sedimentary processes and the sedimentary form of a scour mark were well captured by the simulation.
Castillo and Carrillo [13] studied the scour produced by spillways and outlets of a dam due to the operation of free surface spillways and half-height outlets by using three complementary procedures: obtaining empirical formulas using models and prototypes, using a semiempirical methodology based on the pressure fluctuation-erodibility index, and employing computational fluid dynamics simulations. Jia et al. [14] proposed a local scour scheme for a 3D mobile-bed model and successfully applied the model to the case of a bridge pier under non-uniform sediment conditions. The sediment entrainment, transport, and turbulence around a bridge pier are simulated. The turbulence fluctuations in the water are brought into the sediment layer through strong vertical flows. The ability of the flow to entrain sediment is enhanced by the intruding turbulence fluctuation near the sediment layer. This enhancement could be considered in 3D model by using additional shear stress related to the intruding turbulence fluctuation. However, most 3D models are still limited by the computational efficiency of the model in filed cases.
The sediment entrainment and transport around the scour hole are related to the turbulence fluctuation, characteristics of sediment, bed-load with the gravity projection, and the bed geometry of a slope that have crucial effects on incipient sediment motion. Yang et al. [15] investigated the angle of repose of non-uniform sediment by considering the concept of exposure degree to account for the hidden and exposed mechanism of nonuniform sediment transport. Meng et al. [16] proposed the angle of repose for uniform and non-uniform sediment. The angle of repose is based on the compacted friction relationship between mud and sand particles. Sediment dynamics holds that, when water flows over sediment particles on a loose surface, intersediment collision resists the fluid shear stress. This is consistent with underwater flow and sediment dynamics. Al-Hashemi and Al-Amoudi [17] reviewed the angle of repose of granular materials. This angle is an essential parameter for understanding the microbehavior of granular material and then relating this behavior to the material's macrobehavior.
In general, both an empirical formula and a numerical model can be used to evaluate the equilibrium scour depth around a hydraulic structure. Most of the empirical formulas have limitations to calculate the scouring depth or hole shape changing over time. This phenomenon and mechanism must be analyzed using the 2D or 3D numerical models. Although the empirical formula is relatively convenient to use, the numerical model is often complicated to use in terms of practical applications, especially when it is a 3D mobile-bed model. Therefore, obtaining a balance between these two approaches is crucial.
In this study, an empirical formula for calculating the sediment repose angle [16] and the bed geometry adjustment mechanism are integrated into the CCHE2D numerical model to improve the simulation of the local scour hole around a hydraulic structure. The model is calibrated and validated by simulating an experimental case [18] involving multiple spur dikes. Data are collected on flow depth, bed material grain size, scour depth around structures, and scour hole shape. These basic data are used to calibrate and validate the numerical model. By improving the 2D numerical simulations on local scour hole around spur dike, it will provide practical application value in terms of model efficiency and speed.

Numerical Model
The CCHE2D numerical model is a general surface water flow model. It simulates dynamic processes of flow and sediment transport in natural rivers. The CCHE2D numerical model has been developed at the National Center for Computational Hydroscience and Engineering (NCCHE), the University of Mississippi, for more than 20 years, and the developers continue to improve the model in response to different real problems. The CCHE2D numerical model comprises flow and sediment transport modules. The governing equations of the flow module are 2D depth-integrated Reynolds equations. The free surface elevation of the flow is calculated using the depth-integrated continuity equation. Two methods for calculating the eddy viscosity exist in the model: the depthintegrated parabolic eddy viscosity formula and depth-integrated mixing length eddy viscosity model [19,20]. The sediment transport module comprises the suspended sediment transport and bed load transport simulations for non-uniform sediment. More details of the flow and sediment transport modules can be obtained in the CCHE2D technical report and verification and validation test documentation [2,3].

Flow Module
The governing equations of the flow are 2D depth-integrated Reynolds equations expressed in the Cartesian coordinate system and are as follows: where u and v are the depth-integrated velocity components in the x and y directions, respectively; t is the time; g is the gravitational acceleration; η is the water surface elevation; ρ is the density of water; h is the local water depth; f Cor is the Coriolis parameter; τ xx , τ xy , τ yx , and τ yy are the depth-integrated Reynolds stresses, shown as Figure 2; and τ bx and τ by are the shear stresses on the bed and flow interface. The shear stress terms are not considered at the water surface because the wind shear driven effect is small and not considered in the model. The turbulence Reynolds stresses in the governing equations are approximated according to the Bousinesq's assumption that are related to the main rate of the strains of the depth-averaged flow field with a coefficient of eddy viscosity: Free surface elevation of the flow was calculated using the depth-integrated continuity equation: ∂h ∂t Two methods for calculating eddy viscosity exist in the model: the depth-integrated parabolic eddy viscosity formula and depth-integrated mixing length eddy viscosity model. First, the eddy viscosity coefficient ν t is calculated using the depth integrated parabolic eddy viscosity formula: where C s is the integration constant, u * is shear velocity, κ is the von Karman's constant (0.41) and ς is the relative depth of the flow. A xy is a coefficient to adjust the value of the eddy viscosity. Its default value is set to 1 and it can be adjusted by users from 1-10.
In addition to this approach, a depth integrated mixing length eddy viscosity model is also available: where the depth integrated velocity gradient along vertical coordinate ∂U/∂z is introduced to account for the effect of turbulence generated from the bed surface. The eddy viscosity defined by Equation (7) would be zero in the uniform flow condition without this term. It is determined in the way that eddy viscosity shall be the same as that of the uniform flow in the absence of other terms. Assuming that the flow is of logarithmic profile along the depth of the water, the total velocity U of the vertical gradient should be:

Sediment Transport Module
The sediment transport module performs suspended sediment transport and bed load transport simulations for non-uniform sediment. The module is also designed to simulate a channel bed with large slopes and curved channel secondary flow effects. The capability of simulating the non-equilibrium sediment transport is achieved by solving the following equations.
Suspended load: Bed load: Bed changes: where ε s is the eddy diffusivity of sediment in the vertical z direction; C k is the concentration of the k-th size class, and C *k is the corresponding transport capacity; α s is the adaptation coefficient for suspended load; ω sk is the sediment settling velocity; δ b is the thickness of bed layer; c bk is sediment load of the k-th size class; q b*k , q *k , q bxk and q byk are the bed load transport capacity, the bed load transport rate, and transport rate components in x and y directions, respectively; L t is the adaptation length for bed load that characterizes the distance for a sediment process adjusting from a non-equilibrium state to an equilibrium state, and the model computes the non-equilibrium sediment transport process including the actual and equilibrium sediment load by considering the adaption length and factors; p' is the porosity of bed material, and Z bk is the bed change.
In most of the early sediment transport models, it is usually assumed that the actual sediment transport rate or a saturated concentration is equal to the capacity of flow carrying sediment at equilibrium conditions. Based on this sediment transport state, a sediment transport model is called an equilibrium transport model. It could be a steady exchange of particles between bed material and the sediment load. However, this equilibrium assumption is not realistic in natural alluvial rivers because of variations in flow conditions and bed material properties. Non-equilibrium sediment transport models adapt sediment transport equations to determine the realistic bed load and suspended load transport rates.
In the processes of vertical integration, the water must be assumed to be shallow relative to channel width and the vertical variation in the flow to be negligible; otherwise, the dispersion terms must be maintained to preserve the effect of vertical sediment profiles on changes in the bed form. In the second case, the source term is nonzero but represents the dispersion terms. Computation of the dispersion term is complicated and requires knowledge of the vertical velocity and suspended sediment profiles in the x and y directions. The dispersion terms in suspended load transport equation are usually combined with the diffusion terms. More details can be found in the technical report of the sediment module [2]. Figure 3 displays the modules and flowchart of the calculation process in this study. First, the flow module calculates the flow conditions, e.g., flow field, velocity, and water depth. After calculating the flow conditions, the sediment transport module calculates the sediment transport process, e.g., sediment transport capacity, bed material sorting, and bed changes. Finally, local scour holes around hydraulic structures are adjusted using the bed geometry adjustment mechanism proposed in this study. The flow velocity profile in the vertical direction is generally non-uniform. The predicted local scour depth and volume are often underestimated due to the hypothesis of the depth-averaged velocity in the 2D model. In this study, an empirical formula for calculating the sediment repose angle [16] and an algorithm for modifying the bed geometry in the local scour hole are integrated into the 2D mobile-bed model to improve the model's ability to simulate local scours around the structures. Experiments were conducted by Meng et al. [16] by using cohesionless heavy sediment to obtain the angle of repose by using the natural accumulation method and angle of surface friction by using the tilting method. The experimental results revealed that the angle of repose over the water is bigger than that in the submarine condition. Moreover, the angle of repose does not increase monotonously with an increase in the sorting coefficient but increases with a decrease in the asymmetry coefficient S of non-uniform sediment for the same median diameter D 50 .

Bed Geometry Adjustment Mechanism for a Local Scour Hole
The empirical formula (expressed in Equations (13)-(15)) proposed by Meng et al. [16] for uniform and non-uniform sediment is employed in this study. The procedure of the bed geometry adjustment mechanism for local scour holes can be explained as follows.
First, the angle of repose is calculated according to the sorting coefficient and bed material (e.g., uniform sand, non-uniform sand, gravel, or cobble). Then, the bed elevation around the hydraulic structures in the stream is calculated and modified using trigonometry with different bed slope trends (e.g., downhill or uphill slope). Finally, the adjusted bed geometry of the local scour hole around the hydraulic structures is obtained.
The angle of repose for uniform and non-uniform sediment proposed by Meng et al. [16] can be expressed as follows: S 0 = D 75 /D 25 (15) where φ i is the angle of repose at computational node i and D i is the mean diameter of the bed material at computational node i. Moreover, S 0 is the sorting coefficient, in which D 75 and D 25 are the particle sizes at 75% and 25% weight of the bed materials. The average angle of repose between computational nodes i and i + 1 can be expressed as follows: Depending on whether there is a downhill or uphill bed slope around the scour hole (Figure 4), the adjusted bed geometry can be expressed as follows: where Z i is the bed elevation at computational node i that is modified by the angle of repose and ∆x is the distance between computational nodes i and i + 1. In the computational nodes of bed elevation in the 2D model, where the set of i are the nodes along the longitudinal direction (x-direction) and the set of j are those along the lateral direction (y-direction). Figure 5 shows a flowchart of the bed geometry adjustment mechanism for a local scour hole. According to the characteristics of sediment classes, the bed geometry of local scour hole could be adjusted in the 2D model by considering the adjustment mechanism mentioned above.

Initial and Boundary Conditions for the Experimental Case
An experimental case of multiple spur dikes is used for model calibration and validation. Figure 6 Table 2 shows the parameters selected in the model calibration. In the computational mesh, the total grids are 21,231 where the minimum lateral and longitudinal mesh sizes around the spur dikes are 0.5 cm. In the flow module, the computational time step is 0.1 s, and the equilibrium total simulation time is about 18,000 s. The mixing length eddy viscosity model is selected. In the case of a spur dike, the distance to the wall should be used as the length scale instead of that to the bed. Otherwise, the depth integrated coefficients for the eddy viscosity would be too large when the interior nodes are close to the wall. According to the CCHE2D technical report, the turbulence viscosity coefficient is a coefficient to adjust the value of the eddy viscosity. Its default value is set to 1 and it can be adjusted by users from 1 to 10. Another important issue regarding the mixing length model is the wall effect. Very close to the wall, the distance to the wall should be used as the length scale. Otherwise, the depth integrated coefficients for the eddy viscosity would be too large when the interior nodes are close to the wall. In CCHE2D model, the normal distance from a node to the wall is used to calculate the mixing length. It is also used to calculate the parabolic profile. This approach avoids the prediction of very large eddy viscosity near the wall. More details of the turbulence model and wall effect modification can be obtained in the CCHE2D technical report [2]. In the sediment transport module, the sediment transport formula of Wu et al. [21] is selected for the model calibration. It related the bed-load transport rate to the nondimensional excess grain shear stress with critical shear stress for incipient motion and grain shear stress. The variation of bed material gradation in the mixing layer can be described by the partial differential equation, e.g., Equations (10) and (11), while in other layers under the mixing layer, the bed material gradation can be determined by using the mass conservation law. The bed material is often divided into multiple layers, and the top layer is the mixing layer. The variation of bed material gradation is subject to exchange with those moving with flow. The mixing layer thickness is related to sand dune height or particle diameters of bed material. The minimum mixing layer thickness is set to 0.01 m in this study.

Parameters for the Case Study
To calculate the non-equilibrium sediment transport, the adaptation length for bed load and adaptation factor for suspended load are considered, respectively, in the model. In general, the adaptation length characterizes the distance for sediment to adjust from a non-equilibrium state to an equilibrium state. It is a length scale for the river bed to respond to the disturbance of the environment, such as hydraulic structures, channel geometry changes and incoming sediment variation. The adaptation factor for suspended load is a use-specified parameter which is related to water depth, bed shear velocity, and settling velocity of sediment particles. In this study, the adaptation length for bed load transport is set to 0.4 m, and the adaptation factor for suspended load is set to 0.5 which have the best results. Figure 8 shows the bed material grain size distribution of the experimental case. The median diameter is 0.42 m, and the specific gravity is 2.5. Three bed material size classes are selected in the model based on the measured data of experimental case.

Flowchart of Model Calibration and Validation
The selection of flow and sediment parameters in the model calibration and validation are based on the sensitivity test. The sensitive parameters would be selected for the model calibration. Figure 9 shows the flowchart of model calibration and validation using the original model. First, it is needed to setup the model parameters, initial bed elevation, initial water depth, and boundary conditions. Then, the steady flow conditions on fixed bed are calculated and compared with the measured water stages. Under the standard of error analysis, the time step, turbulence model, turbulence viscosity coefficient, and wall slipness coefficient of flow module should be adjusted and calibrated by comparing the simulated results with the measured data. After the calibration of flow module, the sediment transport module on mobile-bed would be computed and compared with the measured bed changes. The sediment parameters (e.g., sediment transport formula, minimum mixing layer thickness, adaptation length for bed load transport, adaptation factor for suspended load) should be adjusted and calibrated. Finally, to complete the model validation, another set of the same experimental case is carried out based on the same model parameters. Figures 10 and 11 show an overview of the simulated flow fields and velocity magnitude around the multiple spur dikes. The whole velocity magnitude of case G5a is higher than case G4a. Due to the setup of spur dikes in the experiment, there are contraction effects near the region, and these effects cause the flow velocity to increase. Moreover, eddy areas are present behind the dikes. The reattachment point shows that the divided streamline reattaches to the dike. When the flow passes through multiple spur dikes, a uniform flow is gradually obtained because of the increased width of the channel. Parts of the water flow into the spur dike field then make the circulation. The transported sediment causes deposition in this region around the dikes. Figures 12 and 13 show the simulated and measured non-dimensional bed changes around multiple spur dikes in the original model when the equilibrium time is about 18,000 s. These changes are divided on the basis of water depth to compare the simulated and measured non-dimensional results. There are deposition bars occurring in the fields of multiple spur dikes, and the erosion occurs around the head of the first dike and the last dike. There is general erosion occurring along the longitudinal direction of the experimental flume. The location and magnitude of maximum erosion depth are predicted accurately in both experimental simulations. The calculated depths are 2.03 and 4.03 cm for the G4a and G5a cases. Moreover, the measured values are 2.00 and 4.19 cm. However, the angle of repose and modified bed geometry for the scour hole are not yet considered, and the predicted scour location is biased.

Simulated Bed Changes in the Original and Adjusted Model
Based on the parameters of model calibration and validation, the improved bed geometry adjustment is considered and compared with the original model. Figure 14 displays an overview of the simulated non-dimensional bed changes around the multiple spur dikes when considering the angle of repose and bed geometry adjustment mechanism proposed in this study. According to the simulation results obtained using the adjusted model, significant local scouring is observed around the head of the first spur dike. The scour hole shape around the first spur dike could be simulated more obviously.   6. Discussion and Conclusions 6.1. Discussion Figure 16 displays a sketch of a local scour hole around a spur dike; the sketch is used to compare the measured and simulated local scour hole shape [22]. Here, L is the distance of a vertical spur dike to a flume wall, and L sa , L sb , and L sc are the distance from the head of the spur dike to the upstream, lateral, and downstream sides of the scour hole. The nondimensional distance is calculated to compare the modeling and measured scour hole shape results. Table 3 presents a comparison of the non-dimensional distance of the scour hole around the spur dike. The simulated results reveal that the modified model makes 21.43% more accurate predictions than the original model, especially for the non-dimensional distance of L sa /L in case G4a. This is discovered because the sediment repose angle and bed geometry adjustment mechanism are considered in the model. The simulated nondimensional distances from the head of the spur dike to the lateral and downstream sides of the scour hole are in agreement with the measured data. The proposed mechanism provides an improved modification for 2D simulation of the area around hydraulic structures. Figure 16. Sketch of a local scour hole around a spur dike. Modified from a study by Copeland [22]. The longitudinal bed profiles in the experimental flume from A-Line to E-Line are compared with the original and adjusted modeling results to investigate the applicability and accuracy of the proposed bed geometry adjustment mechanism. It shows the relative scour trends around the spur dikes. The Mean Relative Error (MRE) and Root Mean Square Error (RMSE) are used for error analysis to evaluate the improvement obtained using the adjusted model. The MRE and RMSE can be expressed as follows: where Z mea i and Z simu i are the measured and simulated bed changes at computational node i, respectively, and n is the total number of nodes measured along the longitudinal bed profile. Figure 17 displays the measured and simulated bed elevation profiles along the longitudinal channel. The simulated results and measured data exhibit similar trends. There is a local scour near the head and body of the spur dikes. Due to the contraction effects on the hydraulic structures, general scour holes occur in the middle reach of the channel. The sediment particle size of the scour hole is gradually increased during the scouring process. Compared with the original model, the adjusted model determines a more accurate local scour hole. The improvement scope of adjusted model is decreasing from the head of spur dike to the other side of channel.  Table 4 compares the MREs and RMSEs along the longitudinal bed profiles. In case G4a, the simulated bed elevations obtained using the adjusted model are in favorable agreement with the measured data along the A-line and B-line. Because of the increase in the inflow discharge of experiment, the scale of local scouring is more obvious in case G5a ( Figure 18). The adjusted model provides a superior simulation of scour hole shape with more reasonable accuracy and reliability from the A-line to the D-line. The E-line has exceeded the scope of influence of local scour hole so that the original and adjusted model have the same trends.  The empirical formulas proposed by Liu et al., Gill,, which are three common empirical formulas, are selected to calculate the head scouring depth of the first spur dike and compare the measured results with the modeling results. Table 5 presents a comparison of the empirical formula results and modeling results. The measured head scouring depth of the first spur dike is approximately 2 cm, and the depth in the modified model is approximately 2.25 cm. These values are in agreement with the measured data values. The original model calculates an underestimated head scouring depth of approximately 0.41 cm due to the simplification of 2D simulation around the structure. The calculation result obtained using each scour depth empirical formula has a large error compared with the measured scour depth. The proposed adjusted model can thus predict local scour around spur dikes with more reasonable accuracy.

Conclusions
In the research of hydraulics and river dynamics, the understanding of the local scour phenomenon by the streams or the flood around hydraulic structures is important to reduce the risk of foundation structure failure. The empirical formulas are usually convenient to use for the calculation of maximum scour depth in the upstream or downstream of the structures. However, most of the empirical formulas have limitations to calculate the scouring depth or hole shape changing over time. This phenomenon and mechanism must be analyzed using the 2D or 3D numerical models.
In this study, a repose angle formula and bed geometry adjustment mechanism are integrated into a 2D mobile-bed model to improve the model's simulation accuracy for local scours around structures. The sediment classes of uniform sand, non-uniform sand, gravel, and cobble are considered in the proposed model. A multiple spur dike experimental case is simulated for the model calibration and validation. By adjusting the bed geometry of a scour hole depending on whether it is on a downhill or uphill slope, scour depth and geometry of local scour hole can be obtained that are in agreement with the measured situation. The Mean Relative Error (MRE) and Root Mean Square Error (RMSE) were calculated to compare the model results. In the validation case (G5a), the simulated results of the bed profile through the spur dike (A-Line) by adjusted model have a significant improvement effect. The improvement scope is decreasing from the head of spur dike to the side of channel.
According to the simulated results of the multiple spur dike experimental case, the sediment particle size of the scour hole is gradually increased during the scouring process. The scour hole side slope is steep at upstream of the spur dike, it is most affected by the sediment repose angle. In the comparison of the scour hole shape with the measured data, the modified model is 21.43% more accurate than the original model. A comparison of the calculated bed changes obtained using an empirical formula with measured bed change data reveals that the proposed numerical model, which includes the sediment repose angle and bed geometry adjustment mechanism, can predict local scours around spur dikes with more reasonable accuracy. Acknowledgments: The writers would like to acknowledge the Ministry of Science and Technology, Taiwan, that has supported this study. The Water Resources Agency, Ministry of Economic Affairs, Taiwan, has provided valuable field data and technical assistance.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
A xy turbulence viscosity coefficient B channel width C depth-integrated sediment concentration C k concentration of the k-th size class C *k concentration of the k-th size class corresponding transport capacity C s integration constant c bk sediment load of the k-th size class D 25 particle size at which 25% by weight of the bed materials is finer D 50 mean diameter D 75 particle size at which 75% by weight of the bed materials is finer D i mean diameter of bed material at computational node i E scour depth f Cor Coriolis parameter F r Froude number g gravitational acceleration h water depth i computational node k number of sediment size class L distance of spur dike vertical to flume wall l depth integrated mixing length L p length of dike L sa distance from the head of spur dike to the upstream side of scour hole L sb distance from the head of spur dike to the lateral side of scour hole L sc distance from the head of spur dike to the downstream side of scour hole L t adaptation length for bed load n total number of nodes along the longitudinal bed profile p' porosity of bed material q discharge per unit width q b*k bed load transport capacity q *k bed load transport rate q bxk bed load transport rate components in x direction q byk bed load transport rate components in y direction Q discharge S downhill or uphill bed slope between two nodes calculated by the bed geometry t time in governing equations U total velocity u depth-integrated velocity component in x direction u * shear velocity v depth-integrated velocity component in y direction V averaged velocity x direction along the longitudinal bed profile y direction along the lateral bed profile Z bk bed change