**Guillermo J. Acuña 1,2,\*, Humberto Ávila <sup>1</sup> and Fausto A. Canales <sup>2</sup>**


Received: 25 May 2019; Accepted: 29 June 2019; Published: 5 July 2019

**Abstract:** Numerical models are important tools for analyzing and solving water resources problems; however, a model's reliability heavily depends on its calibration. This paper presents a method based on Design of Experiments theory for calibrating numerical models of rivers by considering the interaction between different calibration parameters, identifying the most sensitive parameters and finding a value or a range of values for which the calibration parameters produces an adequate performance of the model in terms of accuracy. The method consists of a systematic process for assessing the qualitative and quantitative performance of a hydromorphological numeric model. A 75 km reach of the Meta River, in Colombia, was used as case study for validating the method. The modeling was conducted by using the software package MIKE-21C, a two-dimensional flow model. The calibration is assessed by means of an Overall Weighted Indicator, based on the coefficient of determination of the calibration parameters and within a range from 0 to 1. For the case study, the most significant calibration parameters were the sediment transport equation, the riverbed load factor and the suspended load factor. The optimal calibration produced an Overall Weighted Indicator equal to 0.857. The method can be applied to any type of morphological models.

**Keywords:** calibration; river modeling; design of experiments; MIKE-21C model; Meta River

#### **1. Introduction**

Numerical models have become an essential tool for researching and developing engineering solutions related to water resources problems [1]. Models enable complex underlying processes to be captured, and facilitate the analysis of interrelationships between variables in cases with limited data [2]. Modeling has major applications in fields such as hydrology, maritime and coastal studies [3,4], river hydraulics [5] and water quality [6].

Calibration is one of the most important activities within the modeling process, because the model's credibility strongly depends on it [7]. The calibration process can be defined as adjusting the parameter values of a model in order to reproduce the real-world response within an accuracy range defined in the performance criteria (i.e., an acceptable level of adjustment between model and reality) [8,9]. The importance and impact of calibration on hydrological and hydraulics models has been assessed and confirmed by several authors [10–13].

According to Troy et al. [14], the process for calibrating the parameters required in a model can be classified in three categories: trial and error (manual) calibration, automatic optimization (generally through computer programs) and multistep automatic methods that take advantage of combining the strengths of manual and automatic calibration. The trial and error approaches provide more

control to the modeler; however, manual adjustment may become too complex when the parameters to calibrate are numerous and correlated. The automatic optimization methods are based on three elements: an objective function, an optimization algorithm (that usually includes a set of constraints) and a convergence criterion. This type of calibration has been widely used in recent works related to hydrological and hydrodynamic modeling [10,11,15–17]. However, this type of calibration faces a numerical problem, the equifinality of the problem (i.e., that there might exist a set of different combinations of parameters for which the objective function returns the same value) [18]. It is also possible that the combination of values obtained in the calibration process is inconsistent with real world phenomena that the model tries to represent.

In general, river models consist of a combination of two or three of the following components: hydrodynamics, sediments and morphology. A river model calibration is usually focused on the hydrodynamic component, whose usual indicators are water surface levels and discharge under steady flow conditions [19]. The most common indicator for sedimentological calibration is the sediment transport rate [20,21]; however, this information is not always available and, furthermore, it may include high levels of uncertainty depending on the measuring or estimation techniques employed to obtain the data. Morphology calibration typically involves comparing field data on bathymetry as well as erosion and sediment transport rates, with results obtained from simulation. This type of calibration requires a set of parameters that allows sensitivity analysis to be performed on the simulated river form, in order to emulate field data conditions from an initial state to a final state [22]. Morphology calibration is complex and rarely performed, because topobathymetric data is usually unavailable [3]. Many numerical models are unable to effectively reproduce the physical processes related to the morphological evolution of the river channel [1], and understanding the relationship between the river channel morphology and erosion processes is the subject of many recent studies [23–25] aiming for a better understanding and forecasting of the river evolution.

Because of the difficulties commonly caused by the lack of sufficient and accurate data, as well as by poor model performance under certain conditions, the calibration of a hydromorphological model is usually complicated due to the interaction of its calibration parameters; consequently, the sequential adjustment of each parameter value might not be the most efficient calibrating method [26]. This may be even more relevant for rivers whose morphology is heavily affected by hydrological and sedimentological variations, because the slightest modification in the parameters of one of the three components (hydrodynamics, sediments or morphology) would possibly affect the adjustment of the other two. For example, a significant change in the river cross-section could influence flow velocities, water surface levels and shear stress, which in turn would impact on erosion and sediment transport processes. In such cases, a common calibration approach used by most consultants and researchers is to follow recommendations given by software developers and more experienced users. Nevertheless, these recommendations do not represent a structured method, and they are also hard to replicate when different conditions arise [27].

