Trace Metal Modelling of a Complex River Basin Using the Suite of Models Integrated in the OpenMI Platform

Modelling trace metal dynamics is essential in an integrated modelling framework as trace metals have the potential to be fatal, even when present at low concentrations. Since the degree of bioavailability of a metal depends on its presence in the dissolved phase, it is necessary to keep track of both the dissolved and particulate phase of metals. In general, the well-known partitioning coefficient approach is widely used for trace metal speciation. As such, we applied a parametric approach to relate the partitioning coefficient to several environmental variables. These environmental variables are made available by a suite of physically based models (a hydrologic and diffuse pollution model, Soil and Water Assessment Tool (SWAT); a hydraulic model, Storm Water Management Model (SWMM); a stream temperature model; an in-stream water quality conversion model; and a sediment transport model) integrated using the Open Modelling Interface (OpenMI). For trace metal speciation, two regression techniques, (a) the multi-linear regression (MLR) and (b) the principle component regression (PCR), were used. It is then tested in the Zenne river basin, Belgium, to simulate four trace metals (copper, cadmium, zinc and lead) dynamics. We demonstrated the usefulness of the OpenMI platform to link different model components for integrated trace metal transport modelling of a complex river basin. It was found that the integrated model simulated different metals with ‘satisfactory’ accuracy. The MLRand PCR-based model results were also comparable. From a management perspective, the river is not heavily polluted except for the level of dissolved zinc. We believe that the availability of such a model will allow for a better understanding of trace metal dynamics, which could be utilized to improve the present condition of the river.


Introduction
The European Union (EU) Water Framework Directive (WFD) established in 2000 [1] advocated the requirement of a good ecological status for all surface waters in all member states by 2015.In order to attain this required status, the directive has specified Environmental Quality Standards (EQS) for several priority substances, including trace metals [2].In response to the WFD's call for enhancing the status of aquatic ecosystems and reducing pollution from priority substances, large investments were made for the management of the wastewaters draining to the Zenne River, Belgium, resulting in an increase of the quality of the river [3].Despite these investments, the river still receives high loads of pollutants, especially considering its low discharge.Consequently, the water quality downstream from Brussels does not comply with the requirements set by the WFD.It is in this context that an interuniversity, multidisciplinary research project titled, 'Good Ecological Status of the river Zenne' (GESZ), was launched to evaluate the effects of the wastewater management plans on the ecological functioning of the river.It aims to provide advice related to sewage management processes in Brussels through the development of various tools that can compute the evolution of water quantity and quality in the river by measuring multiple variables, including trace metals.
The need for integrated models to gain insight into environmental problems is well appreciated in recent times [4].Modelling of trace metal dynamics is essential in an integrated framework because some trace metals, such as mercury (Hg), lead (Pb) and cadmium (Cd), have the potential to be fatal even if they are present at low concentrations.This fact is reflected in the WFD, which specifies Cd as a priority hazardous substance and Pb as hazardous substance (Table 1).The WFD requires that these substances be progressively phased out.Note that not all of the trace metals are toxic; some of them are essential as micro-nutrients and are present in high concentrations [5].Since the degree of bioavailability and thus, the toxicity, of a metal depends on its presence in the dissolved and particulate phases, it is necessary to keep track of both forms of metals in the water column.[7]; d Aquatic Life Criteria US EPA [8], criterion continuous concentration.$ Criterion maximum concentration; # Basic quality standard, and yearly average value for rivers and oceans VLAREM [9].
The Zenne river basin, Belgium, can be considered a unique and complex river basin.From its upstream predominantly agricultural sub-basins, the river receives mostly the non-point sources pollution.The river is made to exchange water fluxes with a canal at multiple locations to prevent the city of the Brussels from flooding.The proportion of flows from Brussel's sewer network can be more than 50% of the river's flow, especially during dry weather flow conditions [3].Furthermore, the river flow is influenced by the tide of the North Sea at its mouth.Hence, a single model such as the Soil and Water Assessment Tool (SWAT) [10], a hydrologic model, is not suitable as the model cannot simulate either complex interactions of sewer systems or the back water effects of the tide.While hydraulic models such as the Storm Water Management Model (SWMM) [11] can simulate the complex sewer systems and the tidal effect, they are not suited for simulating the non-point sources of pollution.Therefore, in order to have flexibility during the selection of the simulators for the different processes occurring within the basin (i.e., hydrology of the river basin, hydraulics in the river and sewers, erosion and sediment transport, fecal bacteria transport, the fate of trace metals and the carbon-nitrogen-phosphorus cycle), we opted to use the Open Modelling Interface (OpenMI) [12,13].The OpenMI is an interface which allows for the dynamic linking of different models that are compliant to its pre-established standard.The OpenMI has been successfully and widely applied for the purposes of model integration in such areas as urban hydrology and hydraulics [14,15], sediment transport modeling [16], in-stream water quality modelling [17] and fecal bacteria modelling [18].Additionally, ample support for using OpenMI has come from multiple different groups and organizations [19][20][21].
Furthermore, there have been growing calls to follow the so-called 'open source route' in environmental modelling discipline so as to create a 'functional collaborative community' [4,21,22].As such, use of such platform (the OpenMI) which promotes the use of existing open-source simulators would enhance model re-usability, thereby saving huge amount of time and efforts in developing a new one [23].
In this paper, we used the OpenMI to integrate several physically based models: a hydrological model, the Soil and Water Assessment Tool (SWAT) [10]; a hydraulic model, the Storm Water Management Model (SWMM) [11]; a stream water temperature model; a sediment transport model; and a water quality model.The integrated model would then provide different environmental variables for the newly developed trace metal transport model which then predicts both dissolved and particulate phase of metals using different regression techniques.In order to be able to integrate the models on the OpenMI platform, SWAT and SWMM were migrated to the OpenMI platform, (refer to Shrestha et al. [24] for information on SWMM OpenMI migration).The migration process involved a re-structuring of the engine core of the simulator in such a way that allows it to inherit the running sequences of OpenMI [12].The other simulators were developed in such a way that they directly became OpenMI compliant.The integrated model was then used to simulate the trace metal dynamics of the Zenne river, Belgium.We believe that with the availability of such an integrated trace metal model, the trace metal dynamics can be better understood and such information could easily be translated into improving the present condition of the river.

