Calibration and prediction uncertainty analysis of a hydraulic-water quality coupling model using a modified moth-flame optimizer

Parameter auto-calibration is a hot topic in the realm of environmental modeling in this decade due to its high efficiency compared to manual calibration. As a result of good global optimum locating ability for high-dimensional parameter space and good parallelizability, heuristic algorithms have become one of the main streams in auto-calibration tasks. To date, the successful application of some classic heuristic algorithms in environmental model calibration has been reported by many studies. For example, the Shuffled Complex Evolution-University of Arizona (SCE-UA) proposed by Duan et al. (1992) has been widely adopted in hydrologic models (Yang et al. 2020), water quality models, and groundwater model calibrations. It was also encapsulated by the Uncertainty Quantification Python Laboratory (Wang et al. 2016) for more extensive use; Particle Swarm Optimation (PSO) algorithm (Kennedy & Eberhart 1995) is one of the most popular calibration algorithms for the Soil and Water Assessment Tool (SWAT) model (Liang et al. 2021). Although with such wide implementations, there are still two main problems that can be improved to pursue a more efficient calibration process and a more accurate and robust calibration result.

The first problem is the lag between the cutting edge of machine intelligence and its applications in environmental model calibration, especially for the model associated with complex modeling software. In this decade, several novel heuristic algorithms, for example, the grey wolf optimizer (GWO) (Mirjalili et al. 2014), moth-flame optimizer (MFO) (Mirjalili 2015), and sparrow search algorithm (SSA) (Xue & Shen 2020) have been proposed. These algorithms usually outperform the classic algorithms, such as PSO and genetic algorithm (GA), on both the testing functions and real-world problems due to their special mechanisms for avoiding falling into local optimums. Besides, these newly proposed algorithms usually have fewer hyper-parameters, which it is claimed can significantly affect the final performance in a calibration task (Reddy & Kumar 2020), thus, they can have better robustness and usability than the traditional ones. However, the applications of the aforementioned techniques in environmental model calibration are still rare (Reddy & Kumar 2020). With the increasing need for models with high accuracy in refined management, a closer connection between them will be required.

Being one of the representatives of those state-of-the-art algorithms, the MFO is a paradigm inspired by the moth's movement in the night. It combines a population-based algorithm and local strategy for both global exploration and local exploitation, which leads to an increase in its efficiency. Besides, its other main advantages over other algorithms can include simplicity, flexibility, and hybridizability (Shehab et al. 2020). Due to its superiority, it has been utilized to solve a variety of problems in different realms (Jangir et al. 2016), including medical diagnosis (Wang et al. 2017), machine learning (Li et al. 2016), and image processing (Zhou et al. 2018), and it is worth being introduced to the field of environmental modeling. However, since there is a ‘no free lunch’ theorem for optimization (Lynn & Suganthan 2017), every algorithm has its limitations. For the MFO, weaknesses were claimed to be a relaxed convergence (Jain & Saxena 2019) and the risk of being stuck at local optima in some specific applications (Lai et al. 2018). As a result, further enhancement of the original MFO is necessary to fit the potential requirements of the complicated environmental model calibration tasks. One feasible option is the coupling with several algorithm improvement techniques, including chaotic mapping (Gandomi et al. 2013), Lévy flight (Houssein et al. 2020), Opposition-Based Learning (OBL) (Tizhoosh 2005), and Elite-OBL (EOBL) (Tubishat et al. 2019). These improvements can be made at different stages of the optimization and result in hybrid mechanisms. For instance, chaotic mapping can be embedded into the algorithm for random number generations (Saremi et al. 2014), while the EOBL can be employed at the end of each iteration of the algorithm to improve the directionality for better convergence (Ma et al. 2021). Among them, EOBL is a new technology in intelligence computation. Its main idea is to calculate the opposite solution of a feasible solution and replace it if the calculated position is better (Zhou et al. 2016). By doing this, the algorithm can focus more on searching in a feasible area, and therefore both the convergence and the diversity can be improved. Considering the ‘relax convergence’ characteristic of the MFO due to its ‘local exploitation’ behavior after the initial stage, EOBL will be quite fit for it. The crucial thing to be considered when integrating it with MFO is how to place this mechanism appropriately to keep its advantages of ‘local exploitation’ for a better ergodic, and also in case of the occurrence of premature convergence. To do this, a comprehensive understanding of both the MFO and EOBL is required.

