Combined Use of High-Resolution Numerical Schemes to Reduce Numerical Diffusion in Coupled Nonhydrostatic Hydrodynamic and Solute Transport Model

Abstract

In three-dimensional simulations of free-surface flow where the vertical velocities are relevant, such as in lakes, estuaries, reservoirs, and coastal zones, a nonhydrostatic hydrodynamic approach may be necessary. Although the nonhydrostatic hydrodynamic approach improves the physical representation of pressure, acceleration and velocity fields, it is not free of numerical diffusion. This numerical issue stems from the numerical solution employed in the advection and diffusion terms of the Reynolds-averaged Navier–Stokes (RANS) and solute transport equations. The combined use of high-resolution schemes in coupled nonhydrostatic hydrodynamic and solute transport models is a promising alternative to minimize these numerical issues and determine the relationship between numerical diffusion in the two solutions. We evaluated the numerical diffusion in three numerical experiments, for different purposes: The first two experiments evaluated the potential for reducing numerical diffusion in a nonhydrostatic hydrodynamic solution, by applying a quadratic interpolator over a Bilinear, applied in the Eulerian–Lagrangian method (ELM) step-ii interpolation, and the capability of representing the propagation of complex waves. The third experiment evaluated the effect on numerical diffusion of using flux-limiter schemes over a first-order Upwind in solute transport solution, combined with the interpolation methods applied in a coupled hydrodynamic and solute transport model. The high-resolution methods were able to substantially reduce the numerical diffusion in a solute transport problem. This exercise showed that the numerical diffusion of a nonhydrostatic hydrodynamic solution has a major influence on the ability of the model to simulate stratified internal waves, indicating that high-resolution methods must be implemented in the numerical solution to properly simulate real situations.

1. Introduction

When the ratio of vertical to horizontal motion scales is not small (e.g., flows over abruptly changing bottom topography, orbital movements in short-wave motions, or intensive vertical circulation), a nonhydrostatic approach may be necessary to accurately simulate three-dimensional free-surface flows in large aquatic ecosystems such as lakes, estuaries, reservoirs, and coastal zones [1,2,3]. There is a trade-off between method efficiency and the computational costs when hydrostatic and nonhydrostatic approaches are compared. The nonhydrostatic approach, in general, requires more complex interactive numeric methods [4,5], which increases the computation cost compared to the hydrostatic approach using the same computation grid. The nonhydrostatic approach also improves the physical representation of the phenomenon and requires fewer vertical layers to obtain physically satisfactory results than does the hydrostatic approach. Therefore, to obtain similar results, the hydrostatic approach requires more vertical layers, necessitating a higher computation cost than the nonhydrostatic approach [6].

Whatever the approach (hydrostatic or nonhydrostatic), the numerical diffusion issue remains critical, and many studies have attempted to reduce it by adopting diverse approaches. Higher-resolution techniques (e.g., quadratic, cubic and spline interpolators, and flux limiters) have proven to be a promising alternative in order to minimize these numerical issues [7,8,9,10,11]. Despite previous studies showing that higher-resolution schemes might not maintain some important characteristics of a numerical solution (e.g., spurious oscillations of internal waves) [12], they may better represent the nonlinear behavior of a velocity field than low-resolution schemes (upwind solute transport solver, bilinear interpolation). In general, linear schemes might be more stable, monotonic, and easily implemented, but they still might give unsatisfactory results due to higher numerical diffusion. Generally, higher-resolution techniques are applied to evaluate the reduction of numerical diffusion separately in RANS-based (Reynolds-averaged Navier–Stokes) or solute transport models, without determining the interrelationship with the numerical diffusion caused by both solutions.

The Eulerian–Lagrangian Method (ELM) is one of the most popular and accurate techniques used to numerically solve the advection terms in RANS and solute transport equations [7,8,13,14,15,16,17,18,19,20]. The purpose of ELM is to combine the simplicity of the fixed Eulerian computational grid with a stable and accurate Lagrangian approach. In summary, the ELM has two major steps [21]: (1) a particle-tracking method, which consists of obtaining the location of the departure point of a fluid particle in the Lagrangian step by integrating a characteristic equation backward in time using a certain number of sub-time steps, and (2) a repeated interpolation of the advected field in each sub-time step required to estimate the velocity components at the departure point. As the result of the interpolation errors in the two steps, the ELM may introduce a substantial numerical diffusion when low-resolution interpolation schemes are used, or spurious oscillation when high-resolution interpolation schemes are used in the hydrodynamic and mass transport solutions [7,8,22].

The interpolation technique used in the second step has a major impact on the total numerical diffusion produced in ELM, especially when low-resolution interpolation schemes are used, as both particle tracking and velocity retrieval are substantially affected by the interpolation process [21]. Accordingly, Hodges et al. [16] proposed a stable non-conciliatory quadratic Lagrange interpolation for a three-dimensional mesh, where 27 grid points are used to estimate the velocity values during the particle-tracking process. This technique has not yet been formally analyzed (Hodges et al. [23], without any posterior record of formal analyses), although previous studies successfully modeled hydrodynamic and solute transport simulations, strongly indicating that the solution satisfactorily represented internal gravity waves and may improve the ability of ELM to solve the free surface motion [16,23,24,25,26,27,28].

The solute transport solution can also be obtained using other numerical schemes besides ELM, such as high-order Total Variation Diminishing (TVD) schemes, Random Walk Particle Tracking (RWPT), and Smoothed Particle Hydrodynamics (SPH), which have been thoroughly discussed by Boso et al. [29]. One simple scheme to discretize the mass conservation equation, with a low cost of implementation, is the first-order Upwind, to which the flux-limiter technique can be applied in order to reduce the well-known numerical diffusion of this discretization [30,31]. The flux-limiter scheme proved to be one of the most effective approaches to constructing a nonlinear high-resolution scheme [9]. These schemes are simple functions that define a convection scheme based on a ratio of local gradients in the solution field, with a limited behavior in order to maintain the monotonicity, usually defined by a Total Variation Diminishing (TVD) condition [9,32]. This method has been widely used to retrieve less diffusive and stable solutions, even in critical situations, such as in sharp concentration gradients in regions with a small Courant number [10,11,17,19,30,31,33,34].

For solute transport models, it is usually assumed that the numerical solution is free of hydrodynamic numerical diffusion and that any error is leads to the mass transport solution, and therefore the formal evaluation of the effect of hydrodynamic numerical diffusion on the solute transport solution may not have been properly addressed [7,8,13]. Therefore, coupling of hydrodynamic and solute transport models may be useful in order to determine the interrelationship between numerical diffusion in the two solutions, which is still underexplored. Usually, studies using a coupled model propose improvements (high-resolution methods or higher order equations) of the hydrodynamic or transport module to investigate numerical diffusion in one of the solutions, separately; Without accounting for the possibility that the use, or not, of high-resolution methods in one solution may affect the other [10,11,17,19,30,31,35]. Thus, further effort is needed to understand how numerical diffusion, produced by low-resolution schemes in the hydrodynamic solution, may be transferred to solute transport and the interrelationship between solutions.

To explore this question, we applied high-resolution methods to hydrodynamic and solute transport solutions using a coupled model. We also investigated the individual performance of each method in reducing numerical diffusion in one of the solutions, and the performance of the combined use of the methods. Thus, demonstrating the impact that numerical diffusion in one solution can have on another solution. Specifically, we tested a nonlinear quadratic interpolation in contrast to a linear one in ELM step-ii, combined with several flux-limiter schemes applied in the solute transport module, to reduce numerical diffusion, and evaluated the effect of the combined use of these methods on the solute transport numerical diffusion. We used a coupled hydrodynamic solute transport model with nonhydrostatic and flux-limiter approaches to evaluate the numerical diffusion in three numerical experiments, for different purposes: (a) the first two experiments evaluated the potential for reducing the numerical diffusion in the hydrodynamic solution, by applying a nonlinear interpolator rather than a linear one, and the capability of representing the propagation of complex waves; (b) the third experiment evaluated the effect of using, or not, high-resolution schemes on the numerical diffusion of combined hydrodynamic and solute transport models.

2. Mathematical Considerations

2.1. Governing Equations

The RANS equations are used to describe three-dimensional free-surface flows. These equations express the physical principle of volume, mass, and momentum conservation. The momentum equations for an incompressible fluid have the following form,

