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.

Figure 1 shows the input data needed by IberWQ and its relations with the other modules of the Iber package. IberWQ relies on the results of water depth, depth-averaged velocity and turbulent viscosity of the hydrodynamic module in order to solve the convection–diffusion equation for each species. Since all modules are integrated in the same windows graphical interface and they share the same unstructured finite volume mesh to solve their respective transport equations, the coupling between modules is handled automatically.
Figure 1

Integration of the module IberWQ in the software package Iber, and input data needed to run each module.

Figure 1

Integration of the module IberWQ in the software package Iber, and input data needed to run each module.

The water quality parameters which need to be defined by the user are the reaction constants and settling velocities defined in Figure 2. In addition, IberWQ allows the user to introduce the radiation stresses due to wave forcing in coastal areas as an external input field obtained from a wave propagation model, as for example the SWAN model (Holthuijsen 2007).
Figure 2

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.

Figure 2

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 numerical solver of IberWQ is coupled to the hydrodynamic module of Iber, which solves the 2D shallow water equations written in conservative form as: 
formula
1
 
formula
 
formula
where is the bed elevation, h is the water depth, are the two components of the unit discharge, is the modulus of the depth-averaged velocity, are the two components of the depth-averaged velocity, n is the Manning coefficient, g is the gravity acceleration and is the turbulent viscosity, which is computed from a depth-averaged turbulence model (Cea et al. 2007).

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

The following generic convection–diffusion equation is solved for each of the species considered in the model: 
formula
2
where C is the depth-averaged concentration of the species, 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).
The physical processes and biochemical reactions considered for each species are modelled by means of the source term in Equation (2), and represented in Figure 2. All the interactions between species are modelled as first order reactions (Bowie et al. 1985; Chapra 1997; Kannel et al. 2011). The expressions of the source/sink terms for each species are the following: 
formula
3
 
formula
4
 
formula
5
 
formula
6
 
formula
7
 
formula
8
 
formula
9
with: 
formula
10
where 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 reaction constants - and the settling velocities - must be introduced by the user as model parameters. Typical ranges of variation for these constants are given in Bowie et al. (1985) and Eaton et al. (2005). On the other hand, the reaeration constant is automatically evaluated as: 
formula
11
where accounts for the effect of water depth and velocity on reaeration and for the effect of wind velocity. The value of is computed as detailed in Covar (1976) as: 
formula
12
where h is the water depth in m and is the modulus of the depth-averaged water velocity in m/s. The effect of wind velocity on is considered in , which is evaluated following Banks & Herrera (1977) as: 
formula
13
where is the wind velocity in m/s at 10 m above the water surface level.

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).

In the water temperature equation the time series of net incident radiation () must be introduced by the user, while the other heat losses components (, and ) are computed as (Chapra 1997): 
formula
14
 
formula
15
 
formula
16
where all the heat losses are expressed in , is the water radiation emissivity coefficient, is the Stefan–Boltzmann constant, is the wind velocity at 7 m above the water surface, is the air temperature in Celsius and is the atmospheric relative humidity.
By default, the coupling between the hydrodynamics module and IberWQ is just in one direction, i.e., the flow field is necessary to solve the species transport equations, but it is not influenced by the results of the water quality variables. However, if the fields of salinity or temperature are computed, it is possible to include in the hydrodynamic equations the effects of baroclinic pressure, in order to account for the horizontal gradients of water density. In this case, the depth-integrated baroclinic force terms per unit mass are evaluated as: 
formula
17
where is the water density, which is computed from the salinity and temperature fields using the international one atmosphere equation of state of seawater (Millero & Poisson 1981): 
formula
18
with: 
formula
 
formula
 
formula
 
formula
where is the water temperature in degrees Kelvin, and S the water salinity in . The depth-integrated baroclinic force terms are introduced, respectively, in the x- and y-momentum conservation equations of Equation (1).

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