Another problem is the conflict between ‘deterministic optimization’ and the so-called parameter ‘equifinality’ as well as different sources of uncertainties in environmental models. Typically, heuristic algorithms are based on the idea that there is a unique global optimum in the search space, while the global optimum of the parameter of an environmental model may be sometimes ambiguous due to the ‘equifinality’. As a result, the prediction based on a specific parameter set would not be robust enough. To make up for this problem, researchers tried to couple statistical methods with the heuristic algorithm to estimate the parameter uncertainty. Blasone et al. (2008) coupled the SCE-Metropolis-UA (SCEM-UA) algorithm with the Generalized Likelihood Uncertainty Estimation (GLUE) framework for model calibration, parameter uncertainty analysis, and prediction uncertainty analysis. Their results suggested that the SCEM-UA can significantly improve the sampling efficiency in the high probability density region (HPDR). The well-known calibration software SWAT-Calibration and Uncertainty Program (SWAT-CUP) combined the PSO with the frequency statistic-based method to quantify the prediction uncertainty. Its performance on uncertainty analysis has been tested and discussed in different types of research areas (Zhang et al. 2015; Liang et al. 2021; Tang et al. 2021). The principle of these methodologies is utilizing the search history of the heuristic algorithm as a kind of sampling technique for the subsequent statistical methods. Since the searching behavior of the heuristic algorithm is affected by the posterior information, sampling using heuristic algorithms can generate statistically sufficient behavioral samples within fewer total samplings. Most importantly, the ensemble modeling based on a set of behavioral samples makes the prediction more robust. When using recently proposed state-of-the-art algorithms for environmental model calibration, the abovementioned conflict still exists, since these algorithms share similar core ideas. Therefore, it is necessary to combine these algorithms with uncertainty quantification techniques.

Inspired by these contexts, this study attempted to fill in the gap between the cutting-edge development of the heuristic algorithm and the application of them in the complex environment models' auto-calibration. Based on our knowledge of MFO and the EOBL theory, a novel EOBL-modified MFO was proposed to fix the limitations of the original algorithm. Meanwhile, the limitation of ‘deterministic optimization’ was also fixed by coupling the algorithm with the frequency statistic-based uncertainty estimation method. The capability of the proposed scheme was examined on the calibration of a Storm Water Management Model (SWMM)-MIKE 11 coupling model, which is a complex hydrologic/hydraulic-water quality coupling model with two software executables. Besides, the global sensitivity analysis and behavioral sampling result analysis were performed to characterize the parameter space of the model and the sampling behavior of the algorithm. The proposed modified algorithm and calibration methodology are expected to be good references for other models' calibration and application, which can include hydrologic models, hydraulic models, deep learning models, and even models in different realms.

For flow generation and hydraulic simulation, the model parameters were pre-calibrated (the coefficient of Nash–Sutcliffe efficiency (NSE) > 0.8 and NSE > 0.7 for calibration and validation, respectively) using the same algorithm EOMFO. The calibration and validation processes will be demonstrated in our other papers. When doing the non-point source pollutant build-up and wash-off simulation, we adopted the exponential build-up model and exponential wash-off model, the governing equations can refer to Tu & Smith (2018). At the stage of the river water quality simulation, the MIKE 11 Advection-Dispersion module uses the 1-dimensional advection-dispersion equation as its governing equation (DHI 2003).

