Understanding the complexity of the siltation process and sediment resuspension in shallow reservoirs is vital for maintaining the reservoir functionality and implementing sustainable sediment management strategies. The geometry of reservoirs plays an indispensable role in the appearance of various flow structures inside the basin and, consequently, the pattern of the morphological evolution. In this study, a three-dimensional numerical model, coupled with optimization algorithms, is used to investigate the morphological bed changes in two symmetric shallow reservoirs having hexagon and lozenge shapes. This work aims to evaluate the applicability, efficiency, and accuracy of the automatic calibration routine, which can be a suitable replacement for the time-consuming and subjective method of manual model calibration. In this regard, two sensitive parameters (i.e., roughness height and sediment active layer thickness) are assessed. The goodness-of-fit between the calculated bed levels and the measured topography from physical models are presented by different statistical metrics. From the results, it can be concluded that the automatically calibrated models are in reasonable agreement with the observations. Employing a suitable optimization algorithm, which finds the best possible combination of investigated parameters, can considerably reduce the model calibration time and user intervention.

  • The flow structure and sedimentation pattern in symmetric expansions are numerically studied.

  • Local (GML) and a global (BB–BC) optimization algorithms are used to calibrate the numerical models of two symmetric shallow reservoirs.

  • GML outperforms BB–BC considering the convergence speed (efficiency) without trapping in local minima by having the same predicted parameter values (robustness).

As general multi-purpose hydraulic structures, shallow reservoirs are used in different fields, serving as sediment and pollutant trapping tanks, retention ponds, water storage basins, and also for aquafarming (Dewals et al. 2020). Depending on several factors, such as geometrical aspects (i.e., expansion ratio: reservoir width/inlet width; and aspect ratio: reservoir length/reservoir width), hydraulic and boundary conditions, and sediment properties, the velocity field inside a reservoir emerges with various configurations, influencing the sedimentation along the basin. In turn, deposited sediments can further modify the flow pattern (Kantoush 2008; Dufresne et al. 2012). The flow behavior in shallow reservoirs, such that the horizontal dimensions are considerably larger than the water depth, can be considered a condition in which the flow passes through a narrow inlet channel into a sudden or moderate expanded basin. Thus, the flow field involves an entering jet associated with large-scale 2D horizontal coherent turbulent structures and recirculation zones, controlling the mass and momentum exchange. In addition, the flow field in symmetric expansions is specified by symmetry-breaking bifurcation under certain conditions. The importance of this behavior in engineering applications has been a subject of interest among many researchers in the field of fluid dynamics (Fani et al. 2012).

A comprehensive experimental study of laminar flow through a planar symmetric sudden expansion dates back to Durst et al. (1974), who found that the Reynolds number (Re) strictly governs the flow structure. They defined a critical Re value, beyond which the flow tends to become asymmetric. Cherdron et al. (1978) performed a similar study and confirmed that the reason for asymmetric flow formation lies in disturbances generated at the edge of the expansion, which amplify in the shear layers. Sobey & Drazin (1986) investigated the instabilities and bifurcation of flow in slight expansions by asymptotic analysis, numerical methods, and laboratory experiments. Shapira et al. (1990) carried out a linear stability analysis of symmetric flow considering different angles at the expansion section. A numerical and experimental investigation of the flow field in a symmetric expansion, indicating a symmetry-breaking (pitchfork) bifurcation point, was published by Fearn et al. (1990). They introduced a slight asymmetry in the inlet channel in their numerical model to reproduce the imperfections, which are inevitable in physical models. This small perturbation was also considered by Hawa & Rusak (2000). The effect of different expansion ratios on the occurrence of non-symmetrical flow regimes in symmetric sudden expansions was studied by Drikakis (1997). Many other researchers investigated the discussed phenomenon (i.e., the transition of the symmetric flow into the asymmetric state) in sudden expansions (e.g., Mizushima et al. 1996; Sarma et al. 2000; Quaini et al. 2016). The main focus of these studies is on determining a threshold for the critical Reynolds number as the bifurcation initiating factor and the effect of geometrical aspects on the flow structure.

In contrast to extensive investigations of laminar flow behavior downstream of sudden expansions, relatively limited studies exist about turbulent flow configurations (Escudier et al. 2002). As one of the earliest studies, Abbott & Kline (1962) experimentally investigated the turbulent flow through a symmetric plane sudden expansion. They noticed two same-length recirculation zones with a predominant central plane jet after the expansion. However, when the expansion ratio was increased, the asymmetric pattern was observed with two disproportionate recirculation regions by the deflection of the main jet toward one of the walls. They claimed that the reattachment length highly depends on the expansion ratio and is not sensitive to the Reynolds number value (turbulence intensity). Similar results were obtained by Mehta (1981), arguing the insensitivity of the asymmetric pattern for a range of high Reynolds numbers. Among other related studies on turbulent flow through symmetric expansions, the laboratory experiments of Casarsa & Giannattasio (2008) and the numerical simulations by De Zilwa et al. (2000) are worth mentioning. Most of these studies report the formation of the asymmetric flow pattern as a function of the expansion ratio.

Narrowing down the discussed phenomenon to shallow reservoirs, Kantoush (2008) and Kantoush et al. (2008a, 2008b) performed a comprehensive series of systematic investigations of the flow field with suspended sediments inside different geometries with suddenly expanded regions experimentally and numerically. They studied the effect of reservoir geometry alteration on the flow field, the silting process pattern, and the reservoir trap efficiency. The primary conclusion of these works is that an asymmetric flow pattern appears in symmetric geometries under certain conditions, affecting the process of suspended sediment transport and adjusting the spatial distribution of sediment deposits. Furthermore, the accumulation of deposited sediments due to the additional sediment supply can, in turn, alter the flow pattern. Regarding the related experimental studies in this field, Dufresne et al. (2010) classified flow patterns in shallow rectangular reservoirs into four stable categories. They found a symmetric flow pattern containing a straight jet from the inlet to the outlet with two equal-size recirculation zones at the sides for short-length reservoirs (called S0 type) and three asymmetric patterns with one, two, or three reattachment points (A1, A2, and A3) with unequal recirculation regions, depending on the length of the basin, for long-length reservoirs. In reservoirs with intermediate length, an unstable flow field was also identified (S0/A1), in which the flow randomly oscillates between the symmetric (S0) and asymmetric (A1) patterns. They investigated the effects of dimensionless length and flow depth, lateral expansion ratio, and the Froude number on the median reattachment lengths and defined a shape factor to predict the flow behavior based on the geometric aspects. Similar experimental results to those of Kantoush (2008) were obtained by Camnasio et al. (2011), who categorized the flow field inside shallow rectangular reservoirs into a channel-like (CH-L) type for reservoirs with very short width, two symmetric (S0 and S1), and two asymmetric (A1 and A2) stable types, based on the reservoir expansion and aspect ratios.

