Abstract

In this study, a 2D shallow water flow solver integrated with a water quality model is presented. The interaction between the main water quality constituents included is based on the Water Quality Analysis Simulation Program. Efficiency is achieved by computing with a combination of a Central Processing Unit (CPU) and a Graphics Processing Unit (GPU) device. This technique is intended to provide robust and accurate simulations with high computation speedups with respect to a single-core CPU in real events. The proposed numerical model is evaluated in cases that include the transport and reaction of water quality components over irregular bed topography and dry–wet fronts, verifying that the numerical solution in these situations conserves the required properties (C-property and positivity). The model can operate in any steady or unsteady form allowing an efficient assessment of the environmental impact of water flows. The field data from an unsteady river reach test case are used to show that the model is capable of predicting the measured temporal distribution of dissolved oxygen and water temperature, proving the robustness and computational efficiency of the model, even in the presence of noisy signals such as wind speed.

HIGHLIGHTS

  • Implementation of a two-dimensional water quality model with high-performance computing.

  • Conservative and well-balanced scheme in the presence of bed and friction source terms.

  • Robust scheme in the presence of unsteady flows and pollutants.

  • Real-time prediction of contaminants even in the presence of wet–dry fronts.

INTRODUCTION

Many of today's water pollution problems require continuous and detailed attention. This alteration of the environment, directly or indirectly caused by a project or activity in a given area in the past, has been characterized by models that evaluate the environmental impact in a watershed. The complex interaction of species present in surface water bodies opened the discussion of which processes should be incorporated into these water quality models.

A wide range of water quality simulation models can be found in the literature. Some of them such as AQUATOX (Park & Clough 2014), SIMulation of CATchments (SIMCAT) (Warn 1987), Temporal Overall Model for Catchments (TOMCAT) (Jamieson & Fedra 1996), QUAL2Kw (Wells 2019), QUASAR (Whitehead et al. 1997), Simplified Method Program for Multiple Discharge (Multi-SMP) (Hosseinipour & Martin 1992), among others, are able to predict the evolution of some chemical species in one-dimensional (1D) steady-state flows. More ambitious models evaluate the concentration of water quality components in the framework of one-dimensional (1D) and two-dimensional (2D) unsteady flows. Examples of these in 1D are CE-QUAL-RIV1 (Martin et al. 2002), MIKE-21 (Warren & Bach 1992) and the 1D unsteady river flow model recently proposed by Gordillo et al. (2019) based on Water Quality Analysis Simulation Program (WASP) (Wool et al. 2006). IberWQ (Cea et al. 2016) is an example of 2D hydraulic model that includes a water quality model. Both Cea et al. (2016) in 2D and Gordillo et al. (2019) in 1D are based on a coupled discretization of the hydraulic and water quality parts that ensure positivity by considering the effects of bed slope and friction source terms. On this subject, different techniques have been developed to control the two main basic properties that any numerical scheme with solute transport must comply with: C-property and positivity (Cea & Vázquez-Cendón 2012; Murillo et al. 2012; Liu 2019).

On the other hand, high-performance computing (HPC) has recently emerged as a necessity when trying to develop flood-hazard maps and evaluate the environmental impacts produced by floods. The reason is that complex models require excessively long computational times in practical cases to provide predictions on large domains. In Vacondio et al. (2014), the simulation of real flood events was carried out by means of an explicit second-order finite volume scheme resolved on a Cartesian grid and the speedup in the computation was achieved using Graphics Processing Units (GPU). A 2D shallow water model was also presented in Echeverribar et al. (2019) and Lacasta et al. (2014) to predict real-time flooding events using GPU on triangular unstructured grids. The explicit first-order scheme demonstrated absolute control over the numerical instabilities present in scenarios with complex geometries. Smith & Liang (2013) also reported a significant computational acceleration for realistic cases using a second-order finite volume Godunov-type scheme implemented on GPU.

There is a lack of water quality models that take advantage of the efficiency of GPU devices. In this context, the novelty of this work consists in developing a GPU implementation of a 2D numerical model that not only provides flood-hazard maps but also allows the water quality characterization. To achieve this goal, a water quality model based on 10 transported variables (WASP) that take into account the main interactions that affect dissolved oxygen levels and temperature of a given water volume is coupled to the 2D hydraulic flow model, enabling the prediction of water quality in flood events with sufficient accuracy in a moderate computational time. The system of non-linear hyperbolic partial differential equations governing the water flow and the chemical species transport is discretized by a finite volume Godunov-type scheme. The resulting set of computations is efficiently carried out using a GPU device.