Trace Metal Dynamics in a Riverine System
Trace metals enter the riverine system via two main sources: natural (e.g., weathering, erosion and volcanic actions) and anthropogenic (e.g., domestic and industrial waste water, mining and dredging) [25].In the river, they are partitioned between the dissolved phase and the particulate phase.Some major processes of metal partitioning are adsorption-desorption, complexation-decomplexation, precipitation-dissolution and biological uptake-organic matter degradation [26][27][28][29].Since the degree of bioavailability and thus, the toxicity, of a metal depends on its presence in the dissolved and particulate phases, it is necessary to keep track of both forms of metals in the water column.It should also be noted that the bed sediments constitute a reservoir for dissolved trace metals [30][31][32].
As per current modelling practices, the use of metal partitioning coefficients-also known as distribution coefficients-to partition the metals between the dissolved and the particulate form is widespread [33].The partitioning coefficient (K d in L/kg) is defined as the ratio of the particulate metal concentration (MeP in µg/kg of sorbed material) to the dissolved metal concentration (MeD in µg/L of water).In most cases, K d is considered a constant.However, there is growing realization that K d varies in accordance with changes in environmental variables, such as pH and dissolved oxygen (DO) [34,35].
Different techniques that may be used to simulate K d are: (a) parametric models; (b) isotherm adsorption models and (c) mechanistic adsorption models [36].A parametric model describes the variation of the distribution coefficient through empirically-derived relationships between K d and sediment and water quality variables.Such models are applicable only within the range of the values of the variables used to derive the empirical relationship.Furthermore, owing to the statistical nature of the equations, they do not provide an explanation of how metals associate with particulate matter [36].An adsorption isotherm model describes the association of the metals with the solids for constant values of variables such as pH and temperature [37].This kind of model provides information on the adsorption mechanism of the metals.The particulate metal concentration is evaluated based on the dissolved metal concentration and constant parameters.This simplifies the dependence of the metal concentration on variables such as organic matter [36].Isotherm models are available in a variety of forms, including the Languir model [38], the Freundlich model [39] and the Dubinin-Radushkevich equation [40].A mechanistic model takes into account the dependency of the distribution coefficient on a number of variables such as the metal concentrations, the concentrations of competing ions and the surface charge of the adsorbent [36].Although these models may be expected to provide better results than the parametric models or the isotherm models [41], they are complex and require a large number of input variables [36].
Once a trace metal is separated into the dissolved and the particulate forms, it needs to be routed.In general, its fate is governed by a 3-D advective-diffusion equation [42].However, several simplifications of the equation exist for various situations.For example, for a river that is long enough when compared to the distance of complete mixing over the entire x-section, a 1-D approach is appropriate [43].A further simplification for a situation in which a river is divided into several river segments and within which perfect mixing is assumed, the 1-D equation further reduces to an ordinary differential equation which is commonly known as the completely stirred tank reactor (CSTR).

The Study Area
The Zenne River is a part of the international Scheldt basin.It drains an area of 1162 km 2 and runs through the three administrative regions of Belgium: the Walloon Region (574 km 2 ), the Brussels Capital Region (162 km 2 ) and the Flemish Region (426 km 2 ).After a length of approximately 100 km, it meets with the Dijle River, where it is subject to the tidal influence of the Scheldt River.The river basin is crossed by the Brussels-Charleroi canal and the Brussels-Scheldt Sea Canal.In the Walloon Region, the canal is fed by former tributaries of the river (Figure 1).In the area between the borders of the Walloon and the Brussels Region, different tributaries discharge into the river.The canal also receives water from the river through several overflow structures in order to prevent flooding in the Brussels region.Excess water in the canal can be discharged back into the river downstream from Brussels.The basin has a population density of about 1260 inhabitants/km 2 , which makes it the most densely populated basin of the Scheldt basin [44].The land-use in the catchment is dominated by agriculture (51%), primarily in the upstream part of the basin and by urban regions (19%) in the downstream section.There are several waste water treatment plants (WWTPs) in the basin.In the Brussels area, the river receives the sewage waters from the WWTP Brussels South (360,000 inhabitant-equivalents, in operation since 2000) and the WWTP Brussels North (1,000,000 inhabitant-equivalents, in operation since 2007).
The Zenne River is also a fine example of how a river system's ecological and hydrologic continuum can be modified by anthropogenic pressures.The prosperous textile craft industry at the beginning of 11th century along the banks of the Zenne River led to a significant population increase.As a result of industrialization and associated urbanization and sanitation needs of the inhabitants living in the catchment, the river had already been heavily polluted [45,46].The drawn discharge from the Zenne to feed a canal which ran parallel to the river also contributed to worsening of the situation [47].Thus, the river was required to be covered for sanitation purposes.Once the river was covered, water quality preservation was no longer a priority for Brussels' authorities [3] until the introduction of enforcements imposed by the Urban Waste Water Treatment Directive (UWWT) [48] in 1991.