$$\begin{array}{ccc}\hfill \frac{\partial \mathit{u}}{\partial t}+\mathit{u}\frac{\partial \mathit{u}}{\partial x}+\mathit{v}\frac{\partial \mathit{u}}{\partial y}+\mathit{w}\frac{\partial \mathit{u}}{\partial z}-f\mathit{v}& =& -\frac{\partial {\mathit{p}}_a}{\partial x}-g\frac{\partial \eta }{\partial x}-g\frac{\partial }{\partial x}\left[{\int }_z^{\eta }\frac{\rho -{\rho }_0}{{\rho }_0}d\zeta \right]-\frac{\partial \mathit{q}}{\partial x}\hfill \\ \hfill +{\nu }^h\left(\frac{{\partial }^2\mathit{u}}{\partial x^2}+\frac{{\partial }^2\mathit{u}}{\partial y^2}\right)+\frac{\partial }{\partial z}\left({\nu }^v\frac{\partial \mathit{u}}{\partial z}\right)\end{array}$$

(1)

$$\begin{array}{ccc}\hfill \frac{\partial \mathit{v}}{\partial t}+\mathit{u}\frac{\partial \mathit{v}}{\partial x}+\mathit{v}\frac{\partial \mathit{v}}{\partial y}+\mathit{w}\frac{\partial \mathit{v}}{\partial z}+f\mathit{u}& =& -\frac{\partial {\mathit{p}}_a}{\partial y}-g\frac{\partial \eta }{\partial y}-g\frac{\partial }{\partial y}\left[{\int }_z^{\eta }\frac{\rho -{\rho }_0}{{\rho }_0}d\zeta \right]-\frac{\partial \mathit{q}}{\partial y}\hfill \\ \hfill +{\nu }^h\left(\frac{{\partial }^2\mathit{v}}{\partial x^2}+\frac{{\partial }^2\mathit{v}}{\partial y^2}\right)+\frac{\partial }{\partial z}\left({\nu }^v\frac{\partial \mathit{v}}{\partial z}\right)\end{array}$$

(2)

$$\begin{array}{c}\hfill \frac{\partial \mathit{w}}{\partial t}+\mathit{u}\frac{\partial \mathit{w}}{\partial x}+\mathit{v}\frac{\partial \mathit{w}}{\partial y}+\mathit{w}\frac{\partial \mathit{w}}{\partial z}=-\frac{\partial \mathit{q}}{\partial z}+{\nu }^h\left(\frac{{\partial }^2\mathit{w}}{\partial x^2}+\frac{{\partial }^2\mathit{w}}{\partial y^2}\right)+\frac{\partial }{\partial z}\left({\nu }^v\frac{\partial \mathit{w}}{\partial z}\right)\end{array}$$

(3)

where $\mathit{u}(x,y,z,t)$, $\mathit{v}(x,y,z,t)$, and $\mathit{w}(x,y,z,t)$ are the velocity components in the horizontal (x and y) and vertical (z) directions, respectively; ${\nu }^h$ e ${\nu }^v$ are the horizontal and vertical turbulent eddy viscosity coefficients, respectively; t is the time; $p_a(x,y,z,t)$ is the atmospheric pressure; and $\eta$ is the free-surface elevation from a water-level reference. The second and third terms on the right-hand side of Equations (1) and (2) represent the barotropic and the baroclinic contributions to the hydrostatic pressure; $q(x,y,z,t)$ denotes the nonhydrostatic pressure component; f is the Coriolis parameter; and g is the gravitational acceleration.

When a simple hydrostatic approach is considered, Equation (3) is neglected and q is assumed to be equal to zero in Equations (1) and (2). In this case, it is assumed that the vertical acceleration does not have a significant effect on the velocity field in comparison with the horizontal acceleration, which is the assumption usually applied for simulation of shallow waters (see, e.g., in [36,37,38,39]).

The volume conservation is expressed by the incompressibility condition and the continuity equation, given by

$$\frac{\partial \mathit{u}}{\partial x}+\frac{\partial \mathit{v}}{\partial y}+\frac{\partial \mathit{w}}{\partial z}=0$$

(4)

Integrating Equation (4) over depth leads to the following equation,

$${\int }_{-h}^{\eta }\left[\frac{\partial \mathit{u}}{\partial x}+\frac{\partial \mathit{v}}{\partial y}+\frac{\partial \mathit{w}}{\partial z}\right]dz={\int }_{-h}^{\eta }\frac{\partial \mathit{u}}{\partial x}dz+{\int }_{-h}^{\eta }\frac{\partial \mathit{v}}{\partial y}dz+{\int }_{-h}^{\eta }\frac{\partial \mathit{w}}{\partial z}dz=0$$

(5)

where h is the bathymetry measured from the theoretical undisturbed water surface. Using the Leibniz integration rule in each direction, in Equation (5), and using a kinematic condition at the free surface, leads to the following free-surface equation,

$$\frac{\partial \eta }{\partial t}+\frac{\partial }{\partial x}{\int }_{-h}^{\eta }\mathit{u}dz+\frac{\partial }{\partial y}{\int }_{-h}^{\eta }\mathit{v}dz=0$$

(6)

Finally, the mass conservation of a conservative scalar variable may be expressed by the following equation,

$$\frac{\partial C}{\partial t}+\frac{\partial \left(\mathit{u}C\right)}{\partial x}+\frac{\partial \left(\mathit{v}C\right)}{\partial y}+\frac{\partial \left(\mathit{w}C\right)}{\partial z}=\frac{\partial }{\partial x}\left(K^h\frac{\partial C}{\partial x}\right)+\frac{\partial }{\partial y}\left(K^h\frac{\partial C}{\partial y}\right)+\frac{\partial }{\partial z}\left(K^v\frac{\partial C}{\partial z}\right)$$

(7)

where C is the concentration of a conservative substance being transported (e.g., salinity), and $K^h$ and $K^v$ are the horizontal and vertical turbulent eddy diffusivity coefficients, respectively.

For both the velocity field and scalar transport solutions, the boundary conditions are implemented under the assumption of “free-slip” boundaries. The Dirichlet and Neumann conditions were assigned to represent the normal and tangential velocities in the solid boundaries, respectively. For the scalar solution, a no-flux boundary condition is assumed in solid boundaries.

The tangential stress boundary conditions for the momentum equations at the free-surface are specified by the prescribed wind stresses, which can be approximated as

$${\nu }^v\frac{\partial \mathit{u}}{\partial z}={\gamma }_T({\mathit{u}}_a-\mathit{u}),{\nu }^v\frac{\partial \mathit{v}}{\partial z}={\gamma }_T({\mathit{v}}_a-\mathit{v});atz=\eta$$

(8)

where ${\mathit{u}}_a$ and ${\mathit{v}}_a$ are the horizontal wind velocity components and ${\gamma }_T$ is a non-negative wind stress coefficient. The bottom friction is specified by

$${\nu }^v\frac{\partial \mathit{u}}{\partial z}={\gamma }_{B}u,{\nu }^v\frac{\partial \mathit{v}}{\partial z}={\gamma }_{B}v;atz=-h$$

(9)

where ${\gamma }_B$ is a non-negative bottom friction coefficient, which is typically represented by means of a Manning or Chezy coefficient.

2.2. Grid and Variable Locations

The computation grid can be described as a generic unstructured orthogonal grid, having $N_p$ elements, each having an arbitrary number of sides $S_i\ge 3,i=1,2,\dots ,N_p$ (Figure 1).

Let $N_s$ be the total number of sides in the grid. The length of each side is ${\lambda }_j,j=1,2,\dots ,N_p$. The vertical faces of the i-th element are identified by an index $j_{(i,l)},l=1,2,\dots ,S_i$, so that $1\le j_{(i,l)}\le N_s$. Similarly, the two polygons that share the j-th vertical face of the grid are identified by the indices $i_{(j,1)}$ and $i_{(j,2)}$, so that $1\le i_{(j,1)}\le N_p$ and $1\le i_{(j,2)}\le N_p$. The non-zero distance between centers of two adjacent polygons that share the j-th side is denoted by ${\delta }_j$.

Along the vertical direction, a simple finite difference discretization, not necessarily uniform, is adopted. By denoting a given level surface as $\mathsf{\Delta }{z}_{k+\frac{1}{2}}$ the vertical discretization step is defined by

$$\mathsf{\Delta }{z}_k=\mathsf{\Delta }{z}_{k+\frac{1}{2}}-\mathsf{\Delta }{z}_{k-\frac{1}{2}}k=1,2,\cdots ,N_s$$