The paper is organized as follows: section ‘Governing equations’ is a short introduction to the governing equations; in particular, ‘2D hydraulic model’ provides information about the shallow water equations (SWE) while ‘2D water quality model’ gives a description of the transport equation with the rate processes included in the model. The numerical scheme used for both the water flow equations and chemical species transport equations is presented in the section ‘Numerical scheme and HPC implementation’. Next, section ‘HPC implementation of the water quality model’ contains the HPC implementation for the sequential and intensive parts of the model. In section ‘Numerical experiments’, the scheme is verified with several academic tests and one real case. Finally, summary and concluding remarks are presented in ‘Conclusions’.

GOVERNING EQUATIONS

2D hydraulic model

The 2D SWE allows to describe free surface flows (Fernández-Pato & García-Navarro 2018). This system of equations represents depth-averaged mass and momentum conservation along the main x- and y-directions of the flow:
formula
(1)
formula
(2)
formula
(3)
where h represents the water depth, and are the unit discharges, with u and v the depth-integrated components of the velocity vector in the x and y directions, respectively, g stands for the acceleration due to gravity, are the bed slopes and are the friction slopes here formulated using Manning's law (Chow 1964).

2D water quality model

Depth-averaged transport equations for 10 scalar variables are used in the present work to evaluate the physical, chemical and biological processes of interest in rivers. The equations for the WASP water quality state variables considered – ammonium nitrogen (), nitrate nitrogen (), inorganic (ortho) phosphorus (IP), phytoplankton (A), Carbonaceous Biological Oxygen Demand (C-BOD), dissolved oxygen (DO), organic nitrogen (ON), organic phosphorus (OP), temperature (T) and total coliform bacteria (TC) can be written in conservative form as:
formula
(4)
where is the depth-averaged concentration of the s state variable (). and are the diffusion–dispersion coefficients in the x and y directions, respectively, is the rate of change in the substance concentration due to the physical, chemical and biological processes as a function of each concentration and the p model parameters. The description of these processes of degradation of organic material, growth and death of algae, nitrification, hydrolysis of organic nitrogen and phosphorus, re-aeration, sedimentation of algae, phosphorus and nitrogen, sediment uptake of oxygen are detailed in Thomann & Mueller (1987); Cox (2003) and Gordillo et al. (2019). The set of state variables with their conversion processes form a matrix of size where is the number of processes and is the number of chemical species. The formulation of the individual processes for each rate parameter related to each component is described in Gordillo et al. (2019).

Finally, are point and non-point sources of each species. It is worth noting that, since most of the kinetic coefficients depend on the temperature, the proposed model calculates the water temperature among the set of water quality variables (). The equation for heat exchange is defined in Gordillo et al. (2019). This expression uses atmospheric radiation and sensitive heat fluxes in the heat balance (Edinger et al. 1968).

NUMERICAL SCHEME AND HPC IMPLEMENTATION

This section describes the process followed to guarantee the conservation property and to avoid non-physical results of solute transport as well as its HPC implementation. The details of the numerical scheme can be found in Murillo & García-Navarro (2012) and Murillo et al. (2009). Therefore, they are not described again in the present work. However, the process to achieve a fully conservative method with the implementation in GPU of the transport of chemical species especially in real-scale scenarios (hundreds of ) is described below.

As the interest of this work is focused on the evolution of the concentrations of each chemical species within complex hydrodynamic systems, the system (1)–(3) with active solute (4) can be formulated in a single system of equations:
formula
(5)
with:
formula
(6)
where accounts for the conserved variables, are the fluxes, are the source terms, are the solute process rates, is the diffusion term and is an empirical dispersion-diffusion matrix.

Following previous work (Morales-Hernández et al. 2019; Murillo et al. 2009), the numerical solution of the hydraulic equations and that of the transport equations are performed sequentially.

Discretization of the hydrodynamic system

First, the hydrodynamic variables are solved by means of a finite volume method based on an augmented Roe's solver.

To update the flow variables in cell i of area , from to , an explicit first-order upwind scheme is formulated by García-Navarro et al. (2019):
formula
(7)
where is the number of edges in cell i, k denotes the cell edges, is the edge length and m corresponds to the eigenvalues counter. The expression (7) is written in terms of the incoming part of the Jacobian eigenvalues,
formula
(8)
the Jacobian eigenvectors , and the projection of fluxes and source terms on that basis, .