The Trace Metal Simulator
As stated in Section 2.1, there exists different techniques to simulate the partitioning coefficient (K d ).We opted for a parametric model because various input data required for other techniques (e.g., isotherm adsorption models and mechanistic adsorption) were not available.The parametric model for K d was based on readily available (monitored) or computable environmental variables.As we intended to use the trace metal model as a component model in an integrated modelling chain, a simple yet robust model was needed to limit the computation time which is considered to be a severe bottleneck of an integrated model.
In this study, we applied two regression techniques: a multilinear regression (MLR) and a principal component regression (PCR).In an MLR, the dependent variables are determined directly from the independent variables while in a PCR, new components are defined to be used as the independent variables.For both the MLR and the PCR, we first determined an empirical relationship for the calculation of K d based on physicochemical variables and total metal concentrations (MeT).The considered physicochemical variables were: electrical conductivity (EC), dissolved oxygen concentration (DO), the pH, the suspended particulate matter (SPM) concentration and the temperature (T s ).Interaction terms were introduced and selected based on the assessment of the variance inflation factor (VIF) [49].For an MLR, interaction terms with excessive multi-collinearity (VIF > 10,000) were eliminated.This retained the terms with smaller VIF values which indicated lessened multi-collinearity.Additionally, the Bayesian information criterion (BIC) [50] was used to determine the model parameters that provided a good compromise between goodness of fit (R 2 ) and predictive power reflected by the predicted residual sums of squares (PRESS) statistic [51].The model with the lowest BIC was selected as the final model.
It should be noted that we transformed all of the variables before they were subjected to regression using a Box-Cox (BC) transformation [52] to ensure a normal distribution of the variables.The parameter (λ-value) of the BC-transformation was subjected to calibration.We found the parameters to be −0.208,−0.403, −0.065 and −0.204, respectively, for MeT-Cu, MeT-Cd, MeT-Pb and MeT-Zn.Similarly, λ-values of 2.057, 1.101, 0.396, −0.632 and −0.359 were found for T s , pH, DO, EC and SPM, respectively.The assumption of normality was checked using the Shapiro-Wilk test [53].
Once a function for the partitioning coefficient, K d (L/kg), is defined, MeD (µg/L) and MeP (µg/g) can be derived from the total metal concentration, MeT (µg/L) and the suspended particulate matter, SPM (mg/L), concentration as: (2) Furthermore, we adopted the principles of CSTR to route the metal species.The metal species were considered not to undergo any conversion process within the time frame of relevance.
An accurate information of water fluxes was required in order to route the metal constituents.Owing to the nature and complexity of the river basin under consideration, we selected the SWAT model for the simulation of the river basin processes from the predominantly rural sub-catchments in the upstream part of the river basin.SWAT is a physically-based, semi-distributed, hydrologic simulator and one of the most widely used models to simulate the basin processes [10].The SWAT, being a hydrological model, cannot simulate complex hydraulic processes such as back water effects, nor two way interactions between, for instance, a sewer system and a river.We used a hydraulic simulator, the SWMM, for the downstream reach of the river to address all of the different components of the system (e.g., the sewer systems with their auxiliary structures, the canal and the river in tidal reach).The SWMM is, primarily, a dynamic rainfall runoff simulator for computing the runoff quantity and quality from urban areas [11].The SWAT-SWMM combination was therefore deemed necessary to simulate water fluxes for a complex basin such as the Zenne.
An additional simulator needed to be developed for the sediment as we intended to have a trace metal model which would track both dissolved and particulate forms (see Equations ( 1) and ( 2)) using the partitioning coefficient.The sediment transport simulator we developed is based on Shield's approach for the initiation of motion based on critical shear stress and critical particle diameter in conjunction with Velikanov's energy equation for the limitation of the carrying capacity.A thorough description of the simulator is presented in Shrestha et al. [54].
Since the partitioning coefficient depends on various physiochemical variables, a separate simulator for the water quality conversion processes needed to be developed.We used river water quality model No. 1 (RWQM1) principles for the water quality conversion processes.A thorough description of the simulator is presented in Shrestha et al. [17].
Additionally, since many kinematics of the sediment transport simulator and the water quality simulator depended on the stream water temperature, a temperature simulator also needed to be developed.The temperature simulator we developed is based on a non-linear regression fit between the air and the stream water temperatures, as suggested by Mohseni et al. [55].Further information of the simulator can be found in Shrestha et al. [56].Finally, a trace metal transport simulator, described in more detail in following sections, was developed to simulate trace metal dynamics and speciation.