Regarding the numerical assessment of the flow behavior in rectangular shallow basins, 2D models using the shallow water equations were employed by Dewals et al. (2008) and Dufresne et al. (2011). On the other hand, the use of 3D numerical models, employing the Reynolds-averaged Navier–Stokes equations, can be found in Esmaeili et al. (2016) and Lakzian et al. (2020). Numerical models of symmetric expanding channels are complicated cases where the entering jet can randomly follow one side of the reservoir. It means that for such cases, multiple solutions for the Navier–Stokes equations can exist, where the presence of an obstacle (deposited sediments), the grid resolution, discretization schemes, and the turbulence model may cause the convergence of the equations with a different flow direction. Hence, in contrast to the extensively studied rectangular geometries, we investigated the performance of a fully 3D numerical model (SSIIM) coupled with an automatic calibration tool (PEST) to evaluate the morphological bed changes in two shallow reservoirs having hexagon and lozenge shapes. The reservoirs are selected based on available high-resolution measurements for calibration and to test the automatic model calibration for cases with complicated flow behavior. Since PEST uses a gradient-based local optimization algorithm, to confirm its ability in finding the global optimum point over the search space and avoid the local minima, the models are also calibrated by a global optimization algorithm.

Due to the fact that hydro-morphodynamic models involve a series of physically unmeasurable parameters, their accuracy and reliability highly depend on the calibration process, which means the adjustment of uncertain input parameters to minimize the misfit between simulation results and corresponding physical measurements. Generally, the calibration of such models has been manually carried out through trial-and-error parameter adjustments based on the user's understanding of the model structure and features of the environmental system until a satisfactory agreement between simulated and measured values can be achieved. By having multiple parameters for manual calibration, the typical way is to consider each parameter separately for tuning by keeping the others constant. The procedure can then be repeated for the rest of the parameters one at a time. However, in reality, the combination of the best values of every single parameter may not result in the overall best fit. Thus, having a complicated model with several uncertain input parameters, the manual calibration method becomes cost- and time-consuming, involving a high degree of subjectivity. Nevertheless, employing optimization algorithms is an ingenious approach to the model fitting process. The application of automatic model calibration in different fields of environmental studies, such as hydrologic or groundwater models, has gained popularity over the last four decades; however, there is a considerable gap in applying automatic calibration in hydro-morphodynamic studies. This work aims to evaluate the efficiency and accuracy of automatic model calibration based on mathematical optimization, which can be an innovative practice to overcome the time-consuming and subjective nature of manual model calibration.

Experimental data

The experimental data obtained by Kantoush (2008) are used to set up the numerical models and calibrate them against measured bed levels. As an outline of the experimental work, a reference rectangular shallow reservoir with 4.0 m of width, 6.0 m of length, and 0.3 m of depth was constructed at the Laboratory of Hydraulic Constructions (LCH) of the Swiss Federal Institute of Technology (EPFL). There were up- and downstream rectangular channels centrally connected to the reservoir with 0.25 m of width and 1.0 m of length. The walls were made of movable PVC plates, which could be adjusted to create different geometrical shapes. The water-sediment mixture was supplied from a mixing tank into the basin by gravity. The thickness of deposited sediments was measured by a mini echo sounder. Fine ground non-uniform walnut shells were used as sediment particles with a median diameter of d50 = 50 μm, a geometric standard deviation of particle size gradation of , a density of = 1.5 g/cm3, and an average entrance concentration of C = 3.0 g/l. Further details regarding the sediment characteristics can be found in Table 1. The discharge rate was constant during different experiments with a value of Q = 7 l/s and a corresponding water depth of h = 0.2 m, which was adjusted by a flap gate at the outlet.

Table 1

Sediment characteristics used for numerical modeling

Sediment size classes
Size (mm) 0.025 0.03 0.05 0.06 0.09 0.125 0.28 
Proportion (%) 20 10 20 10 20 10 10 
Cumulative proportion (%) 20 30 50 60 80 90 100 
Fall velocity (mm/s) 0.2 0.25 0.7 4.2 20 
Sediment size classes
Size (mm) 0.025 0.03 0.05 0.06 0.09 0.125 0.28 
Proportion (%) 20 10 20 10 20 10 10 
Cumulative proportion (%) 20 30 50 60 80 90 100 
Fall velocity (mm/s) 0.2 0.25 0.7 4.2 20 

Among the several investigated reservoir shapes in the laboratory, we used two configurations, lozenge and hexagon (Figure 1), for numerical modeling. The sediment feeding was done in four steps with different periods (3 × 1.5 hrs + 1 × 3 hrs for the lozenge-shaped reservoir, and 3 × 1.5 hrs + 1 × 4.5 hrs for the hexagonal reservoir). Furthermore, for both cases, the Froude number was small (0.1), and the Reynolds number was high enough (28,000) at the inlet to ensure subcritical, fully developed turbulent flow conditions.
Figure 1

Photos of the physical models and their sketched plan views for (a and b) hexagonal and (c and d) lozenge-shaped reservoirs.

Figure 1

Photos of the physical models and their sketched plan views for (a and b) hexagonal and (c and d) lozenge-shaped reservoirs.

Close modal

Numerical modeling