Being (7) an explicit scheme, the time step size is limited for stability reasons by the cell size and the eigenvalues of the system. The allowed time step is dynamically calculated every time step as detailed in Morales-Hernández et al. (2019). Additionally, the proposed scheme may produce non-physical values on dry–wet fronts, thus an alternative is to reduce the time step size. Another option is to control the intermediate states of the conserved variables of the flow in order to guarantee a non-negative solution. The complete formulation can be found in Murillo & García-Navarro (2010). These techniques ensure that the scheme is conservative and maintains water depth positivity on irregular beds.

Discretization of the transport equations

Following the same philosophy, a first-order upwind scheme can be formulated for the transport of conservative and non-conservative solutes. Previous studies (Murillo et al. 2009; Cea & Vázquez-Cendón 2012) have shown that the scheme is able to guarantee conservation. However, it may give rise to unrealistic values when dealing with problems that involve bed variations and transient flows over wet/dry fronts, leading to unbounded values in the final solute concentration (Morales-Hernández et al. 2019). To overcome this problem, the Courant–Friedrichs–Lewy (CFL) condition is frequently implemented, severely restricting the time step size and causing computational inefficiency.

To avoid the inconveniences described above, a scheme with a redistribution in solute fluxes is presented. This technique consists of defining a numerical flux directly related to the Roe's linearization. This allows the solution to decouple the solute transport equation (4) of the hydrodynamic system (1)–(3). The present work is concerned with transport processes dominated by the convection and reaction terms. The influence of diffusion will be neglected and therefore from now on.

The process of updating the set of conserved variables is carried out in three steps:

  • Step (I). Compute the updating value of the hydrodynamic variables , through Equation (7).

  • Step (II). Compute the updating value of the advection part in the set of transported variables as:
    formula
    (9)
    with
    formula
    (10)
    and the numerical solute concentration at edge k separating cells i and :
    formula
    (11)
  • Step (III). Compute the updating value of the reaction R and source f terms as:
    formula
    (12)

With this redistribution of fluxes considering both the hydraulic variations and the biological processes of the constituent, a perfect preservation of the solute positivity and a non-oscillatory solution are guaranteed.

HPC implementation of the water quality model

Generally, the resolution of 2D models has a high computational burden that depends on the physical processes that are considered, on the extension of the spatial domain and on the duration of the simulated event. There are multiple alternatives that accelerate the calculation of the operations present in the simulation process. One of them is the use of GPU devices. Thus, using Compute Unified Device Architecture (CUDA), a kernel containing all the reaction processes is implemented to run in parallel within the GPU.

The simulation tool presented is structured on a heterogeneous environment, that is, the model contains an architecture in the Central Processing Unit (CPU) and another in the GPU. The sequential algorithm is solved in the CPU, while the intensive calculation part is sent to the GPU. The first step of the sequential part is to process all the reading functions of the initial conditions, boundary conditions and coefficients/parameters of the hydrodynamic and water quality part. The next step is to generate the mesh. Then, memory is allocated for all the necessary structures and arrays. Finally, the solute variables and reaction matrix are vectorized and data for both the CPU and the GPU are initialized.

The vectorization of both the solute variables and the reaction matrix is explained here. Assuming, for instance, that the total number of solutes is and the number of processes is , this results in a 5 × 2 reaction matrix. In addition, if the number of cells is , the solute variables vector will be built in groups of for each solute (see Figure 1).

Figure 1

Solute vectorization.

Figure 1

Solute vectorization.

On the other hand, the reaction matrix vector is constituted by groups of reaction processes (5 in this example) for each solute. Figure 2 illustrates the implementation chosen here to transform the reaction matrix into a vector.

Figure 2

Vectorization of the reaction matrix.

Figure 2

Vectorization of the reaction matrix.

A first step is required for transferring the information from the CPU to the GPU. The operations related to the numerical scheme run on the GPU. The computation of both the hydrodynamic variables and solute variables at all grid cell requires long periods of time to find the final solution. For this reason, it is of interest to perform this set of operations on the most efficient computational framework. The details of the syntax used in the functions on the GPU device, the aspects to be taken into account in the implementation of the numerical models and the general structure of the main code loop can be found in Lacasta et al. (2014). In order to represent the results, the information of the updated variables must be transferred from the GPU to the CPU. With the information already in CPU, the results are finally written.

NUMERICAL EXPERIMENTS

This section presents some numerical examples that show the strengths of the proposed tool. The first two cases have been selected to demonstrate the accuracy and robustness of the scheme for the pure advective transport of a single solute in cases including wet/dry fronts in an analytical geometry. The third academic case involves the transport of nine solute variables considering their reaction/decomposition processes over the same idealized geometry. The last case aims to check the computational speed in predicting water flow and water quality evolution in a real river reach event.

