Efficient number of calibrated cross sections bottom levels on a hydrodynamic model using the SCE-UA algorithm . Case study : Madeira River

Hydrodynamic models are important tools for simulating river water level and flow. A considerable fraction of the hydrodynamic model errors are related to parameters uncertainties. As cross sections bottom levels considerably affect water level simulation, this parameter has to be well estimated for flood studies. Automatic calibration performance and processing time depend on the search space dimension, which is related to the number of calibrated parameters. This paper shows the application of the Shuffled Complex Evolution (SCE-UA) optimization algorithm to assess the number of cross sections bottom levels used in calibration. Also was evaluated the extent of algorithm exploration regarding computational processing time and accuracy. It was tested the calibration of 2, 4, 7 and 10 cross sections bottom levels (2PAR, 4PAR, 7PAR and 10PAR calibration configurations) of a 1,100 km reach of the Madeira River. 7PAR and 10PAR representation had better fitness (lower objective function value) on cross sections used for calibration; however, the error on other cross sections (2 validation gauging stations) was higher than 2PAR and 4PAR calibration. The short number (5) of gauging stations used in calibration has limited the number of calibrated parameters to represent adequately the river level profile. Finally, this paper shows a contribution for the parsimonious selection of parameters regarding the spatial distribution of observation sites used in calibration.


INTRODUCTION
Hydrodynamic models are important tools for understanding river dynamics (PAIVA et al., 2013;PONTES et al., 2015); discharge, water level and flood forecasting (FAN et al., 2014a) and simulating pollutant dispersion (CAVALCANTE; MENDES, 2014).For large rivers, these models are often based on one-dimensional momentum and mass conservation equations, named as the Saint Venant equations.Saint Venant full resolution (PAIVA et al., 2013) and some simplified models, such as the 1D Inertial model (YAMAZAKI; ALMEIDA; BATES, 2013;FAN et al., 2014b), can successfully account river water level and adequately represent backwater effects, storage volume and river-floodplain water exchanges on river basin.
The hydrodynamic model performances are considerably dependent to the chosen set of parameters values.Cross section wetted area, given by its width and depth, bottom slope and channel roughness control the modeled flow state variables.For a large-scale basin, data from remote sensing can be helpful for parameter estimation, as field measurements are insufficient to embrace the whole river.Aerial images can facilitate defining river width (PAVELSKY;SMITH, 2008), and Digital Elevation Models (DEM) contributes identifying the drainage system (SIQUEIRA et al., 2016a).On the other hand, parameters such as the river bottom level carry out higher uncertainty as no remote sensing techniques for deep rivers are available other than local measurements (p.e. using ADCP).Therefore, the hydrodynamic model usually is submitted to a calibration process, which is a parameter adjustment to improve model fitness to observation data.
Although most of model applications focus on river roughness coefficient calibration -since it is a subjective coefficient that simplifies forces related to head losses - Paiva et al. (2013) had shown that errors in river bottom level could affect significantly simulations of water surface elevation, flood wave travel time and floodplain area extent.In addition, Wood et al. (2016) calibrated the LISFLOOD-FP model using SAR (synthetic aperture radar) images and observed that simulated flood extent is more sensible to variation on bathymetry than roughness coefficient.Indeed, in large rivers, cross section information (e.g.bottom levels) is sparse or missing and then could be incorporated to the calibration process as a set of parameters.
The calibration process can be manual or automatic.Manual calibration depends on the model user expertise, which selects the next set of parameters values based on the previous results mainly by visual analysis (BOYLE; GUPTA; SOROOSHIAN, 2000).As a negative point manual calibration can be exhaustive and time-consuming, especially in high-complexity models, and still could not reach true optima parameter values.On the other hand, automatic calibration algorithms evaluate model results through a mathematical criteria and search for better fitted parameters sets by optimization methods, hence eliminating subjectivity.
Automatic calibration algorithms can either be of local search or global search.Local search techniques are designed to find local optimum thus becomes deeply dependent of the startup parameters sets.Global search techniques commonly use a probabilistic approach in order to improve exploration over the search space, consequently increasing their chances to find the global optimum instead of converging to the closest optimum (region of attraction) as local search techniques (DUAN, 2003).
In 1992, Duan, Sorooshian and Gupta presented the Shuffled Complex Evolution -University of Arizona algorithm (SCE-UA) which was proved to be a robust global search calibration method.This algorithm was introduced in a CRR (Conceptual Rainfall-Runoff) model calibration context, being initially tested for a six-parameters model (DUAN; SOROOSHIAN;GUPTA, 1992).Subsequent SCE-UA applications had kept a small number of model parameters to be calibrated, for example: six in Siqueira et al. (2016b), eleven in Blasone, Madsen and Rosbjerg (2007), two to thirteen in Duan, Sorooshian and Gupta (1994), sixteen in Zhang et al. (2009) and eighteen in Eckhardt and Arnold (2001).
However, the number of parameters used to represent a hydrological system can be much higher than a couple of dozens, especially due to the growing complexity of current models.High number of parameters increases the search space dimension and consequently the quantity of model simulations needed in the calibration process.Thus complex models calibration demands longer computational processing time, which might be unfeasible.Eventually calibration accuracy becomes limited by computational capacity.
To bypass this situation, correlated parameters can be adjusted in a fixed ratio in order to reduce the search space dimension (ECKHARDT; ARNOLD, 2001;BLASONE;MADSEN;ROSBJERG, 2007).Decreasing number of parameters might also help lowering costs since less field measurements for parameters estimation would be required.
Previous researches have evaluated the impact of the number of calibrated parameters on the model accuracy of distributed, lumped or statistical hydrological models (HER; CHAUBEY, 2015; SCHOUPS; VAN DE GIESEN;SAVENIJE, 2008).However it is necessary to relate these factors to computational time demand and extend the subject to include hydraulic models likewise.
For this purpose, the present paper evaluates calibration effectiveness by investigating how many cross sections bottom levels should be calibrated on a 1D hydraulic model aiming to balance between processing time and model accuracy.The studied area comprises an 1,100 km reach of the Madeira River, in the Amazon River basin, where four optimization configurations with different numbers of calibrated cross sections bottom level (2, 4, 7 and 10) were assessed.