(10)

The three-dimensional space discretization consists of elements whose horizontal faces are the polygons of a given orthogonal grid, represented by the layers at $k+\frac{1}{2}$ (upper face) or $k-\frac{1}{2}$ (bottom face), and whose height, for each layer, is $\mathsf{\Delta }{z}_k$. The water surface elevation ($\eta$) is located at the barycenter of the upper horizontal face for each i-th element. The velocity component normal to each horizontal face is assumed to be constant over the face of each computation cell, which is defined at the point of intersection between the face and the segment joining the centers of the two prisms that share the face. The nonhydrostatic pressure component $q_{i,k}^n$ and the concentrations $C_{i,k}^n$ are located at the center of the i-th computation cell, halfway between $\mathsf{\Delta }{z}_{k+\frac{1}{2}}$ and $\mathsf{\Delta }{z}_{k-\frac{1}{2}}$. Finally, the water depth $h_j$ is specified and assumed constant on each vertical face of an element.

2.3. Numerical Approach

2.3.1. Rans Equations

We used a semi-implicit method ($\theta$-method in [40]) in a finite-difference finite-volume model, with a Eulerian–Lagrangian Method [14] to solve the convective and viscous terms of the RANS equations, and a fractional-step framework [5] to solve the pressure component by splitting the pressure into hydrostatic and nonhydrostatic parts. The first step of the hydrodynamic solution is to compute the provisional water velocity and surface elevation, neglecting the implicit contribution of the nonhydrostatic pressure. In the second fractional step, the provisional velocity ($\tilde{\mathit{u}}$, $\tilde{\mathit{v}}$, and $\tilde{\mathit{w}}$) and provisional surface elevation ($\tilde{\eta }$) are corrected by nonhydrostatic pressure terms. A complete description of the numerical solution for nonhydrostatic flows was provided by Casulli and Lang [5], and is here named CL04. Here, we describe in detail only the Eulerian–Lagrangian method to clarify where the high-resolution methods were applied.

One difficulty in the numerical treatment of RANS equations arises from the discretization of the convective and viscous terms of the convection-diffusion equation in three dimensions (Equations (1)–(3)). The advection–diffusion equation without Coriolis and pressure terms can be described as

$$\frac{dc}{dt}=\frac{\partial c}{\partial t}+\mathit{u}\frac{\partial c}{\partial x}+\mathit{v}\frac{\partial c}{\partial y}+\mathit{w}\frac{\partial c}{\partial z}=\mu \left(\frac{{\partial }^{2}c}{\partial x^2}+\frac{{\partial }^{2}c}{\partial y^2}\right)+\frac{\partial }{\partial z}\left(\nu \frac{\partial c}{\partial z}\right)$$

(11)

where c is a generic variable (e.g., velocity components $\mathit{u}$, $\mathit{v}$, and $\mathit{w}$) [14]. To simply, solve, and improve the stability and accuracy of an explicit finite-difference method, consider Equation (11) in the Lagrangian form:

$$\frac{dc}{dt}=\mu \left(\frac{{\partial }^{2}c}{\partial x^2}+\frac{{\partial }^{2}c}{\partial y^2}\right)+\frac{\partial }{\partial z}\left(\nu \frac{\partial c}{\partial z}\right)$$

(12)

where the substantial derivative $\frac{d}{dt}$ indicates that the rate of change over time, which is calculated along the stream line, is defined by

$$\frac{dx}{dt}=\mathit{u},\frac{dy}{dt}=\mathit{v},\frac{dz}{dt}=\mathit{w}$$

(13)

An explicit discretization of Equation (13) can be given by

$$\frac{(C_{j,k}^{n+1}-C_{j,k}^{*})}{\mathsf{\Delta }t}=\frac{{\nu }_{k+\frac{1}{2}}\frac{C_{j,k+1}^{n+1}-C_{j,k}^{n+1}}{\mathsf{\Delta }{z}_{j,k+\frac{1}{2}}}-{\nu }_{k-\frac{1}{2}}\frac{C_{j,k}^{n+1}-C_{j,k-1}^{n+1}}{\mathsf{\Delta }{z}_{j,k-\frac{1}{2}}}}{\mathsf{\Delta }{z}_{j,k}}+\mu {\mathsf{\Delta }}_h{C}_{j,k}^{*}$$

(14)

where $C_{i,j,k}^{*}$ denotes a generic variable at the j-th side of a grid at vertical level and time step n.

So as to apply the Lagrangian discretization in a Eulerian grid, a backward stream line is required to estimate the departure velocity at time “$\mathit{n}$” (Lagrangian approach) required to reach the final position “$\mathit{j},\mathit{k}$” at time “$\mathit{n}+1$” (Eulerian approach). To estimate the particle departure point, the ELM step-i interpolation is applied, which requires a particle-tracking method. A multistep backward Euler (MSE) was applied, but other methods can be satisfactorily applied as well (e.g., fifth-order Runge–Kutta) [22]. However, the departure position is not a grid point, and an interpolation formula using known node points must be used to define $C_{j,k}^{*}$ (ELM step-ii).

A consistent semi-implicit finite difference discretization is used to calculate the provisional horizontal velocity component from momentum Equations (1) and (2) and takes the following form (CL04),

$$\begin{array}{c}{\tilde{\mathit{u}}}_{j,k}^{n+1}=\mathbf{F}{\mathit{u}}_{j,k}^n-(1-\theta )\frac{\mathsf{\Delta }t}{{\delta }_j}\left[g({\eta }_{i(j,r)}^n-{\eta }_{i(j,l)}^n)+(q_{i(j,r),k}^n-q_{i(j,l),k}^n)\right]-\hfill \\ \theta g\frac{\mathsf{\Delta }t}{{\delta }_j}({\eta }_{i(j,r)}^{n+1}-{\eta }_{i(j,l)}^{n+1})+\frac{\mathsf{\Delta }t}{\mathsf{\Delta }{Z}_{j,k}^n}\left[{\nu }_{j,k+\frac{1}{2}}^v\frac{{\tilde{\mathit{u}}}_{j,k+1}^{n+1}-{\tilde{\mathit{u}}}_{j,k}^{n+1}}{\mathsf{\Delta }{Z}_{j,k+\frac{1}{2}}^n}-{\nu }_{j,k-\frac{1}{2}}^v\frac{{\tilde{\mathit{u}}}_{j,k}^{n+1}-{\tilde{\mathit{u}}}_{j,k-1}^{n+1}}{\mathsf{\Delta }{Z}_{j,k-\frac{1}{2}}^n}\right]\hfill \\ k=m_j,m_{j+1},\cdots ,M_j^n\hfill \end{array}$$

(15)

where ${\mathit{u}}_{j,k}^n$ denotes the horizontal velocity component normal to the j-th side of the grid at vertical level k, and time step n; ${\eta }_{i(j,r)}^n$ is the free-surface level at the right neighbor i-th element from the j-th side and ${\eta }_{i(j,l)}^n$ from left neighbor; $q_{i(j,r),k}^n$ denotes the nonhydrostatic pressure component; and $\mathbf{F}$ is an explicit finite-difference operator, which accounts for the contributions from discretization of the air pressure, Coriolis, baroclinic pressure, advection, and horizontal friction terms.

Analogously to Equation (15), a semi-implicit finite-difference discretization for the provisional vertical velocity component (neglecting the nonhydrostatic pressure contribution) at the top of each computation cell is derived from Equation (3):

$$\begin{array}{c}\hfill {\tilde{\mathit{w}}}_{i,k}^{n+1}=\mathbf{F}{\mathit{w}}_{i,k+\frac{1}{2}}^n-(1-\theta )\frac{\mathsf{\Delta }t}{\mathsf{\Delta }{Z}_{i,k+\frac{1}{2}}^n}\left(q_{i,k+1}^n-q_{i,k}^n\right)+\frac{\mathsf{\Delta }t}{\mathsf{\Delta }{Z}_{i,k+\frac{1}{2}}^n}\\ \hfill \left[{\nu }_{i,k+1}^v\frac{{\tilde{\mathit{w}}}_{i,k+\frac{3}{2}}^{n+1}-{\tilde{\mathit{w}}}_{i,k+\frac{1}{2}}^{n+1}}{\mathsf{\Delta }{Z}_{i,k+1}^n}-{\nu }_{i,k}^v\frac{{\tilde{\mathit{w}}}_{i,k+\frac{1}{2}}^{n+1}-{\tilde{\mathit{w}}}_{i,k-\frac{1}{2}}^{n+1}}{\mathsf{\Delta }{Z}_{i,k}^n}\right]k=m_i,m_{i+1},\cdots ,M_i^n-1\end{array}$$

