WELL-BALANCED POSITIVITY PRESERVING ADAPTIVE MOVING MESH CENTRAL-UPWIND SCHEMES FOR THE SAINT-VENANT SYSTEM

. We extend the adaptive moving mesh (AMM) central-upwind schemes recently proposed in Kurganov et al . [ Commun. Appl. Math. Comput. 3 (2021) 445–479] in the context of one-(1-D) and two-dimensional (2-D) Euler equations of gas dynamics and granular hydrodynamics, to the 1-D and 2-D Saint-Venant system of shallow water equations. When the bottom topography is nonflat, these equations form hyperbolic systems of balance laws, for which a good numerical method should be capable of preserving a delicate balance between the flux and source terms as well as preserving the nonnegativity of water depth even in the presence of dry or almost dry regions. Therefore, in order to extend the AMM central-upwind schemes to the Saint-Venant systems, we develop special positivity preserving reconstruction and evolution steps of the AMM algorithms as well as special corrections of the solution projection step in (almost) dry areas. At the same time, we enforce the moving mesh to be structured even in the case of complicated 2-D computational domains. We test the designed method on a number of 1-D and 2-D examples that demonstrate robustness and high resolution of the proposed numerical approach.


Introduction
We consider the Saint-Venant system of shallow water equations, which was first introduced in [10] and is widely used to model water flow in rivers, canals and coastal areas as well as in atmospheric sciences and oceanography.In the one-dimensional (1-D) case the studied system reads as where ℎ(, ) is the water depth, (, ) is the velocity, () is the bottom topography, and  is the constant gravitational acceleration.The two-dimensional (2-D) Saint-Venant system is where ℎ(, , ) is the water depth, (, , ) and (, , ) are the -and -velocities, and (, ) is the bottom topography.Development of accurate, efficient and robust numerical methods for the systems (1.1) and (1.2) is an important and challenging problem due to several reasons.First, these systems admit nonsmooth solution, which, in the case of discontinuous bottom topography , may not be unique.Second, it is very important to preserve a delicate balance between the flux and source terms since many practically relevant solutions are, in fact, small perturbations of the so-called "lake at rest" steady states: where  is an equilibrium water surface variable (we say that a numerical method is well-balanced if it is capable of exactly preserving the "lake at rest" states).Third, in many practically important situations, one may need to deal with dry or almost dry areas and then it is crucial for the developed numerical method to preserve the nonnegativity of ℎ.
In the past decades, many well-balanced and positivity preserving numerical methods for the Saint-Venant systems have been developed; see, e.g., the review papers [16,24] and references therein.In this paper, we focus on semi-discrete central-upwind schemes for (1.1) and (1.2), which were designed on a variety of fixed grids: uniform Cartesian [5,17,18], unstructured triangular [7], quadrilateral [21], and cell-vertex polygonal [3] ones.We follow the lines of [19] and construct an adaptive moving mesh (AMM) central-upwind scheme for the Saint-Venant systems (1.1) and (1.2).The AMM central-upwind schemes, which have been recently introduced in [19] for 1-D and 2-D hyperbolic systems of PDEs, are based on structural meshes, which are evolved in time according to the movimg mesh differential equations; we refer the reader to [1,8,9,13] for several examples of existing AMM algoritms.Our goal is to develop several techniques required to ensure that the resulting scheme is well-balanced and positivity preserving.In order to achieve this goal, we first generalize the 1-D well-balanced positivity preserving semi-discrete central-upwind scheme from [17] for the 1-D nonuniform grids and modify the 2-D well-balanced positivity preserving semi-discrete central-upwind scheme from [21] to evolve the solution over the structured quadrilateral meshes.In these schemes, second-order well-balanced quadratures for the geometric source terms are developed to ensure the well-balanced property.In order to preserve the positivity of the computed water depth, several measures are taken.First, we either make sure that the reconstructed values of the water surface stay above the corresponding values of the bottom topography or use the positivity preserving reconstructions for the water depth.Second, we use a draining time-step technique, originally proposed in [4], to ensure that the water depth remains positive during the evolution step.Third, we propose special corrections of the solution projection step in (almost) dry areas.We stress that the positivity preserving technique we develop here does not allow us to ensure the well-balanced property of the proposed AMM central-upwind schemes in the presence of dry areas.An alternative positivity preserving approach was proposed in [25,26], where a GRP AMM method on unstructured triangular meshes was introduced.In addition to preserving positivity of the water depth, this method is capable of ensuring the total water conservation, which is enforced by redistributing the conservation error.Such redistribution, however, may lead to a purely artificial appearance of water in dry areas and we prefer not to implement this technique in our AMM method.Finally, we stress that development of a well-balanced AMM method that can accurately handle wetting/drying interfaces and preserve both "lake at rest" and "dry lake" steady states still remains an open problem.
The paper is organized as follows.The proposed AMM algorithms are presented in Section 2. First, in Section 2.1, we introduce the second-order well-balanced and positivity preserving central-upwind scheme on 1-D nonuniform grids for the system (1.1).Then, in Section 2.2, we discuss the 1-D moving mesh equation and its discretization.We complete the description of the 1-D AMM algorithm in Section 2.3, where we introduce the 1-D conservative positivity preserving projection.In Section 2.4, we present the second-order well-balanced and positivity preserving central-upwind scheme on 2-D structured quadrilateral meshes for the system (1.2).We then design the 2-D AMM algorithm with the help of the 2-D moving mesh equation, which is discussed together with its discretization in Section 2.5, and a special 2-D conservative positivity preserving projection, which is introduced in Section 2.6.Finally, in Section 3, we demonstrate the performance of the proposed AMM algorithms through several 1-D and 2-D numerical examples.