IberWQ is integrated in the pre- and post-processing windows graphical environment of the software package Iber. The pre-processing tools allow the user to build an unstructured mesh adapted to the site topography. Several meshing algorithms are available, including triangular irregular network (TIN), right-triangulated irregular network (Evans et al. 2001) and structured meshes (Figure 3). The user can control the mesh size manually or automatically as a function of the topography or the location of the sewage discharges. The post-processing interface allows the user to visualize the results as concentration fields, time series or longitudinal profiles, and to export them in various formats as images and graphical videos, text files or georeferenced rasters. Specific tools to operate with the variables and to perform spatial integration are also integrated in the post-processing interface.
Figure 3

Unstructured TIN mesh in the pre-processing interface (left) and image of E. coli concentration results in the post-processing interface (right).

Figure 3

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.

The spatial and temporal evolution of the three species considered, as well as the maximum concentration levels, are correctly predicted by the model, although there is a certain lag between the analytical and numerical solutions (Figure 4). This lag is due to the numerical errors introduced by the discretization of the advective flux, which are reduced as the mesh size is refined. The mean absolute error (MAE) obtained with three different mesh sizes in combination with the first order and the Gamma schemes is shown in Table 1. Although the Gamma scheme reduces the MAE with respect to the first order scheme, in this case the effect of the mesh size clearly dominates the error on model output.
Table 1

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 
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 
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 
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 
Gamma 0.0031 0.0233 0.0089 
2.5 Gamma 0.0023 0.0179 0.0069 
Figure 4

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.

Figure 4

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

This example analyses the evolution of CBOD, ammonia, nitrates and DO concentrations along a river reach. It is based on a sample problem from Thomann & Mueller (1987), which is summarized in Figure 5. A wastewater treatment plant discharges CBOD and ammonia into a river, affecting its distribution of DO. The observed in-stream concentrations of DO, , -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).
Figure 5

Second example. Definition of the case (from Thomann & Mueller (1987)) and longitudinal concentration profiles of DO, CBOD, NH3-N and NO3-N.

Figure 5

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

In this example, IberWQ is used to model faecal coliform contamination in Carmarthen Bay, an extensive shallow inlet that encompasses the estuaries of the river Loughor, Taf, Tywi and Gwendraeth, in the north-western boundary of the Bristol Channel (Figure 6). The model domain, which covers the main part of the Bristol Channel and a large part of the Severn Estuary, was discretized using an unstructured triangular mesh of approximately 80,000 elements, with finer mesh elements near the location of the sewage discharges.
Figure 6

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.

Figure 6

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.

Available field discharge and concentration data, obtained as described in Schnauder et al. (2007), were used to define the river and sewage treatment works (STW) inputs to the model. A total of 59 locations of river and stream catchment outlets and 15 STWs were considered. Inflows and E. coli concentrations measured in the rivers Taf, Tywi, Gwendraeth Fawr, Gwendraeth Fach and Loughor were introduced in the model as inlet boundary conditions. The remaining stream and STWs inputs were included in the model as point discharges within the domain. The discharge and E. coli concentration for each river input varied within 0.01–14 m3/s and 29–38,000 cfu/100 mL, respectively, while the STW discharges varied within 0.004 and 0.36 m3/s, with E. coli concentrations up to 34 × 106 cfu/100 mL. At the open sea boundary, the model was forced with a time series of tidal water levels and with a zero concentration of E. coli. The formulation proposed by Mancini (1978) was used to compute the bacteria decay rate, based on measured temperature, salinity, turbidity and solar radiation levels. The model thus calculates a spatially and temporally variable decay rate, as shown in Figure 7. For a precise definition of the time varying conditions, the IberWQ data files are available to download at www.iberaula.es.
Figure 7

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 .

Figure 7

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.

Figure 8 shows the predicted DO levels in the river reach for the maximum and minimum discharge rates. DO concentration increases gradually downstream of the dam due to surface reaeration. For the 5 m3/s discharge, oxygen concentrations at the end of the reach are close to saturation, exceeding 8 mg/L. In contrast, DO values remain below 5 mg/L for the whole reach during the peak discharge of 120 m3/s.
Figure 8

Fourth example. DO levels downstream of the reservoir. Computed DO concentrations for minimum (left) and peak (right) discharges.

Figure 8

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).

