Stability analysis and convergence rate of a two-step predictor-corrector approach for shallow water equations with source terms

1.   Introduction

Most open-channel flows of interest in the physical, hydrological, biological, engineering and social sciences are unsteady and can be considered to be one-dimensional ($1$D). In this paper, we are interested in the numerical solutions of one-dimensional shallow water equations with source terms introduced in ^[1,2] and still widely used in modeling flows in rivers, lakes and coastal areas as well as atmospheric and oceanic flows in some regimes. In the case of a prismatic channel, the shallow water equations with source terms read as follows

where the bottom (or bed) slope $(S_{0})$ and the friction slope $(S_{f})$ (see ^[1]) are defined as

where $v = v(t, x)$ is the lateral inflow per unit length along the channel, $T_{1}$ represents the time interval length, $L$ denotes the rod interval length, $A = A(t, x)$ and $Q = Q(t, x)$ are the cross-section and the discharge, respectively, $g$ denotes the acceleration of gravity. $T = T(t, x)$ represents the top width assumed to be constant, $\overline{\tau}$ designates the average shear stress on the water from the channel boundary, $\rho$ is the fluid density and $P = P(x, y(t, x))$ denotes the wetted perimeter (i.e., the length of the boundary of the cross-section that is underwater for a given height of water $(y)$). As in ^[1], the conveyance for a compact channel is defined as

where $R = A/P$ (see for example, ^[1,3]) denotes the hydraulic radius and $n_{1}$ represents the manning's roughness coefficient.

Since the parameters $g$ and $T$ are positive constants, it is easy to see that $g\frac{A}{T}\frac{\partial A}{\partial x} = \frac{\partial }{\partial x}\left(\frac{gA^{2}}{2T}\right).$ Using this equation together with equality $\frac{\partial }{\partial x}\left(\frac{gA^{2}}{2T}\right)+\frac{\partial}{\partial x}\left(\frac{Q^{2}}{A}\right) = \frac{\partial}{\partial x}\left(\frac{gA^{2}}{2T}+\frac{Q^{2}}{A}\right)$, the second equation of the system $(1.1)$ is equivalent to

This equation combined with the first equation in $(1.1)$ results in

which can be rewritten as

where $W = (A, Q)^{T}, \; F = (Q, \frac{gA^{2}}{2T}+\frac{Q^{2}}{A})^{T},$ and $S = (v, gA(S_{0}-S_{f}))^{T}.$ Equation $(1.4)$ emphasizes the conservative character of system $(1.1)$.

The one-dimensional shallow water equations with source terms $(1.4)$ are highly nonlinear and therefore do not have global analytical solutions ^[1]. When solving the system of balance laws $(1.4)$ numerically, one typically faces several difficulties. One difficulty stems from the fact that many physically relevant solutions of $(1.4)$ are small perturbations of steady-state solutions. So, using a wrong balance between the flux and geometric source term in Eq $(1.4),$ the solution may develop spurious waves of a magnitude that can become larger than the exact solution. Another drawback occurs when the cross-section is very small. In that case, even small numerical oscillations in the computed solution can result in a very large discharge, which is not only physically irrelevant, but cause the numerical scheme to break down. To overcome these numerical challenges, one needs to use a numerical scheme that is both well-balanced and positivity-preserving.

A number of well-balanced and positivity-preserving numerical methods frequently used in the models based on the shallow water equations have been proposed in the literature ^[4,5], or on the boussinesq equations, which are reduced to shallow water equations, in order to simulate breaking waves ^[6,7,8,9]. Although the MacCormack scheme is less accurate than the more recent methods, it is commonly used for engineering problems due to its greater simplicity. So, we have to approximate the exact solution of problem $(1.4)$ by a numerical method based on a two-level predictor-corrector scheme. This algorithm lies in the class of higher order finite difference methods (temporal second-order convergent and fourth-order accurate in space) which provide an effective way of joining the spectral method for accuracy and robust characteristics of finite difference schemes. For example, to compute unsteady flow specifically in the presence of discontinuity, inherent dissipation and stability, one such widely used method is the MacCormack method ^[10]. This technique has been used successfully to provide a time-accurate solution for fluid flow and aeroacoustics problems. The applications of this technique to $1$D shock tubes and $2$D acoustic scattering problems provide good results when compared with the analytical solution. The original MacCormack introduces a simpler variation of the Lax-Wendroff scheme which is basically a two-step scheme with second order Taylor series expansion in time and fourth order in spatial accuracy ^[10,11]. This algorithm is computationally efficient and easy to implement which can be appropriate to obtain reliable results. By using this scheme with two nodes, the flow field can be simulated for unsteady flows especially for shallow water problems in the presence of discontinuity and strict gradient conditions. Furthermore, to capture the fluid flow in transition over long periods of time and distance, numerical spatial derivative are required to be determined in a few grid points while error-controlled can be accurately computed. The authors ^{[12,13,14,15,16,17]} extended the MacCormack scheme ^[10] to an implicit-explicit scheme (by coupling the original MacCormack approach with Crank-Nicolson method), an implicit compact differencing scheme and a three-level time-split MacCormack procedure (by splitting the derivative operator of the scheme into one-sided forward and backward difference formulations). The one-sided nature of the time-split MacCormack approach is an essential advantage especially when severe gradients are present.

