Setting up a water quality ensemble forecast for coastal ecosystems: a case study of the southern North Sea

Water quality is a crucial factor for both coastal ecosystems and human societies. Phytoplankton blooms in the North Sea may cause mortality of mussels and other benthic organisms. Furthermore, fisheries and aquacultures are influenced by the algal primary production since it is the base of the food web. Consequently, accurate real-time phytoplankton concentration prediction is required for ecosystem and economic benefits. Timely information regarding water quality allows for early warning and adequate response such as mitigation measures and targeted monitoring.

Existing hybrid ecosystem models are powerful tools for modelling water quality; however, their reliability depends highly on the uncertainty stemming from different sources. Uncertainty originating from external forcings is further propagated and complicated by the non-linear ecological processes, with numerous intercorrelated parameters incorporated in the water quality model. Considering the high level of uncertainty in the coastal water quality forecasting process, it is assumed that a single-valued deterministic forecast may not be sufficiently reliable for decision making. Thus, a strong need arises for an ensemble forecasting system that could potentially account for the uncertainty associated with the driving forces, the model simplifications or the parameterization.

Ensemble forecasting of water quality might require considerably more effort than a deterministic approach, though in return gives added value to the forecast. Deterministic forecasts only provide a point estimate, whereas probabilistic forecasts determine the probability density function (PDF) of the predictand. By performing a statistical analysis of the acquired PDF, various measures can be derived such as the mean and standard deviation of the model output distribution, or the likelihood of a specific output value. Through these measures the uncertainty in the model output can be quantified, which serves decision support functions.

The application of ensemble prediction systems is well known in various fields, especially in numerical weather prediction, yet few examples (Fiechter 2012; Trolle et al. 2014; Huang & Gao 2017) show their applicability in water quality forecasting. A recent study by Ahn et al. (2016) investigated the possibility of developing a water quality forecasting system for riverine ecosystems with an ensemble streamflow prediction method. In that study a historical rainfall and temperature ensemble was applied as a forcing condition in order to produce a probabilistic water quality prediction. Based on those findings it can be concluded that implementing an ensemble system for water quality forecasts by addressing the uncertainty in the input variables is advantageous. Another research in ensemble forecast accuracy in the southern North Sea was conducted by Herrmann (2015). The ensemble forecasting system was set up using perturbed meteorological forcings and the Delft3D-FLOW model. Despite expectations, the results showed that the impact of the meteorological ensemble input is not sufficient to estimate the total uncertainty in the forecasted parameters (salinity and temperature). In order to complement those findings this paper suggests the application of an extended input ensemble with the aim to provide further uncertainty estimation.

This paper aims to set up the framework for an operational water quality ensemble forecasting system, even though the proposed method is only applied in a hindcasting case study and further steps are required to reach the operational stage. The model output of interest is the Chlorophyll-a concentration which is a commonly used indicator for water quality. The simulation is carried out using the Delft3D Generic Ecological Model (GEM) with the advanced algal speciation module-BLOOM which is part of the Delft3D-WAQ software package.

The paper starts with a description of the GEM, its application to the North Sea, and the different data sources used for the forecast verification. The research methodology is presented in the following sections including the significant parameter selection, the Latin Hypercube Sampling with Dependence (LHSD) steps, and the applied forecast verification method. The performance analysis of the ensemble forecasting system applied to the southern North Sea is given together with the validation and spatial results. The conclusions and recommendations for further research are presented in the final section of the paper.

The modelling instrument is a comprehensive hybrid ecological model combining a three-dimensional hydrodynamic model (Delft3D-FLOW) and the GEM biogeochemical model. The GEM includes an array of modules reproducing water quality processes that are then combined with the transport model. Most importantly, the model computes primary production and Chlorophyll-a concentration while integrating dynamic process modules for dissolved oxygen and nutrient concentration calculation. Furthermore, the GEM includes a phytoplankton module (BLOOM) that simulates the growth, respiration and mortality of phytoplankton. Using this module, the species competition and their adaptation to limiting nutrients or light can be simulated. The model offers flexibility in the processes selection and provides general applicability in diverse case studies (Blauw et al. 2009). In recent years the GEM has already been applied to the southern North Sea. A three-dimensional GEM application was presented by Los et al. (2008), while a two-dimensional application was completed by Salacinska et al. (2010).

