Next Article in Journal
Variability of the Initial Abstraction Ratio in an Urban and an Agroforested Catchment
Next Article in Special Issue
Investigation of Morphological Changes in the Tamsui River Estuary Using an Integrated Coastal and Estuarine Processes Model
Previous Article in Journal
Rethinking the Framework of Smart Water System: A Review
Previous Article in Special Issue
Comparative Analysis of Four Baseflow Separation Methods in the South Atlantic-Gulf Region of the U.S.
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

IberWQ: A GPU Accelerated Tool for 2D Water Quality Modeling in Rivers and Estuaries

1
Environmental Physics Laboratory (EPhysLab), CIM-UVIGO, Universidade de Vigo, Campus As Lagoas s/n, 32004 Ourense, Spain
2
Environmental and Water Engineering Group, Department of Civil Engineering, Universidade da Coruña, Elviña, 15071 A Coruña, Spain
*
Author to whom correspondence should be addressed.
Water 2020, 12(2), 413; https://doi.org/10.3390/w12020413
Submission received: 18 December 2019 / Revised: 30 January 2020 / Accepted: 31 January 2020 / Published: 4 February 2020

Abstract

:
Numerical models are useful tools to analyze water quality by computing the concentration of physical, chemical and biological parameters. The present work introduces a two-dimensional depth-averaged model that computes the most relevant and frequent parameters used to evaluate water quality. High performance computing (HPC) techniques based on graphic processing unit (GPU) parallelization have been applied to improve the efficiency of the package, providing speed-ups of two orders of magnitude in a standard PC. Several test cases were analyzed to show the capabilities and efficiency of the model to evaluate the environmental status of rivers and non-stratified estuaries. IberWQ will be freely available through the package Iber.

1. Introduction

The application of Environmental Quality Standards (EQS) has become a widely used technique to assess the status of surface water bodies. These standards are based on the spatial and temporal evolution of certain variables and pollutants that might pose a significant risk to the environment (temperature, salinity, fecal contamination, nutrients, dissolved oxygen, pH, etc.). The concentration of these variables, which depends on the receiving waters, must be within a range of values to protect the environment and the human health. More generally, not only the exceedance of a given threshold should be controlled, but also the frequency and duration of those exceedance events. Concentration-duration-frequency curves can be used to define the thresholds that should not be exceeded to guarantee a correct environmental status [1,2,3].
In this context, numerical models have become valuable tools to help decision makers to evaluate alternative measures and solutions for the control of sewage spills. Some of the most known and used models for this purpose are the Water Quality Analysis Simulation Program (WASP) originally developed by [4,5], the QUAL2K model [6] and the CE-QUAL-W2 model [7,8]. Since these models assume different hydrodynamic approximations, their suitability is strongly dependent on the case under study. QUAL2K is a 1D steady flow model for rivers that assumes that the flow is well-mixed in the whole cross section. CE-QUAL-W2 is a 2D laterally averaged model that solves the hydrodynamics in the longitudinal and vertical direction, being therefore appropriate for rivers, reservoirs and estuaries well-mixed in the lateral direction. WASP can solve 1D, 2D and 3D problems with a large variety of pollutant types. However, this model needs to be linked with a hydrodynamic model to provide accurate flow velocities and water elevation. Many other codes have been used in the scientific literature to model different water quality parameters in rivers and estuaries [9,10,11,12,13].
The model IberWQ was presented in [14] as a numerical tool for 2D water quality modeling in rivers and non-stratified estuaries. Its use is especially adequate to simulate the mixing of effluents when the flow is well mixed in the vertical direction. This model is integrated into the software package Iber and is freely available at http://www.iberaula.com.
IberWQ is not a suitable tool to deal with flows with vertical stratification or with relevant vertical gradients of velocity and concentration. In this context, it might be used as a complementary tool to 3D near-field models as VisualPLUMES [15] or CORMIX [16]. It could also be integrated in a hydrological-hydraulic modeling cascade approach [17] as a complementary tool of a basin-scale water quality hydrological model as SWAT [18].
In its original version, IberWQ could be used to compute the spatial and temporal evolution of: Escherichia coli, dissolved oxygen (DO), carbonaceous biochemical oxygen demand (CBOD), organic nitrogen (org-N), ammoniacal nitrogen (NH3-N), nitrate-nitrite nitrogen (NO3-N), water temperature (T) and salinity (S). The model solved one 2D depth-averaged advection-diffusion-reaction equation for each water quality variable. The diffusion and advection by the mean flow were computed using the velocity and water depth fields obtained from the hydrodynamic module of Iber, which solves the 2D shallow water equations using an explicit finite volume solver [19].
An important limitation of IberWQ was the fact that it did not simulate variables as phosphorus, phytoplankton and pH. However, these variables are relevant parameters that are often used to evaluate the environmental status of water bodies. Phosphorus is a necessary nutrient for algae and phytoplankton, which can be a limiting factor for plant growth. Neglecting phosphorus and phytoplankton prevented the application of the model to evaluate eutrophication due to the excessive concentration of nutrients (mainly nitrates and phosphates). pH is also an important water quality indicator that should be inside a certain range to preserve the ecosystems. This parameter determines the solubility of chemicals (e.g., metals are more soluble, and thus toxic, at lower pH values) and the amount of nutrients that can be utilized by aquatic life [20].
Another practical limitation of IberWQ was the computational burden required to compute real cases. Most studies imply the simulation of several water quality parameters over large spatial domains and rather long time intervals. The implementation of water quality early warning systems also requires very fast computational codes to immediately and effectively respond to hazardous events. As IberWQ solves the conservation and water quality transport equations on a fine numerical grid, it requires a high computational time (hours or days depending on the specific case). A natural way to overcome this limitation is the implementation of HPC (high performance computing) techniques. This approach had already been implemented for the hydrodynamic module of Iber [21] using GPU (graphical processing unit) parallelization techniques. Due to these developments, the new implementation is able to achieve speed-ups up to nearly two orders of magnitude.
This work presents an improved version of IberWQ model that addresses the previously mentioned limitations. The paper is organized as follows. First, the implementation of new water quality parameters and the parallelization of the code are briefly described in Section 2. New variables include organic phosphorus (org-P), phosphates (PO4-P), phytoplankton, inorganic carbon (inorg-C), alkalinity and pH. Several examples showing the capabilities of the code and the speed-ups achieved in problems with different computational burden are presented in Section 3. Finally, the conclusions are reported in Section 4.

2. Model Structure and Equations

2.1. Model Structure

IberWQ consists of a set of routines that solve a series of unsteady 2D depth-averaged transport equations for different water quality parameters. This includes advection by the mean flow, diffusion due to concentration gradients and reaction with other variables. Those routines are fully integrated within the shallow water model Iber [19], which solves the 2D Saint-Venant equations including turbulence. Thus, the resulting model for water quality and hydrodynamics is fully coupled.
Figure 1 shows the different species and biochemical reactions considered in the extension of the IberWQ module presented in this paper. The new components considered in relation to [14] are organic and inorganic phosphorus, phytoplankton, alkalinity and inorganic carbon. In addition, the pH is computed as a function of inorganic carbon and alkalinity. Phytoplankton is represented as the concentration of Chlorophyll-A (Chl-A), as commonly done in other water quality models [22]. Thus, in the case of computing phytoplankton the user must introduce the ratios between nitrogen, phosphorus, carbon and Chl-A in order to evaluate the corresponding biochemical reactions.
Inorganic carbon includes carbonates ( CO 3 2 ), bicarbonates ( HCO 3 ) and carbon dioxide ( CO 2 ). The three species were modelled together as total inorganic carbon dissolved in water. A posteriori estimation of the concentration of each species was done as a function of pH (specific equations are presented in Appendix B).
The whole model included therefore 14 water quality parameters that are commonly used to evaluate the environmental status of rivers and estuaries (the 13 variables shown in Figure 1 plus pH). The different model components can be activated or deactivated depending on the problem to be solved. The kinetic constants of most of the biochemical reactions depend on water temperature and salinity. In 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.

2.2. Model Equations

