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.

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

*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

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

*et al.*1985; Chapra 1997; Kannel

*et al.*2011). The expressions of the source/sink terms for each species are the following:

*et al.*(1985) and Eaton

*et al.*(2005). On the other hand, the reaeration constant is automatically evaluated as: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: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: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).

*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

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

## 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 m^{2}/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.

. | . | MAE (mg/L) . | ||
---|---|---|---|---|

Δx (m)
. | Numerical scheme . | org-N . | NH_{3}-N
. | NH_{3}-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 . | NH_{3}-N
. | NH_{3}-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 |

### Second example: wastewater discharge of CBOD and ammonia in a river

### Third example: faecal contamination in Carmarthen Bay, UK

*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 m

^{3}/s and 29–38,000 cfu/100 mL, respectively, while the STW discharges varied within 0.004 and 0.36 m

^{3}/s, with

*E. coli*concentrations up to 34 × 10

^{6}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.

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 m^{3}/s and the second one to a minimum environmental flow of 5 m^{3}/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 m^{2}.

^{3}/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 m

^{3}/s.

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