Generalized Likelihood Uncertainty Estimation Method in Uncertainty Analysis of Numerical Eutrophication Models : Take BLOOM as an Example

Uncertainty analysis is of great importance to assess and quantify a model’s reliability, which can improve decision making based on model results. Eutrophication and algal bloom are nowadays serious problems occurring on a worldwide scale. Numerical models offer an effective way to algal bloom prediction and management. Due to the complex processes of aquatic ecosystem, such numerical models usually contain a large number of parameters, which may lead to important uncertainty in the model results. This research investigates the applicability of generalized likelihood uncertainty estimation (GLUE) to analyze the uncertainty of numerical eutrophication models that have a large number of intercorrelated parameters. The 3-dimensional primary production model BLOOM, which has been broadly used in algal bloom simulations for both fresh and coastal waters, is used.


Introduction
Eutrophication and algal bloom are serious problems occurring on a worldwide scale, which deteriorate the water qualities in many aspects, including oxygen depletion, bad smell, and production of scums and toxins.Accurate and reliable predictions of algal blooms are essential for early warning and risk mitigating.Numerical eutrophication models offer an effective way to algal blooms prediction and management.There exist several well-developed eutrophication models, such as CE-QUAL-ICM [1,2], EUTRO5 [3][4][5], BLOOM [6][7][8][9][10], CAEDYM [11,12], and Pamolare [13,14].The choice of the most appropriate model may depend on the specific research objectives and data availability.
Due to the complexity of algal bloom processes, these numerical models usually have a large number of parameters, which inevitably brings uncertainty to model results.Modeling practice typically includes model development, calibration, validation, and application, while uncertainty analysis is often neglected.Uncertainty analysis is essential in the assessment and quantification of the reliability of models.Prior to the use of model results, information about model accuracy and confidence levels should be provided to guarantee that results are in accordance with measurements [15] and that the model is appropriate for its prospective application [16].There are three major sources of uncertainty in modeling systems: parameter estimation, input data, and model structure [17][18][19][20][21]. Understanding and evaluating these various sources of uncertainty in eutrophication models are of importance for algal bloom management and aquatic ecosystem restoration.
Several methods for parameter uncertainty analysis are available, for example, probability theory method, Monte Carlo analysis, Bayesian method, and generalized likelihood uncertainty estimation (GLUE) method.Probability theory method employs probability theory of moments of linear combinations of random variables to define means and variances of random functions.It is straightforward for simple linear models, while it does not apply to nonlinear systems [22].The Monte Carlo analysis computes output statistics by repeating simulations with randomly sampled input variables complying with probability density functions.It is easily implemented and generally applicable, but the results gained from Monte Carlo analysis are not in an analytical form and the joint distributions of correlated variables are often unknown or difficult to derive.[18,20,23].Bayesian methods quantify uncertainty by calculating probabilistic predictions.Determining the prior probability distribution of model parameters is the key step in Bayesian methods [24].GLUE is a statistical method for simultaneously calibrating the input parameters and estimating the uncertainty of predictive models [25].GLUE is based on the concept of equifinality, which means that different sets of input parameter may result in equally good and acceptable model outputs for a chosen model [26].It searches for parameter sets that would give reliable simulations for a range of model inputs instead of searching for an optimum parameter set that would give the best simulation results [27].Furthermore, model performance in GLUE is mainly dependent on parameter sets rather than individual parameters, whence interaction between parameters is implicitly accounted for.
Beven and Freer [28] pointed out that in complex dynamic models that contain a large number of highly intercorrelated parameters, many different combinations of parameters can give equivalently accurate predictions.In consideration of equivalence of parameter sets, the GLUE method is particularly appropriate for the uncertainty assessment of numerical eutrophication models, which are an example of complex dynamic models with highly intercorrelated parameters [28][29][30].
The GLUE method has already been adopted for uncertainty assessments in a variety of environmental modeling applications, including rainfall-runoff models [25,31,32], soil carbon models of forest ecosystems [33], agricultural nonpoint source (NPS) pollution models [34], groundwater flow models [35], urban stormwater quality models [36,37], crop growth models [38], and wheat canopy models [39].The popularity of the GLUE method can be attributed to its simplicity and wide applicability, especially when dealing with nonlinear and nonmonotonic ecological dynamic models.
The objectives of this paper are to make use of the broadly used eutrophication and algal bloom model BLOOM [6,7], in order to investigate the applicability of the GLUE method to analyze and quantify the uncertainty in numerical eutrophication models that have a large number of intercorrelated parameters and to provide a reference for method selection when conducting uncertainty analysis for similar types of models.