In ^[18,19,20] the authors compared the Lax-Wendroff scheme to many numerical methods of high order accuracy, such as, the linear Central Weighed Essential Non-Oscillatory (CWENO) scheme which is superior to full nonlinear CWENO method, to high-resolution TVD conservative procedures along with high order Central Schemes for hyperbolic systems of conservative Laws ^[3,21] and to Central-upwind schemes for the shallow water system ^[22]. In a search for stable and more accurate shock capturing numerical approach, they observed that the Lax-Wendroff technique is one of the most frequently encountered in the literature related to classical shock-capturing schemes. However, difficulties have been reported when trying to include source terms in the discretization and to keep the second order accuracy at the same time ^[23]. The proposed two-level solver which can be considered a predictor-corrector version of the Lax-Wendroff algorithm provides a reasonably good result at discontinuities. The developed method is more easier to apply than the Lax-Wendroff method because the Jacobian does not appear. For linear equations, both algorithms provide similar amplification factors and stability constraints (for instance, see ^[24], P. $202$–$206$). It is worth noticing to mention that the solutions obtained for a given problem at the same courant number are different from those obtained using the Lax-Wendroff scheme. This situation is due to the switched differencing in the predictor, corrector and the nonlinear nature of the governing equations. It should be noted that reversing the differencing in the predictor and corrector steps leads to quite different results. Furthermore, the explicit MacCormack time discretization for nonlinear Burgers equations (which can serve as model equations for a wide set of nonlinear PDEs: Navier-Stokes equations, Stokes-Darcy equations, Parabolized Navier-Stokes equations, ...) give a suitable stability restriction which should be used with an appropriate safety factor (see ^[24], P. $227$–$228$). As regards the proposed predictor-corrector formulation, like other explicit approaches, it requires a time-step limitation. In general, the maximum time-step allowable in the natural MacCormack scheme applied to linear hyperbolic equations is limited by the CFL condition, as for all explicit finite difference methods. However, the considered overland flow equations are nonlinear and a rigorous stability analysis of numerical techniques is exceedingly difficult. The source terms place additional and problem-dependent restrictions on the maximum admissible time-step for stability. Therefore, the CFL condition can only be considered as a general guideline here, and the maximum allowable time-step for any particular problem will be less than predicted by the CFL condition and determined by numerical experimentation (see ^[25], page $223$).

Most recently, several researchers have deeply analyzed implicit/explicit (for example: Crank-Nicolson/MacCormack) approach, compact finite difference methods, two-level factored Crank-Nicolson scheme, three-level time-split MacCormack algorithms, central upwing schemes and Lax-Wendroff technique in the approximate solutions of mixed Stokes-Darcy's model, evolutionary reaction-diffusion equation, nonstationary convection-diffusion equation, time-dependent convection-diffusion-reaction problems, unsteady coupled Burgers'equations, Navier-Stokes equations, time fractional equations and shallow water problems. For more details, we refer the readers to ^{[5,9,12,13,17,22,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41]}. In this paper, we are interested in a computed solution of a nonlinear system $(1.1)$, but in the sense of a linear stability condition and convergence rate of the numerical scheme. In particular, we consider the case where the channel is prismatic and the interesting result is that the algorithm is second-order accurate in time and spatial fourth-order convergent, whereas the stability limitation is not similar to the CFL condition widely studied in the literature for hyperbolic PDEs. However, while the stability requirement is highly unusual, the result has a potential positive implication since the stability constraint controls the CFL condition. Indeed the nice feature is that, as required in a stability context, we normally find a linear stability requirement which can be considered as a necessary condition. In addition, it follows from this analysis that instability occurs when $\Delta t$ is greater than $\Delta t_{\max} = (\Delta t)_{CFL}.$ More specifically, the attention is focused in the following three items:

(i1) Full description of a two-step predictor-corrector method for solving one-dimensional shallow water equations with source terms.

(i2) Stability analysis of the proposed algorithm: this item together with item (i1) represent our original contributions and they improve some works studied in the literature ^[3,12,42,43].

(i3) Numerical experiments which considers the convergence rate of the method, the simulation of the numerical solution together with the analytical ones, and regarding the effectiveness of the proposed algorithm according to the theoretical analysis given in the first two items.

The paper is organized as follows. Section $2$ deals with a full description of a two-step predictor-corrector formulation for one-dimensional shallow water equations with source terms. A necessary condition of stability of the proposed numerical scheme is deeply analyzed in Section $3$. In Section $4,$ some numerical experiments which consider the convergence rate of the developed approach and some simulations are presented and discussed. We draw the general conclusion and present the future works in Section $5$.