Model calibration and validation of the three-dimensional GEM for the southern North Sea was carried out in previous studies by Los et al. (2008) using an in-situ dataset for the year 1989, in Keetels et al. (2012) for the year 2007, and in Arentz et al. (2012) for years 2009–2010 also using the in-situ dataset. In general, all the reports mention acceptable Chlorophyll-a prediction accuracy, however, the results are not homogeneous across the stations and are highly dependent on proper description of the suspended particulate matter (SPM) field. Furthermore, it should be noted that in the Dutch coastal zone, the observed gradients of algal biomass are very steep and there is considerable natural variability in the Chlorophyll-a concentration. The peak spring bloom could be shifted ±1 month and its magnitude could vary ±80% (Keetels et al. 2012).

• is concentration of the state variables

• ⁠, v, w are velocity vector components

• ⁠, ⁠, dispersion tensor components

• ⁠, y, z are coordinates

• is the source and sink term of mass due to loads and boundaries

• is the source and sink term of mass due to processes

• is time ⁠.

The concentrations of 30 state variables including algae concentrations, nutrients, and salinity are calculated through various ecological processes. The most important processes are related to nutrient cycles, oxygen dynamics, energy availability and phytoplankton processes. The GEM ecological processes are considered to be moderately complex by Los et al. (2008) based on the fact that the GEM excludes microbial loop, explicit grazing and higher trophic levels; also the benthic processes are relatively simple. On the other hand, full nutrient cycles and complex phytoplankton kinetics were implemented. Figure 1 contains the schematic overview of all possible variables and processes in the GEM. The reader may refer to Blauw et al. (2009) and Los et al. (2008) for further description of the ecological processes.

The GEM makes use of more than 400 parameters for the process calculations. These parameters are related to algae's characteristics such as nutrient-to-carbon ratios and growth rates. Some describe light availability through extinction coefficients and settling velocities and others are connected to the nutrient cycles. The GEM is part of an integrated modelling system which contains separate modules for hydrodynamics, waves and sediment transport calculation. The water quality simulation therefore requires external forcings such as meteorological conditions, hydrodynamics, SPM concentration field, as well as nutrient loadings from atmospheric deposition and riverine loads.

The domain decomposition of the southern North Sea is a three-dimensional curvilinear grid with finer resolution along the coast as shown in Figure 2. The grid contains 12 sigma layers with unequal thicknesses (see Table 1). The distribution of the layers was designed to provide higher resolution on the top and the bottom of the water column to enable more accurate suspended matter calculation for light attenuation on the top and more detailed resuspension calculation near the bed. The resolution of the coarse segments varies from 6-by-5 km to 20-by-30 km and the resolution of the fine grid ranges from 1-by-2 km to 2.5-by-3 km.

In order to validate the water quality model in this paper three types of data sources are used: a set of historical surface in-situ measurements, and remote sensing images together with their gap-filled reconstruction (see Figure 3). Fifteen stations were selected with available Chlorophyll-a in-situ measurements for the given years 2007 and 2009 in the focus area provided by Rijkswaterstaat (RWS, Dutch Ministry of Infrastructure and Environment). It should be noted that the samples are taken close to the water surface (usually the upper 3–5 metres of the water column) and therefore they can only be used to validate the first layer of the model.

The second type of dataset is remote sensing which gives a more comprehensive view of the model's success in capturing seasonal variability since the spatial and temporal resolution of the in-situ dataset is low. The remotely sensed sea colour images applied in this study are retrieved from the MEdium Resolution Imaging Spectrometer instrument (MERIS), which was on board the European Space Agency ENVironmental SATellite (ENVISAT) spacecraft. Since the ENVISAT ended its mission in 2012, the application of MERIS data in this paper is considered as a preparation for the new Copernicus Sentinel dataset usage which has recently stepped into the fully operational stage.

