ADER-WAF Schemes for the Homogeneous One-Dimensional Shallow Water Equations

Pavlos Stampolidis, Maria Ch. Gousidou-Koutita

Department of Mathematics, Faculty of Sciences, Aristotle University, Thessaloniki, Greece.

DOI: 10.4236/am.2025.161004                

Abstract

ADER-WAF methods were first introduced by researchers E.F. Toro and V.A. Titarev. The linear stability criterion for the model equation for the ADER-WAF schemes is $C_{\text{CFL}}\le 1$ , where $C_{\text{CFL}}$ denotes the Courant-Friedrichs-Lewy (CFL) coefficient. Toro and Titarev employed $C_{\text{CFL}}=0.95$ for their experiments. Nonetheless, we noted that the experiments conducted in this study with $C_{\text{CFL}}=0.95$ produced solutions exhibiting spurious oscillations, particularly in the high-order ADER-WAF schemes. The homogeneous one-dimensional (1D) non-linear Shallow Water Equations (SWEs) are the subject of these experiments, specifically the solution of the Riemann Problem (RP) associated with the SWEs. The investigation was conducted on four test problems to evaluate the ADER-WAF schemes of second, third, fourth, and fifth order of accuracy. Each test problem constitutes a RP characterized by different wave patterns in its solution. This research has two primary objectives. We begin by illustrating the procedure for implementing the ADER-WAF schemes for the SWEs, providing the required relations. Afterward, following comprehensive testing, we present the range for the CFL coefficient for each test that yields solutions with diminished or eliminated spurious oscillations.

Keywords

1D Shallow Water Equations, ADER-WAF Schemes, Finite Volume Schemes, Riemann Problem

1. Introduction

The SWEs, often known as the Saint-Venant system, are widely utilized by researchers to model and simulate various physical phenomena related to water flows with a free surface influenced by gravity. The SWEs are employed in numerous applications, including open channel and river flows, storm surges, dam-break waves, tidal fluctuations, tsunami waves, atmospheric flows, and more [1] [2]. As stated in [3], the SWEs serve as an approximation to the full free surface problems and can be obtained from the Navier-Stokes equations. Further information concerning the SWEs is available in the textbooks [3]-[5] and the associated references.

The resulting SWEs constitute a time-dependent system of nonlinear Partial Differential Equations (PDEs). Usually, the governing equations are hyperbolic and fall within the category of Hyperbolic Conservation Laws (HCLs), which (in one space dimension) can be expressed as a system of hyperbolic PDEs in the following form

${\partial }_{t}U+{\partial }_{x}F\left(U\right)=0,$ (1)

where $U\left(x,t\right)={\left[u_1\left(x,t\right),u_2\left(x,t\right),\cdots ,u_n\left(x,t\right)\right]}^{\top }$ is called the vector of conserved variables and $F\left(U\right)={\left[f_1,f_2,\cdots ,f_n\right]}^{\top }$ is the vector of fluxes with $f_i=f_i\left(u_1,u_2,\cdots ,u_n\right)$ , $i=1,\cdots ,n$ .

HCLs describe a diverse range of phenomena across various fields of physics, including quantum mechanics, continuum mechanics, and gravitational physics. In addition to SWEs, notable examples include Maxwell’s equations of electromagnetism, Einstein’s equations of gravitation, Navier’s equations of elasticity, and the Euler equations of gas dynamics.

The exact solutions of HCLs exist only in a limited number of special cases, thereby necessitating the development and evaluation of numerical methods for their approximate solution. The HCLs admit solutions that contain discontinuities, even for smooth Initial Conditions (ICs), due to hyperbolicity and non-linearity [6] [7]. The primary difficulty of numerical methods is the resolution of problems in discontinuous regions.

In recent decades, extensive research has focused on designing and implementing computational methods capable of accurately approximating discontinuous solutions. The numerical methods range over finite difference methods, finite volume methods, finite element methods, discontinuous Galerkin methods, spectral methods, and others. A powerful class of numerical methods is referred to as high-resolution methods after Harten [8]. These methods provide numerical solutions that are free of spurious oscillations and have a high resolution of discontinuities while achieving second or higher order accuracy in the smooth regions of the solution. Numerous high-resolution methods have been devised, including Total Variation Diminishing (TVD) schemes [8] [9], Essentially Non-Oscillatory (ENO) schemes [10]-[14], Weighted Essentially Non-Oscillatory (WENO) schemes [15]-[17], Arbitrary accuracy DErivative Riemann problem (ADER) methods [18]-[25], and others.

The ADER family includes fully discrete, one-step methods with arbitrary spatial and temporal accuracy, designed for solving HCLs and hyperbolic equations with source terms. These methods are a generalization of the Godunov method [3] [4] [26] and operate within the framework of finite volumes and discontinuous Galerkin methods [27] [28]. The initial publication of ADER methods was in [18], which addressed linear problems on Cartesian meshes. Subsequently, in [19]-[24], they were expanded to encompass non-linear problems. In their paper [23], Toro and Titarev introduced a new version of the ADER approach for HCLs titled the flux expansion version of ADER, as opposed to its predecessor, the state expansion version of ADER [20] (for HCLs). Furthermore, they propose the substitution of the first-order upwind Godunov flux [3] [4] [26] with the second-order upwind TVD flux of the Weighted Average Flux (WAF) method [3] [4] [29] [30], resulting in the development of a TVD variant of ADER methods. The scheme is referred to as the ADER-WAF scheme. Additional advancements of ADER schemes are also documented in [25].

The present study employs the flux expansion version of the ADER-WAF scheme of Toro and Titarev [23], which operates within the framework of finite volume methods. To develop ADER-WAF schemes of arbitrary spatial and temporal accuracy $k$ , the authors in [23] formulate and solve (approximately) a Generalised Riemann Problem (GRP) with ICs comprising two polynomials of degree $k-1$ separated by a discontinuity. The initial phase of the methodology entails piecewise polynomial data reconstruction, utilising nonlinear ENO or WENO reconstruction procedures to reduce or eliminate spurious oscillations near discontinuities and steep gradients of the solution. Following the reconstruction step, a GRP is formed at each cell interface. In the second phase of the methodology for $k$ -th order accuracy, the challenge of addressing the GRP is simplified to solving a series of $k$ conventional RPs: one nonlinear RP for the ICs and $k-1$ linear RPs for the $\nu$ -th order spatial derivatives of the ICs, where $\nu =1,\cdots ,k-1$ . Further details are provided in Section 3.4. It is noted that for the $k$ -th order accurate scheme (in both time and space), the reconstruction polynomials must be of degree $k-1$ . The current paper employed the characteristic-wise finite volume WENO procedure [31] [32] for the reconstruction step, with the $k$ -th order ADER-WAF scheme referred to as ADERk-WAF. Unless otherwise specified, the term “ADERk-WAF” implies that $k$ is equal to 2, 3, 4, or 5.

Linear ADER methods utilized for the linear advection equation with constant coefficient have the stability criterion $C_{\text{CFL}}\le 1$ [18], where $C_{\text{CFL}}$ denotes the Courant-Friedrichs-Lewy (CFL) coefficient. This is also the case for the ADER-WAF schemes, as stated in [23]. Additional information concerning the analysis of ADER and ADER-WAF schemes is available in [33], which investigates the stability properties and truncation errors of the finite volume ADER and ADER-WAF schemes on structured meshes, applied to the linear advection equation with constant coefficients in one, two, and three spatial dimensions. On the other hand, the CFL coefficient is typically applied empirically to nonlinear systems. For the experiments conducted by the authors in [20] and [23], the CFL coefficient was set at 0.95 for the ADER and ADER-WAF schemes.

The present research investigates the numerical behaviour of the ADER-WAF schemes as applied to RPs for the 1D SWEs. We have noted that employing $C_{\text{CFL}}=0.95$ for the problems discussed in this study results in spurious oscillations in the numerical solution in certain instances. Regrettably, these oscillations persist despite the refinement of the mesh. We observed that decreasing the CFL coefficient diminishes or eliminates the artificial oscillations. The main objective of this study is to determine the range of the CFL coefficient for each problem that produces solutions with decreased or entirely eliminated oscillations. The second objective of this work is to outline the procedure for applying the ADER-WAF method for the SWEs while providing the required equations to facilitate its coding.

The manuscript is organized as follows. In Section 2, we present the mathematical elements required to describe the SWEs, as well as the terminology and basic notation that will be used throughout the paper. In Section 3, we describe the finite volume discretization and provide the ADERk-WAF procedure. In Section 4, we assess the numerical efficiency of the ADERk-WAF method employing two distinct CFL coefficient values for solving RPs of SWEs. Concluding remarks are given in Section 5. In an effort to guarantee the comprehensiveness of this investigation, we have introduced supplementary information in the Appendices, some of which, to our knowledge, are not available in literature.

Throughout this paper, vectors and matrices are represented using bold face letters, while scalars are represented using plain letters. All vector-related operations are carried out component-wise.

2. General Background and Notation

2.1. Governing Equations

The 1D SWEs are derived by supposing an incompressible fluid within a channel of unit width, characterized by negligible vertical velocity and nearly constant horizontal velocity throughout any cross-section of the channel. This assumption holds true when examining small-amplitude waves in a fluid that is comparatively shallow relative to its wavelength (see [5]). In differential form, the conservative formulation of the homogeneous 1D SWEs is

${\partial }_{t}U+{\partial }_{x}F\left(U\right)=0,$ (2)

where $U$ and $F\left(U\right)$ are the vectors of conserved variables and fluxes, given respectively by

$U=\left[\begin{array}{l}{u}_1\hfill \\ u_2\hfill \end{array}\right]=\left[\begin{array}{c}h\\ hu\end{array}\right],F\left(U\right)=\left[\begin{array}{c}{f}_1\left(u_1,u_2\right)\\ f_2\left(u_1,u_2\right)\end{array}\right]=\left[\begin{array}{c}{u}_2\\ \frac{u_2^2}{u_1}+\frac{1}{2}g{u}_1^2\end{array}\right]=\left[\begin{array}{c}hu\\ h{u}^2+\frac{1}{2}g{h}^2\end{array}\right],$ (3)

where $u_i=u_i\left(x,t\right)$ , $i=1,2$ , are conserved variables, while $h=h\left(x,t\right)$ and $u=u\left(x,t\right)$ are primitive variables. These variables are functions of space $x$ and time $t$ . Equations (2), (3) express the physical laws of conservation of mass and momentum, with $h$ representing water depth, $u$ indicating the fluid’s horizontal velocity, and $g$ signifying the gravitational constant. In our computations, the gravitational constant is set as 9.81 m/s².

2.2. Eigen Structure of the System

A different formulation for the SWEs (2)-(3) in quasi-linear form is possible for a smooth solution

${\partial }_{t}U+A\left(U\right){\partial }_{x}U=0,$

where the coefficient matrix $A\left(U\right)$ is the Jacobian matrix

$\begin{array}{c}A\left(U\right)=\frac{\partial F}{\partial U}=\left[\begin{array}{cc}\partial f_1/\partial u_1& \partial f_1/\partial u_2\\ \partial f_2/\partial u_1& \partial f_2/\partial u_2\end{array}\right]=\left[\begin{array}{cc}0& 1\\ -{\left(\frac{u_2}{u_1}\right)}^2+g{u}_1& \frac{2u_2}{u_1}\end{array}\right]\\ =\left[\begin{array}{cc}0& 1\\ -u^2+gh& 2u\end{array}\right]=\left[\begin{array}{cc}0& 1\\ -u^2+c^2& 2u\end{array}\right],\end{array}$ (4)

where $c=\sqrt{gh}$ represents the celerity. The eigenvalues of $A\left(U\right)$ are

${\lambda }_1=\frac{u_2}{u_1}-\sqrt{g{u}_1}=u-c,{\lambda }_2=\frac{u_2}{u_1}+\sqrt{g{u}_1}=u+c,$ (5)

with the corresponding right eigenvectors

$R^{\left(1\right)}=\left[\begin{array}{c}1\\ {\lambda }_1\end{array}\right]=\left[\begin{array}{c}1\\ u-c\end{array}\right],R^{\left(2\right)}=\left[\begin{array}{c}1\\ {\lambda }_2\end{array}\right]=\left[\begin{array}{c}1\\ u+c\end{array}\right].$ (6)

The matrix with columns that are eigenvectors (6) is represented as

$R=R\left(u,c\right)=\left[\begin{array}{cc}1& 1\\ u-c& u+c\end{array}\right].$ (7)

Then, the rows of the subsequent matrix

$L=R^{-1}\left(u,c\right)=\left[\begin{array}{cc}\frac{u+c}{2c}& \frac{-1}{2c}\\ \frac{-u+c}{2c}& \frac{1}{2c}\end{array}\right],$ (8)

are left eigenvectors of $A\left(U\right)$ .

The SWEs (2)-(3) are hyperbolic, and when the celerity $c$ stays positive, the equations are strictly hyperbolic [4].

2.3. The Riemann Problem for the SWEs

The following notations agree with the ones found in [3]. The RP for (2)-(3) is the initial value problem for (2)-(3) with ICs

$U\left(x,0\right)=\{\begin{array}{ll}{U}_L={\left[u_{1L},u_{2L}\right]}^{\top }\hfill & \text{if}x<x_0,\hfill \\ U_R={\left[u_{1R},u_{2R}\right]}^{\top }\hfill & \text{if}x>x_0,\hfill \end{array}$ (9)

and is indicated by the notation RP ($U_L$ , $U_R$ ). The subscripts L and R denote the left and right states, respectively, of the piecewise constant data exhibiting a singular jump discontinuity at a specific point, namely $x=x_0$ . The initial value problem (2), (3), (9) can be solved exactly.

The solution is a similar solution, meaning it is a function only dependent on the ratio $x/t$ . Two waves, each associated with an eigenvalue ${\lambda }_i$ for $i=1,2$ , divide the solution into three constant states. The waves emanate from the initial discontinuity with constant velocities. The solution to the left of the ${\lambda }_1$ -wave corresponds to the initial data $U_L$ , while the solution to the right of the ${\lambda }_2$ -wave corresponds to $U_R$ . The area enclosed by the ${\lambda }_1$ and ${\lambda }_2$ waves is usually referred to as the Star Region, with $U^{\text{*}}$ symbolising the solution inside that domain.

Two distinct kinds of waves exist: shock waves and rarefaction waves. From this point forward, shock waves and rarefaction waves will be termed as shocks and rarefactions, respectively. There are four possible wave formations which are depicted in Figure 1. At any given time moment, the solution varies continuously with $x$ through rarefactions, while it exhibits discontinuous jumps through shocks.

When solving the RP (2), (3), (9), the velocities of the head and tail of a left or right rarefaction will be denoted as $S_{HL}$ , $S_{TL}$ , $S_{HR}$ , and $S_{TR}$ , correspondingly. Whereas the velocities of a left or right shock will be designated as $S_{SL}$ and $S_{SR}$ , respectively. The formulas required to calculate these velocities are provided in

Figure 1. Potential wave patterns in the solution of the RP for the 1D SWE: (a) left rarefaction, right shock (b) left shock, right rarefaction (c) left rarefaction, right rarefaction (d) left shock, right shock.

[4]. The solution within the left or right rarefaction, represented as $U_{Lfan}$ and $U_{Rfan}$ , respectively, when rarefactions are present, can be obtained by

$U_{Lfan}\left(x,t\right)={\left[\frac{{\left(2c_L+u_L-S\right)}^2}{9g},\frac{\left(u_L+2S+2c_L\right){\left(2c_L+u_L-S\right)}^2}{27g}\right]}^{\top },$ (10)

$U_{Rfan}\left(x,t\right)={\left[\frac{{\left(2c_R-u_R+S\right)}^2}{9g},\frac{\left(u_R+2S-2c_R\right){\left(2c_R-u_R+S\right)}^2}{27g}\right]}^{\top },$ (11)

where $S=\left(x-x_0\right)/t$ , $u_K=u_{2K}/u_{1K}$ , $c_K=\sqrt{g{u}_{1K}}$ , $K=L,R$ . Refer to [4]-[6] for detailed information concerning the RP solutions. The exact solution to the RP for SWEs can be found in [4] [5].

3. ADER-WAF Scheme

Three crucial subjects must be covered before outlining the ADER-WAF scheme: the formulation of finite volume methods for HCLs, the determination of the time step $\text{Δ}t$ , and the reconstruction procedure.

3.1. The Finite Volume Framework

Consider the 1D HCLs (1). The computational domain is discretized into grid cells, or finite volumes $I_i=\left[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}\right]$ . In the $x-t$ plane, assume a control volume $V=I_i\times \left[t^n,t^{n+1}\right]$ with dimensions $\text{Δ}x=x_{i+\frac{1}{2}}-x_{i-\frac{1}{2}}$ and $\text{Δ}t=t^{n+1}-t^n$ .

In essence, two options exist for discretizing (1). One approach is to formulate entirely discrete one-step schemes, while the alternative is to maintain the continuity of the time variable $t$ and contemplate semi-discrete schemes. The former approach is employed to develop the ADER schemes.

(1) is integrated over volume $V$ with respect to $x$ and $t$ , producing the exact relation