Adaptive Moving Mesh central-upwind schemes
In this section, we present the AMM central-upwind schemes for the systems (1.1) and (1.2).The schemes are based on the 1-D and 2-D AMM central-upwind schemes introduced in [19].They are constructed in two steps.Given the solution at a certain time level, it is first evolved to the new time level on a given mesh (see Sects.2.1 and 2.4), which is updated at the end of the evolution step by solving the moving mesh differential equation (see Sects.2.2 and 2.5).The solution is then projected in a conservative manner to the new finitevolume mesh (see Sects.2.3 and 2.6).
We use the same notation as in [19], which is briefly reviewed in this paper for the sake of completeness.

One-dimensional semi-discrete scheme
We first rewrite the system (1.1) in terms of the water depth ℎ and discharge  := ℎ: This system can be put into the vector form with )︃ , ( , ) := We then apply the 1-D AMM central-upwind scheme from [19] to the system (2.1), (2.2).As mentioned in Section 1, we will need to develop the scheme that preserves the positivity of ℎ and is well-balanced in the sense that it is capable of exactly preserving the "lake at rest" steady states ( ≡ Const,  ≡ 0).
Assume that the computational domain is covered with nonuniform cells , and that at a certain time , the cell averages of the computed solution, are available.They are then evolved from time level  to  + ∆ using the semi-discrete central-upwind scheme on nonuniform grids from Section 2.1 of [19]: where the numerical fluxes are given by Here,  ± ), and the second component of the source term is approximated using the well-balanced quadrature developed in [17]: In (2.4)Note that in (2.4), (2.5), all of the indexed quantities depend on , but from now on, we will omit this dependence for the sake of brevity.

Positivity preserving reconstruction
We begin by following Section 2.1 of [19] and obtain the generalized minmod piecewise linear reconstruction of the water surface  := ℎ +  and discharge  := ℎ, which are used to evaluate  ± + 1 2 and  ± + 1

2
. Unfortunately, for some  the obtained point values  ± + 1 2 may be smaller than the corresponding value  + 1  2 , which would lead to negative point values of the water depth ℎ ± 2 .We therefore follow [17] and reconstruct ℎ instead of  in potentially dry areas.
In this paper, the cell   is called "dry" if at least one of the following inequalities is satisfied: where  > 0 is a small parameter (we take  = 10 −16 in all of the numerical experiments reported in Sect.3) and   = 1 Δ ∫︀  () , which should be numerically computed using a quadrature (in our numerical examples, we have used either the trapezoidal or Simpson's rule).All other cells are called "wet".
In "dry" cells, we use the piecewise linear positivity preserving reconstruction (described in [19], Sect.2.1) of ℎ based on the cell averages of water depth ℎ  :=   −  .After computing ℎ − +  + 1 2 there.Remark 2.1.In practice, one may replace the condition (2.6) with its simplified version: which is much easier to check.We stress that the use of the simplified condition (2.8) does not guarantee positivity of the projected values of ℎ (see the proof of Theorem 2.5 below), but negative ℎ may appear only in a small number of cells located near the wet/dry interfaces.If this occurs, we reclassify these particular cells from "wet" to "dry" and reconstruct ℎ instead of  in them.

Desingularization and one-sided local speeds
In order to avoid division by zero (or by a very small number), we follow [18] and desingularize the computation of the velocity point values needed to evaluate numerical fluxes in (2.4) by setting and then for consistency we modify the corresponding values of the discharge by recalculating In (2.9),  is a small desingularization parameter, which we take to be equal to  in (2.7).
Equipped with the values of  ± , we estimate the one-sided local speeds of propagation needed in (2.4) as follows: , 0 }︁ .