Materials and Methodology
2.1.Study Area.The Meiliang Bay (31 ∘ 27  N/120 ∘ 10  E), which locates at the north of Taihu Lake in China (Figure 1), is chosen as the study area.Taihu Lake has high level of eutrophication, and algal blooms that cause enormous damage to drinking water safety, tourisms, and fish farming frequently break out in summer and autumn.
The Meiliang Bay has a length of 16.6 km from south to north, a width of 10 km from east to west, and an average depth of 1.95 m.There are two main rivers: the Zhihu Gang that flows into Taihu Lake and the Liangxi River that flows out of the lake.The exchange of substance between Meiliang Bay and the main body of Taihu Lake is taken into account in the study.The monthly observed data from four monitoring sites in the Meiliang Bay were collected during 2009 to 2011 for model calibration.These data include river discharge, water level, irradiance, temperature, concentrations of ammonia, nitrate, nitrite, phosphate, and biomass concentration of bluegreen algae, green algae, and diatom.
In the composition of the algae blooms, blue-green algae is the dominant species and has the highest percentage of total biomass.Therefore it is selected to be the output variable of BLOOM on which the uncertainty analysis is performed.

BLOOM Model.
BLOOM is a generic hydroenvironmental numerical model that can be applied to calculate primary production, chlorophyll-a concentration, and phytoplankton species composition [6][7][8][9][10].Fifteen algae species can be modeled, including blue-green algae, green algae, and diatoms.Each algae species has up to three types, the N-type, P-type, and E-type, which correspond to nitrogen limiting conditions, phosphorus limiting conditions, and energy limiting conditions, respectively.Algae biomass in BLOOM mainly depends on primary production and transport.
The transport of dissolved or suspended matter in the water body is modeled by solving the advection-diffusion equation numerically: where : concentration (kg⋅m −3 );   ,   : dispersion coefficient in and -direction respectively (m 2 ⋅s −1 ); S: source terms;  (, ): reaction terms; : time (s).
Primary production is mainly dependent on the specific rates of growth, mortality, and maintenance respiration, which are modulated according to the temperature: where   : potential specific growth rate of the fastest growing type of algae species (d −1 ); Proalg 0  : growth rate at 0 ∘ C (d −1 ); TcPalg  : temperature coefficient for growth (-);   : specific mortality rate (d −1 ); Moralg 0  : mortality rate at 0 ∘ C (d −1 ); TcMalg  : temperature coefficient for mortality (-);   : specific maintenance respiration rate (d −1 ); Resalg 0  : maintenance respiration rate at 0 ∘ C (d −1 ); TcRalg  : temperature coefficient for respiration (-).More details of BLOOM can be found in Delft Hydraulics [6,7].[25] is based upon a large number of model runs performed with different sets of input parameter, sampled randomly from prior specified parameter distributions.The simulation result corresponding to each parameter set is evaluated by means of its likelihood value, which quantifies how well the model output conforms to the observed values.The higher the likelihood value, the better the correspondence between the model simulation and observations.Simulations with a likelihood value larger than a user-defined acceptability threshold will be retained to determine the uncertainty bounds of the model outputs [33,40].The major procedures for performing GLUE include determining the ranges and prior distributions of input parameters, generating random parameter sets, defining the generalized likelihood function, defining threshold value for behavioral parameter sets, and calculating the model output cumulative distribution function.