The recorded optical reflectance of the water surface is transformed into Chlorophyll-a or SPM concentration using the HYDROPT algorithm given in van der Woerd & Pasterkamp (2008). The algorithm not only computes the concentration of the substances but also provides a measure of error in the estimates (El Serafy et al. 2007). The retrieved pixel data are then interpolated onto the southern North Sea Domain Decomposition (ZUNO-DD) grid taking into account the provided error estimate, and resulting in the so-called gridded data (Blaas 2013). It should be noted that the raw pixel data were first interpolated using nearest-to-the-centre interpolation onto the ZUNO-COARSE gridding and then it was subsequently transposed onto the ZUNO-DD domain gridding.

Most forecast verification metrics are based on the differences between the prediction and observation; however, these differences can only be calculated if both values are available for the same time-step at the same location. Considering the incomplete temporal and spatial coverage of the samples from the MERIS measurements due to clouds, the number of matchups is fairly limited. As an attempt to tackle this issue, a third dataset is used for the model validation, which is the gap-filled version of the gridded MERIS images. The gap filling is done using a data interpolating algorithm called DINEOF. The DINEOF algorithm determines the major spatial and temporal patterns of variation in the MERIS dataset and produces the gap-filled reconstruction at all segments and all time-steps (Blaas 2013). The gap-filled data applied in this study was constructed using the MERIS data only from the years of interest: 2007 and 2009.

The main objective of the paper is to set up a viable and generally applicable water quality ensemble forecasting system for coastal ecosystems. The proposed methodology (see Figure 4) includes the following important steps:

• Selection of the significant model process parameters and external forcings responsible for the uncertainty in the Chlorophyll-a concentration prediction.

• Statistically efficient sampling of the input ensemble (stratified sample with parameter dependency) using the selected significant parameters and forcings.

• Post-processing of the ensemble model run by fitting PDFs to the model outputs at every time-step and all locations.

• Forecast verification to assess and compare the performance of the deterministic and ensemble predictions.

The previously introduced ecological processes of the GEM make use of more than 400 model process parameters in total. Incorporating all these parameters into the ensemble generation is not feasible and would not yield a better uncertainty estimate, since most of them do not affect the Chlorophyll-a concentration.

A previous study investigated the sensitivity of the Chlorophyll-a concentration to process parameters in the North Sea application of the GEM (Salacinska et al. 2010) and identified 20 parameters that obtained the highest rank. These parameters consist of growth rates, extinction coefficients and nutrient-to-carbon ratios (e.g. P:C, N:C) of the different algal species. Furthermore, according to Li et al. (2013) primary production mainly depends on the specific rates of growth, mortality, and maintenance respiration as well as the temperature coefficient because these rates are calculated as a function of the temperature. In addition, Daggers (2013) conducted a sensitivity analysis of the nutrient concentrations influenced by the different processes such as denitrification, nitrification, mineralization and burial. Those findings suggest that in coastal regions nitrate is mainly sensitive to the denitrification rate in the top sediment layer (S1), while phosphate and all other nutrients are mainly sensitive to the burial rate. This extensive list of significant parameters was narrowed down to six parameters (see Table 2) considering the dominating algal type (Diatoms type E) and the coastal environment in the southern North Sea. Further description of the significant parameters can be found below.

The algal biomass concentration is calculated using an optimization technique (linear programming) that solves a set of linear equations and constraints. The applied constraints are the energy-, growth-, nutrient- and mortality constraints. This paper focuses on the peak algae bloom prediction which is not limited by the mortality constraint which for this reason is not considered. Therefore, the selected model process parameters can be divided into three main groups depending on whether they affect the energy-, growth- or nutrient constraints as shown in Table 2.

In addition to the model process parameters, certain model forcings are selected to be perturbed based on previous research findings focusing on the southern North Sea coastal ecosystem.