Case 1: Quiescent flow on irregular bottom

To demonstrate that the proposed scheme guarantees the C-property and solute positivity in the presence of wet/dry interfaces, the following example is presented. A rectangular domain 500 m long in the x direction and 300 m long in the y direction is considered. The bed roughness is characterized by a uniform Manning coefficient . The bed topography includes three conical islands centered at (, ), (, ) and (, ), respectively. Therefore, the bed level z is defined as:
formula
(13)
being and , some geometrical constants. In this example, and for the first two islands and and for the third island. Figure 3 illustrates the topography configuration and the triangular unstructured mesh used.
Figure 3

Topography and computational mesh for cases 1, 2 and 3.

Figure 3

Topography and computational mesh for cases 1, 2 and 3.

The initial flow conditions are a uniform water level and zero velocity:
formula
The initial condition for the single and non-reactive solute is a uniform concentration:
formula
(14)

As for the boundary conditions, solid wall is assumed on the four sides of the domain. The computational domain is discretized in 23,306 triangular cells and CFL = 0.95. Figure 4 shows the water surface elevation at both as a surface (a) and as a longitudinal profile along the central line of the x axis (b). Figure 4(b) also shows the concentration profile on the same line. As observed, the numerical results of both the water surface elevation and the solute concentration remain constant. No oscillations appear even in the presence of partially submerged and completely submerged islands. The method is able to reproduce the correct values on both the wet and dry sides.

Figure 4

3D view of the quiescent water surface (a) and water level, bed level and concentration profile of the solution along the central axis (b) at . Case 1.

Figure 4

3D view of the quiescent water surface (a) and water level, bed level and concentration profile of the solution along the central axis (b) at . Case 1.

Case 2: Dambreak with scalar transport on dry irregular bottom

To check the robustness of the method in the presence of unsteady flow and source terms, a dambreak with conservative solute transport is simulated. The bed topography (13) and the bed roughness are the same as in Case 1.

The initial conditions for the flow are zero velocity and discontinuous water surface level at . An initial water level is assumed on the left. Dry bed conditions are assumed on the right part of the discontinuity:
formula
(15)
formula
The initial conditions for the single non-reactive solute are given by:
formula
(16)

Solid wall boundary conditions are imposed at all domain boundary edges. Figure 5(a) and 5(c) shows the water depth distribution at and , respectively. As shown, the solution is free of oscillations on wet/dry fronts in any of the simulation times presented, even in the presence of the irregular bed level. Figure 5(b) and 5(d) shows the distribution of concentration at the same times. The proposed method keeps not only the positivity of the concentration but also its uniform value perfectly when it is transported in a body of water with irregular topography. The concentration is not compromised even in the presence of the wet/dry interface.

Figure 5

Contour of water depth (left) and concentration (right) at (upper) and (lower). Case 2.

Figure 5

Contour of water depth (left) and concentration (right) at (upper) and (lower). Case 2.

Case 3: Dambreak with reactive solute transport on irregular bottom

The aim of this academic case is to evaluate the performance of the method when including reactive solutes in a dambreak flow over irregular topography. This scenario is carried out with the same bed level function as the previous cases. The model in this case is formed by the hydrodynamic equations and the water quality transport equations (4) with .

The initial conditions for the water flow are zero velocity and discontinuous water surface level at where is assumed on the left and is assumed on the right:
formula
(17)
formula

Table 1 shows the initial conditions for the solute variables expressed in and those of the water temperature expressed in °C.

Table 1

Initial conditions of non-conserved variables (Case 3)

Variableif if
 6.0 
 0.1 0.05 
 5.0 0.02 
 0.065 0.025 
 30.0 2.0 
 6.0 9.25 
 6.0 
 0.1 0.01 
 23.0 12.0 
Variableif if
 6.0 
 0.1 0.05 
 5.0 0.02 
 0.065 0.025 
 30.0 2.0 
 6.0 9.25 
 6.0 
 0.1 0.01 
 23.0 12.0 

All the water quality model parameters are chosen according to their default value following Bowie et al. (1985). In this case, real records of the meteorological variables are used in order to consider the temporal variations throughout the day, thus avoiding raising the water temperature to non-real values (Figure 12). The domain is assumed closed as in Cases 1 and 2.