Positivity preserving evolution
In order to guarantee the positivity of ℎ during the evolution step, a draining time-step technique introduced in [4] is employed.
We first consider the forward Euler discretization of (2.3), the first component of which reads as where ∆ is the time step constrained by the CFL condition: We then denote by and replace (2.10) by where the time step ∆ + 1 2 is defined as: The corresponding forward Euler step for   is where the advective and gravitational parts of the fluxes are defined by It has been shown in [5] that in the case of uniform grid the resulting central-upwind scheme is both positivity preserving and well-balanced.The proof can be directly extended to the case of a nonuniform grid and the three-stage third-order strong stability preserving (SSP) Runge-Kutta method (see, e.g., [11,12]) as it can be written as a convex combination of forward Euler steps.
Remark 2.2.An alternative way to ensure the positivity of ℎ during the time evolution is to use a more restrictive CFL condition, which we derive in Appendix C.1.This approach, however, will significantly affect the efficiency of the overall method as the size of time steps is to be reduced by a factor of about 3. We have carefully compared the numerical results obtained by both approaches and realized that using much smaller CFL number does not lead to any improvement in the quality of the computed solutions.Therefore, in all of the numerical examples reported in Section 3, we have used the draining time-step technique.

One-dimensional moving mesh equation
In this section, we briefly describe the 1-D moving mesh equation and its numerical solution algorithm; see Section 3.1 of [19] for details.
In addition to the computational domain [, ] covered by the nonuniform mesh

}︁
, we introduce the uniform logical mesh  + 1 2 = ∆,  = 0, . . ., , with ∆ = 1/ .Let us denote the one-to-one coordinate transformation from the logical domain to the computational one by The mesh is distributed according to the following moving mesh equation (see, e.g., [13] for a detailed derivation): (  )  = 0, ( where  is a monitor function and  is a differential operator (see, e.g., [2,13,23]).In this paper, we use (for some component of  ), which is approximated using the second-order centered difference: The function  in (2.11) is a smoothing filter designed as follows.We first compute  0  = |  |, and then smooth  0  out by averaging over the neighboring cells for each  for a prescribed number of iterations, that is, we introduce and then set ((| |))  :=    , which is used in (2.11).In our numerical experiments, we have taken  = 4. Finally,  in (2.11) is the intensity parameter employed to control the mesh concentration: the use of larger values of  leads to the higher concentration of grid points in the "rough" areas.We follow [14] and choose  to be where  ∈ (0, 1) is a prescribed fraction of mesh points to be concentrated in the "rough" areas of the computed solution.
Equipped with the monitor function , we discretize the moving mesh equation (2.11) using the centered difference approximation, which results in the following linear algebraic system for the mesh points locations: + 1 2 = .We numerically solve this system using the Jacobi iterations combined with the mesh relaxation procedure needed to ensure that the mesh does not get distracted during the iterations.Denoting by   + 1   2   the grid nodes in the beginning of the ( + 1)-th iteration step (with the initial guess  0 + 1 2 being the grid nodes from the previous evolution step), we take one Jacobi sweep where

}︁
, which is then relaxed by setting which implies 2 and thus  +1 )︁ , which means that the logical structure of the mesh does not change.
Remark 2.3.One can set up a stopping criterion for the iteration (2.12), for instance, the iterations should stop as soon as max }︁ <  for a given  > 0. In practice, however, the mesh does not typically change much in one evolution time step.We therefore improve the efficiency of the resulting AMM method by stopping the iteration process after several iterations.In all of our numerical experiments, the upper bound on the number of iterations has been set to 4.
Remark 2.4.In order to avoid the appearances of excessively small cells, which would lead to severe time step restrictions, we modify the mesh movement as follows.If , where ∆ min is the smallest allowed cell size, then we locally freeze the movement of the mesh by setting , and in the next iteration of moving mesh process, we set  0 −1 = 0,  0  = 0 and  0 +1 = 0 to reduce the mesh concentration nearby.In all of the numerical examples, we have set ∆ min = 0.1∆ unif , where ∆ unif is the size of cells of the corresponding uniform mesh, which contains precisely the same total number of cells as used in the AMM computations.

