Pollutant transport by shallow water equations on unstructured meshes: Hyperbolization of the model and numerical solution via a novel flux splitting scheme

Abstract

The purpose of this paper is twofold. First, using the Cattaneo's relaxation approach, we reformulate the system of governing equations for the pollutant transport by shallow water flows over non-flat topography and anisotropic diffusion as hyperbolic balance laws with stiff source terms. The proposed relaxation system circumvents the infinite wave speed paradox which is inherent in standard advection–diffusion models. This turns out to give a larger stability range for the choice of the time step. Second, following a flux splitting approach, we derive a novel numerical method to discretise the resulting problem. In particular, we propose a new flux splitting and study the associated two systems of differential equations, called the “hydrodynamic” and the “relaxed diffusive” system, respectively. For the presented splitting we analyse the resulting two systems of differential equations and propose two discretisation schemes of the Godunov-type. These schemes are simple to implement, robust, accurate and fast when compared with existing methods. The resulting method is implemented on unstructured meshes and is systematically assessed for accuracy, robustness and efficiency on a carefully selected suite of test problems including non-flat topography and wetting and drying problems. Formal second order accuracy is assessed through convergence rates studies.

Introduction

Numerical simulations of the advection–diffusion of a passive tracer (e.g. nutrients, pollutants, temperature, etc.) in natural water bodies are nowadays fundamental tools for assessing water pollution, designing measures to improve water quality and evaluating possible effects on the ecosystems (e.g. [1], [2], [3]). The shallow water equations are the usual model adopted for describing flows in environments characterised by horizontal dimensions much larger than the vertical extent. If the passive tracer is vertically well-mixed, their dynamics can be represented by a depth-averaged advection–diffusion equation. The coupling of the shallow water equations with the advection–diffusion equation results in a system of partial differential equations (PDE) of parabolic type. They govern different physical processes having different scales: in the advective process small waves propagate at a finite velocity while in the diffusive process small perturbations of either the flow field or concentration of the passive tracer propagate with infinite velocity. This aspect together with the presence of domains characterised by complex morphologies in natural water bodies, where also wetting and drying regions may appear, makes this problem particularly challenging from the numerical point of view. Over the past decade, numerical schemes for scalar transport problems have attracted much attention with particular focus on the case of scalar advection [4], [5], [6], the more complex problem taking into account both advection and diffusion [7], [8], [9] and the case of advection—diffusion of a reactive tracer [10], [11]. Many other related studies can be found in the literature.

The aim of this paper is to present a novel, robust and accurate numerical approach for solving the shallow water advection diffusion equations over complex topographies. First, we reformulate the original advection–diffusion problem via a relaxation approach proposed by Cattaneo [12], [13] and Vernotte [14], transforming the initial parabolic system of PDE in a hyperbolic system containing stiff source terms [15]. This approach consists of augmenting the Fourier law with a transient term involving a parameter named “relaxation time”: when such parameter is small, the original Fourier law is recovered. This relaxation technique is appropriate for the construction of numerical schemes that transform advection–diffusion equations into hyperbolic systems. This technique i) allows a large enough relaxation parameter suitable for numerical implementations; ii) allows relaxation of the spatial gradients of the passive tracer concentration removing second order derivatives present in the original model with only minor changes in the mathematical formulation of the problem and maintaining the physical meaning of the relaxed model; iii) allows to compute accurate and efficient solutions by implementing numerical schemes able to handle stiff source terms.

The first numerical applications of the relaxation in the sense of Cattaneo were implemented by Nishikawa [16], who computed steady-state solutions of model advection diffusion equations. The Cattaneo's relaxation approach was put on firm grounds in terms of stiff hyperbolic conservation laws for time dependent problems by Toro and Montecinos [15] and it has been applied to solve different kinds of advection–diffusion time marching problems. For example, Gómez et al. [17] proposed a finite element solution of the convection-diffusion problems and Montecinos et al. [18] relaxed the one-dimensional problem for the solution of blood flow in viscoelastic collapsible tubes. Here, using the Cattaneo's relaxation approach (sensu Toro and Montecinos [15]) we reformulate the system of PDE governing the advection–diffusion processes of a passive scalar transported by shallow flows. The resulting hyperbolic relaxation system contains stiff source terms and is solved numerically using a Godunov-type approach proposing a new flux splitting method (for background and reviews see for example [19], [20]). A distinct feature of the new scheme is that it separates the hydrodynamics from the relaxed diffusive problem. This provides the possibility of taking advantage of the broad availability of Riemann solvers (exact or approximate) for the numerical solution of the hydrodynamic part [e.g. 21]. It turns out that the relaxed diffusive system is hyperbolic and a simple solution of the associated Riemann problem proves to be simple and robust and the flux evaluation is immediate.

The stiff source terms are handled by means of a locally implicit ADER scheme proposed by Montecinos and Toro [22] and topographical source terms are taking into account adopting the robust approach proposed by Duran et al. [23]. Friction terms are included using a standard splitting technique.

Finally, the proposed model is applied to well-established test problems with the aim of demonstrating that it is able i) to correctly treat the propagation of wetting and drying fronts preserving mass which is challenging when treating unsteady flows over complex topographies; ii) to preserve particular types of steady states, thus correctly balancing fluxes and source terms (well-balanced property); this contributes to improve the solution of steady and transient flows at any flow condition [24]; iii) to achieve second order of accuracy in space and time.

The rest of this paper is structured as follows: Section 2 briefly reviews the governing equations and in Section 3 we reformulate the full problem analysing the properties of the newly developed relaxation system. In section 4 we present the new developed flux splitting scheme. In Section 5 and 6 we describe how the source terms are discretised and how second order is achieved, while in Section 7 we discuss about the best choice for the relaxation parameter. Finally in Section 8 we present numerical results for a range of carefully selected test problems to assess both the robustness and accuracy of the schemes proposed in this paper. Conclusions are drawn in Section 9.