Under these conditions, a sensitivity test of the results to the mesh is carried out with a simulation time of one day. The first mesh used to simulate this example consists of 23,306 triangular cells (referred to in this study as Mesh 1). The second mesh consists of 46,618 cells (Mesh 2), and the third mesh includes 93,109 cells (Mesh 3). Figure 6 shows the temporal evolution of the DO concentration and water temperature at point P (x = 200 m, y = 150 m) computed on the three meshes. There are significant differences between Mesh 1 and Mesh 2; however, between Mesh 2 and Mesh 3, the results are very similar. It is verified that, in two-dimensional cases, an appropriate refinement of the computational mesh is required so that the results do not present severe variations. The computational times with the GPU device for 1-day simulation time for the three meshes studied were 0.85, 1.48 and 3.58 h, respectively (see Table 2).

Table 2

Case 3: summary of meshes and GPU times

MeshCellsGPU time
23,306 0.85 h 
46,618 1.48 h 
91,101 3.58 h 
MeshCellsGPU time
23,306 0.85 h 
46,618 1.48 h 
91,101 3.58 h 
Figure 6

Time evolution of DO concentration (a) and water temperature (b) at P(200,150) with different meshes. Case 3.

Figure 6

Time evolution of DO concentration (a) and water temperature (b) at P(200,150) with different meshes. Case 3.

Once the mesh convergence analysis has been carried out, Mesh 2 is chosen to simulate the complete scenario of five days. Figure 7 represents the variables and at days. The following features can be observed:

  • The concentration distribution of the reactive solutes is nil when , showing a non-oscillatory behavior even in the presence of bed irregularity, verifying correctly the C-property.

  • Under these conditions, the active decomposition zone is maintained in the first 100 m.

  • Nutrients do not pose a problem in the bloom of plants (Figure 7(b)). The discharge of C-BOD has important effects on the concentration level of DO (Figure 7(d)). The water temperature slightly increases around the third island due to the minimal water depth in that area (Figure 7(e)).

Figure 7

Water depth (a) and concentrations of algae (b), C-BOD (c), DO (d) and water temperature (e) at . Case 3.

Figure 7

Water depth (a) and concentrations of algae (b), C-BOD (c), DO (d) and water temperature (e) at . Case 3.

Representing the most relevant variables after day (Figure 8), the following can be seen:

  • The decomposition zone is located mainly near the boundaries.

  • The algae growth is controlled (Figure 8(a)).

  • The concentration of C-BOD decreases significantly, causing a sag in DO concentration. A value of approximately can be observed in some areas (Figure 8(c)).

  • The maximum temperature is concentrated in the decomposition zone with values between 15 and 18.5 °C approximately (Figure 8(d)).

Figure 9 plots the numerical results of some state variables present in the eutrophication model after 5 days of simulation. Figure 9(a) again reflects that algae blooming is controlled. Regarding concentration of C-BOD (see Figure 9(b)), it reaches an approximate value of 3.5 g = m3 which, together with the slight increase in temperature, generates a considerable consumption of DO during this period of time (Figure 9(c)).

Figure 8

Concentrations of non-conservative solutes of algae (a), C-BOD (b), DO (c) and temperature (d) at day. Case 3.

Figure 8

Concentrations of non-conservative solutes of algae (a), C-BOD (b), DO (c) and temperature (d) at day. Case 3.

Figure 9

Concentrations of non-conservative solutes of algae (a), C-BOD (b), DO (c) and temperature (d) at day. Case 3.

Figure 9

Concentrations of non-conservative solutes of algae (a), C-BOD (b), DO (c) and temperature (d) at day. Case 3.

Case 4: Ebro river reach

A reach of the Ebro River is considered in the study. It is located between the cities of Alagón and Zaragoza with an approximate length of 40 km (see Figure 10). This river reach is used to test the simulation model involving the three hydrodynamic equations and the transport equations corresponding to the ammonium nitrogen (), nitrate nitrogen (), carbonaceous biological oxygen demand (C-BOD), dissolved oxygen (DO), organic nitrogen (ON) and temperature (T).

Figure 10

Topography mesh for Case 4.

Figure 10

Topography mesh for Case 4.

Mesh construction

A Digital Terrain Model (DTM) of resolution 5 × 5 m is used as input data to define the topography of the study region. However, since this information uses LIDAR technology, the bathymetry of the river bottom is not completely captured due to the water surface. To overcome this problem, additional cross-sections provided by the river basin authorities (Confederación Hidrográfica del Ebro, CHE) were used.

The roughness coefficient is estimated from the land use (Sistema de Información de Ocupación del Suelo de España (SIOSE)) at a reference scale of 1:25,000. Soil uses and other factors such as the irregularity of the land are analyzed to obtain a map of roughness values for the different types of soil.