One-dimensional conservative positivity preserving projection
After obtaining the new mesh, we project the solution from the cells ]︁ to the new cells ,  +1 + . The conservative solution projection step from [22] is given by where and are the point values reconstructed over the grid    as described in Section 2.1.1.In particular, this projection step for the water surface variable  reads as (2.15) We note that this projection step preserves "lake at rest" steady states, for which =  for all , and thus (2.15) reduces to We also note that the total amount of  after completion of this projection step is conserved, but the total amount of ℎ is not conserved, since, in general, where are approximations of the cell averages of , computed by a proper numerical quadrature.As one can expect, the use of a higher-order quadrature may help to reduce the conservation error.In the numerical experiments reported in Section 3, we have compared the results obtained using the trapezoidal and Simpson's rules and clearly seen the advantage of the higher-order quadrature.The projection step (2.15), however, cannot be used in the "dry" areas, since one of the following may happen: -appearance of negative water depth, that is, for some ; -artificial propagation of water into dry areas, that is, for  located in the totally dry area (far away from the wetting/drying interface), where  +1  must be equal to  +1  .
We therefore replace the projection step (2.15) with where Here,   + 1 2 and ℎ  + 1 2 are determined using (2.14) and is an approximation of the average of  over the interval ]︁ depending on whether . Finally, the "dry" cell correction term  corr  is defined as We note that the magnitude of  corr  is determined by the accuracy of the quadratures used to approximate the integrals in (2.16) and (2.18).Indeed, assuming that these integrals are evaluated using a quadrature of order , we have

}︁
. We now state the positivity preserving property of the described projection step.
Theorem 2.5.The projection step (2.17)-(2.19)preserves the positivity of ℎ  , namely, provided    ≥    , ∀ and the integrals in (2.16) and (2.18) are evaluated exactly.The proof of this important theorem is quite technical and provided in Appendix A.
Remark 2.6.The proof of Theorem 2.5 is still true if the integrals in (2.16) are computed by a numerical quadrature with nonnegative weights (for example, the trapezoidal and Simpson's rules, which have been used in our numerical experiments).We only need to assume that the integral in (2.18) is computed exactly.However, in our numerical experiments, we have used either the trapezoidal or Simpson's rule and  +1  −  +1  always remained positive.

Two-dimensional semi-discrete scheme
Similarly to the 1-D case, we rewrite the system (1.2) in terms of water depth ℎ and discharges   := ℎ and   := ℎ: This system can be put into the vector form with  = (ℎ,   ,   ) ⊤ and (2.22) Figure 1.A typical quadrilateral cell  , with its four neighbors.
We apply the 2-D AMM central-upwind scheme from [19] to the system (2.21), (2.22).As in the 1-D case, we will need to develop the scheme that preserves the positivity of ℎ and is well-balanced in the sense that it is capable of exactly preserving the "lake at rest" steady states ( = Const,   ≡   ≡ 0).
Assume that the computational domain is covered with a structured irregular quadrilateral mesh consisting of cells  , of size | , |, and use the following notations (see Fig. 1): 2 )︁)︁ : the unit outer normal vector to the edge )︁)︁ : the unit outer normal vector to the edge Assume that at a certain time , we have computed an approximate solution, realized in terms of its cell averages: They are then evolved from time level  to  + ∆ using the semi-discrete central-upwind scheme on structured quadrilateral grids from Section 2.2 of [19]: where the numerical fluxes along the boundaries of  , are given by and the second and third components of the cell averages of the source term  , are to be approximated in the well-balance manner; see Section 2.4.

Well-balanced source term quadratures
In order to derive the desired second-order well-balanced quadrature for the source term  , , we follow [3,7,21] and first use Green's formula together with the mean-value theorem to obtain ]︁ . (2.25) The line integrals in the terms I , and II , are then approximated using the midpoint rule, which results in (2.26) and Finally, the cell average of the bottom topography, (2.28) should be numerically computed using a quadrature (in our numerical examples, we have split the quadrilateral  , into four triangles and used 7-point Gaussian quadrature in each of the triangles; see Appendix B for details).
The second term in the numerical source is and it is equal to zero since Thus, we have proved that d d (2) Similarly, one can show that d d , = 0 at "lake at rest" steady states, which completes the proof of the theorem.

Positivity preserving reconstruction
As in the 1-D case, we begin with applying the generalized minmod piecewise linear reconstruction, which has been extended to structured quadrilateral meshes in Section 2.2 of [19], to evaluate  ±  .In order to preserve the positivity of the water depth, we implement the same strategy as in the 1-D case and reconstruct ℎ in the potentially dry areas.
Similarly to the 1-D case, the cell  , is called "dry" if the amount of water there is very small, namely, if at least one of the following inequalities is satisfied: (2.33) where ,  > 0 is the same small parameter as in (2.7), and  , is given by (2.28).All other cells are called "wet".
As in the 1-D case, in "dry" cells, we use the piecewise linear positivity preserving reconstruction (described in [19], Sect.2.2) of ℎ based on the cell averages of water depth ℎ , :=  , − , .After computing ℎ − +  ,+ 1 2 there.Remark 2.8.In practice, one may replace the condition (2.33) with its simplified version: which is much easier to check.As in the 1-D casse, the use of the simplified condition (2.35) does not guarantee positivity of the projected values of ℎ, but negative ℎ may appear only in a small number of cells located near the wet/dry interfaces.If this occurs, we reclassify these particular cells from "wet" to "dry" and reconstruct ℎ instead of  in them.