The water quality series from 1 May 2021 to 31 July 2021 was used for calibration, and the series from 1 August 2021 to 1 November 2021 was for validation. During the calibration, 90 moth search agents and 10% EOBL calculations in each iteration were set, and the iteration number was set as 30. As a result the algorithm was designed to minimize the objective function value, 1 – NSE used in this study. The search ranges of the pollutant build-up and wash-off parameters were determined according to the pre-calibration analysis and relevant studies (Chow et al. 2015; Zeng et al. 2019; Taghizadeh et al. 2021), and the value ranges are given in Table 1. Since these parameter values have significant regional variability, the search ranges are not identical to the literature. In addition to the build-up and wash-off parameters, there were still other parameters that needed to be calibrated, including the sweeping efficiency of the roads (or plazas), the pollutant concentration in the rain, and the decay coefficients in pipelines and reaches. The search ranges of those parameters are given in Table 2. The total dimension of the parameter space was 32 for both the NH₃-N and TP simulations.

To comprehensively evaluate the superiority of the proposed EOMFO, PSO and the original MFO were adopted as two benchmarks. For the benchmark algorithms, the search ranges of parameters were set to be identical to EOMFO. The population size and the number of iterations were 100 and 30, respectively. To be noticed, the actual calculations of EOMFO in each iteration were less than 100 due to the shrinking of the flame population. Therefore, the total calculations of EOMFO were less than PSO and MFO, resulting in a more rigorous requirement for the convergence of EOMFO to obtain better results.

The prediction uncertainty estimation in this study was based on a frequency statistic method, which is also the method adopted in SWAT-CUP for considering interval prediction based on behavioral samples and widely used in hydrologic simulations (Liang et al. 2021). By using it, a band of 95% of prediction uncertainty (95PPU) was generated, and as a result, the problem of ‘equifinality for different parameters’ faced by heuristic algorithms was overcome. It was carried out on PSO, MFO, and EOMFO for a crosswise comparison. The performance of the uncertainty band also reflects the diversity and representativeness of the behavioral samples obtained by the algorithm. The specific steps of the prediction uncertainty estimation method include:

• i.

  Select the behavioral samples according to the threshold objective function value, while those non-behavioral samples are discarded.

• ii.

  Sort the simulation results of the behavioral samples at each time step according to the prediction value.

• iii.

  Determine the number of intervals (typically 10) according to the extent of the prediction value ranges at each time step.

• iv.

  Calculate the frequency distribution of prediction values at different intervals.

• v.

  Calculate the 2.5 and 97.5% percentile according to the frequency distribution.

According to the p-values obtained from the global sensitivity analysis (Figure 5(b)), for the NH₃-N simulation, ⁠, ⁠, ⁠, ⁠, ⁠, ⁠, ⁠, ⁠, ⁠, and are parameters that have no significant effect on the simulation performance. While for the TP simulation, the non-sensitive parameters include ⁠, ⁠, ⁠, ⁠, ⁠, ⁠, ⁠, ⁠, ⁠, ⁠, and ⁠. Interestingly, unlike the ⁠, the decay coefficient for pipelines is not sensitive, which is because the components have longer detention time in the reach. In addition, the only accounts for the pollutant concentration during rainfall events, which only endures very short periods and minor effect on the pollutant concentration, while the affects the pollutant concentration in the whole simulation period.

Except for a direct comparison between sampling behaviors and sampling qualities of different algorithms, there is still a detail to be noticed. According to the response surfaces of obtained by all the algorithms for NH₃-N simulation, a clear decreasing trend of the objective function value with the increasing parameter value can be observed. This result suggests that the optimal value of the decay rate may not be covered by the search range. However, the parameter out of the given ranges would not be practical and physically meaningful. Since the performance response surface is under the synthesized impact of the parameter, the calibration data, the model structure, and other unconsidered factors, an excessive parameter range will overestimate the effect of parameters and result in the overfitting problem. Therefore, prior knowledge about the parameter is important for a reasonable calibration result and should be carefully evaluated; otherwise, the algorithm would focus on searching around impractical value ranges to compensate for the prediction errors.