${\overline{U}}_i^{n+1}={\overline{U}}_i^n-\frac{\text{Δ}t}{\text{Δ}x}\left({\overline{F}}_{i+\frac{1}{2}}-{\overline{F}}_{i-\frac{1}{2}}\right),$

where

${\overline{U}}_i^n=\frac{1}{\text{Δ}x}{\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}U\left(x,t^n\right)\text{d}x}$(12)

represents the spatial-integral average of the solution $U$ in the cell $I_i$ at time $t^n$ , commonly known as the cell average, and

${\overline{F}}_{i+\frac{1}{2}}=\frac{1}{\text{Δ}t}{\displaystyle {\int }_{t^n}^{t^{n+1}}F}\left(U\left(x_{i+\frac{1}{2}},t\right)\right)\text{d}t$(13)

denotes the time-integral average of the physical flux at the cell interface $x=x_{i+\frac{1}{2}}$ . By providing specific approximations $U_i^n$ and ${\widehat{f}}_{i+\frac{1}{2}}$ to the integrals ${\overline{U}}_i^n$ and ${\overline{F}}_{i+\frac{1}{2}}$ , respectively, in (12) and (13), one can construct a specific finite volume method. The conservative finite volume numerical scheme is expressed as

$U_i^{n+1}=U_i^n-\frac{\text{Δ}t}{\text{Δ}x}\left({\widehat{f}}_{i+\frac{1}{2}}-{\widehat{f}}_{i-\frac{1}{2}}\right),$(14)

where ${\widehat{f}}_{i+\frac{1}{2}}$ is referred to as the numerical flux. The numerical flux at the cell interface $x_{i+\frac{1}{2}}$ is defined as a monotone function of two values,

${\widehat{f}}_{i+\frac{1}{2}}=H\left(U_{i+\frac{1}{2}}^{-},U_{i+\frac{1}{2}}^{+}\right)$ (15)

which are referred to as the left and right boundary extrapolated values denoted as $U_{i+\frac{1}{2}}^{-}$ and $U_{i+\frac{1}{2}}^{+}$ , respectively. These extrapolated values $U_{i+\frac{1}{2}}^{-}$ and $U_{i+\frac{1}{2}}^{+}$ are approximations to $U\left(x_{i+\frac{1}{2}},t\right)$ from the left and right, respectively. They are derived from cell averages via a high-order polynomial reconstruction that is essentially non-oscillatory.

The evaluation of the cell averages $U_i^0$ at the initial time $t^0=0$ should be conducted using the exact integration (12) of the suitable ICs, that is

$U_i^0={\overline{U}}_i^0=\frac{1}{\text{Δ}x}{\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}U}\left(x,0\right)\text{d}x,$(16)

where $U\left(x,0\right)$ are the ICs.

3.2. Time Step

The following condition is implemented for the calculation of the time step $\text{Δ}t$

$\text{Δ}t=C_{\text{CFL}}\frac{\text{Δ}x}{S_{\max }^n},$ (17)

where $C_{\text{CFL}}\in \left(0,1\right]$ , $\text{Δ}x$ denotes the spatial step and $S_{\max }^n$ indicates the maximum wave velocity present in the domain at time level $t=t^n$ . The calculation of $S_{\max }^n$ in the present research is carried out using the formula

$S_{\max }^n=\underset{i}{\text{max}}\left\{\left|S_{i+\frac{1}{2}}^L\right|,\left|S_{i+\frac{1}{2}}^R\right|\right\},$ (18)

where $S_{i+\frac{1}{2}}^L$ and $S_{i+\frac{1}{2}}^R$ represent the wave velocities of the left and right nonlinear waves, respectively, in the solution of the RP($U_i^n,U_{i+1}^n$ ), as determined by a Riemann solver. The velocity of the head is selected for rarefactions, whereas the shock velocity is selected for shocks. Consult [3] [4] for additional information on the computation of time steps and Riemann solvers.

3.3. The Reconstruction Procedure

For an in-depth comprehension of the subsequent subjects, refer to [31] [32]. Based on the cell averages $\left\{{\overline{u}}_j\right\}$ of a piecewise smooth function $u\left(x\right)$ , specify $k$ ($k=1,2,\cdots$ ) reconstruction polynomials $p_i^{\left(k,s\right)}\left(x\right)$ over the stencils $T_s\left(i\right)=\left\{I_{i-s},\cdots ,I_{i-s+k-1}\right\}$ , $s=0,1,\cdots ,k-1$ , of degree at most $k-1$ at each cell $I_i$ that fulfil the following criterion

$\frac{1}{\text{Δ}x}{\displaystyle {\int }_{I_j}{p}_i^{\left(k,s\right)}\left(\xi \right)\text{d}\xi }={\overline{u}}_j,j=i-s,\cdots ,i-s+k-1,s=0,1,\cdots ,k-1.$

Subsequently, a convex combination of the $k$ reconstruction polynomials $p_i^{\left(k,s\right)}\left(x\right)$ is utilized to obtain $\left(2k-1\right)$ -th order approximations of the extrapolated values $u_{i+\frac{1}{2}}^{\pm }$ as per

$u_{i+\frac{1}{2}}^{-}={\displaystyle \sum }_{s=0}^{k-1}{\omega }_s{p}_i^{\left(k,s\right)}\left(x_{i+\frac{1}{2}}\right),$

$u_{i-\frac{1}{2}}^{+}={\displaystyle \sum }_{s=0}^{k-1}{\tilde{\omega }}_s{p}_i^{\left(k,s\right)}\left(x_{i-\frac{1}{2}}\right)\Rightarrow u_{i+\frac{1}{2}}^{+}={\displaystyle \sum }_{s=0}^{k-1}{\tilde{\omega }}_s{p}_{i+1}^{\left(k,s\right)}\left(x_{i+\frac{1}{2}}\right),$

where the coefficients ${\omega }_s$ and ${\tilde{\omega }}_s$ are commonly known as nonlinear weights, with their computation detailed in Procedure 1.

Following that, put them into a numerical flux $\widehat{f}\left(u_{i+\frac{1}{2}}^{-},u_{i+\frac{1}{2}}^{+}\right)$ to construct a scheme that corresponds to (14) in the scalar case.

For illustration purposes, when $k=3$ , each cell $I_i$ provides three reconstruction quadratic polynomials, $p_i^{\left(3,0\right)}\left(x\right)$ , $p_i^{\left(3,1\right)}\left(x\right)$ , and $p_i^{\left(3,2\right)}\left(x\right)$ , over the stencils $T_0\left(i\right)=\left\{I_i,I_{i+1},I_{i+2}\right\}$ , $T_1\left(i\right)=\left\{I_{i-1},I_i,I_{i+1}\right\}$ , and $T_2\left(i\right)=\left\{I_{i-2},I_{i-1},I_i\right\}$ , respectively. These polynomials must meet the following criteria

$\frac{1}{\text{Δ}x}{\displaystyle {\int }_{I_j}{p}_i^{\left(3,s\right)}\left(\xi \right)\text{d}\xi }={\overline{u}}_j,j=i-s,\cdots ,i-s+2,s=0,1,2.$

The extrapolated values are then acquired by

$u_{i+\frac{1}{2}}^{-}={\displaystyle \sum }_{s=0}^2{\omega }_s{p}_i^{\left(3,s\right)}\left(x_{i+\frac{1}{2}}\right),u_{i+\frac{1}{2}}^{+}={\displaystyle \sum }_{s=0}^2{\tilde{\omega }}_s{p}_{i+1}^{\left(3,s\right)}\left(x_{i+\frac{1}{2}}\right).$

A key component of the WENO reconstruction procedure is the calculation of what are typically referred to as smoothness indicators. Represented as ${\beta }_s^k$ , these indicators are defined as [16]

${\beta }_s^k={\displaystyle \sum }_{l=1}^{k-1}{\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}\text{Δ}{x}^{2l-1}}{\left(\frac{{\text{d}}^l{p}_i^{\left(k,s\right)}\left(x\right)}{\text{d}{x}^l}\right)}^2\text{d}x,s=0,1,\cdots ,k-1.$ (19)

They are expressly located in [31] for $k=2,3$ , and in [17] for $k=4,5$ . Additionally, they can be found in Appendix A.

The reconstruction for systems can be carried out by applying the WENO reconstruction procedure to each distinct component of the vector $U$ . Nonetheless, spurious oscillations may arise from the component-wise WENO reconstruction. The reconstruction is performed utilising local characteristic variables to avoid this issue. The characteristic reconstruction procedure involves the initial transformation of conservative variables into characteristic variables, succeeded by the implementation of the WENO reconstruction technique on each component of these variables. The final values are obtained by reverting to conservative variables. This study employs reconstruction in local characteristic variables.

3.4. An Overview of the Scheme

This section provides a concise overview of the ADER-WAF scheme developed by Toro and Titarev [23]. Procedure 1 contains the comprehensive calculations. For a deeper understanding, refer to [23] and its accompanying references. In the following, we employ the convention $x^{\left(0\right)}=x$ to any quantity $x$ .

As previously stated, there are two variations of the ADER methods: the flux expansion and the state expansion. In this study, we utilized the flux expansion approach.

In [23], the authors formulate the Taylor series expansion in time of the physical flux at $x_{i+\frac{1}{2}}$ , specifically

$F\left(x_{i+\frac{1}{2}},\tau \right)=F\left(x_{i+\frac{1}{2}},0+\right)+{\displaystyle \sum }_{\nu =1}^{k-1}\frac{{\tau }^{\nu }}{\nu !}{\partial }_t^{\nu }F\left(x_{i+\frac{1}{2}},0+\right),$

where

${\partial }_t^{\nu }F\left(x,t\right):=\frac{{\partial }^{\nu }}{\partial t^{\nu }}F\left(x,t\right),0+:=\underset{\tau \to 0+}{\text{lim}}\tau ,\tau =t-t^n.$

The numerical flux is determined by the combination of the aforementioned equation and (13)

${\widehat{f}}_{i+\frac{1}{2}}=F_{i+\frac{1}{2}}^{\left(0\right)}+{\displaystyle \sum }_{\nu =1}^{k-1}\frac{\text{Δ}{t}^{\nu }}{\left(\nu +1\right)!}{F}_{t,i+\frac{1}{2}}^{\left(\nu \right)},$(20)

where $F_{i+\frac{1}{2}}^{\left(0\right)}$ and $F_{t,i+\frac{1}{2}}^{\left(\nu \right)}$ are approximations to $F\left(x_{i+\frac{1}{2}},0+\right)$ and ${\partial }_t^{\nu }F\left(x_{i+\frac{1}{2}},0+\right)$ , respectively. Therefore, to compute ${\widehat{f}}_{i+\frac{1}{2}}$ , it is required to evaluate $F_{i+\frac{1}{2}}^{\left(0\right)}$ and $F_{t,i+\frac{1}{2}}^{\left(\nu \right)}$ .

The calculation of the leading term $F_{i+\frac{1}{2}}^{\left(0\right)}$ involves the solution of the (nonlinear) RP

${\partial }_{t}U+{\partial }_{x}F\left(U\right)=0,$