Kantoush (2008) indicated three-dimensional flow characteristics in shallow reservoirs (i.e., the presence of secondary currents and 3D stretching vortices). As long as the effect of secondary currents, generated by the streamline curvature in recirculation zones, could be excluded (such as studying the flow field without considering sediments), using 2D numerical models is reasonable. However, such models cannot directly simulate the 3D effect of secondary currents and their contribution to sedimentation. Although such flows are weaker compared to the primary flow, their role becomes significant in morphological studies with suspended sediment transport and the presence of bedforms. Accordingly, 3D numerical models can provide a more precise assessment of morphological processes in shallow reservoirs (Esmaeili et al. 2017). In this work, the fully three-dimensional numerical model SSIIM 2 (Sediment Simulation In Intakes with Multiblock option) (Olsen 2014) is used for hydro-morphodynamic simulations, which has been proven to yield reliable results in the field of reservoir sedimentation/flushing studies (Haun & Olsen 2012; Hillebrand et al. 2017; Mohammad et al. 2020). This software solves the Reynolds-averaged Navier–Stokes (RANS) equations together with the continuity equation (Equations (1) and (2)). SSIIM 2 generates an adaptive, unstructured, three-dimensional, non-orthogonal grid. The adaptive grid refinement in SSIIM 2 is based on the calculated free water surface and bed level changes resulting from the implemented wetting/drying algorithm. This algorithm calculates the number of cells that can be generated in the vertical direction after each time step as the function of the water depth so that the computational domain changes spatiotemporally and can be adjusted for the next time step. The finite-volume approach is used as the spatial discretization scheme, while an implicit scheme is employed for the temporal discretization.
(1)
(2)
where U represents the averaged velocity over the time t, x is the space coordinate, is the water density, P is the dynamic pressure, denotes the Kronecker delta (equal to 1 if , or 0 if ), and indicates the turbulent Reynolds stress term. This term is calculated according to Boussinesq's approximation with the concept of eddy viscosity and by applying the standard k–ε model for turbulence closure (Equations (3) and (4)).
(3)
(4)
Turbulent eddy-viscosity can be determined by turbulent kinetic energy k, turbulent dissipation , and a dimensionless constant equal to 0.09. Regarding the other terms of the RANS equations, the pressure term is handled by the semi-implicit method for pressure-linked equations, and the convective term is modeled by the second-order upwind scheme. The free water surface is calculated based on the computed pressure gradient by using the Bernoulli equation. The flow discharge, turbulence parameters, and inflow sediment concentration are prescribed at the inlet as a Dirichlet boundary condition. At the outlet, the Neumann type zero gradient boundary condition is used for all variables. Wall laws for rough boundaries introduced by Schlichting (1979) are used in the simulations (Equation (5)).
(5)
where U is the flow velocity, is the shear velocity, is the von Kármán constant equal to 0.41, is the distance between the wall and the center of the closest cell, and is the equivalent sediment roughness height.
Concerning the grid resolution, both models (lozenge and hexagon) are composed of 120 × 80 × 5 cells in the streamwise, lateral, and vertical directions, respectively. The suspended sediment transport computation is carried out by solving the advection-diffusion equation (Equation (6)), considering van Rijn's (1984a) formula for the equilibrium near-bed concentration as a boundary condition (Equation (7)).
(6)
where is the sediment concentration of size, w is the sediment fall velocity, and is the turbulent diffusion coefficient, where is the Schmidt number. The value for the Schmidt number is taken to be equal to unity by assuming no deviation between the eddy viscosity and the turbulent diffusivity.
Suspended sediment transport calculation requires the specification of the near-bed sediment concentration, considering a reference level above the bed, where sediment resuspension occurs (reference concentration). The empirical formula of van Rijn is used to calculate the equilibrium suspended sediment concentration in the cells close to the bed as the bottom boundary condition for solving the advection-diffusion equation. The near-bed volumetric sediment concentration for the th fraction is calculated as follows:
(7)
where is the diameter of the th fraction, a represents the roughness height as the reference level, is the bed shear stress, is the critical shear stress for the movement of the th fraction according to the Shields curve, and are the density of sediment and water, respectively, g is the gravity acceleration, and v represents the kinematic viscosity. The sediment resuspension from the bed can be calculated by converting the sediment concentration into the entrainment rate.
The bedload transport is also calculated by the empirical formula of van Rijn (1984b) (Equation (8)), which represents the transport rate of the ith fraction per unit width.
(8)

The sediment continuity equation (Exner equation) is used to calculate the bed elevation changes with respect to the continuity defect in cells close to the bed. The difference between the sediment inflow and outflow in a cell is multiplied by the time step and divided by the horizontal plane area of the cell. The bed is then raised or lowered after each time step, and the grid is regenerated.

Regarding the velocity distribution at the inlet, a minor linear perturbation is introduced to the velocity profile, described by Dewals et al. (2008), to provide an initiating factor and a condition for the genesis of the asymmetric flow in symmetric numerical models. This small perturbation can be considered as an initial boundary condition that is unavoidable in experimental set-ups, and can be further damped and become a straight jet regarding the reservoirs with stable-symmetric flow structure. In other words, a perfect symmetric numerical model needs a symmetry-breaking condition to reproduce an asymmetric flow field.

Model calibration

Model calibration is regarded as the initial stage of appraising the performance of a computational model that represents the complex behavior of a real-world system. This inverse process depends on how sensible the uncertain affecting parameters are adjusted according to the misfit minimization between calculated and corresponding measured values. The traditional manual trial-and-error method of model calibration has nowadays been overshadowed by the concept of mathematical optimization, which paves the way for automatic model calibration owing to its interesting features such as objective-based judgment of goodness-of-fit rather than relying on the user subjectivity and being less time-consuming (Evangelista et al. 2017). The principal elements of an automatic calibration routine are: an objective function to evaluate the model performance, an optimization algorithm to explore the parameter space by repetitive adjustment of the uncertain parameters, and a termination criterion to stop the search when the convergence of the objective function or the maximum allowable number of iterations is satisfied (Vidal et al. 2007).

In this study, the model-independent nonlinear Parameter ESTimation and predictive analysis package PEST (Doherty 2016), which employs the gradient-based Gauss–Marquardt–Levenberg (GML) local optimization algorithm, is used to calibrate the numerical models. This tool has demonstrated promising results for sensitivity analysis and automatic calibration of numerical models in different environmental studies (Shoarinezhad et al. 2020a, 2020b). Further, a global optimization algorithm, Big Bang Big Crunch (BB–BC) (Erol & Eksin 2006), is applied to validate PEST performance in finding the global optimum over the search space.