Desingularization and directional local speeds
Similarly to the 1-D case, we follow [18] and desingularize the computation of the velocity point values needed to evaluate numerical fluxes by setting )︀ and then for consistency we modify the corresponding values of the discharges by recalculating In (2.36),  is a small desingularization parameter, which we take to be equal to  in (2.34).
Equipped with the values ℎ ± , ,  ± , , and  ± , , we estimate the directional local speeds of propagation needed in (2.24) as follows: where  is the velocity in the direction normal to the corresponding side of  , , namely, )︀ . (2.37)

Positivity preserving evolution
As in the 1-D case, we enforce the positivity of ℎ during the evolution step by implementing the draining time-step technique from [4].
We first consider the forward Euler discretization of (2.23), the first component of which reads as where ∆ is the time step constrained by the CFL condition: We then denote by ,+ 2 )︁ , and replace (2.38) with where ∆ + 1 2 , and ∆ ,+ 1 2 are defined as follows: 2 The corresponding forward Euler steps for   and   are , , where the advective and gravitational parts of the fluxes can be obtained by Remark 2.9.Similarly to the 1-D case, one can alternatively use a more restrictive CFL condition (derived in Appendix C.2) to ensure the positivity of ℎ during the time evolution, but this will slow down the computations by a factor of about 2. At the same time, as in the 1-D case the use of smaller CFL number does not lead to any improvement in the quality of the computed solutions.In all of the numerical examples reported in Section 3, we have used the draining time-step technique.

Two-dimensional moving mesh equation
In this section, we briefly describe the 2-D MMPDE and its numerical solution algorithm; see Section 3.2 of [19] for details.
Assume that the computational domain [, ]×[, ] is covered by the nonuniform mesh

}︁
. We introduce the uniform rectangular logical mesh where ∆ = 1/ and ∆ = 1/ are the spatial scales in the -and -directions, respectively.Let us denote the one-to-one coordinate transformation from the logical domain to the computational one by ).We assume that (0, ) =  and (1, ) =  for all  as well as (, 0) =  and (, 1) =  for all .
The mesh is distributed according to the following MMPDE: where  := (, ),  is a monitor function, and  is a differential operator.In this paper, we use  =  (for some component of  ), which is approximated using the second-order centered differences: The function  in (2.39) is a smoothing filter designed as follows.We first compute  0 , := | , | and then smooth  0 , out by averaging the values over the neighboring cells for each ,  for a prescribed number of iterations, that is, we introduce )︀ , ℓ = 0, . . .,  − 1, and then set ((‖ ‖)) , :=   , , which is used in (2.39).In our numerical experiments, we have taken  = 4. Finally,  in (2.39) is an intensity parameter needed to control the mesh concentration.In our computation, we choose  to be where  ∈ (0, 1) is the prescribed fraction of mesh points to be concentrated at the "rough" areas of the computed solution.
Equipped with the monitor function , we discretize the MMPDEs (2.39) using the centered difference approximation, which results in the following linear algebraic system for the mesh points locations: where  ,+ 1 2 := ( , +  ,+1 )/2 and  + 1 2 , := ( , +  +1, )/2.Similarly to the 1-D case, we numerically solve this system using the Jacobi iterations combined with the mesh relaxation procedure to avoid rapid change of mesh.Denoting by   + 1  2 ,+ 1 2 the grid nodes in the beginning of the ( + 1)th iteration step (with the initial guess  0 being the grid nodes from the previous evolution step), we take one Jacobi sweep where   ,+ )︁}︁ , which is then relaxed by setting )︁ ,  = 1, . . .,  − 1,  = 1, . . .,  − 1. (2.40) Remark 2.10.Similarly to 1-D case, we stop the iteration process in (2.40) after a fixed number of iterations.In all of our numerical experiments, the upper bound on the number of iterations has been set to 4.