(16)

where $Fw$ accounts only for contributions from the discretization of the advection and horizontal friction terms. The Eulerian–Lagrangian method can be applied to discretize the finite difference operator $\mathbf{F}$ in Equations (15) and (16), given by

$$\mathbf{F}{\mathit{u}}_{j,k}^n={\mathit{u}}_{j,k}^{*}+{\nu }^h\mathsf{\Delta }t{\mathsf{\Delta }}_h{\mathit{u}}_{j,k}^{*}+f\mathsf{\Delta }t{\mathit{v}}_{j,k}^{*}-g\frac{\mathsf{\Delta }t}{{\delta }_j{\rho }_0}\sum _{\ell =k}^{M_j}{\omega }_{\ell }\left[{\rho }_{i(j,2),\ell }^n-{\rho }_{i(j,1),\ell }^n\right]\mathsf{\Delta }{z}_{j,\ell }^n$$

(17)

$$\mathbf{F}{\mathit{w}}_{i,k+\frac{1}{2}}^n={\mathit{w}}_{i,k+\frac{1}{2}}^{*}+{\nu }^h\mathsf{\Delta }t{\mathsf{\Delta }}_h{\mathit{w}}_{i,k+\frac{1}{2}}^{*}$$

(18)

where ${\mathit{u}}_{j,k}^{*}$ denotes the horizontal velocity component normal to the j-th side of the grid interpolated at time $t^n$ at the end of the Lagrangian trajectory; ${\mathit{v}}_{j,k}^{*}$ denotes the horizontal velocity component orthogonal to $u_{j,k}^{*}$; ${\mathit{w}}_{i,k+\frac{1}{2}}^{*}$ denotes the vertical velocity component normal to the upper horizontal face of element “i” at layer “$k+\frac{1}{2}$”. ${\omega }_k=\frac{1}{2}$ and ${\omega }_{\ell }=1$, for $\ell \ne k$; $\mathsf{\Delta }h$ is the discretization of the horizontal Laplacian. The Lagrangian trajectory is calculated by integrating the velocity backward in time from a face’s $(\mathit{j},\mathit{k})$ barycenter at $t^{n+1}$ to its location at time $t^n$. The second step of the interpolation technique is discussed further below.

2.3.2. Solute Transport Equation

We used a conservative finite-volume scheme, with a semi-implicit approach based on CL04 with time-accurate local time stepping based on that in [41], to discretize the solute transport equation. The CL04 approach yields a conservative solution, also respecting the max-min property.

Assume that $S_k^{+}$ represents the set of vertical faces belonging to the computation cell $(\mathit{i},\mathit{k})$, through which water is leaving the respective computation cell; $S_k^{-}$ represents the set of vertical faces, through which water is entering the same computation cell; and $p(\mathit{i},\mathit{j})$ is the neighbor of the computation cell $(\mathit{i},\mathit{k})$ that shares the vertical face ($\mathit{i},\mathit{j}$). For each computational cell, Equation (7) for scalar transport is discretized as follows,

$$\begin{array}{c}{P}_i\mathsf{\Delta }{z}_{i,k}^{n+1}{C}_{i,k}^{n+1}=P_i\mathsf{\Delta }{z}_{i,k}^n{C}_{i,k}^n-\mathsf{\Delta }t\left[\sum _{j\in S_{i,k}^{+}}|Q_{j,k}^{n+\theta }|C_{i,k}^n-\sum _{j\in S_{i,k}^{-}}\left|Q_{j,k}^{n+\theta }\right|C_{r(i,j),k}^n\right]\hfill \\ -\mathsf{\Delta }t\left[|Q_{i,k+\frac{1}{2}}^{n+\theta }|C_{i,k}^n-|Q_{i,k-\frac{1}{2}}^{n+\theta }|C_{i,k-1}^n\right]+\mathsf{\Delta }t\sum _{j\in S_{i,k}^{+}\cup S_{i,k}^{-}}{d}_{j,k}^{n+\theta }[C_{r(i,j),k}^n-C_{i,k}^n]\hfill \\ -\mathsf{\Delta }t\left[d_{i,k+\frac{1}{2}}^{n+\theta }(C_{i,k+1}^{n+1}-C_{i,k}^{n+1})-d_{i,k-\frac{1}{2}}^{n+\theta }(C_{i,k}^{n+1}-C_{i,k-1}^{n+1})\right]\hfill \\ -\frac{\mathsf{\Delta }t}{2}\left[{\psi }_{i,k+\frac{1}{2}}^n|Q_{i,k+\frac{1}{2}}^{n+\theta }|(C_{i,k+1}^n-C_{i,k}^n)-{\psi }_{i,k-\frac{1}{2}}^n|Q_{i,k-\frac{1}{2}}^{n+\theta }|(C_{i,k}^n-C_{i,k-1}^n)\right]\hfill \\ -\frac{\mathsf{\Delta }t}{2}\sum _{j\in S_{i,k}^{+}\cup S_{i,k}^{-}}{\psi }_{j,k}^n\left|Q_{j,k}^{n+\theta }\right|(C_{r(i,j),k}^n-C_{i,k}^n)\hfill \end{array}$$

(19)

where

$$d_{i,k\pm \frac{1}{2}}^{n+\theta }=\max \left(0,D_{i,k\pm \frac{1}{2}}^n-\frac{1}{2}\left|Q_{i,k\pm \frac{1}{2}}^{n+\theta }\right|\right)$$

(20)

$$d_{j,k}^{n+\theta }=\max \left(0,D_{j,k}^n-\frac{1}{2}\left|Q_{j,k}^{n+\theta }\right|\right)$$

(21)

where $Q_{j,k}^{n+\theta }={\lambda }_j\mathsf{\Delta }{z}_{j,k}^n{\mathit{u}}_{j,k}^{n+\theta }$ and $Q_{i,k\pm \frac{1}{2}}^{n+\theta }=P_i{\mathit{w}}_{j,k\pm \frac{1}{2}}^{n+\theta }$ are the advective flux coefficients, and $D_{j,k}^n={\lambda }_j\mathsf{\Delta }{z}_{j,k}^n\frac{K_{j,k}^h}{{\delta }_j}$ and $D_{i,k\pm \frac{1}{2}}^n=P_i\frac{K_{i,k\pm \frac{1}{2}}^v}{\mathsf{\Delta }{z}_{i,k\pm \frac{1}{2}}^n}$ are the diffusive flux coefficients. In Equation (19) the first term on the right side is the mass in layer k at time n; the second is the horizontal-advection term; the third is the vertical advection, followed by the horizontal and vertical-diffusion terms. The last two terms are the numerical diffusion reduction terms, which depend on the high-resolution scheme to estimate $\psi$.

The $\psi$ represents the flux limiter itself, given by

$${\psi }_{j,k}^n={\mathsf{\Phi }}_{j,k}^n-{\phi }_{j,k}^n$$

(22)

where $\mathsf{\Phi }$ is the partial flux limiter and $\varphi$ is a function that assures the independence between the max-min property and mesh size, given by

$${\phi }_{j,k}^n=min\left(1,\frac{2D_{j,k}^n}{Q_{j,k}^{n+\theta }}\right)$$

(23)

When the horizontal and vertical diffusion is set equal to zero, ${\phi }_{j,k}^n=0$, hence ${\psi }_{j,k}^n={\mathsf{\Phi }}_{j,k}^n$. The upwind scheme may be easily obtained when $\mathsf{\Phi }\left(r\right)=0$.

3. High-Resolution Schemes to Reduce Numerical Diffusion

For the hydrodynamic solution, two different interpolators at ELM step-ii were tested: a simple bilinear interpolation [14] and a quadratic interpolator [16], both applied in a regular structure grid. Regarding the solute transport solution, three different flux limiters were used—MUSCL, Superbee, and Ultimate Quickest—which performed best for coupled hydrodynamic and transport simulations in previous studies [11,16,19,30].

3.1. Flux-Limiter

In order to retrieve some accuracy from the solute transport first order upwind scheme, an additional “antidiffusive” term is used (Equation (22)), the so-called flux limiter function. The high-resolution methods applied to the solute transport solution are described as follows,