$U\left(x,0\right)=\{\begin{array}{ll}{U}_L=p_i\left(x_{i+\frac{1}{2}}\right)\hfill & \text{if}x<x_{i+\frac{1}{2}},\hfill \\ U_R=p_{i+1}\left(x_{i+\frac{1}{2}}\right)\hfill & \text{if}x>x_{i+\frac{1}{2}}.\hfill \end{array}$ (21)

The values represented as $p_i\left(x_{i+\frac{1}{2}}\right)$ and $p_{i+1}\left(x_{i+\frac{1}{2}}\right)$ correspond to the values of the reconstruction polynomials $p_i\left(x\right)$ and $p_{i+1}\left(x\right)$ of degree $k-1$ at $x=x_{i+\frac{1}{2}}$ in cells $I_i$ and $I_{i+1}$ , respectively.

To obtain $F_{t,i+\frac{1}{2}}^{\left(\nu \right)}$ , $\nu =1,2,\cdots ,k-1$ , the following procedure is applied:

• Begin by solving the linearized RP

${\partial }_t\left({\partial }_x^{\nu }U\left(x,t\right)\right)+A_{i+\frac{1}{2}}\left(U_{i+\frac{1}{2}}^{\left(0\right)}\right){\partial }_x\left({\partial }_x^{\nu }U\left(x,t\right)\right)=0,$

${\partial }_x^{\nu }U\left(x,0\right)=\{\begin{array}{l}{U}_L^{\left(\nu \right)}={\text{d}}^{\nu }{p}_i\left(x_{i+\frac{1}{2}}\right)\text{if}x<x_{i+\frac{1}{2}},\hfill \\ U_R^{\left(\nu \right)}={\text{d}}^{\nu }{p}_{i+1}\left(x_{i+\frac{1}{2}}\right)\text{if}x>x_{i+\frac{1}{2}},\hfill \end{array}$ (22)

where $A_{i+\frac{1}{2}}\left(U_{i+\frac{1}{2}}^{\left(0\right)}\right)$ is the Jacobian matrix provided in (4), evaluated at $U=U_{i+\frac{1}{2}}^{\left(0\right)}$ , with the computation of $U_{i+\frac{1}{2}}^{\left(0\right)}$ contingent upon the solution of the RP (21). The values $U_L^{\left(\nu \right)}$ and $U_R^{\left(\nu \right)}$ are obtained by taking the $\nu$ -th derivative of the WENO reconstruction polynomials $p_i\left(x\right)$ and $p_{i+1}\left(x\right)$ in (21) with respect to $x$ , and then evaluating them at $x=x_{i+\frac{1}{2}}$ . The reconstruction polynomials $p_i\left(x\right)$ and all of their derivatives are subject to the same weights and smoothness indicators, as recommended in [20].

• Utilise the Cauchy-Kowalewski method, also referred to as the Lax-Wendroff procedure [34], to express all time derivatives ${\partial }_t^{\nu }U$ as functions of space derivatives ${\partial }_x^{\nu }U$ . See (33)-(42) in Appendix C.1 as a demonstration.

• Derive the time derivatives ${\partial }_t^{\nu }F\left(U\right)$ of the flux $F\left(U\right)$ with respect to the time derivatives ${\partial }_t^{\nu }U={\left[{\partial }_t^{\nu }{u}_1,{\partial }_t^{\nu }{u}_2\right]}^{\top }$ of the conserved variables $u_1$ , $u_2$ of the vector $U$ . Refer to (43)-(46) in Appendix C.2 for an illustration.

As stated previously, the implementation of the ADER-WAF scheme requires the solution to the linear RP. Hence, we provide the solution in the subsequent proposition.

Proposition 1. Consider the strictly hyperbolic system of PDEs with constant coefficients

${\partial }_{t}U+\overline{A}{\partial }_{x}U=0,U\left(x,0\right)=\{\begin{array}{l}{U}_L\text{if}x<x_0,\hfill \\ U_R\text{if}x>x_0,\hfill \end{array}\left(x,t\right)\in ℝ\times \left(0,+\infty \right),$(23)

where

$U=\left[\begin{array}{c}{u}_1\\ u_2\end{array}\right],\overline{A}=\left[\begin{array}{cc}0& 1\\ -{\left(\frac{{\overline{u}}_2}{{\overline{u}}_1}\right)}^2+g{\overline{u}}_1& \frac{2{\overline{u}}_2}{{\overline{u}}_1}\end{array}\right],U_K=\left[\begin{array}{c}{u}_{1K}\\ u_{2K}\end{array}\right],K=L,R.$(24)

The exact solution of the RP (23)-(24) at $\left(x,t\right)$ is

$U\left(x,t\right)=\{\begin{array}{l}{U}_L\text{if}{\overline{\lambda }}_1>\left(x-x_0\right)/t,\hfill \\ U^{\text{*}}\text{if}{\overline{\lambda }}_1<\left(x-x_0\right)/t<{\overline{\lambda }}_2,\hfill \\ U_R\text{if}{\overline{\lambda }}_2<\left(x-x_0\right)/t,\hfill \end{array}$

where ${\overline{\lambda }}_1$ , ${\overline{\lambda }}_2$ are the eigenvalues of $\overline{A}$ provided in (5) and

$U^{\text{*}}=\left[\begin{array}{c}{u}_1^{\text{*}}\\ u_2^{\text{*}}\end{array}\right]=\frac{1}{2\overline{c}}\left[\begin{array}{c}-{\overline{\lambda }}_1{u}_{1L}+{\overline{\lambda }}_2{u}_{1R}+u_{2L}-u_{2R}\\ {\overline{\lambda }}_2{u}_{2L}-{\overline{\lambda }}_1\left[{\overline{\lambda }}_2\left(u_{1L}-u_{1R}\right)+u_{2R}\right]\end{array}\right],\overline{c}=\sqrt{g{\overline{u}}_1}.$

Proof. For information regarding the solution’s structure, refer to Section 2.3. Given that the eigenvectors ${\overline{R}}^{\left(1\right)}$ and ${\overline{R}}^{\left(2\right)}$ (provided in (6)) are linearly independent, we can express the data $U_L$ and $U_R$ as linear combinations of the set ${\overline{R}}^{\left(1\right)}$ and ${\overline{R}}^{\left(2\right)}$ , as follows

$U_L=a_1{\overline{R}}^{\left(1\right)}+a_2{\overline{R}}^{\left(2\right)},U_R=b_1{\overline{R}}^{\left(1\right)}+b_2{\overline{R}}^{\left(2\right)},a_i,b_i\in ℝ,i=1,2.$(25)

The first equation in (25) yields

$\left[\begin{array}{c}{u}_{1L}\\ u_{2L}\end{array}\right]=a_1\left[\begin{array}{c}1\\ {\overline{\lambda }}_1\end{array}\right]+a_2\left[\begin{array}{c}1\\ {\overline{\lambda }}_2\end{array}\right].$

Solving for the unknown coefficients $a_1$ and $a_2$ , we get

$a_1=\frac{{\overline{\lambda }}_2{u}_{1L}-u_{2L}}{2\overline{c}},a_2=\frac{u_{2L}-{\overline{\lambda }}_1{u}_{1L}}{2\overline{c}}.$

Likewise, by employing the second equation in (25) and solving for the coefficients $b_1$ and $b_2$ , we acquire

$b_1=\frac{{\overline{\lambda }}_2{u}_{1R}-u_{2R}}{2\overline{c}},b_2=\frac{u_{2R}-{\overline{\lambda }}_1{u}_{1R}}{2\overline{c}}.$

It is established (see e.g. [3]) that the solution, represented as $U^{\text{*}}$ , within the star region bounded by the ${\overline{\lambda }}_1$ and ${\overline{\lambda }}_2$ waves, i.e. ${\overline{\lambda }}_1<\left(x-x_0\right)/t<{\overline{\lambda }}_2$ , is

$U^{\text{*}}=b_1{\overline{R}}^{\left(1\right)}+a_2{\overline{R}}^{\left(2\right)}.$

Thus, subsequent to a series of algebraic operations, we derive

$\begin{array}{c}{u}_1^{\text{*}}=\frac{1}{2\overline{c}}\left({\overline{\lambda }}_2{u}_{1R}-{\overline{\lambda }}_1{u}_{1L}+u_{2L}-u_{2R}\right)\\ =\frac{1}{2}\left(u_{1L}+u_{1R}\right)+\frac{1}{2\overline{c}}\left[\overline{u}\left(u_{1R}-u_{1L}\right)+u_{2L}-u_{2R}\right],\overline{u}=\frac{{\overline{u}}_2}{{\overline{u}}_1},\end{array}$

$\begin{array}{c}{u}_2^{\text{*}}=\frac{1}{2\overline{c}}\left[{\overline{\lambda }}_2{u}_{2L}-{\overline{\lambda }}_1\left({\overline{\lambda }}_2\left(u_{1L}-u_{1R}\right)+u_{2R}\right)\right]\\ =\frac{\overline{u}}{2\overline{c}}\left[\overline{u}\left(u_{1R}-u_{1L}\right)+u_{2L}-u_{2R}\right]+\frac{1}{2}\left[\overline{c}\left(u_{1L}-u_{1R}\right)+u_{2L}+u_{2R}\right].\end{array}$

The solution to the left of the ${\overline{\lambda }}_1$ -wave, namely when $\left(x-x_0\right)/t<{\overline{\lambda }}_1$ , is $U_L$ , whereas the solution to the right of the ${\overline{\lambda }}_2$ -wave, namely when $\left(x-x_0\right)/t>{\overline{\lambda }}_2$ , is $U_R$ .

3.5. The Procedure of the Scheme

The spatial domain $\left[0,L\right]$ is discretized into $M$ computing cells $I_i=\left[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}\right]$ with a dimension of $x_{i+\frac{1}{2}}-x_{i-\frac{1}{2}}=\text{Δ}x=\frac{L}{M}$ , with $i=1,\cdots ,M$ .

Procedure 1. ADERk-WAF Scheme, $k=2,3,4,5$ .

Step 1. Initial conditions.

1.1. For $i=1,\cdots ,M$ set the ICs using (16).

1.2. Use the cell averages ${\left\{U_i^n\right\}}_{i=1}^{i=M}$ , $n\ge 0$ , as input.

Step 2. Boundary Conditions (BCs).

BCs are implemented through the application of the periodicity condition. Whenever necessary, we positioned $k+1$ fictitious cells adjacent to each boundary. For the left boundary, the fictitious cells are denoted by $i=0,-1,-2,\cdots$ and for the right boundary they are denoted by $i=M+1,M+2,\cdots$ . Then for $j=0,1,\cdots ,k+1$ apply

$U_{0-j}^n=U_{1+j}^n\left(\text{Left BCs}\right)$

$U_{M+1+j}^n=U_{M-j}^n\left(\text{Right BCs}\right)$

Step 3. Computation of time step.

3.1. For $i=0,\cdots ,M$ transform the vector of conserved variables $U_i^n={\left[u_{1i}^n,u_{2i}^n\right]}^{\top }$ to the vector of primitive variables $U{P}_i^n={\left[h_i^n,u_i^n\right]}^{\top }={\left[u_{1i}^n,u_{2i}^n/u_{1i}^n\right]}^{\top }$ .

3.2. For $i=0,\cdots ,M$ solve RP ($U{P}_i^n,U{P}_{i+1}^n$ ), utilising an exact Riemann solver [4], to compute $S_{\max }^n$ using (18) and then calculate $\text{Δ}t$ according to (17).

Step 4. Computation of the numerical flux.

4.1. For $i=-2,\cdots ,M+2$ , do:

Transformation into characteristic variables

4.1.1. Utilise (7) and (8) to calculate $R=R\left(\tilde{u},\tilde{c}\right)$ and $L=L\left(\tilde{u},\tilde{c}\right)$ , where $\tilde{u}$ and $\tilde{c}$ denote the Roe average [35] between $U_i$ and $U_{i+1}$ which are defined as

$\tilde{u}=\frac{u_i\sqrt{h_i}+u_{i+1}\sqrt{h_{i+1}}}{\sqrt{h_i}+\sqrt{h_{i+1}}},\tilde{c}=\sqrt{\frac{1}{2}\left(c_i^2+c_{i+1}^2\right)},c_i=\sqrt{g{h}_i}.$

4.1.2. Transform into characteristic variables using

$V_j=L{U}_j,j=i-k+1,\cdots ,i+k.$

4.1.3. Estimation of the left extrapolated values $V_{i+\frac{1}{2}}^{\left(\nu \right)-}$ , $\nu =0,1,\cdots ,k-1$ .

4.1.3.1. Compute ${\beta }_s^k={\beta }_s^k\left(V,i\right)$ , $s=0,1,\cdots ,k-1$ , employing (29)-(32) for $k=2,3,4,5$ respectively.

4.1.3.2. Compute ${\alpha }_s^k=\frac{d_{k,s}}{{\left(ϵ+{\beta }_s^k\right)}^2}$ , $s=0,1,\cdots ,k-1$ .

The constant $ϵ$ is introduced to prevent the denominator from becoming zero. For the ADERk-WAF scheme, we take $ϵ={10}^{-24}$ [23]. The values of $d_{k,s}$ can be located in [31] for $k=2,3$ , and in [17] for $k=4,5$ . Table A6 of Appendix B contain them, as well. In the literature, the coefficients $d_{k,s}$ are commonly referred to as ideal, optimal, or linear weights.

4.1.3.3. Compute ${\omega }_s^k=\frac{{\alpha }_s^k}{{{\displaystyle \sum }}_{j=0}^{k-1}{\alpha }_j^k}$ , $s=0,1,\cdots ,k-1$ .

For $\nu =0,1,\cdots ,k-1$ , calculate

4.1.3.4. $p_i^{\left(\nu \right)\left(k,s\right)}\left(x_{i+\frac{1}{2}}\right)={\displaystyle {\sum }_{j=0}^{k-1}{c}_{k,s,j}^{\left(\nu \right)}{V}_{i-s+j}}$ , $s=0,1,\cdots ,k-1$ ,

4.1.3.5. $V_{i+\frac{1}{2}}^{\left(\nu \right)-}={\displaystyle {\sum }_{s=0}^{k-1}{\omega }_s^k{p}_i^{\left(\nu \right)\left(k,s\right)}}\left(x_{i+\frac{1}{2}}\right)$ .

The coefficients $c_{k,s,j}^{\left(0\right)}$ are provided in [31] and in Table A1 of Appendix B. Tables A2-A5 in Appendix B provide the coefficients $c_{k,s,j}^{\left(\nu \right)}$ , $\nu =1,2,3,4$ .

4.1.4. Estimation of the right extrapolated values $V_{i+\frac{1}{2}}^{\left(\nu \right)+}$ , $\nu =0,1,\cdots ,k-1$ .

4.1.4.1. Compute ${\beta }_s^k={\beta }_s^k\left(V,i+1\right)$ , $s=0,1,\cdots ,k-1$ , employing (29)-(32) for $k=2,3,4,5$ respectively.

4.1.4.2. Compute ${\alpha }_s^k=\frac{d_{k-1-s}}{{\left(ϵ+{\beta }_s^k\right)}^2}$ , $s=0,1,\cdots ,k-1$ .

4.1.4.3. Compute ${\omega }_s^k=\frac{{\alpha }_s^k}{{{\displaystyle \sum }}_{j=0}^{k-1}{\alpha }_j^k}$ , $s=0,1,\cdots ,k-1$ .

For $\nu =0,1,\cdots ,k-1$ , calculate

4.1.4.4. $p_{i+1}^{\left(\nu \right)\left(k,s\right)}\left(x_{i+\frac{1}{2}}\right)={\displaystyle {\sum }_{j=0}^{k-1}{c}_{k,s-1,j}^{\left(\nu \right)}{V}_{i+1-s+j}}$ , $s=0,1,\cdots ,k-1$ ,

4.1.4.5. $V_{i+\frac{1}{2}}^{\left(\nu \right)+}={\displaystyle {\sum }_{s=0}^{k-1}{\omega }_s^k{p}_{i+1}^{\left(\nu \right)\left(k,s\right)}}\left(x_{i+\frac{1}{2}}\right)$ .

Transformation back into conservative variables

4.1.5. Calculate $U_{i+\frac{1}{2}}^{\left(\nu \right)\pm }=R{V}_{i+\frac{1}{2}}^{\left(\nu \right)\pm }={\left[u_{1i+\frac{1}{2}}^{\left(\nu \right)\pm },u_{2i+\frac{1}{2}}^{\left(\nu \right)\pm }\right]}^{\top }$ , $\nu =0,1,\cdots ,k-1$ , and store the values $U_{i+\frac{1}{2}}^{\left(\nu \right)\pm }$ , $\nu =1,\cdots ,k-1$ .

Conversion to primitive variables

4.1.6. Compute $U{P}_{i+\frac{1}{2}}^{\left(0\right)\pm }={\left[u_{1i+\frac{1}{2}}^{\left(0\right)\pm },u_{2i+\frac{1}{2}}^{\left(0\right)\pm }/u_{1i+\frac{1}{2}}^{\left(0\right)\pm }\right]}^{\top }={\left[h_{i+\frac{1}{2}}^{\left(0\right)\pm },u_{i+\frac{1}{2}}^{\left(0\right)\pm }\right]}^{\top }$ .

Solution of the nonlinear RP (21)

4.1.7. Solve the RP ($U{P}_{i+\frac{1}{2}}^{\left(0\right)-},U{P}_{i+\frac{1}{2}}^{\left(0\right)+}$ ) utilising a Harten, Lax and van Leer (HLL) Riemann solver [4] [36] and store the subsequent values

$c_{1,i+\frac{1}{2}}^{\left(0\right)}=S_{i+\frac{1}{2}}^L\frac{\text{Δ}t}{\text{Δ}x},c_{2,i+\frac{1}{2}}^{\left(0\right)}=S_{i+\frac{1}{2}}^R\frac{\text{Δ}t}{\text{Δ}x},$

$\text{Δ}{h}_{1,i+\frac{1}{2}}^{\left(0\right)}={\widehat{h}}^{\left(0\right)\text{*}}-h_{i+\frac{1}{2}}^{\left(0\right)-},\text{Δ}{h}_{2,i+\frac{1}{2}}^{\left(0\right)}=h_{i+\frac{1}{2}}^{\left(0\right)+}-{\widehat{h}}^{\left(0\right)\text{*}},$

$\text{Δ}{F}_{1,i+\frac{1}{2}}={\widehat{F}}^{\text{*}}-F_{i+\frac{1}{2}}^{-},\text{Δ}{F}_{2,i+\frac{1}{2}}=F_{i+\frac{1}{2}}^{+}-{\widehat{F}}^{\text{*}},$

$S{F}_{i+\frac{1}{2}}=F_{i+\frac{1}{2}}^{-}+F_{i+\frac{1}{2}}^{+},S{U}_{i+\frac{1}{2}}^{\left(0\right)}=U_{i+\frac{1}{2}}^{\left(0\right)-}+U_{i+\frac{1}{2}}^{\left(0\right)+},$

$\text{Δ}{U}_{1,i+\frac{1}{2}}^{\left(0\right)}=U^{\left(0\right)\text{*}}-U_{i+\frac{1}{2}}^{\left(0\right)-},\text{Δ}{U}_{2,i+\frac{1}{2}}^{\left(0\right)}=U_{i+\frac{1}{2}}^{\left(0\right)+}-U^{\left(0\right)\text{*}},$

where $F_{i+\frac{1}{2}}^{-}=F\left(U_{i+\frac{1}{2}}^{\left(0\right)-}\right)$ , $F_{i+\frac{1}{2}}^{+}=F\left(U_{i+\frac{1}{2}}^{\left(0\right)+}\right)$ , and $U^{\left(0\right)\text{*}}={\left[{\widehat{h}}^{\left(0\right)\text{*}},{\widehat{h}}^{\left(0\right)\text{*}}{\widehat{u}}^{\left(0\right)\text{*}}\right]}^{\top }$. The values ${\widehat{h}}^{\left(0\right)\text{*}}$, ${\widehat{u}}^{\left(0\right)\text{*}}$, and ${\widehat{F}}^{\text{*}}$ are estimates of the exact solutions for $h^{\text{*}}$ , $u^{\text{*}}$ , and $F\left(U^{\text{*}}\right)$ in the star region, respectively. The two-rarefaction Riemann solver [4] is employed to calculate ${\widehat{h}}^{\left(0\right)\text{*}}$ and ${\widehat{u}}^{\left(0\right)\text{*}}$. ${\widehat{F}}^{\text{*}}$, $S_{i+\frac{1}{2}}^L$ , and $S_{i+\frac{1}{2}}^R$ are estimated from the HLL Riemann solver.

End do

Computation of $U_{i+\frac{1}{2}}^{\left(0\right)}$ in (22) and solution of the linear RP (22)

4.2. For $i=-1,\cdots ,M+1$ , do:

Computation of $U_{i+\frac{1}{2}}^{\left(0\right)}$ .

For $\nu =0$ , calculate and store

4.2.1. $r_{j,i+\frac{1}{2}}^{\left(\nu \right)}=\{\begin{array}{l}\frac{\text{Δ}{h}_{j,i-\frac{1}{2}}^{\left(\nu \right)}}{\text{Δ}{h}_{j,i+\frac{1}{2}}^{\left(\nu \right)}}\text{if}{c}_{j,i+\frac{1}{2}}^{\left(\nu \right)}>0,\hfill \\ \frac{\text{Δ}{h}_{j,i+\frac{3}{2}}^{\left(\nu \right)}}{\text{Δ}{h}_{j,i+\frac{1}{2}}^{\left(\nu \right)}}\text{if}{c}_{j,i+\frac{1}{2}}^{\left(\nu \right)}<0,\hfill \end{array}j=1,2.$ (26)

4.2.2. $U_{i+\frac{1}{2}}^{\left(\nu \right)}=\frac{1}{2}S{U}_{i+\frac{1}{2}}^{\left(\nu \right)}-\frac{1}{2}{\displaystyle \sum }_{j=1}^2\text{sign}\left(c_{j,i+\frac{1}{2}}^{\left(\nu \right)}\right){\varphi }_{i+\frac{1}{2}}\left(r_{j,i+\frac{1}{2}}^{\left(\nu \right)},c_{j,i+\frac{1}{2}}^{\left(\nu \right)}\right)\text{Δ}{U}_{j,i+\frac{1}{2}}^{\left(\nu \right)},$ (27)

where $\varphi =\varphi \left(r,c\right)$ is a WAF limiter function that corresponds to the well-known flux limiter SUPERBEE of Roe [37] and is defined as

$\varphi \left(r,c\right)=\{\begin{array}{ll}1\hfill & \text{if}r\le 0,\hfill \\ 1-2\left(1-\left|c\right|\right)r\hfill & \text{if}0<r<\frac{1}{2},\hfill \\ \left|c\right|\hfill & \text{if}\frac{1}{2}\le r\le 1,\hfill \\ 1-\left(1-\left|c\right|\right)r\hfill & \text{if}1<r<2,\hfill \\ 2\left|c\right|-1\hfill & \text{if}r\ge 2.\hfill \end{array}$ (28)

Solution of the linear RP (22)

4.2.3. Compute the value of $A_{i+\frac{1}{2}}=A_{i+\frac{1}{2}}\left(U_{i+\frac{1}{2}}^{\left(0\right)}\right)$ , which is the Jacobian matrix $A\left(U\right)$ provided in (4) and evaluated at $U=U_{i+\frac{1}{2}}^{\left(0\right)}$ , as obtained from (27).

4.2.4. Solve the linear ${\text{RP }}^{\nu }=\text{RP}\left(U_{i+\frac{1}{2}}^{\left(\nu \right)-},U_{i+\frac{1}{2}}^{\left(\nu \right)+}\right)$ for $\nu =1,\cdots ,k-1$ using Proposition 1 with

$\overline{A}=A_{i+\frac{1}{2}},U_L=U_{i+\frac{1}{2}}^{\left(\nu \right)-},U_R=U_{i+\frac{1}{2}}^{\left(\nu \right)+},$

$U=U^{\left(\nu \right)}={\left[u_1^{\left(\nu \right)},u_2^{\left(\nu \right)}\right]}^{\top },U^{\text{*}}=U^{\left(\nu \right)\text{*}}={\left[u_1^{\left(\nu \right)\text{*}},u_2^{\left(\nu \right)\text{*}}\right]}^{\top },$

and maintain the subsequent values for each $\nu$

$c_{j,i+\frac{1}{2}}^{\left(\nu \right)}={\lambda }_{j,i+\frac{1}{2}}\frac{\text{Δ}t}{\text{Δ}x},j=1,2,\text{Δ}{h}_{1,i+\frac{1}{2}}^{\left(\nu \right)}=h^{\left(\nu \right)\text{*}}-h_{i+\frac{1}{2}}^{\left(\nu \right)-},$

$\text{Δ}{h}_{2,i+\frac{1}{2}}^{\left(\nu \right)}=h_{i+\frac{1}{2}}^{\left(\nu \right)+}-h^{\left(\nu \right)\text{*}},S{U}_{i+\frac{1}{2}}^{\left(\nu \right)}=U_{i+\frac{1}{2}}^{\left(\nu \right)-}+U_{i+\frac{1}{2}}^{\left(\nu \right)+},$

$\text{Δ}{U}_{1,i+\frac{1}{2}}^{\left(\nu \right)}=U^{\left(\nu \right)\text{*}}-U_{i+\frac{1}{2}}^{\left(\nu \right)-},\text{Δ}{U}_{2,i+\frac{1}{2}}^{\left(\nu \right)}=U_{i+\frac{1}{2}}^{\left(\nu \right)+}-U^{\left(\nu \right)\text{*}},$

where ${\lambda }_{j,i+\frac{1}{2}}$ , $j=1,2$ are the eigenvalues of the matrix $A_{i+\frac{1}{2}}$ (the same for all $\nu$ ), $h_{i+\frac{1}{2}}^{\left(\nu \right)\pm }=u_{1i+\frac{1}{2}}^{\left(\nu \right)\pm }$ and $h^{\left(\nu \right)\text{*}}=u_1^{\left(\nu \right)\text{*}}$ . Based on linear theory, the eigenvalues ${\lambda }_{1,i+\frac{1}{2}}$ and ${\lambda }_{2,i+\frac{1}{2}}$ , respectively, are the speeds of the left and right linear waves present in the solution of the linear ${\text{RP }}^{\nu }$ (see, e.g., [3]).

End do

Computation of $F_{i+\frac{1}{2}}^{\left(0\right)}$ , $F_{t,i+\frac{1}{2}}^{\left(\nu \right)}$ , and ${\widehat{f}}_{i+\frac{1}{2}}$ in (20)

4.3. For $i=0,\cdots ,M$ , do:

Computation of $F_{i+\frac{1}{2}}^{\left(0\right)}$

4.3.1. Compute $r_{j,i+\frac{1}{2}}^{\left(\nu \right)}$ , $\nu =1,\cdots ,k-1$ using (26),

4.3.2. Compute $F_{i+\frac{1}{2}}^{\left(0\right)}=\frac{1}{2}S{F}_{i+\frac{1}{2}}-\frac{1}{2}{\displaystyle \sum }_{j=1}^2\text{sign}\left(c_{j,i+\frac{1}{2}}^{\left(0\right)}\right){\varphi }_{i+\frac{1}{2}}\left(r_{j,i+\frac{1}{2}}^{\left(0\right)},c_{j,i+\frac{1}{2}}^{\left(0\right)}\right)\text{Δ}{F}_{j,i+\frac{1}{2}}$ ,

Computation of $F_{t,i+\frac{1}{2}}^{\left(\nu \right)}$ , $\nu =1,\cdots ,k-1$

4.3.3. Compute $U_{i+\frac{1}{2}}^{\left(\nu \right)}$ , $\nu =1,\cdots ,k-1$ applying (27).

4.3.4. Perform the following operations in the specified sequence, depending on the value of $k$ . It should be mentioned that the subsequent quantity is determined by considering the preceding ones.

4.3.4.1. For $k=2$ , set $U_{i+\frac{1}{2}}^{\left(0\right)}={\left[u_1,u_2\right]}^{\top }$ and $U_{i+\frac{1}{2}}^{\left(1\right)}={\left[u_{1x},u_{2x}\right]}^{\top }$ , then produce ${\partial }_{t}U={\left[u_{1t},u_{2t}\right]}^{\top }$ and ${\partial }_{t}F$ using (33) and (43), respectively. Set ${\partial }_{t}F=F_{t,i+\frac{1}{2}}^{\left(1\right)}$ after conducting the calculations.

4.3.4.2. For $k=3$ , assign $U_{i+\frac{1}{2}}^{\left(2\right)}={\left[u_{1xx},u_{2xx}\right]}^{\top }$ . Then, apply (34), (35), and (44) to derive ${\partial }_x{\partial }_{t}U={\left[u_{1xt},u_{2xt}\right]}^{\top }$ , ${\partial }_t^{2}U={\left[u_{1tt},u_{2tt}\right]}^{\top }$ , and ${\partial }_t^{2}F$ , respectively. Set ${\partial }_t^{2}F=F_{t,i+\frac{1}{2}}^{\left(2\right)}$ following the calculations.

4.3.4.3. For $k=4$ and setting $U_{i+\frac{1}{2}}^{\left(3\right)}={\left[u_{1xxx},u_{2xxx}\right]}^{\top }$ , we can calculate ${\partial }_t{\partial }_x^{2}U={\left[u_{1xxt},u_{2xxt}\right]}^{\top }$ , ${\partial }_t^2{\partial }_{x}U={\left[u_{1xtt},u_{2xtt}\right]}^{\top }$ , ${\partial }_t^{3}U={\left[u_{1ttt},u_{2ttt}\right]}^{\top }$ , and ${\partial }_t^{3}F$ employing (36)-(38), and (45), respectively. Upon performing computations, assign ${\partial }_t^{3}F=F_{t,i+\frac{1}{2}}^{\left(3\right)}$ .

4.3.4.4. For $k=5$ , set $U_{i+\frac{1}{2}}^{\left(4\right)}={\left[u_{1xxxx},u_{2xxxx}\right]}^{\top }$ and subsequently determine ${\partial }_t{\partial }_x^{3}U={\left[u_{1xxxt},u_{2xxxt}\right]}^{\top }$ , ${\partial }_t^2{\partial }_x^{2}U={\left[u_{1xxtt},u_{2xxtt}\right]}^{\top }$ , ${\partial }_t^3{\partial }_{x}U={\left[u_{1xttt},u_{2xttt}\right]}^{\top }$ , ${\partial }_t^{4}U={\left[u_{1tttt},u_{2tttt}\right]}^{\top }$ , and ${\partial }_t^{4}F$ utilising (39)-(42), and (46), respectively. After completing the calculations, put ${\partial }_t^{4}F=F_{t,i+\frac{1}{2}}^{\left(4\right)}$ .

Computation of ${\widehat{f}}_{i+\frac{1}{2}}$

4.3.5. Finally, derive and store ${\widehat{f}}_{i+\frac{1}{2}}=F_{i+\frac{1}{2}}^{\left(0\right)}+{\displaystyle \sum }_{\nu =1}^{k-1}\frac{{\left(\text{Δ}t\right)}^{\nu }}{\left(\nu +1\right)!}{F}_{t,i+\frac{1}{2}}^{\left(\nu \right)}$.

End do

Step 5. Updating of solution and advancing to the next time level.

Updating of solution

5.1. For $i=1,\cdots ,M$ apply the conservative equation (14) to update the conserved variables and obtain $U_i^{n+1}$ .

Advancing to the next time level

5.2. Utilising ${\left\{U_i^{n+1}\right\}}_{i=1}^{i=M}$ as initial values, proceed to Step 1.2 to advance to the next time level.

4. Numerical Results

This section compares the numerical performance of the ADER-WAF scheme applied to test problems employing two different values of the CFL coefficient. We established one at 0.95, as referenced in [23], and determined the other value subsequent to conducting numerical experiments. The test problems are RPs for the 1D SWEs, and each one corresponds to a certain wave pattern, which is depicted in Figure 1.

Tests 1 and 2 were obtained from [4], while Test 3 came from [5]. To complete the four wave patterns, we developed Test 4. To normalise the given spatial domain in the original Tests 1, 2, and 3 to the interval [0, 1], we transformed the initial discontinuity’s location and the output time. For each test problem, the interval [0, 1] is discretised by using M = 100, 200, 400, and 800 cells. The exact solutions to the test problems are obtained by the use of an exact Riemann solver.

In order to quantify the magnitude of the error $E$ , we utilized the $L_1$ norm [5]

${\Vert E\Vert }_1=\text{Δ}x{\displaystyle \sum }_{i=1}^M{\Vert {\overline{U}}_i^n-U_i^n\Vert }_1.$

Here, ${\overline{U}}_i^n$ is the cell average of $U\left(x,t\right)$ over the interval $I_i$ at time $t^n$ , $U_i^n$ is an approximation to the ${\overline{U}}_i^n$ provided by the numerical method, and ${\Vert U\Vert }_1={\Vert {\left[u_1,u_2\right]}^{\top }\Vert }_1=\left|u_1\right|+\left|u_2\right|$ is the $L_1$ vector norm on ${ℝ}^2$ .

The CFL coefficient utilized for each ADERk-WAF scheme was determined from our numerical experiments. We implemented an investigation that employed the values of $C_{\text{CFL}}=0.04+0.01i$ , $i=1,\cdots ,91$ for each Test. The decision for the selection of the $C_{\text{CFL}}$ values was based on the $L_1$ norm error and numerical oscillation. To compare the performance of the ADERk-WAF scheme between the two CFL coefficient values, we provide for each test graphs for either the water depth or the velocity and tabulate the $L_1$ norm error of cell averages of the solution for both CFL coefficient values. We also employ the notation $a:b:c$ , where $a,b,c\in ℝ$ , indicating the range from $a$ to $c$ in increments of $b$ . As an illustration, the values $0.\text{21}+0.0\text{1}i$ , where $i=0,1,\cdots ,9$ , are implied by the notation 0.21:0.01:0.30.

4.1. Computer-Related Issues

It is essential to provide some observations concerning computer-related issues. The first is about the solution to the RP (2), (3), (9). In particular, in order to obtain $h\left(x,t\right)$ in the star region, it is necessary to solve an algebraic nonlinear equation. To the best of our knowledge, no closed-form solution exists for this equation. Thus, we employ the Newton-Raphson iteration method [38] to numerically solve the equation to a tolerance of $TOL={10}^{-6}$ .

The second pertains to the cessation of the algorithms’ procedures upon attaining the specified output time, denoted as $t_{out}$ . The algorithm stops whenever $\left|{\displaystyle {\sum }_{n=1}^m\text{Δ}{t}^{\left(n\right)}}-t_{out}\right|\le TOL$ for some $m$ , where $\text{Δ}{t}^{\left(n\right)}$ represents the time step at time level $t^n$ .

The third issue pertains to the computation of the ratios $r_{j,i+\frac{1}{2}}^{\left(\nu \right)}$ in (26), specifically when the denominator is of small magnitude. Following [23], we set

$\text{Δ}{h}_{j,i+l}^{\left(\nu \right)}=TOL\cdot \text{sign}\left(1,\text{Δ}{h}_{j,i+l}^{\left(\nu \right)}\right)\text{if}\left|\text{Δ}{h}_{j,i+l}^{\left(\nu \right)}\right|\le TOL,l=\pm \frac{1}{2},\frac{3}{2},j=1,2.$

4.2. Test 1: Left Rarefaction and Right Shock

Consider the RP (2), (3), (9) with $U_L={\left[1,2.5\right]}^{\top }$ , $U_R={\left[0.1,0\right]}^{\top }$ , $\left(x,t\right)\in \left[0,1\right]\times \left(0,+\infty \right)$ and $x_0=0.2$ . The numerical solution is computed at the output time of $t_{out}=0.14$ (the original test [4] had $x\in \left(0,50\right)$ , with $x_0=10$ and $t_{out}=7$ ). The initial data of this test generate a right shock and a left rarefaction. Figure 10, left, depicts the structure of the solution to the RP in the $x-t$ plane. The exact solution is

$U\left(x,t\right)=\{\begin{array}{ll}{U}_L,\hfill & \text{if}S\le S_{HL},\hfill \\ U_{Lfan}\left(x,t\right),\hfill & \text{if}{S}_{HL}<S<S_{TL},\hfill \\ U^{\text{*}},\hfill & \text{if}{S}_{TL}\le S<S_{SR},\hfill \\ U_R,\hfill & \text{if}{S}_{SR}<S,\hfill \end{array}S=\frac{x-0.2}{t},$

where $S_{HL}\approx -0.632092$ , $S_{TL}\approx 1.415611$ , $S_{SR}\approx 4.620578$ , $U^{\text{*}}\approx {\left[0.611638,2.364063\right]}^{\top }$ , and $U_{Lfan}$ is determined by (10).

Upon examination of the exact solution, it is obvious that, when $t_{out}<0.8/S_{SR}\approx 0.173$ , the solution has a rich structure consisting of four distinct states that can be defined by smooth regions, derivative discontinuities (the corners at the endpoints of the rarefaction), and a jump discontinuity (shock). The numerical results indicate that the main difficulty encountered by the ADERk-WAF scheme in this test is effectively handling the corner at the tail of the rarefaction. Spurious oscillations near the shock and within the constant area between the two waves (which resemble a sinusoidal wave in the numerical solution) provide a challenge; however, they can be mitigated by decreasing the CFL coefficient value.

Our numerical experiments, not included herein, demonstrated that employing the values 0.36:0.01:0.42 for the CFL coefficient in ADER2-WAF and ADER3-WAF, and 0.37:0.01:0.42 in ADER4-WAF and ADER5-WAF, effectively prevents numerical oscillation (to the eyes) in the aforementioned areas. These CFL coefficient values are valid for spatial step sizes of up to $\text{Δ}x=0.00125$ . Employing a larger or smaller spatial step size may expand or contract the range of the CFL coefficient values, respectively. We selected a uniform value of 0.37 for $C_{\text{CFL}}$ in ADERk-WAF, without implying that it is the optimal choice for each ADERk-WAF individually.

Figures 2-9 present the computed water depth $h\left(x,t\right)$ with $C_{\text{CFL}}=0.37$ against the solution obtained with $C_{\text{CFL}}=0.95$ . The calculated water depth is plotted at $t_{out}=0.14$ with $M=100$ (Figures 2-5) and $M=400$ (Figures 6-9). The $L_1$ norm error of cell averages for both CFL coefficient values is tabulated in Table 1.

Figures 2-9 demonstrate that the ADERk-WAF exhibits superior performance with $C_{\text{CFL}}=0.37$ , particularly in the constant region between the two waves (star region) and in proximity to the jump discontinuity. The shock is captured without numerical oscillation, unlike the case with $C_{\text{CFL}}=0.95$ , which exhibits overshoot accompanied by spurious oscillations. Additionally, observe the significant reduction and/or elimination of numerical oscillations in the star region that results from the application of $C_{\text{CFL}}=0.37$ . It is also worth noticing that for $M=400$ and $C_{\text{CFL}}=0.95$ , the ADER2-WAF avoids the initiation of oscillations in the star region (Figure 6(c)) that commence with ADER3-WAF (Figure 7(c)) and become more severe with ADER4-WAF and ADER5-WAF, ultimately resembling a sinusoidal wave (Figure 8(c) and Figure 9(c)).

The resolution of the rarefaction is slightly improved with $C_{\text{CFL}}=0.37$ ; however, it is also satisfactory with $C_{\text{CFL}}=0.95$ . The tail of the rarefaction exhibits a more accurate approximation at its corner when utilising $C_{\text{CFL}}=0.37$ , in contrast to $C_{\text{CFL}}=0.95$ , with the exception of ADER4-WAF (Figure 4(b) and Figure 8(b)). Regrettably, fluctuations accompanying the rarefaction’s tail are inevitable. Nevertheless, they are eventually successfully damped, unlike the case with $C_{\text{CFL}}=0.95$ .

Upon examining Table 1, it is evident that for $M=100$ , the $L_1$ error is greater with $C_{\text{CFL}}=0.37$ compared to $C_{\text{CFL}}=0.95$ , with the exception of ADER5-WAF. This discrepancy is puzzling, considering that Figures 2(b)-4(b) visually validate the superiority of the solution with $C_{\text{CFL}}=0.37$ over $C_{\text{CFL}}=0.95$ . This also holds true for $M=200$ , but not for $M=400$ and $M=800$ . The reason for this is that the shock is better resolved when $C_{\text{CFL}}=0.95$ is employed, with the exception of ADER5-WAF. For $M=100$ and $M=200$ , the error produced

Figure 2. Numerical solution of Test 1 on 100 cells by ADER2-WAF scheme and its partially enlarged view.

Figure 3. Numerical solution of Test 1 on 100 cells by ADER3-WAF scheme and its partially enlarged view.

Figure 4. Numerical solution of Test 1 on 100 cells by ADER4-WAF scheme and its partial enlarged view.

Figure 5. Numerical solution of Test 1 on 100 cells by ADER5-WAF scheme and its partially enlarged view.

Figure 6. Numerical solution of Test 1 on 400 cells by ADER2-WAF scheme and its partially enlarged views.

Figure 7. Numerical solution of Test 1 on 400 cells by ADER3-WAF scheme and its partially enlarged views.

Figure 8. Numerical solution of Test 1 on 400 cells by ADER4-WAF scheme and its partially enlarged views.

Figure 9. Numerical solution of Test 1 on 400 cells by ADER5-WAF scheme and its partially enlarged views.

Table 1. Test 1. Comparison study of the $L_1$ norm error of the ADERk-WAF schemes at $t=0.\text{14}$ for two distinct values of the CFL coefficient.

$M$	$\Delta x$	ADER2-WAF		ADER3-WAF		ADER4-WAF		ADER5-WAF
		$L_1$ norm error
		0.37	0.95	0.37	0.95	0.37	0.95	0.37	0.95
100	0.01	4.4722E−03	3.9278E−03	3.9204E−03	3.1390E−03	3.7874E−03	3.4525E−03	3.3348E−03	1.7687E−02
200	0.005	2.7008E−03	2.5659E−03	2.3161E−03	1.8487E−03	2.0468E−03	2.2699E−03	1.9981E−03	1.7860E−02
400	0.0025	1.2829E−03	1.4190E−03	1.2017E−03	1.3926E−03	1.1567E−03	1.5628E−03	9.7840E−04	1.7831E−02
800	0.00125	5.4169E−04	6.3721E−04	5.0980E−04	5.2160E−04	4.6025E−04	8.7931E−04	6.1114E−04	1.7908E−02

Table 2. Test 1. $L_1$ norm error of the ADER2-WAF in three regions at $t=0.\text{14}$ with $C_{\text{CFL}}=0.37$ and $C_{\text{CFL}}=0.97$ .

$M$	up to the shock, $x\in \left[0,0.833\right)$		shock, $x\in \left[0.833,0.857\right]$		sum, $x\in \left[0,0.857\right]$
	0.37	0.95	0.37	0.95	0.37	0.95
100	1.1561E−03	2.4874E−03	3.3156E−03	1.4404E−03	4.4717E−03	3.9278E−03
200	5.2670E−04	1.3041E−03	2.1741E−03	1.2618E−03	2.7008E−03	2.5659E−03
400	2.5735E−04	6.4539E−04	1.0255E−03	7.7361E−04	1.2829E−03	1.4190E−03
800	1.3030E−04	3.2895E−04	4.1139E−04	3.0826E−04	5.4169E−04	6.3721E−04

in the shock transition zone with $C_{\text{CFL}}=0.37$ is dominant over the error caused by numerical oscillations with $C_{\text{CFL}}=0.95$ . However, for $M=400$ and $M=800$ , the error caused by numerical oscillations with $C_{\text{CFL}}=0.95$ is significantly greater than the error produced in the shock transition zone when $C_{\text{CFL}}=0.37$ is utilized, resulting in a smaller $L_1$ error with $C_{\text{CFL}}=0.37$ than with $C_{\text{CFL}}=0.95$ . As an illustration, we have included the $L_1$ error for ADER2-WAF generated in the region up to the shock’s position, in the shock transition zone, and their sum in Table 2. Observe that for $M=100$ and $M=200$ , the discrepancy in errors between the two CFL coefficient values in the region preceding the shock’s position is less than the discrepancy in errors within the shock transition zone, resulting in a larger $L_1$ error with $C_{\text{CFL}}=0.37$ .

4.3. Test 2: Left and Right Rarefaction

Consider the RP (2), (3), (9) with $U_L={\left[1,-5\right]}^{\top }$ , $U_R={\left[1,5\right]}^{\top }$ , $\left(x,t\right)\in \left[0,1\right]\times \left(0,+\infty \right)$ and $x_0=0.5$ . The numerical solution is calculated at the output time of $t_{out}=0.05$ (the original test [4] had $x\in \left(0,50\right)$ , with $x_0=25$ and $t_{out}=2.5$ ). The initial data of this test produces two symmetric rarefactions. The two waves exhibit symmetry with respect to the line $x=0.5$ . Figure 10, right, illustrates the structure of the solution to the RP in the $x-t$ plane. The exact solution is

$U\left(x,t\right)=\{\begin{array}{ll}{U}_L,\hfill & \text{if}S\le S_{HL},\hfill \\ U_{Lfan}\left(x,t\right),\hfill & \text{if}{S}_{HL}<S<S_{TL},\hfill \\ U^{\text{*}},\hfill & \text{if}{S}_{TL}\le S\le S_{TR},\hfill \\ U_{Rfan}\left(x,t\right),\hfill & \text{if}{S}_{TR}<S<S_{HR},\hfill \\ U_R,\hfill & \text{if}{S}_{HR}\le S,\hfill \end{array}S=\frac{x-0.5}{t},$

where $S_{HL}\approx -8.132092$ , $S_{TL}\approx -0.632092$ , $S_{TR}\approx 0.632092$ , $S_{HR}\approx 8.132092$ , $U^{\text{*}}\approx {\left[0.040728,0\right]}^T$ . The values of $U_{Lfan}$ and $U_{Rfan}$ are obtained by (10) and (11) accordingly.

Upon examining the exact solution, it is evident that for $t_{out}<0.5/S_{HR}\approx 0.061$ , the solution exhibits a rich structure with five distinct states characterized by smooth parts and derivative discontinuities. The numerical findings demonstrate that the primary challenge faced by the ADERk-WAF scheme in this test is the appropriate management of the tails of the two rarefactions and the constant state (star region) between these waves.

Given that the star region was the greatest challenge for the ADER-WAF scheme, we investigated which CFL coefficient values yielded the minimal error in that region. The answer to this problem is more complex than that of Test 1. The range of CFL coefficient values is contingent upon the order $k$ of the ADERk-WAF scheme and the spatial step size $\text{Δ}x$ . In certain instances, the range also depends on whether we seek the optimal approximation in the star region for $h\left(x,t\right)$ or $u\left(x,t\right)$ . Table 3 and Table 4 illustrate the variations, presenting three CFL coefficient values that yielded the least $L_1$ error in the star region for $u\left(x,t\right)$ and $h\left(x,t\right)$ associated with each ADERk-WAF scheme and spatial step size $\text{Δ}x$ . These values are arranged in ascending order of the $L_1$ error obtained when utilized in ADERk-WAF. Regrettably, a uniform value for $C_{\text{CFL}}$ does not exist for the ADERk-WAF schemes in this Test. The numerical solutions for the velocity are plotted in Figures 11-14 for $M=100$ and in Figures 15-18 for $M=400$ . Utilizing CFL coefficient values from Table 3 and Table 4, the ADER-WAF scheme yields numerical results comparable to those depicted in the aforementioned figures. The $L_1$ errors are provided in Table 5.

Figures 11(b)-14(b) demonstrate that, for $C_{\text{CFL}}=0.95$ , the ADERk-WAF failed to fully realize the constant state on the 100-cell mesh. It is completely absent. The star region was partially captured by ADER2-WAF and ADER3-WAF (for $M=100$ ) for the selected values of $C_{\text{CFL}}$ from Table 3 (see Figure 11(b) and Figure 12(b)). Unexpectedly, ADER4-WAF and ADER5-WAF were unable to reproduce the star region for any of the investigated $C_{\text{CFL}}$ values in this Test. The values $C_{\text{CFL}}=0.61$ and $C_{\text{CFL}}=0.58$ , displayed in Figure 13(b) and Figure 14(b), respectively, produced the minimal $L_1$ error in the star region. Nevertheless, as illustrated in Figure 13(b) and Figure 14(b), these values of $C_{\text{CFL}}$ yielded an improved approximation to the rarefactions in comparison to $C_{\text{CFL}}=0.95$ .

The mesh refinement improves the realization of the constant space, but it causes overshoots and undershoots near the tails of the left and right rarefaction, as shown in Figures 15(b)-18(b) for $M=400$ . The magnitude of these overshoots and

Table 3. Test 2. CFL coefficient values for ADER2-WAF and ADER3-WAF.

$M$	ADER2-WAF		ADER3-WAF
	$u\left(x,t\right)$	$h\left(x,t\right)$	$u\left(x,t\right)$	$h\left(x,t\right)$
100	0.71, 0.68, 0.69	0.46, 0.45, 0.44	0.57, 0.56, 0.58	0.68, 0.67, 0.66
200	0.58, 0.57, 0.59	0.42, 0.41, 0.40	0.65, 0.66, 0.64	0.65, 0.66, 0.64
400, 800	0.42, 0.41, 0.40	0.42, 0.41, 0.40	0.65, 0.66, 0.64	0.65, 0.66, 0.64

Table 4. Test 2. CFL coefficient values for ADER4-WAF and ADER5-WAF.

$M$	ADER4-WAF		ADER5-WAF
	$u\left(x,t\right)$	$h\left(x,t\right)$	$u\left(x,t\right)$	$h\left(x,t\right)$
100	0.61, 0.60, 0.59	0.17, 0.19, 0.18	0.58, 0.59, 0.57	0.20, 0.21, 0.22
200	0.22, 0.21, 0.20	0.22, 0.21, 0.20	0.25, 0.27, 0.26	0.25, 0.27, 0.26
400, 800	0.19, 0.21, 0.20	0.19, 0.21, 0.20	0.25, 0.24, 0.26	0.25, 0.24, 0.26

Table 5. Test 2. Comparison study of the $L_1$ norm error of the ADERk-WAF schemes at $t=0.0\text{5}$ for two distinct values of the CFL coefficient.

$M$	$\Delta x$	ADER2-WAF		ADER3-WAF		ADER4-WAF		ADER5-WAF
		$L_1$ norm error
		0.42	0.95	0.65	0.95	0.19	0.95	0.25	0.95
200	0.005	1.6380E−02	2.6017E−02	1.4446E−02	1.2111E−02	1.4045E−02	1.5620E−02	1.3898E−02	2.7788E−02
400	0.0025	8.2063E−03	1.3231E−02	7.2296E−03	6.0804E−03	6.9899E−03	7.8202E−03	6.9420E−03	2.0350E−02
800	0.00125	4.1070E−03	6.6624E−03	3.6146E−03	3.0413E−03	3.4835E−03	3.9106E−03	3.4827E−03	1.6592E−02

Figure 10. The structure of the solution of the RP in the $x-t$ plane for Test 1 (left) and Test 2 (right).

Figure 11. Numerical solution of Test 2 on 100 cells by ADER2-WAF scheme and its partial enlarged view.

Figure 12. Numerical solution of Test 2 on 100 cells by ADER3-WAF scheme and its partially enlarged view.

Figure 13. Numerical solution of Test 2 on 100 cells by ADER4-WAF scheme and its partially enlarged view.

Figure 14. Numerical solution of Test 2 on 100 cells by ADER5-WAF scheme and its partially enlarged view.

Figure 15. Numerical solution of Test 2 on 100 cells by ADER2-WAF scheme and its partially enlarged view.

Figure 16. Numerical solution of Test 2 on 100 cells by ADER3-WAF scheme and its partially enlarged view.

Figure 17. Numerical solution of Test 2 on 100 cells by ADER4-WAF scheme and its partially enlarged view.

Figure 18. Numerical solution of Test 2 on 100 cells by ADER5-WAF scheme and its partially enlarged view.

undershoots (both in terms of height and width) vary depending on the order $k$ of the ADERk-WAF scheme and the value of $C_{\text{CFL}}$ . Figures 15(b)-18(b) indicate that, generally, the constant state is more accurately approximated as the order $k$ increases for both values of $C_{\text{CFL}}$ . Nonetheless, as illustrated in Figures 15(b)-18(b), utilizing the chosen values of $C_{\text{CFL}}$ yields a significantly improved approximation of the rarefactions, their tails’ positions, and the star region.

4.4. Test 3: Left and Right Shock

We solve the RP (2), (3), (9) with $U_L={\left[1,0.5\right]}^{\top }$ , $U_R={\left[1,-0.5\right]}^{\top }$ , $\left(x,t\right)\in \left[0,1\right]\times \left(0,+\infty \right)$ and $x_0=0.5$ . The numerical solution is computed at the output time of $t_{out}=0.1$ . The initial data of this test produces two symmetric shocks that propagate in opposite directions. The two waves exhibit symmetry with respect to the line $x=0.5$ . Figure 27, left, shows the structure of the solution to the RP in the $x-t$ plane. The exact solution is

$U\left(x,t\right)=\{\begin{array}{ll}{U}_L,\hfill & \text{if}S<S_{SL},\hfill \\ U^{\text{*}},\hfill & \text{if}{S}_{SL}<S<S_{SR},\hfill \\ U_R,\hfill & \text{if}{S}_{SR}<S,\hfill \end{array}S=\frac{x-0.5}{t},$

where $S_{SL}\approx -3.018779$ , $S_{SR}\approx 3.018779$ , $U^{\text{*}}\approx {\left[1.165630,0\right]}^{\top }$ .

Analysis of the exact solution reveals that, for $t_{out}<0.5/S_{SR}\approx 0.166$ , the solution consists of three distinct constant states characterized by jump discontinuities. The numerical results reveal that the principal problem encountered by the ADERk-WAF was producing a solution devoid of spurious oscillations around the shocks and in the constant state between them (star region). Decreasing the CFL coefficient value can alleviate these oscillations; nevertheless, it will lead to slightly increased errors in the regions surrounding the shocks.

Our numerical studies, which are not presented here, indicated that utilising the values 0.20, 0.24, 0.30, 0.40, and 0.60 for $C_{\text{CFL}}$ in ADER2-WAF, as well as 0.60 and 0.61 in ADERk-WAF for $k=\text{3,4}$ , and 5, successfully restricts numerical oscillation in the specified regions. These CFL coefficient values are valid for spatial step sizes of up to $\text{Δ}x=0.00125$ . Employing a larger or smaller spatial step size may expand or contract the range of the CFL coefficient values, respectively. In contrast to Test 2, a uniform value of $C_{\text{CFL}}$ , specifically $C_{\text{CFL}}=0.60$ , is applicable to all ADERk-WAF schemes and spatial step sizes $\text{Δ}x$ .

Figures 19-26 display the calculated water depth $h\left(x,t\right)$ with $C_{\text{CFL}}=0.60$ against the solution obtained with $C_{\text{CFL}}=0.95$ . The computed water depth is plotted at $t_{out}=0.1$ with $M=100$ (Figures 19-22) and $M=400$ (Figures 23-26). The $L_1$ norm error of cell averages for both CFL coefficient values is provided in Table 6.

Figures 19-26 illustrate that the ADERk-WAF demonstrates greater efficiency with $C_{\text{CFL}}=0.60$ in comparison to $C_{\text{CFL}}=0.95$ concerning numerical oscillations. The distinction is significant with ADER4-WAF and even more apparent with ADER5-WAF, particularly when a smaller spatial step size is employed; refer to Figure 21(b) and Figure 22(b) for $M=100$ and Figure 25 and Figure 26 for $M=400$ .

The shocks are resolved without numerical oscillation, in contrast to the case when $C_{\text{CFL}}=0.95$ is employed. Nevertheless, the numerical solution derived with reduced values for $C_{\text{CFL}}$ exhibits a more rounded shape in the “corners” caused by the shocks; refer to Figure 19(b) and Figure 20(b). This resulted in increased errors in those regions when $C_{\text{CFL}}=0.60$ was utilized, overshadowing the errors induced by numerical oscillations, as observed in Test 1. This explains why the $L_1$ error for ADER2-WAF ($M$ = 100, 200, 400) and ADER3-WAF ($M=100,200$ ) is greater when $C_{\text{CFL}}=0.60$ is used in comparison to $C_{\text{CFL}}=0.95$ ; refer to Table 6. In the remaining cases presented in Table 6, numerical oscillations predominated, leading to a higher $L_1$ error when $C_{\text{CFL}}=0.95$ was applied.

It is also noteworthy that for $M=400$ and $C_{\text{CFL}}=0.95$ , the oscillations in the numerical solution obtained with ADER2-WAF are scarcely discernible in the star region away from the shocks’ location; refer to Figure 23(c). The oscillations

Table 6. Test 3. Comparison study of the $L_1$ norm error of the ADERk-WAF schemes at $t=0.\text{1}$ for two distinct values of the CFL coefficient.

$M$	$\Delta x$	ADER2-WAF		ADER3-WAF		ADER4-WAF		ADER5-WAF
		$L_1$ norm error
		0.60	0.95	0.60	0.95	0.60	0.95	0.60	0.95
100	0.01	4.3884E−03	2.1926E−03	3.6632E−03	2.8176E−03	3.1206E−03	3.2583E−03	3.1098E−03	6.3580E−03
200	0.005	1.4067E−03	8.3955E−04	1.0449E−03	8.8885E−04	8.1664E−04	1.4837E−03	1.1136E−03	5.5693E−03
400	0.0025	8.1357E−04	7.9129E−04	7.0994E−04	8.1941E−04	6.2994E−04	1.4933E−03	7.8303E−04	6.6263E−03
800	0.00125	2.7523E−04	3.9518E−04	1.9629E−04	3.9686E−04	1.6054E−04	1.0646E−03	6.7788E−04	1.0534E−02

Figure 19. Numerical solution of Test 3 on 100 cells by ADER2-WAF scheme and its partially enlarged view.

Figure 20. Numerical solution of Test 3 on 100 cells by ADER3-WAF scheme and its partially enlarged view.

Figure 21. Numerical solution of Test 3 on 100 cells by ADER4-WAF scheme and its partially enlarged view.

Figure 22. Numerical solution of Test 3 on 100 cells by ADER5-WAF scheme and its partially enlarged view.

Figure 23. Numerical solution of Test 3 on 400 cells by ADER2-WAF scheme and its partially enlarged views.

Figure 24. Numerical solution of Test 3 on 400 cells by ADER3-WAF scheme and its partially enlarged views.

Figure 25. Numerical solution of Test 3 on 400 cells by ADER4-WAF scheme and its partially enlarged views.

Figure 26. Numerical solution of Test 3 on 400 cells by ADER5-WAF scheme and its partially enlarged views.

are becoming apparent with ADER3-WAF (Figure 24(c)) and intensify with ADER4-WAF and ADER5-WAF (Figure 25(c) and Figure 26(c)).

4.5. Test 4: Left Shock and Right Rarefaction

We solve the RP (2), (3), (9) with $U_L={\left[2,3.5\right]}^{\top }$ , $U_R={\left[3,3\right]}^{\top }$ , $\left(x,t\right)\in \left[0,1\right]\times \left(0,+\infty \right)$ and $x_0=0.5$ . The numerical solution is computed at the output time of $t_{out}=0.05$ . The initial data for Test 4 was chosen to induce a left shock and a narrowly defined right rarefaction in the solution. Figure 27, right, displays the structure of the solution to the RP in the $x-t$ plane. The exact solution is

$U\left(x,t\right)=\{\begin{array}{ll}{U}_L,\hfill & \text{if}S<S_{SL},\hfill \\ U^{\text{*}},\hfill & \text{if}{S}_{SL}<S\le S_{TR},\hfill \\ U_{Rfan}\left(x/t\right),\hfill & \text{if}{S}_{TR}<S<S_{HR},\hfill \\ U_R,\hfill & \text{if}{S}_{HR}\le S,\hfill \end{array}S=\frac{x-0.5}{t},$

Figure 27. The structure of the solution of the RP in the $x-t$ plane for Test 3 (left) and Test 4 (right).

where $S_{SL}\approx -3.770040$ , $S_{TR}\approx 5.486301$ , $S_{HR}\approx 6.424942$ , $U^{\text{*}}\approx {\left[2.663932,0.996948\right]}^{\top }$ , and $U_{Rfan}$ is determined by (11).

Investigating the exact solution indicates that, when $t_{out}<0.5/S_{HR}\approx 0.077$ , the solution displays a rich structure comprising four distinct states distinguished by smooth regions, derivative discontinuities, and a jump discontinuity. The primary issue faced by the ADERk-WAF method was producing a solution that addressed the shock, the star region, and the tail of the rarefaction without inducing oscillations in those regions. The solution features a thin zone with steep gradients, attributable to the narrowness of the rarefaction, complicating its accurate representation on a 100-cell mesh.

A multitude of values for $C_{\text{CFL}}$ can be utilized with the ADER2-WAF to achieve a solution devoid of oscillations. The values are 0.23:0.01:0.30, 0.32:0.01:0.38, 0.42, 0.43, 0.56, and 0.57. The options for $C_{\text{CFL}}$ are significantly limited for ADERk-WAF, where $k=\text{3},\text{4},\text{5}$ . For $C_{\text{CFL}}$ with ADER3-WAF, we may utilize the range of 0.50:0.01:0.59; for $C_{\text{CFL}}$ with ADER4-WAF, the range of 0.54:0.01:0.59; and for $C_{\text{CFL}}$ with ADER5-WAF, the specific values of 0.55, 0.56, 0.58, and 0.59. These CFL coefficient values are valid for spatial step sizes of up to $\text{Δ}x=0.00125$ . Employing a larger or smaller spatial step size may expand or contract the range of the CFL coefficient values, respectively. To employ a uniform value for $C_{\text{CFL}}$ in ADERk-WAF, we opted for a value of 0.56, without suggesting that it is the most suitable option for each ADERk-WAF independently.

Figures 28-35 depict the computed water depth $h\left(x,t\right)$ with $C_{\text{CFL}}=0.56$ against the solution derived with $C_{\text{CFL}}=0.95$ . The calculated water depth is plotted at $t_{out}=0.05$ with $M=100$ (Figures 28-31) and $M=400$ (Figures 32-35). The $L_1$ norm error of cell averages for both CFL coefficient values is given in Table 7.

Figures 28-35 indicate that the ADERk-WAF exhibits better performance with $C_{\text{CFL}}=0.56$ compared to $C_{\text{CFL}}=0.95$ regarding numerical oscillations in the star region and subsequent to the shock. However, an undershoot accompanied by oscillations prior to the rarefaction’s tail is inevitable for both values of $C_{\text{CFL}}$ . This effect becomes more noticeable as the order $k$ of the ADERk-WAF increases. Even so, the undershoot is more constrained when $C_{\text{CFL}}=0.56$ is

Table 7. Test 4. Comparison study of the $L_1$ norm error of the ADERk-WAF schemes at $t=0.0\text{5}$ for two distinct values of the CFL coefficient.

$M$	$\Delta x$	ADER2-WAF		ADER3-WAF		ADER4-WAF		ADER5-WAF
		$L_1$ norm error
		0.56	0.95	0.56	0.95	0.56	0.95	0.56	0.95
100	0.01	1.9990E−02	2.0060E−02	1.8604E−02	1.5927E−02	1.9180E−02	1.8132E−02	1.8415E−02	2.7788E−02
200	0.005	8.5406E−03	9.8360E−03	7.7080E−03	7.4430E−03	7.9815E−03	8.8488E−03	8.2340E−03	2.3091E−02
400	0.0025	4.2135E−03	5.1817E−03	3.5843E−03	4.4731E−03	3.6053E−03	5.4612E−03	4.1629E−03	2.2517E−02
800	0.00125	2.4498E−03	3.3087E−03	2.2051E−03	3.2944E−03	2.2089E−03	3.9794E−03	3.0156E−03	2.5087E−02

Figure 28. Numerical solution of Test 4 on 100 cells by ADER2-WAF scheme and its partially enlarged view.

Figure 29. Numerical solution of Test 4 on 100 cells by ADER3-WAF scheme and its partially enlarged view.

Figure 30. Numerical solution of Test 4 on 100 cells by ADER4-WAF scheme and its partially enlarged view.

Figure 31. Numerical solution of Test 4 on 100 cells by ADER5-WAF scheme and its partially enlarged view.

Figure 32. Numerical solution of Test 4 on 400 cells by ADER2-WAF scheme and its partially enlarged views.

Figure 33. Numerical solution of Test 4 on 400 cells by ADER3-WAF scheme and its partially enlarged views.

Figure 34. Numerical solution of Test 4 on 400 cells by ADER4-WAF scheme and its partially enlarged views.

Figure 35. Numerical solution of Test 4 on 400 cells by ADER5-WAF scheme and its partially enlarged views.

used, with the exception of the ADER4-WAF on a 100-cell mesh; refer to Figure 30(b).

Utilising $C_{\text{CFL}}=0.56$ resolves the shock without oscillations in contrast to $C_{\text{CFL}}=0.95$ . The rarefaction is more accurately represented (not displayed here) when $C_{\text{CFL}}=0.56$ is applied with ADER2-WAF and ADER5-WAF, and when $C=0.95$ is used with ADER3-WAF and ADER4-WAF. Analogous conclusions to those in Test 1 and Test 3 can be drawn for the discrepancies in $L_1$ error between the two values of $C_{\text{CFL}}$ in Table 7 and the numerical performance of ADER2-WAF on a 400-cell mesh.

5. Conclusions and Discussion

In this study, we explore the ADERk-WAF method for solving the RP for homogeneous 1D SWEs. Our motivation for the study was to demonstrate the significance of cautiously selecting the CFL coefficient value while employing the ADER-WAF method. This fact was illustrated by several scenarios from the SWEs, indicating that the approach’s poor performance (especially for the fourth and fifth order) was attributable to the high CFL coefficient value rather than the method itself. Consequently, the objective was to identify an appropriate range of values for the CFL coefficient to ensure that its application with the ADERk-WAF scheme yields solutions with diminished or eradicated spurious oscillations for the problems examined in this paper. Nonetheless, it is not always feasible to eliminate oscillations by reducing the CFL coefficient value of a scheme. In our study [39], we utilized numerical methods from the WENO family for identical experiments, but oscillations produced by high-order WENO schemes could not be mitigated by reducing the value of $C_{\text{CFL}}$ .

Initially, our goal was to outline the ADER-WAF procedure for the SWEs and supply the requisite equations to facilitate its coding. As far as we know, certain ones of these have not been reported in the literature. These include the exact expressions for 1) the values of the $\nu$ -th derivative ($\nu =\text{1},\text{2},\text{3},\text{4}$ ) of the WENO reconstruction polynomials $p_i^{\left(k,s\right)}\left(x\right)$ at the cell interfaces $x=x_{i+\frac{1}{2}}$ , 2) the partial derivatives of $U\left(x,t\right)$ as functions of the partial derivatives of the conserved variables $u_i\left(x,t\right)$ , $i=1,2$ , for the SWEs, 3) the time derivatives of the flux $F\left(U\left(x,t\right)\right)$ as functions of the time derivatives of the conserved variables $u_i\left(x,t\right)$ , $i=1,2$ for the SWEs, and 4) the solution to the linear RP for the SWEs in terms of the conserved variables $u_i\left(x,t\right)$ , $i=1,2$ .

After conducting exhaustive testing, our objective was to provide a suitable range of values for the CFL coefficient for each ADERk-WAF scheme. Our numerical results demonstrate that adjusting $C_{\text{CFL}}$ to an appropriate value for the ADERk-WAF scheme in the examined tests effectively eliminates spurious oscillations around shocks and in the star region, distant from the rarefactions’ tail position. By eliminating oscillations, we refer to their invisibility at a certain magnification level. This is particularly apparent with the ADER4-WAF and ADER5-WAF schemes. Nonetheless, certain oscillations toward the tail of the rarefactions were inevitable. While reducing the $C_{\text{CFL}}$ value did not eradicate oscillations near the tail of the rarefactions, it did manage to control them.

Prior to concluding, we wish to present two final remarks. Initially, it was observed that on a refined mesh, there was a slight discrepancy between the two values of $C_{\text{CFL}}$ when utilized with the ADER2-WAF. Therefore, the ADER2-WAF scheme can be employed with $C_{\text{CFL}}=0.95$ for very long time evolution problems to decrease CPU time. Nonetheless, with ADER4-WAF and particularly with ADER5-WAF, we highly recommend reducing the value of $C_{\text{CFL}}$ . Secondly, it was noted that employing $C_{\text{CFL}}=0.95$ yielded improved shock resolution and, in certain instances, a more accurate approximation of rarefaction, with the exception of the ADER5-WAF. If the resolution of shock or rarefaction is critical, one may employ the value $C_{\text{CFL}}=0.95$ , albeit at the expense of oscillations in specific regions.

The ADER-WAF approach described in this research can be applied to any hyperbolic conservation laws with certain code adjustments. Our future plans include the implementation of the ADER-WAF scheme to address more intricate models and conduct a comparative analysis with methods from the WENO family.

Acknowledgements

Research of Pavlos Stampolidis has been financially supported by General Secretariat for Research and Technology (GSRT) and the Hellenic Foundation for Research and Innovation (HFRI). The authors would also like to thank the anonymous referees for their insightful comments and recommendations.

Appendices

Appendix A. The Smoothness Indicators

For $k=2$ , (19) yields [31]

${\beta }_0^2\left(u,i\right)={\left(u_{i+1}-u_i\right)}^2,{\beta }_1^2\left(u,i\right)={\left(u_i-u_{i-1}\right)}^2.$ (29)

For $k=3$ , (19) yields [31]

$\begin{array}{l}{\beta }_0^3\left(u,i\right)=\frac{13}{12}{\left(u_i-2u_{i+1}+u_{i+2}\right)}^2+\frac{1}{4}{\left(3u_i-4u_{i+1}+u_{i+2}\right)}^2,\\ {\beta }_1^3\left(u,i\right)=\frac{13}{12}{\left(u_{i-1}-2u_i+u_{i+1}\right)}^2+\frac{1}{4}{\left(u_{i-1}-u_{i+1}\right)}^2,\\ {\beta }_2^3\left(u,i\right)=\frac{13}{12}{\left(u_{i-2}-2u_{i-1}+u_i\right)}^2+\frac{1}{4}{\left(u_{i-2}-4u_{i-1}+3u_i\right)}^2.\end{array}$ (30)

For $k=4$ , (19) yields [17]

$\begin{array}{l}{\beta }_0^4\left(u,i\right)=\frac{1}{240}[u_i\left(2107u_i-9402u_{i+1}+7042u_{i+2}-1854u_{i+3}\right)+u_{i+1}(11003u_{i+1}\\ -17246u_{i+2}+4642u_{i+3})+u_{i+2}\left(7043u_{i+2}-3882u_{i+3}\right)+547u_{i+3}^2],\\ {\beta }_1^4\left(u,i\right)=\frac{1}{240}[u_{i-1}\left(547u_{i-1}-2522u_i+1922u_{i+1}-494u_{i+2}\right)+u_i(3443u_i\\ -5966u_{i+1}+1602u_{i+2})+u_{i+1}\left(2843u_{i+1}-1642u_{i+2}\right)+267u_{i+2}^2],\\ {\beta }_2^4\left(u,i\right)=\frac{1}{240}[u_{i-2}\left(267u_{i-2}-1642u_{i-1}+1602u_i-494u_{i+1}\right)+u_{i-1}(2843u_{i-1}\\ -5966u_i+1922u_{i+1})+u_i\left(3443u_i-2522u_{i+1}\right)+547u_{i+1}^2],\\ {\beta }_3^4\left(u,i\right)=\frac{1}{240}[u_{i-3}\left(547u_{i-3}-3882u_{i-2}+4642u_{i-1}-1854u_i\right)+u_{i-2}(7043u_{i-2}\\ -17246u_{i-1}+7042u_i)+u_{i-1}\left(11003u_{i-1}-9402u_i\right)+2107u_i^2],\end{array}$ (31)

where ${\beta }_i^4$ is equal to $\frac{1}{240}{\text{IS}}_{3-i}^4$ , in [17], $i=0,1,2,3$ .

For $k=5$ , (19) yields [17]

$\begin{array}{c}{\beta }_0^5\left(u,i\right)=\frac{1}{5040}[u_i\left(107918u_i-649501u_{i+1}+758823u_{i+2}-411487u_{i+3}+86329u_{i+4}\right)\\ +u_{i+1}\left(1020563u_{i+1}-2462076u_{i+2}+1358458u_{i+3}-288007u_{i+4}\right)\\ +u_{i+2}\left(1521393u_{i+2}-1704396u_{i+3}+364863u_{i+4}\right)\\ +u_{i+3}\left(482963u_{i+3}-208501u_{i+4}\right)+22658u_{i+4}^2],\end{array}$

$\begin{array}{c}{\beta }_1^5\left(u,i\right)=\frac{1}{5040}[u_{i-1}\left(22658u_{i-1}-140251u_i+165153u_{i+1}-88297u_{i+2}+18079u_{i+3}\right)\\ +u_i\left(242723u_i-611976u_{i+1}+337018u_{i+2}-70237u_{i+3}\right)\\ +u_{i+1}\left(406293u_{i+1}-464976u_{i+2}+99213u_{i+3}\right)\\ +u_{i+2}\left(138563u_{i+2}-60871u_{i+3}\right)+6908u_{i+3}^2],\end{array}$

$\begin{array}{c}{\beta }_2^5\left(u,i\right)=\frac{1}{5040}[u_{i-2}\left(6908u_{i-2}-51001u_{i-1}+67923u_i-38947u_{i+1}+8209u_{i+2}\right)\\ +u_{i-1}\left(104963u_{i-1}-299076u_i+179098u_{i+1}-38947u_{i+2}\right)\end{array}$

$\begin{array}{c}+u_i\left(231153u_i-299076u_{i+1}+67923u_{i+2}\right)\\ +u_{i+1}\left(104963u_{i+1}-51001u_{i+2}\right)+6908u_{i+2}^2],\end{array}$ (32)

$\begin{array}{c}{\beta }_3^5\left(u,i\right)=\frac{1}{5040}[u_{i-3}\left(6908u_{i-3}-60871u_{i-2}+99213u_{i-1}-70237u_i+18079u_{i+1}\right)\\ +u_{i-2}\left(138563u_{i-2}-464976u_{i-1}+337018u_i-88297u_{i+1}\right)\\ +u_{i-1}(406293u_{i-1}-611976u_i+165153u_{i+1})\\ +u_i\left(242723u_i-140251u_{i+1}\right)+22658u_{i+1}^2],\end{array}$

$\begin{array}{c}{\beta }_4^5\left(u,i\right)=\frac{1}{5040}[u_{i-4}\left(22658u_{i-4}-208501u_{i-3}+364863u_{i-2}-288007u_{i-1}+86329u_i\right)\\ +u_{i-3}\left(482963u_{i-3}-1704396u_{i-2}+1358458u_{i-1}-411487u_i\right)\\ +u_{i-2}(1521393u_{i-2}-2462076u_{i-1}+758823u_i)\\ +u_{i-1}\left(1020563u_{i-1}-649501u_i\right)+107918u_i^2],\end{array}$

where ${\beta }_i^5$ is equal to $\frac{1}{5040}{\text{IS}}_{4-i}^5$ , in [17], $i=0,\cdots ,4$ .

Appendix B. Tables of Coefficients

Table A1. the coefficients $c_{k,s,j}^{\left(0\right)}$ [31] of Procedure 1.

$k$	$s$	$j=0$	$j=1$	$j=2$	$j=3$	$j=4$
2	−1	3/2	−1/2
	0	1/2	1/2
	1	−1/2	3/2
3	−1	11/6	−7/6	1/3
	0	1/3	5/6	−1/6
	1	−1/6	5/6	1/3
	2	1/3	−7/6	11/6
4	−1	25/12	−23/12	13/12	−1/4
	0	1/4	13/12	−5/12	1/12
	1	−1/12	7/12	7/12	−1/12
	2	1/12	−5/12	13/12	1/4
	3	−1/4	13/12	−23/12	25/12
5	−1	137/60	−163/60	137/60	−21/20	1/5
	0	1/5	77/60	−43/60	17/60	−1/20
	1	−1/20	9/20	47/60	−13/60	1/30
	2	1/30	−13/60	47/60	9/20	−1/20
	3	−1/20	17/60	−43/60	77/60	1/5
	4	1/5	−21/20	137/60	−163/60	137/60

Table A2. The coefficients $c_{k,s,j}^{\left(1\right)}$ of Procedure 1.

$k$	$s$	$j=0$	$j=1$	$j=2$	$j=3$	$j=4$
2	−1, 0, 1	−1/Δx	1/Δx
3	−1	−2/Δx	3/Δx	−1/Δx
	0	−1/Δx	1/Δx	0
	1	0	−1/Δx	1/Δx
	2	1/Δx	−3/Δx	2/Δx
4	−1	−35/12Δx	23/4Δx	−15/4Δx	11/12Δx
	0	−11/12Δx	3/4Δx	1/4Δx	−1/12Δx
	1	1/12Δx	−5/4Δx	5/4Δx	−1/12Δx
	2	1/12Δx	−1/4Δx	−3/4Δx	11/12Δx
	3	−11/12Δx	15/4Δx	−23/4Δx	35/12Δx
5	−1	−15/4Δx	109/12Δx	−35/4Δx	17/4Δx	−5/6Δx
	0	−5/6Δx	5/12Δx	3/4Δx	−5/12Δx	1/12Δx
	1	1/12Δx	−5/4Δx	5/4Δx	−1/12Δx	0
	2	0	1/12Δx	−5/4Δx	5/4Δx	−1/12Δx
	3	−1/12Δx	5/12Δx	−3/4Δx	−5/12Δx	5/6Δx
	4	5/6Δx	−17/4Δx	35/4Δx	−109/12Δx	15/4Δx

Table A3. the coefficients $c_{k,s,j}^{\left(2\right)}$ of Procedure 1.

$k$	$s$	$j=0$	$j=1$	$j=2$	$j=3$	$j=4$
3	−1, ..., 2	1/Δx2	−2/Δx2	1/Δx2
4	−1	5/2Δx2	−13/2Δx2	11/2Δx2	−3/2Δx2
	0	3/2Δx2	−7/2Δx2	5/2Δx2	−1/2Δx2
	1	1/2Δx2	−1/2Δx2	−1/2Δx2	1/2Δx2
	2	−1/2Δx2	5/2Δx2	−7/2Δx2	3/2Δx2
	3	−3/2Δx2	11/2Δx2	−13/2Δx2	5/2Δx2
5	−1	17/4Δx2	−27/2Δx2	16/Δx2	−17/2Δx2	7/4Δx2
	0	7/4Δx2	−9/2Δx2	4/Δx2	−3/2Δx2	1/4Δx2
	1	1/4Δx2	1/2Δx2	−2/Δx2	3/2Δx2	−1/4Δx2
	2	−1/4Δx2	3/2Δx2	−2/Δx2	1/2Δx2	1/4Δx2
	3	1/4Δx2	−3/2Δx2	4/Δx2	−9/2Δx2	7/4Δx2
	4	7/4Δx2	−17/2Δx2	16/Δx2	−27/2Δx2	17/4Δx2

Table A4. The coefficients $c_{k,s,j}^{\left(3\right)}$ of Procedure 1.

$k$	$s$	$j=0$	$j=1$	$j=2$	$j=3$	$j=4$
4	−1, ..., 3	−1/Δx3	3/Δx3	−3/Δx3	1/Δx3
5	−1	−3/Δx3	11/Δx3	−15/Δx3	9/Δx3	−2/Δx3
	0	−2/Δx3	7/Δx3	−9/Δx3	5/Δx3	−1/Δx3
	1	−1/Δx3	3/Δx3	−3/Δx3	1/Δx3	0
	2	0	−1/Δx3	3/Δx3	−3/Δx3	1/Δx3
	3	1/Δx3	−5/Δx3	9/Δx3	−7/Δx3	2/Δx3
	4	2/Δx3	−9/Δx3	15/Δx3	−11/Δx3	3/Δx3

Table A5. The coefficients $c_{k,s,j}^{\left(4\right)}$ of Procedure 1.

$k$	$s$	$j=0$	$j=1$	$j=2$	$j=3$	$j=4$
5	−1, ..., 4	1/Δx4	−4/Δx4	6/Δx4	−4/Δx4	1/Δx4

Table A6. Optimal weights $d_{k,s}$ [31, 17] of Procedure 1.

	$s=0$	$s=1$	$s=2$	$s=3$	$s=4$
$k=2$	2/3	1/3
$k=3$	3/10	3/5	1/10
$k=4$	4/35	18/35	12/35	1/35
$k=5$	5/126	20/63	10/21	10/23	1/126

Appendix C. ADERk-WAF Scheme

In this appendix, the time derivatives of $U\left(x,t\right)$ are obtained through the Cauchy-Kowalewski procedure. The time derivatives of the flux $F\left(U\right)$ are given as well. The following notations are employed for the sake of simplicity

$\frac{{\partial }^{\nu +\mu }}{{\partial }_x^{\nu }{\partial }_t^{\mu }}U={\partial }_x^{\nu }{\partial }_t^{\mu }U,\frac{{\partial }^{\nu +\mu }}{{\partial }_x^{\nu }{\partial }_t^{\mu }}{u}_i=u_{i\underset{\nu \text{times}}{\underbrace{x\cdots x}}\underset{\mu \text{times}}{\underbrace{t\cdots t}}},i=1,2.$

C1. Cauchy-Kowalewski Procedure

From (2), applying the chain rule, we obtain

${\partial }_{t}U=-{\partial }_{x}F\left(U\right)=-\frac{\partial F}{\partial U}{\partial }_{x}U={\left[-u_{2x},-\frac{1}{u_1^2}\left(-u_2^2{u}_{1x}+2u_1{u}_2{u}_{2x}\right)-g{u}_1{u}_{1x}\right]}^{\top }$ (33)

Differentiating (33) with respect to (w.r.t.) $x$ and w.r.t. $t$ , we have

$\begin{array}{l}{\partial }_x{\partial }_{t}U=[-u_{2xx},-\frac{2}{u_1^3}{u}_2^2{u}_{1x}^2+\frac{1}{u_1^2}\left(4u_2{u}_{1x}{u}_{2x}+u_2^2{u}_{1xx}\right)\\ -{\frac{2}{u_1}\left(u_{2x}^2+u_2{u}_{2xx}\right)-g\left(u_{1x}^2+u_1{u}_{1xx}\right)]}^{\top }\end{array}$ (34)

$\begin{array}{l}{\partial }_t^{2}U=[-u_{2xt},-\frac{2}{u_1^3}{u}_2^2{u}_{1x}{u}_{1t}+\frac{1}{u_1^2}\left(2u_2{u}_{1t}{u}_{2x}+2u_2{u}_{1x}{u}_{2t}+u_2^2{u}_{1xt}\right)\\ -{\frac{2}{u_1}\left(u_{2x}{u}_{2t}+u_2{u}_{2xt}\right)-g\left(u_{1x}{u}_{1t}+u_1{u}_{1xt}\right)]}^{\top }\end{array}$ (35)

Differentiating (34) w.r.t. $x$ and w.r.t. $t$ , and (35) w.r.t. $t$ , yields

$\begin{array}{c}{\partial }_x^2{\partial }_{t}U=[-u_{2xxx},\frac{6}{u_1^4}{u}_2^2{u}_{1x}^3-\frac{6}{u_1^3}\left(2u_2{u}_{1x}^2{u}_{2x}+u_2^2{u}_{1x}{u}_{1xx}\right)+\frac{1}{u_1^2}(6u_{1x}{u}_{2x}^2+u_2^2{u}_{1xxx}\\ +6u_2{u}_{2x}{u}_{1xx}+6u_2{u}_{1x}{u}_{2xx})-{\frac{2}{u_1}\left(3u_{2x}{u}_{2xx}+u_2{u}_{2xxx}\right)-g\left(3u_{1x}{u}_{1xx}+u_1{u}_{1xxx}\right)]}^{\top }\end{array}$ (36)

$\begin{array}{c}{\partial }_t^2{\partial }_{x}U=[-u_{2xxt},\frac{6}{u_1^4}{u}_2^2{u}_{1x}^2{u}_{1t}-\frac{2}{u_1^3}(2u_2{u}_{1x}^2{u}_{2t}+2u_2^2{u}_{1x}{u}_{1xt}+u_2^2{u}_{1t}{u}_{1xx}\\ +4u_2{u}_{1x}{u}_{2x}{u}_{1t})+\frac{1}{u_1^2}(2u_{2x}^2{u}_{1t}+2u_2{u}_{2t}{u}_{1xx}+2u_2{u}_{1t}{u}_{2xx}+u_2^2{u}_{1xxt}\\ +4u_2{u}_{1x}{u}_{2xt}+4u_{1x}{u}_{2x}{u}_{2t}+4u_2{u}_{2x}{u}_{1xt})-\frac{2}{u_1}(2u_{2x}{u}_{2xt}+u_{2t}{u}_{2xx}\\ {\stackrel{}{}+u_2{u}_{2xxt})-g\left(2u_{1x}{u}_{1xt}+u_{1t}{u}_{1xx}+u_1{u}_{1xxt}\right)]}^{\top }\end{array}$(37)

$\begin{array}{l}{\partial }_t^{3}U=[-u_{2xtt},\frac{6}{u_1^4}{u}_2^2{u}_{1x}{u}_{1t}^2-\frac{2}{u_1^3}\left(2u_2{u}_{1t}^2{u}_{2x}+4u_2{u}_{1x}{u}_{1t}{u}_{2t}+2u_2^2{u}_{1t}{u}_{1xt}+u_2^2{u}_{1x}{u}_{1tt}\right)\\ +\frac{1}{u_1^2}(4u_{1t}{u}_{2x}{u}_{2t}+4u_2{u}_{1t}{u}_{2xt}+2u_2{u}_{2x}{u}_{1tt}+2u_{1x}{u}_{2t}^2+4u_2{u}_{2t}{u}_{1xt}+2u_2{u}_{1x}{u}_{2tt}\\ +u_2^2{u}_{1xtt})-{\frac{2}{u_1}\left(2u_{2t}{u}_{2xt}+u_{2x}{u}_{2tt}+u_2{u}_{2xtt}\right)-g\left(u_{1x}{u}_{1tt}+2u_{1t}{u}_{1xt}+u_1{u}_{1xtt}\right)]}^{\top }\end{array}$ (38)

Differentiating (36) w.r.t. $x$ and w.r.t. $t$ , (37) w.r.t. $t$ , and (38) w.r.t. $t$ , we derive

$\begin{array}{c}{\partial }_x^3{\partial }_{t}U=[-u_{2xxxx},-\frac{24}{u_1^5}{u}_2^2{u}_{1x}^4+\frac{12}{u_1^4}\left(4u_2{u}_{1x}^3{u}_{2x}+3u_2^2{u}_{1x}^2{u}_{1xx}\right)-\frac{2}{u_1^3}(12u_{1x}^2{u}_{2x}^2\\ +24u_2{u}_{1x}{u}_{2x}{u}_{1xx}+12u_2{u}_{1x}^2{u}_{2xx}+3u_2^2{u}_{1xx}^2+4u_2^2{u}_{1x}{u}_{1xxx})+\frac{1}{u_1^2}(12u_{2x}^2{u}_{1xx}\\ +24u_{1x}{u}_{2x}{u}_{2xx}+8u_2{u}_{2x}{u}_{1xxx}+u_2^2{u}_{1xxxx}+12u_2{u}_{1xx}{u}_{2xx}+8u_2{u}_{1x}{u}_{2xxx})\\ -{\frac{1}{u_1}\left(6u_{2xx}^2+8u_{2x}{u}_{2xxx}+2u_2{u}_{2xxxx}\right)-g\left(3u_{1xx}^2+4u_{1x}{u}_{1xxx}+u_1{u}_{1xxxx}\right)]}^{\top }\end{array}$ (39)

$\begin{array}{c}{\partial }_t^2{\partial }_x^{2}U=[-u_{2xxxt},-\frac{24}{u_1^5}{u}_2^2{u}_{1x}^3{u}_{1t}+\frac{6}{u_1^4}(6u_2{u}_{1x}^2{u}_{2x}{u}_{1t}+3u_2^2{u}_{1x}{u}_{1t}{u}_{1xx}+3u_2^2{u}_{1x}^2{u}_{1xt}\\ +2u_2{u}_{1x}^3{u}_{2t})-\frac{2}{u_1^3}(6u_{1x}^2{u}_{2x}{u}_{2t}+6u_2{u}_{1x}{u}_{2t}{u}_{1xx}+6u_2{u}_{1x}^2{u}_{2xt}+12u_2{u}_{1x}{u}_{2x}{u}_{1xt}\\ +3u_2^2{u}_{1xt}{u}_{1xx}+3u_2^2{u}_{1x}{u}_{1xxt}+6u_2{u}_{2x}{u}_{1t}{u}_{1xx}+u_2^2{u}_{1t}{u}_{1xxx}+6u_{1x}{u}_{2x}^2{u}_{1t}\end{array}$

$\begin{array}{c}+6u_2{u}_{1x}{u}_{1t}{u}_{2xx})+\frac{1}{u_1^2}(6u_{2x}{u}_{1t}{u}_{2xx}+6u_{2x}^2{u}_{1xt}+6u_{2x}{u}_{2t}{u}_{1xx}+6u_2{u}_{2xt}{u}_{1xx}\\ +2u_2{u}_{2t}{u}_{1xxx}+6u_2{u}_{1xt}{u}_{2xx}+2u_2{u}_{1t}{u}_{2xxx}+6u_2{u}_{2x}{u}_{1xxt}+u_2^2{u}_{1xxxt}\\ +12u_{1x}{u}_{2x}{u}_{2xt}+6u_2{u}_{1x}{u}_{2xxt}+6u_{1x}{u}_{2t}{u}_{2xx})-\frac{2}{u_1}(3u_{2xt}{u}_{2xx}+3u_{2x}{u}_{2xxt}\\ {\stackrel{}{}+u_{2t}{u}_{2xxx}+u_2{u}_{2xxxt})-g\left(3u_{1xt}{u}_{1xx}+3u_{1x}{u}_{1xxt}+u_{1t}{u}_{1xxx}+u_1{u}_{1xxxt}\right)]}^{\top }\end{array}$(40)

$\begin{array}{c}{\partial }_t^3{\partial }_{x}U=[-u_{2xxtt},-\frac{24}{u_1^5}{u}_2^2{u}_{1x}^2{u}_{1t}^2+\frac{6}{u_1^4}(4u_2{u}_{1x}{u}_{2x}{u}_{1t}^2+u_2^2{u}_{1t}^2{u}_{1xx}+4u_2^2{u}_{1x}{u}_{1t}{u}_{1xt}\\ +4u_2{u}_{1x}^2{u}_{1t}{u}_{2t}+u_2^2{u}_{1x}^2{u}_{1tt})-\frac{2}{u_1^3}(2u_{2x}^2{u}_{1t}^2+8u_2{u}_{2x}{u}_{1t}{u}_{1xt}+2u_2{u}_{1t}^2{u}_{2xx}\\ +8u_{1x}{u}_{2x}{u}_{1t}{u}_{2t}+4u_2{u}_{1t}{u}_{2t}{u}_{1xx}+8u_2{u}_{1x}{u}_{2t}{u}_{1xt}+8u_2{u}_{1x}{u}_{1t}{u}_{2xt}+2u_2^2{u}_{1xt}^2\\ +2u_2^2{u}_{1t}{u}_{1xxt}+4u_2{u}_{1x}{u}_{2x}{u}_{1tt}+u_2^2{u}_{1xx}{u}_{1tt}+2u_2^2{u}_{1x}{u}_{1xtt}+2u_{1x}^2{u}_{2t}^2\\ +2u_2{u}_{1x}^2{u}_{2tt})+\frac{1}{u_1^2}(8u_{2x}{u}_{2t}{u}_{1xt}+4u_{1t}{u}_{2t}{u}_{2xx}+8u_{2x}{u}_{1t}{u}_{2xt}+8u_2{u}_{1xt}{u}_{2xt}\\ +4u_2{u}_{1t}{u}_{2xxt}+2u_{2x}^2{u}_{1tt}+2u_2{u}_{2xx}{u}_{1tt}+4u_2{u}_{2x}{u}_{1xtt}+2u_{2t}^2{u}_{1xx}\\ +8u_{1x}{u}_{2t}{u}_{2xt}+4u_2{u}_{2t}{u}_{1xxt}+4u_{1x}{u}_{2x}{u}_{2tt}+2u_2{u}_{1xx}{u}_{2tt}+4u_2{u}_{1x}{u}_{2xtt}\\ +u_2^2{u}_{1xxtt})-\frac{2}{u_1}(2u_{2xt}^2+2u_{2t}{u}_{2xxt}+u_{2xx}{u}_{2tt}+2u_{2x}{u}_{2xtt}+u_2{u}_{2xxtt})\\ {\begin{array}{c}\\ \end{array}-g\left(u_{1xx}{u}_{1tt}+2u_{1x}{u}_{1xtt}+2u_{1xt}^2+2u_{1t}{u}_{1xxt}+u_1{u}_{1xxtt}\right)]}^{\top }\end{array}$ (41)

$\begin{array}{c}{\partial }_t^{4}U=[-u_{2xttt},-\frac{24}{u_1^5}{u}_2^2{u}_{1x}{u}_{1t}^3+\frac{6}{u_1^4}(6u_2{u}_{1x}{u}_{1t}^2{u}_{2t}+3u_2^2{u}_{1t}^2{u}_{1xt}+3u_2^2{u}_{1x}{u}_{1t}{u}_{1tt}\\ +2u_2{u}_{1t}^3{u}_{2x})-\frac{2}{u_1^3}(6u_{1t}^2{u}_{2x}{u}_{2t}+6u_2{u}_{1t}{u}_{2x}{u}_{1tt}+6u_2{u}_{1t}^2{u}_{2xt}+6u_{1x}{u}_{1t}{u}_{2t}^2\\ +12u_2{u}_{1t}{u}_{2t}{u}_{1xt}+6u_2{u}_{1x}{u}_{2t}{u}_{1tt}+6u_2{u}_{1x}{u}_{1t}{u}_{2tt}+3u_2^2{u}_{1xt}{u}_{1tt}+3u_2^2{u}_{1t}{u}_{1xtt}\\ +u_2^2{u}_{1x}{u}_{1ttt})+\frac{1}{u_1^2}(6u_{2x}{u}_{2t}{u}_{1tt}+12u_{1t}{u}_{2t}{u}_{2xt}+6u_{1t}{u}_{2x}{u}_{2tt}+6u_2{u}_{2xt}{u}_{1tt}\\ +6u_2{u}_{1t}{u}_{2xtt}+6u_{2t}^2{u}_{1xt}+6u_{1x}{u}_{2t}{u}_{2tt}+6u_2{u}_{1xt}{u}_{2tt}+6u_2{u}_{2t}{u}_{1xtt}+2u_2{u}_{2x}{u}_{1ttt}\\ +2u_2{u}_{1x}{u}_{2ttt}+u_2^2{u}_{1xttt})-\frac{2}{u_1}(3u_{2xt}{u}_{2tt}+3u_{2t}{u}_{2xtt}+u_{2x}{u}_{2ttt}+u_2{u}_{2xttt})\\ {\begin{array}{c}\\ \end{array}-g\left(3u_{1xt}{u}_{1tt}+3u_{1t}{u}_{1xtt}+u_{1x}{u}_{1ttt}+u_1{u}_{1xttt}\right)]}^{\top }\end{array}$ (42)

C2. Time Derivatives of the Flux

Firstly, we differentiate the vector of fluxes $F\left(U\right)$ , as provided in (3), w.r.t. $t$ . To get the subsequent equation, we differentiate the obtained one again w.r.t. $t$ .

${\partial }_{t}F={\left[u_{2t},\frac{1}{u_1^2}\left(-u_2^2{u}_{1t}+2u_1{u}_2{u}_{2t}\right)+g{u}_1{u}_{1t}\right]}^{\top }$ (43)

${\partial }_t^{2}F={\left[u_{2tt},\frac{2}{u_1^3}{u}_2^2{u}_{1t}^2-\frac{1}{u_1^2}\left(4u_2{u}_{1t}{u}_{2t}+u_2^2{u}_{1tt}\right)+\frac{2}{u_1}\left(u_{2t}^2+u_2{u}_{2tt}\right)+g\left(u_{1t}^2+u_1{u}_{1tt}\right)\right]}^{\top }$ (44)

$\begin{array}{l}{\partial }_t^{3}F=[u_{2ttt},-\frac{6}{u_1^4}{u}_2^2{u}_{1t}^3+\frac{6}{u_1^3}\left(2u_2{u}_{1t}^2{u}_{2t}+u_2^2{u}_{1t}{u}_{1tt}\right)-\frac{1}{u_1^2}(6u_{1t}{u}_{2t}^2+6u_2{u}_{1t}{u}_{2tt}\\ +6u_2{u}_{2t}{u}_{1tt}+u_2^2{u}_{1ttt})+{\frac{2}{u_1}\left(3u_{2t}{u}_{2tt}+u_2{u}_{2ttt}\right)+g\left(3u_{1t}{u}_{1tt}+u_1{u}_{1ttt}\right)]}^{\top }\end{array}$ (45)

$\begin{array}{l}{\partial }_t^{4}F=[u_{2tttt},\frac{24}{u_1^5}{u}_2^2{u}_{1t}^4-\frac{12}{u_1^4}\left(4u_2{u}_{1t}^3{u}_{2t}+3u_2^2{u}_{1t}^2{u}_{1tt}\right)+\frac{2}{u_1^3}(12u_{1t}^2{u}_{2t}^2+24u_2{u}_{1t}{u}_{2t}{u}_{1tt}\\ +12u_2{u}_{1t}^2{u}_{2tt}+4u_2^2{u}_{1t}{u}_{1ttt}+3u_2^2{u}_{1tt}^2)-\frac{1}{u_1^2}(12u_{2t}^2{u}_{1tt}+24u_{1t}{u}_{2t}{u}_{2tt}\\ +12u_2{u}_{1tt}{u}_{2tt}+8u_2{u}_{1t}{u}_{2ttt}+8u_2{u}_{2t}{u}_{1ttt}+u_2^2{u}_{1tttt})\\ +{\frac{2}{u_1}\left(3u_{2tt}^2+4u_{2t}{u}_{2ttt}+u_2{u}_{2tttt}\right)+g\left(3u_{1tt}^2+4u_{1t}{u}_{1ttt}+u_1{u}_{1tttt}\right)]}^{\top }\end{array}$ (46)

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Kinnmark, I. (1986) The Shallow Water Wave Equations: Formulation, Analysis and Application. Springer-Verlag. https://doi.org/10.1007/978-3-642-82646-7
[2]	Stoker, J.J. (1957) Water Waves: The Mathematical Theory with Applications. Inter-science Publishers. https://doi.org/10.1007/978-3-642-82646-7
[3]	Toro, E.F. (2009) Riemann Solvers and Numerical Methods for Fluid Dynamics—A Practical Introduction. 3rd Edition, Springer-Verlag. https://doi.org/10.1007/b79761
[4]	Toro, E.F. (2001) Shock-Capturing Methods for Free-Surface Shallow Flows. Wiley and Sons Ltd.
[5]	LeVeque, R.J. (2002) Finite Volume Methods for Hyperbolic Problems. Cambridge University Press. https://doi.org/10.1017/cbo9780511791253
[6]	Thomas, J.W. (1999) Numerical Partial Differential Equations: Conservation Laws and Elliptic Equations. Springer.
[7]	Gousidou-Koutita, M. (2008) Numerical Methods with Applications to Ordinary and Partial Differential Equations. Lecture Notes for Postgraduate Studies, Aristotle University of Thessaloniki.
[8]	Harten, A. (1983) High Resolution Schemes for Hyperbolic Conservation Laws. Journal of Computational Physics, 49, 357-393. https://doi.org/10.1016/0021-9991(83)90136-5
[9]	Harten, A. (1984) On a Class of High Resolution Total-Variation-Stable Finite-Difference Schemes. SIAM Journal on Numerical Analysis, 21, 1-23. https://doi.org/10.1137/0721001
[10]	Harten, A., Osher, S., Engquist, B. and Chakravarthy, S.R. (1986) Some Results on Uniformly High-Order Accurate Essentially Nonoscillatory Schemes. Applied Numerical Mathematics, 2, 347-377. https://doi.org/10.1016/0168-9274(86)90039-5
[11]	Harten, A. and Osher, S. (1987) Uniformly High-Order Accurate Nonoscillatory Schemes. I. SIAM Journal on Numerical Analysis, 24, 279-309. https://doi.org/10.1137/0724022
[12]	Harten, A., Engquist, B., Osher, S. and Chakravarthy, S.R. (1987) Uniformly High Order Accurate Essentially Non-Oscillatory Schemes, III. Journal of Computational Physics, 71, 231-303. https://doi.org/10.1016/0021-9991(87)90031-3
[13]	Shu, C. and Osher, S. (1988) Efficient Implementation of Essentially Non-Oscillatory Shock-Capturing Schemes. Journal of Computational Physics, 77, 439-471. https://doi.org/10.1016/0021-9991(88)90177-5
[14]	Shu, C. and Osher, S. (1989) Efficient Implementation of Essentially Non-Oscillatory Shock-Capturing Schemes, II. Journal of Computational Physics, 83, 32-78. https://doi.org/10.1016/0021-9991(89)90222-2
[15]	Liu, X., Osher, S. and Chan, T. (1994) Weighted Essentially Non-Oscillatory Schemes. Journal of Computational Physics, 115, 200-212. https://doi.org/10.1006/jcph.1994.1187
[16]	Jiang, G. and Shu, C. (1996) Efficient Implementation of Weighted ENO Schemes. Journal of Computational Physics, 126, 202-228. https://doi.org/10.1006/jcph.1996.0130
[17]	Balsara, D.S. and Shu, C. (2000) Monotonicity Preserving Weighted Essentially Non-Oscillatory Schemes with Increasingly High Order of Accuracy. Journal of Computational Physics, 160, 405-452. https://doi.org/10.1006/jcph.2000.6443
[18]	Toro, E.F., Millington, R.C. and Nejad, L.A.M. (2001) Towards Very High Order Godunov Schemes. In: Toro, E.F., Ed., Godunov Methods, Springer, 907-940. https://doi.org/10.1007/978-1-4615-0663-8_87
[19]	Toro, E.F. and Titarev, V.A. (2002) Solution of the Generalized Riemann Problem for Advection-Reaction Equations. Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 458, 271-281. https://doi.org/10.1098/rspa.2001.0926
[20]	Titarev, V.A. and Toro, E.F. (2002) ADER: Arbitrary High Order Godunov Approach. Journal of Scientific Computing, 17, 609-618. https://doi.org/10.1023/a:1015126814947
[21]	Toro, E.F. and Titarev, V.A. (2005) ADER Schemes for Scalar Non-Linear Hyperbolic Conservation Laws with Source Terms in Three-Space Dimensions. Journal of Computational Physics, 202, 196-215. https://doi.org/10.1016/j.jcp.2004.06.014
[22]	Titarev, V.A. and Toro, E.F. (2005) ADER Schemes for Three-Dimensional Non-Linear Hyperbolic Systems. Journal of Computational Physics, 204, 715-736. https://doi.org/10.1016/j.jcp.2004.10.028
[23]	Toro, E.F. and Titarev, V.A. (2005) TVD Fluxes for the High-Order ADER Schemes. Journal of Scientific Computing, 24, 285-309. https://doi.org/10.1007/s10915-004-4790-8
[24]	Toro, E.F. and Titarev, V.A. (2006) Derivative Riemann Solvers for Systems of Conservation Laws and ADER Methods. Journal of Computational Physics, 212, 150-165. https://doi.org/10.1016/j.jcp.2005.06.018
[25]	Busto, S., Chiocchetti, S., Dumbser, M., Gaburro, E. and Peshkov, I. (2020) High Order ADER Schemes for Continuum Mechanics. Frontiers in Physics, 8, Article 32. https://doi.org/10.3389/fphy.2020.00032
[26]	Godunov, S.K. (1959) Finite Difference Methods for the Computation of Discontinu-ous Solutions of the Equations of Fluid Dynamics. Matematicheskii Sbornik, 47, 271-306.
[27]	Dumbser, M. and Munz, C. (2005) Building Blocks for Arbitrary High Order Discontinuous Galerkin Schemes. Journal of Scientific Computing, 27, 215-230. https://doi.org/10.1007/s10915-005-9025-0
[28]	Fambri, F., Dumbser, M. and Zanotti, O. (2017) Space-Time Adaptive ADER-DG Schemes for Dissipative Flows: Compressible Navier-Stokes and Resistive MHD Equations. Computer Physics Communications, 220, 297-318. https://doi.org/10.1016/j.cpc.2017.08.001
[29]	Toro, E.F. (1989) A Weighted Average Flux Method for Hyperbolic Conservation Laws. Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, 423, 401-418. https://doi.org/10.1098/rspa.1989.0062
[30]	Toro, E.F. (1992) The Weighted Average Flux Method Applied to the Euler Equations. Philosophical Transactions. Physical Sciences and Engineering, 341, 499-530. https://doi.org/10.1098/rsta.1992.0113
[31]	Shu, C. (1998) Essentially Non-Oscillatory and Weighted Essentially Non-Oscillatory Schemes for Hyperbolic Conservation Laws. In: Quarteroni, A., Ed., Advanced Numerical Approximation of Nonlinear Hyperbolic Equations, Springer, 325-432. https://doi.org/10.1007/bfb0096355
[32]	Shu, C. (2020) Essentially Non-Oscillatory and Weighted Essentially Non-Oscillatory Schemes. Acta Numerica, 29, 701-762. https://doi.org/10.1017/s0962492920000057
[33]	Titarev, V.A. and Toro, E.F. (2006) Analysis of ADER and ADER-WAF Schemes. IMA Journal of Numerical Analysis, 27, 616-630. https://doi.org/10.1093/imanum/drl033
[34]	Lax, P. and Wendroff, B. (1960) Systems of Conservation Laws. Communications on Pure and Applied Mathematics, 13, 217-237. https://doi.org/10.1002/cpa.3160130205
[35]	Roe, P.-L. and Pike, J. (1984) Efficient Construction and Utilisation of Approximate Riemann Solutions. Proceedings of the Sixth International Symposium on Computing Methods in Applied Sciences and Engineering, Versailles, 12-16 December 1983, 499-518.
[36]	Harten, A., Lax, P.D. and van Leer, B. (1983) On Upstream Differencing and Godunov-Type Schemes for Hyperbolic Conservation Laws. SIAM Review, 25, 35-61. https://doi.org/10.1137/1025002
[37]	Roe, P.-L. (1985) Some Contributions to the Modelling of Discontinuous Flows. Proceedings of the Fifteenth Summer Seminar on Applied Mathematics, La Jolla, 27 June-8 July 1983, 163-193. https://ui.adsabs.harvard.edu/abs/1985ams..conf..163R
[38]	Burden, R.L. and Faires, J.D. (2010) Numerical Analysis. 9th Edition, Brooks/Cole, Cengage Learning.
[39]	Stampolidis, P. and Gousidou-Koutita, M.C. (2024) A Numerical Study of Riemann Problem Solutions for the Homogeneous One-Dimensional Shallow Water Equations. Applied Mathematics, 15, 765-817. https://doi.org/10.4236/am.2024.1511044