Within this context, the main contribution of this paper is to present a calibration method based on Design of Experiments (DOE) theory that allows: calibration of the model considering the interplay between the calibration parameters; identification of the most sensitive parameters within the calibration process; and determination of a value or a range of values for which the calibration parameters produce an adequate performance of the model in terms of accuracy. The method consists of a systematic process for assessing the qualitative and quantitative performance of a hydromorphological numeric model based on calibration parameters (riverbed level changes, velocity vectors, sediment transport rates, etc.).

The combination of sensitivity analysis and optimization techniques for a better model calibration has been described by van Waveren et al. [27] as a good modeling practices. The DOE approach presented in this paper incorporates these good practices while also allowing the modelers to get a better understanding of the effect of the calibration parameters on the adjustment indicators, which is a useful feature, especially for beginners in modelling. Using DOE to define the number of simulations might even reduce the computational times when compared to other heuristic approaches. Because it is based on DOE theory, the method aims at easy implementation and adjustment.

To validate the method, this paper presents a case study evaluating a 75 km-long reach of the Meta River, in Colombia. The modeling was conducted with MIKE-21C, a two-dimensional flow model which can analyze spatial and temporal variation in depth, bed level, velocity, and shear stress during extended time intervals [28]. Nevertheless, the method described in this paper can be coupled with similar hydromorphological models, many of them briefly described in [29].

The remainder of this paper has the following structure: Section 2 reports the set of parameters required for calibrating the model, and descriptions of the method and case study; Section 3 presents the model setup and the main results obtained by applying the method, followed by the corresponding discussion. Section 5 lists the main conclusions.

#### **2. Materials and Methods**

The method is summarized as follows: (i) definition of the model objective and possible simplifications; (ii) selection the calibration parameters and indicators used for hydrodynamic, sedimentologic and morphologic components; (iii) experiment design and calibration. After presenting the method, the river reach taken as the case study is described.

#### *2.1. Modeling Objective and Simplifications*

Some of the typical modeling purposes include: understanding the hydrodynamics, sediment transport or morphological development of a specific segment of a river; aiding in the design of ports, bridges and other hydraulic and navigation structures; assessing river restoration and management plans; analyzing aquatic ecosystems; etc. Defining the model objective is perhaps the most important step of the process. Once the objective is defined, and before calibrating the model, it is necessary to reduce the number of calibration parameters, in order to have a better understanding and control over the process. Depending on the spatiotemporal resolution of the simulations, the use that will be given to the results and the characteristics of the river in consideration, some parameters may become less significant. For example, the river hydrodynamic conditions may be such that the eddy viscosity calibration becomes unnecessary, possibly because it is not a sensitive parameter in the hydrodynamic adjustment process [30]. Another example is the existence of hydraulic structures or marginal protections providing riverbank erosion control [31], which in turn make the erosion rate parameter expendable.

#### *2.2. Parameters Used as Hydrodynamic, Sedimentologic and Morphological Indicators*

For every parameter considered during the calibration process, and in order to determine the goodness-of-fit of the model, a set of qualitative and quantitative indicators is defined for comparing measurements with the results from simulations. Qualitative indicators refer to variables difficult to compare at a specific time or place, or variables whose interpretation requires modeler experience. For this method, the goodness-of-fit of the calibration based on the quantitative indicators is measured by means of the coefficient of determination, also known as R-squared (R2), a statistical measure of how close the data fit the regression line, in this case, comparing whether the simulations match real-data. The method admits other statistical measures for goodness-of-fit, like the mean squared error (MSE) and the mean absolute error (MAE).

Based on suggestions by Matte et al. [32] and some common outputs from river models, a set of quantitative indicators were chosen to compose an Overall Weighted Indicator (OWI) for assessing calibration. The OWI is expressed as follows:

$$\text{OVI} = \beta\_1 \cdot \text{WL} + \beta\_2 \cdot \text{QL} + \beta\_3 \cdot \text{FD} + \beta\_4 \cdot \text{MP} + \beta\_5 \cdot \text{SST} + \beta\_6 \cdot \text{SBT} + \beta\_7 \cdot \text{ST} + \beta\_8 \cdot \text{BL} + \beta\_9 \cdot \text{BE},\tag{1}$$