METHODOLOGY Hydrodynamic modeling
The hydrodynamic model used in this study is based on the so-called inertial model (BATES; HORRITT; FEWTRELL, 2010) which is basically a simplification of the shallow water equations that ignores the advective inertia term from the conservation of momentum.Despite of being introduced recently, the inertial model had already been applied successfully in several studies (YAMAZAKI; ALMEIDA; BATES, 2013;FAN et al., 2014b;PONTES et al., 2015).The 1D inertial model can be summarized in two equations: Where Q is the discharge; A is the cross section area; lat q is lateral inflow by unit of length; g is gravity; y is the water surface level (given by depth + bottom level, which is the calibrated parameter); n is Manning roughness; R is the hydraulic radius; x and t are the independent variables that represent longitudinal distance and time, respectively.
The numerical approach to solve the 1D inertial model equation system is described by Fan et al. (2014b).It is an explicit scheme that calculates flow from the momentum equation and water level/depth from the continuity equation and uses an adaptive time step based on the Courant-Friedrichs-Levy condition.The modeled cross section is rectangular (wetted area = depth × width) and the water volume that overflows riverbanks behaves as water storage.
The case study was an 1,100 km reach of the Madeira river, starting at the city of Porto Velho and ending at the river's outlet on the Amazon River (Figure 1).The Madeira River basin drainage area is about 1.4 million km 2 , which accounts for 20% of the Amazon basin.The portion of the basin upstream of the modeled river reach corresponds to ¾ of the total drainage area.This river reach was selected for this case study because it has no significant hydraulic structures or point lateral inflows that would require complex modelling.
The modeled river reach was discretized in 211 reaches of approximately 5 km length.Each reach width was determined using a river mask obtained from near infrared images (band 5) of Landsat 8.The adopted Manning roughness coefficient was 0.030.As boundary conditions, it was used a water level time series at the outlet reach and a flow time series at the first reach upstream, both data obtained from Brazilian National Water Agency (ANA) gauging station -15940000 and 15400000, respectively.
A simplified method was used to estimate lateral inflow through a sinusoidal function proportional to the incremental drainage area of the respective reach, obtained by an approximation of the difference between two ANA flow gauging stations distant 600 km from each other (15400000 and 15700000 -Figure 1).The lateral inflow ( i lat q ) was calculated by the following equations and results are presented in Figure 2: ( ) Where t is the simulation time; max is the maximum and minimum discharge difference between the two gauging stations (15400000 and 15700000), i IDA is the reach i incremental drainage area; DA is the accumulated drainage area of the gauging stations; i L is the distance between the river reach and the gauge station 15400000; and 10 days relates to the delay in which the two gauging stations discharge series present the highest correlation.