Generalized Likelihood Uncertainty Estimation. The GLUE methodology
BLOOM contains hundreds of parameters.Ideally, all the parameters should be regarded stochastically and included in the uncertainty analysis.However, a more practical and typical manner to conduct uncertainty analysis is to focus on a few key parameters [25,28,41].In eutrophication models, algal biomass is most closely related to the growth, mortality, and respiration processes, resulting in the selection of seven key parameters about blue-green algae according to (2).Table 1 summarizes the main characteristics of these seven parameters.The initial ranges of the parameters are obtained by model calibration and a uniform prior distribution reported in the literature [28] is considered for all parameters.
Latin Hypercube Sampling (LHS), which is a type of stratified Monte Carlo sampling, is employed in this study to generate random parameter sets from the prior parameter distributions.In total, 60,000 parameter combinations are generated for the model runs.
The GLUE method requires the definition of a likelihood function in order to quantify how well simulation results conform to observed data.The likelihood measure should increase monotonically with increasing conformity between simulation results and observations [25].Various likelihood functions have been proposed and evaluated in the literature [35,37,38,42].Keesman and van Straten [43] defined the likelihood function based on the maximum absolute residual; Beven and Binley [25] defined the likelihood function based on the inverse error variance with a shape factor ; Romanowicz et al. [44] defined the likelihood function based on an autocorrelated Gaussian error model; Freer et al. [45] defined the likelihood function based on the Nash-Sutliffe efficiency criterion with shape factor , as well as the exponential transformation of the error variance with shaping factor ; Wang et al. [41] defined the likelihood function based on minimum mean square error.In this study, the likelihood function (  | ) of the model run corresponding to the th set of input parameters (  ) and observations  is defined based on the exponential transformation of the error variance  2  and the observation variance  2 0 with shape factor  [37,45]: where  2  = ∑ ( sim −  obs ) 2 ;  2  = ∑( obs −  obs ) 2 ;  sim is the simulated blue-green algae biomass;  obs is the observed blue-green algae biomass;  obs is the average value of  obs .
The sensitivity of the choice of the shape factor  will be analyzed and discussed.If the likelihood value of a simulation result is larger than a user-defined threshold, the model simulation is considered "behavioral" and retained for the subsequent analysis.Otherwise, the model simulation is considered "nonbehavioral, " and removed from further analysis.There are two main methods for defining the threshold value for behavioral parameter sets: one is to allow a certain deviation from the highest likelihood value in the sample, and the other is to use a fixed percentage of the total number of simulations [46].The latter is used in this study, and the acceptable sample rate (ASR) is defined as 60%.The sensitivity of the choice of the threshold in the form of the acceptable sample rate (ASR) will be analyzed and discussed.
The likelihood function is then normalized, such that the cumulative likelihood of all model runs equals 1: where   (  ) is the normalized likelihood for theth set of input parameters (  ).The uncertainty analysis is performed by calculating the cumulative distribution function (CDF) of the normalized likelihood together with prediction quantiles.The GLUE-derived 90% confidence intervals for the biomass of blue green are then obtained by reading 5% and 95% percentiles of the cumulative distribution functions.

BLOOM Model Results.
The calibration result for bluegreen algae is shown in Figure 2, and the calibration  2. The mean values of the bluegreen algae biomass for the four sample sites are similar.Therefore, in order to reduce the sampling uncertainties, the average of the four sampling sites has been retained as dependent variable in the present study.The biomass of blue-green has a yearly cycle (Figure 2), with low values during spring, followed by a rapid increase towards peak values in summer or autumn.The growth periodicity of blue-green algae is mainly attributed to the periodic variation of temperature and algae dormancy.Taihu Lake experiences a subtropical monsoon climate, with four distinct seasons.The lowest temperature is about 2.8 ∘ C in average and appears in January, and the highest temperature is about 29.4 ∘ C in average and usually appears in August.The suitable temperature range for growth of blue-green algae is 25∼35 ∘ C. As a result, the biomass of blue-green algae is low in spring.When temperature increases in summer, it is appropriate for blue-green algae breeding, leading to the sharp increase in biomass and the occurrence of the peak value around August.
The modes capture satisfactorily the observed evolution of the blue-green algae biomass, which indicates further analyses on model uncertainty are meaningful.The coefficient of determination (CoD), which is given by ( 5), is 0.85: where   : the simulated biomass of blue-green algae at time step ;   : the mean value of observed data;   : the observed value of blue-green at time step .