With the previous information, an unstructured triangular mesh was created using 169,918 triangular cells of variable size, more refined in regions requiring a higher level of detail such as the river bed or thin levees. The length of the edges of the cells ranged from 13.0 m (near the river bed) to 159.0 m (limit of the domain).

Initial and boundary conditions

To reproduce and evaluate a transient event, information on hydrodynamic and quality parameters is required. This information is imposed at the beginning of the calculation () as initial condition and at the domain boundaries. Through the monitoring stations managed by the CHE, information is obtained with a periodicity of 15 min for the discharge flow rate, water temperature and dissolved oxygen concentration (the latter measures less frequently on some occasions). In the present case, a 7-day event corresponding to a 2012 flood event (December 18–25) is simulated. The measured discharge hydrograph as well as the measured DO concentrations and water temperatures at Alagón were used as time variable inlet boundary conditions. They are plotted in Figure 11. For the rest of the water quality components, the inlet boundary value was assumed constant. On the other hand, external heat contributions, such as air temperature, relative humidity, wind speed and solar radiation, were obtained from the Sociedad Aragonesa de Gestión Agroambiental (SARGA). Meteorological data are recorded every hour so that when there are no records of the meteorological and chemical variables a linear interpolation is performed. Figure 12 reports the meteorological variables measured in Zaragoza (nearest measurement station). Finally, the rest of the water quality parameters required by the WASP model has been taken from the existing literature (Bowie et al. 1985).

Figure 11

Upstream boundary condition of discharge Q (left), DO concentration and water temperature (right). Case 4.

Figure 11

Upstream boundary condition of discharge Q (left), DO concentration and water temperature (right). Case 4.

Figure 12

Meteorological data for the water temperature equation. Case 4.

Figure 12

Meteorological data for the water temperature equation. Case 4.

The initial conditions for the simulation were set as the steady state computed in a previous run of the model using steady inlet boundary conditions corresponding to the first value of the inlet variables in the temporal series.

Results

Analyzing the variable h in t = 25 h and t = 168 h (see Figure 13(a) and 13(b)) in a section of the river, we can see areas that are slightly flooded when transporting flows around . This event involves the transport of some chemical species that cause a decrease in the concentration of DO. In this case, for example, from t = 25 h to t = 168 h, the DO concentration varies from 11.07 to . As follows from the analysis of the numerical results presented in Figure 13(c) at t = 25 h, there are DO concentrations below due to the wet–dry front. These values are not significant in relation to the reach of the river studied, as they appear in a punctual way on the banks of the domain. At t = 168 h, the DO concentration levels are more uniform, an average value of is observed at the end of the simulation (Figure 13(d)).

Figure 13

Snapshots of water elevation (a,b), DO concentration (c,d) and water temperature (e and f) at t = 25 h (left) and t = 168 h (right). Case 4.

Figure 13

Snapshots of water elevation (a,b), DO concentration (c,d) and water temperature (e and f) at t = 25 h (left) and t = 168 h (right). Case 4.

Regarding the variable water temperature, a variation from upstream to downstream is observed from 8.3 to 7.87 °C at (Figure 13(e)). Finally, at , we can see the distribution of temperatures that go from 9.49 to 9.28 °C in the downstream part of the reach (Figure 13(f)).

Figure 14 illustrates the comparison between the temporal evolution of the simulated flow rate with the measurements at Zaragoza. The simulated values of discharge in certain periods of time are superior to the recorded data. This could be due to the Manning coefficient estimation process since, although it is evaluated according to the land use, there exists a lack of accuracy on the soil measurements and on the estimation process itself.

Figure 14

Temporal evolution of the flow rate at the Zaragoza monitoring station. Case 4.

Figure 14

Temporal evolution of the flow rate at the Zaragoza monitoring station. Case 4.

Figure 15(a) and 15(b) compare the measured and simulated DO concentrations and water temperatures recorded during the 7 days. Marked differences are found during the first day, which matches when the flow rate transported by the river bed is around .

Figure 15

Temporal evolution of DO concentration (a) and water temperature (b) at the Zaragoza monitoring station. Case 4.

Figure 15

Temporal evolution of DO concentration (a) and water temperature (b) at the Zaragoza monitoring station. Case 4.

Except for these differences, the results between the values measured and simulated with the proposed model follow the same trends. In addition, the calculated distributions do not show any instability even with high noise meteorological signals such as wind speed.