The numerical models are calibrated in accordance with the measured bed levels in different longitudinal and cross-sections, taking 8,600 and 16,500 points into account for the lozenge- and hexagon-shaped reservoirs, respectively. The residual sum of squares between measured and calculated bed levels, , is used as the optimization objective function for the automatic calibration of the models. The measured points were located in every 10 cm of cross-sections (61 cross-sections) with a lateral distance of about 1 cm from each other. Among the four phases of sediment feeding in the experiments, the measured bed levels from the first 1.5 hrs (considered a warm-up period) are used as the initial bed levels to simulate the following three phases. This is performed to provide an initial condition for the roughness calibration considering the surface friction development from the hydraulically smooth bed surface.

Gauss–Marquardt–Levenberg algorithm of PEST

PEST iteratively adjusts the uncertain parameters within a pre-defined space with the aid of the GML search algorithm until reaching the minimum deviation between the measured and calculated values according to the residual sum of squares as the objective function. Using the GML algorithm, as a combination of the gradient descent method and Gauss–Newton algorithm, PEST runs the model and linearizes the relationship between input parameters and model outputs by Taylor's expansion of the actual parameter set. The number of adjustable parameters, subject to calibration, determines the number of model runs in a single PEST iteration. During each iteration, a Jacobian matrix of partial derivatives of model outputs is generated, followed by an upgrade vector (Equation (9)) to alter the parameters for the next iteration up to reaching either the minimum of the objective function or termination criteria.
(9)
where u is the parameter upgrade vector, represents the transpose of the Jacobian matrix J, Q is the diagonal weight matrix, is the Marquardt lambda acting as a damping factor, I is an identity matrix, and r is the vector of residuals.

Big Bang Big Crunch algorithm

BB–BC is a metaheuristic optimization algorithm inspired by a theory for the evolution of the universe. The initial population is uniformly generated by spreading random candidate solutions over the entire search space in the Big Bang phase. The main concept of this phase is based on the energy dissipation in nature which creates chaos and randomness in the population. Then, the fitness value for each candidate solution is calculated (i.e., the ‘mass’ of each particle). During the Big Crunch phase, the randomly distributed population shrinks to a single point , which is called the ‘center of mass’, and calculated according to Equation (10).
(10)

Here, is the position vector of the th candidate in an n-dimensional search space, N is the population size, and denotes the fitness value of the th candidate. In the next Big Bang, new individuals are mainly generated around the former-calculated center of mass according to a normal distribution, where the standard deviation of the normal distribution decreases as the optimization proceeds. This step is followed by a contraction according to the recalculation of the new center of mass. In order to ensure the global convergence of the method, the algorithm always generates a number of new solutions far from the center of mass with a diminishing probability as iterations go forward. Eventually, with regard to a defined termination criterion, this successive two-phase scheme (explosion-contraction) converges to the optimum point (Kaveh & Bakhshpoori 2019).

Parameters selection

A series of sensitivity analyses are performed to find out the significance of uncertain input parameters (e.g., roughness coefficient, the porosity of compacted bed sediments, the coefficient for the bedform smoothing algorithm, the thickness of the upper active sediment layer, and the angle of repose for sediments) on the system behavior as the initial stage prior to the model calibration. SENSAN, a model-independent local SENSitivity ANalyzer as one of the PEST utilities, is used for the sensitivity analysis. SENSAN conducts several runs according to the pre-defined sets of parameter values and records the model output sensitivities to parameter changes. Among the assessed parameters, the two most sensitive ones are selected for auto-calibration as follows:

  • Nikuradse equivalent roughness height (), which is generally considered to be proportionately related to the representative grain size, such that . There are several values for and various sediment sizes as the representative grain size in the literature (e.g., see Dey (2014)). Here, the range of this parameter is set to be d50 ≤ 10d90.

  • Active layer thickness (ALT), as the superficial exchange layer depth, engaged in the entrainment and deposition of sediment particles. The sorting mechanism occurs in the active layer, where the sediment continuity equation is computed separately for each size fraction inside a cell during each time step. Depending on the transport regime, ALT scales with the representative grain size of the sediment mixture or is defined as a fraction of bedforms height. The range for this parameter is selected to be d50ALT ≤ 5dmax (Malcherek 2007).

Since the GML algorithm is a gradient-based approach, it might likely find the local optimum point on the objective function surface rather than the desired global one. Therefore, it is worthwhile to reassess the calibration procedure by using different random initial values within the parameter space in the case of using local optimization algorithms. Here, two different pre-defined starting values are considered for the investigated parameters:

  • GML#1: = d90 = 0.013 cm ALT1 = dmax = 0.028 cm

  • GML#2: = 3d90 = 0.039 cm ALT2 = 3dmax = 0.084 cm

The geometry of a reservoir plays a vital role in the flow field development and, consequently, the sedimentation pattern. Regarding this fact, different relations in the literature can be found, the so-called shape factors, which correlate geometrical aspects with the appearance of various flow field categories (e.g., Kantoush (2008); Dufresne et al. (2011)). Above a critical value for the shape factor, the symmetric flow inside a symmetric reservoir evolves toward the asymmetric pattern due to the Coanda effect, where the flow deviates to one side of the reservoir according to a tiny imperfection in the physical symmetric configuration.

Figure 2 shows the calculated final surface velocity magnitude for hexagonal (a, b, c) and lozenge-shaped (d, e, f) reservoirs in different periods (i.e., phases 1, 2, and 3, up to bottom). Concerning the hexagonal configuration, the main jet enters the reservoir and keeps its straight path up to the outlet with recirculation zones on both sides. It is evident that the deposited sediments during different periods do not affect the stability of the velocity field and do not interrupt the flow symmetricity (Figure 2(a)–2(c)). Nevertheless, there is a significant deviation in the flow trajectory in the lozenge-shaped reservoir, where the entering jet reattaches to one side, resulting in a single large recirculation zone. Here, the flow field is much more complex, and the unstable nature of the flow gives rise to a shift in the velocity field from one side to the opposite side of the reservoir due to the accumulated sediments. Although there is a gradual shift of the flow pattern from clockwise in phase 1 to counterclockwise in phase 2, the velocity field keeps its counterclockwise route during phase 3 (Figure 2(d)–2(f)).
Figure 2