The importance of the concentration of suspended inorganic matter forcing, often referenced as SPM, was previously demonstrated in Equation (4). In the continental coastal waters, 25 to 75% of the light extinction is caused by SPM (Los & Blaas 2010). Furthermore, Perez (2015) investigated the effect of using different types of SPM sources on the Chlorophyll-a concentration in the southern North Sea and the results suggested considerable impact. El Serafy et al. (2011) also indicated that the SPM data sources in the North Sea, such as the field measurements and the Delft3D-WAQ-SPM model results, are rather uncertain. Based on these findings, it was decided to include the SPM concentration field in the ensemble generation. By reason of simplicity, the spatial correlation in the SPM concentration field is not considered, and only a simple error range was applied to the input segment function that contains the SPM concentration for each segment at all time-steps. Spatial correlation, however, was previously found to be important in Kamel et al. (2014) since the seasonal variation in fine sediment dynamics forms several ellipsoidal-shaped sedimentation traps along the Dutch coast.

Further model forcing affecting the Chlorophyll-a concentration, especially in coastal ecosystems, is the riverine nutrient load. Nutrient loads from rivers provide a significant portion of the total nitrogen load (12–17%) and total phosphorus load (8–11%) in the North Sea (Daggers 2013). Since the information about the river discharges and nutrient loads is often estimated, and the daily discharge values are interpolated from the less frequently available flow data, the riverine nutrient input is also included in the ensemble generation. For practical reasons, instead of directly perturbing the nutrient concentrations, the river discharges are altered, which linearly affects the nutrient concentrations. In this paper only those nine river discharges are perturbed that are identified to be influential in the study area based on the nutrient composition matrix derived by Los et al. (2014).

A simplified sensitivity analysis was carried out in order to comprehend the Chlorophyll-a concentration's response to the parameter value modification and to verify if previous findings can be confirmed in the present case. The sensitivity analysis was a one-factor-at-a-time method. This method allows us to gain information about the model's sensitivity to a specific parameter; however, it does not account for the dependencies between variables since the parameters interact in a non-linear way.

In the sensitivity test, model runs were executed using the baseline, lower-, and upper bound values of all identified parameters. The changes in residuals between the baseline and the scenarios were analysed at one station over one year with particular attention given to the peak concentration and timing. A brief summary of the observed correlation patterns between the selected parameters and the Chlorophyll-a concentration can be found in Table 3. The table demonstrates the type and strength of the correlation between the investigated parameters and the output parameter, also it indicates whether the time of the peak concentration was affected or not. While the results confirm previous findings, it should be noted that the denitrification rates were found to be less significant compared to other parameters, hence, these might be negligible in deeper coastal waters.

After selecting the significant parameters, it is crucial to identify the parameter statistics, such as the possible value range, mean and standard deviations, and to specify the desired probability distributions to obtain through the sampling. The parameter statistics are given in Table 4. The mean values are assumed to be the calibrated baseline values of the deterministic model setup. The parameter ranges for process parameters were mainly found in Salacinska et al. (2010) and Blauw et al. (2009), whereas for the model inputs value ranges are set according to the estimated uncertainty level. The uncertainty level in the river discharges was assumed as 25% based on a study conducted by Lloyd et al. (2016) on discharge and nutrient uncertainty in streams in the UK, which reported observational discharge uncertainties ranging from ±2 to 25%. The standard deviations for four process parameters are used as described in Jiayuan (2015). For the remaining process parameters, and for the model inputs, the standard deviations are assumed in such a way that the resulting sampled values fall between the specified value ranges, given that for normal distribution 99.7% of the data are within three standard deviations of the mean. Information about the parameter dependency was available for four process parameters from Jiayuan (2015), while having no information about the rest of the parameters they are assumed to be independent.