where *<sup>n</sup> <sup>i</sup>*=<sup>1</sup> β*<sup>i</sup>* = 1, and:


The different values of β*<sup>i</sup>* correspond to the relative weight of each indicator on the OWI. These weights are usually chosen by consensus among the modeling team based on experience or based on the uncertainty related to the variables of each indicator. For beginners, the authors of this paper recommend using the same value for all β*i*, as explained by Chaves and Alipaz [33] based on Shannon's principle of maximum entropy, which warrants that the probabilities of underestimating or overestimation the parameters would be the same, considering that the indicators and the parameters are random variables (with also random distributions).

It is worth mentioning that the OWI can be subdivided in terms of the components of the river model. Elements WL, QL, FD and MP are calibration parameters of the hydrodynamic component (OWIHD); SST, SBT and ST correspond to the sedimentological component (OWIST); while BL and BE are calibration parameters accounting for the morphological component (OWIMF).

The input parameters can be categorized into two types: those that can be determined from field measurements or obtained from literature or other reliable sources; and those that must be determined through calibration [8]. Table 1 shows some of the most common parameters required for modeling the hydrodynamic, sedimentological and morphological characteristics of a river. It is worth noticing that the calibration parameters listed in Table 1 emphasize the ones used by MIKE-21C, the computer model used as the main software tool in the present work.


**Table 1.** Typical calibration parameters. (Source: Adapted from [28]).

Depending on the river characteristics or the chosen computational model (or numerical approach), the model might require more or less parameters. For instance, MIKE-21C was employed by de Villiers [35] for studying cohesive sediment transport in a shallow reservoir, including within this analysis some additional calibration parameters such as Critical Shear Stress for deposition (τ*de*), Critical Shear Stress for erosion (τ*er*), Erosion constant (*E*0) and Exponent of the erosion. Beck and Basson [36] assessed the hydraulic and sedimentological behavior of the Klein river estuary by means of MIKE-21C, including as calibration parameters: Hydrodynamic time step, Morphological time step, Flooding Depth, Drying depth, Median grain diameter, Mass density of sediment and Porosity, as well as the calibration parameters shown in Table 1. Therefore, if additional calibration parameters are required for a more accurate representation of the process, the authors of the present paper recommend to the reader to conduct a literature review in order to define the range of feasible values for each calibration parameter. For example, the critical shear stress for deposition (τ*de*) and the critical shear stress for erosion (τ*er*) are very variable parameters which have a great impact on sediment transport mechanisms, but some authors [25,35] have performed sensitivity analysis on these calibration parameters, that could be used as baseline for setting the range of feasible values.

In this paper, two types of visual analysis are defined: (i) a comparison between simulated and measured (ADCP) velocity vectors, and (ii) an assessment of the simulated and historical morphological behavior of the river. Low performance model setups are easily detected by visual analysis before any numerical analysis.

#### *2.3. Experiment Design and Calibration*

The aim of DOE is to maximize the information obtained from a minimum number of experiments. It also helps to determine which factors might affect the performance of the model [37].

Once the calibration parameters and goodness-of-fit statistics have been chosen, the proposed method considers a screening design of experiment (definitive screening type) that allows identifying which calibration parameters and interactions have a significant effect on the overall accuracy of the model [38,39]. For this screening design, authors propose using a 2k factorial type of experiment, which evaluates *k* calibration parameters, each of them having two alternatives or levels that must be defined according to its corresponding physical or numerical meaning [37,38,40]. In consequence, goodness-of-fit and OWI calculations must be performed for each possible combination defined in the 2<sup>k</sup> experiment.

A second DOE aims at optimizing the calibration parameters whose effects (or the interactions with other parameter) were found as significant from the screening DOE. The present authors recommend a 3<sup>k</sup> type of experiment, that means, evaluating these parameters for at least three different representative levels from their previously defined range of possible values [38]. However, following this recommendation depends on available resources. As in the previous DOE, goodness-of-fit and OWI calculations must be performed for each possible combination defined in the 3k experiment.

Depending on the results from the second DOE, it is possible to obtain a regression equation expressing the relationship between calibration parameters and the model response, as well as the corresponding goodness-of-fit statistic, allowing optimizing a specific component of the OWI. The calibration process proposed in this paper is represented, for a better understanding, in the flowchart shown in Figure 1.