SCE-UA
The SCE-UA (DUAN; SOROOSHIAN; GUPTA, 1992) algorithm is a global search technique that combines the concepts of controlled random search, complex shuffling and competitive evolution to minimize an objective-function (OF).This algorithm is able to overcome the main calibration problems such as: multiple regions of attraction, roughness surface due to discontinuous derivatives, parameter interdependency and varying sensitivity Efficient number of calibrated cross sections bottom levels on a hydrodynamic model using the SCE-UA algorithm.Case study: Madeira River (DUAN; SOROOSHIAN; GUPTA, 1994).A summary of this method is given on the following paragraphs.Firstly, a population of s points (set of parameters values) is randomly generated.Each point is evaluated by an OF which typically is related to the difference between the model results and the observed values.Then the population is ranked based on the OF values and divided into p complexes of m points.Each complex evolves alone within a generation, according to a Competitive Complex Evolution (CCE) procedure.After a generation, the complexes are grouped and shuffled to form new complexes in the next generation.This process continues until a convergence criterion is attached or after a pre-defined number of generations.
The CCE procedure is based on the Simplex optimization method (NELDER; MEAD, 1965).Initially q points from the complex are selected to form a subcomplex.The selection considers a trapezoidal probability distribution in which the most adapted points (lower OF values) have higher chances to be selected.Then a procedure similar to Simplex takes place.Firstly the algorithm calculates the centroid of the subcomplex points except the less adapted one (higher OF value), then new points are evaluated through a reflection step, a contraction step or a mutation step.This process is repeated β times followed by a complex shuffled for the next generation (Figure 3).For a detailed explanation of the SCE-UA algorithm the reader is kindly redirected to Duan, Sorooshian and Gupta (1992).
For the purpose of this paper most of the parameters of the SCE-UA algorithm was set to default values as in Duan, Sorooshian and Gupta (1994): m = 2n + 1, q = n + 1 and β = 2n + 1, where n is the number of model parameters to be calibrated.Only the number of complexes ( p = 2, 4 or 6), which is an important parameter related to the exploration of the search space, was previously tested to assess the algorithm convergence.It was considered that the calibration process reaches the convergence criteria when the relative difference between the population median and lower OF values is less than 0.5%.

Experiment design
The SCE-UA algorithm was applied to calibrate the bottom level of the cross sections of the hydrodynamic model.Four different configurations were tested during calibration of the cross section bottom levels:  , 21, 42, 63, 101, 139, 176, 188, 200 and 211 -Figure 4) had their bottom level calibrated.
The bottom levels of the remaining internal reaches were defined through linear interpolation of the closest cross sections selected for calibration.SCE-UA was tested 10 times for each value of p (2, 4 and 6) in order to verify the algorithm performance.
To compare results and to define the search space, an initial guess (IG) of the river bottom level was acquired reducing a maximum depth from the river surface level obtained from a 500 m resolution digital elevation model (DEM) of the Shuttle Radar Topography Mission (SRTM -FARR et al., 2007).The maximum depth was estimated by a geomorphologic equation, presented on Paiva et al. (2013), which assumes the maximum depth (h) as a function of the drainage area (DA): The search space was defined as 15 m above and 5 m bellow of the bottom levels initial guess (except for 2 parameters configuration, which was 15 m above and 10 m below).The search space limits is shown on the Figure 5.
The objective function of the calibration process was the root mean squared error (RMSE) based on data from five water level gauging stations from ANA -15400000, 15630000, 15700000, 15850000 and 15900000 (Figure 1).Every water level gauging station used in this paper was georeferenced to EGM08 geoid according to Moreira (2016).Data from these gauging stations were compared to the model results at the closest cross section, which were 211, 163, 91, 63 and 34 respectively.
The calibration process refers to a three years simulation period with 2 warm up months (Jan and Feb of 2002) and the remaining time was used for RMSE calculation (Mar/2002 to Dec/2004).Data from Jan/2005 to Jun/2010 were used to validate the results obtained from calibration (temporal validation).Additionally two other water level gauging stations (15490000 and 15860000) were used in order to find out if the model represents other reaches of the river satisfactorily (spatial validation).The gauging stations water level time series are presented on box plots at Figure 6.