• UpWind

  $$\mathsf{\Phi }\left(r_{\mathrm{f}}\right)=0$$

  (24)
• MUSCL [42]

  $$\mathsf{\Phi }\left(r_{\mathrm{f}}\right)=\max [{\varphi }_{j,k}^n;\min (2r_{\mathrm{f}};0,5+0,5r_{\mathrm{f}};2)]$$

  (25)
• Superbee [43]

  $$\mathsf{\Phi }\left(r_{\mathrm{f}}\right)=\max [{\varphi }_{j,k}^n;\min (1;2r_{\mathrm{f}});\min (2,r_{\mathrm{f}})]$$

  (26)
• Ultimate Quickest [44]

  $$\mathsf{\Phi }\left(r_{\mathrm{f}}\right)=\max \left[{\varphi }_{j,k}^n;\min \left(\frac{1}{2}(1+r_{\mathrm{f}})+\frac{1}{6}(1-r_{\mathrm{f}})(1-2|C\left|\right);\frac{2}{1-\left|C\right|};\frac{2r_{\mathrm{f}}}{\left|C\right|}\right)\right]$$

  (27)

  where $C_r$ here is the Courant number.

The consecutive gradients $\left(r_{\mathrm{factor}}\right)$ play an important role in the numerical diffusion of the mass transport equation. Several ways to estimate the $r_{\mathrm{factor}}$ have been proposed [10,11,30,33,45,46]. The $r_{\mathrm{factor}}$ proposed in [10,11] gave the best results among all existing methods Ye et al. [11]. For the sake of simplicity, we adopted the $r_{\mathrm{factor}}$ proposed in [10], which we found to be easier to implement, as the $r_{\mathrm{factor}}$ is calculated only based on the concentration of neighbor cells.

The $r_{\mathrm{factor}}$ is calculated only for horizontal and vertical water-leaving faces. For vertical faces (horizontal flux), the consecutive gradient is given by

$$r_{j,k}^n=\frac{C_{i,k}^n-C_{r(i,j_o),k}^n}{C_{r(i,j),k}^n-C_{i,k}^n}$$

(28)

where ${\mathit{j}}_o$ represents the opposite face from face $\mathit{j}$. For horizontal faces (vertical flux), the $r_{\mathrm{factor}}$ is given by:

$$r_{j,k-\frac{1}{2}}^n=\frac{C_{i,k}^n-C_{i,k+1}^n}{C_{i,k-1}^n-C_{i,k}^n}$$

(29)

$$r_{j,k+\frac{1}{2}}^n=\frac{C_{i,k}^n-C_{i,k-1}^n}{C_{i,k+1}^n-C_{i,k}^n}$$

(30)

3.2. Bilinear Interpolator

For a structure rectangular grid, the bilinear interpolator [14] uses eight node points within a cell. The interpolated velocity component at $t^n$ can be found by

$$\begin{array}{c}{c}_{i-a,j-b,k-d}^n=\left(1-d\right)\left(1-a\right)\left[\left(1-b\right){\mathit{u}}_1^n+b{\mathit{u}}_2^n\right]+a\left[\left(1-b\right){\mathit{u}}_4^n+b{\mathit{u}}_3^n\right]+d\left(1-a\right)\hfill \\ \left[\left(1-b\right){\mathit{u}}_5^n+b{\mathit{u}}_6^n\right]+a\left[\left(1-b\right){\mathit{u}}_8^n+b{\mathit{u}}_7^n\right]\hfill \end{array}$$

(31)

where a, b, and d are the distances of the particle position at the end of a sub-time step $(x_p,y_p,z_p)$ normalized by the position of the nodes in each direction, where $u_1$ was set as the initial position $(x_{u1},y_{u1},z_{u1})$. For the following example (Figure 2), the normalized distances are set as

$$a=\frac{|x_p-x_{u1}|}{|x_{u1}-x_{u2}|};b=\frac{|y-y_{u1}|}{|y_{u1}-y_{u4}|};d=\frac{|z_p-z_{u5}|}{|z_{u1}-z_{u2}|}$$

(32)

We used two different procedures to select the eight node points used for bilinear interpolation of $\mathit{Fu}$ and $\mathit{Fw}$ at the end of the Lagrangian trajectory: Regarding $\mathit{Fu}$ (Figure 2 left), the interpolation process follows the steps:

(1)

For each time step, the particle starts at the barycenter of the j-th face at layer k, where the multistep backward Euler stream line is defined by linear interpolations to find the particle position at the end of a sub-time step (ELM step-i);

(2)

For each sub-time step, the particle position is analyzed in relation to the initial position (up-right, up-left, down-right, or down-left);

(3)

Eight node points are selected in the cell, with four points in layer k (the same layer where the particle stops) and four points in layer $\mathit{k}\pm 1$, depending on the particle position (up or down). The used points are always 4-face barycenters and 4 edge centers, except for the top and bottom cells, which uses 2-face barycenters, 4 edge centers, and 2 nodes;

(4)

The node indices are defined anticlockwise, where the first node is the initial position of the particle at time n+1;

(5)

Equation (31) is used to calculate the velocity component in the particle position in the sub-time step;

(6)

Steps (1–5) are repeated until the end of the Lagrangian trajectory.

For Fw (Figure 2 right), a similar procedure is adopted:

(1)

For each time step, the particle starts at the horizontal-face barycenter of the i-th element at layer $(\mathit{k}+\frac{1}{2})$, where the multistep backward Euler stream line is defined by linear interpolations to find the particle position at the end of a sub-time step (ELM step-i)

(2)

For each sub-time step, the particle position is analyzed in relation to the initial position (up or down and in the direction of one of the four quadrants)

(3)

Eight node points are selected in the cell, with four points in layer $(\mathit{k}+\frac{1}{2})$ (the same layer where the particle starts) and four points in layer $(k-\frac{1}{2})$, if the particle goes down, or $(\mathit{k}+\frac{3}{2})$, if the particle goes up. The used point always has a two-horizontal face barycenter, 4 edge centers, and 2 nodes;

(4)

The node indices are defined anticlockwise, where the first node is the initial position of the particle at time $\mathit{n}+1$.

(5)

Equation (31) is used to calculate the velocity component in the particle position in the sub-time step.

(6)

Steps (1–5) are repeated until the end of the Lagrangian trajectory.

3.3. Quadratic Interpolator

The Quadratic Lagrangian interpolation was adapted from Hodges et al. [16], who extended the 8-point upwind bilinear stencil to a 27-grid point stencil using at least eight computation cells. We proposed to use 27 node points inside a single computation cell, using the calculated velocities at the face barycenters and interpolated velocities at the nodes, edge centers, face barycenters, and element center (Figure 3).

The quadratic interpolation is more generic and uses the same procedure for the horizontal and vertical faces. First, for each time step, the particle starts at the face barycenter, where the backward stream line is defined in 10 sub-time steps by linear interpolation to find the departure position at time n (ELM step-i). Second, for each sub-time step, the interpolation follows three major steps, illustrated in Figure 3, to determine the Lagrangian polynomial coefficients (see, e.g., in [16]). The Lagrangian polynomial coefficients are given by

$${\mathcal{L}}^{\beta }=\prod _{\alpha =0\alpha \ne \beta }^2\frac{z^p-z_{n\pm \alpha }}{z_{n\pm \beta }-z_{n\pm \alpha }}$$

(33)

where each interpolation has three $\mathcal{L}$ coefficients, one for each node used in the interpolation line. $z^p$ is the particle position after the backtrack, $\beta$ is the analyzed node position in the Upwind direction, and $\alpha$ is the other two points, both assuming values of 1, 2, or 3.

In summary, the interpolation process follows these steps:

(1)

First, 9 z-direction interpolations are carried out (black lines in Figure 3). Each vertical interpolation generates a velocity component in the horizontal plane that passes through the z-position of the particle, given by

$${\tilde{\mathit{u}}}_{l+\gamma ,m+\psi }={\mathcal{L}}_{l+\gamma ,m+\psi }^1{\mathit{u}}_{l+\gamma ,m+\psi ,n}+{\mathcal{L}}_{l+\gamma ,m+\psi }^2{\mathit{u}}_{l+\gamma ,m+\psi ,n+1}+{\mathcal{L}}_{l+\gamma ,m+\psi }^3{\mathit{u}}_{l+\gamma ,m+\psi ,n+2}$$

(34)