The process begins with hydrodynamic calibration. For this component, the main calibration parameter is the roughness coefficient, since water depth and flow velocity are function of this parameter, and it also has influence on flow distribution between the islands. If the river topology presents a cross-section contraction or sudden change of direction, it is likely that the eddy viscosity coefficient will be another significant calibration parameter [41].

**Figure 1.** Flowchart describing the calibration process.

The present method is suitable for calibrating and modeling rivers with a significant relationship between hydrodynamics and morphology. Therefore, it is necessary to include the effect of river morphological changes when conducting hydrodynamic calibration. For the method explained in this document, it is proposed the use of a morphological seed or preliminary values for morphological calibration parameters. This would allow a preliminary calculation of the morphological component of the model, during hydrodynamics calibration. The selection of the morphological seed can be done through a conceptual analysis, by using reference values available in the literature or based on known values used in other sections of the same river or other rivers with similar characteristics.

Once the hydrodynamic calibration is complete, the process continues with the sedimentological and morphological calibration. As explained before, such calibration is performed by an experimental design which can be used to define the significant parameters, the interactions between them, and the range of values for the calibrating parameters that provide a good fit between the results of the model and the observed data. When assessing the calibration of the morphological (OWIMF) and sedimentological (OWIST) components of the Overall Weighted Indicator, it is important to maintain the calibration of the hydrodynamic component and its corresponding parameters.

It is worth mentioning that the method can be used for unidimensional or multidimensional models (2D or 3D). A greater quantity of calibration parameters and indicators increases the required computational processing time, but it does not guarantee more accuracy.

#### *2.4. Case Study: A Reach of the Meta River, Colombia*

The Meta River is one of the major tributaries of the Orinoco river. The methods and results presented in this paper were produced from a project called "Update of studies and designs for navigation between Cabuyaro (K804) and Puerto Carreño (K0)" [42], conducted by Universidad del Norte (UNINORTE) in 2013, for the Colombian National Roads Institute (INVIAS) and the Colombian Ministry of Transport. From the field measurements made, the length of the Meta River was calculated as 1002 km, from its source near the town of Guamal, in the Meta Department, to its mouth in the Orinoco river. The Meta River watershed is approximately 99,500 km2 in area.

Recovering and improving navigability in Colombian rivers has been one of the main goals of the Colombian government over the last years [43]. For this reason, the aforementioned project was commissioned in order to analyze and model the Meta River, aiming at identifying the actions required to improve its navigability conditions.

The main data for this study consists of water depth, velocity and flow rate measurements, taken at the Orinoco river mouth (K0) and near the town of Cabuyaro (K796). The measurements were taken in two field campaigns that occurred between August 2012 and January 2013. Besides these on-site measurements, hydrological records from the Institute for Hydrology, Meteorology and Environmental Studies of Colombia (IDEAM) and Landsat (NASA) satellite imagery were other sources of information used for this case study.

The navigable channel of the Meta River is very unstable, and the flow and water depth variations through the year are significant and often abrupt, which causes some boats to become stranded. At the same time, the river presents high erosion and sedimentation rates [44]. From historical records between 1983 and 2010 at the IDEAM gauging station "Aceitico" (located 127 km upstream the mouth of the Orinoco), it was found that the average hydrograph presents maximum flow rates of around 10,323 m3/s, and minimum of around 755 m3/s. This ratio of 13.67 between maximum and minimum average indicates wide variability in flow rate in any given year. Furthermore, a unimodal trend was found in all hydrometric records; with high flow rates between May and August, and low flow rates between early December and late February.

Sedimentological measurements at the IDEAM stations in Puerto Texas (K669), Aguaverde (K360) and Aceitico (K127) show that sand and a small fraction of silt are the main materials transported by the river. The average grain size in the bed (D50-bed) is approximately 0.35 mm and the average grain size of the suspended particles (D50-susp) ranges between 0.057 and 0.17 mm. It is also worth noting that the flow rate in the river presents strong correlation with the sediment transport rate. According to the IDEAM records, this transport rate ranges between 419,680 and 4057 ton/day, based on data from Aceitico station for maximum and minimum average flow rates.

The results from the geomorphological studies carried out during the project showed that the Meta River morphology could be described as: 34% tabular, 26% sinuous / meandering, 24% straight, 10% braided and 6% anastomosed, based on satellite images from August and September 2012. Although the Meta River does not hold a particular morphological form, the selected satellite images, as well as aerial photographs taken between years 1986 and 2012 evidenced a high lateral and frontal mobility of the channel. From the analysis of the satellite images, it was also found that the floodplain of the river ranges from 3 to 9 km wide.

