The mR scheme to the shallow water equation with horizontal density gradients in one and two dimensions

1.   Introduction

Recently, geophysical fluid dynamics modeling has grown in prominence due to the importance of forecasting and understanding the time development of a wide range of atmospheric and oceanic flows. The characteristics of the flows and the corresponding scales are used to construct models for simulating geophysical flows. A significant horizontal length scale in proportion to depth distinguishes several geophysical and atmospheric processes. That is, they are shallow, and shallow water equations are suited for explaining their time evolution in those conditions ^[1]. The one-dimensional shallow water model in Eulerian coordinates is considered ^[2,3].

One way to depict fluid flow is as a multi-layered flow where one layer flows over another to generate a new layer. The shallow water is a narrow layer of constant density fluid in hydrostatic equilibrium, bounded from below by a rigid surface and from above by a free surface. The shallow water model is one of the most vital models of the hyperbolic systems. In applied science and new physics, nonlinear hyperbolic systems of conservation laws are crucial for creating mathematical representations of a variety of natural processes ^[4,5,6,7,8]. It is the most fundamental layer example of an incompressible fluid moving across a free surface ^[9]. Utilizing this approach, channels, hydraulic jumps, tsunamis and reservoirs may all be simulated. Because the topography may be irregular and the geometry may be complicated, replicating these natural flow problems is difficult. Some bottom topographies need moving meshes in Eulerian coordinates that are stationary meshes in mass Lagrangian coordinates ^[2].

Various mathematical methods for finding analytical and numerical solutions to nonlinear partial differential equations (NPDEs) have been developed over the years ^{[10,11,12,13,14,15]}. In the ongoing work, we take into account a system of flows for shallow water where horizontal density fluctuations are considered. The Euler equations are vertically averaged to produce this model. We develop the mR scheme to solve the shallow water flows via horizontal density fluctuations. This approach has predictor and corrector stages. The Riemann invariants and limiters principles serve as the foundation for the control parameter for numerical diffusion through the predictor stage. The balance conservation equation is recovered in the corrector stage, see ^{[16,17,18,19,20,21]}. In most common schemes, the numerical flux was estimated using the Riemann solution. In contrast to earlier schemes, the intriguing property of the modified Rusanov (mR) scheme is to evaluate the numerical flux in the absence of the Riemann solution ^[22]. As long as the prerequisite for the canonical Courant-Friedrichs-Lewy $(CFL)$ is obeyed, this scheme is linearly stable.

The remainder of this article's framework is constructed as follows. Section 2 introduces the mathematical model for shallow water equation with horizontal density gradients. Section 3 provides the 1D mR method to solve the proposed model. Section 4 shows various 1D test cases for investigating the development mechanisms of constructed waves. Section 5 depicts the 2D shallow water with horizontal density gradients. Section 6 offers various numerical test cases to check the accuracy and performance of the proposed scheme in 2D. Conclusions and observations about the current results are presented in Section 7.

2.   Mathematical model

We consider the shallow water model with influences of horizontal density gradients, which given as follows ^[1]:

$h$ is the height of the water above the bottom, $g$ is the gravity acceleration, $u$ is the water velocity, $Z$ is the function that describes the bottom topography and $\rho$ is the vertical averaged density. The above system can write in the conservation form as the following

The vector-valued functions $W$; $F(W)$ and $Q(W)$ in $R^3$ are

It is clear that system of Eq (2.2) is strictly hyperbolic, and the associated eigenvalues are

3.   Modified Rusanov scheme

The integration of the (2.2) on the domain $[x_{i-1/2}, x_{i+1/2}] [t_n, t_{n+1}]$, which gives the finite volume scheme, which equivalent the corrector stage and we can be written as the following

$W_i^n$ is the average value of the solution $W$ over $[x_{i-1/2}, x_{i+1/2}]$ at time $t_n$ i.e.,

and $F(W_{i\pm1/2}^n)$ represents the numerical flux at the time $t_n$ and place $x = x_{i\pm1/2}$. Generally, in the finite volume scheme (3.1), the Riemann solution at the cell interfaces $x_{i\pm1/2}$ is required for the building of the numerical fluxes $F(W_{i\pm1/2}^n)$. The objective is to create the intermediate states, which are denoted by $W_{i\pm 1/2}^{n}$. To achieve this, we must integrate the Eq (2.2) through the interval $[t_n, t_{n}+\theta^n_{i+1/2}]\times[x^-, x^+]$. We achieve the intermediate state, which has the following written form:

whereas, the distance measurements, $\Delta x^-$ and $\Delta x^+$ can be stated as follows.

$Q^n_{i+\frac{1}{2}}$ refers to an approximation of the source term $Q(W)$

When, we substitute $x^-$ by $x_i$ and $x^+$ by $x_{i+1}$ in the Eq (3.2), we get the following equation.