A generic 2D depth-averaged advection-diffusion equation with reaction terms is solved for each of the species considered in the model. In order to compute the advection terms, the water quality module is linked to a hydrodynamic module that solves the 2D shallow water equations, given by:
h t + q x x + q y y = 0
q x t + x ( q x 2 h ) + y ( q x q y h ) = g h ( h + z b ) x g n 2 h 7 3 | q | q x 1 2 g h 2 ρ ρ x
q y t + x ( q x q y h ) + y ( q y 2 h ) = g h ( h + z b ) y g n 2 h 7 3 | q | q y 1 2 g h 2 ρ ρ y
t ( h   C i ) + x ( q x C i ) + y ( q y C i ) = x ( h   ν e C i x ) + y ( h   ν e C i y ) + S i h
where h is the water depth, z b is the bed elevation (topography), ( q x ,   q y ) are the two components of the unit discharge, | q | is the modulus of the unit discharge, C i is the depth-averaged concentration of the species i (the available species are shown in Figure 1), S i is a generic reaction term for the species i , ρ is the density of water, g is the gravity acceleration, n is the Manning coefficient and ν e is the effective viscosity computed using a depth-averaged turbulence model. The available turbulence models include the k ε model, a mixing-length model and a parabolic model, and are described in detail in [23]. If the baroclinic pressure is activated in the simulation, the water density is computed from the salinity and temperature fields. For this, the International One Atmosphere Equation of State of Seawater [24] is employed. The total number of partial differential equations solved ( n e q ) depends on the number of species considered on each specific test case ( n i ), and it is equal to n e q = 3 + n i .
The biochemical reactions considered for each species were modeled with the source terms S i . Each solid arrow in Figure 1 represents a biochemical reaction and was modeled as a first order reaction [22,25]. The formulations used to compute the different source terms are basically the same as those used in well-known and extendedly used models such as WASP, QUAL-2K and CE-QUAL-W2, which follow [22]. The details of the equations implemented in IberWQ are given in Appendix A and Appendix B, while the model constants to be introduced by the user are detailed in Appendix C.

2.3. Numerical Solver

The 2D depth-averaged transport equations for each species were solved with an unstructured finite volume solver, using a computational grid formed by triangular and quadrilateral elements. The same finite volume mesh was used to solve the transport equations of the water quality species and the hydrodynamic equations. Thus, the water depth ( h ) and the unit discharges ( q x , q y ) needed to compute the advective fluxes in the scalar transport equations were obtained from the hydrodynamic module of Iber. This module solved the 2D Saint Venant equations using an upwind Godunov scheme with the approximate Riemann solver of Roe. In order to transfer the water velocity and depth from the hydraulic module to the water quality advection–diffusion equations, the mass conservative scheme detailed in [26] was used.
The advective terms can be discretized with either a first order upwind scheme [27] or a second order upwind scheme. In particular, the Gamma scheme proposed in [28] was implemented in the solver to obtain second order accuracy in space. Using this approach only implies an increase on the CPU (Central Processing Unit) time of approximately 5–10%.

2.4. Parallelization

In order to address the limitations in terms of efficiency of the Iber model, a new implementation was developed. The new implementation, Iber+ onwards, was first presented in [21]. It was developed in C++ using the object-oriented paradigm to improve the modularity and the maintenance of the source code. Iber+ was parallelized for shared memory systems using OpenMP and for GPUs using the Nvidia CUDA (Compute Unified Device Architecture) [29] platform. The GPUs are a cost-effective solution to accelerate numerical models without the necessity of expensive HPC facilities. GPUs are a solution available for servers, workstations and laptops. GPU computing is attractive not only to reduce the computational time of the simulations but also to improve the time spent on design and test cases. There are numerous cases of success in the literature for both meshless [30] and mesh-based [31,32,33,34] codes. Iber+ features both CPU and GPU parallelization, however the focus of the implementation is on GPU computing, so in this study only the GPU implementation is analyzed. Iber+ is included in the official package of Iber since version 2.5 and can be downloaded without any cost at http://www.iberaula.com. The new developments presented in this paper will be included in the next release.
Originally developed exclusively for computer graphics, the GPUs feature a highly parallel architecture that provides several TFLOPS (1012 floating point operations per second) of computing power in a single card [35]. This high-throughput hardware can be employed for general applications due to the GPGPU (general processing graphics processing unit) APIs (application programming interface), being Nvidia CUDA one of the most popular in scientific applications. CUDA exposes the GPUs capabilities to programmers due to an extension to common programming languages like C, C++ or Fortran. These extensions provide three abstractions to programmers: the thread hierarchy, shared memory and barrier synchronization [29].
The Nvidia GPUs feature a SIMT (single instruction multiple thread) architecture that issues threads in groups of 32 named warps. Each warp is executed in parallel in a SM (stream multiprocessor). All the threads in a warp execute the same instruction. If a branch instruction is processed, a divergence starts in this case, the threads that follow one of the paths are stalled while the others are executed and vice versa. Once the branch is finished, all the threads converge again. It is essential to consider this peculiarity of the GPUs in order to achieve a high throughput. The algorithms should be revised to avoid unnecessary branching and reorganize data to avoid divergence.
Discrete GPUs have their own memory, independent of the system memory. Data transfers should be made from the system memory to the GPU’s memory, usually through the PCI-Express bus. These data transfers can be a potential bottleneck due to the limited bandwidth and higher latency compared with the system memory accesses. Figure 2 shows the flowchart of Iber+, the memory transfers have been reduced as much as possible. The problem data is transferred before the simulation starts and when it is necessary to write it to the hard disk. Only the current timestep of the simulation is transferred on every loop iteration of the simulation. The write of the simulation state to the hard disk can suppose an important part of the total run time, depending on the case and system configuration. To reduce these times, data is written to disk in background by a thread separated from the main execution thread. On the other hand, the data locality is essential in modern architectures to reduce the cache miss rate and thus improve performance. For this, the elements of the mesh are reordered using a space-filling curve, more specifically the Hilbert curve [36] was used for this purpose.
Global synchronization in GPUs is expensive, making some algorithms like reductions more complicated than their CPU counterpart. In order to implement reductions efficiently, Iber+ employs the library Nvidia CUB (CUDA unbound) [37].
Modern GPUs have significant performance differences when working with single or double precision. The performance ratio single/double precision may vary from 1/2 for the Nvidia Tesla V100 [35] targeted to the HPC market up to 1/32 for the Nvidia RTX 2080ti [38] targeted for the consumer market. For this reason, most of the computations are performed in single precision. The use of double precision is restricted to potentially sensitive variables like the simulation time or the computation of distances.
For the implementation of the water quality module, an object-oriented approach was used. In this way, new species can be easily added by implementing an abstract class without significantly affecting the existing code.

3. Test Cases

In this section, four test cases were analyzed to illustrate the capabilities and performance of IberWQ and its GPU implementation. For a detailed validation of the test cases the reader is addressed to previous publications from the authors [14]. More information about the data sources used for the cases analyzed can be found in Appendix F. All the cases were run in a workstation equipped with an AMD Ryzen 7 2700X processor, 32 GB of RAM (random access memory) and an Nvidia RTX 2080ti graphics card. The run time measurements were carried out using the execution time corresponding to the simulation loop.

3.1. Faecal Contamination in a Coastal Estuary

3.1.1. Description

The first test case simulated the concentration of fecal contamination in a coastal estuary. This is a typical application of 2D depth-averaged water quality models [9,13,39,40,41,42,43] in which only one species, E. coli, needs to be computed.
The model was applied to compute the concentration of E. coli in the coastal estuary of Ferrol, located in the NW of Spain (Figure 3). The estuary extends over 32 km2 with a relatively narrow and elongated shape. The total length of the estuary is 17 km, while its width varies from 400 m to a couple of kilometers in its widest part. Tidal ranges vary from 0.9 (neap tides) to 4.6 m (spring tides), being freshwater inflows from rivers negligible. As we are using a 2D depth-averaged model, we assumed well-mixed conditions over the water depth, which are characteristic of the autumn, winter and spring seasons in this estuary [39].
Approximately six spring semi-diurnal tidal cycles were simulated (3 days), with an approximate tidal range of 4.2 m. Two continuous sewage spills at different locations were considered. The first one (V2 in Figure 3) is located in the narrowest region of the estuary, where the high water velocities activate a rapid mixing with the receiving waters. The second one (V1 in Figure 3) is located in the inner part of the estuary, where mixing and dispersion were slower due to the lower water velocity. Both discharges were characterized by an E. coli concentration of 107 cfu/100 mL and a sewage discharge of 0.1 m3/s, both values were constant in time. Even though the model included the possibility of using the formulation of Mancini to evaluate the degradation of E. coli, in this case we had fixed the value of the degradation constant to 5.5 days−1, following previous studies performed in this estuary [44].
To perform the computations, the estuary was discretized with an unstructured computational grid with 143,993 elements (with an average element size of 222 m2). The levels of E. coli were sampled every 600 s with Iber and Iber+ at control points P1, P2 and P3.