Contour maps of calculated velocity fields for the two reservoirs in different periods: 1.5 hrs of (a and d) phase 1 and (b and e) phase 2 for both reservoirs; (c) 4.5 hrs and (f) 3 hrs of phase 3 for the hexagonal and lozenge-shaped reservoirs, respectively.

Figure 2

Contour maps of calculated velocity fields for the two reservoirs in different periods: 1.5 hrs of (a and d) phase 1 and (b and e) phase 2 for both reservoirs; (c) 4.5 hrs and (f) 3 hrs of phase 3 for the hexagonal and lozenge-shaped reservoirs, respectively.

Close modal
Figure 3 depicts the middle cross-sectional profile (x = 3 m) of (a) hexagonal and (b) lozenge-shaped reservoirs at the end of the final phase. The color contour maps show the calculated streamwise velocity (Ux), where the vectors are the resultant of the lateral (Uy) and vertical (Uz) velocity components.
Figure 3

Central cross-sectional view of the final velocity field for (a) hexagonal and (b) lozenge-shaped reservoirs.

Figure 3

Central cross-sectional view of the final velocity field for (a) hexagonal and (b) lozenge-shaped reservoirs.

Close modal

The auto-calibration process is based on the pairwise comparison of the calculated and measured bed levels. Table 2 shows the calibration results for both reservoirs using GML (with two different starting values for the selected parameters) and BB–BC (with randomly sampled initial values) algorithms. It can be seen that the calibrated values for the investigated parameters (calibrated and ALT calibrated) are very similar for GML#1 and GML#2 regarding both reservoirs, indicating the repeatability and robustness of the auto-calibration procedure based on the GML algorithm. Moreover, these values are almost the same as the results of the BB–BC algorithm. It means the gradient-based GML algorithm is not affected by local minima, and PEST can reliably calibrate the numerical models in a global manner. Furthermore, among innumerable parameter combinations, PEST calibrates the models with 30–40 runs depending on the starting values, which shows the method's efficiency compared to the BB–BC algorithm.

Table 2

Initial values of the investigated parameters and calibration results (final calibrated values and the number of model runs)

Reservoir shapeAlgorithmInitial values
Calibration results
ksinitialALT initialkscalibratedALT calibratedModel runs
Lozenge GML#1 0.013 0.028 0.0212 0.0308 32 
GML#2 0.039 0.084 0.0214 0.0311 38 
BB–BC Randomly sampled 0.0212 0.0309 471 
Hexagon GML#1 0.013 0.028 0.0238 0.0382 31 
GML#2 0.039 0.084 0.0235 0.0377 35 
BB–BC Randomly sampled 0.0237 0.0380 438 
Reservoir shapeAlgorithmInitial values
Calibration results
ksinitialALT initialkscalibratedALT calibratedModel runs
Lozenge GML#1 0.013 0.028 0.0212 0.0308 32 
GML#2 0.039 0.084 0.0214 0.0311 38 
BB–BC Randomly sampled 0.0212 0.0309 471 
Hexagon GML#1 0.013 0.028 0.0238 0.0382 31 
GML#2 0.039 0.084 0.0235 0.0377 35 
BB–BC Randomly sampled 0.0237 0.0380 438 

Note: , roughness height at the bed; ALT, active layer thickness. Units are in centimeters.

Considering d90 = 0.013 cm as the representative grain size, the calibrated values can be rewritten as: ks ≈ 1.64d90 and ALT ≈ 2.4d90 for the lozenge-shaped reservoir; and ks ≈ 1.84d90 and ALT ≈ 2.94d90 regarding the hexagonal reservoir.

Figures 4 and 5 illustrate plan views of the measured and simulated bed topography at the end of each period (i.e., phases 1, 2, and 3, up to bottom) for hexagonal and lozenge-shaped reservoirs, respectively. In the hexagonal reservoir (Figure 4), the flow field has a stable behavior during different periods with a continuous straight jet. Hence, the main part of sediment particles is settled along this mid-longitudinal section. In addition to the main central flow, two large recirculation zones are in charge of the lateral depositions. In phase 1 (Figure 4(a) and 4(d)), the magnitude of these lateral eddies in the experimental set-up and the numerical model are almost identical. However, in phase 2 (Figure 4(b) and 4(e)), there is a stronger vortex at the right side of the flow direction in the laboratory model. Here, more sediments are transported and settled at the right side of the reservoir in the experimental set-up compared to the uniform pattern of the numerical model. According to the final bed levels for the hexagonal reservoir (Figure 4(c) and 4(f)), a scour hole occurs in the immediate upstream of the experiment, which could not be simulated in the numerical model. The overall patterns of final bed levels are similar, with the maximum deposition along the first half of the central-longitudinal section.
Figure 4

Comparison of bed level changes in the experiment (left) with the simulation results (right) in (a and d) phase 1, (b and e) phase 2, and (c and f) phase 3 for the hexagonal reservoir.

Figure 4

Comparison of bed level changes in the experiment (left) with the simulation results (right) in (a and d) phase 1, (b and e) phase 2, and (c and f) phase 3 for the hexagonal reservoir.

Close modal
Figure 5

Comparison of bed level changes in the experiment (left) with the simulation results (right) in (a and d) phase 1, (b and e) phase 2, and (c and f) phase 3 for the lozenge-shaped reservoir.

Figure 5

Comparison of bed level changes in the experiment (left) with the simulation results (right) in (a and d) phase 1, (b and e) phase 2, and (c and f) phase 3 for the lozenge-shaped reservoir.

Close modal

Despite the stable, straight, and symmetric nature of the flow in the hexagonal reservoir, which results in a central-longitudinal pattern of sediment depositions, the flow structure in the lozenge-shaped reservoir is unstable and fluctuating with main sediment deposits at the sides of the reservoir (Figure 5). During the first phase, the flow direction is clockwise, and sediments are mostly settled at the left part, resulting in an asymmetric bed topography (Figure 5(a) and 5(d)). The shift in the flow path during the second phase, due to the deposited sediments that can modify the unstable and sensitive flow structure in the lozenge-shaped reservoir, changes the bed topography to a semi-symmetric pattern (Figure 5(b) and 5(e)). Here, the numerical model underestimates the bed levels at the second half of the right side of the reservoir. Considering the final bed levels, sediments are almost symmetrically deposited along the left and right parts of the basin in the experiment, whereas the maximum deposition is along the first half of the right side in the numerical model (Figure 5(c) and 5(f)). Furthermore, in all phases, an erosion area can be seen just in front of the inlet in the experiment, which is deeper and larger in size compared to the results of the numerical model.

