A GPU-based 2D shallow water quality model

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

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

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.

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.

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

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.

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:

  

  (9)

  with

  

  (10)

  and the numerical solute concentration at edge k separating cells i and ⁠:

  

  (11)
• Step (III). Compute the updating value of the reaction R and source f terms as:

  

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

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

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.

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.

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.

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.

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.

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.

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

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

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

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

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

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

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

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

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.

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

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

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 [email protected] 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.

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.

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.

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

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