This paper presents a new freeware simulation tool (IberWQ) for 2D water quality modelling in rivers and non-stratified estuaries. The model computes the spatial and temporal evolution of several species and variables which are relevant for the evaluation of the environmental status of rivers and estuaries, including: Escherichia coli, dissolved oxygen, carbonaceous biochemical oxygen demand, organic nitrogen, ammoniacal nitrogen, nitrate–nitrite nitrogen, water temperature and salinity. A depth-averaged transport equation is solved for each variable with a mass conservative unstructured finite volume solver. IberWQ is fully coupled to the hydrodynamic module of the software Iber, a freeware simulation tool for solving the 2D shallow water equations. Both models are integrated in the same windows graphical environment, including several tools which allow the user to generate unstructured meshes adapted to the site topography, define spatially variable input parameters and visualize model outputs. We present four application examples to illustrate the possibilities of the software for the dynamic simulation of water quality in rivers and estuaries.
SOFTWARE AVAILABILITY
- •
Name of software: IberWQ
- •
Developed by: Iberaula (http://www.iberaula.com/)
- •
First available year: 2015
- •
Software requirements: Windows operating system
- •
Hardware requirements: IBM-compatible PC
- •
Programming language: Fortran (numerical solver); C ++ (pre- and post-processing operations) and Tcl/Tk (GUI)
- •
Availability: IberWQ is integrated in the software package Iber, which can be downloaded for free at www.iberaula.com
NOTATION AND UNITS
- C
depth-averaged concentration (
)
- CBOD
depth-averaged concentration of carbonaceous biochemical oxygen demand (
)
- DO
depth-averaged concentration of dissolved oxygen (
)
depth-averaged saturation concentration of dissolved oxygen (
)
depth-integrated baroclinic force per unit mass in the x direction (
)
depth-integrated baroclinic force per unit mass in the y direction (
)
depth-integrated diffusive/dispersive mass flux in the x direction (kg/m/s)
depth-integrated diffusive/dispersive mass flux in the y direction (kg/m/s)
attenuation factor on CBOD degradation at low oxygen concentration (−)
attenuation factor on nitrification at low oxygen concentration (−)
increase factor on denitrification at low oxygen concentration (−)
- g
gravity acceleration (
)
- h
water depth (m)
reaeration rate (
)
CBOD degradation rate (
)
sediment oxygen demand rate (
)
organic nitrogen hydrolysis rate (
)
nitrification rate (
)
denitrification rate (
)
E. coli decay rate (
)
depth-averaged concentration of ammoniacal nitrogen (
)
depth-averaged concentration of nitrate-nitrite nitrogen (
)
depth-averaged concentration of organic nitrogen (
)
net incident radiation (short-wave and atmospheric radiation) (
)
outgoing long-wave back radiation flux (
)
heat loss by conduction to atmosphere (
)
heat loss by evaporation (
)
- S
depth-averaged concentration of salinity (
)
- T
depth-averaged water temperature (Celsius)
air temperature (Celsius)
depth-averaged water temperature (Kelvin)
depth-averaged velocity in the x direction (m/s)
depth-averaged velocity in the y direction (m/s)
modulus of the depth-averaged velocity (m/s)
wind velocity at 7 m above the water surface (m/s)
wind velocity at 10 m above the water surface (m/s)
settling velocity of CBOD (m/s)
settling velocity of org-N (m/s)
water radiation emissivity coefficient (−)
water density (
)
Stefan–Boltzmann constant (
)
atmospheric relative humidity (−)
INTRODUCTION
The European and national water-related legislations published in the last years, and more specifically the Water Framework Directive 2000/60/EC (EC 2000) and the Bathing Water Directive 2006/7/EC (EC 2006), require the assessment of the environmental, physical and chemical status of surface water bodies according to various criteria as temperature, oxygenation, nutrient conditions, faecal contamination and other specific pollutants. Environmental quality standards (EQS) have been proposed by the European Commission as a way to assess the status of water bodies by quantifying the presence of certain pollutants that pose a significant risk to the environment. The EQS specify maximum concentrations for specific water pollutants which should not be exceeded in order to protect human health and the environment. Specific implementation of these standards is the responsibility of local policy-makers and water resources managers. In this context, the Urban Pollution Management Manual (FWR 1998) defines concentration–duration–frequency thresholds which depend on the ecosystem and pollutant being analysed, and which should not be violated to reach the acceptable ecological status of the receiving water bodies.
Numerical models capable of simulating the concentration of pollutants in the aquatic environment are valuable and increasingly used tools for better implementation of the water quality regulations. In recent years, several codes have implemented 1D, 2D and 3D water quality modules which allow decision-makers to evaluate alternative measures for the control of sewage spills by modelling the temporal and spatial evolution of pollutants and other environmental variables in the receiving waters.
This paper presents IberWQ, a new software tool for 2D water quality modelling in rivers and non-stratified estuaries which computes the spatial and temporal evolution of several species and variables including: Escherichia coli, dissolved oxygen (DO), carbonaceous biochemical oxygen demand (CBOD), organic nitrogen (org-N), ammoniacal nitrogen (-N), nitrate–nitrite nitrogen (
-N), water temperature (T) and salinity (S). The present version of the model does not simulate phosphorus, which might be a water quality limiting factor in some fresh surface water bodies. IberWQ is integrated in the freeware package Iber (www.iberaula.com), which includes a hydrodynamic module for solving the 2D shallow water equations (Bladé et al. 2014b). The structure of the model and several examples showing its computational capabilities are presented.
SOFTWARE DESCRIPTION
Software structure
IberWQ is a new water quality model integrated in the freeware package Iber (www.iberaula.com). This software package includes a hydrodynamic module which solves the 2D shallow water equations (Bladé et al. 2014b), and a sediment transport module which solves the 2D Exner equation including bed and suspended load. All the modules share the same pre- and post-processing windows interface, which is based on the software GiD (www.gidhome.com) and includes several tools for the automatic generation of unstructured meshes adapted to the site topography, the definition of spatially variable input parameters and the visualization of results in different formats.
Integration of the module IberWQ in the software package Iber, and input data needed to run each module.
Integration of the module IberWQ in the software package Iber, and input data needed to run each module.
Structure of IberWQ. The variables inside black boxes are computed by the model. Solid arrows indicate interactions between variables. Dashed lines indicate dependence of reaction kinetic constants on model variables.
Structure of IberWQ. The variables inside black boxes are computed by the model. Solid arrows indicate interactions between variables. Dashed lines indicate dependence of reaction kinetic constants on model variables.
The structure of the water quality module and the biochemical reactions between the different species considered are shown in Figure 2. Most of the reaction kinetics depend on the water temperature and salinity, which might be computed by the model or introduced as input data by the user. In the case that the temperature field is computed by the model, additional atmospheric input data (time series of net solar and atmospheric radiation, air temperature, atmospheric relative humidity and wind velocity) must be introduced by the user. The different species considered in the model can be activated or deactivated according to the specific problem being solved. The only restriction is that the three nitrogen species must be always solved together. The concentration of the species which are not activated is assumed to be zero except for the DO, which in the case of not being computed is assumed to be equal to its saturation value, and the water temperature, which is set to 20 °C if its value is not specified by the user.
Hydrodynamic equations
The shallow water equations are solved using the finite volume method and the numerical scheme of Roe (1986). The reader is referred to Bladé et al. (2014b) and the references therein for a detailed description and experimental validation of the numerical schemes used to solve the shallow water equations, which are not detailed here for the sake of conciseness. The hydraulic module has been applied to several hydraulic studies in the past, including river inundation modelling (Bladé et al. 2014a), overland flow (Cea & Bladé 2015), evaluation of gully restoration measures (Castillo et al. 2014) and wood transport in rivers (Ruiz-Villanueva et al. 2014).
Water quality equations
is a generic source/sink term which depends on the species considered, h is the water depth and are the two components of the depth-averaged velocity. The terms are the two components of the depth-integrated diffusive/dispersive mass flux. If wave field data are entered by the user, a non-isotropic dispersion coefficient due to short waves can be computed automatically using the formulation proposed by Law (2000), which depends on wave height, period and direction of propagation. Implementation details can be found in Cea et al. (2011). Turbulent diffusion is assumed to be isotropic in the two horizontal directions and proportional to the eddy viscosity coefficient computed by a turbulence model. The user can choose between several depth-averaged turbulence models, including a mixing-length model (Cea et al. 2007) and the model of Rastogi & Rodi (1978).
are the kinetic constants for the processes depicted in Figure 2, and are, respectively, the settling velocities of CBOD and organic nitrogen, T is the water temperature in Celsius, is the concentration of DO at saturation, is the net incoming radiation through the water surface (short-wave solar radiation and long-wave atmospheric radiation), is the outgoing long-wave back radiation flux, is the heat loss by conduction from the water surface to the atmosphere and is the heat loss by evaporation. Notice that, since salinity is considered a conservative variable, no source terms are included in its transport equation.
The concentration of DO at saturation () is computed from the formulation presented in Eaton et al. (2005), which depends on the water temperature, salinity and altitude above sea level. Anoxia is considered in the biodegradation and nitrification reactions by reducing the oxidation rates at low oxygen levels (Chapra et al. 2008) by means of the reduction factors
and
, while the function
accounts for the increase of denitrification rate at low oxygen concentrations (Chapra et al. 2008).
Regarding the degradation of E. coli, the user might introduce a decay rate calibrated from field data ( in Figure 2 and in Equation (8)), or choose between the standard formulations of Canteras et al. (1995) and Mancini (1978), which take into account the effects of water temperature, salinity and solar radiation in the E. coli decay rate. The latter is considered one of the most complete models of the faecal coliform decay process (Manache et al. 2007).
Numerical solver
The convection–diffusion equation for each species is solved in an unstructured finite volume mesh formed by triangular and quadrilateral elements. In order to solve Equation (2), the water depth (h) and water fluxes ( and
) are obtained from the hydrodynamic module of Iber. The most important characteristic of the numerical scheme used to solve Equation (2) is that discretization of the advective flux is coupled to the mass flux discretization in the shallow water equations. This is a necessary requirement for a mass conservative discretization of a depth-averaged scalar transport equation (Audusse & Bristeau 2003; Benkhaldoun et al. 2007; Murillo et al. 2008; Cea & Vázquez-Cendón 2012). Specific details of the numerical implementation in IberWQ can be found in Cea & Vázquez-Cendón (2012). The most relevant implication is that the discretization of the transport equations solved in IberWQ is mass conservative.
A first order upwind discretization (Versteeg & Malalasekera 1995), as well as the second order Gamma scheme proposed in Jasak et al. (1999), is implemented in the solver. The Gamma scheme uses second order central differencing wherever this scheme satisfies a boundedness criterion in the normalized variable diagram (Leonard 1988), and first order upwind differencing otherwise. In terms of numerical accuracy, the Gamma scheme introduces less numerical diffusion in the solution than the first order scheme, which implies a sharper resolution of the spatial gradients of concentration. Using the Gamma scheme implies an increase on the CPU time of approximately 5%.
The CPU time required to solve the transport equation for each species is, depending on the species considered, between 10% and 15% of the CPU time required by the hydrodynamic module to solve the 2D shallow water equations. This means that the CPU time required to solve the eight transport equations considered in IberWQ is approximately the same as the time required to solve the shallow water equations (i.e., the total CPU time doubles when solving the complete water quality model). On the other hand, if only a few species are considered, the percentage increase in CPU time can be estimated between 0.10 and 0.15 times the number of transport equations being solved.
User interface
Unstructured TIN mesh in the pre-processing interface (left) and image of E. coli concentration results in the post-processing interface (right).
Unstructured TIN mesh in the pre-processing interface (left) and image of E. coli concentration results in the post-processing interface (right).
MODEL VALIDATION AND APPLICATION
This section presents four application examples of IberWQ which show some of the model's capabilities and limitations. In the first example the model results are compared with an analytical solution which considers the transformation of the three species of nitrogen in a simple 1D channel geometry. The second example considers organic matter (CBOD) and DO in addition to nitrogen. In this case the results of the model are compared to field measurements presented in Thomann & Mueller (1987). The third example shows the application of the model to the simulation of E. coli concentration in a large estuary. Model results are compared to field measurements, showing the predictive capabilities that can be expected in this kind of application. In the last application example the recovering rate of DO concentration is modelled in a river reach located downstream of a reservoir.
The original data files of the examples can be downloaded from www.iberaula.es and used directly in IberWQ, allowing users to modify the input data at their convenience and verify the model results.
First example: discharge of organic nitrogen in a 1D channel
The purpose of the first example is to compare the results of the model with an analytical solution in a 1D channel, including the different transformations of nitrogen considered in IberWQ (hydrolysis, nitrification and denitrification). The test case consists of a discharge of organic nitrogen in a 1,200 m channel with a mean longitudinal velocity of 1 m/s and a dispersion coefficient of 5 m2/s. A 500-second pulse type injection of organic nitrogen with a concentration of 1 mg/L is imposed at the upstream boundary. The three-member decay chain of nitrogen is considered (organic nitrogen > ammonia > nitrates), with values of 0.004, 0.001 and 0.002 s−1 for the first order kinetic coefficients of hydrolysis, nitrification and denitrification (,
and
in Figure 2). In order to validate the numerical results with an analytical solution, the organic nitrogen settling velocity is neglected and the correction factors that incorporate the influence of oxygen concentration on the nitrification and denitrification rates are fixed to
in Equations (6) and (7). The analytical solution for this problem can be found in Genuchten et al. (2013).
The 1D channel was discretized in IberWQ with a longitudinal mesh made of quadrilaterals. Since the problem is 1D, only one element was used in the transversal direction. In the longitudinal direction, three different mesh sizes were tested, with mesh spacings equal to 2.5 m, 5 m and 10 m.
First example: MAE computed from the concentration profiles 1,000 s after the beginning of the pulse for different mesh sizes (Δx) and numerical schemes
| MAE (mg/L) | ||||
|---|---|---|---|---|
| Δx (m) | Numerical scheme | org-N | NH3-N | NH3-N |
| 10 | First order | 0.0052 | 0.037 | 0.0145 |
| 5 | First order | 0.0033 | 0.0242 | 0.0095 |
| 2.5 | First order | 0.0025 | 0.0182 | 0.0070 |
| 10 | Gamma | 0.0047 | 0.0349 | 0.0134 |
| 5 | Gamma | 0.0031 | 0.0233 | 0.0089 |
| 2.5 | Gamma | 0.0023 | 0.0179 | 0.0069 |
| MAE (mg/L) | ||||
|---|---|---|---|---|
| Δx (m) | Numerical scheme | org-N | NH3-N | NH3-N |
| 10 | First order | 0.0052 | 0.037 | 0.0145 |
| 5 | First order | 0.0033 | 0.0242 | 0.0095 |
| 2.5 | First order | 0.0025 | 0.0182 | 0.0070 |
| 10 | Gamma | 0.0047 | 0.0349 | 0.0134 |
| 5 | Gamma | 0.0031 | 0.0233 | 0.0089 |
| 2.5 | Gamma | 0.0023 | 0.0179 | 0.0069 |
First example. Analytical and numerical concentration profiles for the three nitrogen species 1,000 s after the beginning of the discharge (left) and at the location x = 400 m (right). The numerical results were computed with longitudinal mesh sizes of 2.5 m and 10 m, and with the second order Gamma scheme.
First example. Analytical and numerical concentration profiles for the three nitrogen species 1,000 s after the beginning of the discharge (left) and at the location x = 400 m (right). The numerical results were computed with longitudinal mesh sizes of 2.5 m and 10 m, and with the second order Gamma scheme.
Second example: wastewater discharge of CBOD and ammonia in a river
, -N and -N downstream the discharge presented in Thomann & Mueller (1987) provide a benchmark comparison for the numerical results obtained with IberWQ. The spatial domain covers a total distance of 50 km, and was discretized with a longitudinal mesh size of 25 m. Since the length of the river reach is three orders of magnitude larger than the river width, the computational mesh has a single element in the transverse direction. A constant river slope of 0.0001 and a Manning coefficient of 0.10 were used in order to satisfy the depth and velocity values given in Thomann & Mueller (1987). The water quality variables considered in the model are DO, CBOD, ammonia and nitrates. In order to compare the computed values of total CBOD with the field data of , a constant ratio of CBOD/ = 2 was considered. The first order kinetic coefficients for organic matter degradation and nitrification were obtained from Thomann & Mueller (1987) with the following values: and . The discharge of CBOD and ammonia induces a decay in the concentration of DO along the first 10 km downstream the discharge (Figure 5) due to the oxygen demand of the processes of biodegradation and nitrification. Further downstream the biodegradation and nitrification rates diminish due to the descent on the concentration of CBOD and ammonia, and the concentration of DO increases progressively by surface reaeration, but it takes more than 30 km to recover the original DO levels of the river. The numerical model reproduces correctly the observed spatial evolution of the four species considered (Figure 5).Second example. Definition of the case (from Thomann & Mueller (1987)) and longitudinal concentration profiles of DO, CBOD, NH3-N and NO3-N.
Second example. Definition of the case (from Thomann & Mueller (1987)) and longitudinal concentration profiles of DO, CBOD, NH3-N and NO3-N.
Third example: faecal contamination in Carmarthen Bay, UK
Faecal contamination in Carmarthen Bay. Spatial domain included in the numerical model and location of control points and STW discharges. The colours represent the depth of the seabed below mean sea level.
Faecal contamination in Carmarthen Bay. Spatial domain included in the numerical model and location of control points and STW discharges. The colours represent the depth of the seabed below mean sea level.
Faecal contamination in Carmarthen Bay. Comparison of measured and predicted E. coli concentrations and evolution of computed values at control points 1 (a) and (b) and 2 (c) and (d). is the time in hours required to reach a 90% reduction in E. coli concentration, and is related to the decay rate as .
Faecal contamination in Carmarthen Bay. Comparison of measured and predicted E. coli concentrations and evolution of computed values at control points 1 (a) and (b) and 2 (c) and (d). is the time in hours required to reach a 90% reduction in E. coli concentration, and is related to the decay rate as .
Results from a field monitoring campaign were used to evaluate the performance of the model in predicting the E. coli concentration (Figure 7). The comparison between simulated and observed E. coli levels at two control points in the estuary shows that the model captures the order of magnitude of the observed concentrations, the maximum values reaching 100 cfu/100 mL at the first control point and 10 cfu/100 mL at the second. Nonetheless, the exact details of the temporal dynamics are not properly captured. While the simulated concentrations of E. coli show a repeated pattern over the tidal cycle, the measured concentrations exhibit a more random behaviour, which is typically observed in marine environments (Quilliam et al. 2011; Bedri et al. 2013). The high levels of uncertainty involved in the numerical modelling and field measurement of pathogens makes it difficult to obtain accurate predictions of the temporal patterns of E. coli concentration in coastal waters, even with more complex 3D models (Chan et al. 2013; Bedri et al. 2014). Even so, 2D models like IberWQ can serve as a useful tool to determine the factors controlling the E. coli concentrations and to evaluate the impact of water management strategies in non-stratified systems (De Brauwere et al. 2011; Chan et al. 2013; Boye et al. 2015).
Fourth example: DO levels downstream of a reservoir
In this hypothetical example, IberWQ is used to analyse the influence of river discharge on the DO levels in a 10-km-long river reach located downstream of an hydroelectric power station. The turbine intakes draw water from the deeper level of the reservoir, characterized by low DO levels. Two different upstream boundary conditions were modelled. The first one corresponds to a peak discharge of 120 m3/s and the second one to a minimum environmental flow of 5 m3/s. At the downstream boundary the normal depth at the boundary cross section was imposed. The DO concentration in the water released through the dam during generation was set to 3 mg/L in both cases. Organic matter (CBOD) and nitrogen concentration are assumed to be zero in this example. Water temperature and salinity, which are needed to compute the saturated DO concentration, were fixed to 20 °C and zero, respectively. The only species considered in the model is therefore DO. The DO concentration does not need to be specified at outlet boundaries and was therefore not imposed at the downstream cross section. The river reach was discretized with a structured mesh of 10,428 elements, with an average element size of 22 m2.
Fourth example. DO levels downstream of the reservoir. Computed DO concentrations for minimum (left) and peak (right) discharges.
Fourth example. DO levels downstream of the reservoir. Computed DO concentrations for minimum (left) and peak (right) discharges.
CONCLUSIONS
IberWQ is a 2D depth-averaged water quality module which allows the user to simulate different scenarios regarding the mixing, transport and degradation of several species related to the environmental status of rivers and estuaries. The numerical solver, as well as the graphical user interface, is fully integrated in the freeware package Iber, and coupled to its hydrodynamic module, which solves the 2D shallow water equations. In this context, IberWQ is a useful numerical tool to assist the environmental management of surface water bodies via the application of EQS as those required in European water related directives. The windows pre- and post-processing environment allows the user to set up a simulation and analyse the spatial and temporal evolution of salinity, temperature, DO, CBOD, organic nitrogen, ammoniacal nitrogen, nitrate–nitrite nitrogen and E. coli. The CPU time required to run the water quality model varies between 10% and 100% of the CPU time required by the hydrodynamic module to solve the Saint Venant equations, depending on the number of species considered. Future versions of the model will include other species of interest to define the environmental status of rivers and estuaries, such as phosphorous, which might be a water quality limiting factor in fresh surface water bodies.
ACKNOWLEDGEMENTS
The author Maria Bermudez gratefully acknowledges financial support from the Spanish Regional Government of Galicia, Plan I2C (postdoctoral grant reference ED481B 2014/156-0).
%20cropped.png?versionId=7142)