Two-dimensional conservative positivity preserving projection
After obtaining the new mesh, we need to project the solution from the cells   , , whose vertices are   ± , respectively, and introduce the following quantities measuring the mesh shift (for their precise geometric meaning; see [19], Sect.3.2): )︁]︁ .
The conservative solution projection step from [22] is given by where are the point values reconstructed over the grid   , as described in Section 2.4.2.In particular, this projection step for the water surface variable  reads as Again, by the same reason as in the 1-D case, this projection step cannot be used in the "dry" areas, and therefore we replace the projection step (2.42) with where 2 , ≤ 0 or  +1, is "wet" and 2 , , otherwise, and (, ) d d,  into the four corresponding triangles and then using in each of the triangles the 7-point Gaussian quadrature, described in Appendix B. Finally, the "dry" cell correction term  corr , is defined as Remark 2.11.For long time simulations, the mesh grid distribution may not be balanced around the shocks, and as a result, more grid points can be concentrated on one side of the shock curve than on the other side.This can be fixed by increasing the number of iterations of moving mesh process, since within 1-4 iterations only few grid points can cross the shock curve.In this case, however, the resulting method will be extremely inefficient.
In order to improve the efficiency of the AMM algorithm, we re-project the solution back onto uniform mesh from time to time and then re-start the moving mesh adjustment.This significantly improves the distribution of resulting meshes.In all of our numerical experiments, we re-project the solution based on its second-order reconstruction.

Numerical examples
In this section, we test the developed AMM central-upwind schemes for the 1-D and 2-D Saint-Venant systems and compare the obtained results with the ones computed by the central-upwind scheme from [18] implemented on uniform meshes.
In all of the examples, the minmod parameter in the piecewise linear reconstructions is taken to be 1.3; see [19] for details.
The ODE systems (2.3) and (2.23) are numerically solved by the three-stage third-order SSP Runge-Kutta method.Each forward Euler stage of the SSP solver is described in Sections 2.1.3and 2.4.4 for the 1-D and 2-D case, respectively.
Remark 3.1.The parameter  will be taken different in different examples as the fraction of mesh points that needs to be concentrated in the "rough" areas of the solution depends on the structure of the solution, in particular, on the numbers and distribution of shocks.One can also select  dynamically in order to keep the smallest cell and/or time step sizes within a tolerable range.This is, however, beyond the scope of the current paper.

Example 1 -One-dimensional case
In this example, we solve the 1-D Saint-Venant system of shallow water equations on the computational domain [−, ] with the bottom topography consisting of one hump, the following initial conditions: (, 0) ≡ 1, (, 0) = −sign() cos .
and the (reflective) solid wall boundary conditions at both  = ±.We initially split the computational domain into  = 200 uniform finite-volume cells and compute the solution by the proposed AMM central-upwind scheme using the monitor function in (2.11) with  =   and the mesh concentration parameter  = 0.8.
In Figure 2, we plot the water surface  together with the bottom topography  at times  = 0, 0.5 and 2. As one can see, the AMM results are sharper than the ones obtained using the uniform mesh with  = 400 and are comparible with the ones computed using the uniform mesh with  = 1600.In Figure 3 (left), we present the corresponding time-space distribution of mesh cells, which clearly shows that the AMM is able to capture and follow the solution discontinuities.Finally, in Figure 3 (right), we compare the conservation errors in the computation of ℎ, which is measured by As we have explained in Section 2.3, this conservation error is caused by the use of a quadrature in the approximation of cell averages of the bottom topography in (2.16).We have tested the trapezoidal and Simpson's rules and, as one can see in Figure 3 (right), the use of Simpson's rule leads to much smaller conservation errors.

Example 2 -Accuracy test
In this example, we experimentally check the order of accuracy of the proposed AMM central-upwind scheme.
We take  = 0.5 and compute a sequence of numerical solutions using  = 50, 100, 200 and 400 (with the 400 × 400 solution being used as the reference solution) until the final time  = 0.07, by which the solution remains smooth.In Table 1, we show the  1 -,  2 -and  ∞ -norms of the errors in the computation of ℎ together with the experimental rates of convergence.As one can see, the rates are close to 2 as expected.

Example 3 -Waves in a water tank
In this example, we consider the initial-boundary value problem (IBVP) for the 2-D Saint-Venant system in the computational domain [−1, 1] × [−1, 1] with the (reflective) solid wall boundary conditions.The bottom topography function is and the initial conditions are given by We initialy split the computational domain into 200 × 200 uniform finite-volume cells and solve the studied IBVP using the designed AMM central-upwind scheme.We choose the monitor function with  = ∆ in (2.11) and  = 0.7.The contour plots of  and the corresponding meshes at times  = 0.125, 0.25, 0.375 and 0.5  The re-projection (see Rem. 2.11) onto the uniform mesh has been performed at times  = 0.025, 0.05, . . ., 0.475, 0.5, that is, every 0.025 s.This re-projection process leads to additional conservation errors, but it can significantly improve the quality of the mesh distribution.The time evolution of the conservation errors in ℎ and , as well as the comparison between meshes with/without this re-projection process at time  = 0.15 are shown in Figure 5.As one can see, without the re-projection procedure, a large portion of the grid points get trapped by the circular shocks, which leads to the lack of cells in the central part of the domain and thus prevents sharp resolution of shocks.

