{'fraction_non_alphanumeric': 0.04723422042636199, 'fraction_numerical': 0.01978542566531977, 'mean_word_length': 4.415770609318996, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
High-order accurate entropy stable finite difference schemes for the shallow water magnetohydrodynamics

This paper develops the high-order accurate entropy stable (ES) finite difference schemes for the shallow water magnetohydrodynamic (SWMHD) equations.They are built on the numerical approximation of the modified SWMHD equations with the Janhunen source term. First, the second-order accurate well-balanced semi-discrete entropy conservative (EC) schemes are constructed, satisfying the entropy identity for the given convex entropy function and preserving the steady states of the lake at rest (with zero magnetic field). The key is to match both discretizations for the fluxes and the non-flat river bed bottom and Janhunen source terms, and to find the affordable EC fluxes of the second-order EC schemes. Next, by using the second-order EC schemes as building block, high-order accurate well-balanced semi-discrete EC schemes are proposed. Then, the high-order accurate well-balanced semi-discrete ES schemes %satisfying the entropy inequality are derived by adding a suitable dissipation term to the EC scheme with the WENO reconstruction of the scaled entropy variables in order to suppress the numerical oscillations of the EC schemes. After that, the semi-discrete schemes are integrated in time by using the high-order strong stability preserving explicit Runge-Kutta schemes to obtain the fully-discrete high-order well-balanced schemes. The ES property of the Lax-Friedrichs flux is also proved and then the positivity-preserving ES schemes are studied by using the positivity-preserving flux limiter. Finally, extensive numerical tests are conducted to validate the accuracy, the well-balanced, ES and positivity-preserving properties, and the ability to capture discontinuities of our schemes.


Introduction
The shallow water equations are widely used in atmospheric flows, tides, storm surges, river and coastal flows, lake flows, tsunamis, etc. They describe the flow with free surface under the influence of gravity and the bottom topology, where the vertical dimension is much smaller than any typical horizontal scale. Numerical simulation is an effective tool to solve them and a great variety of numerical methods are available in the literature, e.g. [4,31,37,47,48,52,53,54] and the references therein.
Here we are concerned with numerical methods for the shallow water magnetohydrodynamic (SWMHD) equations, which take into account the effect of the magnetic field, originally proposed in [23] for studying the global dynamics of the solar tachocline. The two-dimensional (2D) SWMHD equations with non-flat bottom topography [23,43] read the following quasi-linear hyperbolic balance laws with the divergence-free condition ∇ · (hB) = 0, (1.2) where U is the conservative variables vector, and F 1 and Existing numerical studies for the SWMHD equations include the evolution Galerkin scheme [30], space-time conservation element solution element (CESE) method [42], centralupwind schemes [56], Roe-type schemes [29], second-order entropy stable (ES) finite volume scheme (satisfying the semi-discrete entropy inequality) [50], high-order CESE scheme up to 4th-order [1], etc. Some have dealt with the non-flat topography [50,56], and are wellbalanced in the sense that the schemes can preserve the lake at rest. For the numerical solutions of the SWMHD equations, we need to deal carefully with the divergence-free constraint (1.2). In the ideal magnetohydrodynamic (MHD) case, many works have focused on this issue, for example, the projection method [8], the constrained transport method and its variants [2,17,36,44], the eight-wave formulation of the MHD equations [41], the hyperbolic divergence cleaning method [13], the locally divergence-free DG method [34], the "exactly" divergence-free central DG method [35], and so on. However, if the solutions contain discontinuities, then the above identity does not hold.

Definition 2 (Entropy solution).
A weak solution U of (1.1) is called an entropy solution or a physically relevant solution if for all entropy pairs (η, q ), the entropy inequality or (1.6) holds in the sense of distributions.
The entropy conditions are of great importance in the well-posedness of hyperbolic conservation laws or balance laws, e.g. (1.1), and may improve the robustness of the numerical schemes, thus it is meaningful to seek their numerical schemes satisfying the discrete entropy inequality. For the scalar conservation laws, the conservative monotone schemes are nonlinearly stable and satisfy the discrete entropy conditions so that their solutions can converge to the entropy solution [11,24]; the E-schemes satisfy the semi-discrete entropy condition for any convex entropy [39,40]. However, those schemes are only first-order accurate. Generally, it is very hard to prove that the high-order accurate schemes of the scalar conservation laws and the schemes for the system of hyperbolic conservation laws satisfy the entropy inequality for any convex entropy function. Two relative works are presented in [7] and [27]. The former is second-order accurate and not in the standard finite volume form, while the latter approximates the entropy variables and needs solving nonlinear equations at each time step. In view of this, the researchers usually try to study the high-order accurate entropy conservative (EC) (resp. entropy stable (ES)) schemes, which satisfy the discrete entropy identity (resp. inequality) for a given entropy pair. The second-order EC schemes were built in [45,46], and their higher-order extension was studied in [32]. It is known that EC schemes may become oscillatory near the shock waves, thus additional dissipation terms have to be added into the EC schemes to obtain the ES schemes. Combining the EC flux of the EC schemes with the "sign" property of the ENO reconstructions, the' arbitrary high-order ES schemes were constructed by using high-order dissipation terms [20]. Some ES schemes based on the discontinuous Galerkin (DG) framework were also studied, e.g. the ES space-time DG formulation [3,25] and the ES nodal DG schemes using suitable quadrature rules [10]. The ES schemes based on summation-by-parts (SBP) operators were developed for the Navier-Stokes equations [18]. As a base of those works, the construction of the affordable two-point EC flux is one of the key parts, and has been extended to the shallow water equations [19,22], the MHD equations [9,49], the relativistic (magneto-)hydrodynamic equations [16,15,51], and so on. Those ES MHD schemes are built on the modified MHD equations with the non-conservative source terms, e.g. the Powell source terms [9,14,41] or the Janhunen source terms [28,49]. The Powell source term can be used to obtain the symmetrizable MHD equations, but the solution to the Riemann problem of Powell's MHD equations with initial positive gas pressures may contain a nonphysical negative gas pressure [28], while the Janhunen source term can preserve the conservation of the momentum and energy and be used to restore the positivity of gas pressure. Due to the adding source terms, the sufficient condition proposed in [45] for a finite difference scheme satisfying an entropy identity should be modified, see [9]. One should also take care of the discretization of the source terms, to match it to the discretization of the flux gradients and ensure that the final schemes satisfy the semi-discrete entropy inequality.
The paper aims at constructing the high-order accurate ES finite difference schemes for the SWMHD equations (1.1) for a given entropy pair. With suitable discretization of the non-flat bottom topography and Janhunen source terms in the modified SWMHD equations, a two-point EC flux is derived for constructing the semi-discrete second-order accurate well-balanced EC schemes satisfying the entropy identity. Our discretization of the source terms is essential to achieve both high-order accuracy and well-balance, and does not meet the computational issue in [50] when the magnetic field is zero. The high-order well-balanced EC schemes are constructed by using the above two-point EC schemes as building blocks. In order to avoid the numerical oscillation produced by the EC schemes around the discontinuities, some suitable dissipation terms utilizing the weighted essentially non-oscillatory (WENO) reconstruction in the scaled entropy variables are added to the EC fluxes to get the high-order accurate well-balanced ES schemes satisfying the semi-discrete entropy inequality. The above semi-discrete EC and ES schemes are integrated in time by using the high-order accurate explicit strong-stability preserving Runge-Kutta schemes to obtain the fully-discrete high-order accurate well-balanced schemes. The ES property of the Lax-Friedrichs flux is also proved and then the positivity-preserving ES schemes are developed by using the positivity-preserving flux limiter.
The rest of the paper is organized as follows. Section 2 presents the modified SWMHD equations, the entropy pair and necessary notations. Section 3 constructs the affordable twopoint EC flux and the semi-discrete EC and ES schemes for the 1D SWMHD equations, and proves the well-balanced properties of the EC and ES schemes. The positivity-preserving ES scheme is also studied. Section 4 gives the 2D well-balanced EC and ES schemes. Extensive numerical tests are conducted in Section 5 to validate the effectiveness and performance of our schemes. Section 6 gives some conclusions.

Modified SWMHD equations
This section gives an entropy analysis of the SWMHD equations with the non-flat bottom topography.
If the solutions are smooth, then the SWMHD equations (1.1)-(1.3) can be equivalently cast into the primitive variable form Defining the mathematical entropy as the total energy [19,50] η(U , b) : and using (2.1) gives which means that under the constraint ∇ · (hB) = 0, the following quantities where the vector V = (∂η/∂U ) T is explicitly given by It is not difficult to verify that the matrix ∂U /∂V is symmetric and positive definite.
Similar to the Janhunen source term for the ideal MHD equations [28], one needs to add some non-conservative source terms to get a modified SWMHD system for (1.1)-(1.3) where Ψ = (0, 0, 0, v T ), and satisfies ΨV = Φ(V ) := v · B. Note that the modified SWMHD equations without bottom topography have been discussed in the literature, see e.g. [50].
Taking the dot product of V with (2.4) yields where we have used the identity Notice that the identity (2.5) is obtained without using the divergence-free condition, and will be useful in constructing an entropy stable (ES) scheme because the numerical divergence of ∇ · (hB) may not be zero. Moreover, we define the "entropy potential" ψ from the given which makes the following identity true The "entropy potential" plays an important role in obtaining the sufficient condition for the two-point entropy conservative (EC) fluxes.
Remark 2.1. Notice that the entropy identity (2.5) is slightly different from (1.5), since the source terms in the SWMHD equations have special structure, and then V T G(U ) in (1.5) can be absorbed into the left-hand side by using In other words, the difference between the entropy pair (η, q ) in (2.5) and that in (1.5) is (ghb, ghbv ).

Second-order EC schemes
Let us consider a uniform mesh in x: x 1 < x 2 < · · · < x Nx , with the spatial step size ∆x = x i − x i−1 , i = 2, · · · , N x and the semi-discrete conservative finite difference scheme for (3.1) as follows , and the second-order central difference is used to approximate ∂(hB 1 )/∂x and ∂b/∂x in the source terms. satisfies a semi-discrete entropy identity for some numerical entropy flux q i+ 1 2 consistent with the given physical entropy flux q 1 .
The following lemma gives a sufficient condition for the semi-discrete scheme (3.2) to be EC, with the discretization of the source terms.
is used in (3.2), where a and { {a} } denote the jump and mean of a, respectively, then the semi-discrete scheme (3.2) is second-order accurate and EC, with the numerical entropy flux has been used in the third equality. Thus the scheme (3.2) with F i+ 1 is EC in the sense of The discretization of the source terms is second-order accurate since the second-order central difference is used, and the results in [45] show that the discretization of the flux gradient is second-order accurate, therefore the the scheme (3.2) with F i+ 1 Remark 3.1. Since the central difference is used for approximating the non-flat river bed bottom and Janhunen source terms, the sufficient condition (3.4) is different from that in [19,50]. Moreover, such discretization of the source terms is essential to achieve high-order accuracy and well-balance, see the subsection 3.2.
Below we present such EC flux satisfying (3.4).
is an EC flux, consistent with the physical flux Proof. The key is to use the identity and rewrite the jumps of the entropy variables V , the potential ψ 1 , Φ, hv 1 and hbv 1 as some linear combinations of the jumps of a specially chosen parameter vector. For simplicity in derivation, the parameter vector is chosen as (h, b, v 1 , v 2 , B 1 , B 2 ), then Substituting them into (3.4), and equating the coefficients of the same jump terms gives the numerical flux in (3.6).
If letting U i = U i+1 , it is easy to see the consistency of the numerical flux in (3.6).

Remark 3.2.
For the y-split system of (2.4), the rotational invariance may be used to get the EC fluxes consistent to F 2 (U ) defined in (1.3). The EC flux (3.6) is the same as the one in [50].
Remark 3.3. The present discretization of the source terms is different from that in [50], and does not meet the computational issue in [50] Remark 3.4. If the magnetic field B ≡ 0, then the SWMHD equations reduce to the shallow water equations (SWEs) and the above SWMHD scheme defined in (3.6) reduces to the well-balanced EC scheme for the SWEs with the EC flux which is the same as the EC flux in [19], except for the second component due to the different approximation of the source terms.
Lemma 3.1 and Theorem 3.2 tell us that the SWMHD scheme is second-order accurate and EC. Moreover, we can show that it is well-balanced.
Theorem 3.3. The scheme (3.2) with the EC flux (3.6) is well-balanced, in the sense that when the magnetic field is zero, it preserves the lake at rest, that is to say, for the given initial data the solutions of (3.2) satisfy Proof. Under the hypotheses, one can verify that the scheme (3.2) satisfies Therefore the scheme (3.2) preserves the lake at rest.
Remark 3.5. The second-order EC scheme (3.2) satisfies the 1D moving equilibrium state [19] In fact, it is easy to verify that

High-order EC schemes
To develop the high-order well-balanced EC schemes, our task is to get the high-order numerical fluxes and conduct the matched high-order discretization of the source term related to the non-flat river bed bottom and the Janhunen source term in (3.1).
Following the way in [32], the EC flux of the 2pth-order (p ∈ N + ) accurate scheme can be obtained by using the linear combinations of the "second-order accurate" EC flux (3.6) as follows F 2pth The readers are referred to [20,32] for more details on constructing the "high-order accurate" EC flux.
To make the resulting schemes high-order accurate, well-balanced and EC, it is essential that the high-order finite difference approximations of the spatial derivatives ∂(hB 1 )/∂x and ∂b/∂x in the source terms should match the "high-order accurate" EC flux (3.8). To this end, based on the observation that the second-order central differences for the source terms used in the second-order EC scheme (3.2), have the same form as the discretization of the flux gradient, using those second-order central differences as a building block can obtain the high-order approximations of the source terms as follows ( hB 1 ) 2pth where the linear combination coefficients are the same as those in the "high-order accurate EC flux" (3.8). Similarly, it is not difficult to verify Such treatment can also be found in [14].
In summary, by approximating (3.1), we obtain the following 2pth order semi-discrete EC scheme which satisfies the entropy identity (3.3) with the numerical entropy flux ( q) 2pth It is a linear combination of the two-point numerical entropy flux (3.5). For example, when p = 3, the expression of the "6th-order accurate" EC flux is explicitly given as follows F 6th It can also be verified that the scheme (3.9) is well-balanced in the sense of Theorem 3.3, since the numerical fluxes and the numerical source terms in (3.9) are formed by the linear combinations of the fluxes and the source terms in the second-order scheme (3.2) with the same coefficients, specifically, the second equation in (3.9) is written as follows For the 1D moving equilibrium states discussed in Remark 3.5, one needs to impose very restrictive conditions

ES schemes
It is known that for the quasi-linear hyperbolic conservation laws or balance laws, the entropy identity is available only if the solution is smooth. In other words, the entropy is not conserved if the discontinuities such as the shock waves appear in the solution. Moreover, the EC scheme may produce serious unphysical oscillations near the discontinuities. Those motivate us to develop the ES scheme in this section by adding a suitable dissipation term to the EC scheme to avoid the unphysical oscillations produced by the EC scheme and to satisfy the entropy inequality for the given entropy pair.
Following [45], adding a dissipation term to the EC flux F i+ 1 2 gives the ES flux where D i+ 1 2 is a symmetric positive semi-definite matrix. It is easy to prove that the scheme Motivated by the Cholesky decomposition and the dissipation term in the (local) Lax- in (3.11) can be chosen as Here and RR T is the Cholesky decomposi- .
To obtain the arbitrary high-order accurate ES scheme, the dissipation term in (3.11) has to be improved. , and then define Combining such reconstructed jump with the "2pth-order EC flux"F 2pth and the 2pth-order discretization of the source terms gives the kth order ES scheme as follows where p = k/2 for even k and p = (k + 1)/2 for odd k, and F kth (3.14) The semi-discrete numerical schemes (3.13) is ES if the reconstruction satisfies the "sign" property [20] sign( w i+ 1 2 ) = sign( w i+ 1 2 ).
which does hold for the ENO reconstructions [21].
Moreover, one can also obtain higher-order accurate ES scheme with the WENO reconstruction instead of the ENO reconstruction, if the same number of candidate points values are used. In view of that a general WENO reconstruction may not satisfy the "sign" property, following [5], the dissipation term in (3.11) may be modified as follows is a switch function defined by Remark 3.6. At the steady state, the entropy variables V T become (gC, 0, 0, 0, 0), so that the low-order dissipation term with V and the high-order dissipation term with w = R V all vanish. Thus the constructed ES schemes are well-balanced.

Time discretization
This paper uses the following third-order accurate strong stability preserving explicit Runge-Kutta (RK3) time discretization to their right-hand side.

Positivity-preserving ES schemes
This section restricts to the flat bottom topography. Generally, the high-order ES scheme where h( F ) denotes the first component of the vector F .
Before that, we first prove that the scheme (3.2) with the local Lax-Friedrichs flux is positivity-preserving. is positivity-preserving under the CFL condition Proof. The first component of the Lax-Friedrichs scheme can be split as Thus h ±,LF i > 0 under the CFL condition (3.19), and then h n+1 i > 0.
The following lemma shows that the Lax-Friedrichs flux (3.18) is ES even if hB 1 or hB 2 is not constant. Its proof is direct without the assumption that the 1D exact Riemann solution of x-split system is ES. In fact, when there are jumps in hB 1 or hB 2 at the cell interface in the 2D case, whether such assumption is available needs further investigation.

Lemma 3.5. The Lax-Friedrichs flux is ES when the bottom topography is flat.
Proof. Substituting the Lax-Friedrichs flux into the inequality (3.12) and using identity (3.7) gives If hB 1 0, then B 0, because Similarly, C 0. Therefore the Lax-Friedrichs flux is ES.
Based on those discussions, the high-order positivity-preserving ES schemes can be constructed by using the Lax-Friedrichs flux and the positivity-preserving limiter [26]. The limited numerical flux F kth,PP i+ 1 2 is given by is the scaling factor corresponding to the two neighboring grid points, which share the same flux F kth i+ 1 2 , and It is worth noting that the discretization of the source terms should also be replaced by the corresponding convex combinations as follows It is easy to verify that the water height updated by the high-order positivity-preserving ES schemes satisfies

Two-dimensional schemes
This section extends the 1D high-order EC and ES schemes developed in Section 3 to the 2D SWMHD system (1.1). For convenience, the notation (x 1 , x 2 ) is replaced with (x, y). Our attention is limited to a uniform Cartesian mesh {(x i , y j ), i = 1, · · · , N x , j = 1, · · · , N y } with the spatial step sizes ∆x, ∆y so that the extension of the 1D schemes to (1.1) can be done by approximating (2.4) in a dimension by dimension fashion. To avoid repetition, the detailed extension is not described below.
Similarly, using the EC scheme (4.1) as building block can give a 2pth-order well-balanced semi-discrete EC scheme for the 2D SWMHD system (1.1) as follows Then adding a suitable dissipation term to (4.2) gives a kth-order well-balanced semi-discrete ES scheme for the 2D SWMHD system (1.1) as follows where p = k/2 for even k and p = (k + 1)/2 for odd k,

Numerical results
This section conducts some numerical experiments to validate the performance of our EC and ES schemes for the SWMHD equations (1.1) and its x-split system. Unless otherwise stated, all computations take the CFL number µ as 0.5, and the 5th-order schemes with the fifth-order WENO reconstruction in [6]. For the accuracy tests, the time stepsize ∆t is taken as µ∆x 6/3 (resp. µ∆x 5/3 ) for the 6th order EC schemes (resp. the 5th order ES schemes) to make the spatial error dominant.

One-dimensional case
Example 5.1 (Accuracy test). This example is used to verify the accuracy. The computational domain is [0, 1] with periodic boundary conditions, and g = 1. The constructed exact solution is given as follows   The surface level h + b and the bottom b are shown in Figure 5.1, and the errors in h and v 1 are given in Table 5.2. It can be seen that the errors are at the level of round-off errors for the double precision, and the well-balanced property is verified.   [56]). This example is adapted from the problem in [55] and used to check the dissipative and dispersive errors in the ES scheme.
The computational domain and the bottom topography are the same as the last problem,  g = 9.812, and the initial data are The results obtained by using the ES scheme with N x = 50, 100 are shown in Figure 5.2, and the reference solutions are obtained by using the ES scheme with N x = 1000. We can see that the accurate solutions can be obtained even with the coarse mesh N x = 50.     The initial discontinuity will be decomposed into two magnetogravity waves and two Alfvén waves propagating away in two directions as the time increases. The problem is solved until t = 0.4 with g = 1.   with v max = 0.2, B max = 0.1, r = x 2 + y 2 . Using the Galilean transformation x = x − t, y = y − t, t = t can give a time-dependent exact solution