2.   Full description of a two-level explicit predictor-corrector method

In this section, we give a detailed description of a two-step explicit predictor-corrector algorithm for the system of nonlinear Eq $(1.4).$ The proposed scheme uses the forward difference in the predictor step while the corrector step considers the backward difference. Since the aim of this section is to analyze the linear stability of the method, without loss of generality we should use a constant time step $\Delta t$ and mesh size $\Delta x$. However, this assumption does not compromise the result on the stability. Let $K_{1}$ and $M_{1}$ be two positive integers. Set $x_{r} = r\Delta x, \; t^{i} = i\Delta t$ be the discrete points and let the superscript and subscript be the time level and space level, respectively, of the approximation. We denote $W_{r}^{i} = (A_{r}^{i}, Q_{r}^{i})^{T}$ be the approximate solution of equations $(1.4),$ provided by the constructed two-step MacCormack technique and $W(t^{i}, x_{r}) = (A(t^{i}, x_{r}), Q(t^{i}, x_{r}))^{T}$ be the analytical ones. Furthermore, the domain $\Omega = (0, L)$ is partitioned into $M_{1}+1$ grid points $\{x_{r}:\text{ }r = 0, 1, ..., M_{1}\},$ whereas the time interval $(0, T)$ is subdivided into $K_{1}+1$ discrete points $\{t^{i}:\text{ }i = 0, 1, ..., K_{1}\}.$ Using this, a full description of the proposed method for solving the system $(1.4)$ reads: Given $W^{i}_{r}$, find an approximate solution $w^{i+1}_{r},$ for $0\leq i\leq K_{1}-1$ and $0\leq r\leq M_{1},$ satisfying

Predictor step: solve Eq $(2.1)$ for predicted value

Corrector step: use the predicted value obtained in Eq $(2.1)$ to compute the corrected one

Lemma 2.1. Let $K_{1}$ and $M_{1}$ be two positive integers. Setting $t^{i} = i\Delta t$, for $i = 0, 1, ..., K_{1},$ and $x_{r} = r\Delta x,$ for $r = 0, 1, ..., M_{1},$ where $\Delta t$ and $\Delta x$ are time step and mesh size, respectively. The proposed numerical scheme for solving the system of nonlinear Eq $(1.4)$ is given by: for $i = 0, 1, ..., K_{1}-1,$ and $r = 0, 1, ..., M_{1}-1,$

Predictor step: solve Eqs $(2.3)$ and $(2.4)$ for predicted values $A_{r}^{\overline{i+1}}$ and $Q_{r}^{\overline{i+1}}$

Corrector step: use $A_{r}^{\overline{i+1}}$ and $Q_{r}^{\overline{i+1}}$ obtained in Eqs $(2.3)$ and $(2.4)$ to compute $A_{r}^{i+1}$ and $Q_{r}^{i+1}$

Approximations $(2.3)$–$(2.6)$ are subjects to appropriate initial and boundary conditions.

Proof. Since $W = (A, Q)^{T}$, combining Eqs $(2.1)$ and $(2.2)$ to get

Substituting Eq $(1.3)$ into relation $(1.2)$ yields

In a similar way, substituting $(1.2)$ and $(2.11)$ into Eqs $(2.8)$ and $(2.10)$, respectively, results in

and

A combination of Eqs $(2.7), \; (2.9), \; (2.12)$ and $(2.13)$ provides a full description of the developed two-step predictor-corrector scheme, that is,

Here, the terms $A^{\overline{i+1}}$ and $Q^{\overline{i+1}}$ are "predicted" values of $A$ and $Q,$ respectively, at the time level $i+1.$ Assuming that the superscript $\overline{i+1}$ is a time-level, it is easy to see that the proposed MacCormack algorithm is a three-level method, so the initial conditions $A^{0}$ and $Q^{0}$ are needed to begin the algorithm. However, appropriate initial and boundary conditions must be specified. Further, the presence a cross-section in the denominator of several terms disallows zero cross sections. Therefore, a minimum cross section should be assigned to each node that is ponded. It is primarily the term $A^{10/3}$ in the denominator of the friction slope term given by relation $(2.11)$ that limits the magnitude of the minimum cross section and discharge. When the cross sections are very small, the friction slope is very large compared with the other terms in the second equation of system $(1.1)$. As cross sections increase rapidly during the early stages of flow development, the friction slope term magnitude changes much faster than the other terms. This phenomenon renders the second equation in system $(1.1)$ stiff and severely limits the maximum admissible time step for stability. Indeed, this phenomenon likely forced many researchers to consider the smallest time-steps relative to their mesh size (courant number $\ll1$) and keep lateral inflows and initial cross sections large ^[25].

3.   Stability analysis of the proposed numerical scheme