Example 4 -Asymmetric dam break
In this example, we test the designed scheme on a benchmark taken from [20].The computational domain is We initially split the computational domain into the uniform Cartesian cells of size ∆ = ∆ = 1 and solve the studied IBVP using the proposed AMM central-upwind scheme.During the moving mesh process, the vertices at the boundary are kept fixed and all of the other grid points initially located on the boundary are allowed to move according to the moving mesh equation, but only along the boundaries.We choose the monitor function with  = ∆ in (2.11) and  = 0.7.The contour plots of  and the corresponding meshes at times  = 1.5, 3 and 7.5 are presented in Figure 6.As one can see, the mesh is concentrated at the discontinuous and other rough parts of the solution, which is captured with an extremely high resolusion.
The mesh is initialized as follows.We split the rectangular domain [0, 3] × [−0.5, 0.5] into 300 × 100 uniform finite-volume cells whose vertices are denoted by  * )︁ , and project them onto the computational domain using 2 We then solve the IBVP using the designed AMM central-upwind scheme.During the moving mesh process, the grid points initially located on the upper and lower boundaries are evolved in time by changing their - coordinates according to the moving mesh equation, while their -coordinates are computed using the boundary function  =   () in order to keep these grid points on the corresponding parts of the boundary.We choose the monitor function with  = ∆ in (2.11) and  = 0.2.The contour plots of  and the corresponding meshes at times  = 0.5, 1.5 and 5 are presented in Figure 7.As one can see, the solution, including the steady state reached by time  = 5, has been accurately captured by the proposed AMM central-upwind scheme and the mesh is automatically adjusted to the computed solution structure.
We initially split the computational domain into 300 × 200 uniform Cartesian cells and solve the studied IBVP using the proposed AMM central-upwind scheme.We choose the monitor function with  = ∆ℎ in (2.11) and  = 0.3.In Figure 8, we show the time evolution of the dam-break wave propagating over the initially dry bed together with the corresponding meshes.It can be observed that as in previous examples, the designed scheme sharply captures a complicated wave structure and, in addition, the scheme performs well when handling both wetting and drying processes.
Appendix A. Proof of Theorem 2.5 We consider several different cases depending on where the cell   is located in.
Case 1 (  as well as  −1 and  +1 are "wet").In this case, equation (2.17) reduces to (2.15).Moreover, since we reconstruct  (not ℎ) in the cell   , we have )︁ , which after being substituted into (2.15)gives Due to the relaxation step (2.12), we have )︂ (), where the latter inequality is guaranteed by (2.13).Thus, the inequality (2.20) immediately follows from (A.2) and (A.3).Case 2 (  as well as  −1 and  +1 are "dry").In this case, we reconstruct ℎ (not ) in the cell   , and thus we have ℎ . Using this and the fact that    = ℎ   +    , we rewrite (2.17) as follows: , where we have used the following notation: Using (2.14), this equation can be rewritten as (compare with (A.1)): )︁  , which are nonnegative since the reconstruction of ℎ is positivity preserving.Therefore, ℎ +1 ,1 ≥ 0 and thus (A.2) is satisfied.Case 3 (  is "wet", while  +1 or  −1 or both are "dry").In this case, the cell   is located near the wetting/drying interface and the projection procedure depends on the type of the neighbouring cells and also on whether the cell interface moves into the "wet" or "dry" area.We will first consider the situation with only one "dry" neighbor, say,  +1 (Case 3a) and then we will consider the case when   is an isolated "wet" cell (Case 3b).Case 3a ( −1 and   are "wet", while  +1 is "dry").There are two possible situations: either the interface   > 0 )︁ into the "dry" area.In the former situation, the proof is identical to Case 1.We therefore consider the case  + , and then rewrite it as .We also note that since the reconstruction of ℎ is positivity preserving, ℎ + and that due to the relaxatoin step (2.12), we have ≥ 0. The desired bound (2.20) is then obtained by estimating the right-hand side (RHS) of (A.4) as follows: .
Case 3b (  is "wet", while both  −1 and  +1 are "dry").The only situation in which the projection step is different from the one that have been already investigated, is when both interfaces of the cell   propagate into the "dry" areas, that is, when In this case, we once again use (2.14) so that (2.17 ≥ 0, namely: Case 4 (  is "dry", while  +1 or  −1 or both are "wet").As in Case 3, the cell   is located at the wetting/drying interface and the projection procedure depends on the type of the neighbouring cells and also on whether the cell interface moves into the "wet" or "dry" area.We will first consider the situation with only one "wet" neighbor, say,  +1 (Case 4a) and then we will consider the case when   is an isolated "dry" cell (Case 4b).Case 4a ( −1 and   are "dry", while  +1 is "wet").There are two possible situations: either the interface   > 0 )︁ into the "wet" area.In the former situation, the proof is identical to Case 2. We therefore consider the case Taking into account (2.14), the last equation can be rewritten as: Since the piecewise linear reconstruction of  satisfies the right inequality in (A.2) and the cell  +1 is "wet", we have )︂ () Next, due to the relaxation step (2.12), we have ≥ 0, which together with the identity and positivity preserving property of the reconstruction of ℎ imply that 2 ≥ 0 and thus (2.20) is satisfied.Case 4b (  is "dry", while both  −1 and  +1 are "wet").The only situation in which the projection step is different from the one that have been already investigated, is when both interfaces of the cell   propagate into the "wet" areas, that is, when  )︁ + , where we have used the following notation: Similarly to Case 3a, we note that since piecewise linear reconstruction of  satisfies the right inequality in (A.2) and the cells  −1 and  +1 are "wet", we have . Therefore, we conclude that ℎ +1 ,3 ≥ 0 and the inequality (2.20) immediately follows.

