Comparison of Local and Global Optimization Methods for Calibration of a 3D Morphodynamic Model of a Curved Channel

Abstract

In curved channels, the flow characteristics, sediment transport mechanisms, and bed evolution are more complex than in straight channels, owing to the interaction between the centrifugal force and the pressure gradient, which results in the formation of secondary currents. Therefore, using an appropriate numerical model that considers this fully three-dimensional effect, and subsequently, the model calibration are substantial tasks for achieving reliable simulation results. The calibration of numerical models as a subjective approach can become challenging and highly time-consuming, especially for inexperienced modelers, due to dealing with a large number of input parameters with respect to hydraulics and sediment transport. Using optimization methods can notably facilitate and expedite the calibration procedure by reducing the user intervention, which results in a more objective selection of parameters. This study focuses on the application of four different optimization algorithms for calibration of a 3D morphodynamic numerical model of a curved channel. The performance of a local gradient-based method is compared with three global optimization algorithms in terms of accuracy and computational time (model runs). The outputs of the optimization methods demonstrate similar sets of calibrated parameters and almost the same degree of accuracy according to the achieved minimum of the objective function. Accordingly, the most efficient method concerning the number of model runs (i.e., local optimization method) is selected for further investigation by setting up additional numerical models using different sediment transport formulae and various discharge rates. The comparisons of bed topography changes in several longitudinal and cross-sections between the measured data and the results of the calibrated numerical models are presented. The outcomes show an acceptable degree of accuracy for the automatically calibrated models.

1. Introduction

Numerical models have become useful behavioral representatives of complex environmental systems. In water-related domains, the growth of knowledge about the underlying processes and recent developments regarding numerical solvers and computational techniques have increased the performance and speed of simulations. However, a significant challenge remains for modelers to minimize the misfit between simulation outputs and corresponding physical observations, which is necessary to obtain a reliable predictive model [1,2]. Hydro-morphodynamic models are characterized by a large number of input parameters dealing with a considerable amount of uncertainty. This uncertainty arises from the sophisticated behavior of environmental fluid systems, the simplified structure of models, implemented empirical equations, unknown boundary conditions, and imprecise input data. Some of the physical input parameters, which are subject to calibration, can only be quantified at specific locations over a limited period, and for some it is not physically feasible to measure them.

Modeling practice can be defined as a four-step task, involving model setup, calibration, validation, and application [3]. Consequently, the model accuracy depends not only on the model structure and the quality of input data but also on the model calibration and validation. Calibration of numerical models is the foremost step for quantifying and accrediting computational simulations, involving two major stages: parameter specification (i.e., selection of sensitive parameters that are subject to adjustment) and parameter estimation (i.e., determination of optimal or quasi-optimal values of specified parameters). Calibration is an inverse, multistep problem, with the aim of uncertainty diminution by updating the model through subsequent comparisons between observations and computational results, in order to achieve a good agreement with a reasonable tolerance [4]. In other words, calibration is the process of model adaption by adjusting independent variables to gain accordance between computed and measured distributions of dependent variables, so that the model output deviates marginally from the real data specified in the performance criteria [5,6]. Nevertheless, even calibrated models potentially involve a certain amount of uncertainty and seldom address the underlying deficiencies because of both model errors (modeling assumptions, simplifications, and approximations) and data errors (the lack of error-free measurements) [7,8]. Hence, a calibrated model needs to be validated before applying it to practical problems to ensure reliable predictions [9,10,11]. On the other hand, according to Refsgaard and Henriksen [12], a validated model operates in a precise manner with regards to site-specific applications and predefined accuracy criteria, and thus the validity of a model is always limited in terms of time, space, boundary conditions, and the type of application.

Model calibration can be accomplished manually, automatically, or by using a multistep method combining the two approaches [13,14,15]. The most widely used trial-and-error method is not only highly dependent on users’ knowledge of the model structure and their level of expertise, but also on their understanding of the environmental system characteristics and the properties of measured data [16]. The trial-and-error method follows the simple approach of manual adjustment of uncertain or unknown parameters and the comparison of predictions with measured values, involving human judgment to attain the best fit for model parameters [17]. As the model structure becomes complicated, the manual calibration procedure becomes laborious, time- and cost-consuming, and for models with a large number of uncertain or unknown variables (e.g., morphological models), sometimes impractical. The demanding practice related to the manual method has been persuading modelers to improve the complex inverse calibration technique based on optimization algorithms (e.g., deterministic, metaheuristic, stochastic, or uncertainty-based). This can be done by coupling the model with an optimization engine, with the aim of speeding up the calibration process and establishing an objective scheme to address the “user-subjectivity” issue [18,19,20]. As an example of the user influence on model calibration results, Botterweg [21] compared the outcomes of a hydrological model that was calibrated independently by two users with an identical set of measured data and reported how different sets of calibrated parameter values could yield reasonable results.