This section deals with a linear stability analysis of a three-level MacCormack procedure in a numerical solution of the system of nonlinear PDEs $(1.1)$. In the literature, a wide set of multi-level schemes for solving shallow water equations and unsteady transport problems have been used to advance the solution in time. The most popular of these techniques is the original explicit MacCormack approach, which was widely applied to obtain a time-accurate solution for fluid flow and aeroacoustic problems. Compared to a broad range of explicit methods such as, Lax-Wendroff formulation, the natural MacCormack scheme is computationally efficient, easy to implement and it provides less computing time. In fact, this algorithm lies in a class of higher order finite difference methods for solving unsteady flow specifically in the presence of discontinuity, inherent dissipative, dispersion and stability ^[10], whereas the Lax-Wendroff method indicates some difficulties to preserve the stability and the temporal second-order accuracy at the same time ^[23]. However, for linear problems both numerical schemes have provided almost the same amplification factor and stability restriction (for instance, see ^[24], P: $227$–$228$). To analyze the stability of the proposed approach, we apply the Fourier stability method to the difference Eqs $(2.3)$–$(2.6)$, by computing the amplification factor to obtain an algebraic criterion. Following the Von Neumann stability analysis for necessary condition of stability, we assume that the approximate solution $W_{r}^{i} = (A_{r}^{i}, Q_{r}^{i})$ can be expressed in the form of Fourier series as

where $b = b_{1}+\widehat{i}b_{2}$, $a = a_{1}+\widehat{i}a_{2},$ with $a_{j}, b_{j}\in\mathbb{R}$, $1\leq j\leq2, \; \phi = k\Delta x$ is the phase angle, $k$ denotes the wave number, $\Delta x$ represents the mesh grid and $\widehat{i}$ indicates the imaginary unit. For the sake of readability, we assume in the following that $a_{2} = b_{2}$ and the ratios $|v||A|^{-1}$ and $|v||Q|^{-1}$ are very small.

The analysis of the stability of the proposed method requires the following remark.

Remark 3.1. Since the considered problem is nonlinear, the solution may contain discontinuity even if the initial conditions are smooth enough. To overcome this numerical challenge, we should assume that the phase angle $\phi = k\Delta x$ satisfies $|\phi| < < \pi.$ In fact, the method could be generally stabilized by adding additional dissipation to the scheme without affecting the order of accuracy.

The following Lemma gives the "temporary" linear stability constraint of the two-step predictor-corrector algorithm described in section $2.$

Lemma 3.1. A necessary condition of stability for the numerical scheme $(2.5)$ is given by

where $0\leq \epsilon < 1$, $\mu = \underset{0\leq i\leq K_{1}}{\max}\mu^{i}$, with $\mu^{i} = \frac{|Q^{i}|}{|A^{i}|}$, $\widehat{A} = \underset{0\leq i\leq K_{1}}{\min}|A^{i}| = \underset{0\leq i\leq K_{1}}{\min}|e^{a_{1}t^{i}}|\neq0$ (according to relation $(3.1)$), for any $a_{1}\in \mathbb{R}$, $\Gamma_{0} = \frac{gn_{1}^{2}}{1.49^{2}}P^{\frac{4}{3}}$ and $\phi = k\Delta x,$ where $k\neq0$ is the wave number.

Proof. Firstly, we should provide a simple expression of $A^{i+1}.$ Combining relations $(2.3)$ and $(2.4),$ simple computations give

and

Substituting Eqs $(3.3)$ and $(3.4)$ into $(2.5)$, this results in

Neglecting the last term in $(3.5),$ we obtain

where

Indeed, since the terms $|v(x, t)||A(x, t)|^{-1}$ and $|v(x, t)||Q(x, t)|^{-1}, \; \forall \; (x, t)\in[0, L]\times[0, T_{1}],$ are too small, the tracking of the last term in $(3.5)$ does not compromise the result.

Substituting Eq $(3.1)$ into relation $(3.6)$ to get

Multiplying both sides of relation $(3.8)$ by $e^{-at^{i}}e^{-\widehat{i}kx_{r}}$ to obtain

Using the identities: $e^{\widehat{i}\phi}-e^{-\widehat{i}\phi} = 2\widehat{i}\sin(\phi), \; e^{2\widehat{i}\phi}-2+e^{-2\widehat{i}\phi} = -4\sin^{2}(\phi)$ and $1-e^{\frac{4}{3}\widehat{i}\phi} = 2\sin^{2}(\frac{2}{3}\phi)+\widehat{i}\sin(\frac{4}{3}\phi),$ Eq $(3.9)$ becomes

where

Setting

A combination of Eqs $(3.12)$ and $(3.10)$ yields

Of course the aim of this report is to give a general picture of the necessary condition of stability. Since the formulae can become quite heavy, for the convenient of writing, the proof of our results will consider Remark $3.1$. However, this leads to a linear stability condition which can be observed as a necessary condition of stability.

We start the analysis with some extreme cases: $|\phi| = \pi$ and $\phi = 0$.

$\bullet$ Case $1.$ For $|\phi| = \pi,$ it comes from Eq $(3.13)$ that the amplification factor becomes

The squared modulus of the amplification factor gives

But there exist values of $\alpha_{2}$ for which $\sin\alpha_{2} = 1.$ So $|e^{a\Delta t}|^{2} > 1.$ Thus, the scheme is unconditionally unstable.

$\bullet$ Case $\phi = 0.$ In that case, the amplification factor given by Eq $(3.13)$ provides $e^{a\Delta t} = 1$ which implies $|e^{a\Delta t}| = 1.$ Thus, the numerical scheme is stable. This suggests that the proposed technique is not dissipative in the sense of Kreiss ^[44] and when applied to the system of nonlinear Eq $(1.1).$ That is, the computations should become unstable in certain circumstances. This instability is entirely due to the non-linearity of the equations, since the same scheme applied to linear shallow water equations without source terms does not diverge, although strong oscillations are generated (for example, see ^[10,17,45]).

$\bullet$ Case where $0 < |\phi| < < \pi.$ Using the Taylor expansion around $\phi = 0,$ and neglecting high-order terms, the squared modulus of the amplification factor given by $(3.13)$ is approximated as

For $|e^{a\Delta t}|^{2}$ to be less than one, the quantity $\left(1-\frac{2\Delta t^{2}}{3\Delta x}|\gamma_{2}|\phi\sin\alpha_{2}\right)^{2}+ \left(\frac{\Delta t}{\Delta x}|\gamma_{1}|+\frac{2\Delta t^{2}}{3\Delta x}|\gamma_{2}|\cos\alpha_{2}\right)^{2}\phi^{2}$ must be less than one. This implies

for any value of $\alpha_{2}$, where $0\leq\epsilon < 1$. Using simple calculations, it is not hard to see that estimates given by $(3.15)$ are equivalent to

for every $\alpha_{2}\in\mathbb{R}$. But, the second inequality in $(3.16)$ implies $\frac{3}{2}(\epsilon-1)\leq\frac{\Delta t}{\Delta x}\left(\frac{3}{2}|\gamma_{1}|+\Delta t|\gamma_{2}|\right)|\phi|\leq \frac{3}{2}(1-\epsilon)$. Since $\frac{\Delta t^{2}}{\Delta x}|\gamma_{2}|\phi\sin\alpha_{2}\leq\frac{\Delta t^{2}}{\Delta x}|\gamma_{2}||\phi|$, these facts combined with the first estimate in $(3.16)$ show that the numerical scheme $(2.5)$ is stable if

Multiplying these inequalities side by side and rearranging terms to get

Since $|\mu^{i}| = \left|\frac{Q_{j}^{i}}{A_{j}^{i}}\right| = e^{(b_{1}-a_{1})t^{i}},$ it comes from relation $(3.12)$ that

But $\underset{0\leq i\leq K_{1}}{\max}|A^{i}|^{-1} = \left[\underset{0\leq i\leq K_{1}}{\min}|A^{i}|\right]^{-1}$, taking the maximum in both sides of estimates $(3.17)$, for $i = 0, 1, ..., K_{1},$ the proof of Lemma $2$ is completed thank to relations $(3.18)$ and equality $|A^{i}| = e^{a_{1}t^{i}}$.

Lemma 3.2. The numerical scheme $(2.6)$ is linearly stable if the following estimate holds:

for $i = 0, 1, ..., K_{1}$, where $\beta^{*}\in(0, 1)$ and

with

Estimate $(3.19)$ represents a necessary condition of stability for the numerical scheme $(2.6).$

Proof. Using the Taylor series expansion, the terms: $\frac{1}{2}\left(Q^{i}_{r}+Q_{r}^{\overline{i+1}}\right)$, $\frac{g}{4T}\left[(A_{r}^{\overline{i+1}})^{2}- (A_{r-1}^{\overline{i+1}})^{2}\right]$, $\frac{\left(Q_{r}^{\overline{i+1}}\right)^{2}}{A_{r}^{\overline{i+1}}}-\frac{\left(Q_{r-1}^{\overline{i+1}}\right)^{2}} {A_{r-1}^{\overline{i+1}}}$ and $\frac{Q_{r}^{\overline{i+1}}|Q_{r}^{\overline{i+1}}|}{\left(A_{r}^{\overline{i+1}}\right)^{\frac{7}{3}}}$ can be approximated as:

and

where

where

where the functions $H_{1}$, $H_{2}$, $K_{1}^{2}$, $K_{2}^{2}$ and $K_{1}K_{2}$ are given by

and

$C_{ls}, \; \overline{C}_{ls}, \; l, s = 1, 2,$ are given by $(3.25)$ and $\alpha_{2}$ and $\alpha_{3}$ follow from relation $(3.26).$ Furthermore

Where the functions $C_{1l}$ and $\widehat{C}_{1l}, \; l = 1, 2,$ are defined by relations $(3.25)$ and $(3.26), \; C_{3l}$ and $\overline{C}_{3l}$, $l = 1, 2,$ are given by

Combining Eqs $(3.23)$, $(3.24)$, $(3.27)$ and $(3.34)$, the amplification factor of the numerical scheme $(2.6)$ is approximated as

Where $\hat{i}$ is the unit complex number,

$H_{1}$, $H_{2}$, $K_{1}^{2}$, $K_{2}^{2}$ and $K_{1}K_{2}$ are given by Eqs $(3.28)$–$(3.32)$, respectively; $C_{ls}, \; \overline{C}_{ls}, \; s, l = 1, 2;$ come from $(3.25)$; $\alpha_{2}$ and $\alpha_{3}$ follow from Eq $(3.26)$ and $C_{3l}, \; \overline{C}_{3l}, \; l = 1, 2,$ are defined by relation $(3.35).$ The functions $\theta_{1}$ and $\theta_{2}$ are given implicitly by relations

and

Taking the squared modulus on both sides of approximation $(3.36)$ and performing straightforward computations to get estimates $(3.19)$–$(3.22)$ given in Lemma $3.2$.

Using the above results (namely, Lemmas $3.1$ and $3.2$) we are ready to give a necessary stability restriction of the two-step predictor-corrector method $(2.3)$–$(2.6)$ and to compare it with what is available in the literature (for example, Courant-Friedrich-Lewy condition for linear hyperbolic PDEs).

Theorem 3.1. The two-step predictor-corrector scheme $(2.3)$–$(2.6)$ applied to nonlinear shallow water equations with source terms $(1.4)$ is linearly stable if

where $|\phi| = |k\Delta x| < < \pi$, $\mu = \underset{0\leq i\leq K_{1}}{\max}\mu^{i}$, with $\mu^{i} = \frac{|Q^{i}|}{|A^{i}|}$, $\widehat{A} = \underset{0\leq i\leq K_{1}}{\min}|A^{i}| = \underset{0\leq i\leq K_{1}}{\min}|e^{a_{1}t^{i}}|\neq0$, for any $a_{1}\in \mathbb{R}$, $\Gamma_{0} = \frac{gn_{1}^{2}}{1.49^{2}}P^{\frac{4}{3}}$, $0 < \epsilon, \mathit{\text{}}\beta^{*} < 1, \; W_{1}(\Delta t, \Delta x)$ and $W_{2}(\Delta t, \Delta x)$ are given by relations $(3.20)$ and $(3.21)$, respectively.

Proof. The proof of Theorem $3.1$ is obvious according to Lemma $3.1$ and $3.2.$

The Von Neumann stability approach, based on a Fourier analysis in the space domain has been developed for nonlinear one-dimensional shallow water equations with source terms. Although the stability condition has not to be derived analytically, we have analyzed the properties of the amplification factor numerically (by use of Taylor series expansion), which contain information on the dispersion and diffusion errors of the considered numerical scheme. It is worth noticing that we used a local, linearized stability analysis to obtain an estimate $(3.38),$ which must be considered as a necessary condition of stability for the numerical scheme $(2.3)$–$(2.6).$

Some important remarks on stability analysis

This section considers some useful remarks based on the stability restrictions and compares them with what is known in the literature: CFL condition.

(i1) The linear stability condition $(3.38)$ suggests that a small space step $\Delta x$ cannot force the time step $\Delta t$ to be more potentially small (indeed: $|\phi|^{2} = k^{2}\Delta x^{2}$). This improves the convergence speed of the proposed numerical scheme. Moreover, because consistency requires that $\frac{\Delta t^{i}}{\Delta x} \; (i\geq1)$ approaches zero as $\Delta t$ and $\Delta x$ tend zero, an acceptable time step (not too small) should be applied to guarantee the stability condition $(3.38)$. For this reason, the developed three-level MacCormack method is suitable for the calculation of steady solutions (where time accuracy is unimportant) and the unsteady ones.

(i2) The developed approach $(2.3)$–$(2.6)$ for one-dimensional surface water equations has a linear stability limitation $(3.38)$ that limits the maximum time step. This stability requirement does not coincide with the CFL condition obtained for linear hyperbolic PDEs (for example: linear advection equation, wave equation, linearized burgers equations, etc...) because the considered algorithm is applied to complex time-dependent PDEs. As a discussion on the stability restrictions one can refer to the stability analysis of the two-step Lax-Wendroff method and the MacCormack scheme applied to complete burgers equations (for example, see ^[24], P. $245$–$247$). The linear stability condition $(3.38)$ is highly unusual. Since we normally find this condition from a Fourier stability analysis, it follows from inequality $(3.38)$ that instability occurs when $|\Delta t|$ is greater than $|\Delta t|_{\max} = (\Delta t)_{CFL}.$ However, it comes from the linear stability restriction $(3.38)$ that the empirical formula

can be used with an appropriate safety factor. It should be remembered that the "heuristic" stability analysis, i.e., estimates $(3.38)$ can only provide a necessary condition for stability. Thus, for some finite difference algorithms, only partial information about the complete stability bound is obtained and for others (such as algorithms for the heat equation, wave equation and linearized Burgers equations) a more complete theory must be employed.

4.   Numerical experiments and Convergence rate

This section simulates the two-step explicit predictor-corrector scheme $(2.3)$–$(2.6)$ for solving the one-dimensional shallow water Eq $(1.4)$. We consider two examples described in ^[46] to demonstrate the efficiency and robustness of the proposed technique. A practical application of a shallow water flow deals with the Benoué river. The river is located in Cameroon with a $7000m$ long reach of the upstream part (altitude = $174.22$ m). Furthermore, the characteristics of the flow consider a rather steep part in the first kilometers together with strong irregularities in the cross section and a low base discharge $(708m^{3}/s)$, altogether, produce a high velocity basic flow, transcritical in some parts. Specifically, the floods problem observed in this river in $2012$ is discussed since it is a classical example of unsteady nonlinear flow with shocks to expect floods and to test conservation in numerical schemes. In addition, the considered problem is assumed to be generated by a one-dimensional shallow water equation for the ideal case of a flat and frictionless channel with a prismatic cross-section wetted perimeter ($P = 366, 4$ m) and top width ($T = 348$ m). The initial data are provided by Eq $(4.1).$

Let $w\in L^{2}(0, T_{1};L^{2}(0, L))$, $k = \Delta t$ and $\tau = \Delta x$, we introduce the following discrete norm

Denote $W_{r}^{i} = (A_{r}^{i}, Q_{r}^{i})$, be the computed solution provided by the two-step approach $(2.3)$–$(2.6)$ and $w(t^{i}, x_{r}) = \left(A(t^{i}, x_{r}), Q(t^{i}, x_{r})\right)$, be the exact one at the discrete point $(t^{i}, x_{r})$, the error at time $t^{i}$ and position $x_{r}$ is given by $e(t^{i}, x_{r}) = w(t^{i}, x_{r})-W_{r}^{i} = \left(A(t^{i}, x_{r})-A_{r}^{i}, Q(t^{i}, x_{r})-Q_{r}^{i}\right)$. Thus,

$\bullet$ Problem 1 (Dam break on a dry domain with friction). The analytical solution is computed by Dressler's dam break with friction ^[46]. In the literature different approaches are deeply analyzed. Dressler's analyzed Chézy friction law using a perturbation method in the Ritter's scheme, i.e., both velocity: $u(m/s)$ and height: $h(m)$, of the water are expanded as power series in the friction coefficient $C_{f} = 1/C^{2}$. The initial conditions are defined as

We assume that $C = 40m^{1/2}/s$ (Chézy coefficient), $h_{l} = 5\times 10^{-3}$m, $x_{0} = L/2, \; T_{1} = 1$s, and $L = 1$m. Dressler's first order developments for the flow resistance give the following corrected height and velocity

where

and

Following Dressler's technique, four regions are considered: from upstream to downstream (a steady state region $(h_{l}, 10^{-1})$ for $x\leq x_{1}(t)$); a corrected region ($(h_{c}, u_{c})$ for $x_{1}(t)\leq x\leq x_{2}(t)$); the tip region (for $x_{2}(t)\leq x\leq x_{3}(t)$) and the dry region ($(0, 0)$ for $x_{3}(t)\leq x\leq L$). In the tip region, the friction term is preponderant thus $(4.2)$ is no more valid. The velocity increases in the corrected region with $x$, thus Dressler assumes that the velocity reaches the maximum of $u_{c}$ at $x_{2}(t)$ and it is constant in space in the tip region.

Armed with these assumptions together with Eqs $(4.2)$–$(4.5)$, the exact solution of the considered shallow water problem is defined as

and

where $\alpha_{1}$ and $\alpha_{2}$ are given by Eqs $(4.3)$ and $(4.4),$ respectively, $x_{1}(t) = x_{0}-t\sqrt{gh_{l}}, \; x_{3}(t) = x_{0}+2t\sqrt{gh_{l}}$ and $x_{2}(t)\in[x_{1}(t), x_{3}(t)]$ is the point where the velocity $u_{c}(t, x)$ attains its maximum.

$\bullet$ Problem 2 (Dam break on a dry domain without friction). We consider the exact solution introduced by Dressler's dam break with friction ^[46]. We also consider Ritter's solution corresponding to an ideal dam break (case of a reservoir with a constant height $h_{l}$) on a dry region. Similar to the Stoker's solution, the dam break is instantaneous, the bottom is flat and there is no friction. Using this fact, the initial condition (Riemann problem) is given by equation $(4.1)$. In addition, we assume that $C = 40m^{1/2}/s$ (Chézy coefficient), $h_{l} = 5\times10^{-3}$m, $x_{0} = L/2, \; T_{1} = 4$s, and $L = 8$m. At time $t > 0,$ the free surface is the constant water height $(h_{l})$ at rest connected to a dry zone $(h_{r})$ by a parabola. This parabola is limited upstream (resp. downstream) by the abscissa $x_{A}(t)$ (resp. $x_{B}(t)$). The analytical solution is defined as

and

where $x_{A}(t) = x_{0}-t\sqrt{gh_{l}}$ and $x_{B}(t) = x_{0}+2t\sqrt{gh_{l}}$.

The analysis considers the case where the channel is prismatic with a constant top width $(T)$ and average velocity $(u)$ defined as: $u(t, x) = Q(t, x)/A(t, x)$. Thus, the following formulas are satisfied

Since the water height is constant in the tip region, it follows from $(4.10)$ that the cross section $(A)$ and discharge $(Q)$ are not modified in that region. A combination of Eqs $(4.6)$–$(4.10)$ results in

$\bullet$ Problem 1 (Dam break on a dry domain with friction).

The table suggests that the proposed approach is second-order accurate in time and fourth-order convergent in space.

$\bullet$ Problem 2 (Dam break on a dry domain without friction).

and

The following values are used in the simulations: shear stress $\overline{\tau} = 1.329N/m^{2};$ Top width $T = 348m,$ wetter perimeter $P = 366, 4m;$ wavelength $K_{\lambda} = 2\pi\simeq6.28m,$ manning's number $n_{1} = 0.025s/m^{1/3},$ the acceleration of gravity $g = 10m/s^{2},$ the rainfall intensity is described as

We observe from this table that the proposed method is temporal second-order convergent and spatial fourth-order accurate.

The mathematical model for this ideal overland flow is the following: we consider a uniform plane catchment whose overall length in the direction of flow is $L(m)\in\{1, 8\}$. The surface roughness and shear stress are assumed invariant in space and time. It comes from Eq $(4.11)$ that the constant rainfall excess is defined as

The space step $\Delta x\in\{2^{-l}, \; l = 2, \cdots, 7\}$, while the time step $\Delta t$ varies in the range: $2^{-l}, \; l = 5, \cdots, 12. \; I$ is the rainfall intensity defined by relation $(4.11), \; t_{0} = 0s$ and $T_{1}(s)\in\{1, 4\}$, are initial and final times, respectively, and $L$ is the rod interval length. The numerical solutions given by Eqs $(2.3)$ and $(2.6)$ are displayed in Figures 1 and 2. Before $3$ iterations are encountered, the discharge wave propagates with almost a perfectly constant value at different positions (Figures 1 and 2). Furthermore, after these iterations, the discharge wave approaches zero at different times (Figures 1 and 2). So, the graphs suggest that the computed solution cannot grow with time and should satisfy the necessary condition $(3.38)$. Similar remarks are observed for the cross section. In addition, setting $k = \tau^{2}: = h_{1}^{2}$, Tables 1 and 2 indicate that the developed numerical technique is second-order accurate in time and spatial fourth-order. Furthermore, the graphs show that the numerical solutions start to destroy after a fixed time. Thus, physical insight must be used when the stability restriction $(3.38)$ of the constructed two-step predictor-corrector method $(2.3)$–$(2.6)$ is investigated. Finally, both Tables 1 and 2 and Figures 1 and 2 indicate that the approximate solutions do not increase with time and converge to the analytical one. Specifically, they show that stability for the proposed numerical scheme is subtle. It is not unconditionally unstable, but stability depends on the parameters $\Delta x$ and $\Delta t.$

5.   General conclusions and future works

In this paper, we have developed a two-step explicit predictor-corrector approach $(2.3)$–$(2.6)$ for solving the evolutionary nonlinear problem $(1.4)$. A necessary condition of stability of the proposed numerical scheme has been deeply analyzed using the Von Neumann stability approach whereas the convergence rate of the algorithm is numerically computed using the $L^{2}$-norm. The graphs (Figures 1 and 2) show that the considered method is both stable and convergent while Tables 1 and 2 suggest that the algorithm is temporal second-order convergent and fourth-order accurate in space. After a few number of iterations, the figures indicate that the numerical solutions strongly converge with the analytical ones. One should observe from Figures 1 and 2 that the only case where the exact solutions do not tend to zero corresponds to the initial condition. This follows from assumptions made by Dressler when constructing the analytical solutions. Our future investigations will develop a two-step explicit predictor-corrector scheme for the two-dimensional shallow water model.

Acknowledgment

The authors appreciate the valuable comments and remarks of anonymous referees which helped to greatly improve the quality of the paper.

Conflict of interest

The authors declare that they have no conflicts of interest.