In order to quantitatively evaluate the overall performance of the automatically calibrated numerical models, various statistical metrics are used to compare the final measured and calculated bed levels, as shown in Table 3. The results of uncalibrated models using the initial parameter values of GML#1 are also presented to see the deviation between the results of the calibrated models and our initial guess. Mean Bias Error (MBE) is applied as a bias indicator, describing the degree of underprediction (negative values) or overprediction (positive values) of the model. Root Mean Squared Error (RMSE) and Mean Absolute Error (MAE) both reflect the average magnitude of the error. While MAE indicates a linear behavior of individual errors, RMSE gives more importance to big errors by giving higher weights to them. Pearson correlation coefficient (R) shows the linear correlation between estimated and measured bed levels. As a multi-component goodness-of-fit, Kling-Gupta efficiency (KG) is used, which combines Pearson correlation coefficient, bias, and variability within a single objective function.

Table 3

Statistical performance of the automatic model calibration

Reservoir shapeGoodness-of-fit
MBE (cm)RMSE (cm)MAE (cm)R (-)KG (-)
Lozenge Calibrated −0.024 0.054 0.044 0.80 0.77 
Initial guess −0.061 0.288 0.212 0.56 0.52 
Hexagon Calibrated −0.013 0.055 0.039 0.86 0.85 
Initial guess −0.039 0.252 0.208 0.70 0.68 
Reservoir shapeGoodness-of-fit
MBE (cm)RMSE (cm)MAE (cm)R (-)KG (-)
Lozenge Calibrated −0.024 0.054 0.044 0.80 0.77 
Initial guess −0.061 0.288 0.212 0.56 0.52 
Hexagon Calibrated −0.013 0.055 0.039 0.86 0.85 
Initial guess −0.039 0.252 0.208 0.70 0.68 

Note: MBE, Mean Bias Error; RMSE, Root Mean Squared Error; MAE, Mean Absolute Error; R, Pearson correlation coefficient; KG, Kling-Gupta efficiency.

According to the negative values of MBE, the numerical models underestimate the final bed levels in both cases. Although having a similar calibrated RMSE, the lower MAE value for the hexagonal model represents its more precise calibration (the best value is 0). This can further be confirmed by higher values of correlation and Kling-Gupta efficiency metrics for the hexagonal model (the best value is 1). Nevertheless, the statistical performance of the calibrated lozenge-shaped model also agrees reasonably with the measured data.

In this study, we applied an automatic calibration tool (PEST), which uses the gradient-based Gauss–Marquardt–Levenberg optimization algorithm, to calibrate the 3D morphodynamic numerical models of two reservoirs (hexagonal and lozenge-shaped configurations) against the experimental data. The two most affecting input parameters (roughness height and active layer thickness) in the numerical models are selected for calibration through a sensitivity analysis. In order to verify the ability of the gradient-based optimization algorithm to find the global optimum values of the parameters over the search space rather than sticking to a local optimum point (evaluating the robustness and convergence of the algorithm), in addition to using the global optimization algorithm Big Bang Big Crunch, we tested the gradient-based method with two different pre-defined initial values with a reasonable range based on the literature. Achieving an almost identical set of calibrated values confirmed the accuracy and reliability of the optimization procedure. The parameters are calibrated with 30–40 model runs by PEST, which shows its efficiency and superiority over the trial-and-error manual calibration, suggesting its potential use for hydro-morphodynamic models calibration.

So far, most related research works have been on rectangular reservoirs with a 90° expansion angle. In this work, the effect of lower expansion angles on the flow field development and the sedimentation pattern is assessed numerically. Keeping the maximum inner dimensions constant, the lower expansion angle gives a higher chance of asymmetric flow appearance. Regarding the hexagonal reservoir, which has a stable flow structure, the main part of sediment particles is settled along the mid-longitudinal section. However, the unstable flow pattern inside the lozenge-shaped reservoir causes sediments to be distributed at the sides of the basin. The calculated bed levels are compared with the measured topography of the physical models in different time steps. Considering the primary features and specifications of both reservoirs regarding the flow field, jet direction, recirculation zones, and bed topography in the experimental work, the calibrated numerical models can reasonably reproduce similar patterns. The achieved results regarding the flow fields and erosion/deposition patterns are used as a criterion to judge the ability of the Gauss–Marquardt–Levenberg search algorithm to calibrate the numerical models (optimization performance). A difference between the bed levels in the physical set-up and simulation results is related to the scour hole in front of the inlet, likely due to boundary effects, which cannot be predicted by the numerical model.

The quality and performance of the calibrated numerical models are also investigated by different statistical metrics, comparing the predicted and measured bed levels. The mean bias error shows an underestimation of bed levels in both reservoirs. Nevertheless, according to the low values of the root mean squared error and the mean absolute error; and high values of the correlation coefficient and Kling-Gupta efficiency, the overall performance of the automatic calibration procedure is reasonable.

The manual trial-and-error model calibration approach with just one parameter subject to alteration may be reasonable and sufficient in most cases. However, if there are several uncertain parameters, their innumerable combinations cause the manual calibration method to become much more complex, time- and cost-consuming, and impractical. What is more, since the issue of subjectivity is involved in manual model calibration, there is no guarantee that the best possible combination of parameters can be achieved. The overall outcome of this study is that using suitable optimization algorithms for hydro-morphodynamic models calibration, which is not a common practice among researchers in this field, can considerably reduce the calibration time and user intervention/subjectivity and concurrently increase the precision of the process.

The first author received a Ph.D. scholarship from the German Federal Ministry of Education and Research (BMBF) through the Sustainable Water Management-NaWaM program of the German Academic Exchange Service (DAAD). The last author is indebted to the Baden-Württemberg Stiftung for the financial support from the Elite program for Postdocs. The authors also acknowledge Prof. Nils Reidar Bøe Olsen, the developer of SSIIM software, for his valuable suggestion regarding numerical modeling.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.