REFERENCES

REFERENCES
Audusse
E.
Bristeau
M. O.
2003
Transport of pollutant in shallow water: a two time steps kinetic method
.
Mathematical Modelling and Numerical Analysis
37
(
2
),
389
416
.
Banks
R. B.
Herrera
F. F.
1977
Effect of wind and rain on surface reaeration
.
Journal of the Environmental Engineering Division
103
(
3
),
489
504
.
Bedri
Z.
Corkery
A.
O'Sullivan
J. J.
Alvarez
M. X.
Erichsen
A. C.
Deering
L. A.
Demeter
K.
O'Hare
G. M.
Meijer
W. G.
Masterson
B.
2014
An integrated catchment-coastal modelling system for real-time water quality forecasts
.
Environmental Modelling & Software
61
,
458
476
.
Bladé
E.
Cea
L.
Corestein
G.
2014a
Modelización numérica de inundaciones fluviales
[Numerical modelling of river inundations]
.
Ingeniería del Agua
18
(
1
),
71
82
.
Bladé
E.
Cea
L.
Corestein
G.
Escolano
E.
Puertas
J.
Vázquez-Cendón
M. E.
Dolz
J.
Coll
A.
2014b
Iber: herramienta de simulación numérica del flujo en ríos
[Iber: tool numerical simulation of flow in rivers]
.
Revista Internacional de Métodos Numéricos para Cálculo y Diseño en Ingeniería
30
(
1
),
1
10
.
Bowie
G. L.
Mills
W. B.
Porcella
D. B.
Campbell
C. L.
Pagenkopf
J. R.
Rupp
G. L.
Johnson
K. M.
Chan
P. W. H.
Gherini
S. A.
Chamberlin
C. E.
1985
Rates, Constants, and Kinetic Formulations in Surface Water Quality Modeling
.
US Environmental Protection Agency, ORD
,
Athens, GA
,
ERL, EPA/600/3-85/040
.
Boye
B. A.
Falconer
R. A.
Akande
K.
2015
Integrated water quality modelling: application to the Ribble Basin, UK
.
Journal of Hydroenvironment Research
9
,
187
199
.
Castillo
C.
Pérez Alcántara
R.
Gómez Calero
J. A.
2014
A conceptual model of check dam hydraulics for gully control: efficiency, optimal spacing and relation with step-pools
.
Hydrology and Earth System Sciences
18
,
1705
1721
.
Cea
L.
Puertas
J.
Vázquez-Cendón
M. E.
2007
Depth-averaged modelling of turbulent shallow water flow with wet-dry fronts
.
Archives of Computational Methods in Engineering (ARCME)
14
(
3
),
303
341
.
Chapra
S. C.
1997
Surface Water-Quality Modeling
.
Waveland Press
,
Long Grove, IL
,
USA
.
Chapra
S. C.
Pelletier
G. J.
Tao
H.
2008
QUAL2K: A Modeling Framework for Simulating River and Stream Water Quality, Version 2.11: Documentation and Users Manual
.
Civil and Environmental Engineering Department, Tufts University
,
Medford, MA
,
USA
.
Covar
A. P.
1976
Selecting the proper reaeration coefficient for use in water quality models
. In:
U.S. EPA Conference on Environmental Simulation and Modeling
.
Cincinnati, OH
,
USA
.
De Brauwere
A.
De Brye
B.
Servais
P.
Passerat
J.
Deleersnijder
E.
2011
Modelling Escherichia coli concentrations in the tidal Scheldt river and estuary
.
Water Research
45
(
9
),
2724
2738
.
Eaton
A. D.
,
APHA, AWWA & WEF
2005
Standard Methods for the Examination of Water and Wastewater
.
American Public Health Association (APHA)
,
Washington, DC
,
USA
.
EC
2000
Directive 2000/60/EC of the European Parliament and of the Council of 23 October 2000 establishing a framework for Community action in the field of water policy. Official Journal L327
.
EC
2006
Council Directive 2006/7/EC of 15 February 2006 Concerning the Management of Bathing Water Quality and Repelling Directive 76/160/EEC. Official Journal L064
.
Evans
W.
Kirkpatrick
D.
Townsend
G.
2001
Right-triangulated irregular networks
.
Algorithmica
30
(
2
),
264
286
.
FWR
1998
Urban Pollution Management Manual: a Planing Guide for the Management of Urban Wastewater Discharges During Wet Weather
.
Foundation of Water Research
,
Marlow
,
UK
.
Genuchten
M. T. V.
Leij
F. J.
Skaggs
T. H.
Toride
N.
Bradford
S. A.
Pontedeiro
E. M.
2013
Exact analytical solutions for contaminant transport in rivers 2. Transient storage and decay chain solutions
.
Journal of Hydrology and Hydromechanics
61
(
3
),
250
259
.
Holthuijsen
L. H.
2007
Waves in Oceanic and Coastal Waters
.
Cambridge University Press
,
Cambridge
,
UK
.
Jasak
H.
Weller
H. G.
Gosman
A. D.
1999
High resolution NVD differencing scheme for arbitrarily unstructured meshes
.
International Journal for Numerical Methods in Fluids
31
,
431
449
.
Kannel
P. R.
Kanel
S. R.
Lee
S.
Lee
Y.-S.
Gan
T. Y.
2011
A review of public domain water quality models for simulating dissolved oxygen in rivers and streams
.
Environmental Modeling & Assessment
16
(
2
),
183
204
.
Law
A. W. K.
2000
Taylor dispersion of contaminants due to surface waves
.
Journal of Hydraulic Research
38
(
1
),
41
48
.
Leonard
B. P.
1988
Simple high-accuracy resolution program for convective modelling of discontinuities
.
International Journal for Numerical Methods in Fluids
8
,
1291
1318
.
Manache
G.
Melching
C. S.
Lanyon
R.
2007
Calibration of a continuous simulation fecal coliform model based on historical data analysis
.
Journal of Environmental Engineering
133
(
7
).
681
691
.
Mancini
J. L.
1978
Numerical estimates of coliform mortality rates under various conditions
.
Journal (Water Pollution Control Federation)
50
(
11
),
2477
2484
.
Millero
F. J.
Poisson
A.
1981
International one-atmosphere equation of state of seawater
.
Deep Sea Research Part A. Oceanographic Research Papers
28
(
6
),
625
629
.
Murillo
J.
Burguete
J.
Garcia-Navarro
P.
2008
Analysis of a second order upwind method for the simulation of solute transport in 2D shallow water flow
.
International Journal of Numerical Methods in Fluids
56
,
661
686
.
Quilliam
R.
Clements
K.
Duce
C.
Cottrill
S.
Malham
S.
Jones
D.
2011
Spatial variation of waterborne Escherichia coli – implications for routine water quality monitoring
.
Journal of Water and Health
9
(
4
),
734
737
.
Rastogi
A. K.
Rodi
W.
1978
Predictions of heat and mass transfer in open channels
.
Journal of the Hydraulics Division
(
HY3
),
104
(
3
),
397
420
.
Roe
P. L.
1986
A basis for the upwind differencing of the two-dimensional unsteady Euler equations
.
Numerical Methods for Fluid Dynamics
2
,
55
80
.
Ruiz-Villanueva
V.
Bladé-Castellet
E.
Sánchez-Juny
M.
Martí
B.
Díez Herrero
A.
Bodoque
J. M.
2014
Two dimensional numerical modelling of wood transport
.
Journal of Hydroinformatics
16
(
5
),
1077
1096
.
Schnauder
I.
Bockelmann-Evans
B.
Lin
B.
2007
Modelling faecal bacteria pathways in receiving waters
.
Proceedings of the ICE-Maritime Engineering
160
(
4
),
143
153
.
Thomann
R. V.
Mueller
J. A.
1987
Principles of Surface Water Quality Modeling and Control
.
Harper & Row
,
New York
,
USA
.
Versteeg
H. K.
Malalasekera
W.
1995
An Introduction to Computational Fluid Ddynamics. The Finite Volume Method
.
Addison Wesley Longman Ltd
,
Harlow
,
UK
.