Model Inputs, Build up, Calibration and Validation
Table 2 provides an overview of the types of data that were used for the build-up, the calibration and the validation of the developed models.The following paragraphs elaborate the issues for each model component.Emphasis is placed on the inputs, build-up, calibration and validation of the trace metal component.
The SWAT model is the first model of the integrated modelling chain.Hence, it is built up, calibrated and validated on a standalone basis.The data sources and the data resolutions used for the model are presented in Table 2.We used the years 1998-2008 for calibration and the years 1988-1997 for validation, using daily stream flows recorded at Tubize (Figure 1).We optimized the related parameters (see Leta et al. [15] and Leta et al. [57] for details).For sediment and water quality simulations, we used the years 1998-2008 for calibration and the years 1994-1997 for validation.For the calibration of the sediment module, we used sediment concentrations measured at Quenast (Figure 1).Five water quality variables (DO, BOD, NH 4 , NO 3 and DIS-P) measured at Quenast were used for a manual calibration of the water quality-related parameters.
A SWMM model was built for the section of river between the border of the Walloon-Flemish region and for the mouth of the river (Figure 1).The canal network, including the exchange points with the river, was also modelled using SWMM.The SWMM model included simplified representations of the sewer systems that are connected to the major WWTPs (Figure 1).As SWMM requires rural catchment runoff as a boundary condition and the runoffs were simulated using the SWAT model, an OpenMI-integrated SWAT/SWMM model was implemented.The integrated model was calibrated and validated using hourly stream flows observed at three locations along the river-Lot, Vilvoorde and Eppegem (Figure 1).We used the stream flow data from 2007 to 2008 to calibrate the model and flow data from 2009 to 2010 for the validation.As the next step, we added the temperature model to the SWAT-SWMM integrated model.We used stream water temperature data at three locations along the river, Lot, Vilvoorde and Eppegem (Figure 1) from the years 2007 to 2010, as recorded by the VMM, to validate the model.
Next, the sediment model was added.The quantities of flow and sediment exchanged at different element sets such as links, nodes and sub-basin outlets are shown in Figure 2. Imposed boundary conditions for the integrated model are also shown in Table 2.We calibrated the parameters (see Shrestha et al. [58] for details) of the model using sediment concentrations for Lot, Vilvoorde and Eppegem from 2007 to 2008.The integrated model was validated using sediment concentrations from 2009 to 2010.
For the water quality model, SWAT-SWMM-Temperature-Water Quality components were integrated in OpenMI.Figure 2 shows the data exchange between the different model components that were used for the simulation of the water quality.Imposed boundary conditions for the model are presented in Table 2.The parameters of the water quality model (see Shrestha et al. [17] for details) were determined heuristically, based on measured water quality variables during dry weather flow (DWF) conditions.The model was then validated based on long term simulations, using monthly water quality concentrations measured by the VMM at Lot, Vilvoorde and Eppegem (Figure 1).Finally, an integrated trace metal transport model was implemented which consisted of six model components (Figure 2).We considered of metals: copper (Cu), cadmium (Cd), zinc (Zn) and lead (Pb).According to the WFD, Cd is considered to be a 'priority hazardous substance' and Pb is considered a 'priority substance'.
In this model build up, we used data gathered from the GESZ campaigns.It included longitudinal profiles of various physicochemical variables and MeT concentrations measured during dry weather conditions in different seasons (five campaigns).In addition, two campaigns were organized to gain information about the water quality during storm conditions, including a 24-h campaign.In order to calibrate and validate the derived relationship, we divided the data into two groups.'Group-1' (Table 3) consisted of 113 samples taken from four DWF campaigns and one 24-h set of continuous measurements.This group made up the calibration data set.'Group-2' consisted of 39 samples taken from one DWF campaign, one rain campaign and 10 randomly selected datasets from 'Group-1' in order to obtain wider variability in the data set.'Group-2' was used for validation purposes.The decision to select two sampling campaigns for 'Group-2' was made in consideration of the fact that these measurements showed clustering, as shown in the score plot (Figure S1) where four components F1 to F4, describing 83% of the total information from the data, are plotted.A score plot shows a map of these observations.Data points close to the origin display average properties of the data while data grouped together or in a quadrant could indicate data anomalies.Similarly, outliers could be detected too.As can be seen on both plots (Figure S1), the data corresponding to sampling campaign of December 2009 fall in the third quadrant for both plots.Similarly, a clear cluster corresponding to data of June 2011, a wet weather campaign, is especially visible in Figure S1a.In addition to the rural catchment MeT inputs which were estimated with the derived regression equations, we used constant MeT concentrations at WWTP outlets, at CSO sources and at the outlet of the river based on the GESZ measurements.Figure 3 shows the GESZ measurements of total metal concentration at the outlets of WWTP and at GESZ station Z13, located near the river outlet.Hence, mean observed MeT concentrations were imposed at the river outlet during the sampling campaigns.Although there exists variation among the total metal concentrations, we had to impose an average constant concentration due to a lack of an observed time series.Although we are aware that more accurate boundary conditions needed to be imposed at locations such as the WWTP outlets, the CSO outlets and the river outlets, we did not have another viable option.The integrated trace metal model was validated using MeT and MeD concentrations measured by the VMM measurements at three stations: Lot, Vilvoorde and Eppegem (Figure 1).Particulate metal (MeP) concentrations of metals were not available for the considered data period of 2007-2010.Moreover, several measurements of total and dissolved trace metal concentrations were below the detection limit and, therefore, had to be eliminated from the dataset.

Performance Evaluation of the Model
To evaluate the accuracy of model results, we used the percentage bias (PBIAS), the ratio of the root mean square error to the standard deviation of the measured data (RSR) and the Nash-Sutcliffe efficiency (NSE) [59].The main reason for using these statistics was to assign a qualitative rating to the model results.We used the ranges for these goodness-of-fit statistics as proposed by Moriasi et al. [60] to determine a qualitative assessment of the model (i.e., very good, good, satisfactory, or unsatisfactory) for each of these indicators.A global qualitative assessment of the model was made by assigning a score of 4 to a very good model and assigning a score of 1 to an unsatisfactory model for each of the criteria and then calculating the average score for the three criteria.The widely-used root-mean squared error (RMSE) and mean absolute error (MAE) were used to complement the aforementioned dimensionless goodness-of-fit statistics.

Results and Discussion
For clarity purposes, detailed results of each model component are not presented here.Rather, a qualitative assessment of the different components is provided.For the optimized stream flow related parameters of the SWAT model, we refer toLeta et al.Optimized flow-related parameters for one of the collectors of the sewer systems of Brussels-the Paruck collector-can be found in Nkaika et al. [60].Similarly, water quality parameters of the SWAT model and parameters of the integrated water quality model can be found in Shrestha et al. [17].Optimized parameters of the sediment transport and temperature models can be found in Shrestha et al. [54].The following paragraphs, therefore, focus on the detailed results of the trace metal model.