where each ${\mathcal{L}}^{\beta }$ is multiplied by the respective velocity node and $\gamma$ and $\psi$ are the “l” and “m” displacement, respectively, to indicate which vertical line has been interpolated.

(2)

Three x-direction interpolations are carried out (orange line), using the estimated velocities found in step (1), resulting in three new velocity components (white diamonds in Figure 3), given by

$${\overline{\mathit{u}}}_{m+\psi }={\mathcal{L}}_{l,m+\psi }^1{\tilde{\mathit{u}}}_{l,m+\psi }+{\mathcal{L}}_{l+1,m+\psi }^2{\tilde{\mathit{u}}}_{l+1,m+\psi }+{\mathcal{L}}_{l+2,m+\psi }^3{\tilde{\mathit{u}}}_{l+2,m+\psi }$$

(35)

(3)

One interpolation is made to compute the y-direction displacement and find the final velocity for the particle at time $t^n$

$${\mathit{u}}_{i-a,j-b,k-d}^n={\mathcal{L}}_j^1{\overline{\mathit{u}}}_j+{\mathcal{L}}_{j+1}^2{\overline{\mathit{u}}}_{j+1}+{\mathcal{L}}_{j+2}^3{\overline{\mathit{u}}}_{j+2}$$

(36)

(4)

A new departure velocity is used to define the particle position in the next sub-time step.

(5)

Steps (1–4) are repeated until the end of the Lagrangian trajectory.

4. Numerical Experiments

The proposed numerical approaches were used in three consolidated benchmarks usually use to verification and validation of numerical models [5,6,15,19,47,48,49]. The first two experiments tested the numerical diffusion produced by the hydrodynamic solution, and the last experiment evaluated the numerical diffusion of the coupled hydrodynamic and solute transport solution. Each numerical experiment had a different purpose, as follows:

• Standing waves in a three-dimensional closed basin: this test case verifies the capability of the model to simulate 3D linear waves, including phase and amplitude representation [6,20,48]. The motion in the basin is caused only by the initial condition of the free surface. When the roughness, viscosity and diffusivity coefficients are set equal to zero, the motion of the free surface should not lose energy. However, a wave damping is caused by the numerical diffusion of the hydrodynamic solution. We evaluated the differences in the numerical diffusion between the bilinear and quadratic interpolations, applied in ELM step-ii, as well as the numerical diffusion considering the no-advection scheme.We also compared the mass conservation of the computation domain for each time step, the cumulative mass conservation over the course of the simulation, and the mean computation time of one time-step simulation.
• Wave propagation over a submerged bar: this was an experimental model idealized by Beji and Battjes [50], and has been frequently used to validate numerical models (e.g., [4,47,48,49,51,52,53]). The experiment was used to evaluate the accuracy of representing an irregular wave pattern caused by physical changes at the bottom, by comparing the quadratic and bilinear interpolations used in ELM step-ii.
• Gravity wave test: consists of a finite-amplitude deep-water standing wave in an inviscid fluid in a nonequilibrium situation, where the baroclinic pressure makes a major contribution to promote flow. The experiment evaluated the numerical diffusion in terms of density interface expansion, analyzing the difference between the combined uses of the interpolation techniques used in ELM step-ii and different flux-limiter schemes applied in a solute transport solution. This test case also evaluated the individual effect of each interpolation technique used in the hydrodynamic solution on the solute transport solution, using different flux limiters. We also compared the mean computation time of one time-step simulation.

All experiments were run using an Intel${}^{\circledR }$ Xenon${}^{\circledR }$ CPU-E5-1620 3.7 GHz computer with 32 GB of RAM memory in a Fortran based numerical model. We used Bias, root-mean-squared error (RMSE), Volume Error (%), and the Kling–Gupta (KGE) and Nash–Sutcliffe (NSE) coefficients of efficiency as metrics to compare the performance of the methods in each experiment.

4.1. 3D Standing Waves in a Closed Basin

In this experiment, the numerical diffusion was evaluated in terms of the cumulative wave damping, comparing the waves after 30 s of simulation. In the no-advection scheme, ${\mathit{u}}^{*}$, ${\mathit{v}}^{*}$, and ${\mathit{w}}^{*}$ (Equations (17) and (18)) were set directly equal to the horizontal and vertical face velocities. The numeric experiment set-up [6,20,48] are a closed cubic basin with edge of 10 m and wave amplitude set to $0.1$ m (Figure 4). The spatial domain was discretized using a regular grid with $0.5$ m resolution, resulting in 8000 computation cells. The time step was $0.01$ s and the total simulation time was 30 s. We evaluated the numerical diffusion produced by ELM, comparing the free-surface elevation at $x=y=0.25$ m with the exact solution. We also compared the mass conservation between bilinear and quadratic interpolators, neglecting the convective terms.

The analytical solution of the free-surface water elevation is given by

$$\eta =Acos\left(k_{x}x\right)cos\left(k_{y}y\right)cos\left(2\pi \frac{t}{T}\right)$$

(37)

where t is the time (the initial condition of the free surface may be obtained by setting $t=0$) and T is the wave period equal to $3.1$ s, with the wave number $kx=ky=n/L$ and the total wave number $k=\sqrt{k_x^2+k_y^2}=0.44\frac{rad}{m}$. The analytic solution for each velocity component is described as follows,

$$\mathit{u}=\frac{Ag{k}_x}{\omega }\frac{cosh\left[k_x(h+z)\right]}{cosh\left(k_{x}h\right)}sin\left(k_{x}x\right)cos\left(k_{y}y\right)sin\left(\omega t\right)$$

(38)

$$\mathit{v}=\frac{Ag{k}_y}{\omega }\frac{cosh\left[k_y(h+z)\right]}{cosh\left(k_{y}h\right)}cos\left(k_{x}x\right)sin\left(k_{y}y\right)sin\left(\omega t\right)$$

(39)

$$\mathit{w}=\frac{Ag{k}_x}{\omega }\frac{sinh\left[k_x(h+z)\right]}{cosh\left(k_{x}h\right)}cos\left(k_{x}x\right)cos\left(k_{y}y\right)sin\left(\omega t\right)$$

(40)

where $\omega$ is given by

$$\omega =\sqrt{gKtanh\left(Kh\right)}$$

(41)

The results showed that the numerical diffusion can be reduced when a higher-order interpolation is used (Figure 5). The numerical diffusion of the quadratic interpolator (measured in terms of the amplitude error compared to the analytical solution) is ~10 times smaller than the bilinear interpolator after 30 s of simulation. The bilinear scheme did not show a substantial improvement over the simulation neglecting the convective terms (the numerical diffusion was only 10% lower). The quadratic results show best performance between methods (Bias of 0.18 mm and KGE of 0.87, see

Table 1. Metrics between the analytical and simulated results each method for 30 s of simulation.

Metrics	Quadratic	Bilinear	No-Advection
RMSE (mm)	21.86	20.80	20.70
BIAS (mm)	0.18	0.52	0.56
Error (%)	27.22	25.94	25.79
KGE	0.87	0.66	0.64
NSE	0.90	0.91	0.91

), indicating best agreement with the analytical solution, despite similar RMSE and NSE found in the other methods.

The phase representation showed good agreement with the analytical solution for all tested methods, showing a cumulative phase error of only $0.3$ s at the end of the simulation, which is similar to previous studies [20,48]. As expected, the simulation neglecting the convective terms had a shorter computation time (~$1.65$ s) than the bilinear (~$2.23$ s) and quadratic interpolations (~$2.37$ s). The computation cost using a quadratic interpolation was slightly higher (~6%) than the bilinear interpolation.

The mass conservation analysis showed a more conservative behavior for the quadratic simulation (Figure 6). The simulation neglecting the convective terms showed a substantial decreasing pattern after 5 s of simulation, differently from when the bilinear was used, in which the pattern appeared after 20 s. The quadratic interpolation showed a more conservative pattern than the bilinear interpolation.

4.2. Wave Propagation over a Submerged Bar

A scheme of the experiment of the wave propagation over a submerged bar with an uneven bottom is seen in Figure 7 [50]. At the upward slope of the bar, the shoaling wave becomes nonlinear due to the generation of a higher bound harmonic. At the downward slope, the depth increases rather fast and these harmonics become free, resulting in an irregular wave pattern [47]. The numerical reproduction of this pattern has proven to be very demanding with respect to the accuracy of the computed dispersion frequency [4].