In this case study, the selected process parameters are non-measurable and therefore no prior knowledge was available to obtain empirical distributions. Consequently, the theoretical distributions were assumed according to literature (Jiayuan 2015) and general considerations. Based on the literature, exponential distributions were assumed for those parameters which can only take positive values, whereas normal distributions were applied to others. This approach was adopted as best available knowledge on the parameters' behaviour. Moreover, these distributions were adopted partly due to the fact that in the applied sampling procedure a Gaussian copula model is used, which has the underlying assumption that the joint distributions of the random variables are elliptically symmetric. On the other hand, it is important to note that in future cases, when prior information about the parameters is insufficient, choosing uniform prior distributions might be more appropriate as demonstrated by Freni & Mannina (2010). Applying normal distribution without prior knowledge may result in wrong estimation of uncertainty as it leads to a narrower uncertainty band as compared to uniform distribution (Freni & Mannina 2010; Dotto et al. 2012). This may create excessive confidence in the model results especially in water quality modelling applications where the inherent uncertainties are high.

The generated sample has two specific requirements to fulfil: (1) the dependence structure between the variables must be represented; and (2) a good coverage of the parameter space should be achieved with only a limited sample number (variance reduction). The water quality model makes use of numerous model parameters which are in many cases correlated. If the sampling technique ignores those correlations the simulation could result in unrealistic outputs. One way to describe dependency between ecological process variables is to use copulas (Jiayuan 2015). In addition to the parameter dependency, choosing a sampling technique with variance reduction gives a greater precision of the output random variables for a given number of iterations, this way providing a statistically efficient sampling.

The parameters are transformed to the above-mentioned prior distributions. The sample number n = 160 was based on Jiayuan (2015) who conducted a simplified analysis to identify the optimal sample number. Applying the LHSD, smooth and evenly distributed sample distributions could be achieved (see Figure 5) while maintaining the parameter dependency. Samples that fall outside of the predefined range were removed in order to avoid values with non-physical meaning.

The results of the ensemble simulation runs are transformed into probabilities by fitting PDFs to the ensemble members at each time-step. In the applied distribution fitting method, continuous probability distributions were fitted and ranked in order to obtain the best fit. This fitting procedure does not allow multimodal fitting meaning that the PDF cannot have multiple peaks; this simplification might result in inaccuracies.

With the help of multiple verification metrics, different forecast attributes can be assessed and quantified. In this paper the Mean Absolute Error (MAE), Root Mean Square Error (RMSE) (Harmel & Smith 2007) and Percent Bias (PBIAS) (Gupta et al. 1999) are applied to quantify the error between the observation and forecast. In addition, further accuracy measures such as the Brier Score (BS) (Brier 1950) and the Continuous Ranked Probability Score (CRPS) (Brocker 2012) are included in the verification. The CRPS can be calculated for single-valued and probabilistic predictions, allowing a direct comparison between the two types of forecasts. The forecast reliability is assessed with the False Alarm Ratio (FAR), whereas the discrimination attribute is evaluated with the Probability of Detection (POD) (McGovern et al. 2011). Finally, the model's goodness-of-fit is determined by the Index of Agreement (IoA), Coefficient of determination (R²) and Nash–Sutcliffe (NS) model efficiency (Harmel & Smith 2007).

All the metrics are collected in Table 5 together with their perfect scores. The formulas of the deterministic verification metrics are given in Equations (12)–(19), whereas the probabilistic verification measures are calculated as shown in Equations (20) and (21). The computed verification metrics are then used to assess the improvement in forecast performance. Table 6 shows the relationship between the changes in verification metrics, grouped by forecast attributes, and the corresponding changes in forecast performance.

Some of the above-mentioned metrics such as the POD and FAR use the comparison of the measured and forecast occurrence frequency of a predefined event. In this study an event is defined when the Chlorophyll-a concentration in the top layer of the water column exceeds the elevated assessment levels. The mean elevated assessment level during the growth season (March–September) is 7.5 mg/cm³ in the Dutch coastal waters and 12 mg/cm³ in the Dutch Wadden Sea. These levels were previously specified by Baretta-Bekker et al. (2008) for the application of the OSPAR Comprehensive procedure.