The Total Metal Concentrations at The Outlets of Rural Catchment
Table 4 provides an overview of the regression equations for the different metals.In the case of the MLR, the use of EC, pH and SPM as predictors in various combinations resulted in the lowest BIC for the MeT-Cd (Figure S2).For MeT-Cu, two variables (EC and SPM) were the best set of predictors.For MeT-Pb, four variables (EC, SPM, T s and pH) were selected.For MeT-Zn, three variables (EC, pH and SPM) were the best set that could be determined.DO did not appear in the set of predictors for any of the metals.EC and SPM influenced the MeT concentrations the most negatively and positively, respectively.The results were expected as the loading plot (Figure S3) showed that MeT and SPM lie on the same quadrant, while MeT and EC do not.Indeed, a loading plot showed relationships amongst the variables that were taken into consideration.Variables close to one another exhibit high correlation while those on the opposite side of the origin exhibit negative correlations.Moreover, similar samples are placed together.Hence, distance from the origin to the variable and its quadrant location carries special meaning.
The quality of the regression models varies from 'Unsatisfactory' to 'Very Good,' depending on the considered metal, the analysis technique used and the dataset (Table S3).A graphical representation of the results is shown in Figure 4a.In the figure, most of the data points fall inside the +/−1-standard deviation line (i.e., within 68% of the confidence interval) of the predicted value.For all cases, the predicted MeT were underestimated for higher MeT observations.We found that these higher MeT observations corresponded to one specific sampling campaign, the 24-h campaign (Group-1, Table 3) and at one specific location, Z4.It should be noted that the SPM appeared in the regression equations for all MeT (Table 4) and influenced the MeT predictions positively (i.e., higher SPM resulted in higher MeT-see the regression equation in Table 4).We found that the measured SPM in that particular sampling campaign and at that particular location was lower (29 ± 7 mg/L) than the DWF average value (43 ± 35 mg/L).The underestimation could, thus, be related to this.
Globally, it may be concluded that the results of both analysis techniques are comparable.From the statistical indicators and graphical plots, it can be derived that the PCR technique performed slightly better for Cd and Zn, while MLR performed better for Pb.However, a general trend that can be derived from Table 4 is that the PCR fared better during the calibration period (i.e., had higher R 2 values) than did the MLR.On the other hand, the predictive capability of PCR was generally lower than the MLR's as indicated by higher PRESS statistics.Hence, when the derived relationships were subjected to a different range of independent variables, the multi-linear regression (MLR) was expected to fare better than the PCR.

The Distribution Coefficients
In the loading plot (Figure S3), the EC and the Log K d are lying in opposite quadrants, meaning that the EC has a negative influence on the LogK d .The regression equations (Table 4) confirm this, especially for Log K d -Cd and Log K d -Zn.The quality of the regression models ranged from 'Satisfactory' to 'Very Good' (Table S3). Figure 4b shows predicted and observed Log K d data sets.The performances of the MLR and PCR derived Log K d equations are fairly comparable.For most cases, however, it appears that the PCR-based regressions provided somewhat better fits, as indicated by higher values of R2, lower values of PRESS (Table 4) and higher performance ratings (Table S3).However, the goodness-of-fit statistics for Log K d -Zn show that, for both the MLR and the PCR analysis, the correlation was significantly lower than that for the other metals.While the average tendency was represented fairly well (PBIAS < 4%), the values of RSR and NSE, indeed, pointed to a lower accuracy.The NSE was below zero for the PCR derived Log K d -Zn, even for the calibration period.The negative values for NSE and the higher values for RSR during the validation period indicate that the regression equations for all of the Log K d performed poorly when subjected to different conditions than those used for the calibration.

The Results of the Integrated Trace Metal Model
As already indicated, the predicted metal concentrations are validated using observed concentrations at three stations: Lot, Vilvoorde and Eppegem.For clarity purposes, we present simulated and observed total, dissolved and particulate metal concentrations at the most downstream stations-Eppegem (Figure 5).Table 5 shows the goodness-of-fit statistics of the simulated results.It should hereby be noted that the quality of the particulate metal concentrations simulations could not be evaluated as the latter concentrations were not measured.
As evident in Figure 5, at Eppegem, the PCR-based results matched the observed total copper levels well, while the MLR-based simulations were systematically higher than the PCR-based simulations.On average, the bias for the MLR simulations was ca.12%, compared to ca. 4% for the PCR simulations (Table 5).The error statistics (MAE and RMSE) also pointed to a slightly better quality of PCR simulation.On average, the goodness-of-fit statistics pointed to a 'Satisfactory' quality of both model results.For the distribution coefficients, the PCR-based values were systematically higher than the MLR-based simulations.Accordingly, the MLR-based dissolved Cu levels were higher compared to the PCR simulations.For the dissolved Cu, the PCR-based results matched the (2) observations best.The number of observations was, however, too low to draw straightforward conclusions on this issue.Both models produced comparable results for particulate Cu.
For total cadmium, both models showed comparable results.However, this was not the case for the distribution coefficient.The PCR simulations showed a higher variability (4.9-7.4;average 5.5) compared to the MLR simulations (5.9-6.8;average 6.1).The average and the variability of the PCR simulations was comparable to that of the GESZ measurements (4.3-7.1;average 5.5).In this regard, we can consider the PCR predictions better than the MLR predictions.The overestimation of the MLR-based distribution coefficient could be related to the overestimation of dissolved oxygen by the model for Eppegem (PBIAS ≈ 42%-Table 5), as the dissolved oxygen was positively correlated to the distribution coefficient for Cd (Table 4).Better prediction of the PCR simulations was also reflected in the dissolved Cd simulations (Figure 5), where the PCR simulations captured the two observations while the MLR simulations underestimated them.As expected, the MLR-based particulate Cd concentrations were slightly higher than the PCR-based concentrations.
The total zinc concentrations simulated by both models were fairly comparable (Figure 5).While the goodness-of-fit statistics pointed to the same overall quality ('Satisfactory') of the model results, the MLR-based results had slightly lower PBIAS and better RSR and NSE values.As was also evident in the plot, both models significantly overestimated the low Zn concentrations.As was also the case for Lot and Vilvoorde, the distribution coefficients simulated by both models were different.The MLR-based simulation showed a quasi-constant distribution coefficient (4.3-4.8;average 4.4), which was in contrast to the higher variability and average that was observed during the GESZ campaigns (4.0-7.4;average 4.9).The PCR-based simulation showed a higher variability (1.5-5.16;average 4.3).However, the average was too low when compared to the GESZ measurements.This rather poor quality of the models was somehow expected, as both the MLR-based and the PCR-based regression equations showed low R2 values, high PRESS values (Table 4) and a poor overall quality (Satisfactory; Table S3) compared to other metal.As for the total Zn concentration, the dissolved Zn concentration was overestimated by both models.On average, the overestimation from the PCR-based model (PBIAS ≈ 33%) was higher than that from the MLR-based model (PBIAS ≈ 30%).A slightly better performance of the MLR-based model was evident in the MAE and the RMSE.Both of the models tended to represent the trend of the total Pb observations fairly well (Figure 5).On average, the PCR-based simulations showed a lower bias (ca.3%) and lower error statistics (MAE and RMSE values of ca. 14 and 19 µg/L, respectively) than the MLR-based simulations (PBIAS = 9%, MAE = 16 µg/L and RMSE = 20 µg/L).The goodness-of-fit statistics pointed to a 'Satisfactory' quality of the results for both models.The distribution coefficient results showed that the MLR-based simulations (average: 6.9) were comparable to the PCR simulations (average: 6.8).For the dissolved Pb, the goodness-of-fit statistics pointed to an 'Unsatisfactory' quality of both models, whereby the MLR-based results were slightly better than the MLR-based results (Table 5).
Unlike the situation at Lot where rural catchment runoff dominates, the metal dynamics at Vilvoorde, located just outside Brussels, were dominated by WWTP effluents and the sporadic CSO events because the Brussels region added higher water volume than what was entered into it (Table S1).In Brussels, the main pollution sources were the WWTPs discharge and the CSO, occurring during wet events.Hence, the inability of the integrated model to simulate the trace metal dynamics at Vilvoorde with high accuracy (Table 5) was not only due to the propagation of the errors observed at Lot but also due to assumptions inherent in imposed concentrations at the WWTPs effluent and the CSO outlets.Especially given that the streamflows at Vilvoorde (Table 5) were simulated with 'Very Good' and 'Good' accuracy, this rules out the possibility of over-and under-dilution effects.The only major pollution source located after Vilvoorde was the Woluwe CSO outlet.The WWTP effluent of a small plant at Grimbergen was another pollution source.As with Vilvoorde, the inaccuracy of an integrated model to represent the metal dynamics with high accuracy was due to the propagation of model errors and to the underlying simplifications on the imposed boundary condition at the WWTP effluent and the CSO outlet.