Principally, automated calibration consists of three elements: an objective function to assess the differences between model outputs and observations, an optimization algorithm for sequential adjustment of preselected model parameters with regard to the reduction of the objective function’s value, and a convergence criterion [22,23,24]. Optimization algorithms can be categorized into two classes: global methods based on sampling the proposed values of parameters over the entire space, and local methods based on point estimation by finding the optimum point where no further progress can be achieved in the adjacent space of the parameter. Regarding the local methods, which are computationally efficient and need far fewer model runs, the initial values should be chosen carefully, as model calibration proceeds from this point towards the gradient descent of the objective function. If nonlinearity affects the model, there is no assurance that the inverse problem is unimodal, which means local methods can be simply trapped in local minima points instead of finding the global minimum. Therefore, since the calibration process is strongly dependent on the initial starting point, it is required that the user continuously assesses results, adjusts starting values, and restarts the model [25,26,27]. “Equifinality”, as described by Beven [28,29], should also be considered regarding the automated calibration, which may propose the same model prediction by using different parameter sets (also see Straten and Keesman [30] for the term “equally possible” parameter sets). In other words, uncertainty-based calibration methods such as Generalized Likelihood Uncertainty Estimation (GLUE) [31,32,33] or the Metropolis–Hastings algorithm of Markov Chain Monte Carlo (MCMC) method [34], which sample different parameter combinations to produce a calibrated probability distribution, detect several feasible parameter sets, unlike the other methods, which focus on pinpointing a single optimal solution [35].

The performance of autocalibration techniques has been extensively studied for environmental models in the fields of hydrology and groundwater over several decades. However, in the field of fluvial hydraulics, there are relatively few studies available in the literature regarding the application of automated inverse methods for calibrating hydro-morphodynamic models [36,37,38,39,40,41,42]. In the present study, a 3D hydro-morphodynamic model of a 180° curved channel is developed according to a physical model and calibrated against the measured bed elevations of 54 points along different longitudinal and cross-sections. The model calibration is carried out with a local optimization method (Gauss–Marquardt–Levenberg method), using the Parameter ESTimation (PEST) package [43] and three population-based global optimization algorithms: Particle Swarm Optimization (PSO) [44], Shuffled Complex Evolution (SCE-UA) [45], and Covariance Matrix Adaptation Evolution Strategy (CMA-ES) [46]. The used methodology is further explained in Section 2. This work aims to evaluate the efficiency of the automatic calibration procedure based on the local optimization method as compared to the global methods, with the main objective of predicting a single set of optimal parameter values for the model, which is presented in Section 3.1. Further, in order to test the performance of the most efficient method, different experimental runs with various discharge rates are used to set up and calibrate additional models, presented in Section 3.2. Finally, the results are summarized and conclusions are presented in Section 4.

2. Methodology

2.1. Experimental Data

Yen and Lee [47] conducted a series of experiments to investigate the bed topography evolution and transverse sediment sorting in a 180° curved channel under unsteady flow conditions. The rectangular U-shaped channel measured 1 m in width and had a 4 m central radius. Two straight up- and downstream reaches measuring 11.5 m in length were attached to the bend to ensure uniform and fully turbulent flow. A 20 cm flat layer of nonuniform sand with a median diameter of d₅₀ = 1 mm and a geometric standard deviation of ơ_g = 2.5 was placed on the bed, having a final slope of 2‰. Further details regarding the sediment characteristics can be seen in

Table 1. Sediment characteristics of the experiments that were applied in the numerical model.

Sediment Size Classes
Size (mm)	0.25	0.42	0.84	1.19	2.00	3.36	4.76	8.52
Proportion (%)	6.55	10.56	25.36	15.06	20.11	13.02	4.88	4.46
Cumulative proportion (%)	6.55	17.11	42.47	57.53	77.64	90.66	95.54	100
Fall velocity (m/s)	0.03	0.06	0.11	0.14	0.20	0.26	0.32	0.43