whereas $W_{i+1/2}^n$ is define as the following are additional alternatives for $x^+$ and $x^-$. According to the stability analysis in ^[16] for conservation laws, we selected the control parameter $\theta^n_{j+1/2}$ for completing our scheme, and we are able to pick the parameter $\theta^n_{j+1/2}$ as follows:

The Eq (3.2) shows that the

whereas

and the control parameter is $\alpha^n_{i+1/2}$. The kth eigenvalue of (2.2) is represented by $\lambda^n_{k, i}$ and $k$ is the number of of eigenvalues of the system (2.2). As a result of this discussion, we can choose the control parameter in the following ways:

ⅰ) $\alpha_{i+1/2}^n = 1$, the mR scheme becomes the upwind scheme in the linear case with this decision ^[16].

ⅱ) $\alpha_{i+1/2}^n = \frac{\Delta t}{\Delta x}S_{i+1/2}^n$, by making this decision, the suggested scheme return to the Lax-Wendroff scheme.

ⅲ) $\alpha_{i+1/2}^n = \tilde\alpha_{i+1/2}^n = \frac{S_{i+1/2}^n}{s_{i+1/2}^n}$, with this decision, the scheme is transformed into a first-order scheme.

ⅳ) All calculations in this paper were made using the following formula for

where $\Phi_{i+1/2} = \Phi\left(r_{i+1/2}\right)$ is a function that limits slope, with

At last, we write the propose scheme for shallow water equation with horizontal density gradients as the following

The proposed scheme is well-balanced in the sense of ^[23,24,25], if the source term $Q^n_i$ in the corrector stage is approximated in such a way that the still-water balance ($C$-property) ^[26] is obeyed, if the condition

is satisfied the steady state solutions, then, we said the numerical scheme for the shallow water equation with horizontal density gradients satisfies the $C$-property.

Proof of the exact $C$-property. In (2.1) and (2.2), we assume that $u = 0$, which represents a stationary flow in rest. Then Eq (2.1) can be written as follows

We applied the predictor stage in (3.10) for the system (3.12), we have

During the corrector stage, the solution is updated as

To obtain a stationary solution $W^{n+1}_i = W^n_i$, the sum of the discretized flux gradients and source terms in Eq (3.12) must be equal to zero.

Then, the previous expression (3.15), which is equivalent to

Therefore, if the source term in the corrector stage of Eq (3.16) is discretized as follows

the proposed scheme satisfied the $C$-property.

4.   Numerical results

We present several test cases for numerical simulation of the shallow water equation with horizontal density gradients. We selected the stability condition ^[16] in the manner described below in order to demonstrate the effectiveness and accuracy of the suggested finite volume scheme

where the chosen constant $CFL$ equals 0.5.

4.1.   Test $1$

With the following initial conditions ^[1], this case's solution includes a rarefaction wave; contact discontinuity; shock wave:

This test case is simulated using modified Rusanov and classical Rusanov, when compared to the reference solution produced by applying the classical Rusanov technique with a smaller mesh size of $30000$ uniform cells over the range $[-10, 10]$ and time t = 1 s with gravity ($g = 1$). Figures 1–3 show the results of water height, velocity, density, quantity of movement, pressure and parameter of control. For the another the domain $[0, 12]$ and gravity $9.81$, we simulate this test case with a more refined mesh of $20000$ uniform cells. In this case, Figures 4–6 show the results.

4.2.   Test $2$

According to the initial conditions listed in ^[4], this example has been suggested:

The mesh has $300$ cells and the domain of computation is $[-10, 10]$ with a final time of t = 0.5 s, by applying the mR technique and the classical Rusanov method. Figure 7 displays the numerical results for the water height and velocity, which are contrasted with the reference solution that was produced utilizing the classical Rusanov method and a finer mesh made up of $30000$ uniform cells. We find that the results from the modified Rusanov scheme are more precise than those from the classical and modified Rusanov schemes with alpha constant $(\alpha^n_{i+\frac{1}{2}} = 1.4)$. Figure 8 shows the density and quantity of movement. Figure 9 shows the pressure and $\alpha^n_{i+\frac{1}{2}}$.

4.3.   Test $3$ non homogeneous

We consider the same case in the first example with source term as follows

This test case is simulated using classical Rusanov and modified Rusanov, when compared to the reference solution produced by the classical Rusanov method using a smaller mesh size of $20000$ cells over the interval $[0, 12]$ and time t = 0.4 s, with gravity ($g = 9.81$). Figures 10–12 show the results of water height, velocity, density, quantity of movement, pressure and parameter of control.

4.4.   Test $4$ non homogeneous

We consider the same case in the first example with source term and we can write it as follows

Modified Rusanov and classical Rusanov are used to simulate this test case, in comparison to the reference solution produced by the classical Rusanov method using a smaller grid of $20000$ uniform cells over the interval $[0, 12]$ and time t = 0.4 s, with gravity ($g = 9.81$). Figures 13–15 show the results height, velocity, density, quantity of movement, pressure and parameter of control. Also, Figure 16 shows the effect of parameter of control ($\alpha = 1, 1, 5, 2, 2.5$) on the diffusion of density and quantity of movement.