RESULTS AND DISCUSSION
First, it was assessed whether the selected configuration of the SCE-UA algorithm was able to find the global optimum set of parameters values.For this purpose the SCE-UA algorithm  Efficient number of calibrated cross sections bottom levels on a hydrodynamic model using the SCE-UA algorithm.Case study: Madeira River was tested using different number of complexes ( p = 2, 4, 6) and four parameter configurations (n = 2, 4, 7 and 10).Each test was repeated 10 times to verify the performance.The Figure 7 exhibits the algorithm performance relating the number of hydrodynamic model runs the SCE-UA population fitness represented by the median of the OF values.It can be observed that as the number of calibration parameters increases, more hydrodynamic model runs were necessary for algorithm convergence.The calibration of 10 parameters (10PAR) needed approximately 2, 5 and 15 times more hydrodynamic model runs to converge when compared to the calibration of 7 (7PAR), 4 (4PAR) and 2 (2PAR) parameters, respectively.Regarding the number of complexes used for the SCE-UA algorithm configuration, it seems that the quantity of hydrodynamic model runs needed for convergence increases by a nearly constant value.For instance, when 4 parameters are calibrated the number of hydrodynamic model runs necessary for convergence was increased in approximately 300 for every 2 adding complexes, while in a 10PAR configuration this constant value is around 1000.About effective computational demand, the tests were performed on an 8 GB RAM Intel Core i5 CPU; each model simulation takes about 4 seconds.Thus, while a 2PAR-two-complexes calibration configuration takes about 7 minutes to converge (100 model runs), a 10PAR-six-complexes takes approximately 4.5 hours (4000 model runs).
Table 1 demonstrates the highest and lowest OF values from 10 calibration processes of each configuration.As observed, the OF (RMSE) value decreases as the number of parameters increases.However, only a slight difference in the OF value was observed between the calibration of 2 and 4 parameters, although the 4 parameters calibration needs 3 times more hydrodynamic model runs to converge (similar results when compared 10PAR and 7PAR).On the other hand, the 10PAR and 7PAR configuration calibration resulted in a 20% lower OF value compared to the 2PAR and 4PAR.Thus, increasing the parameters number can lead to either a slight or a high reduction of the OF value, although it consistently raises the number of hydrodynamic model runs.These facts enhance the discussion of which should be the number of calibration parameters since a significant increase on