Two-dimensional case
which describes a vortex moving with a constant speed (1,1).
The computational domain is [−8, 8] 2 with periodic boundary conditions, g = 1, h max = 1, and the output time is t = 16 so that the vortex travels and returns to the original position after a period. Table 5.3 lists the errors in h and corresponding orders of convergence. The results show that our EC and ES schemes achieve the optimal convergence order. Figure 5.6 plots the contours of h and the magnitude of the magnetic field |B| with 40 equally spaced contour lines. It can be seen that our schemes can preserve the shape of the vortex well after a whole period.
To check the positivity-preserving property, let us do a test with which implies that the lowest water height is 10 −6 . The errors and orders of convergence obtained by using the ES scheme with or without the positivity-preserving (PP) limiter are listed in Table 5 with zero velocity and magnetic field. Outflow boundary conditions are used and the problem is solved until t = 0.6 with g = 9.812. which describe a right-going disturbance as it propagates past the hump. It can be seen that the complex small features are well captured without any spurious oscillations, and the results are comparable to those in [54].
Example 5.9 (Orszag-Tang like problem [56]). It is similar to the Orszag-Tang problem for the ideal MHD equations [38]. The computational domain is [0, 2π] 2 with periodic boundary The solution of this problem is smooth initially, but the complicated pattern will arise as the time increases and it has the turbulence behavior. Figure 5.8 presents the results obtained by using the ES scheme with N x = N y = 200 at t = 1 and 2. The solutions are in good agreement with those in [56]. The left plot in Figure 5.10 shows the time evolution of the discrete total entropy i,j η(U i,j )∆x∆y with three spatial resolutions N x = N y = 100, 200 and 400. The total entropy is conserved when the solutions are smooth during an initial period, but it decays when discontinuities arise.  Example 5.10 (Rotor like problem [30]). It is an extension of the classical ideal MHD rotor test problem [30]. The computational domain is [−1, 1] 2 with outflow conditions with g = 1.
Initially hB 1 = 1 and hB 2 = 0, and there is a disk of radius r 0 = 0.1 centered at (0, 0), where fluid with large h is rotating in the anti-clockwise direction. The ambient fluid is homogeneous for r > r 0 , where r = x 2 + y 2 . Specifically, the initial data are  without obvious spurious oscillations comparable to those in [30]. The right plot of Figure   5.10 displays the time evolution of the discrete total entropy with three spatial resolutions N x = N y = 100, 200 and 400. The results show that the total entropy decays as expected, and the fully discrete scheme is also ES.

Conclusion
The paper proposed the high-order accurate entropy stable (ES) finite difference schemes for the one-and two-dimensional shallow water magnetohydrodynamic (SWMHD) equations with non-flat bottom topography. The Janhunen source term was added to the conservative SWMHD equations. For the modified SWMHD equations, the second-order accurate well-balanced semi-discrete entropy conservative (EC) finite difference scheme was first constructed, in the sense that it satisfied the semi-discrete entropy identity for the given entropy pair (the total energy served as the mathematical entropy) and preserved the steady state of lake at rest (with zero magnetic field). The key was to match both discretizations for the fluxes and the source term related to the non-flat river bed bottom and the Janhunen source term, and to find the affordable EC fluxes of the second-order EC schemes. Next, the highorder accurate well-balanced EC schemes were obtained by using the second-order accurate EC schemes as building block. In view of that the EC schemes might become oscillatory near the discontinuities, the appropriate dissipation terms were added to the EC fluxes to develop the semi-discrete well-balanced ES schemes satisfying the semi-discrete entropy inequality.
The WENO reconstruction of the scaled entropy variables and the high-order explicit Runge-Kutta time discretization were implemented to obtain the fully-discrete high-order schemes.
The ES and positivity-preserving properties of the Lax-Friedrichs scheme were also proved without the assumption that the 1D exact Riemann solution of x-split system was ES, and then the high-order positivity-preserving ES schemes were developed by using the positivitypreserving flux limiter. Extensive numerical tests showed that our schemes could achieve the designed accuracy, were well-balanced or positivity-preserving, and could well capture the discontinuities.

Appendix A. SWMHD in polar coordinates
In the polar coordinate system (r, θ), the modified SWMHD equations (2.4) without the bottom topography become ∂h ∂t + ∇ · (hv) = 0, and F r , F θ are the radius and azimuth components of the vector F , respectively.