To evaluate the water quality simulation performance of the optimized solution obtained by the tested algorithms, evaluation metrics including NSE, RMSE, and R have been calculated and provided in Table 3. According to the table, both the NH₃-N and TP simulations obtained satisfactory performance with NSE values higher than 0.5 in the calibration stage for EOMFO. In the validation stage, their NSE values are still around 0.5, indicating that the calibrated parameter sets are robust and are not overfitted. The extents of RMSE values of NH₃-N and TP simulations suggest that the optimal parameter sets have moderate but acceptable deviations. Correlation coefficients R over 0.7 in both the calibration stage and validation stage for EOMFO suggest a good reflection of the water quality variation pattern obtained by the optimal solution of the algorithm. In comparison with the benchmark algorithms, the optimal solutions obtained by EOMFO for both the simulated pollutants outperformed those obtained by PSO and MFO regarding almost all evaluation metrics, although the value of R for NH₃-N simulation in calibration was slightly lower. Overall, the performance of deterministic optimization of EOMFO was evidently better than PSO and MFO in both the calibration and validation stage and for both the simulated pollutants. The water quality simulation results support the conclusions obtained from calibration performance and behavioral sampling analysis.

The characteristics of the prediction uncertainty band, as well as the median prediction performance (NSE_50%, RMSE_50%, and R_50%) using behavioral parameters, were also assessed for quantifying the capability of the algorithm on ensemble modeling. All these evaluation metrics are given in Table 4. Since the PSO failed to generate sufficient behavioral samples for the TP simulation, the metrics of its uncertainty band prediction are not available for this pollutant. From the perspective of the uncertainty band prediction, the performances of the EOMFO median prediction result obtained similar but slightly worse NSE, RMSE, and R values than the best parameter sets for both the NH₃-N and TP simulations. For the p-factor, it can be found that the prediction uncertainty band of the NH₃-N simulation did not bracket sufficient observation points in either the calibration stage (6.7%, p-factor = 0.067) or the validation stage (18.0%, p-factor = 0.18). The performance of the uncertainty band was better in the TP simulation, although only 26.4% of the observation points were bracketed by the uncertainty band in the calibration period, this ratio in the validation period was over 80%. For the r-factors, the values indicate the width of the uncertainty bands of both the NH₃-N and TP simulations remained steady throughout the whole simulation. According to a crosswise comparison, the uncertainty band obtained by EOMFO outperformed PSO and MFO from a median prediction perspective with significantly higher NSE_50%, lower RMSE_50%, and higher R_50% values. Making a comparison between the uncertainty band characteristics, MFO and PSO have similar performances considering their coverages (p-factors) and the bandwidths (r-factors), while EOMFO has a minor to moderate advantage according to its higher p-factor and limited increment of r-factor, especially in the validation stages in both the NH₃-N and TP simulation. These results illustrate that the EOMFO not only has superiority in deterministic optimization but also uncertainty band prediction. The potential reason is the higher quality of the behavioral samples discussed in the previous section.

There are several potential reasons for the relatively worse performance of the NH₃-N simulation in this case study. The first of which is the deviation caused by the model structure, including the wash-off calculation or the lack of consideration of the interactions between different nitrogenous components during their fate and transport processes. The wash-off mechanism of the NH₃-N may not exactly follow the exponential pattern at the end of the rainfall event, resulting in overestimations of the value. Additionally, several complicated nitrogenous component transformation processes, for instance, nitrification and denitrification, which are significant in river water conditions (David et al. 2011), can also affect the accuracy of the MIKE 11 AD module that only considers the first-order reaction kinetics and the advection-dispersion mechanism. The second reason is the error of the automatically measured data for relatively low concentrations in dry periods. According to China's national standard for online automatic monitoring stations for NH₃-N (HJ 101-2019), the observation error will be larger when the concentration is lower than 0.15 mg/L. According to our grab sampling data at about 8 km upstream of the observation point in the same year, the NH₃-N concentration determined by the experiment is around 0.2–0.3 mg/L in dry periods, it is significantly higher than the average level of the data at the automatic station (0.02–0.05 mg/L) although the degradation has been considered. This potential error would partially result in the unreasonable response surface of the parameter for NH₃-N simulation mentioned in the previous section, since a higher value can compensate for the deviation caused by it, and it would also further affect the performance of the algorithm on fitting peak values. A direct comparison between the automatic station and the chemical experiment at the same location should be carried out in the future to verify and quantify the potential error. However, since the low concentration value during the dry period is not of great focus, this potential error could have a limited impact on real-world pollution control. Another probable reason is the synthesized effect of the free outfall setting of the coupling model and the characteristic of the NH₃-N discharge. The free outfall may overestimate the pollutant discharged to the channel when the flow rate is not enough to generate the overflow. This effect could frequently occur at the beginning and end of the rainfall event.