Brêda et al. processing time does not guarantee a considerable improvement on model fitness.
Table 1 also informs the algorithm precision rate, which was defined as the number of times the calibration procedure reached the exact same parameters values as the best calibration parameters set (lowest OF).The 2PAR and 4PAR configuration calibration achieved the same parameters values in all calibration processes, which probably means that the OF surface roughness is low and it has only a main region of attraction.In contrast, the 10PAR calibration has not repeated a single set of parameters.Although the precision result looks disappointing at first sight, there is not much difference between the 10PAR calibration outcomes.The highest and lowest OF values differ only in 0.1%.That happens because the bottom levels of the downstream reaches do not seem to influence significantly the OF value (Figure 8) since there is a boundary condition that controls water elevation at the river mouth, creating a highly rough surface near the global optimum.Regarding to the 7PAR configuration, it can be seen that there is a second region of attraction (Table 1) which accounts for a local optimum with OF value of 0.6815 m; however setting the number of complexes to 4 basically ensures enough exploration to find the global optimum.
As has been noted, the number of complexes selected has made a slight impact on the results of the calibration.However, raising the number of complexes increased significantly the number of model simulations necessary for convergence.2PAR and 4PAR configuration calibration had always converged to the same set of parameters and the 10PAR calibration had reached almost the same OF value.Although a higher number of complexes contributed finding the global optimum of the 7PAR configuration, the local optimum OF value (0.6815 m) reached in 2 of 10 times, was still better than the value obtained with a lower number of parameters.Thus, it becomes questionable whether is necessary to raise exploration by increasing the number of complexes and consequently the chances of finding global optimum or just search for an acceptable local optimum and save computational processing time.In this study case, it seems sufficient to use two complexes for any calibration configuration aiming to reduce computational cost.
The calibrated bottom levels obtained in the four different configurations are presented in Figure 9.It can be observed that all configurations have converged to similar values of bottom levels at the "middle region", limited between 300 and 850 km from the river outlet.It probably means the bottom levels at the middle region that provide the most accurate water level simulation is somewhere near to the calibrated values.Since 10PAR and 7PAR have more degrees of freedom, especially at the extremities of the reach, some "bumps" were created on the final bottom level to improve level estimations at the gauging stations cross sections.10PAR calibration has built a small barrier upstream followed by a steep region, creating an unreal but effective bottom level (Figure 9).
As example of the results, observed and simulated surface water levels from August 2007 to December 2008 in each gauging station are presented in Figure 10.The difference between simulated and observed water levels is about few meters during the high and low seasons.The calibration process has virtually eliminated   Efficient number of calibrated cross sections bottom levels on a hydrodynamic model using the SCE-UA algorithm.Case study: Madeira River the model bias.However, the simulated levels amplitude is clearly overestimated.That is due to other parameters uncertainties such as Manning roughness or cross section width, which could be calibrated as well.Thus, the irregular bottom level obtained from the 10PAR calibration might be an attempt to reduce the simulated level amplitude at the gauging stations.
High values of RMSE from the initial guess bottom (IG) level simulation shows that calibration is essential for acquiring reasonable results as observed in validation outcomes (Table 2).Even a simple calibration configuration (2PAR) can significantly improve simulations of the water level when parameters uncertainties are considerable.Most RMSE values of calibrated models were below 1 m.On the other hand, the IG RMSE was higher, reaching 9 m on gauging stations 15400000 and 15630000, which significantly affects quantifying river-floodplain water exchanges.
Validation results were similar to calibration results at the gauging stations used as observation sites on the calibration process (light grey area of Table 2).The RMSE values have remained low and the calibration processes with more parameters have presented better results.10PAR mean RMSE was 0.2 m lower than 2PAR and 4PAR, which represents a 25% error reduction (Table 2).Therefore, the temporal validation has shown that increasing the number of calibration parameters improves the model accuracy related to the gauging stations used on the calibration process.
The situation described above inverts if it is also considered the gauging stations 15860000 and 15490000, which were not used during the calibration process (dark grey on Table 2), on the mean RMSE calculation (spatial validation).A parsimonious representation of the bottom levels revealed better performance on both gauging stations.2PAR and 4PAR calibrations have maintained almost the same RMSE values (which are around 1 m) while 7PAR and 10PAR have doubled or triplicated the RMSE values compared to other gauging stations (Table 2).Although the gauging stations 15490000 and 15400000 are relatively close (Figure 1), around 75 km distance, there are considerable differences between the RMSE values.These results confirm that the irregular bottom levels created by the 7PAR and 10PAR calibration on the river reach extremities deviate from the real bottom level.Therefore, if the model purpose is to simulate the river level profile or obtain water level series on specific points apart from gauging stations used in calibration, then the 4PAR or 2PAR calibration configurations (more parsimonious) are indicated.

