Next Article in Journal
A Novel Thin-Film Technique to Improve Accuracy of Fluorescence-Based Estimates for Periphytic Biofilms
Next Article in Special Issue
Sediment Transport and Water Flow Resistance in Alluvial River Channels: Modified Model of Transport of Non-Uniform Grain-Size Sediments
Previous Article in Journal
Mobile Genetic Elements Drive the Antibiotic Resistome Alteration in Freshwater Shrimp Aquaculture
Previous Article in Special Issue
Experimental Study on the Influence of the Transition Section in the Middle of a Continuous Bend on the Correlation of Flow Movement in the Front and Back Bends
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
Disaster Prevention & Water Environment Research Center, National Yang Ming Chiao Tung University, 1001 University Road, Hsinchu City 30010, Taiwan
2
National Center for Computational Hydroscience and Engineering, University of Mississippi, 2301 S. Lamar Blvd, Oxford, MS 38655, USA
*
Author to whom correspondence should be addressed.
Water 2021, 13(11), 1462; https://doi.org/10.3390/w13111462
Submission received: 21 January 2021 / Revised: 14 May 2021 / Accepted: 21 May 2021 / Published: 23 May 2021
(This article belongs to the Special Issue Research on Hydraulics and River Dynamics)

Abstract

:
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 efficient 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) flow field 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.

1. 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.
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 non-uniform 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.

2. 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 depth-integrated 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].

2.1. Flow Module

The governing equations of the flow are 2D depth-integrated Reynolds equations expressed in the Cartesian coordinate system and are as follows:
u t + u u x + v u y = g η x + 1 h h τ x x x + h τ x y y τ b x ρ h + f C o r v
v t + u v x + v v y = g η y + 1 h h τ y x x + h τ y y y τ b y ρ h f C o r u
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; fCor 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:
τ i j = u i u j ¯ = ν t ( u i , j + u j , i )
Free surface elevation of the flow was calculated using the depth-integrated continuity equation:
h t + u h x + v h y = 0
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:
ν t = A x y C s κ u h
C s = 0 1 ς ( 1 ς ) d ς = 1 6
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 x y 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:
ν t = l ¯ 2 2 u x 2 + 2 v y 2 + u y + v x 2 + U ¯ z 2
l ¯ = 1 h κ z 1 z h d z = κ h 0 1 ς ( 1 ς ) d ς 0.267 κ h
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:
U ¯ z = u κ z = 1 h U z d z

2.2. 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:
( h C k ) t + ( u h C k ) x + ( v h C k ) y = x ε s h C k x + y ε s h C k y + α s ω s k ( C k C k )
Bed load:
( δ b c ¯ b k ) t + ( α b x q b k x ) x + ( α b y q b k y ) y + 1 L t q b k q b k   = 0
Bed changes:
1 p Z b k t = α ω s k C k C k + q b k q b k / L t
where εs is the eddy diffusivity of sediment in the vertical z direction; Ck 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; cbk is sediment load of the k-th size class; qb*k, q*k, qbxk and qbyk are the bed load transport capacity, the bed load transport rate, and transport rate components in x and y directions, respectively; Lt 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 Zbk 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].

3. Bed Geometry Adjustment Mechanism for a Local Scour Hole

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 D50. 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:
ϕ i = 6.08 D i + 30
ϕ i = ( 3.05 D i + 35.7 ) + ( 3.05 D i + 35.7 ) × [ 0.41 ( S 0 i 1 ) 2 + 0.21 ( S 0 i 1 ) 0.06 ]
S 0 = D 75 / D 25
where ϕ i is the angle of repose at computational node i and Di is the mean diameter of the bed material at computational node i. Moreover, S0 is the sorting coefficient, in which D75 and D25 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:
ϕ = ( ϕ i + ϕ i + 1 ) / 2
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:
Z i = Δ x tan ϕ + Z i + 1
Z i + 1 = Δ x tan ϕ + Z i
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.

4. Case Study

4.1. Initial and Boundary Conditions for the Experimental Case