, which also includes the fall velocity of particles calculated according to the empirical equation of Ferguson and Church [48], by assuming ρ_s = 2650 kg/m³ as the density of sediments. The base flow rate for all experiments was set to Q₀ = 0.02 m³/s, according to the incipient motion of sediment particles, with a flow depth of h₀ = 5.44 cm. During the experiments, the discharge was linearly increased up to a maximum value and then progressively lowered to the base flow discharge (triangular-shaped hydrograph). In total, five different experiments were conducted with different durations and peak discharge values.

2.2. Numerical Model

Due to the three-dimensional nature of the flow field in bends caused by the interaction of centrifugal force and the lateral pressure gradient, which results in a helical flow structure, the use of a 3D numerical model is imperative. In the present study, the fully three-dimensional numerical model SSIIM 2 (Sediment Simulation In Intakes with Multiblock option) [49] is employed to simulate hydraulics and morphological bed changes. This Computational Fluid Dynamics (CFD) code has been successfully used to model several cases with different scales in the fields of hydraulic and sedimentation engineering by many researchers [50,51,52,53,54,55,56,57,58,59,60,61]. The numerical model SSIIM 2 solves the Reynolds-averaged Navier–Stokes (RANS) equations together with the continuity equation (Equations (1) and (2)) on an adaptive three-dimensional non-orthogonal and unstructured grid, by using a finite volume approach for the spatial discretization and an implicit scheme for the temporal discretization to compute the water motion for turbulent flow.

$$\frac{\partial U_i}{\partial t}+U_j\frac{\partial U_i}{\partial x_j}=\frac{1}{\rho }\frac{\partial }{\partial x_j}\left(-P{\delta }_{ij}-\rho \overline{u_i{u}_j}\right)$$

(1)

$$\frac{\partial U_i}{\partial x_i}=0i=1,2,3$$

(2)

where $U$ is the averaged velocity over the time $t; x$ is the space direction; $\rho$ is the water density; $P$ is the dynamic pressure; ${\delta }_{ij}$ is the Kronecker delta, which can be either 1 if $i=j$ or 0 if $i\ne j$ [62]; and $-\rho \overline{u_i{u}_j}$ is the turbulent Reynolds stress term. Regarding the RANS equations, the convective term is modeled by using the second-order upwind (SOU) scheme [49]. The semi-implicit method for pressure-linked equations (SIMPLE) is used to evaluate the pressure term [63] by computing pressure-correction from the water continuity defect in a cell. The Reynolds stress term is computed by the concept of eddy viscosity with the standard k–ε turbulence model [64], using the Boussinesq approximation. The free water surface is calculated according to the computed pressure gradient using the Bernoulli equation [65]. A Dirichlet boundary condition is applied for the inflow and a zero gradient boundary condition is defined for variables at the outflow boundary. Wall laws for rough boundaries are used for the bed and the side walls [66].

The grid size is chosen according to a similar study by Fischer-Antze et al. [67], containing 254 × 20 × 5 cells in streamwise, lateral, and vertical directions, respectively. The number of cells in the adaptive grid is adjusted during the computations according to the implemented wetting–drying algorithm [68]. Three different bedload transport formulae by Engelund-Hansen [69], van Rijn [70], and Wu [71] are used for the sediment transport computation. Using Wu’s correction factor, the hiding and exposure effect of nonuniform sediments is applied for the formulae by Engelund-Hansen and van Rijn, which by default do not consider the interaction between different size fractions in calculating the fractional transport rate of nonuniform sediments (i.e., they take the transport of each size individually into account).

2.3. Calibration Procedure

2.3.1. Local Optimization Method

The model-independent Parameter ESTimation tool (PEST) has been successfully used for calibration, sensitivity, and uncertainty analysis of environmental models in the fields of hydrology [72], groundwater [73], ecology and water quality [74,75,76], and flood modeling [77], providing promising results. PEST uses the gradient-based Gauss–Marquardt–Levenberg (GML) [78,79] algorithm for nonlinear parameter estimation through iterative progress refinement and by employing the weighted sum of the squared residuals as the objective function (Equation (3)).

$$SSR={\displaystyle \sum }_{i=1}^n{w}_i{r}_i^2$$