CONCLUSION
This paper has focused on evaluating computational demand and accuracy of the calibration processes with different number of parameters (2, 4, 7 and 10).It is presented a case in which it is possible to discuss the main problems related to calibration The quantity of cross sections that should be selected for a hydraulic model calibration process is highly dependent of the observation spatial distribution.If there are only a few observation sites, the calibrated cross sections will not represent well the whole river.Indeed the calibration will adjust the parameters to provide the best fitness to the calibration gauging stations and neglect the real river level profile.In such cases, less calibrated parameters results in a more realistic river profile estimation.However if the model purpose is matching simulated to observation water levels in gauging stations, then calibration with more parameters can provide better results.
One way that could reduce level profile representation problem is forcing downstream bottom level to be lower than upstream.This is a parameter limitation that can avoid imaginary bottom levels as in the 10PAR calibration.Another option is to create a correlated bottom level as in Yoon et al. (2012), thus close cross sections bottom levels will not differ considerably.
Concerning this case study, the OF value at the local optima was very close to the OF value at the global optima.Thus reducing the calibration algorithm exploration becomes an option to decrease the number of simulations necessary for algorithm convergence and consequently reduce the computational demand.After all it is user task to balance between accuracy and processing time, and to be careful with overfitting and parsimony.
Although the bottom elevation accounts for the main errors in surface water level, alternative parameters could also be tested.Since the calibrated model simulations presented higher amplitude than observed levels, probably the remaining improvements in model fitness is related to calibration of other parameters such as Manning roughness or cross section width.Therefore, a mixed Manning roughness and bottom level calibration turn out to be a promising choice for achieving better outcomes.
These results enhance the knowledge about calibration settings.Her and Chaubey (2015) revealed that raising the number of calibration parameters on the hydrological model SWAT diminish parameters uncertainty and improve model accuracy, however the authors only tested few observation sites and did not consider parameters spatial distribution.On the other hand, the present study has shown that, regarding river cross sections on hydraulic models, a parsimonious number of parameters can provide more reliable results.
This study introduces the importance of selecting a sufficient number of cross sections to adequately represent the river channel on a hydraulic model.Nonetheless this paper presents only a case study.It is recommended that future researches evaluate rivers with different sizes and shapes to test if the cited conclusions apply to other situations.
In conclusion, this paper has presented contributions regarding parsimonious choices of calibration parameters, specially related to river bottom level.Fewer parameters might not only save processing time but also improves model performance.It is commonly assumed that the model would adjust better as the number of calibrated parameters increases; however, this quantity should rely on the number of observation sites and the parameters spatial correlation.

Figure 1 .
Figure 1.Studied reach and gauging stations position at the Madeira Basin.
a) 2 parameters (2PAR) -only the first and the last cross sections bottom level are calibrated (cross sections 1 and 211, Figure 4); b) 4 parameters (4PAR) -the first, the last and other reaches where most important tributaries join the main river (cross sections 1, 63, 176 and 211, Figure 4) had their bottom level calibrated; c) 7 parameters (7PAR) -the four previous ones and one reach between each of them (cross sections 1, 31, 63, 120, 176, 194 and 211, Figure 4) had their bottom level calibrated; d) 10 parameters (10PAR) -the four cross sections presented in b) and other two equally-spaced reaches between them (cross sections 1

Figure 7 .
Figure 7. Number of Model Runs necessary for convergence related to number of complexes (p) and the number of calibrated parameters.

Figure 6 .
Figure 6.River median level profile and level series box plot of the gauging stations/GS (dark -downstream boundary condition; grey -validation GS; white -calibration/validation GS).

Figure 8 .
Figure 8. Calibrated bottom levels of 10PAR the configuration.

Table 1 .
OF RMSE values and Precision Rate.

Table 2 .
Validation results (RMSE) (dark grey -only validation GS; light grey -calibraion/validation GS).In addition, this study highlights the importance of calibrating river bottom level for water elevation studies.The case study was an 1100 km reach of the Madeira River simulated by a hydrodynamic model.The configurations with 2 and 4 parameters have presented results that are more consistent and less time demanding.