Ecological Status of River Zenne Based on Long Term Simulation Results
The ecological status of the Zenne can be categorized by comparing the long-term (2007-2010) simulation of total, dissolved and particulate metal concentrations against the mandatory levels of different legislations (Table 1).To do so, the simulated instantaneous concentrations were compared against US EPA drinking water quality standard [7] and against the US EPA aquatic life standard [8].Similarly, the annual average concentrations were calculated to compare these against the Flemish standard [9], the EU ecological quality standard [2] and the EU standard for the quality of fresh waters to support fish life [10].
At Lot, only the dissolved Cadmium levels exceeded the US EPA drinking water quality standard.However, the frequency was very rare.In fact, the level was crossed only three times over the entire time span of 2007-2010.Similarly, the dissolved Zinc crossed the US EPA aquatic life standard only in the summer season in the PCR-based simulation.Moreover, the annual average of dissolved Zinc concentration has, on average, crossed the Flemish basic quality standard (Table 6).Despite Vilvoorde's location downstream of major pollution sources (e.g., effluents of major WWTPs and combined sewer flows during wet weather), the situation was not exceptionally alarming.In fact, the ecological status of the river was comparable to that of Lot.However, the levels of dissolved Zinc were almost three times that of the levels at Lot.Hence, the instantaneous concentrations have frequently crossed the US EPA aquatic life standard for all seasons and the annual average concentrations were consistently higher (3.3 times on average compared to Lot concentrations) than the Flemish basic quality standard (Table 6).
At Eppegem, the level of dissolved cadmium has, indeed, was higher by 3.5 times, on average, compared with concentrations at Lot and 1.5 times, on average, when compared with concentrations at Vilvoorde.However, it has seldom crossed the US EPA aquatic life criteria standard (Figure 5).Similarly, the dissolved Zinc concentrations continued to increase, crossing the US EPA aquatic life criteria standard for almost the entirety of the time (Figure 5).The annual average concentrations were always higher than the Flemish basic quality standard.In fact, the dissolved Zinc concentrations at Eppegem increased by 1.5 fold and 3.5 fold when compared with concentrations at Vilvoorde and Lot, respectively (Table 6).Overall, the ecological status of the Zenne River was good across all reaches of the river, except for the dissolved instantaneous Zinc concentrations which tended to cross the US EPA aquatic life standard.The annual average concentrations also tended to be higher than the Flemish basic quality standard.