3.1.2. Results

The E. coli concentration was sampled at three control points (P1, P2 and P3 in Figure 3). Figure 4 compares the spatial distribution of E. coli in the estuary at different time steps for the CPU and GPU implementations. The results obtained with Iber and Iber+ were virtually identical. The time series of concentration at the three control points are shown in Figure 5 where no significant differences could be observed. Furthermore, the time series of concentration obtained with the GPU and CPU implementations were statistically compared in terms of the normalized root mean square deviation (in %):
NRMSD = 100 × RMSD y max y min RMSD = i = 1 N ( x i y i ) 2 N
where xi and yi are the values computed with Iber+ and Iber respectively for any variable of interest.
The average error, 〈 NRMSD 〉 = 0.07%, was calculated taking into account all variables at all control points, supporting the good resemblance between series suggested by the visual comparison.
Table 1 shows the performance measurements for the first test case. This was the best case for GPU implementation because it had the largest mesh (143,993 elements) among the cases under study. GPU computing implied a certain overhead, mainly due to memory transfers and kernel launch latency. The higher the GPU workload the less significant the overhead will be. In this case Iber took more than 16 h to complete the simulation whilst Iber+ could run the same problem in less than 6 min. This implied a speedup of 181.

3.2. Organic Matter Contamination in an Estuary

3.2.1. Description

In the second test case the concentration of dissolved oxygen (DO) and carbonaceous biochemical oxygen demand (CBOD) were computed in the estuary of A Coruña, located in the NW of Spain (Figure 6). The estuary has an area of 26 km2. The outer part of the estuary is relatively deep (with depths of the order of 20 m and maximum values of circa 30 m in the mouth), while the inner part, where the sewage spills are located, is relatively shallow (with maximum water depths of the order of 5 m at high tide and many dry regions at low tide).
The physical time simulated for this case extends over 8 tidal cycles (4 days) in which the tidal range moves from neap tides (with a tidal range of 1.2 m) to spring tides (with a tidal range of 3.5 m). The river discharge at the upstream boundary of the model is equal to 20 m3/s, which corresponds to its annual average value. Four discontinuous sewage spills of CBOD were considered, all of them located in the inner estuary (Figure 6). The mean depth at the spills varied within 3 m and 4 m, which were relatively shallow depths. All the spills were characterized by the same uniform concentration of CBOD and DO (511 mg/L of CBOD and 3.5 mg/L of DO). The water discharges were discontinuous and different from one spill to another, with maximum values ranging from 0.06 to 0.5 m3/s. The CBOD degradation rate at 20 °C was set to k d b o c = 0.23   day 1 .
The domain was discretized using a computational grid of 50,329 elements. The mesh size was variable over the estuary. A finer resolution was used in the inner part of the estuary (element sizes of 10 m) since in this region the spills were located and unsteady wet-dry tidal fronts appeared. The size of the elements in the mouth of the estuary was approximately 180 m. The levels of DO and CBOD were sampled every 300 s using both Iber and Iber+ at control points P1, P2 and P3.

3.2.2. Results

In the second test case, the levels of CBOD and DO were sampled at four control points shown in Figure 6. In particular, the time series of CBOD computed at these control points are shown in Figure 7 with identical results for both implementations. Regarding DO concentration (Figure 8), very small differences could be appreciated at control point P2, where the peak values computed with Iber+ were slightly higher than those computed with Iber. Statistically, the differences were not significant with a negligible average error (〈 NRMSD 〉 = 0.32%, covering all points and variables), although higher than in previous case.
This case employed a mesh with 50,329 elements, about a third of the mesh elements used in the first test case. As shown in Table 2, Iber+ was able to perform the simulation 61 times faster than Iber. In spite of having a lower number of elements, Table 2 shows that the real time employed by Iber+ to compute a simulation time step was longer than in the previous test. This can be explained by two reasons. First, this case had to compute two water quality species instead of one, causing a larger kernel launch overhead. Second, the size of the problem was not big enough to saturate the GPU capacity. Thus, this case was not able to take full advantage of the GPU computing capacity, as it happened in the previous case.

3.3. Combined Sewer Overflows in a River Miño Reach

3.3.1. Description

This test case was extracted from the study presented in [1]. In that work, a 2D water quality model is used as a fundamental part of an integrated modeling approach for the design of the sewer network of the city of Lugo (Spain). Here, the impact of combined sewer overflows (CSO) in the river Miño was evaluated by means of the Environmental Quality Standards (EQS) presented in the Urban Pollution Manual [44]. The application of EQS requires an efficient water quality model since usually a large number of simulations must be run over long periods of time.
The objective variables computed in [44] to apply the EQS are the concentrations of dissolved oxygen and ammonia. The evaluation of these variables requires the computation of the five following species: org-N, NH3-N, NO3-N, DO and CBOD.
In the example presented here, two days of discontinuous sewage discharges at three different locations along the river were modeled. All the spills were characterized by the following concentrations: 102 mg/L of CBOD, 5 mg/L of org-N, 1.5 mg/L of NH3-N and 4 mg/L of DO. The concentration of nitrates (NO3-N) in the spills is negligible. The spill discharges were discontinuous in time and different from each other, with maximum values of discharge ranging from 0.18 to 1.5 m3/s. The ambient river concentration of the water quality species was negligible with the exception of the DO concentration, which was set to 9 mg/L. Those ambient values were imposed as initial and upstream boundary conditions in the model. The river discharge, which was obtained from a Water Quality Automatic Information System (SAICA) located upstream the river reach under study, varied from 40 to 60 m3/s during the two days of computation. The same values of the reaction kinetic constants proposed in [45] were used in the simulations, namely (all values at 20 °C): CBOD degradation rate equal to k d b o c = 0.35   day 1 , org-N hydrolysis rate equal to k h n = 0.20   day 1 , nitrification rate equal to k n i t = 0.50   day 1 and denitrification rate equal to k d e n i t = 0.05   day 1 .
Figure 9 shows the bathymetry of the river reach under study and the location of the sewage spills. The computational domain extends over 8 km of river, with an average width of 70 m and a total extension of 0.57 km2. The domain was meshed with 8763 elements with an average size of 65 m2. The levels of the simulated species were sampled every 600 s with both models at the control points P1, P2, P3 and P4.

3.3.2. Results

In the third case, five species (org-N, NH3-N, NO3-N, DO and CBOD) were sampled at the four control points indicated in Figure 9. For the sake of clarity, only the results obtained at P3 are shown in Figure 10, while the rest of time series are shown in Appendix D. In this simulation, the average error (〈 NRMSD 〉 = 0.04%) was negligible and smaller than in the previous two cases.
In this case, Iber+ completed the simulation 29 times faster than Iber (Table 3). This speedup was considerably lower than the one obtained in the previous two cases due to the relatively small size of the mesh, with just 8763 elements. Compared to the previous case, the time needed to process a single time step was lower whilst the number of cells processed per second was much lower. This indicates the difficulty of taking full advantage of the GPU resources in cases with a small number of elements, such as the present one. Nevertheless, even with a small mesh, Iber+ ran significantly faster (29 times) than the non-parallelized CPU version.

3.4. Effluent Discharge from a Wastewater Treatment Plant

3.4.1. Description