The computation domain has a total length of 30 m, with an initial undisturbed water level of $0.4$ m, and was discretized using a regular grid of $0.025$ m resolution. The time step was $0.005$ s and the total simulation time was 39 s. At the left boundary, a sinusoidal wave condition, with period $T=2$ s and amplitude $A=0.01$ m, was imposed to represent the wave generator of the original experiment. At the right outflow boundary, the experimental absorbed beach was computationally represented by a 5-m sponge layer with a combination of a sponge layer technique [54] and a Sommerfeld-type radiation boundary condition, applied to minimize wave reflection, given by

$${ϵ}_i=\left\{\begin{array}{c}\beta {\left(\frac{x_i-x_{io}}{l_i}\right)}^2\left(\frac{z_m-z}{z_m-z_M}\right)u_{i}if{x}_i\ge x_{io}\hfill \\ 0if{x}_i<0\hfill \end{array}\right.$$

(42)

where ${ϵ}_i$ is the sponge layer coefficient, $x_{io}$ is the initial point, and $l_i$ the total length. This term must be added in the right side of Equations (1) and (2).

The methods performance was compared using a few metrics (RMSE, BIAS, Volume Error, KGE, and NSE), and also through graphic visualization (Figure 8) of the wave pattern at each of the six stations. We evaluated the capability of the model to correctly represent the measured free-surface water elevation at between 33 and 39 s of simulation, comparing the bilinear and quadratic interpolation solutions with the experimental results for the “low-amplitude waves” (LW), generated by the shoaling process, and the “higher-amplitude waves” (HW).

The results indicated that the shoaling process on the upward slope was well described by both interpolation techniques, although the bilinear interpolation poorly represented the LW at station “A”. At the beginning of the downward slope, at station “B”, the bilinear simulation showed higher free-surface elevation, poorly representing the LW minimums (B at $34.7$ and $36.7$ s), which is physically incoherent with the measurements and simulations by other studies (e.g., Stelling and Zijlema [4]), due to an overincreased velocity in the de-shoaling process. The quadratic simulation more accurately represented the de-shoaling process, although it did not accurately represent the maximum height of HW (35 and 37 s). Moreover, the simulation without advection scheme was not able to represent the de-shoaling process and was numerically unstable, due to a larger overincrease in velocity at the same point as the bilinear interpolation overincreased it.

The bilinear overestimation of the free-surface level was propagated to the other stations. For the other four stations, the bilinear case better represented the HW maximum amplitudes than the quadratic simulation, but poorly represented the LW phase and shape (see “D” and “F” LW over the oscillation pattern for the bilinear interpolation).

The quadratic interpolator had the best performance among methods (see

Table 2. Metrics between simulated and experimental results for the six stations for bilinear and quadratic results.

M	Station a: x = 13.5 m		Station d: x = 17.3 m
	Bilinear	Quadratic	Bilinear	Quadratic
RMSE (mm)	3.57	2.17	5.18	4.29
BIAS (mm)	2.01	0.18	0.87	−0.56
Error (%)	47.52	23.65	83.78	60.00
KGE	−0.90	0.82	0.23	0.46
NSE	0.83	0.94	0.56	0.70
	Station b:x= 14.5 m		Station e:x= 19.0 m
RMSE (mm)	4.79	4.25	4.97	3.81
BIAS (mm)	0.14	−0.64	2.49	0.76
Error (%)	64.77	50.22	52.79	45.31
KGE	0.39	0.35	−6.12	−1.18
NSE	0.33	0.47	0.60	0.76
	Station c:x= 15.7 m		Station f:x= 21.0 m
RMSE (mm)	4.66	3.51	4.04	3.67
BIAS (mm)	−0.19	−1.25	1.34	−0.28
Error (%)	56.71	39.63	57.48	43.84
KGE	0.65	−1.07	−4.92	−0.26
NSE	0.67	0.81	0.69	0.75

). It better represented the amplitude and the wave pattern for all stations (NSE between 0.47 and 0.94, and RMSE between 2.17 and 4.29 mm). The results showed a considerable volume error (between 23% and 83%), and the bias indicates that, in general, there was an overestimation of the free-surface level in most cases (except for some methods at stations b, c, d, and f). In summary, the quadratic solution was better in interacting with uneven bottoms to represent complex nonlinear wave simulations. Similar results to the quadratic interpolation solution were reported by Chen [3], Walters [18], Dingemans [47], Cui et al. [53], and Monteiro and Schettini [20].

4.3. Gravity Wave

This benchmark set-up is used to evaluate the numerical diffusion of numerical models under sharp density gradients [19]. The equation that represents the initial density condition is

$$\rho (x,z)=-\frac{\mathsf{\Delta }\rho }{2}tanh\left[\frac{2{tanh}^{-1}\alpha }{k\delta }\left(kz-k\zeta +\frac{k}{2}\right)\right]$$

(43)

where $\alpha =0.99$, $k=\pi$ is the wave number, $ka$ represents the fluid interface inclination, $k\delta$ is the nondimensional interface thickness, and $k\zeta$ is the interface initial condition, given by [55]

$$k\zeta \left(x\right)=ka\left[\left(1-\frac{{\left(ka\right)}^2}{64}\right)cos\left(kx\right)-\frac{{\left(ka\right)}^2}{8}cos\left(3kx\right)\right]$$

(44)

The constant values follow those used in [19]; $ka=0.1$, $k\delta =0.05$, and the density difference between the layers is $\frac{\mathsf{\Delta }\rho }{{\rho }_0}=0.03$, with an analytical wave period of $9.43$ s. In this test, the bottom roughness and diffusion coefficients were both equal to zero. The spatial domain was a 1 m × 1 m, discretized using a regular grid of $0.0125$ m resolution, resulting in 6400 computation cells. The time step was $0.0236$ s and the total simulation time was two wave periods.

Because the vertical and horizontal diffusion were set equal to zero (i.e., pure advection transport), the initial and final shape of the gravity wave must be the same. We evaluated the numerical diffusion produced, neglecting the convective terms, and using both bilinear and quadratic interpolations in ELM step-ii combined with different flux-limiter schemes. The model outcomes (i.e., water density profiles) were compared with the exact solution. So as to evaluate the effects of different schemes and interpolation techniques on the numerical diffusion, we calculated the $L_{\mathrm{error}}$ norm of the density (

Table 3. Relative errors (%) for the different advection schemes after two wave periods for the bilinear and quadratic interpolators.

Flux Limiters	${\mathit{L}}_{\mathit{NoAdvc}}$	${\mathit{L}}_{\mathit{Bilinear}}$	${\mathit{L}}_{\mathit{Quadratic}}$	Fringer et al., 2005
UpWind	5.10	4.98	4.98	3.9
MUSCL	1.66	1.23	0.48	0.6
Ultimate-Quickest	1.40	0.94	0.46	0.4
SuperBee	1.36	0.88	0.30	0.4

), and compared the simulated vertical density profiles at $x=0.0$ m and at $x=0.4$ m with the exact solution.

The $L_{\mathrm{error}}$ was calculated using the following expression,

$$L_{\mathrm{error}}=\frac{{\sum }_{i=1}^{N_i}{\sum }_{k=1}^{N_k}|{\rho }_{i,k}^{2T}-{\tilde{\rho }}_{i,k}|\delta A_{i,k}}{{\sum }_{i=1}^{N_i}{\sum }_{k=1}^{N_k}\left|{\tilde{\rho }}_{i,k}\right|\delta A_{i,k}}$$

(45)

where $\tilde{\rho }$ is the density from the exact solution, ${\rho }^{2T}$ is the simulated density after two wave periods, and $\delta A_{i,k}$ is the area of grid cell $(\mathit{i},\mathit{k})$, which in our case is equal for all computation cells.

The comparison between the simulated density profiles and the exact solution is shown in Figure 9. Our findings indicated that the numerical diffusion produced by the hydrodynamic solution had a major influence on the total numerical diffusion.

The upwind scheme showed a higher numerical diffusion than the high-resolution schemes, considering the two interpolation techniques in ELM step-ii (RMSE ~2.8 times higher,

Table 4. Root-mean-squared errors (${\mathrm{kg}/\mathrm{m}}^3$) for the different flux limiter schemes after two wave periods, for the no-advection scenario, bilinear, and quadratic interpolators.

Flux Limiters	RMSE${}_{\mathbf{NoAdvc}}$	RMSE${}_{\mathbf{Bilinear}}$	RMSE${}_{\mathbf{Quadratic}}$
Upwind	2.08	1.97	1.97
MUSCL	1.16	0.72	0.38
Ultimate-Quickest	1.15	0.69	0.28
Superbee	1.16	0.70	0.24

). For the high-resolution schemes, the numerical diffusion was reduced mainly when combined with the quadratic interpolation in ELM (RMSE ~3 times smaller). A similar numerical diffusion pattern was found when the convective terms were neglected or when a bilinear interpolation was used in ELM step-ii, with a higher diffusive behavior in the middle of the spatial domain and an antidiffusive pattern near the boundary domain for the higher-resolution schemes, as seen in Figure 9 and Figure 10.

The upwind scheme showed a higher relative error ($L_{\mathrm{error}}$ norm ~5%) for all simulations. The results showed that the high-resolution methods substantially reduced the relative errors compared to the upwind scheme, and the Superbee performed best among the higher-resolution schemes. The quadratic and Fringer et al. [19] had similar results to the tested high-resolution methods; however, the hydrodynamic discretization and its influence in the result was not discussed in Fringer et al. [19]. The $L_{\mathrm{error}}$ indicated that the bilinear interpolation substantially improved the representation of the density profile compared to the simulations neglecting the convective terms. Regarding the computation cost, the higher-resolution methods had similar computation times for each time step, differing only with the interpolation technique, with $1.25$, $1.51$, and $1.62$ s for the no-advection, bilinear, and quadratic interpolations, respectively. The upwind scheme had a shorter computation time, with 1.01, 1.15, and 1.61 s for the no-advection, bilinear, and quadratic interpolations, respectively.

5. Discussion

The numerical experiments in the hydrodynamic solution showed that the quadratic interpolation method, using 27 node points in a single computation cell, substantially reduced the numerical diffusion in the hydrodynamic solution, which had a positive effect on the solute transport solution. This is the first study to verify and validate the proposed quadratic interpolation method in ELM step-ii.

The results of the first experiment showed a good amplitude representation of the waves, with a satisfactory phase representation only $0.3$ s slower than the analytical result. This was expected, as the phase representation is more related to the nonhydrostatic pressure and vertical momentum discretization [48,52]. The second experiment successfully validated the proposed interpolation, which proved to be capable of representing complex wave problems. The bilinear interpolator applied at ELM step-ii had a numerical diffusion ~10 times higher than the quadratic interpolation. Moreover, the bilinear interpolation did not yield a substantial improvement in terms of numerical diffusion, compared with the simulations neglecting the convective terms. The quadratic interpolation also substantially improved the mass conservation over the course of the simulation, indicating that a high-resolution method can be applied to find mass conservative solutions in free surface simulations in complex wave problems, including shallow or deep waters. Higher numerical diffusion was also observed when the bilinear interpolator in ELM step-ii was used with different flux-limiter schemes in the mass transport solution [7], indicating that high-resolution schemes can be successfully used to attenuate the numerical diffusion for coupled hydrodynamic and transport solutions. The set of numerical experiments showed that the quadratic interpolation is a powerful and promising method to reduce numerical diffusion, with a slight increase in computation cost related to a bilinear interpolator (7.3% longer), and also can be applied in 1D, 2D, or 3D models. Other key factors related to the numerical solution are also responsible for the numerical diffusion, such as the fractional step error [3], and are caused by splitting the pressure solution into hydrostatic and nonhydrostatic modules, which was appropriately treated here following [56] and the interpolation method used to define the stream line trajectory in ELM step-i [21].

In this study, we used different flux-limiter schemes to efficiently solve the mass transport equation, which also proved to be an accurate alternative to reduce the numerical diffusion when coupled with a quadratic interpolation in ELM. The last numerical experiment showed that the Upwind scheme has serious numerical diffusion problems ($L_{\mathrm{error}}$ norm ~5%), and can be considered inappropriate for modeling stratified basins. The SuperBee flux limiter showed the best results among the high-resolution methods used in this study ($L_{\mathrm{error}}$ norm ~0.3% with quadratic interpolation). Although further investigation is still needed to indicate the best flux limiter, the SuperBee has been used as a default in some models (see, e.g., in [11,16,30,31]), has proven to be capable of simulating stratified basins with or without sharp gradients, and has also been mass conservative. Despite the satisfactory results, the flux limiters used here are only nonlinear second-order schemes [9,19]. Numerical diffusion can be reduced even more by a different linear or nonlinear higher-order discretization, but this usually requires a complex numerical solution, which is computationally expensive and also may be vulnerable to unphysical spatial oscillations (wiggles) under some circumstances [44]. The flux limiter is a simple approach that is easily implemented, especially when Kong’s $r_{\mathrm{factor}}$ is applied.

The results indicate that high-resolution schemes are suitable for reducing numerical diffusion, but combining their use in hydrodynamic and solute transport solutions can improve the overall result even more. When the Upwind scheme was used, the differences between interpolation techniques were overcome by the higher Upwind numerical diffusion, although when flux-limiter methods were also used, the differences were clearer. The combined use of high-resolution methods shows that the interpolation technique in the hydrodynamic solution has a substantial effect on the mass transport solution, and that, despite the improvement over Upwind, applying linear interpolators at the ELM step-ii still yields a relative error that is 2 to 3 times larger than with a nonlinear interpolator, with a similar behavior to the no-advection solution, which may generate unsatisfactory results for the transport solution in cases of sharp stratification.

A higher diffusive behavior in the middle of the spatial domain of gravity wave than in the nearby boundary domain was found for all simulations, as expected. Fringer et al. [19] showed that in regions of the flowfield where the local Courant number is reduced (the middle), the interfacial diffusion is increased, which makes a zone that is sensitive to the velocity field. The results of the gravity-wave experiment were strongly influenced by the diffusion-over-time steps, because, as the interface became more diffuse, the wave period and velocity field changed; that is, the use of a more diffusive method in a hydrodynamic solution may generate cumulative errors, leading to unsatisfactory results. The quadratic interpolation proved to better reduce the numerical diffusion in these critical areas, due to the less diffusive behavior in the hydrodynamic solution, and better predicted the sharp change of the velocity field in the middle region of the wave. Moreover, the results for the quadratic interpolation indicated that it performed better than the bilinear interpolation and similarly to the results reported by Fringer et al. [19] in the transport solution. The bilinear interpolation, although it had a smaller relative error than the no-advection solution (1.5 times smaller), showed the same unsatisfactory diffusive behavior in the middle region.

The proposed combined use of the quadratic interpolation applied at ELM step-ii and the flux-limiter technique substantially reduced the numerical diffusion in solving mass transport problems, showing that high-resolution methods must be implemented in the numerical solution to properly simulate more complex real situations. These methods can be easily implemented, as proposed here, and applied in 1D, 2D, or 3D approaches, which are mass conservative and do not introduce excessive numerical and stability problems, as discussed in Casulli and Lang [5].

These methods can be easily implemented, as proposed here, can be applied in 1D, 2D, or 3D models, and they are mass conservative, with no additional stability problems than the inherent numeric solution, as discussed in Casulli and Lang [5].

6. Conclusions

The combined use of high-resolution methods (quadratic interpolation and flux-limiter functions) proved a suitable alternative to reduced numerical diffusion, and with low cost of implementation in relation to higher order discolorations. The analyses in a coupled hydrodynamic and solute transport numerical model allowed to understand how the numerical diffusion at one solution may affect the other. We found that the numeric diffusion at the hydrodynamic solution promoted by low-resolution methods (low-order interpolation at ELM) may have a substantial impact on the solute transport solution, even if a high-resolution method were applied at the solute transport solution (relative error and RMSE ~3 times higher). On the other hand, when a low-order method is used in the solute transport solution (Upwind scheme), the numerical diffusion differences between methods in the hydrodynamic solution were overcome by the higher Upwind numerical diffusion. Thus, to accurately modeling stratified flows in real situations, the combined use of high-resolution methods is mandatory.

Numerical models that use low-resolution methods usually neglected the diffusive part of the transport equation (Equation (7)), employing diffusivity coefficients equal to zero and assuming the numerical diffusion as the physical diffusion, which may difficult the calibration of the model. The applied high-resolution methods allows to use these diffusivity coefficients as parameters to the model calibration. Therefore, we recommended the application of the implemented methods at real situations to evaluated the model’s performance in the representation of the temperature dynamics of stratified flows in deep lakes or reservoirs.