Limitations of the Study
Models used to represent complex and interrelated environmental systems are always a simplification of reality.The level of simplification depends on the intended use of the model and availability of data to build-up, calibrate and validate the model, as well as on external resources (e.g., time and money).While we tried to represent this complex river system using a combination of different existing models (e.g., SWAT and SWMM) and new models (temperature, sediment transport, water quality conversion and trace metal), there were some underlying simplifications that needed some clarification.Most of these assumptions were forced due to data limitations.Furthermore, those who are involved in environmental modeling, are aware of the difficulty of acquiring high quality data of water quality variables, especially metal species.Generally, most of the data are collected to fulfill certain specific goals like we did in the scope of the GESZ project (refer to Table 3 for the sampling campaigns conducted in the scope of the project).Moreover, data for some extreme events, such as the CSOs, are normally difficult to monitor unless an automatic sampler is installed or a dedicated crew is available on a stand-by basis [61].
The major limitations of this work were the assumptions made on the imposed boundary conditions at important pollution sources (e.g., at WWTPs, CSO and the river outlet) which were kept constant in our current implementation and were obvious simplifications of reality.To our knowledge, the management organizations of the Brussels' sewer systems and thus, the WWTPs, have no plan to monitor trace metal dynamics for either of the WWTPs effluents or for the CSOs any time soon.
As the WWTPs inside Brussels were especially large and the influents during the dry weather conditions were quite constant, larger variations during storm conditions were dampened due to higher residence time (order of magnitude was about a week).At CSO sources, the variations in an urban drainage system such as the one in Brussels, will of course be large as have been observed by various researchers such as Lee et al. [62].The peak concentration of a particular pollutant will depend on factors such as the spatial extent of the system, the complexity of the drainage network, rainfall intensity and antecedent dry days.Hence, researchers have traditionally made use of the Event Mean Concentration (EMC) [63].Since the trace metal concentrations during CSO events were determined based on instantaneous measurements, a true EMC could not be calculated.While we are aware of the importance of CSO events as short-term intense pollution sources, it is evident, in our case, that the CSO's volumetric contribution in an annual scale was not significant (Table S1 and Figure S5).To further investigate this issue, we arbitrarily increased the metal concentrations at the CSO outlets by 20% and 50% relative to currently imposed concentrations and found that the trace metal concentrations at a downstream station (Eppegem) would increase by, on average, 5% and 12%, respectively.Had we missed the CSO concentration by 50%, the effect would still have been within the measurement errors of most water quality variables (~15%, Table S2).Furthermore, the governing legislations such as the Flemish (VLAREM, 2010) and the EU (EU, 2006; EU, 2008) impose annual averaged limits to the trace metal concentrations in the Zenne River.Finally, the constant metal concentrations imposed at the river outlet (Figure 3) only affected the river up to 9 km upstream from the outlet.Indeed, hydraulic simulation results showed backwater effects up to that distance (Figure S6).For an accurate boundary condition at the river outlet, a 2-D model of the Dijle River (the outlet of the river Zenne) would have needed to have been extended to the Scheldt estuary and coupled with our integrated model, which was not within the scope of this work.
Furthermore, we did not conduct a thorough evaluation on the suitability of different integrated modelling platforms, which is one of the limitations of this work.It should be noted that integrated modelling platforms such as The Invisible Modelling Framework (TIME) [64], Framework for Risk Assessment Analysis of Multimedia Environmental System (FRAMES) [65], The object Modelling System (OMS) [66] and so forth, are widely used for model integration.Our choice of the OpenMI as model integration platform was slightly influenced by its growing popularity in Europe and many projects in the region such as SEAMLESS [67], OpenMI-Life [68], SENSOR [69], NITROEUROPE [70], EFORWOOD [71] and so forth, used the OpenMI to integrated different legacy and new models.
Moreover, it should be noted that the EU WFD [1] explicitly specifies 'quality elements for the classification of the ecological status' as: (a) biological elements (composition and abundance of aquatic flora and fauna), (b) hydro-morphological elements supporting the biological elements (hydrological regime, river continuity and morphological conditions) and (c) chemical and physicochemical elements supporting the biological elements (thermal, oxygen, nutrient, acidification and salinity conditions) of waterbodies.The directive also specifies a list of priority substances which needs to be systematically phased out.It is thus clear that the directive gives biological analyses precedence over chemical ones.To this end, it can be argued that bioavailability of trace metals should have been validated by analyzing multiple sentinel organisms as various processes and factors affect the bioavailability of different trace metals to different organisms [72][73][74][75][76][77][78].This is another limitation of this work.