In this case an effluent discharge from a wastewater treatment plant into a river located in the south of Spain was modeled. The water quality species computed in this case are DO, CBOD, org-N, NH3-N, NO3-N, org-P, PO4-P, inorg-C, alkalinity, salinity and pH. Thus, 10 transport equations were solved in addition to the three Saint Venant equations (as mentioned in Section 2, the pH did not have an associated transport equation, instead its value was computed at each mesh element from the values of alkalinity and inorg-C).
The river reach modeled (Figure 11) is 18 km long and approximately 200 m wide, with a total area of 4 km2. The downstream boundary of the reach modeled is located 50 km upstream the river mouth into the ocean and is affected by the ocean tidal elevation, although tide does not propagate all the way up to the upstream boundary. Thus, a tidal wave with a range of 1.0 m was imposed at the downstream boundary while a constant river discharge of 3.6 m3/s was imposed at the upstream boundary. The ambient concentrations of the water quality components imposed at the upstream boundary were the following: 8 mg/L of DO, 0 mg/L of CBOD, 2.8 mg/L of org-N, 1.3 mg/L of NH3-N, 4.2 mg/L of NO3-N, 0.3 mg/L of org-P, 0.2 mg/L of PO4-P, 90 mg/L of inorg-C, 346 mg/L CO3Ca (alkalinity) and 0.8 mg/L of salt. These values correspond to mean winter conditions in this river reach.
A continuous spill from a waste water treatment plant (Figure 11) was modeled, with a discharge of 2.1 m3/s and the following concentration of the water quality variables considered in this case: 5 mg/L of DO, 23.6 mg/L of CBOD, 23.9 mg/L of org-N, 28.4 mg/L of NH3-N, 0.3 mg/L of NO3-N, 0.8 mg/L of org-P, 0.4 mg/L of PO4-P, 97 mg/L of inorg-C, 392 mg/L CO3Ca (alkalinity) and 1 mg/L of salt.
The following kinetic constants were used in the simulations (all values at 20 °C): CBOD degradation rate ( k d b o c = 0.23   day 1 ), org-N hydrolysis rate ( k h n = 0.20   day 1 ), nitrification rate ( k n i t = 0.20   day 1 ), denitrification rate ( k d e n i t = 0.05   day 1 ) and organic phosphorus hydrolysis rate ( k h p = 0.20   day 1 ).
The simulated physical time was 3 days, using an unstructured mesh of 90,406 elements. At the end of the simulation a steady state was achieved for the whole reach. All the simulated species were sampled every 600 s with both Iber and Iber+ at control points P1, P2, P3 and P4.

3.4.2. Results

In the last case, ten different species were simulated (DO, CBOD, org-N, NH3-N, NO3-N, org-P, PO4-P, inorg-C, alkalinity, salinity and pH). Time series of concentrations were taken at Points P1 to P4 as indicated in Figure 11. The zone of the sewage spill is shown in detail in Figure 12. The time series of concentration at control point P3 are shown in Figure 13. For the sake of clarity, the results at the other control points were included in Appendix E. The results provided by Iber and Iber+ were almost identical in all cases, confirmed by the negligible average error (〈 NRMSD 〉 = 0.16%) as obtained for the rest of the cases.
Table 4 shows the performance measurements for this case, where a mesh of 90,406 elements (halfway between Test 1 and Test 2) was considered. Iber+ reached a speedup of 92 with respect to Iber. Comparing with Test 2, the new mesh was almost double the size of the former one and ten species were simulated instead of two. As a result of this, the time required to process a single time step was almost 3.5 times longer. However, the number of mesh cells processed per second was only twice lower. In this case, Iber+ could take more advantage of the GPU processing capacity due to the bigger mesh size.

4. Conclusions

Numerical models are useful tools for the evaluation of the environmental status of water bodies. However, it should be noted that these kinds of tools rely on a set of simplifications and are much simpler than the real systems, whose complexity cannot be fully reproduced by the models. Hence, they can provide wrong results, especially when the managers or the decision makers are not aware of those limitations. The parametrization of the models is then a crucial task in order to obtain precise results, which requires an exhaustive sensitivity analysis. In the particular case of Iber+, its 2D nature makes it only applicable to rivers and non-stratified estuaries.
On the other hand, the main advantage of 2D models like Iber+ is their affordable execution time, which makes them especially valuable for fast response purposes when compared with 3D models. However, its application is hindered by the high spatial resolution required for water quality studies in long river reaches or large estuaries.
This paper presented an improved version of a two-dimensional depth-averaged water quality model, whose efficiency was improved by implementing HPC techniques based on GPU parallelization. The model considered the water quality parameters most commonly used in the environmental assessment of receiving waters, including dissolved oxygen, CBOD, organic nitrogen, ammonia, nitrates, organic phosphorus, phosphates, pH, salinity and temperature. The implementation of the code ran in NVIDIA GPUs that are commonly installed in standard laptop and desktop PCs. In the test cases presented here, speedups in computational time between 29 and 181 were obtained when compared with the non-parallelized implementation, keeping the accuracy of the original model. The code will be integrated in the software package Iber, making it freely available.
In summary, the model presented three key features that make it very attractive and useful for the scientific and engineering community: simulation of the most common water quality parameters, high-performance computing in standard PCs and free availability.

Author Contributions

O.G.-F., L.C. and M.G.-G. conceived the study; O.G.-F. and L.C. developed the software; L.C. and J.M.D. supervised the code development; O.G.-F. and J.G.-C. performed the experiments; O.G.-F., L.C., J.G.-C. and M.G.-G. analyzed the results; O.G.-F., L.C., J.G.-C., J.M.D. and M.G.-G. wrote the manuscript. All authors have read and agreed to the published version of the manuscript.

Funding

This research was partially supported by INTERREG-POCTEP under project RISC_ML (Code: 0034_RISC_ML_6_E) co-funded by European Regional Development Fund (ERDF); and by Xunta de Galicia under Project ED431C 2017/64-GRC “Programa de Consolidación e Estruturación de Unidades de Investigación Competitivas (Grupos de Referencia Competitiva)”. O.G.F. is supported by Xunta de Galicia grant ED481A-2017/314.

Acknowledgments

The aerial pictures used in this work are courtesy of the Spanish IGN (Instituto Geográfico Nacional) and part of the PNOA (Plan Nacional de Ortofotografía Aérea) program.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A. Reaction Terms