Appendix B. Quadrature for general quadrilaterals
In this appendix, we provide a detailed description of the 2-D quadrature over quadrilaterals, which was used in our numerical experiments.Here, we consider a general quadrilateral , which may be convex, concave where the sign in front of each of the integrals on the RHS of (B.1) is "+" if the verticies of the triangle are listed in the counterclockwise order, and "−" otherwise.
We then approximate the integrals over each of the triangles using a 7-point Gaussian quadrature, for example, )︁ /6.Notice that in order to establish the inequality (C.1), one needs to use the inequality   ≥   .In the cells, where the latter inequality is not true, that is, where   <   , or in "dry" cells defined in (2.6), (2.7), we reconstruct ℎ instead of  and this approach gives )︁ .
We now prove the following theorem.
Theorem C.1.Let ℎ is evolved in time using the forward Euler discretization (2.10).If ℎ () ≥ 0 for all , and the time step satisfies the following CFL condition: Remark C.2.The proof of Theorem C.1 can be directly extended to the case when the time discretization is performed using the three-stage third-order SSP Runge-Kutta method as it can be written as a convex combination of forward Euler steps.

C.2. Two-dimensional CFL condition
We begin with proving the following lemma.
Lemma C.3.Assume that in the convex quadrilateral cell  , , the water depth ℎ is approximated using a nonnegative conservative linear function.Then,  , respectively, and the proof of the theorem is completed.

Figure 2 .
Figure 2. Example 1: Initial condition and solution profiles ( and ) computed by the AAM and uniform grid central-upwind schemes at different times.

Figure 3 .
Figure 3. Example 1: Time-space distribution of mesh cells (left) and time evolution of the conservation errors in the computation of ℎ (right).

Figure 4 .
Figure 4. Example 3: Time evolution of  (top row) and the corresponding meshes (bottom row).

Figure 5 .
Figure 5. Example 3: Time evolution of the conservation errors in  and ℎ (left); mesh distributions at  = 0.15 with (middle) and without (right) the re-projection process.

Figure 7 .
Figure 7. Example 5: Time evolution of  (left column) and the corresponding meshes (right column).

Figure 8 .
Figure 8. Example 6: Time evolution of  (plotted together with  in the left column), ℎ (middle column) and the corresponding meshes (right column).

2 > 0 .(A. 4 )
We now denote by  max  := max ∈(  −1 ,  +1 ) () and note that since  −1 and   are "wet" and we use the maximum principle preserving reconstruction of  in these cells,    ≥  max

2 > 0 .
In this case, we use(2.14)and the fact that    = ℎ

, 1 2 2 ,
, .Similarly, the third and fourth, the fifth and sixth, and the seventh and eights terms on the RHS of (C.10) and ∆ ≤ −  ,− 1 2 2 − ,− 1 2 and (2.5),  ± is the values of the monitor function  at the grid point  =    computed using the cell averages {︁ , are the values of the monitor function  at the grid points  =