Usefulness of the OpenMI Platform for Integrated Modelling of a Complex River Basin
In light of the growing realization over the need for an integrated approach to holistic water management of river basins [4,22], it was demonstrated that use of loosely coupled, component-based, 'standard' environmental modelling framework-the OpenMI is feasible in such a complex river basin.Moreover, there is a need of promoting, practicing and sharing 'open source' developments among the practitioners [21] and OpenMI standard interface is nicely tailored for such purposes.Use of such 'standard' would inevitably need some extra effort to conform the design of the (standard's) engine interface.It can be stated that even for a user which is not a software specialist (but has some level of knowledge and experience in programming), there is quite a substantial support in terms of documentations, demonstrations and training materials in assisting for the development of OpenMI-compliant components [21].The main benefits of the OpenMI, as the author believes, are that it allows reusability of legacy models which are the results of substantial investment of time and money and it is an 'open standard' thereby allowing to share developments with wider user communities.While concerns have been raised regarding the OpenMI based data-exchange induced computation time overhead (e.g., Buahin & Horsburgh [79]), in our case, we reported in Shrestha et al. [54] that the overheard related to the run-time communication was found not to be significant (just 14% increased due to OpenMI based data exchange).However, with the advancement in computing technology, one can adopt parallel computing, grid computing and so forth, to overcome the slightest overload on computational time.On the other hand, another big advantage of OpenMI, as we see it, is that it allows to check the accuracy of each component of an integrated model and to more easily locate the possible sources of errors in the model.For instance, in this integrated trace metal model, the accuracy of simulation results of SWAT, SWMM and temperature model components (and their combinations, links 'a' to 'f,' Figure 2) is, in general, either 'Good' or 'Very Good' (Table 5).The accuracy of model simulations tends to deteriorate for the sediment transport and water quality model components (Table 5).As the trace metal model component requires input of different environmental variables (links 'g' and 'i', Figure 2) from the sediment transport and water quality model components, the 'satisfactory' accuracy of the trace metal model can thus also be attributed to the results of these components, in addition to the trace metal model component itself.Hence, in order to improve the accuracy of the trace metal model simulations of the integrated model, besides overcoming all the limitations as stated in Section 3.5, one can focus on improving the robustness (hence the accuracy) of the sediment transport and water quality components or one can simply replace these components in a 'pick and plug' manner with a more robust component.This flexibility as offered by the OpenMI platform, we believe, is very important and significant in the field of environmental modelling.

Conclusions
In this paper, we demonstrated the usefulness of the OpenMI platform to link different model components for integrated trace metal transport modelling.This approach is innovative in the sense that the different model components run in parallel over the entire river system.We found such a manner of model integration to be very useful to represent a complex river system such as the Zenne River in Belgium.The integrated model consisted of a hydrologic and diffuse pollution model (SWAT), a hydraulic model (SWMM), a stream temperature model, an in-stream water quality conversion model, a sediment transport model and a trace metal transport model.For this, the existing simulators, SWAT and SWMM, were migrated into the OpenMI platform while water quality, temperature, sediment and trace metal simulators were conceived in such a way as to become directly OpenMI compliant.We tested the integrated model to simulate the dynamics of four metal species, copper, cadmium, zinc and lead, on the Zenne River.The trace metal simulator employed the well-known partitioning coefficient (Kd) approach for metal speciation between the dissolved and the particulate phase.Unlike most simulators for metal transport in rivers, the partitioning coefficients were not constant but were made dependent on environmental variables such as the suspended particulate matter concentration, pH, water temperature and electric conductivity.To this end, we applied a parametric approach to relate the partitioning coefficients to environmental variables for the Zenne River.Two regression techniques were hereby used: the multi-linear regression technique (MLR) and the principle component regression (PCR) technique.
In contrast to our expectations, the accuracies of both regression techniques (MLR and PCR) were fairly comparable.For the total metals originating from rural catchments, the PCR was more accurate, while the MLR was slightly superior for making predictions.For the distribution coefficients, the PCR-based regression results were slightly better.However, both models fared poorly for the distribution coefficient of Zinc.The simulations using the integrated model showed that the total metal concentrations were represented fairly well by both models.However, very limited data were available to assess the accuracy of the simulated dissolved and particulate metal concentrations.While the MLR-based and the PCR-based simulation results for copper and the lead were comparable, the results for the cadmium and the zinc were not.The PCR-based simulations of the distribution coefficient showed higher variability compared to the MLR-based simulation results, affecting the dissolved and particulate metal results.The global performance ratings, calculated using the values of the bias, the root mean squared error scaled by the standard deviation of the observations and the Nash-Sutcliffe efficiency, also pointed to a 'satisfactory' quality of most of the models.From a management point of view, the river is not heavily polluted.The instantaneous concentration of dissolved cadmium exceeded the US EPA aquatic life standard in only certain extreme events while the instantaneous concentration of dissolved zinc exceeded this standard more frequently, especially downstream from Brussels.Consequently, the annual average concentration of dissolved Zinc exceeded the Flemish basic quality standard at all reaches of the river.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2076-3298/5/4/48/s1.Table S1 Yearly volume (in million cubic meters) of water flowing through the selected stations/points for the years 2007-2010.Table S2: Increment on total metal concentration (µg/L) at Eppegem in response to a given percentage increment on the total metal concentrations at possible CSO outlets.Table S3: Performance evaluation of MLR-and PCR-derived equations on representing the total metal concentration (MeT) and partitioning coefficient (Log based: LogK d ) at rural catchment runoff discharging points.Figure S1: The score plot (a) for the first two components (F1 and F2) and (b) for the third and fourth components (F3 and F4), showing with distinction of data between seven GESZ sampling campaigns.Figure S2: The Bayesian information criterion (BIC) scores corresponding to different independent variables used to derive the MLR based regression equation for total cadmium (MeT-Cd).Figure S3: The loading plot, showing relationships amongst different physio-chemical variables taken into consideration and the total metal concentrations with the corresponding distribution coefficients.Figure S4: Streamflow results at Eppegem for calibration (2007)(2008) and validation (2009-2010) at daily time scale.Figure S5: Total CSO discharge sent to river in the year of 2007 through systems of four collector systems in/around Brussels.Also, shown is the CSO volume as percentage of streamflow at Eppegem at daily time scale.Figure S6: Hourly and daily (averaged) stream flow at GESZ sampling station Z12.From simulation, the tidal effects are observed up to 9 km upstream of the outlet of Zenne river.

Figure 1 .
Figure 1.The Zenne basin with different sub-basins, river, canal and sewer network and locations of the Waste Water Treatment Plants (a); and locations of validation stations of different variables, possible Combined Sewer Overflow (CSO) points, sampling stations, meteorological stations and the river network (b).

Figure 2 .
Figure 2. Model components of the integrated trace metal transport model and the data exchange between the components in Open Modelling Interface (OpenMI).

Figure 3 .
Figure 3. Box plots of the total metal concentrations at the outlet of waste water treatment plants (WWTPs) (top row) and at Good Ecological Status of the river Zenne (GESZ) station Z13, near the river outlet (bottom row), measured during the GESZ sampling campaigns.

Figure 4 .
Figure 4. Predicted versus observed (a) total concentrations: MeT (b) Partitioning coefficients (log based: Log K d ) of Cd, Cu, Pb and Zn based on multi-linear regression (MLR) (circles) and principle component regression (PCR) (triangles).Dotted lines indicate the 1-standard deviation prediction interval.The solid lines correspond to the bisector lines.

Figure 5 .
Figure 5. Multi-linear regression (MLR) and principle component regression (PCR) based simulated (daily aggregated) and observed total (MeT), dissolved (MeD) and particulate (MeP) metal concentration for Cu, Cd, Zn and Pb for the time frame of 2007-2010.

Table 1 .
Concentration limits imposed by different legislation.
* Priority hazardous substance; ** Priority substance; † Depending on hardness as expressed in mgCaCO3/l; EU WFD: European Union Water Framework Directive; a EQS: Environmental Quality Standard, yearly average value EU [2]; Directive 2008/105/EC; b Quality of fresh waters to support fish life, yearly average value Directive 2006/44/EC EU [6]; ˆSalmonoid waters; ˆˆCyprinid waters; US EPA: United States Environmental Protection Agency; c Drinking water quality, based on instantaneous concentration US EPA

Table 2 .
Model inputs for the integrated model.

Table 3 .
The Good Ecological Status of the river Zenne (GESZ) sampling campaigns.

Table 5 .
Goodness-of-fit statistics and final performance rating of different model components of the integrated trace metal transport model.

Table 6 .
Simulated yearly average concentration of different metal species at Lot, Vilvoorde and Eppegem and comparison with Flemish and European Union (EU) mandatory levels.