(3)

where $r$ is the residual, defined as the difference between model-generated bed levels and corresponding measured topographic elevation; and $w$ is the weighting factor, which is considered to be 1 for all observations in this study.

The GML algorithm combines the gradient descent method (when the function is far from the minimum) with the Gauss–Newton algorithm (when the function is around the minimum), resulting in a fast and efficient convergence to the minimum of the objective function. PEST invokes the model by using the initial parameter values prespecified by the user. Through Taylor’s expansion of the actual parameter set, PEST linearizes the relationship between input parameters and model outputs. At the first iteration, according to the number of adjustable parameters, PEST runs the model several times and fills the Jacobian matrix containing the partial derivatives of each model output (i.e., filling the Jacobian matrix requires as many runs as the number of adjustable parameters). Subsequently, the upgrade vector is calculated (Equation (4)) with the aim of defining a damping factor, and parameters are changed within the user-defined bounds. This procedure is repeated for each iteration by using the updated values from the previous one until the minimum of the objective function is found.

$$u={\left(J^{T}Q J+\lambda I\right)}^{-1} J^{T}Qr$$

(4)

where $J$ is the (n × m) Jacobian matrix, which has n columns according to the number of adjustable parameters and m rows following the number of model outputs; the superscript T stands for the transpose operator; Q is an (m × m) diagonal weight matrix; λ is a damping parameter (i.e., Marquardt lambda), which is adjusted during each iteration; I is an (n × n) identity matrix; and r is the vector of residuals for the current parameter set.

A reliable estimation of calibrated parameters in PEST highly depends on the predefined starting values. The reason is that in the case of getting trapped in a local minimum point on the objective function surface, the GML algorithm may not be able to leave this region. Therefore, in this study, the model calibration is done with four different initial parameter values for PEST to ensure its stability in finding the global minimum.

2.3.2. Global Optimization Methods

• The Shuffled Complex Evolution-University of Arizona (SCE-UA)

The SCE-UA method was initially developed for calibration of conceptual rainfall-runoff models [45], which uses a global optimization algorithm. The algorithm is based on the competitive evolution [80], the local direct search of downhill simplex method [81], a controlled random search [82], and the concept of complex shuffling. The optimization starts by sampling a population of randomly distributed points over the feasible parameter space and partitioning the points into several groups (complexes). Each complex contains 2n + 1 points (n is the number of parameters to be optimized) and evolves independently to explore the domain in different directions, using the downhill simplex algorithm. Periodically, the population is shuffled and candidate solutions are reassigned to new complexes to share the information obtained from previous complexes during their evolution phase. The evolution and shuffling are progressively repeated until the convergence criterion is reached (i.e., the entire population moves gradually to converge towards the global minimum point). For more details about the algorithm see [83,84].

• Covariance Matrix Adaptation Evolution Strategy (CMA-ES)

CMA-ES is a derivative-free stochastic method for real-parameter optimization of nonlinear and nonconvex functions, using the general operators of evolution strategies, i.e., recombination, mutation, and selection. Each iteration (generation) is initiated by sampling a new population of individuals or candidate solutions from a multivariate normal search distribution. This is followed by the evaluation of individuals according to their performance measures and by updating the parameters of the multivariate normal distribution to guide the search towards the area with lower objective function values. The sequence of iterations is continued until a user-defined termination criterion is met. CMA-ES adapts the covariance matrix representing the pairwise dependencies between the variables in the distribution. The algorithm follows two principles for the parameter adaptation: (a) the maximum likelihood principle for updating the mean of the normal distribution in a way that the likelihood of previously successful candidate solutions is maximized; (b) using two evolution paths, one for the covariance matrix adaptation and the other one for adjusting the mutation strength (or the step size, which is defined as the standard deviation of the normal distribution). For a detailed description of the algorithm see [85].

• Particle Swarm Optimization (PSO)

The PSO algorithm, inspired by the social behavior of swarming animals (e.g., bird flocks or fish schools), is a metaheuristic stochastic global optimization approach. In this method, the search space is initialized with randomly positioned particles as potential solutions forming a swarm. Particles are composed of two vectors (position and velocity vectors) and move in random directions, keeping track of the best position (solution) found by themselves (having intelligence) and the one achieved by the whole swarm (with social interaction). The fitness of particles is measured by an objective function. In each generation, the velocity vector (reflecting the travel direction of particles across the search trajectory) is updated according to the previous best position of each particle (personal best or cognitive component) and that of the entire swarm population (global best or social component). Finally, the new position is computed with regard to the updated velocity. The process is iterated until a termination condition is satisfied [44,86].