A quantitative analysis of the DO and water temperature variables gives a root mean square error (RMSE) of 0.338 and 0.549, respectively. There is great difficulty in reproducing in more detail the processes that occur in a reach of the river when there is a lack of regularity in the water quality monitoring by the agencies in charge of management. Although there are techniques that allow the recovery of information, it does not become the final solution because all these techniques require more periodic information in the measurement of state variables present in the quality model (Gordillo et al. 2020).

Since the computational load in this case is very high due to 2D hydrodynamic and water quality computation, an efficiency analysis is carried out. When the hydrodynamic and water quality models are both solved with a single CPU processor (Intel (R) Xeon (R) CPU E5-2697 v3@2.6 GHz), the calculation time is 108 h, whereas the same scenario requires 3.75 h when the model is solved on a Tesla K80 GPU device. Therefore, the gain in efficiency when using GPU computing in comparison with CPU is 28.8×. On the other hand, in order to measure the extra time required for the water quality algorithms, an additional simulation has been performed considering only the hydrodynamic model (without any water quality calculations), obtaining a simulation time of 1.76 h. This highlights the great complexity of the solute transport and reaction algorithms.

CONCLUSIONS

The attention paid at water quality values is increasing not only for environmental reasons but also because of its importance as a fundamental human resource. This makes modeling techniques a strong alternative for the prediction of complex phenomena, such as water quality evolution during flood events, and its value as a support tool for decision makers. In this context, an efficient coupling between hydrodynamic and water quality tools is needed to provide accurate and fast results. In this study, a robust and accurate 2D numerical model based on HPC is successfully developed to simulate the SWE with transport of chemical species and their interactions. The proposed numerical scheme is based on a finite volume method that is solved with a first-order upwind scheme for unstructured triangular meshes. A fundamental analysis has been carried out to guarantee the proper performance of the model when it is applied to theoretical cases, such as water at rest or a dambreak. It was observed that the method can be calculated without presenting any non-physical state or instability of the solute (passive and reactive) in both, at rest and with a rapid variation of the flow (dambreak). Additionally, positivity in all the solutes transport is ensured by the numerical method in the presence of wet/dry fronts with variable bed without having to resort to time reduction techniques beyond the condition of CFL. In this work, the model has been applied to a realistic case not only to analyze the computational behavior and speedups but also to compare the obtained results with measured data. A stretch of the Ebro River (Spain) has been selected to reproduce a historical event with the water quality module and match the results with measurements. Despite the lack and low quality of the measurements, a good agreement is achieved in terms of DO concentration and temperature evolution, demonstrating the capability of the model to reproduce real situations. Once the model was demonstrated to be robust and accurate enough, an analysis of computational times was carried out, comparing simulation times with GPU devices. Simulations of the real scenario showed that there is a 28.8× acceleration with a GPU device over CPU. This reflects the importance of implementing efficient techniques that reduce the time in predicting both hydraulic and chemical variables. In this way, the use of HPC techniques provides a noticeable gain in the computational performances when simulating reactive solute transport in real scenarios. This computational strategy allows the development of efficient tools for the environmental impact control in real time. Additionally, it is worth mentioning that every real application requires a proper calibration of the model parameters and this process would also take advantage of the HPC techniques due to the high amount of simulations to be carried out.

ACKNOWLEDGEMENTS

This work was funded by the Spanish Ministry of Science and Innovation under the research project PGC2018-094341-B-I00. Additionally, Mario Morales-Hernández was partially supported by the U.S. Air Force Numerical Weather Modeling Program.

SUPPLEMENTARY DATA

The data to reproduce the cases will be made available on request.

DATA AVAILABILITY STATEMENT

All relevant data are included in the paper or its Supplementary Information.

REFERENCES