Uncertainty Analysis Results
. The confidence interval (CI) is obtained by calculating the cumulative distribution functions of model outputs based on the normalized likelihood (4) with  = 1 and ASR = 60%.Figure 3 presents the 90% confidence interval of blue-green algae biomass, which is estimated from the 5% and 95% quantiles of the cumulative distribution functions and the corresponding observations from January 2009 to December 2011.Table 3 summarizes the width of the 90% CI of each month and whether or not the observations are located within the 90% CI.
The 90% CI is narrow from January to May when the biomass of blue-green algae is low.The width of the 90% CI expands as the biomass of blue-green algae increases during summer and autumn.Among the total of 36 observations, 13 are located within the 90% CI of the simulations.
The subjective choice of the shape factor  in (1) considerably influences the GLUE results, whereas  is commonly taken as 1 [35].Figure 4 displays the 90% CI when ASR = 60%, with shape factors  equal 50 and 100, respectively.The simulated 5% and 95% confidence quantiles and the weighted mean, as well as the corresponding observations of bluegreen algae biomass, are shown.
Comparison of Figures 3 and 4 shows that the increase of shape factor  leads to a narrowing of the 90% CI. Figure 5  illustrates the effect of the shape factor , which can be seen as a weight factor for the likelihood corresponding to each simulation.When  = 1, the magnitudes of the likelihood are similar for each simulation, and there is no clear division between acceptable and unacceptable simulations.As a result, the cumulative distribution functions increase gradually.With increasing  (e.g.,  = 50), the high behavioral simulations have a higher weight, resulting in a larger gradient in the cumulative distribution function and a narrower CI.Theoretically, when  = 0, every simulation has equal likelihood, and the widest CI will be obtained.When  → ∞, the single best simulation will have a normalized likelihood of 1, while all other simulations will get a likelihood of zero, resulting in the collapse of the 5% and 95% quantiles on a single line.This corresponds to the traditional calibration method that omits uncertainty analysis.
Previous studies have shown that the choice of threshold values for the likelihood measures is particularly important for the GLUE method [34,36,47].In order to quantify the effect of threshold values on the uncertainty analyses, a series of acceptable sample rates (ASR) of 0.5%, 1%, 5%, 10%, 30%, 60%, 90%, 95%, 99% is investigated.In this study, average relative interval length (ARIL) and percentage of observations covered by the 90% confidence interval ( 90CI ) are adopted as metrics for the analysis.These metrics are defined as follows: where Limit upper, and Limit lower, are the upper and lower boundary values of the 90% confidence interval;  is the number of time steps;  obs, is the observed biomass of bluegreen algae: where  in is the number of observations located within 90% CI;  obs is the total number of observations.Figures 6 and 7 present the influence of ASR on ARIL and  90CI for  = 1, 50, 100. Figure 6 shows that, for all ASR values, ARIL has the highest value for  = 1 and decreases with increasing , which confirms the results of Figure 4.For a given  value, ARIL increases with ASR.When ASR moves from 0.5% to 99%, the ARIL increases by 73.93%, 41.96% and 5.24% for  = 1, 50, 100, respectively.An increasing ASR, which corresponds to a lower threshold of the accepted likelihood, means that simulations with lower likelihood are considered "behavioral, " which inevitably results in a larger ARIL.
From Figure 7, it is seen that  90CI becomes larger as ASR increases for  = 1 and  = 50, while  90CI keeps constant for  = 100.This is because the increase of ASR results in a larger ARIL, which logically leads to an increase in observations located within the 90% CI.When  = 100, the ARIL is low and  90CI does not increase with ASR because the 90% CI does not widen.
The highest  90CI is obtained for ASR close to 100% and  = 1.Its value of about 50% indicates that about half of the observed data remain outside the 90% CI for the greatest ASR.This can be attributed to other sources of uncertainty, such as the input parameters or the observations.

Discussion and Conclusion
The 90% confidence interval of the simulated results fails to enclose the peaks of the observed values in 2009 and 2011 (Figure 3).Such a feature is not unusual, and several reasons can lead to this result.Firstly, there are inherent uncertainties from inputs, boundaries, and model structure, which are not    taken into account explicitly.Secondly, the observed values used for comparison are space averaged through arithmetic mean other than weighted mean, which could also introduce discrepancy.Finally, the original observations from the four stations contain measurement uncertainties.During summer and autumn when the algal blooms break out, the biomass of blue-green algae is high and shows pronounced daily temporal variations and spatial variations due to changes in irradiance, transport by flow and wind drifting.The measurements taken at a particular time and point cannot fully reflect these fine-scale spatial and temporal dynamics.
The model calibrated in this study is, however, capable of simulating blue-green algae dynamics at large spatial (spatial averages) and temporal (seasonal) scales.Ideally, accurate predictions require that the results are consistent with the observations, while the uncertainty spread of the results, quantified by the 90% CI, is as narrow as possible [46].From Figures 6 and 7, it can be seen that while keeping ASR fixed, the 90% CI is narrowed by increasing the shape factor  at the expense of decreasing the percentage of observations that it covers ( 90CI ).Similarly, while keeping  fixed, the 90% CI is narrowed by reducing ASR, but at the same time also  90CI decreases.As a consequence, it is essential to optimally choose  and ASR in order to find the optimal compromise between the uncertainty spread and its coverage of observations.
As illustrated by the application to the BLOOM model for algal bloom, GLUE is an appropriate method for uncertainty analysis that can cope with equifinality between different parameter sets incurred by high level of model complexity.
In conclusion, the study demonstrates that GLUE is an effective method for uncertainty analysis of complex dynamic ecosystem models, which provides a solid foundation for the use of the model predictions in decision making.

Table 1 :
Selected input parameters and their initial ranges.
parameters are summarized in the last column of Table1.The statistical characteristics of the observed blue-green algae biomass are shown in Table

Table 3 :
Width of 90% confidence interval (CI) for each month and indication whether (Y) or not (N) the observations is located within the 90% CI band.