2.3.3. Investigated Parameters

Among several uncertain input parameters in SSIIM 2 numerical model, the three most sensitive parameters are selected for calibration through a manual sensitivity analysis as follows:

• Roughness height at the bed (k_s): For fixed beds, this parameter is typically assumed to be proportional to the representative grain size d_n (the diameter where n% of sediment grains are finer) (i.e., Nikoradse’s equivalent grain roughness). For movable beds, however, the roughness caused by bedforms has to be added to the grain roughness, which may increase it with a higher factor than the grain roughness itself. A collection of different equations regarding the roughness height can be found in the literature (e.g., [87]). In this study, despite the dynamic nature of the roughness coefficient due to the formation of bedforms, the roughness height is considered to be uniform along the whole domain, because of the focus on the automatic calibration procedure. The range of this parameter is selected to be between d₅₀ and 10d₉₀.
• Active layer thickness (ALT): This parameter is described as a function of the representative grain diameter and the bedform height as a dynamic value that depends on the sediment properties and the flow conditions. ALT determines the maximum depth of erosion during one time-step in the numerical model. For this study, a constant value of ALT is also chosen, with a range between d₅₀ and 5d_max [88,89].
• The volume fraction of compacted sediments (VFS): This parameter describes the proportion of deposited sediments in the bed compared to the water content, which depends on the bulk density as a function of grain size distribution and packing of sediment depositions. This parameter’s range is adjusted between 40% and 60% in the calibration process.

2.4. Work Structure

Among the five experiments conducted by Yen and Lee [47], three different discharge rates—Run#1 (the highest discharge), Run#3 (the middle discharge), and Run#5 (the lowest discharge)—are used in this study for model calibration in order to cover different discharge characteristics. First, Run#5 is selected to test all of the optimization algorithms. Afterwards, the most efficient method is also applied for Run#1 and Run#3 to test the performance of the selected method. It should be mentioned that Run#1 and Run#3 show the formation of bedforms at the outer bend in the experimental setup. However, the employed numerical model in this study cannot precisely simulate the bedforms.

Table 2. Characteristics of the three hydrographs used in this study.

Run#	Peak Flow Discharge (m3/s)	Peak Flow Depth (cm)	Duration (min)
1	0.0750	12.9	180
3	0.0613	11.3	240
5	0.0436	9.10	420

presents the features of the three applied hydrographs.

3. Results and Discussion

3.1. Testing the Performance of Different Optimization Methods

In order to evaluate the performance of the four selected calibration methods, in the first step, the numerical model of Run#5 is calibrated as a test case, using Wu’s bedload transport formula, which considers the fractional transport of nonuniform sediments. The initial values of the investigated parameters are defined for all methods to be: k_s = d₉₀, ALT = d_max, and VFS = 50%. However, for ensuring the survival of the gradient-based method by using PEST from being trapped in a local minimum point, three additional starting values are also proposed as:

• PEST#2 (k_s = d₉₀, ALT = 2d_max, and VFS = 55%)
• PEST#3 (k_s = 2d₉₀, ALT = 4d_max, and VFS = 45%)
• PEST#4 (k_s = 5d₉₀, ALT = 2d_max, and VFS = 50%)

Table 3. Calibration results for the numerical model Run#5 by local and global optimization methods, using Wu’s sediment transport formula.

Calibration Results	Calibration Method
	PEST#1	PEST#2	PEST#3	PEST#4	SCE-UA	CMA-ES	PSO
ks (cm)	0.904	0.894	0.899	0.903	0.898	0.892	0.904
ALT (cm)	1.639	1.644	1.629	1.638	1.635	1.631	1.647
VFS (%)	48.1	47.9	48.4	48	48.1	48.1	47.9
Minimum of SSR (cm)	0.28249	0.28251	0.28240	0.28249	0.28247	0.28247	0.28241
Number of model runs	26	30	28	27	176	440	448

Note: k_s, roughness height at the bed; ALT, active layer thickness; VFS, the volume fraction of compacted sediments; SSR, sum of the squared residuals.