REFERENCES
Bowie
G. L.
Mills
W. B.
Porcella
D. B.
Campbell
C. L.
Pagenkopf
J. R.
Rupp
G. L.
Johnson
K. M.
Chan
P.
Gherini
S. A.
Chamberlin
C. E.
1985
Rates, constants, and kinetics formulations in surface water quality modeling
.
EPA
600
,
3
85
.
Cea
L.
Bermúdez
M.
Puertas
J.
Blade
E.
Corestein
G.
Escolano
E.
Conde
A.
Bockelmann-Evans
B.
Ahmadian
R.
2016
IberWQ: new simulation tool for 2D water quality modelling in rivers and shallow estuaries
.
Journal of Hydroinformatics
18
(
5
),
816
830
.
Chow
V. T.
1964
Handbook of Applied Hydrology: A Compendium of Water-Resources Technology
, Vol.
1
.
McGraw-Hill, New York
.
Echeverribar
I.
Morales-Hernández
M.
Brufau
P.
García-Navarro
P.
2019
2D numerical simulation of unsteady flows for large scale floods prediction in real time
.
Advances in Water Resources
134
,
103444
.
Edinger
J. E.
Duttweiler
D. W.
Geyer
J. C.
1968
The response of water temperatures to meteorological conditions
.
Water Resources Research
4
(
5
),
1137
1143
.
García-Navarro
P.
Murillo
J.
Fernández-Pato
J.
Echeverribar
I.
Morales-Hernández
M.
2019
The shallow water equations and their application to realistic cases
.
Environmental Fluid Mechanics
19
(
5
),
1235
1252
.
Gordillo
G.
Morales-Hernández
M.
García-Navarro
P.
2019
Finite volume model for the simulation of 1D unsteady river flow and water quality based on the WASP
.
Journal of Hydroinformatics
22
(
2
),
327
345
.
Gordillo
G.
Morales-Hernández
M.
García-Navarro
P.
2020
A gradient-descent adjoint method for the reconstruction of boundary conditions in a river flow nitrification model
.
Environmental Science: Processes & Impacts
22
(
2
),
381
397
.
Hosseinipour
E.
Martin
J.
1992
Rivmod, A one-Dimensional Hydrodynamic and Sediment Transport Model, Model Theory and Users Manual
.
AScI Corporation at USEPA
,
Athens
, p.
55
.
Lacasta
A.
Morales-Hernández
M.
Murillo
J.
García-Navarro
P.
2014
An optimized GPU implementation of a 2D free surface simulation model on unstructured meshes
.
Advances in Engineering Software
78
,
1
15
.
Liu
X.
2019
A robust numerical model for shallow water governing solute transport with wet/dry interfaces
.
Computer Methods in Applied Mechanics and Engineering
351
,
85
108
.
Martin
J. L.
Wool
T.
Olson
R.
2002
A Dynamic One-Dimensional Model of Hydrodynamics and Water Quality. EPDRiv1 Version 1
.
Morales-Hernández
M.
Murillo
J.
García-Navarro
P.
2019
Diffusion–dispersion numerical discretization for solute transport in 2D transient shallow flows
.
Environmental Fluid Mechanics
19
(
5
),
1217
1234
.
Murillo
J.
García-Navarro
P.
Burguete
J.
2009
Conservative numerical simulation of multi-component transport in two-dimensional unsteady shallow water flow
.
Journal of Computational Physics
228
(
15
),
5539
5573
.
Murillo
J.
Latorre
B.
García-Navarro
P.
2012
A Riemann solver for unsteady computation of 2d shallow flows with variable density
.
Journal of Computational Physics
231
(
14
),
4775
4807
.
Park
R.
Clough
J.
2014
Aquatox (Release 3.1 Plus), Modeling Environmental Fate and Ecological Effects in Aquatic Ecosystems
, Vol.
2
.
Technical documentation. US Environmental Protection Agency (USEPA)
,
Washington, DC
, p.
354
.
Smith
L. S.
Liang
Q.
2013
Towards a generalised GPU/CPU shallow-flow modelling tool
.
Computers & Fluids
88
,
334
343
.
Thomann
R. V.
Mueller
J. A.
1987
Principles of Surface Water Quality Modeling and Control
.
Harper Row, New York
.
Vacondio
R.
Dal Palù
A.
Mignosa
P.
2014
GPU-enhanced finite volume shallow water solver for fast flood simulations
.
Environmental Modelling & Software
57
,
60
75
.
Warn
A.
1987
Simcat-a Catchment Simulation Model for Planning Investment for River Quality
.
Systems Analysis in Water Quality Management Pergamon Press
,
New York
, pp.
211
218
,
10 ref
.
Warren
I.
Bach
H.
1992
Mike 21: a modelling system for estuaries, coastal waters and seas
.
Environmental Software
7
(
4
),
229
240
.
Wells
S. A.
2019
CE-QUAL-W2: a two-dimensional, laterally averaged, hydrodynamic and water quality model, version 4.2. Available from: http://www.ce.pdx.edu/w2/
.
Whitehead
P.
Williams
R.
Lewis
D.
1997
Quality simulation along river systems (quasar): model theory and development
.
Science of the Total Environment
194
,
447
456
.
Wool
T. A.
Ambrose
R. B.
Martin
J. L.
Comer
E. A.
Tech
T.
2006
Water quality analysis simulation program (WASP). Users Manual, Version 6
.