The following source terms were included in the advection–diffusion equation to model the biochemical reactions between the water quality species considered in the model, as represented in Figure 1. Salinity and alkalinity were considered as conservative variables with no source terms.
Chlorophyll-A
S C h l A = S C h l A p h o t o S C h l A r e s p S C h l A d e a t h S C h l A d e p
S C h l A p h o t o = μ m a x F L m i n ( F N , F P , F C )   1.047 ( T 20 )   C C h l A
S C h l A r e s p = k r p   F o x p   1.047 ( T 20 )   C C h l A
S C h l A d e a t h = k d p   1.047 ( T 20 )   C C h l A
S C h l A d e p = V s A h C C h l A
F N = C N H 3 - N + C N O 3 - N C N H 3 - N + C N O 3 - N + k s n p F P = C P O 4 - P C P O 4 - P + k s p p F C = C i n o r g C C i n o r g C + k s c p
F L = 1 k e   h L n ( k l p + I 0 k l p + I · e k e h )
F o x p = C D O C D O + k s o p
Nitrogen
S o r g N = S o r g N h y d r o S o r g N d e p + r n a   S C h l A d e a t h
S N H 3 - N = S o r g N h y d r o S N H 3 - N n i t + r n a   S C h l A r e s p r n a   P a p   S C h l A p h o t o
S N O 3 - N = S N H 3 - N n i t S N O 3 - N d e n i t r n a   ( 1 P a p )   S C h l A p h o t o
S o r g N h y d r o = k h n   1.047 ( T 20 ) C o r g N  
S o r g N d e p = V s N h C o r g N
S N H 3 - N n i t = k n i t   1.083 ( T 20 ) F n   C N H 3 - N
S N O 3 - N d e n i t = k d e n i t   1.045 ( T 20 )   F d n   C N O 3 - N
P a p = k p a C N H 3 - N k p a C N H 3 - N + ( 1 k p a ) C N O 3 - N
F n = C D O k n 1 / 2 + O D ,   F d n = k d n 1 / 2 k d n 1 / 2 + C D O
Phosphorus
S o r g P = S o r g P h y d r o S o r g P d e p + r p a   S C h l A d e a t h
S o r g P h y d r o = k h p   1.047 ( T 20 )   C o r g P
S o r g P d e p = V s P o h C o r g P
S P O 4 - P = S o r g P h y d r o S P O 4 - P d e p + r p a   S C h l A r e s p r p a   P a p   S C h l A p h o t o
S P O 4 - P d e p = V s P i h C P O 4 - P
CBOD
S C B O D =   S C B O D d e g S C B O D d e p 2.86   S N O 3 - N d e n i t
S C B O D d e g = k C B O D   1.047 ( T 20 )   F o x c   C C B O D
S C B O D d e p = V s C B O D h C C B O D
F o x c = C D O k s o c f + C D O
Dissolved Oxygen
S D O = S D O r e a i r S C B O D d e g r a   S N H 3 - N n i t + r o a   S C h l A p h o t o r o a   S C h l A r e s p S D O D O S
S D O r e a i r = k a i r   1.024 ( T 20 )   ( DO sat C DO )
S D O S O D = k S O D h
k a i r = k a i r h + k a i r w h
k a i r h = { 5.32 U 0.67 h 1.85   i f   h < 0.61   m 3.93 U 0.5 h 1.5   i f   h > 3.45   U 2.5   & h > 0.61   m 5.026 U h 1.67 o t h e r w i s e
k a i r w = 0.728   V w 10 0.5 0.317   V w 10 + 0.0372   V w 10 2
Inorganic Carbon
S i n o r g C = S i n o r g C r e a i r + r c c o S C B O D d e g + r c c a   S C h l A r e s p r c c a   S C h l A p h o t o
S i n o r g C r e a i r = 0.923   k a i r   1.024 ( T 20 )   ( C O 2 s a t α 0   C i n o r g C )
r c c o = 1 12 1 r o c r c c a = 1 12 r c a
E. coli
S E C = k d e c × E C
k d e c = ( 0.8 + 0.02 S ) 1.07 ( T 20 ) + 0.086 I 0 k e h ( 1 e ( k e h ) )
Temperature
S T = q r a d q b r q c o n d q e v a p
q b r = 0.97   σ   T 4 σ = 5.669 × 10 8 W m 2 K 4
q c o n d = 0.47   ( 19 + 0.95   V w 7 2 )   ( T T a i r )
q e v a p = ( 19 + 0.95   V w 7 2 )   ( e w a t e r e a i r )
e w a t e r = 4.596 × e x p ( 17.27   ( T 273.15 ) T 35.7 ) e a i r = R H × 4.596 × e x p ( 17.27   ( T a i r 273.15 ) T a i r 35.7 )

Appendix B. pH model

pH computations are based on the values of alkalinity and inorganic carbon. The pH model proposed by [45] was used. The model was based on the following equilibrium, mass balance and electroneutrality equations.
K 1 = [ H C O 3 ]   [ H + ] [ H 2 C O 3 * ] K 2 = [ C O 3 2 ]   [ H + ] [ H C O 3 ] K w = [ H + ] [ O H ] C i n o r g C = [ H 2 C O 3 * ] + [ H C O 3 ] + [ C O 3 2 ] A l k = [ H C O 3 ] + 2   [ C O 3 2 ] + [ O H ] [ H + ]
where K 1 , K 2 and K w are equilibrium constants, A l k is the alkalinity (eq/L), C i n o r g C is the concentration of inorganic carbon (mol/L), [ H 2 C O 3 * ] is the sum of carbon dioxide and carbonic acid dissolved in the water (mol/L), [ H C O 3 ] is the concentration of bicarbonate ions, [ C O 3 2 ] is the concentration of carbonate ions, [ H + ] is the concentration of hydrogen ions and [ O H ] is the concentration of hydroxyl ions. The equilibrium constants are computed as a function of water temperature as:
log 10 ( K w ) = 4787.3 T 7.1321 log 10 ( T ) 0.010365 × T + 22.8 log 10 ( K 1 ) = 356.3094 0.0609196 × T + 21834.37 T + 126.8339 log 10 ( T ) 1684915 T 2 log 10 ( K 2 ) = 107.887 0.0325285 × T + 5151.79 T + 38.9256 log 10 ( T ) 563713.9 T 2
The values of C i n o r g C and A l k are obtained from the solution of their respective depth-averaged transport equations. The previous system of five equations is solved as proposed in [44], by solving the following algebraic non-linear equation for [ H + ] at each mesh element:
( α 1 + 2 α 2 )   C i n o r c C + K w [ H + ] [ H + ] A l k = 0
with:
α 0 = [ H + ] 2 [ H + ] 2 + K 1   [ H + ] + K 1 K 2 α 1 = K 1   [ H + ]   [ H + ] 2 + K 1   [ H + ] + K 1   K 2 α 2 = K 1   K 2   [ H + ] 2 + K 1   [ H + ] + K 1   K 2
The coefficients α 0 , α 1 and α 2 represent respectively the fraction of inorganic carbon in form of carbon dioxide, bicarbonates and carbonates. Once the previous equation is solved at each mesh element, the pH can be computed as:
p H = log 10 [ H + ]
and the concentration of carbon dioxide, bicarbonates and carbonates as:
C H 2 C O 3 * = α 0 C i n o r g C C H C O 3 = α 1 C i n o r g C C C O 3 2 = α 2 C i n o r g C

Appendix C. Model Constants

ConstantUnitsSuggested ValuesDescription
MinMax
r n a mg/mg0.79.0Ratio of Nitrogen to Chl-A in phytoplankton
r p a mg/mg0.12.0Ratio of Phosphorus to Chl-A in phytoplankton
r o a mg/mg14180Ratio of Oxygen to Chl-A in phytoplankton
r c a mg/mg--Ratio of Carbon to Chl-A in phytoplankton
k r p 1/day0.040.8Phytoplankton respiration rate
μ m a x 1/day1.03.0Maximum photosynthesis rate
k d p 1/day0.050.5Phytoplankton death rate
k p a -0.01.0Phytoplankton preference factor for ammonia
V s A m/day--Phytoplankton settling velocity
k s n p mg/L0.010.3Nitrogen half-saturation constant for photosynthesis attenuation
k s p p mg/L0.0010.05Phosphorus half-saturation constant for photosynthesis attenuation
k s c p mg/L--Carbon half-saturation constant for photosynthesis attenuation
k l p W/m20.050.3Light half-saturation constant for photosynthesis attenuation
k s o p mg/L--Oxygen half-saturation constant for respiration attenuation
k e 1/m--Light extinction coefficient in water
V s P o m/day--Organic phosphorus settling velocity
V s P i m/day--Inorganic phosphorus settling velocity
k h p 1/day0.010.7Organic phosphorus hydrolysis rate at 20 °C
k n i t 1/day0.011.0Nitrification rate at 20 °C
k h n 1/day0.020.4Organic nitrogen hydrolysis rate at 20 °C
k d e n i t 1/day0.0010.1Denitrification rate at 20 °C
V s N m/day0.0010.1Organic nitrogen settling velocity
k n 1 / 2 mg/L Oxygen half-saturation constant for nitrification attenuation
k d n 1 / 2 mg/L Oxygen half-saturation constant for denitrification attenuation
k d b o c 1/day0.023.4CBOD degradation rate at 20 °C
V s D B O C m/day0.010.36CBOD settling velocity
k s o c f mg/L Oxygen half-saturation constant for CBOD degradation attenuation
k s o d kg/m2/day0.00.01Sediment oxygen demand rate
r o c mg/mg--Ratio of oxygen consumed per organic carbon oxidized to inorganic carbon
k d e c 1/dayMancini Degradation constant for E. coli

Appendix D. Time Series for Test 3

Figure A1. Time series of dissolved oxygen for Test 3, sampled at control points P1 (a), P2 (b), P3 (c) and 4 (d).
Figure A1. Time series of dissolved oxygen for Test 3, sampled at control points P1 (a), P2 (b), P3 (c) and 4 (d).
Water 12 00413 g0a1
Figure A2. Time series of CBOD for Test 3, sampled at control points P1 (a), P2 (b), P3 (c) and P4 (d).
Figure A2. Time series of CBOD for Test 3, sampled at control points P1 (a), P2 (b), P3 (c) and P4 (d).
Water 12 00413 g0a2
Figure A3. Time series of organic nitrogen for Test 3, sampled at control points P1 (a), P2 (b), P3 (c) and 4 (d).
Figure A3. Time series of organic nitrogen for Test 3, sampled at control points P1 (a), P2 (b), P3 (c) and 4 (d).
Water 12 00413 g0a3
Figure A4. Time series of ammoniacal nitrogen for Test 3, sampled at control points P1 (a), P2 (b), P3 (c) and P4 (d).
Figure A4. Time series of ammoniacal nitrogen for Test 3, sampled at control points P1 (a), P2 (b), P3 (c) and P4 (d).
Water 12 00413 g0a4
Figure A5. Time series of nitrate-nitrite nitrogen for Test 3, sampled at control points P1 (a), P2 (b), P3 (c) and P4 (d).
Figure A5. Time series of nitrate-nitrite nitrogen for Test 3, sampled at control points P1 (a), P2 (b), P3 (c) and P4 (d).
Water 12 00413 g0a5

Appendix E. Time Series for Test 4

Figure A6. Time series of pH for Test 4, sampled at control points P1 (a), P2 (b), P3 (c) and P4 (d).
Figure A6. Time series of pH for Test 4, sampled at control points P1 (a), P2 (b), P3 (c) and P4 (d).
Water 12 00413 g0a6
Figure A7. Time series of alkalinity for Test 4, sampled at control points P1 (a), P2 (b), P3 (c) and P4 (d).
Figure A7. Time series of alkalinity for Test 4, sampled at control points P1 (a), P2 (b), P3 (c) and P4 (d).
Water 12 00413 g0a7
Figure A8. Time series of inorganic carbon for Test 4, sampled at control points P1 (a), P2 (b), P3 (c) and P4 (d).
Figure A8. Time series of inorganic carbon for Test 4, sampled at control points P1 (a), P2 (b), P3 (c) and P4 (d).
Water 12 00413 g0a8
Figure A9. Time series of inorganic phosphorus for Test 4, sampled at control points P1 (a), P2 (b), P3 (c) and P4 (d).
Figure A9. Time series of inorganic phosphorus for Test 4, sampled at control points P1 (a), P2 (b), P3 (c) and P4 (d).
Water 12 00413 g0a9
Figure A10. Time series of organic phosphorus for Test 4, sampled at control points P1 (a), P2 (b), P3 (c) and P4 (d).
Figure A10. Time series of organic phosphorus for Test 4, sampled at control points P1 (a), P2 (b), P3 (c) and P4 (d).
Water 12 00413 g0a10
Figure A11. Time series of nitrite-nitrate nitrogen for Test 4, sampled at control points P1 (a), P2 (b), P3 (c) and P4 (d).
Figure A11. Time series of nitrite-nitrate nitrogen for Test 4, sampled at control points P1 (a), P2 (b), P3 (c) and P4 (d).
Water 12 00413 g0a11
Figure A12. Time series of ammoniacal nitrogen for Test 4, sampled at control points P1 (a), P2 (b), P3 (c) and P4 (d).
Figure A12. Time series of ammoniacal nitrogen for Test 4, sampled at control points P1 (a), P2 (b), P3 (c) and P4 (d).
Water 12 00413 g0a12
Figure A13. Time series of organic nitrogen for Test 4, sampled at control points P1 (a), P2 (b), P3 (c) and P4 (d).
Figure A13. Time series of organic nitrogen for Test 4, sampled at control points P1 (a), P2 (b), P3 (c) and P4 (d).
Water 12 00413 g0a13
Figure A14. Time series of CBOD for Test 4, sampled at control points P1 (a), P2 (b), P3 (c) and P4 (d).
Figure A14. Time series of CBOD for Test 4, sampled at control points P1 (a), P2 (b), P3 (c) and P4 (d).
Water 12 00413 g0a14
Figure A15. Time series of dissolved oxygen for Test 4, sampled at control points P1 (a), P2 (b), P3 (c) and P4 (d).
Figure A15. Time series of dissolved oxygen for Test 4, sampled at control points P1 (a), P2 (b), P3 (c) and P4 (d).
Water 12 00413 g0a15

Appendix F. Data Sources for the Test Cases

Table A1. Data sources for the Test Case 1.
Table A1. Data sources for the Test Case 1.
BathymetryBathymetric survey carried out for previous studies. Spatial resolution of 30 m.
Effluent Discharge and ConcentrationVirtual
StreamflowAnnual average flow from river Grande de Xubia. Available from the regional Meteorological Agency MeteoGalicia (www.meteogalicia.gal).
TideTidal harmonics obtained from the tidal gauge of Ferrol. Available from Puertos del Estado (www.puertos.es).
Table A2. Data sources for the Test Case 2.
Table A2. Data sources for the Test Case 2.
BathymetryBathymetric survey carried out for previous studies. Spatial resolution of 30 m.
Effluent Discharge and ConcentrationSewer network model carried out in previous studies.
StreamflowAnnual average flow from river Mero. Available from the regional Meteorological Agency MeteoGalicia (www.meteogalicia.gal).
TideTidal harmonics obtained from the tidal gauge of A Coruña. Available from Puertos del Estado (www.puertos.es).
Table A3. Data sources for the Test Case 3.
Table A3. Data sources for the Test Case 3.
BathymetryBathymetric survey carried out in [44].
Effluent Discharge and ConcentrationSewer network model carried out in [46] and [46].
StreamflowRiver discharge obtained from a Water Quality Automatic Information System (SAICA) located upstream the river reach under study. Available from the regional water administration Confederación Hidrográfica del Miño-Sil (www.chminosil.es).
Table A4. Data sources for the Test Case 4.
Table A4. Data sources for the Test Case 4.
BathymetryDigital terrain model at 2 m resolution, obtained from LiDAR data from the Spanish National Plan of Aerophotogrammetry (PNOA), available from the Spanish National Geographic Institute (www.ign.es).
Effluent Discharge and ConcentrationVirtual.
StreamflowAnnual average flow obtained from the regional water administration Confederación Hidrográfica del Guadalquivir (www.chguadalquivir.es).

References

  1. Foundation of Water Research. Urban Pollution Management Manual: A Planning Guide for the Management of Urban Wastewater Discharges During Wet Weather; Foundation of Water Research: Buckinghamshire, UK, 1998. [Google Scholar]
  2. Robinson, R.B.; Roby, J.C. Concentration-duration-frequency curves for pH in a stream in the great smoky mountains. J. Environ. Eng. 2006, 132, 1600–1605. [Google Scholar] [CrossRef]
  3. Schwartz, J.S.; Dahle, M.; Bruce Robinson, R. Concentration-duration-frequency curves for stream turbidity: Possibilities for assessing biological impairment. J. Am. Water Resour. Assoc. 2008, 44, 879–886. [Google Scholar] [CrossRef]
  4. Di Toro, D.; Fitzpatrick, J.; Thomann, R. Documentation for Water Quality Analysis Simulation Program (WASP) and Model Verification Program (MVP); United States Environmental Protection Agency: Washington, DC, USA, 1983.
  5. Ambrose, R. WASP4, a Hydrodynamic and Water Quality Model: Model Theory, User’s Manual and Programmer’s Guide; United States Environmental Protection Agency: Washington, DC, USA, 1988.
  6. Chapra, S.C.; Pelletier, G.J.; Tao, H. QUAL2K: A Modeling Framework for Simulating River and Stream Water Quality; Version 2.11: Documentation and Users Manual; Civil and Environmental Engineering Dept., Tufts University: Medford, MA, USA, 2008. [Google Scholar]
  7. Cole, T.M.; Buchak, E.M. CE-QUAL-W2: A Two-Dimensional, Laterally Averaged, Hydrodynamic and Water Quality Model; Version 2.0: User Manual; US Army Corps of Engineers: Vicksburg, MS, USA, 1995. [Google Scholar]
  8. Cole, T.M.; Wells, S.A. CE-QUAL-W2: A Two-Dimensional, Laterally Averaged, Hydrodynamic and Water Quality Model; Version 3.5; US Army Engineering and Research Development Center: Vicksburg, MS, USA, 2006. [Google Scholar]
  9. Kashefipour, S.M.; Lin, B.; Harris, E.; Falconer, R.A. Hydro-environmental modelling for bathing water compliance of an estuarine basin. Water Res. 2002, 36, 1854–1868. [Google Scholar] [CrossRef]
  10. Kay, D.; Stapleton, C.M.; Wyer, M.D.; McDonald, A.T.; Crowther, J.; Paul, N.; Jones, K.; Francis, C.; Watkins, J.; Wilkinson, J.; et al. Decay of intestinal enterococci concentrations in high-energy estuarine and coastal waters: Towards real-time T 90 values for modelling faecal indicators in recreational waters. Water Res. 2005, 39, 655–667. [Google Scholar] [CrossRef] [PubMed]
  11. Wu, Y.; Falconer, R.; Lin, B. Modelling trace metal concentration distributions in estuarine waters. Estuarine Coast. Shelf Sci. 2005, 64, 699–709. [Google Scholar] [CrossRef]
  12. Gao, G.; Falconer, R.A.; Lin, B. Numerical modelling of sediment–bacteria interaction processes in surface waters. Water Res. 2011, 45, 1951–1960. [Google Scholar] [CrossRef] [Green Version]
  13. Abu-Bakar, A.; Ahmadian, R.; Falconer, R.A. Modelling the transport and decay processes of microbial tracers in a macro-tidal estuary. Water Res. 2017, 123, 802–824. [Google Scholar] [CrossRef]
  14. Cea, L.; Bermúdez, M.; Puertas, J.; Bladé, E.; Corestein, G.; Escolano, E.; Conde, A.; Bockelmann-Evans, B.; Ahmadian, R. IberWQ: New simulation tool for 2D water quality modelling in rivers and shallow estuaries. J. Hydroinf. 2016, 18, 816–830. [Google Scholar] [CrossRef] [Green Version]
  15. Frick, W.E.; Roberts, P.J.W.; Davis, L.R.; Keyes, J.; Baumgartner, D.J.; George, K.P. Dilution Models for Effluent Discharges, 4th ed.; (Visual Plumes); U.S. Environmental Protection Agency, Office of Research and Development: Washington, DC, USA, 2003.
  16. Doneker, R.L.; Jirka, G.H. CORMIX-GI systems for mixing zone analysis of brine wastewater disposal. Desalination 2001, 139, 263–274. [Google Scholar] [CrossRef]
  17. Kiesel, J.; Schmalz, B.; Brown, G.L.; Fohrer, N. Application of a hydrological-hydraulic modelling cascade in lowlands for investigating water and sediment fluxes in catchment, channel and reach. J. Hydrol. Hydromech. 2013, 61, 334–346. [Google Scholar] [CrossRef] [Green Version]
  18. D’Ambrosio, E.; Gentile, F.; de Girolamo, A.M. Assessing the sustainability in water use at the basin scale through water footprint indicators. J. Clean. Prod. 2019, 244, 118847. [Google Scholar] [CrossRef]
  19. Bladé, E.; Cea, L.; Corestein, G.; Escolano, E.; Puertas, J.; Vázquez-Cendón, E.; Dolz, J.; Coll, A. Iber: herramienta de simulación numérica del flujo en ríos. Rev. Int Metodos Numer. Para Calculo Diseno Ing. 2014, 30, 1–10. [Google Scholar] [CrossRef] [Green Version]
  20. Michaud, J.P. A Citizien’s Guide to Understanding and Monitoring Lakes and Streams; Washington State Department of Ecology, Publications Office: Olympia, WA, USA, 1991.
  21. García-Feal, O.; González-Cao, J.; Gómez-Gesteira, M.; Cea, L.; Domínguez, J.M.; Formella, A. An accelerated tool for flood modelling based on Iber. Water (Switzerland) 2018, 10, 1459. [Google Scholar] [CrossRef] [Green Version]
  22. Chapra, S.C. Surface Water-Quality Modeling; McGraw-Hill: New York, NY, USA, 1997; ISBN 1478608307. [Google Scholar]
  23. Cea, L.; Puertas, J.; Vázquez-Cendón, M.E. Depth averaged modelling of turbulent shallow water flow with wet-dry fronts. Arch. Comput. Meth. Eng. 2007, 14, 303–341. [Google Scholar] [CrossRef]
  24. Millero, F.J.; Poisson, A. International one-atmosphere equation of state of seawater. Deep Sea Res. Part A 1981, 28, 625–629. [Google Scholar] [CrossRef]
  25. Bowie, G.; Mills, W.; Porcella, D.; Campbell, C.; Pagenkopf, J.; Rupp, G.; Johnson, K.; Chan, P.; Gherini, S.; Chamberlin, C. Rates, Constants, and Kinetics Formulations in Surface Water Quality Modeling; EPA: Athens, GA, USA, 1985.
  26. Cea, L.; Vázquez-Cendón, M.E. Unstructured finite volume discretisation of bed friction and convective flux in solute transport models linked to the shallow water equations. J. Comput. Phys. 2012, 231, 3317–3339. [Google Scholar] [CrossRef]
  27. Versteeg, H.K.; Malalasekera, W. An Introduction to Computational Fluid Dynamics—The Finite Volume Method; Pearson Education: New York, NY, USA, 2007. [Google Scholar]
  28. Jasak, H.; Weller, H.G.; Gosman, A.D. High resolution NVD differencing scheme for arbitrarily unstructured meshes. Int. J. Numer. Methods Fluids 1999, 31, 431–449. [Google Scholar] [CrossRef]
  29. Nvidia Corporation CUDA C Programming Guide. Available online: https://docs.nvidia.com/cuda/pdf/CUDA_C_Programming_Guide.pdf (accessed on 17 December 2019).
  30. Crespo, A.J.C.; Domínguez, J.M.; Rogers, B.D.; Gómez-Gesteira, M.; Longshaw, S.; Canelas, R.; Vacondio, R.; Barreiro, A.; García-Feal, O. DualSPHysics: Open-source parallel CFD solver based on Smoothed Particle Hydrodynamics (SPH). Comput. Phys. Commun. 2015, 187, 204–216. [Google Scholar] [CrossRef]
  31. Crossley, A.; Lamb, R.; Waller, S. Fast solution of the Shallow Water Equations using GPU technology. In Proceedings of the British Hydrological Society 3rd International Symposium, Newcastle, UK, 13–19 July 2010. [Google Scholar]
  32. Vacondio, R.; Dal Palù, A.; Mignosa, P. GPU-enhanced finite volume shallow water solver for fast flood simulations. Environ. Modell. Softw. 2014, 57, 60–75. [Google Scholar] [CrossRef]
  33. Lacasta, A.; Morales-Hernández, M.; Murillo, J.; García-Navarro, P. An optimized GPU implementation of a 2D free surface simulation model on unstructured meshes. Adv. Eng. Softw. 2014, 78, 1–15. [Google Scholar] [CrossRef]
  34. Liu, Q.; Qin, Y.; Li, G. Fast Simulation of Large-Scale Floods Based on GPU Parallel Computing. Water 2018, 10, 589. [Google Scholar] [CrossRef] [Green Version]
  35. NVIDIA Corporation NVIDIA Tesla V100 GPU Architecture. Available online: https://images.nvidia.com/content/volta-architecture/pdf/volta-architecture-whitepaper.pdf (accessed on 7 December 2019).
  36. Moon, B.; Jagadish, H.V.; Faloutsos, C.; Saltz, J.H. Analysis of the Clustering Properties of Hilbert Space-filling Curve. IEEE Trans. Knowl. Data Eng. 2001, 13, 124–141. [Google Scholar] [CrossRef] [Green Version]
  37. Nvidia CUB. Available online: https://nvlabs.github.io/cub/ (accessed on 7 December 2019).
  38. NVIDIA Corporation NVIDIA Turing GPU Architecture. Available online: https://www.nvidia.com/content/dam/en-zz/Solutions/design-visualization/technologies/turing-architecture/NVIDIA-Turing-Architecture-Whitepaper.pdf (accessed on 2 February 2020).
  39. Cea, L.; Bermúdez, M.; Puertas, J. Uncertainty and sensitivity analysis of a depth-averaged water quality model for evaluation of Escherichia Coli concentration in shallow estuaries. Environ. Modell. Softw. 2011, 26, 1526–1539. [Google Scholar] [CrossRef]
  40. García-Barcina, J.M.; Oteiza, M.; de la Sota, A. Modelling the faecal coliform concentrations in the Bilbao estuary. Hydrobiologia 2002, 475, 213–219. [Google Scholar] [CrossRef]
  41. Kashefipour, S.M.; Lin, B.; Falconer, R.A. Modelling the fate of faecal indicators in a coastal basin. Water Res. 2006, 40, 1413–1425. [Google Scholar] [CrossRef] [PubMed]
  42. Manache, G.; Melching, C.S.; Lanyon, R. Calibration of a continuous simulation fecal coliform model based on historical data analysis. J. Environ. Eng. 2007, 133, 681–691. [Google Scholar] [CrossRef]
  43. Bode, A.; Álvarez-Ossorio, M.T.; González, N.; Lorenzo, J.; Rodríguez, C.; Varela, M.; Varela, M.M. Seasonal variability of plankton blooms in the Ria de Ferrol (NW Spain): II. Plankton abundance, composition and biomass. Estuarine Coast. Shelf Sci. 2005, 63, 285–300. [Google Scholar] [CrossRef]
  44. Anta Álvarez, J.; Bermúdez, M.; Cea, L.; Suárez, J.; Ures, P.; Puertas, J. Modelización de los impactos por DSU en el río Miño (Lugo). Ingeniería del agua 2015, 19, 105–116. [Google Scholar] [CrossRef]
  45. Stumm, W.; Morgan, J. Aquatic Chemistry: Chemical Equilibria and Rates in Natural Waters; Wiley: New York, NY, USA, 1996. [Google Scholar]
  46. Piñeiro, J.; Maestro, I.; Aguirre, F.; Ures, P.; Torres, D.; Anta, J.; Puertas, J.; Suárez, J. Análisis del funcionamiento de un depósito-aliviadero en el sistema de saneamiento unitario en la aglomeración de Lugo. In Proceedings of the Actas de las II Jornadas de Ingeniería del Agua, Barcelona, Spain, 5–6 October 2011. [Google Scholar]
Figure 1. Model structure. Solid arrows indicate interactions between variables due to different biochemical processes.
Figure 1. Model structure. Solid arrows indicate interactions between variables due to different biochemical processes.
Water 12 00413 g001
Figure 2. Flowchart of the Iber+ implementation.
Figure 2. Flowchart of the Iber+ implementation.
Water 12 00413 g002
Figure 3. Bathymetry of the Ferrol estuary (NW, Spain). The sewage discharges are located at V1 and V2, whilst the control points are placed at P1, P2 and P3.
Figure 3. Bathymetry of the Ferrol estuary (NW, Spain). The sewage discharges are located at V1 and V2, whilst the control points are placed at P1, P2 and P3.
Water 12 00413 g003
Figure 4. Comparison of Escherichia coli concentrations at different time steps obtained with Iber (CPU) and Iber+ (GPU).
Figure 4. Comparison of Escherichia coli concentrations at different time steps obtained with Iber (CPU) and Iber+ (GPU).
Water 12 00413 g004
Figure 5. Time series of E. coli concentration at control points P1 (a), P2 (b) and P3 (c), for the first test case.
Figure 5. Time series of E. coli concentration at control points P1 (a), P2 (b) and P3 (c), for the first test case.
Water 12 00413 g005
Figure 6. Bathymetry of the estuary of A Coruña (NW, Spain). The sewage spills are located at V1, V2, V3 and V4. The control points are placed at P1, P2 and P3.
Figure 6. Bathymetry of the estuary of A Coruña (NW, Spain). The sewage spills are located at V1, V2, V3 and V4. The control points are placed at P1, P2 and P3.
Water 12 00413 g006
Figure 7. Time series of carbonaceous biochemical oxygen demand (CBOD) for Test 2 at control points P1 (a), P2 (b) and P3 (c).
Figure 7. Time series of carbonaceous biochemical oxygen demand (CBOD) for Test 2 at control points P1 (a), P2 (b) and P3 (c).
Water 12 00413 g007
Figure 8. Time series of DO for Test 2 at control points P1 (a), P2 (b) and 3 (c).
Figure 8. Time series of DO for Test 2 at control points P1 (a), P2 (b) and 3 (c).
Water 12 00413 g008
Figure 9. Bathymetry of the Miño River as it passes through the city of Lugo (NW, Spain). The sewage spills are located at points V1, V2 and V3. The control points are located at P1, P2, P3 and P4.
Figure 9. Bathymetry of the Miño River as it passes through the city of Lugo (NW, Spain). The sewage spills are located at points V1, V2 and V3. The control points are located at P1, P2, P3 and P4.
Water 12 00413 g009
Figure 10. Time series sampled at control point P3 for dissolved oxygen (a), carbonaceous biochemical oxygen demand (b), organic nitrogen (c), ammoniacal nitrogen (d) and nitrite-nitrate nitrogen (e).
Figure 10. Time series sampled at control point P3 for dissolved oxygen (a), carbonaceous biochemical oxygen demand (b), organic nitrogen (c), ammoniacal nitrogen (d) and nitrite-nitrate nitrogen (e).
Water 12 00413 g010
Figure 11. Bathymetry of the river reach in Test Case 4. The sewage spill of the Wastewater Treatment Plant (WWTP) is located at point V1. Control points are located at P1, P2, P3 and P4.
Figure 11. Bathymetry of the river reach in Test Case 4. The sewage spill of the Wastewater Treatment Plant (WWTP) is located at point V1. Control points are located at P1, P2, P3 and P4.
Water 12 00413 g011
Figure 12. Detail of the concentration of CBOD downstream the sewage spill at t = 60 h using both Iber and Iber+ implementations.
Figure 12. Detail of the concentration of CBOD downstream the sewage spill at t = 60 h using both Iber and Iber+ implementations.
Water 12 00413 g012
Figure 13. Time series for the analyzed species in Test 4, sampled at control point P3.
Figure 13. Time series for the analyzed species in Test 4, sampled at control point P3.
Water 12 00413 g013
Table 1. Performance measurements obtained for the first test case (estuary of Ferrol).
Table 1. Performance measurements obtained for the first test case (estuary of Ferrol).
ModelRun Time (s)Time per Step (ms)Millions of Cells per SecondSpeedup vs. Iber
Iber60,54566.82.21
Iber+ GPU3350.4400.0181
Table 2. Performance measurements obtained for the second test case (estuary of A Coruña).
Table 2. Performance measurements obtained for the second test case (estuary of A Coruña).
ModelRun Time (s)Time per Step (ms)Millions of Cells per SecondSpeedup vs. Iber
Iber31,61528.81.81
Iber+ GPU5220.5107.161
Table 3. Performance measurements obtained for the third test case (Miño River, NW Spain).
Table 3. Performance measurements obtained for the third test case (Miño River, NW Spain).
ModelRun Time (s)Time per Step (ms)Millions of Cells per SecondSpeedup vs. Iber
Iber30548.61.01
Iber+ GPU1050.330.229
Table 4. Performance measurements obtained for the fourth test case.
Table 4. Performance measurements obtained for the fourth test case.
ModelRun Time (s)Time per Step (ms)Millions of Cells per SecondSpeedup vs. Iber
Iber44,367146.50.61
Iber+ GPU4821.656.992

Share and Cite

MDPI and ACS Style

García-Feal, O.; Cea, L.; González-Cao, J.; Domínguez, J.M.; Gómez-Gesteira, M. IberWQ: A GPU Accelerated Tool for 2D Water Quality Modeling in Rivers and Estuaries. Water 2020, 12, 413. https://doi.org/10.3390/w12020413

AMA Style

García-Feal O, Cea L, González-Cao J, Domínguez JM, Gómez-Gesteira M. IberWQ: A GPU Accelerated Tool for 2D Water Quality Modeling in Rivers and Estuaries. Water. 2020; 12(2):413. https://doi.org/10.3390/w12020413

Chicago/Turabian Style

García-Feal, Orlando, Luis Cea, José González-Cao, José Manuel Domínguez, and Moncho Gómez-Gesteira. 2020. "IberWQ: A GPU Accelerated Tool for 2D Water Quality Modeling in Rivers and Estuaries" Water 12, no. 2: 413. https://doi.org/10.3390/w12020413

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop