1. In [1], which is well known in the context of Riemann solvers, Godunov introduced the concept of a monotone finite-difference scheme and showed that there are no monotone high-order accurate schemes among the linear two-layer-in-time ones. Further development of the theory of finite-difference shock-capturing schemes for hyperbolic systems of conservation laws was aimed, to a large degree, at overcoming this Godunov order barrier. As a result, various classes of difference schemes were developed in which a high order of accuracy for smooth solutions and monotonicity (in the case of a linear system and a scalar conservation law) are reached via nonlinear flux correction, which leads to the nonlinearity of these schemes even in the case of the linear transport equation. The basic classes of these schemes, which we called NFC (Nonlinear Flux Correction) schemes, include MUSCL [2], WENO [3], DG [4], and CABARET [5] schemes. The main advantage of these schemes is that they localize shock waves with high accuracy and do not generate considerable spurious oscillations.

It was shown that NFC schemes have at most the first order of both local convergence in regions of influence of shock waves [6, 7] and integral convergence on intervals with one of the boundaries lying in a region of shock wave influence [8–10]. At the same time, some high-order accurate nonmonotone schemes having analytic functions of numerical fluxes and, hence, approximating the Rankine–Hugoniot ε-conditions with higher accuracy preserve the high order of convergence in negative norm in integration over domains containing strong discontinuities [8]. As a result, these nonmonotone schemes, in contrast to NFC ones, preserve the high order of convergence in regions of shock influence despite the noticeable spurious oscillations on their fronts.

In this context, a method was proposed in [11] for constructing combined finite-difference shock-capturing schemes that combine the advantages of NFC and classical nonmonotone schemes. Namely, they localize shock fronts with high accuracy, while preserving the high order of convergence in regions of their influence. A combined finite-difference scheme makes use of a nonmonotone basic scheme having a higher order of convergence in regions of influence of shock waves. The basic scheme is used to construct a difference solution in the entire computational domain. In high-gradient regions, where this solution exhibits spurious oscillations, it is corrected by numerically solving internal initial-boundary value problems with the use of an NFC scheme. In [11] the Rusanov scheme [12] of third order of classical approximation was used as a basic and a monotone modification of CABARET [5] of second-order accuracy for smooth solutions was used as an internal NFC scheme. This modification of the CABARET scheme was studied in [10] and, in what follows, we refer to it as CABARETM.

A potential shortcoming of a combined scheme is that the oscillations arising at the shock front in the nonmonotone basic scheme can propagate over time into smooth parts of the computed exact solution (primarily, into the region of shock wave influence), so the computational domain for the internal NFC scheme gradually expands and the efficiency of the combined scheme degrades. However, an opposite situation actually takes place. It was revealed in [13] in the numerical solution of the classical Shu–Osher problem [14] by applying a DG-based NFC scheme [4] that the numerical solution does not converge locally to the exact solution behind the front of a shock wave propagating against the entropy perturbation background. This result is explained by arising numerical oscillations whose amplitude ceases to reduce with decreasing spatial step, starting with a sufficiently small value of the latter.

In this paper, we show that a similar difficulty arises in computing shock waves by applying other NFC schemes. Specifically, in the dam break problem for the shallow water equations, the numerical solution produced by the nonmonotone Rusanov scheme [12], despite the noticeable oscillations at the shock wave, converges monotonically to the exact solution with the second order in the region of shock influence. At the same time, the numerical solutions of this problem produced by the NFC schemes CABARETM [10] and WENO5 [3] exhibit undamped oscillations in the constant-flow domain between the shock wave and the centered rarefaction wave. As a result, these solutions do not converge locally in the region of shock wave influence.

2. In the case of a rectangular horizontal channel without bottom friction, the system of conservation laws of shallow water theory in the first approximation can be written in vector form as

$${{{\mathbf{u}}}_{t}} + {\mathbf{f}}{{({\mathbf{u}})}_{x}} = 0;$$
(1)

here,

$${\mathbf{u}} = \left( {\begin{array}{*{20}{c}} H \\ q \end{array}} \right),\quad {\mathbf{f}}({\mathbf{u}}) = \left( {\begin{array}{*{20}{c}} q \\ {{{q}^{2}}{\text{/}}H + g{{H}^{2}}{\text{/}}2} \end{array}} \right),$$
(2)

where H(x, t) is the fluid depth, \(q(x,t)\) is the flow rate of the fluid, and g = 9.81 is the acceleration of gravity. For system (1), (2), we consider the dam break problem, i.e., the Riemann problem with the following piecewise constant initial data:

$$H(x,0) = \left\{ \begin{gathered} 5,\quad x \leqslant 0 \hfill \\ 1,\quad x > 0, \hfill \\ \end{gathered} \right.\quad q(x,0) = 0.$$
(3)

The solution of this problem consists of a shock wave propagating at the constant speed D = 6.64 and a centered depression wave with a constant flow region in between. A numerical solution of problem (1)–(3) is constructed on a uniform rectangular grid \({{x}_{j}} = jh\), \({{t}_{n}} = n\tau (h)\) with the time step \(\tau (h)\) determined by the Courant stability condition

$$\tau (h) = \frac{{zh}}{{\mathop {\max }\limits_{x,t} \max \left( {{\text{|}}{{\lambda }_{ + }}({\mathbf{u}}(x,t)){\text{|}},\;{\text{|}}{{\lambda }_{ - }}({\mathbf{u}}(x,t)){\text{|}}} \right)}},$$
(4)

where \({{\lambda }_{ \pm }} = q{\text{/}}H \pm \sqrt {gH} \) are the velocities of the characteristics in system (1), (2); \({\mathbf{u}}(x,t)\) is the exact solution of problem (1)–(3); and z = 0.45 is the safety factor.

Figures 1 and 2 show the numerical results for problem (1)–(3) produced by the Rusanov [12], CABARETM [10], and WENO5 [3] schemes at the time T = 1. In Fig. 1a, the exact solution for the fluid depth is compared with the numerical solution obtained on a grid with the spatial step h = 0.36. It can be seen that the Rusanov nonmonotone scheme exhibits spurious oscillations at the shock front, whereas the NFC schemes CABARETM and WENO5 do not. Moreover, the shock wave and the weak discontinuities at the boundaries of the centered depression wave in CABARETM are smeared significantly less than in the Rusanov and WENO5 schemes.

Fig. 1.
figure 1

(a) Fluid depth and (b) local orders of convergence at the time T = 1 produced by the Rusanov (crosses), CABARETM (solid circles), and WENO5 (open circles) schemes. In panel (a), the solid curve depicts the exact solution, and the dashed line shows the initial value of the fluid depth.

Full size image
Fig. 2.
figure 2

Fluid depth in the constant-flow domain as produced by (a) CABARETM and (b) WENO5 schemes on grids with spatial steps \(h = 0.009\) (squares), \(h = 0.003\) (open circles), and \(h = 0.001\) (solid circles). The exact solution is shown by the horizontal line.

Full size image

Figure 1b presents the orders of local convergence of the difference solutions computed using the Runge formula

$${{\tilde {\rho }}_{j}} = \mathop {\log }\nolimits_{1/3} \frac{{{\text{|}}{{{\mathbf{v}}}_{{h/3}}}({{x}_{j}}(h),T) - {\mathbf{u}}({{x}_{j}}(h),T){\text{|}}}}{{{\text{|}}{{{\mathbf{v}}}_{h}}({{x}_{j}}(h),T) - {\mathbf{u}}({{x}_{j}}(h),T){\text{|}}}}$$
(5)

and corrected with the help of the limiter function

$${{\rho }_{j}} = \left\{ \begin{gathered} {{{\tilde {\rho }}}_{j}},\quad 0 \leqslant {{{\tilde {\rho }}}_{j}} \leqslant 2.5 \hfill \\ 0,\quad {{{\tilde {\rho }}}_{j}} < 0 \hfill \\ 2.5,\quad {{{\tilde {\rho }}}_{j}} > 2.5, \hfill \\ \end{gathered} \right.$$
(6)

where \({{x}_{j}}(h) = jh\), \(T = N\tau (h)\), N is a positive integer, and \({{{\mathbf{v}}}_{h}}\) and \({{{\mathbf{v}}}_{{h/3}}}\) are the numerical solutions obtained on grids with spatial steps h and h/3, respectively. The orders of convergence ρj were computed on the basis grid with \(h = 0.009\) and are shown in Fig. 1b for every 40th spatial grid node \(j = 40i\). Figure 1b shows that all three schemes have the first order of convergence within the centered depression wave. In the constant-flow domain between the shock and the depression wave, the Rusanov scheme has the second order of convergence, while the values of \({{\rho }_{j}}\) based on formulas (5), (6) for CABARETM and WENO5 strongly oscillate, so the order of local convergence of these schemes in the region of shock influence is uncertain.

To explain these results, we performed a series of test computations on a sequence of refined grids. It was found that the difference solution produced by the Rusanov scheme is monotone (with respect to both fluid depth and the flow rate) in the region of shock influence outside some neighborhoods of the shock front and the weak discontinuity at the right boundary of the depression wave, where the difference solution converges with the second order to the exact constant one. At the same time, the difference solutions based on the NFC schemes CABARETM and WENO5 exhibit undamped oscillations in the region of shock influence, and the oscillation structure depends on the value of the parameter z involved in stability condition (4). Figure 2a shows that, for z = 0.45 on the interval [3.2, 4.2], which lies within the region of shock influence, CABARETM exhibits numerical oscillations with roughly identical amplitudes for the spatial steps \({{h}_{1}} = h,\) \({{h}_{2}} = h{\text{/}}3,\) and \({{h}_{3}} = h{\text{/}}9\), where \(h = 0.009\), while the wave length is reduced roughly by a factor of three in the transition from hi to \({{h}_{{i + 1}}}\), i.e., it is proportional to \({{h}_{{i + 1}}}{\text{/}}{{h}_{i}}\). Figure 2b demonstrates a similar result for the numerical solutions produced by WENO5 on the interval [3.6, 4.6]. An analogous behavior of oscillations in the region of shock influence was obtained in [13] in the case of the DG method [4] applied to the Shu–Osher problem [14]. Figure 2 also shows that the amplitude of the oscillations obtained using CABARETM is about tenfold larger than in the case of WENO5, while the length of the oscillation waves is nearly identical for both schemes for a fixed spatial step hi.

Thus, the following general tendency is observed: the difference solutions produced by the NFC schemes may not exhibit local convergence to the exact solution in the regions of shock wave influence. In this case, in view of the Lax–Wendroff theorem [15], the limiting discontinuous solutions of the conservative NFC schemes are only weak solutions of the approximated system of conservation laws in the regions of shock wave influence.