The above-presented ensemble method was tested for hindcasting the Chlorophyll-a concentration in the southern North Sea for different hydrodynamic years, first for year 2009 and then for year 2007. The ensemble method's performance is analysed separately in the Dutch Wadden Sea and Dutch coastal waters using three observation types with varying temporal coverage and accuracy. In-situ measurements are considered as the most reliable data but their temporal coverage is low, only 15 measurements per year on average. In comparison, the remote sensing data has increased temporal coverage, 30 MERIS observations per year and 200 gap-filled MERIS observations per year on average, allowing better verification of the model but introducing more uncertainty.

Ensemble forecasts are often assessed by verifying how well the ensemble band could capture the observations. Figure 6 and Table 7 show the percentage of measurements captured by the ensemble band and the confidence intervals (50%, 75%, 95%) at the in-situ stations. Overall, it can be concluded that the percentage of measurements lying within the ensemble band is relatively low, 60% on average considering all stations and both years. The percentage results for the 95% and 75% confidence intervals are even lower, but the accuracy of the former is still acceptable for the research. Nevertheless, particular differences can be observed in the ensemble forecast's skill once the two areas are analysed separately. In the Wadden Sea only 35 to 57% of the measurements could be captured in the ensemble band while in the coastal waters this percentage is as high as 67 to 80%.

In this paper, processed remote sensing images are also applied as a verification dataset. Given the high degree of uncertainty in this dataset, it is advised to provide an error estimate together with the results. The MERIS data are compared with the in-situ measurements to obtain the error estimate at the stations where a sufficient number of matchups are available. Considering all stations, the MERIS dataset has an average MAE of 41% in 2009 and 48% in 2007, whereas the DINEOF gap-filled MERIS data's average error is 52% and 48% for the same years. Thus, applying averaged error bands for the remote sensing dataset instead of the uncertain single point values might provide a better picture of the ensemble forecast's performance. The reader may refer to Harmel & Smith (2007) for further information about considering measurement uncertainty in the evaluation of goodness-of-fit metrics in water quality modelling.

Taking into account the measurement uncertainty in the remote sensing data, the percentage of captured measurements would increase considerably up to 84% in the Wadden Sea, and up to 89% in the coastal waters, however, 100% coverage could not be achieved. A possible reason for this might be that the coastal zone and the Wadden Sea are shallow, dynamically varying ecosystems with high turbidity and therefore the deterministic simulation, which serves as the base for the ensemble forecast, should be recalibrated here to achieve improved results. This model uncertainty should be considered when evaluating the results. Moreover, the fact that not all uncertainty sources were taken into account may explain the additional uncertainty stemming from the model structure and from the meteorological and hydrodynamic inputs. Finally, the presumed level of uncertainty in the inputs (e.g. SPM concentration field), might be underestimated or their assumed sample distribution may not be appropriate.

The direct comparison of the deterministic and ensemble forecast was achieved through calculating a selection of verification metrics (introduced above) for both types of predictions. In order to be able to compute deterministic scores for the ensemble forecast, the ensemble median (50-th percentile) prediction was selected in this study, although further options are available such as the ensemble mean or any other percentile. The ensemble median (or ensemble mean) usually verifies better than the deterministic forecast by most verification scores, because it presents the most predictable elements of the forecast and smooths out the extreme unpredictable elements. It indicates the future values of the model output variable that can be predicted with confidence, but it will rarely capture the extreme events, and therefore should not be relied upon on its own.

The percentage improvement in the groups of verification metrics if the ensemble median prediction is used instead of the deterministic forecast is presented in Table 8. It is important to note that all metrics are equally weighted within the groups. Nevertheless, it might be possible that in other studies the metrics would be weighted or some specific metrics would be omitted according to the forecasters' need. The group members are shown in Table 6. The verification metrics are averaged values of six stations in the Dutch Wadden Sea and nine stations in the coastal waters. Furthermore, the percentage improvements in the metrics are computed as an average of the results using all measurement types.