Abbott
D. E.
&
Kline
S. J.
1962
Experimental investigation of subsonic turbulent flow over single and double backward facing steps
.
Journal of Basic Engineering
84
,
317
325
.
https://doi.org/10.1115/1.3657313
.
Camnasio
E.
,
Orsi
E.
&
Schleiss
A. J.
2011
Experimental study of velocity fields in rectangular shallow reservoirs
.
Journal of Hydraulic Research
49
,
352
358
.
https://doi.org/10.1080/00221686.2011.574387
.
Casarsa
L.
&
Giannattasio
P.
2008
Three-dimensional features of the turbulent flow through a planar sudden expansion
.
Physics of Fluids
20
,
015103
.
https://doi.org/10.1063/1.2832780
.
Cherdron
W.
,
Durst
F.
&
Whitelaw
J. H.
1978
Asymmetric flows and instabilities in symmetric ducts with sudden expansions
.
Journal of Fluid Mechanics
84
,
13
.
https://doi.org/10.1017/S0022112078000026
.
Dewals
B. J.
,
Kantoush
S. A.
,
Erpicum
S.
,
Pirotton
M.
&
Schleiss
A. J.
2008
Experimental and numerical analysis of flow instabilities in rectangular shallow basins
.
Environmental Fluid Mechanics
8
,
31
54
.
https://doi.org/10.1007/s10652-008-9053-z
.
Dewals
B.
,
Archambeau
P.
,
Bruwier
M.
,
Erpicum
S.
,
Pirotton
M.
,
Adam
T.
,
Delhez
E.
&
Deleersnijder
E.
2020
Age of water particles as a diagnosis of steady-state flows in shallow rectangular reservoirs
.
Water
12
,
2819
.
https://doi.org/10.3390/w12102819
.
Dey
S.
2014
Fluvial Hydrodynamics: Hydrodynamic and Sediment Transport Phenomena
.
Springer
,
Berlin, Heidelberg
.
https://doi.org/10.1007/978-3-642-19062-9.
De Zilwa
S. R. N.
,
Khezzar
L.
&
Whitelaw
J. H.
2000
Flows through plane sudden-expansions
.
International Journal for Numerical Methods in Fluids
32
,
313
329
.
https://doi.org/10.1002/(SICI)1097-0363(20000215)32:3 < 313::AID-FLD940 > 3.0.CO;2-B
.
Doherty
J.
2016
PEST Model-Independent Parameter Estimation User Manual Part I, 6th ed
.
Watermark Numerical Computing
,
Brisbane
,
Australia
.
Drikakis
D.
1997
Bifurcation phenomena in incompressible sudden expansion flows
.
Physics of Fluids
9
,
76
87
.
https://doi.org/10.1063/1.869174
.
Dufresne
M.
,
Dewals
B. J.
,
Erpicum
S.
,
Archambeau
P.
&
Pirotton
M.
2010
Classification of flow patterns in rectangular shallow reservoirs
.
Journal of Hydraulic Research
48
,
197
204
.
https://doi.org/10.1080/00221681003704236
.
Dufresne
M.
,
Dewals
B. J.
,
Erpicum
S.
,
Archambeau
P.
&
Pirotton
M.
2011
Numerical investigation of flow patterns in rectangular shallow reservoirs
.
Engineering Applications of Computational Fluid Mechanics
5
,
247
258
.
https://doi.org/10.1080/19942060.2011.11015368
.
Dufresne
M.
,
Dewals
B.
,
Erpicum
S.
,
Archambeau
P.
&
Pirotton
M.
2012
Flow patterns and sediment deposition in rectangular shallow reservoirs
.
Water and Environment Journal
26
,
504
510
.
https://doi.org/10.1111/j.1747-6593.2012.00310.x
.
Durst
F.
,
Melling
A.
&
Whitelaw
J. H.
1974
Low reynolds number flow over a plane symmetric sudden expansion
.
Journal of Fluid Mechanics
64
,
111
128
.
https://doi.org/10.1017/S0022112074002035
.
Erol
O. K.
&
Eksin
I.
2006
A new optimization method: Big Bang–Big Crunch
.
Advances in Engineering Software
37
,
106
111
.
https://doi.org/10.1016/j.advengsoft.2005.04.005
.
Escudier
M. P.
,
Oliveira
P. J.
&
Poole
R. J.
2002
Turbulent flow through a plane sudden expansion of modest aspect ratio
.
Physics of Fluids
14
,
3641
3654
.
https://doi.org/10.1063/1.1504711
.
Esmaeili
T.
,
Sumi
T.
,
Kantoush
S. A.
,
Haun
S.
&
Rüther
N.
2016
Three-dimensional numerical modelling of flow field in shallow reservoirs
.
Proceedings of the Institution of Civil Engineers-Water Management
169
,
229
244
.
https://doi.org/10.1680/wama.15.00011
.
Esmaeili
T.
,
Sumi
T.
,
Kantoush
S.
,
Kubota
Y.
,
Haun
S.
&
Rüther
N.
2017
Three-dimensional numerical study of free-flow sediment flushing to increase the flushing efficiency: a case-study reservoir in Japan
.
Water
9
,
900
.
https://doi.org/10.3390/w9110900
.
Evangelista
S.
,
Giovinco
G.
&
Kocaman
S.
2017
A multi-parameter calibration method for the numerical simulation of morphodynamic problems
.
Journal of Hydrology and Hydromechanics
65
,
175
182
.
https://doi.org/10.1515/johh-2017-0014
.
Fani
A.
,
Camarri
S.
&
Salvetti
M. V.
2012
Stability analysis and control of the flow in a symmetric channel with a sudden expansion
.
Physics of Fluids
24
.
https://doi.org/10.1063/1.4745190.
Fearn
R. M.
,
Mullin
T.
&
Cliffe
K. A.
1990
Nonlinear flow phenomena in a symmetric sudden expansion
.
Journal of Fluid Mechanics
211
,
595
608
.
https://doi.org/10.1017/S0022112090001707
.
Haun
S.
&
Olsen
N. R. B.
2012
Three-dimensional numerical modelling of reservoir flushing in a prototype scale
.
International Journal of River Basin Management
10
,
341
349
.
https://doi.org/10.1080/15715124.2012.736388
.
Hawa
T.
&
Rusak
Z.
2000
Viscous flow in a slightly asymmetric channel with a sudden expansion
.
Physics of Fluids
12
,
2257
2267
.
https://doi.org/10.1063/1.1287610
.
Hillebrand
G.
,
Klassen
I.
&
Olsen
N. R. B.
2017
3D CFD modelling of velocities and sediment transport in the Iffezheim hydropower reservoir
.
Hydrology Research
48
,
147
159
.
https://doi.org/10.2166/nh.2016.197
.
Kantoush
S. A.
2008
Experimental Study on the Influence of the Geometry of Shallow Reservoirs on Flow Patterns and Sedimentation by Suspended Sediments
.
PhD Thesis
,
Écolé Polytechnique Fédérale de Lausanne (EPFL)
,
Lausanne
,
Switzerland
.
Kantoush
S. A.
,
Bollaert
E.
&
Schleiss
A. J.
2008a
Experimental and numerical modelling of sedimentation in a rectangular shallow basin
.
International Journal of Sediment Research
23
,
212
232
.
https://doi.org/10.1016/S1001-6279(08)60020-7
.
Kantoush
S. A.
,
De Cesare
G.
,
Boillat
J. L.
&
Schleiss
A. J.
2008b
Flow field investigation in a rectangular shallow reservoir using UVP, LSPIV and numerical modelling
.
Flow Measurement and Instrumentation
19
,
139
144
.
https://doi.org/10.1016/j.flowmeasinst.2007.09.005
.
Kaveh
A.
&
Bakhshpoori
T.
2019
Big Bang-Big Crunch Algorithm
. In:
Metaheuristics: Outlines, MATLAB Codes and Examples
.
Springer International Publishing
,
Cham
, pp.
31
40
.
https://doi.org/10.1007/978-3-030-04067-3_4.
Lakzian
E.
,
Saghi
H.
&
Kooshki
O.
2020
Numerical simulation of sediment deposition and trapping efficiency estimation in settling basins, considering secondary flows
.
International Journal of Sediment Research
35
,
347
354
.
https://doi.org/10.1016/j.ijsrc.2020.02.001
.
Malcherek
A.
2007
Sedimenttransport und Morphodynamik
.
Scriptum Institut für Wasserwesen, Bundeswehr University Munich
,
Munich
,
Germany
.
Mehta
P. R.
1981
Separated flow through large sudden expansions
.
Journal of the Hydraulics Division
107
,
451
460
.
https://doi.org/10.1061/JYCEAJ.0005648
.
Mizushima
J.
,
Okamoto
H.
&
Yamaguchi
H.
1996
Stability of flow in a channel with a suddenly expanded part
.
Physics of Fluids
8
,
2933
2942
.
https://doi.org/10.1063/1.869072
.
Mohammad
M. E.
,
Al-Ansari
N.
,
Knutsson
S.
&
Laue
J.
2020
A computational fluid dynamics simulation model of sediment deposition in a storage reservoir subject to water withdrawal
.
Water
12
,
959
.
https://doi.org/10.3390/w12040959
.
Olsen
N. R. B.
2014
A Three Dimensional Numerical Model for Simulation of Sediment Movement in Water Intakes with Multiblock Option
.
Department of Hydraulic and Environmental Engineering, The Norwegian University of Science and Technology
,
Trondheim
,
Norway
.
Quaini
A.
,
Glowinski
R.
&
Čanić
S.
2016
Symmetry breaking and preliminary results about a Hopf bifurcation for incompressible viscous flow in an expansion channel
.
International Journal of Computational Fluid Dynamics
30
,
7
19
.
https://doi.org/10.1080/10618562.2016.1144877
.
Sarma
A. S. R.
,
Sundararajan
T.
&
Ramjee
V.
2000
Numerical simulation of confined laminar jet flows
.
International Journal for Numerical Methods in Fluids
33
,
609
626
.
https://doi.org/10.1002/1097-0363(20000715)33:5 < 609::AID-FLD745 > 3.0.CO;2-%23
.
Schlichting
H.
1979
Boundary Layer Theory
.
McGraw-Hill
,
New York
,
USA
.
Shapira
M.
,
Degani
D.
&
Weihs
D.
1990
Stability and existence of multiple solutions for viscous flow in suddenly enlarged channels
.
Computers & Fluids
18
,
239
258
.
https://doi.org/10.1016/0045-7930(90)90009-M
.
Shoarinezhad
V.
,
Wieprecht
S.
,
Haun
S.
,
2020a
Automatic calibration of a 3D morphodynamic numerical model for simulating bed changes in a 180° channel bend
. In:
Recent Trends in Environmental Hydraulics, GeoPlanet: Earth and Planetary Sciences
(
Kalinowska
M. B.
,
Mrokowska
M. M.
&
Rowiński
P. M.
, eds).
Springer International Publishing
,
Cham
, pp.
253
262
.
https://doi.org/10.1007/978-3-030-37105-0_22
Shoarinezhad
V.
,
Wieprecht
S.
&
Haun
S.
2020b
Comparison of local and global optimization methods for calibration of a 3D morphodynamic model of a curved channel
.
Water
12
,
1333
.
https://doi.org/10.3390/w12051333
.
Sobey
I. J.
&
Drazin
P. G.
1986
Bifurcations of two-dimensional channel flows
.
Journal of Fluid Mechanics
171
,
263
287
.
https://doi.org/10.1017/S0022112086001441
.
van Rijn
L. C.
1984a
Sediment transport, part II: suspended load transport
.
Journal of Hydraulic Engineering
110
,
1613
1641
.
https://doi.org/10.1061/(ASCE)0733-9429(1984)110:11(1613)
.
van Rijn
L. C.
1984b
Sediment transport, part I: bed load transport
.
Journal of Hydraulic Engineering
110
,
1431
1456
.
https://doi.org/10.1061/(ASCE)0733-9429(1984)110:10(1431)
.
Vidal
J.-P.
,
Moisan
S.
,
Faure
J.-B.
&
Dartus
D.
2007
River model calibration, from guidelines to operational support tools
.
Environmental Modelling & Software
22
,
1628
1640
.
https://doi.org/10.1016/j.envsoft.2006.12.003
.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY 4.0), which permits copying, adaptation and redistribution, provided the original work is properly cited (http://creativecommons.org/licenses/by/4.0/).