An experimental case of multiple spur dikes is used for model calibration and validation. Figure 6 displays an overview of the experimental flume and spur dikes. The flume is 10 m long and 0.4 m wide with five spur dikes in the middle reach. The spur dikes are 0.1 m long and 0.02 m wide and are placed at a mutual distance of 0.2 m. The bed elevation and morphological change are obtained using the laser rangefinder system established on the flume. The terrain is measured at every 1-cm in the longitudinal and lateral directions of the experimental domain.
Table 1 lists the initial and boundary conditions of the experimental cases used for model calibration (Case No. G4a) and model validation (Case No. G5a). The inflow is steady with clear water. Figure 7 illustrates the computational domain of the experimental case. The grid numbers along the x-direction and y-direction are 337 and 63, respectively. The output profiles of the longitudinal bed elevation, e.g., A-line (Y = 0.344 m), B-line (Y = 0.252 m), C-line (Y = 0.196 m), D-line (Y = 0.146 m), and E-line (Y = 0.084 m) are presented. Three size classes ranging from 0.149 to 0.84 mm are selected as the bed material to compute the associated sediment transport.

4.2. Parameters for the Case Study

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 non-dimensional 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.
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.

4.3. 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.

5. Results

5.1. Simulation Results for the Experimental Case, Obtained Using the Original Model

Figure 10 and Figure 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.
Figure 12 and Figure 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.

5.2. 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.
Figure 15 presents a local view and comparison of the simulated results obtained using the original and adjusted models. In the original model, head scouring is found to only occur near the first spur dike. The local scour holes around the spur dikes and the erosion depth are predicted more accurately when using the proposed mechanism.

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 Lsa, Lsb, and Lsc are the distance from the head of the spur dike to the upstream, lateral, and downstream sides of the scour hole. The non-dimensional 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 Lsa/L in case G4a. This is discovered because the sediment repose angle and bed geometry adjustment mechanism are considered in the model. The simulated non-dimensional 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.
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:
M R E = 1 n i = 1 n Z i s i m u Z i m e a Z i m e a
R M S E = 1 n i = 1 n ( Z i s i m u Z i m e a ) 2
where Z i m e a and Z i s i m u 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, and Neil [23,24,25], 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.

6.2. 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.

Author Contributions

Conceptualization, K.-C.Y. and C.-T.L.; methodology, Y.-C.L. and C.-T.L.; numerical simulation, Y.-C.L. and C.-T.L.; experiment, R.-K.J.; writing—original draft preparation, C.-T.L.; writing—review and editing, K.-C.Y., C.-T.L. and Y.J. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

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

Axyturbulence viscosity coefficient
Bchannel width
Cdepth-integrated sediment concentration
Ckconcentration of the k-th size class
C*kconcentration of the k-th size class corresponding transport capacity
Csintegration constant
cbksediment load of the k-th size class
D25particle size at which 25% by weight of the bed materials is finer
D50mean diameter
D75particle size at which 75% by weight of the bed materials is finer
Dimean diameter of bed material at computational node i
Escour depth
fCorCoriolis parameter
FrFroude number
ggravitational acceleration
hwater depth
icomputational node
knumber of sediment size class
Ldistance of spur dike vertical to flume wall
l ¯ depth integrated mixing length
Lplength of dike
Lsadistance from the head of spur dike to the upstream side of scour hole
Lsbdistance from the head of spur dike to the lateral side of scour hole
Lscdistance from the head of spur dike to the downstream side of scour hole
Ltadaptation length for bed load
ntotal number of nodes along the longitudinal bed profile
pporosity of bed material
qdischarge per unit width
qb*kbed load transport capacity
q*kbed load transport rate
qbxkbed load transport rate components in x direction
qbykbed load transport rate components in y direction
Qdischarge
Sdownhill or uphill bed slope between two nodes calculated by the bed geometry
ttime in governing equations
Utotal velocity
udepth-integrated velocity component in x direction
u*shear velocity
vdepth-integrated velocity component in y direction
Vaveraged velocity
xdirection along the longitudinal bed profile
ydirection along the lateral bed profile
Zbkbed change
Z i bed elevation at computational node i modified by the angle of repose
Z i m e a measured bed changes at computational node i
Z i s i m u simulated bed changes at computational node i
α coefficient of Neill (1980) [22]
αsadaptation coefficient for suspended load
βcoefficient to convert the turbulence eddy viscosity to eddy diffusivity for sediment
Δxdistance from computational node i to i + 1
δbthickness of bed layer
εseddy diffusivity of sediment
κvon Karman’s constant
ρdensity of water
ηwater surface elevation
θangle calculated based on the mesh nodes and the elevation of bed geometry
ς relative depth of the flow
υteddy viscosity coefficient
τxxdepth integrated Reynolds stress in x direction at a plane perpendicular to x direction
τxydepth integrated Reynolds stress in y direction at a plane perpendicular to x direction
τyxdepth integrated Reynolds stress in x direction at a plane perpendicular to y direction
τyydepth integrated Reynolds stress in y direction at a plane perpendicular to y direction
τbxshear stresses on the bed and flow interface in x direction
τbyshear stresses on the bed and flow interface in y direction
ϕ i angle of repose at computational node i
ωsksettling velocity of k-th size class sediment particle