The goodness-of-fit and accuracy metrics show moderate improvement in both years and both areas but the magnitude of the improvement is markedly higher in the coastal waters nonetheless. The Brier Score is only calculated for the ensemble forecast and thus it cannot be compared to the deterministic prediction. Considering all measurement types the average Brier Score in the Wadden Sea is 0.13 in 2009 and 0.12 in 2007, whereas at the Dutch Wadden Sea stations the Brier Scores for the same years are 0.05 and 0.09. Reliability and discrimination measures only experience a minor improvement in the coastal waters and no improvement at all in the Wadden Sea. The verification results confirm the previous finding that the ensemble forecast performs better in the Dutch coastal waters than in the Wadden Sea. In the problematic Wadden Sea area, the ensemble method's efficiency is less convincing; on the other hand, in the coastal waters the preliminary results are promising.

Making use of the three-dimensional water quality model the forecast can be visualized as a Chlorophyll-a concentration map. As mentioned above, the ensemble median prediction should not be relied upon on its own but with the indication of the degree of certitude. The ensemble spread can provide a measure of the level of uncertainty in the output parameter; hence, it is appropriate to complement the ensemble median. Figure 7 illustrates the ensemble median prediction map at a specific time-step together with the ensemble spread map calculated as the standard deviation of the Chlorophyll-a concentration.

Observing the ensemble spread map allows us to draw conclusions on the prediction's spatial uncertainty. In late March, during the peak spring algal bloom, the prediction in most off-shore areas has low standard deviation, while significant spread can be observed in the near shore, shallow areas of the North Sea and in the Dutch Wadden Sea. This might be partly due to the direct effect of perturbations on the river loads in the focus area and the high spread may indicate the rivers' zone of reach. Nonetheless, it could be also explained with the specific regional system dynamics, since in the near shore shallower zones the water is low-dynamic and the primary production levels are already elevated.

This paper can contribute to existing knowledge in probabilistic forecasting by providing an application to water quality prediction with a three-dimensional ecosystem model. The promising results indicate that ensemble prediction techniques can produce enhanced forecast of water quality indicators due to their ability to account for the error in the model output variable. Therefore, it is suggested that in the future, by further developing the proposed forecast set-up, the implementation of operational ensemble forecasting systems for coastal water quality prediction is viable. Nevertheless, the potential of the ensemble forecasting system ultimately depends on the input ensemble. Moreover, the sampling technique that is used to generate the input ensemble should be tailored to the specific requirements of the application.

In the presented case study, the ensemble forecast moderately outperforms the deterministic forecast at the coastal stations, which might be advantageous in decision making if the underlying baseline (deterministic) model is sufficiently well calibrated and validated. Moreover, the degree of certitude and the likelihood of an event, which could be expressed through ensemble forecasting, provide opportunity to set risk-based criteria for the response measures. In other words, the uncertainty estimate produced by the proposed ensemble forecast promotes rational decision making, and offers potential for additional ecosystem and economic benefits.

Due to their fundamental importance, the identified significant parameters should be revised based on a comprehensive sensitivity analysis. In addition, the parameter value ranges and prior distributions should be re-evaluated for future applications. Furthermore, the uncertainty from the meteorological and hydrodynamic inputs could also be included in the ensemble generation in order to better estimate the model input uncertainty. The spatial correlation of the SPM concentration field should be also considered to achieve a more realistic SPM perturbation.

Further necessary advancement of the presented ensemble prediction system would be to develop the method into an operational forecasting system. This could be achieved by implementing the proposed technique in a platform (e.g. Delft-FEWS) which helps to organize the integrated modelling system elements and the various types of observations from data importing and processing through the ensemble simulation with the GEM to the model output post-processing. Finally, a statistical post-processing technique could complement the ensemble forecast in order to estimate the total predictive uncertainty. This approach includes statistical post-processing uncertainty estimation for each ensemble member independently for each lead time. The product is a separated PDF for each ensemble member which can then be averaged to obtain the total predictive uncertainty.

This study was conducted within the framework of the ECOPOTENTIAL project. This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No. 641762. The authors would like to express their gratitude for the contribution that was provided by Alexander Ziemba.