5.   Two-dimensional shallow water with horizontal density gradients

The shallow water with horizontal density gradients can be expressed as follows in two dimensions:

with

and

where $\rho, h, u, v$ are the density, water height and velocity respectively

5.1.   Two-dimensional problems with modified Rusanov scheme

In order to deducing the modified Rusnov scheme in two dimension, we integrate the Eq (5.1) on a generic control volume $c_i$ as shown in Figure 17, which are produced when the domain is divided into several control volumes, presents

where $A_i$ refers to the cell's area $c_i$, $K(i)$ is the index set of adjacent triangles that share an edge with the cell $c_i$, $\vec{n}_{ij}$ denotes the unit outward normal vector to the surface that surround the control volume, $\int_{\gamma_{ij}}\mathcal{F}(W, \vec{n}_{ij})d\sigma = \phi(W^n_i, W^n_j, \vec{n}_{ij})mes(\gamma_{ij})$ is the numerical flux, $mes(\gamma_{ij})$ is the interface between cells $c_i$ and $c_j$, with $W^n_i$ and $W^n_j$ are the left and right values of $W$ at cell $c_i$ and $c_j$.

It is possible to express the modified Rusanov scheme as

where $S^n_{ij}$ is the local Rusanov speed represented by

and the local control parameter $\alpha^n_{ij}$ is chosen to meet the stability condition; for more information, see ^[16]. The eigenvalues of the Jacobien matrix (5.1) are$\lambda^n_{p, i}$ and $\lambda^n_{p, j}$.

5.2.   Calculation of $W^n_{ij}$

Here, we aim to rewrite the numerical flux to resemble the 1D situation. To achieve this, we write $\phi(W^n_i, W^n_j, \vec{n}_{ij}) = \mathcal{F}(W^n_{ij}, \vec{n}_{ij})$, $W^n_{ij}$ is determined at the predictor stage. To determine the $W^n_{ij}$, we project the system (5.1) on the local cell's outward normal $\mu$ and tangential $\mu^{\bot}$ as shown in Figure 17. Then, we have

where $U = \vec{V}\cdot\eta = un_x+vn_y$ represent the normal speed and $V = \vec{V}\cdot\mu^{\bot} = -un_y+vn_x$ represents the tangential velocity. The predictor phase (5.5) expressed as

with $\Delta_{ij}(X) = X_j-X_i$ and $M_{ij}(X) = \frac{1}{2}(X_i+X_j)$. Where $\tilde{ \rho h}_{ij}$ is the linear interpolation between $(\rho h)_i$ and $(\rho h)_j$ given by $\tilde{\rho h}_{ij} = \frac{A_i (\rho h)_i+A_j (\rho h)_j}{A_i+A_j}$. The solution $W^n_{ij}$ is recovered via the transformation

and

6.   Numerical results in two dimensions

To test the precision and effectiveness of the suggested scheme in 2D, we provide several test cases. Due to presented CFL condition ^[16], the theoretical maximum stable time step $\Delta t$ is specified.

the constant Currant number $Cr$ is chosen to be less than unity, in this case $Cr = 0.2$ in all simulations, and the time step $\Delta t$ is changed using Eq (6.1), with $M = \underset{i, j}{max}(S^n_{ij})$, $m = \underset{i, j}{min}(s^n_{ij})$, where $S^n_{ij}$ is the local Rusanov speed, $A_i$ refers to the cell's area $c_i$, the control parameter for the suggested scheme is $\alpha^n_{ij} = 1.8$, and $|\delta C_i|$ denotes the perimeter of the cell $C_i$.

We consider the this test case with and without source term. The initial conditions are

At the end, we simulate this test case using the mR scheme on the domain $[0, 12]\times[0, 1]$ meshgrid $150\times 15$ equivalent 1903 cells and final time $t = 0.4$. Figures 18–20 display water height, velocity and density with out and with source term.

7.   Conclusions

The mR scheme was proposed and applied to numerically simulate the shallow water equation with horizontal density gradients. This approach can properly represent discontinuance profiles and prevents numerical diffusion in the solution. In contrast to earlier schemes, the mR scheme has the intriguing ability to compute the numerical flow in the absence of the Riemann solution. This method, in fact, can be utilized as a box solver for a wide range of other conservation law models. Numerical results in one and two dimensions were introduced to depict the efficiency of the presented scheme. The Rusanov approach and analytical solutions were compared with the mR scheme. Several test cases were introduced in 1D and 2D. The numerical results display the high resolution of the proposed mR scheme and show that it can provide precise and effective simulations for the shallow water equation with horizontal density gradients.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors extend their appreciation to the Deputyship for Research & Innovation, Ministry of Education in Saudi Arabia for funding this research work through the project number 445-9-164.

Conflict of interest

The authors declare that they have no competing interests.