summarizes the values for the calibrated parameters obtained from the four PEST runs (four different starting conditions) and those obtained by the three global optimization methods. In addition, the minimum of the achieved objective function (sum of the squared residuals, SSR) and the number of models runs are presented in the table.

According to the results presented in Table 3, PEST returns very similar outputs by using different initial values. This reveals that the software explored the whole parameter space and did not stick to a local minimum point. By comparing the results of PEST with the three global optimization algorithms, a good agreement can be found between them regarding the calibrated parameter values, which proves the ability of PEST to calibrate the model globally by applying a local optimization algorithm. Nevertheless, by taking the number of model runs into account, considerably fewer runs are used by PEST, representing an advantage for the calibration of morphological models, for which long simulation time occurs.

3.2. Application of the Selected Calibration Method for Additional Numerical Setups

Based on its efficient performance in model calibration, PEST is selected to calibrate additional numerical models (i.e., Run#1, Run#3, and Run#5, with the sediment transport formulae of Wu, van Rijn, and Engelund-Hansen; nine model calibrations in total). From

Table 4. Calibrated values of the three selected parameters for nine simulations.

Calibration Parameters	Sediment Transport Formula
	Wu			Van Rijn			Engelund-Hansen
	Run#1	Run#3	Run#5	Run#1	Run#3	Run#5	Run#1	Run#3	Run#5
ks (cm)	1.52	1.34	0.90	0.63	0.61	0.37	0.48	0.31	0.25
ALT (cm)	2.04	1.85	1.64	2.20	1.94	1.24	1.31	1.12	0.74
VFS (%)	49	51	48	60	60	60	52	51	53

, which presents the calibration results for nine models, it can be seen that the active layer thickness (ALT), as well as the roughness height (k_s), gain higher values with increase of the flow discharge. It can be argued that the formation and evolution of bedforms along the channel bottom account for this effect.

The volume fraction of sediments (VFS) compared to the water volume is quite constant for the simulations using the formulae of Wu and Engelund-Hansen (around 50%). However, a value of 60% is found by employing the sediment transport formula of van Rijn, which corresponds to the upper predefined boundary for the search algorithm, indicating that the best solution for the optimization problem may lie outside the physically realistic bounds. As already mentioned, the selection of threshold values, between which the optimization algorithm adjusts parameters, should be taken according to a reasonable range for each parameter.

Figure 1, Figure 2 and Figure 3 show longitudinal sections of the initial and final bed levels after Run#1 (the highest discharge), Run#3 (the mid discharge), and Run#5 (the lowest discharge), respectively. The results are presented along the sections: (a) 10 cm from the inner convex bank, (b) the central line, and (c) 10 cm from the outer concave bank.

In Run#1, according to the measurements, the maximum deposition occurs near the inner wall at the 75° cross-section, which is well predicted in terms of position and magnitude by the numerical model using the formula of Wu. The model using Engelund-Hansen’s formula also predicts the location of this point at the same cross-section (75°), with just a minor discrepancy. Nevertheless, the model of van Rijn shows higher bed levels compared to the measurements. In this case, the position of the point bar is at the central cross-section (90°), having a slightly higher magnitude (Figure 1a). Along the central longitudinal-section (Figure 1b), the predicted patterns of bed elevations by all of the numerical models have a similar trend compared to the measurements, with a modest overestimation of the deposition for the last 40° range of the bend (140°–180°). Near the outer bend, the maximum scour depth in the experiment can be seen at the 165° cross-section, whereas the numerical models give this point at the 180° cross-section. In the second half of the bend (90°–180°), measurements show larger fluctuations caused by the development of bedforms compared to the smoother patterns obtained from the simulations (Figure 1c).

Regarding Run#3, near the inner bank, the maximum calculated deposition point is slightly shifted to the upstream compared to the experiment. The location of this point is at around 60°–70° according to the numerical models, although measurements show it at about 75° (Figure 2a). Additionally, in the range of 90°–160°, higher bed levels can be seen in the numerical models compared to the experiment. Along the central line (Figure 2b), numerical models marginally underpredict the erosion part up to the 60° cross-section. However, the overall patterns of the bed levels are similar to the experiment. Although having almost identical patterns near the outer bend (Figure 2c), a contrast can be found in the range of the last 30° of the bend between the simulations and measurements, especially the underestimation of the maximum erosion depth at the 180° cross-section.

In Run#5, similar to Run#3, the position of the maximum bed level near the inner wall (Figure 3a) is shifted further to the upstream in the numerical models. Here, the results of the model using Engelund-Hansen’s formula diverge largely from the measured data. Although showing some degree of underprediction regarding the erosion part up to the 60° cross-section, the simulated bed levels along the central line and near the concave bank are in good agreement with the experimental data (Figure 3b,c).

In Figure 4a–f, the cross-sectional bed profiles for all three runs are plotted at the locations of the maximum deposition height (Run#1: 75°; Run#3: 75°; Run#5: 90°) and the maximum scour depth (Run#1: 165°; Run#3: 180°; Run#5: 180°), according to the experiments. It can be seen that the results of the numerical models are in good agreement with the experimental data regarding the pattern and magnitude of the bed levels. However, the models using the formula of Engelund-Hansen are not successful in accurately predicting the steepness of the scour flank and its depth in the downstream part while increasing the discharge rate (compare Figure 4b,d,f).

According to

Table 5. Goodness of fit for the calibrated models with different sediment transport formulae.

Goodness of Fit	Sediment Transport Formula
	Wu			Van Rijn			Engelund-Hansen
	Run#1	Run#3	Run#5	Run#1	Run#3	Run#5	Run#1	Run#3	Run#5
R2 (-)	0.90	0.89	0.95	0.90	0.89	0.94	0.88	0.83	0.90
RMSE (cm)	1.49	1.04	0.67	1.57	1.08	0.69	1.63	1.33	0.81

Note: R², the coefficient of determination; RMSE, the root mean squared error.

, which presents the overall performance of the models based on the coefficient of determination (R²) and the root mean squared error (RMSE) statistics, all three calibrated models using Wu’s formula have the best agreement with the observations, with R² values ranging between 0.89 and 0.95 and RMSE values ranging between 0.67 cm and 1.49 cm. The simulation results with van Rijn’s formula show similar statistical performance as those obtained with Wu’s formula. Nevertheless, the results for the models using Engelund-Hansen’s formula have the lowest correlation and the highest error compared to the measured data.

In Figure 5, the bed level changes normalized by the base water depth (ΔZ/h₀, where h₀ = 5.44 cm) are presented for all runs using Wu’s formula (right side), and are compared with the experimental data (left side).

4. Summary and Conclusions

In the present study, four different optimization methods (GML, SCE-UA, CMA-ES, and PSO) are employed for automatic calibration of a 3D morphodynamic numerical model of a 180° bend. Three input parameters, namely roughness height, active layer thickness, and the sediment volume compared to water at the bed, are selected for calibration, using their initial values and a reasonable range based on the literature. The results of the model calibration using all optimization methods (a local method and three global methods) show very similar calibrated values. All of the methods facilitate the calibration process by reducing the user intervention. However, considering the convergence speed, using the GML algorithm of PEST software is considerably more efficient compared to the other methods, as it needs far fewer model runs. Accordingly, PEST is selected for further investigation (In total, nine models are calibrated by using three different discharge rates and the sediment transport formulae of van Rijn, Wu, and Engelund-Hansen). Nonetheless, as PEST uses a gradient-based algorithm, it potentially involves a high risk of being trapped in a local minimum point over the search space rather than finding the global minimum. Therefore, the initial parameter values should be taken with care. It is also worthwhile reassessing the calibration procedure with different starting values to ensure that PEST finds the global minimum of the objective function. Moreover, characters and discontinuous values, which are defined in models (such as the selection of the sediment transport formula), cannot be processed by PEST and have to be adjusted manually.

The final bed levels predicted by the calibrated numerical models are compared with the measurements at various longitudinal and cross-sections. The overall capabilities of the numerical models are evaluated using the coefficient of determination and the root mean squared error statistics. It is concluded that the calibrated models using Wu’s formula can predict the general characteristics of bed deformations (regions of deposition, erosion, and their magnitudes) better than the other two formulae by having the highest correlation and the lowest error between simulations and experimental data. The formula of van Rijn also gives reasonably acceptable results. However, the models which use the formula of Engelund-Hansen have the highest disagreement with the observations. It can also be mentioned that the bedforms, developed at the outer bend by increasing the discharge rate, cannot be accurately simulated by the numerical model. This effect is independent of the selection of the sediment transport formula and the calibration routine.