The applicability of the proposed algorithm, as well as the presented scheme, is not limited to the specific type of research cases. Except for the mechanism-based water quality model used in this case, data-driven models, such as artificial neural networks (ANNs), support vector machines, and random forests, have also been widely and successfully adopted for different types of water quality modeling tasks (Alizadeh et al. 2018; Shamshirband et al. 2019; Hadjisolomou et al. 2021; Kouadri et al. 2021; AlDahoul et al. 2022). Although these models typically do not have a parameter or a calculation step to represent a particular physical, chemical, or biological process, their parameters or hyper-parameters can still significantly affect their final performance and are needed to be carefully tuned. Trying to make full use of the advantages of heuristic algorithms, some studies started to integrate them into the calibration of data-driven models. For instance, Deng et al. (2021) integrated PSO and GA in the calibration of the ANN model and reported that the best performance occurred when combining the PSO with the Levenberg-Marquardt method. The outperformance of PSO over GA in their results, and also similar effects observed by others (Mirjalili 2015), indicated that the advantage of an algorithm is usually not case-specific and can be reflected in various applications. Therefore, the potential to further integrate the proposed EOMFO with the calibration of data-driven models is profound. On the other hand, most of the methods used in the abovementioned data-driven modeling studies did not consider the equifinality of parameters and uncertainty or only had a primary consideration based on the best model. Interval prediction enabled by a prediction uncertainty estimation based on a set of behavioral samples can avoid the estimation based on an overfitted model and can cover more peak value situations than only using the best parameter set, hence it can make data-driven prediction more robust and reliable for decision-makers.

This study presented a water quality model calibration work using an EOBL-modified MFO coupled with the frequency-based prediction uncertainty estimation method for ensemble modeling. In addition, the global sensitivity analysis and behavioral sampling result analysis were performed to characterize the properties of the parameter space as well as the sampling behavior. The main findings of this article include:

• i.

  Build-up parameters and wash-off coefficient of the urban area (roads and plazas) and roofs, as well as the decay coefficient in reach, are sensitive parameters for both the NH₃-N and TP simulation. Among them, the build-up rate of the urban area is the most sensitive parameter due to the continuous wet season simulation.

• ii.

  EOMFO has a better convergence compared with PSO and MFO from the perspective of global optimization. It can also obtain behavioral samples with higher quality near the true response surface of the parameter space, although the behavioral sampling rate at the early stage of the calibration process is still relatively slow.

• iii.

  EOMFO outperformed the benchmarks in both the NH₃-N simulation and TP simulation. The coupled uncertainty estimation method can generate robust uncertainty bands for ensemble modeling, although its performance for NH₃-N simulation was probably hindered by several potential factors.

Future works can focus on improving the limitations of the calibration scheme and the constructed model. To improve the slow behavioral sampling rate at the early stage, hybridization with other algorithms can be considered. Besides, the algorithm has a very high ratio of behavioral samples. This may sometimes result in an oversized behavioral sample set, increasing the time needed for validation and ensemble modeling. To resolve this problem, a bootstrap resampling method for shrinking the sample set can be adopted. Alternatively, stop criteria can be introduced to terminate the algorithm when the expected convergence has been achieved. For the constructed model, more details including the CSO system, the outfall type, and the fate and transport of nitrogenous components can be considered when the supporting data is sufficient. Additionally, the property of the water quality data at automatic stations can be further investigated to better support the relevant applications.

Thanks for the water quality monitoring data provided by the Environmental Knowledge Service System developed by the Chinese Research Academy of Environmental Sciences.

Data cannot be made publicly available; readers should contact the corresponding author for details.

The authors declare there is no conflict.