Based on the previous characterization and due to the technical and financial unfeasibility to model the whole river, a representative reach 75 km long was selected, between abscissas K235 (6◦06 18.81" N/69◦12 55.04" W) and K310 (6◦02 50.38" N/69◦44 40.11" W) (See Figure 2). The selection of this analysis section was based on:


**Figure 2.** Representative section of the river, between K235–K310.

2.4.1. Indicators for the Meta River Model

As described in the previous section, the calibration process begins by selecting the available parameters that can be used to describe the river behavior. The rating curves for Meta River were obtained from the unidimensional model presented by UNINORTE [42], and this information was used as boundary conditions for hydrodynamic calibration.

The following comparable elements, based on the measurements obtained in the aforementioned study, were used as indicators for the calibration process of the Meta River model:


$$Q\_S = 0.0096 \cdot Q\_L^{1.774} \text{\AA} \tag{2}$$

where *QS* is suspended sediment transport rate (ton/day) and *QL* is the river discharge (m3/s). To calculate R2, comparison was made between the *QS* values obtained by using discharge measurements and the ones calculated by using the results from the model.


For this case study, the OWI described in Equation (1) was divided into two components: The weighted indicator of hydrodynamic adjustment (OWIHD) and the weighted indicator of sedimentological and morphological adjustment (OWIST, MF). The β*<sup>i</sup>* coefficients were defined by the team of modelers participating of this project, based on agreed confidence levels for each variable measured during the field campaigns. Lower β*<sup>i</sup>* coefficients were assigned to variables with higher uncertainty levels.

$$\text{OVI}\_{\text{HD}} = 0.60 \text{ WL} + 0.40 \text{ FD} \tag{3}$$

$$\text{OWI}\_{\text{ST}, \text{MF}} = 0.30 \text{ SST} + 0.70 \text{ BL} \tag{4}$$

By considering that both components have the same weight, the global OWI can be expressed as:

$$\text{OVI} = 0.\text{30 WL} + 0.20 \text{ FD} + 0.15 \text{ SST} + 0.35 \text{ BL} \tag{5}$$

2.4.2. Parameters for Modeling the Meta River

For this case study, the modeling objective was to assess navigation capabilities in the Meta River. Based on this objective and the river characteristics, the modelers defined that:

• For the river reach under consideration, the Chézy coefficient is used for quantifying bed resistance. In MIKE-21C, this value is estimated as a function of water depth, following the approach explained by Talmon [46]:

$$\text{Chézy} = \mathbb{C} \cdot h^{0.17} \tag{6}$$

where Chézy is given in m1/2/s and h is the water depth in meters. Coefficient *C* (in m1/3/s) is the parameter to calibrate, and corresponds to the reciprocal of the Manning's roughness coefficient [47], considered constant by the authors for this study. According to Talmon [46], this approach allows to incorporate small bathymetric irregularities (due to dunes, ripples, etc.) as bed resistance. The adequate calibration of the bed resistance directly impacts on the accuracy of the model to correctly estimate flow distribution and morphologic evolution. For example, if the resistance over a shallow region is too high, too much flow will be deflected and there will be a greater tendency towards developing sandbars [35].


with vegetation, which favors stability. The erosion rate at the riverbank and the corresponding coefficients were not considered as calibration parameters.

• Due to absence of significant patterns and/or phenomena affecting the morphology, the helical flow coefficient HL was simplified to its default value of 1.00.

Based on the previous considerations, the calibration parameters for this case study can be summarized as:


### **3. Results**

This section describes the main results of the calibration process.

#### *3.1. Hydrodynamic Calibration*

To initiate the hydrodynamic calibration, it was necessary to define a sedimentological and morphological seed given as:


The hydrodynamic calibration was focused on the adjustment of coefficient *C* in the roughness equation (Equation (6)). As seen in Tables 2 and 3, three coefficient *C* values were evaluated for this purpose: 50, 55 and 60. As described in the previous section, their accuracy was assessed through the measured water levels and percentages of flow distribution through river branches as indicators. The best fit was obtained using *C* = 55, which results in an OWIHD = 0.9062, as shown in Table 3.

**Table 2.** Flow distribution through river branches around islands of the river reach.


**Table 3.** Hydrodynamic calibration results.


Once the hydrodynamic calibration has been conducted, the behavior of the velocity vectors is visually assessed, basically to verify if the model adequately represents this feature. An example of the model results for these visual comparisons is presented in Figure 3.

**Figure 3.** Velocity vectors from the model for hydrodynamic calibration.