References

  1. Guney, M.S.; Aksoy, A.O.; Bombar, G. Experimental study of local scour versus time around circular bridge pier. In Proceedings of the 6th International Advanced Technologies Symposium (IATS’11), Elazığ, Turkey, 16–18 May 2011; Firat University: Elazver, Turkey, 2011. [Google Scholar]
  2. Jia, Y.; Wang Sam, S.Y. CCHE2D: Two-Dimensional Hydrodynamic and Sediment Transport Model for Unsteady Open Channel Flows Over Loose Bed; Technical Report No. NCCHE-TR-2001-1; National Center for Computational Hydroscience and Engineering, The University of Mississippi: Oxford City, MS, USA, 2001. [Google Scholar]
  3. Jia, Y.; Wang Sam, S.Y. CCHE2D Verification and Validation Tests Documentation; Technical Report No. NCCHE-TR-2001-2; National Center for Computational Hydroscience and Engineering, The University of Mississippi: Oxford City, MS, USA, 2001. [Google Scholar]
  4. Lai, Y.G.; Greimann, B.P.; Wu, K. Soft bedrock erosion modeling with a two-dimensional depth-averaged model. J. Hydraul. Eng. 2011, 137, 804–814. [Google Scholar] [CrossRef]
  5. Liao, C.T.; Yeh, K.C.; Huang, M.W. Development and application of 2-D mobile-bed model with bedrock river evolution mechanism. J. Hydro-Environ. Res. 2014, 8, 210–222. [Google Scholar] [CrossRef]
  6. Juez, C.; Murillo, J.; García-Navarro, P. Numerical assessment of bed-load discharge formulations for transient flow in 1D and 2D situations. J. Hydroinform. 2013, 15, 1234–1257. [Google Scholar] [CrossRef]
  7. Juez, C.; Murillo, J.; García-Navarro, P. A 2D weakly-coupled and efficient numerical model for transient shallow flow and movable bed. Adv. Water Resour. 2014, 55, 93–109. [Google Scholar] [CrossRef]
  8. Juez, C.; Soares-Frazao, S.; Murillo, J.; García-Navarro, P. Experimental and numerical simulation of bed load transport over steep slopes. J. Hydraul. Res. 2017, 55, 455–469. [Google Scholar] [CrossRef]
  9. Horvat, Z.; Isic, M.; Spasojevic, M. Two dimensional river flow and sediment transport model. Environ. Fluid Mech. 2015, 15, 595–625. [Google Scholar] [CrossRef]
  10. Jia, Y.; Wang, S.S.Y. Simulation of horse-shoe vortex around a bridge pier. In Proceedings of the ASCE International Conference on Water Resource Engineering, Seattle, WA, USA, 8–12 August 1999. [Google Scholar]
  11. Nagata, N.; Hosoda, T.; Nakato, T.; Muramoto, Y. Three-dimensional numerical model for flow and bed deformation around river hydraulic structures. J. Hydraul. Eng. 2005, 131, 1074–1087. [Google Scholar] [CrossRef]
  12. Burkow, M.; Griebel, M. A full three dimensional numerical simulation of the sediment transport and the scouring at a rectangular obstacle. Comput. Fluids 2016, 125, 1–10. [Google Scholar] [CrossRef]
  13. Castillo, L.G.; Carrillo, J.M. Scour, velocities and pressures evaluations produced by spillway and outlets of dam. Water 2016, 8, 68. [Google Scholar] [CrossRef] [Green Version]
  14. Jia, Y.; Altinakar, M.; Guney, M.S. Three-dimensional numerical simulations of local scouring around bridge piers. J. Hydraul. Res. 2017. [CrossRef]
  15. Yang, F.G.; Liu, X.N.; Yang, K.J.; Cao, S.Y. Study on the angle of repose of nonuniform sediment. J. Hydrodyn. 2009, 21, 685–691. [Google Scholar] [CrossRef]
  16. Meng, Z.; Wang, H.; Yang, W. Experiment on angle of repose and angle of surface friction of cohesionless sediment. J. Tianjin Univ. Sci. Technol. 2015, 48, 11. [Google Scholar]
  17. Al-Hashemi, H.M.B.; Al-Amoudi, O.S.B. A review on the angle of repose of granular materials. Powder Technol. 2018, 330, 297–417. [Google Scholar] [CrossRef]
  18. Jhong, R.K. Study on Local Scour and Sediment Exchange around Spur Dikes. Ph.D. Thesis, Department of Civil Engineering, National Chiao Tung University, Hsinchu City, Taiwan, 2017. [Google Scholar]
  19. Jia, Y.; Wang, S.S.Y. Numerical model for channel flow and morphological change studies. J. Hydraul. Eng. ASCE 1999, 125, 924–933. [Google Scholar] [CrossRef]
  20. Jia, Y.; Wang, S.Y.Y.; Xu, Y. Validation and application of a 2D model to channels with complex geometry. Int. J. Comput. Eng. Sci. 2002, 3, 57–71. [Google Scholar] [CrossRef]
  21. Wu, W.M. CCHE2D Sediment Transport Model; Technical Report No. NCCHE-TR-2001-3; National Center for Computational Hydroscience and Engineering, The University of Mississippi: Oxford City, MS, USA, 2001. [Google Scholar]
  22. Copeland, R.R. Bank Protection Techniques Using Spur Dikes; Hydraulics Laboratory, U.S. Army Engineer Waterways Experiment Station: Vicksburg, MS, USA, 1983.
  23. Liu, H.K.; Chang, F.M.; Skinner, M.M. Effect of bridge construction on scour and backwater. In CER 60 HKL 22; Colorado State University Civil and Environmental Engineering Section: Fort Collins, CO, USA, 1961. [Google Scholar]
  24. Gill, M.A. Erosion of sand beds around spur dikes. J. Hydraul. Div. 1972, 98, 1587–1602. [Google Scholar] [CrossRef]
  25. Neill, C.R. River-Bed Scour—A Review for Engineers; Technical Publication No. 23; Canadian Good Roads Association: Ottawa, ON, Canada, 1964. [Google Scholar]
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.
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.
Water 13 01462 g001
Figure 2. The shear stress in two dimensions.
Figure 2. The shear stress in two dimensions.
Water 13 01462 g002
Figure 3. Flowchart of the computational modules used in this study.
Figure 3. Flowchart of the computational modules used in this study.
Water 13 01462 g003
Figure 4. Different situations of bed elevation in the scour hole: (a) downhill bed slope (S > 0) and (b) uphill bed slope (S < 0).
Figure 4. Different situations of bed elevation in the scour hole: (a) downhill bed slope (S > 0) and (b) uphill bed slope (S < 0).
Water 13 01462 g004
Figure 5. Flowchart of the bed geometry adjustment mechanism of a local scour hole.
Figure 5. Flowchart of the bed geometry adjustment mechanism of a local scour hole.
Water 13 01462 g005
Figure 6. Overview of the multiple spur dike experimental channel: (a) top view, (b) side view, (c) top view of the morphological changes after conducting the experiment, (d) side view of the morphological changes after conducting the experiment, and (e) laser rangefinder system.
Figure 6. Overview of the multiple spur dike experimental channel: (a) top view, (b) side view, (c) top view of the morphological changes after conducting the experiment, (d) side view of the morphological changes after conducting the experiment, and (e) laser rangefinder system.
Water 13 01462 g006
Figure 7. Computational domain of the experimental case: (a) computational mesh (63 × 337 grid) and (b) initial bed elevation and layout of multiple spur dikes.
Figure 7. Computational domain of the experimental case: (a) computational mesh (63 × 337 grid) and (b) initial bed elevation and layout of multiple spur dikes.
Water 13 01462 g007
Figure 8. Bed material grain size distribution of the experimental case.
Figure 8. Bed material grain size distribution of the experimental case.
Water 13 01462 g008
Figure 9. Flowchart of model calibration and validation using the original model.
Figure 9. Flowchart of model calibration and validation using the original model.
Water 13 01462 g009
Figure 10. Overview of the simulated flow fields around multiple spur dikes in the experimental case: (a) case G4a and (b) case G5a.
Figure 10. Overview of the simulated flow fields around multiple spur dikes in the experimental case: (a) case G4a and (b) case G5a.
Water 13 01462 g010
Figure 11. Overview of the simulated velocity magnitude around multiple spur dikes in the experimental case: (a) case G4a and (b) case G5a.
Figure 11. Overview of the simulated velocity magnitude around multiple spur dikes in the experimental case: (a) case G4a and (b) case G5a.
Water 13 01462 g011
Figure 12. Overview of the simulated non-dimensional bed changes around multiple spur dikes: (a) case G4a and (b) case G5a. E is the scour depth, and h is the average water depth.
Figure 12. Overview of the simulated non-dimensional bed changes around multiple spur dikes: (a) case G4a and (b) case G5a. E is the scour depth, and h is the average water depth.
Water 13 01462 g012
Figure 13. Overview of the measured non-dimensional bed changes around multiple spur dikes: (a) case G4a and (b) case G5a. E is the scour depth, and h is the average water depth.
Figure 13. Overview of the measured non-dimensional bed changes around multiple spur dikes: (a) case G4a and (b) case G5a. E is the scour depth, and h is the average water depth.
Water 13 01462 g013
Figure 14. Overview of the simulated non-dimensional bed changes around multiple spur dikes when considering the angle of repose and adjusted bed geometry: (a) case G4a and (b) case G5a.
Figure 14. Overview of the simulated non-dimensional bed changes around multiple spur dikes when considering the angle of repose and adjusted bed geometry: (a) case G4a and (b) case G5a.
Water 13 01462 g014
Figure 15. Local view of the simulated non-dimensional bed changes around the spur dike: original and adjusted model in (a) case G4a and (b) case G5a.
Figure 15. Local view of the simulated non-dimensional bed changes around the spur dike: original and adjusted model in (a) case G4a and (b) case G5a.
Water 13 01462 g015
Figure 16. Sketch of a local scour hole around a spur dike. Modified from a study by Copeland [22].
Figure 16. Sketch of a local scour hole around a spur dike. Modified from a study by Copeland [22].
Water 13 01462 g016
Figure 17. Comparison of the measured and simulated bed profiles from the A-line to the E-line: (a) case G4a and (b) case G5a.
Figure 17. Comparison of the measured and simulated bed profiles from the A-line to the E-line: (a) case G4a and (b) case G5a.
Water 13 01462 g017
Figure 18. Local view of the comparison of the measured and simulated bed profiles at the A-line and E-line: (a) case G4a and (b) case G5a.
Figure 18. Local view of the comparison of the measured and simulated bed profiles at the A-line and E-line: (a) case G4a and (b) case G5a.
Water 13 01462 g018
Table 1. Initial and boundary conditions for the experimental cases used in model calibration and validation.
Table 1. Initial and boundary conditions for the experimental cases used in model calibration and validation.
Case No.Q (m3/s)V (m/s)h (m)Fr
G4a0.00180.1660.0270.32
G5a0.00230.170.0330.30
Note: Q = discharge, V = averaged velocity, h = averaged water depth, Fr = Froude number.
Table 2. Parameters used in model calibration.
Table 2. Parameters used in model calibration.
ParametersValues
Computational meshTotal grids (lateral × longitudinal)63 × 337
Minimum lateral mesh size (cm)0.5
Maximum lateral mesh size (cm)2.0
Minimum longitudinal mesh size (cm)0.5
Maximum longitudinal mesh size (cm)2.0
Flow moduleTime step (s)0.1
Total simulation time (s)18,000
Turbulence modelMixing length model
Turbulence viscosity coefficient2
Wall slipness coefficient1
Sediment transport moduleSediment transport formulaWu et al. [21]
Minimum mixing layer thickness (m)0.01
Adaptation length for bed load transport (m), Lt0.4
Adaptation factor for suspended load, αs0.5
Bed material size class 1 (mm)0.149
Bed material size class 2 (mm)0.42
Bed material size class 3 (mm)0.84
Table 3. Comparison of the measured and simulated non-dimensional distance of the scour hole around the spur dike.
Table 3. Comparison of the measured and simulated non-dimensional distance of the scour hole around the spur dike.
Case No.Non-Dimensional Distance
Lsa/LLsb/LLsc/L
G4aMeasured 0.981.172.21
Original Model0.410.961.90
Adjusted Model0.621.151.92
Improved (%)21.4316.240.90
G5aMeasured 1.341.543.12
Original Model0.781.022.16
Adjusted Model1.021.322.36
Improved (%)17.9119.486.41
Table 4. Comparison of the Mean Relative Errors (MREs) and Root Mean Square Errors (RMSEs) along the bed profiles from the A-line to the E-line.
Table 4. Comparison of the Mean Relative Errors (MREs) and Root Mean Square Errors (RMSEs) along the bed profiles from the A-line to the E-line.
Case No.Bed Profiles
A-LineB-LineC-LineD-LineE-Line
G4aMRE (cm)-Original0.2470.3350.0560.0920.190
MRE (cm)-Adjusted0.2260.3190.0560.0920.190
G4aRMSE (cm)-Original0.2550.4220.1230.1250.240
RMSE (cm)-Adjusted0.2440.4030.1230.1250.240
G5aMRE (cm)-Original0.4360.3470.3060.2700.292
MRE (cm)-Adjusted0.1950.2960.2840.2680.292
G5aRMSE (cm)-Original0.7760.3700.3680.3790.322
RMSE (cm)-Adjusted0.3470.3510.3420.3130.322
Table 5. Comparison of the experiment, numerical model, and empirical formula results.
Table 5. Comparison of the experiment, numerical model, and empirical formula results.
Measurement, Numerical Model, and Empirical Formula for Calculating the Head Scouring DepthScour Depth (cm)
Experimental Measurement-2.00
Original Model-0.41
Adjusted Model-2.25
Liu et al. [23] (1961) E h = 2.15 ( L p h ) 0.4 F r 0.33 7.54
Gill [24] (1972) E + h h = 8.375 ( D 50 / h ) 0.25 ( 1 L p / B ) 0.857 12.02
Neill [25] (1964) E = α 2.5 q 2 D 50 0.318 0.333 7.06–10.07
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Liao, C.-T.; Yeh, K.-C.; Lan, Y.-C.; Jhong, R.-K.; Jia, Y. Improving the 2D Numerical Simulations on Local Scour Hole around Spur Dikes. Water 2021, 13, 1462. https://doi.org/10.3390/w13111462

AMA Style

Liao C-T, Yeh K-C, Lan Y-C, Jhong R-K, Jia Y. Improving the 2D Numerical Simulations on Local Scour Hole around Spur Dikes. Water. 2021; 13(11):1462. https://doi.org/10.3390/w13111462

Chicago/Turabian Style

Liao, Chung-Ta, Keh-Chia Yeh, Yin-Chi Lan, Ren-Kai Jhong, and Yafei Jia. 2021. "Improving the 2D Numerical Simulations on Local Scour Hole around Spur Dikes" Water 13, no. 11: 1462. https://doi.org/10.3390/w13111462

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop