ADAPTIVE PHYSICAL-CONSTRAINTS-PRESERVING UNSTAGGERED CENTRAL SCHEMES FOR SHALLOW WATER EQUATIONS ON QUADRILATERAL MESHES

. A well-balanced and positivity-preserving adaptive unstaggered central scheme for two-dimensional shallow water equations with nonflat bottom topography on irregular quadrangles is presented. The irregular quadrilateral mesh adds to the difficulty of designing unstaggered central schemes. In particular, the integral of the source term needs to subtly be dealt with. A new method of discretizing the source term for the well-balanced property is proposed, which is one of the main contributions of this work. The spacial second-order accuracy is obtained by constructing piecewise bilinear functions. Another novelty is that we introduce a strong-stability-preserving Unstaggered-Runge–Kutta method to improve the accuracy in time integration. Adaptive moving mesh strategies are introduced to couple with the current unstaggered central scheme. The well-balanced property is still valid. The positivity-preserving property can be proved when the cells shrink. We prove that the current adaptive unstaggered central scheme can obtain the stationary solution (“lake at rest” steady solutions) and guarantee the water depth to be nonnegative. Several classical problems of shallow water equations are shown to demonstrate the properties of the current numerical scheme.


Introduction
We consider two-dimensional shallow water equations which read as follows: where ℎ(, , ) represents the water depth, (, , ) and (, , ) denote the depth-averaged fluid velocity in -and -direction respectively,  is a constant related to the acceleration due to the gravity, (, ) denotes the bottom topography, which is time-independent.We first rewrite the system (1.1) as a compact form: where W := (ℎ, ℎ, ℎ) ⊤ is the conservative vector, and fluxes are and a source term is (1.4) Stationary solutions of the system (1.1) satisfies in which contains moving steady states (see [1][2][3][4]) and a "lake at rest" stationary solution which takes the following form  := ℎ +  = Const.,  =  = 0. ( It is indispensable to preserve the stationary solution at a discrete level to avoid the numerical storm [5].Numerical schemes that can preserve the stationary solution at the discrete level are called "well-balanced" scheme.We refer the reader to [1,[3][4][5][6][7][8][9][10][11][12] and references therein.Since the calculation of the acoustic speed √ ℎ is needed for determining the time step ∆ under the CFL condition, the numerical schemes need to guarantee the water depth ℎ to be nonnegative.The positivity-preserving solutions are also physically relevant.It is worth mentioning that the positivity-preserving property of the water depth is crucial in preserving the conservation of the water mass and is indispensable for successfully simulating physically relevant phenomena.For several positivity-preserving numerical schemes, we refer the reader to [6-9, 11, 12] and references therein. Central schemes have attracted tremendous attention in the past decade because these schemes are Riemannproblem-solver-free.In [13], they proposed "Lax-Friedrichs" (Lxf) schemes for hyperbolic conservation laws.Nevertheless, the Lxf scheme produces large numerical diffusions that smear discontinuous solutions (shocks, contact discontinuous etc.).In [14], they proposed a second-order central scheme termed as the "Nessyahu-Tadmor" (NT) scheme by replacing the piecewise constant function with a piecewise linear function in the MUSCL framework.They also proved that the NT scheme satisfies an entropy condition using homotopy methods.The NT scheme has subtle difficulty in preserving synchronization problems and the boundary conditions, which become more complex for the multi-dimensional system due to the application of staggered grids.Thus, unstaggered central schemes were introduced in [15].They inherited the advantages of the NT scheme and solved their difficulties.In [16], they proposed "Central-Upwind" scheme based on discretizing the system on exactly Riemann fans to reducing the numerical diffusions of the NT scheme.They also obtained semi-discrete central-upwind schemes, which can be discretized solvers of higher-order ODEs.The central-upwind scheme was extended for solving shallow water equations in [6].Unstaggered central schemes for shallow water equations can be found in [12,17].In [18], they extended the unstaggered central scheme for the Euler equation with a gravitation.In [9], they proposed unstaggered central schemes for the open channel flow, which are well-balanced and positivity-preserving.
Capturing the complex wave structures is a challenging work, especially for computing discontinuous weak solutions.In order to correctly reflect the complex wave patterns, finer computational meshes are needed.
Nevertheless, it makes the computational effort to be expensive.In the smooth region, we need a coarse mesh to compute the numerical solution, in the other region, we need a finer mesh to obtain a higher resolution.Adaptive moving mesh (AMM) methods can perfectly overcome these difficulties.The AMM methods can adaptively shift old mesh points that towards to the region where the large variation can be constructed.The new mesh points are constructed based on a mesh partial differential equation [19,20].In [21], they proposed "central-upwind" schemes for the hyperbolic system based on AMM methods using irregular quadrangles.In [22], they proposed a WENO reconstructions for compressible flows together with the AMM method on triangles.In [23], numerical schemes for two-dimensional ideal magnetohydrodynamics based on AMM methods were discussed.
In this paper, we propose well-balanced and positivity-preserving adaptive unstaggered central schemes on irregular quadrilateral meshes for 2-D shallow water equations with nonflat bottom topography.To the best of our knowledge, this paper is the first attempt designing adaptive unstaggered central schemes on irregular quadrangles.Compared with the adaptive approaches in [24,25], the current adaptive unstaggered central scheme is Riemann-problem-solver-free and uncoupled with the PDE solvers.Our method has less restrictions on the quality of the meshes.A key difficulty is to discrete the integral of the source term for the well-balanced property.Notice that the discretized source term should be consistent with the homogeneous shallow water equation for correctly reflecting the wave structure of shallow water equations without bottom topography.Motivating by the method in [11], we use the divergence theorem based on a new auxiliary vector (2.13) which results in well-balanced and consistent discretization of the source term.In particular, the adaptive unstaggered central schemes cannot preserve the "lake at rest" steady states when the computational domain contains wet-dry fronts.
It is nontrivial to prove the positivity-preserving property of the current unstaggered central schemes based on irregular quadrangles because of the complex geometrical structure of irregular quadrangles.The formulae (B.1) used to compute the cell average of the unknown variable is very complex.It is not easy to find the relation between the cell average of the water depth ℎ  + 1  2 ,+ 1   2 and the numerical flux (2.7) in proving the positivity of the water depth ℎ +1 + 1  2 ,+ 1   2 .It is remarkable that the positivity-preserving property of the forward step and the backward step can be addressed using the method [21].Even though the calculations are complex, we still found that the cell average of the water depth on the staggered cell has same terms based on the formulae (2.24), which makes it easy to prove the positivity-preserving property.We prove that the current scheme can guarantee the water depth to be nonnegative without using the "draining" time-step technique as in [12].
We would like to point out that, when moving the meshes, the well-balanced property is still valid and the positivity-preserving property can be proved when the cells shrink.We fix the points of the associated with the mesh when the adaptive current scheme cannot guarantee the water depth to be nonnegative.
It is remarkable that, in order to improve the accuracy of the time integration of the central scheme, the usual method is to apply the Taylor expansion or the Central Runge-Kutta (CRK) method proposed in [26] to compute the intermediate values.The first unstaggered central scheme [15] uses classical Runge-Kutta solvers to obtain the intermediate value for improving the accuracy of the time integration.The intermediate values were again computed using the Taylor expansions on the unstaggered cell.This methodology was applied to design well-balanced unstaggered central schemes for shallow water equations by the authors [18] and for the Euler equations by the authors [27].Despite the CRK method obtains high-order accuracy in time, the positivity of the physically relevant unknown variables (e.g.density, water depth, energy et al. ) is an open problem.In particular, the strong-stability-preserving (SSP) property is also an open problem.We follow the method proposed in [28] to discretize the ODEs (2.6) in time integration.This method was termed as SSP Unstaggered-Runge-Kutta method, which inherits the property of the SSP Runge-Kutta method proposed in [29,30].
We summarize the main contributions of the proposed adaptive unstaggered central scheme for twodimensional shallow water equations as: -The scheme is Riemann-problem-solver-free; -The scheme can adaptively capture complex wave structures; -The scheme has second-order accuracy both in time and in space; -The scheme can preserve the "lake at rest" stationary solution; -The scheme can guarantee the water depth to be nonnegative and is robust when the computational domain contains wet-dry fronts.
The paper is organised as follows.In Section 2, we introduce the unstaggered central scheme for 2-D shallow water equations with nonflat bottom topography.The discretized source term is investigated in Section 2.2.1.In Section 2.4, we introduce the strong-stability-preserving Unstaggered-Runge-Kutta method.The well-balanced and positivity-preserving properties are proved in Section 3. Adaptive moving mesh strategies are introduced in Section 4. We show several computed results obtained by the current scheme for 2-D shallow water equations in Section 5. A conclusion follows in Section 6.

Numerical schemes
In this section, we investigate a well-balanced and positivity-preserving unstaggered central scheme for twodimensional shallow water equations with nonflat bottom topography on quadrilateral meshes.The numerical scheme consists of three steps: the forward step, the corrector step, the backward step.Let us remark that the forward step is used to compute the cell average of unknowns on staggered cells; the corrector step is used to evolve the cell average of unknowns on staggered cells; the backward step is used to recover solutions on unstaggered cells.
It is worth mentioning that each step needs to be capable of preserving the steady states and guarantee the positivity of the water depth.Once the previous step cannot preserve the steady states, the errors will be accumulated.In particular, the velocity maybe vary large to break down the code.The well-balanced property of the forward step and the backward step can be obtained by using the surface gradient method (e.g.[5,9,12,31]) in the MUSCL framework.It is a challenging work to guarantee the well-balanced and positivity-preserving properties of the corrector step, since the discretized source term of the system (1.1) does not necessarily balance the difference of the flux.In particular, when the computational domain is discretized by using irregular meshes, guaranteeing the well-balanced property is more difficult.

The forward step
This section computes the cell average of the solutions on the staggered cell using the constructed piecewise bilinear polynomial, which is the approximate function of the exact solution on the unstaggered cell.We first introduce notations used in this paper and illustrate the method in constructing approximate numerical derivatives of the solution.We close this subsection in computing the cell average of solutions on staggered cells.

The discretization of the computational domain
The computational domain is divided into quadrangles.It is worth noticing that, the quadrangle is not necessarily to be uniform rectangles.We denote by  , = (︁ ), the unstaggered and staggered cell, respectively.We use x , := (  ,   ) and )︁ to denote the barycenter of the unstaggered cell  , and the staggered cell  + 1 2 ,+ 1 2 , respectively.Since the computational mesh is staggered, it is nontrivial to introduce the notations that we will use in showing the details of constructing numerical schemes.As showed in Figure 1, we first introduce the following notations associated with the information of the side:  ,+1 denotes the side connecting points x , and x +1, . ,+1 denotes the side connecting points x , and x ,+1 . ,−1 denotes the side connecting points x , and x −1, . ,−1 denotes the side connecting points x , and x ,−1 . 24 denotes the side connecting points  +  Next, we introduce the following notations associated with the area: Remark 2.1.We would like to point out that the discretized quadrangles have the uniform topological structure in this paper.Even so, the present unstaggered central scheme can be applied to the unstructured quadrangles.The numerical results showed in Section 5 are based on structured irregular quadrangles, which are obtained by disturbing uniform rectangles.

Second-order extensions
This section introduces the method for constructing piecewise bilinear polynomials using the cell average of the solutions on unstaggered cells.Since the mesh is not uniform, the construction based on dimension-by-dimension method is not valid any more.The piecewise constant function approximates the exact solution, the resulting numerical schemes obtain the first-order accuracy in space.Although the first-order scheme is easily constructed, it produces more numerical diffusions than higher order schemes.IN particular, the computed shocks may be seriously smeared.It is subtle in constructing higher order schemes, even for the second-order scheme.The constructed piecewise bilinear function should share the essentially non-oscillatory property.In [11], they used the information of neighboring cells to construct piecewise bilinear polynomials on triangles in constructing wellbalanced Central-Upwind schemes.Using the similar methodology, the authors [21] introduced robust Central-Upwind schemes on irregular quadrangles for compressible Euler equations and granular hydrodynamics.We use the method proposed in [21] to construct the piecewise bilinear polynomials on the irregular quadrangles.
For the convenience of the presentations, in what follows, we temporarily omit the superscript time notation  when there is no ambiguity.
Next, we begin with the following notations used in this section to give the details of constructions, )︁ ⊤ denotes the numerical derivative of the bilinear polynomial  1 , obtained by constructing a bilinear function which passes three points, x , , x ,−1 and x +1, .
denotes the numerical derivative of the bilinear polynomial  2 , obtained by constructing a bilinear function which passes three points, x , , x ,+1 and x +1, .
denotes the numerical derivative of the bilinear polynomial  3 , obtained by constructing a bilinear function which passes three points, x , , x −1, and x ,+1 .
denotes the numerical derivative of the bilinear polynomial  4 , obtained by constructing a bilinear function which passes three points, x , , x −1, and x ,−1 .
These notations are shown in Figure 2. The calculations of the above numerical derivative can be obtained in Appendix A.
In order to minimize the spurious oscillations and preserve the second-order accuracy, we compute the numerical derivative ∇W , in the MUSCL framework, i.e. the numerical derivative ∇W , is obtained by using the minmod function (e.g.[6,9,15,16]) based on the computed neighbouring four numerical derivatives.The same methodology was also considered in [21].Noting that, the constructed bilinear polynomial can preserve the conservation of the water mass.
The numerical derivative of the solution ∇W , on the unstaggered cell is computed by Remark 2.2.We would like to point out that, a direct application of the constructed bilinear function can not guarantee the positivity of the computed water depth due to the applied gradient limiter.The negative water depth is not physical and maybe break down the code.The positivity of the constructed water depth ̃︀ ℎ , (, ) can be guaranteed using the method proposed in [21] through introducing a positivity-preserving limiter, which will be discussed in latter.

Cell averages of solutions on staggered cells
In this section, we compute the cell average of the solution on the staggered cell  + 1 2 ,+ 1 2 using the constructed piecewise bilinear functions on the unstaggered cell  , .Comparing with the computational method showed in [12] whose computational domain is divided by uniform rectangles, one can not use midpoint formulae to compute the cell average of the solutions on distinct unstaggered cells.This case is more complicated due to the used irregular quadrangles.Even so, we can use the analogous methodology to compute the cell average of the solutions on the staggered cell.That is, we first divide the staggered cell into eight triangles, and then compute the cell average of the solutions on these triangles using the constructed bilinear functions on the unstaggered cell.
As showed in Figure 3, the cell average of the solution W on the staggered cell  + 1 2 ,+ 1 2 is defined as Since the constructed function is piecewise bilinear function, we divide the staggered cell into eight triangles.The calculations of the cell average of the solution on staggered cells can be obtained in Appendix B.
It is worth noticing that the conservative variables and the water surface are computed in the componentwise manner.Since we do not construct the bottom topography directly, then the cell average of the bottom topography is computed using the following method, via the equation (1.5).The cell average of the bottom topography depends on the temporal variable due to the cell average of the water surface and water depth depending on the temporal variable.
Figure 3. Schematic: the staggered cell of the computational domain.

The corrector step
This section evolves the cell average of solutions on staggered cells.We adopt the method of lines to discretize the system (1.2).This methodology allows us to convert the partial differential equations into ordinary differential equations (ODEs).Through only integrating the system in the spatial direction, we can obtain an ordinary differential equation with respect to the temporal direction.Then, one can use robust higher-order ODEs solvers to compute the numerical solutions.
In order to give a clear presentation, we first introduce following notations: 2 ,1 denotes the side connecting points x , and x +1, . + 1 2 ,+ 1 2 ,2 denotes the side connecting points x +1, and x +1,+1 . + 1 2 ,+ 1 2 ,3 denotes the side connecting points x +1,+1 and x ,+1 . + 1 2 ,+ 1 2 ,4 denotes the side connecting points x ,+1 and x , .We only integrate the system (1.2) on the staggered cell  + 1 2 ,+ 1 2 which results in the following form, d d Noting that, the flux is continuous with respect to the temporal variable, since the constructed bilinear function (2.4) on the unstaggered cell is continuous.This property is crucial in defining the numerical fluxes.Unlike the Godunov-type upwind scheme, in which the definition of the numerical fluxes is based on exact or approximate Riemann solvers emerged from the cell interface.The Godunov-type upwind scheme is timeconsuming in analyzing the characteristic fields of the Jacobian matrix of the flux.The central scheme does not need this information in defining numerical fluxes.This property is one of the key ingredients of central schemes, we refer the reader to [12,14,32] about the central schemes.
Using the divergence theorem, we can compute the integral of the flux, with the outer normal vector n of the corresponding side and  + 1 2 ,+ 1 2 denotes the boundary of the staggered cell  + 1 2 ,+ 1 2 .Using the trapezoidal quadrature formula, the right hand side of the equation (2.7) can be computed by 2 , is the angle of the outer normal vector of the -th side  + 1 2 ,+ 1 2 , of the staggered cell, and equipped with F , := F(W , ).We can compute the G + 1 2 ,+ 1 2 , using the similar method.Then, we can obtain a semi-discrete scheme for the system (1.2), here and S + 1 2 ,+ 1 2 is the approximate integral of the source term It is greatly nontrivial to design the discretization of the source term, because the discretization of the source term should be well balanced by the difference of the numerical flux when the system admits the steady states (1.5).
The author [12] designed well-balanced discretized source term based on the hydrostatic reconstruction methods on the uniform meshes.In [17], they designed well-balanced discretized source term using the constant water surface on uniform meshes.However, these methodology are both not extended to the irregular quadrangles.
We discuss the details of the discretization of the source term in the latter section.

The well-balanced discretization of the source term
A direct computation of the integral of the source term resulted in the numerical scheme that cannot preserve the stationary solution and fails to correctly capture its small perturbations.In [12], they proposed well-balanced discretized source terms based on hydrostatic reconstructions.Nevertheless, the methodology of the construction of the discretized source term discussed in [12] cannot be applied due to the irregular quadrangles.The discretized source term can be obtained by introducing an auxiliary vector.
As discussed in [33], the discretization of the source term can be obtained using the first Green's formula, here, the auxiliary vector is defined as As an example, we can compute the discretization of the source term associated with the derivative of the bottom topography along the  direction, ,+ 1 2 where,   and   are the numerical derivative of the water surface and the bottom topography, respectively.Notice that, the discretization of the source term (2.12) is not consistent with the homogeneous shallow water equation, i.e. the discretization of the source term is not vanish when the bottom topography to be a constant.Especially, even the bottom topography is zero, the discretization of the source term (2.12) doest not vanish.
In order to overcome this issue, we introduce a new auxiliary vector to discretize the integral of the source term.We define the vector We compute the integral, According to the first Green's formula, we obtain here,   is the -component of the outer normal vector n.Using the above formula, we can compute the discretization of the second component of the source term, which is computed as ,+ 1 2 with introduced notations Similarly, we can compute the discretization of the third component of the source term, which is computed as with introduced notations where, we compute )︀ , )︀ , It is worth mentioning that the constructed bilinear function (2.4) is discontinuous at the cell interface.Thus, the point values W + 1 2 ,+ 1 2 , can not be directly computed.The equation (2.18) is used to define the arithmetic mean value ̃︁ W(, ) at the cell interface.We notice that the water surface is a constant when the system (1.1) admits steady states (1.5) and does not contain the wet-dry fronts.In this case, the approximate water surface and velocity functions are continuous.
We observe that, the discretization of the source term is only consistent with the homogeneous shallow water equation when the bottom topography to be zero, but this property is not shared by the discretized methodology of [11] where the discretized source term does not vanish even the bottom topography become zero.The discretized source term of the current unstaggered central scheme is carried out based on the new auxiliary variable (2.13).We will prove the well-balanced property of the current unstaggered central scheme in Section 3.
Finally, we summarize the discretized source terms as It is worth mentioning that the discretized source terms (2.16) and (2.17) may be discontinuous.
Remark 2.4.The discretized source term (2.16) and (2.17) is consistent with the homogeneous shallow water equation since the  −  terms have been replaced by  − ℎ terms.Then, the water surface equals to the water depth when the bottom topography vanishes.
Remark 2.5.We would like to stress that the discretized source term (2.16) and (2.17) is modified to be zero when the cell average of the water depth on the staggered cell vanishes, i.e. ℎ + 1 2 ,+ 1 2 = 0 for ,  ∈ Z.

The backward step
This section recovers solutions on the unstaggered cell using the evolved solutions on staggered cells.It is worth mentioning that this step is crucial in avoiding the alternate mesh at each time step.The key difference between the unstaggered central scheme and staggered central scheme is that the unstaggered central scheme only integrates the system (1.1) on the mesh  + 1  2 at each time steps, but the staggered central scheme integrates the system (1.1) on the mesh  , and  + 1  2 at even time steps and odd time steps, respectively.The staggered central scheme makes the boundary conditions that are more complex than the unstaggered central schemes.We refer the reader to [15] for more details.
For the sake of brevity, we first omit the time mark  + 1.The piecewise bilinear function is denoted by ︁ W(, ) on the staggered cell, which can be constructed using the analogous method in the forward step, then here ( 2 can be obtained using the similar method discussed in Section 2.1.2.
The cell average of the solution on the unstaggered cell can be obtained in Appendix C. It is worth noticing that, we construct the water surface, the discharge along the -direction and the discharge along the -direction for guaranteeing the well-balanced property that can be maintained in the backward step.This methodology is related to the surface gradient method discussed in [12,34].Thus, the cell average of the water depth is computed as on unstaggered cells.
Remark 2.6.The water depth obtained by the equation (2.22) cannot be guaranteed positivity since the constructed water surface maybe lower than the bottom topography, especially near dry areas.The negative water depth is not physically relevant and make the simulation within the risk of break.For the positivitypreserving property of this step, we replace the formula (2.22) as follows that is, we construct the water depth instead of the water surface when the water depth (2.22) becomes negative.
Remark 2.7.We would like to point out that the positivity of the constructed water depth ̂︀ ℎ + 1 2 ,+ 1 2 (, ) can be guaranteed using the similar method discussed in Section 3.2.
Remark 2.8.The velocities are solved through using the following formula, , otherwise, ( here  = 10 −9 is a parameter to switch the formula.This method is used to avoid the cancellations near dry areas.We refer the reader to [6,9,11] for more details.

The discretization of the system in the time integration
This sections solves the time integration of the ODEs (2.9) for improving the accuracy in time.It is worth mentioning that, the first unstaggered central scheme [15] uses classical Runge-Kutta solvers to obtain the intermediate value for improving the accuracy in time.The intermediate values are computed using the Taylor expansions on the unstaggered cell.This methodology was applied to design well-balanced unstaggered central schemes for shallow water equations by the authors [18] and for the Euler equations by the authors [27].Based on this methodology, the authors [26] proposed Central Runge-Kutta (CRK) method for obtaining a high-order accuracy in time.
Although the CRK method obtains high-order accuracy in time, the positivity of the physically relevant unknown variables, such as, density, water depth, energy, etc., is an open problem.Especially, the strongstability-preserving (SSP) property is also an open problem.
We follow the method proposed in [28] to discretize the ODEs (2.9) in time integration.This method was termed as Unstaggered-Runge-Kutta (URK) method, which inherits the property of the SSP Runge-Kutta method proposed in [29,30].In this paper, we use the second-order URK method, which is defined as ; , 2 2 ; , where the residuum is defined as (2.26) The intermediate values ̃︁ W , are obtained at the backward step based on the values W (1) . We can compute the intermediate cell average ̃︁ W (1) 2 at the forward step based on the values ̃︁ W , .We can obtain the values W +1 , using the backward step based on the values W +1 + 1  2 ,+ 1 2 .Remark 2.9.The URK scheme can obtain the strong-stability-preserving property together with the TVD reconstruction of the forward step and the backward step.The updated value W +1 + 1  2 ,+ 1

Well-balanced and positivity-preserving properties
This section proves that the unstaggered central scheme can preserve the stationary solution and guarantee the water depth to be nonnegative.

Well-balanced properties
When the system (1.1) admits the stationary solution, the gradient of the physical flux can be balanced by the source term in the smooth sense.We do not consider the case where the computational domain contains wet-dry fronts when the water is at rest.Therefore, the numerical scheme should satisfy that the difference of the numerical flux should be balanced by the discretized source term, that is, the designed numerical scheme is well-balanced.We begin with the proof of well-balanced properties of the current scheme.
Proof.Since the first equation of the system does not contain the source term, we only consider the second and third component of the numerical scheme.Thanks to the central numerical flux, we can obtain d d ℎ + 1 2 ,+ 1 2 () = 0 due to the velocity vanishes.Indeed, we only need to prove that d d We can use the analogous method to obtain d d () = 0. Thus, for the stationary solution (1.5), we can obtain here,  is a constant.Then, using the equation (2.7), we compute the numerical flux Using the equation (2.16), we can compute the discretized source term In the last line, we have used   = 0 according to the condition (3.2).Notice that, we have Next, we discuss the positivity-preserving property of the unstaggered central scheme.The positivitypreserving property of the water depth is crucial in preserving the conservation of the water mass and is indispensable for successfully simulating physically relevant phenomena.In this paper, it is nontrivial to prove the positivity-preserving property, since the present unstaggered central scheme is designed on the irregular quadrangles.The formulae (B.1) used to compute the cell average of the unknown variable is very complex.It is difficult to couple with the numerical flux (2.7) in proving the positivity of the water depth ℎ 2 .It is remarkable that, the positivity-preserving property of the forward step and the backward step can be addressed using the method [21].Even though the calculations are complex, we still found that, the cell average of the water depth on the staggered cell has same terms based on the formulae (2.24), which makes it easy to prove the positivity-preserving property.

Positivity-preserving modifications
This section redefines the bilinear polynomial associated with the water depth to guarantee that the water depth is nonnegative at eight points.As showed in Figure 1, we redefine the polynomial associated with the water depth when one of these point values is negative.
We first give the constructed polynomial of the water depth Notice that, we can compute the cell average of the water depth as According to the conservation of the water mass, we obtain the following lemma.and Proof.According to the definition of the cell average of the water depth ℎ , , we can directly compute here After several calculations give here It is not difficulty to verify that  = ∑︁    , for  = 1, . . ., 8. (3.8) Then, the cell average of the water depth ℎ , is a convex combination of the point values ̃︀ ℎ , (  ).We complete the proof.
Thanks to the gradient limiter, the computed point values of the water depth ℎ  does not necessary to be nonnegative.We use the method proposed in [35] to redefine a bilinear polynomial to approximate the water depth function while inherits the second-order accuracy in space.The reconstructed bilinear polynomial is given by here Obviously, the conservation of the water mass is still valid based on the reconstructed polynomial (3.9).We can obtain the following theorem.
Theorem 3.3.Assuming the cell average of the water depth ℎ , ≥ 0, then the point values of the water depth satisfy providing the reconstructed polynomial (3.9).
In order to prove the positivity-preserving property of the present unstaggered central scheme clearly, we first use the first-order forward Euler method to discretize the scheme (2.9) in time integration.The notations   and  +1 :=   + ∆ are used to denote the corresponding time level.Theorem 3.4.Considering the two-dimensional shallow water equation with nonflat bottom topography (1.1).The unstaggered central scheme (2.9) endows with the numerical flux (2.7) and the discretized source term (2.16) and (2.17).Assuming the water depth is nonnegative at points   ,  = 1, . . ., 8 and ( + 1 2 ,+ 1   2 ), and providing the following CFL condition here, and Proof.In the following process, for the sake of clarity, we first omit the time notation .Using the equation (B.1), we can compute the cell average of the water depth on the staggered cell The equation (3.17) is complex, we can not see the relations with the numerical fluxes (2.7) directly.Notice that, the water depths computed at nine points,  1 , . . .,  8 and the point )︁ are nonnegative, it is a key ingredient of the current scheme in proving the positivity-preserving property.
After several simple calculations give We observe that, the numerical fluxes includes the information of four points, which are points  2 ,  4 ,  6 ,  8 .As showed in Figure 3.This observation helps us rewrite the expression of the water depth and the numerical fluxes.According to the equation (2.7), we can compute the first component of the numerical flux as )︁(︁(︁ ̃︀ ℎ +1,+1 ( 4 )̃︀  +1,+1 )︁ ( 4 ) + ̃︀ ℎ ,+1 ( 6 )̃︀  ,+1 ( 6 ) A direct calculation gives the equivalent form, Let us define the following notations for the ease of presentation, According to the first equation of (2.9), and omitting the time mark , we obtain )︂ )︂ )︂ Since the water depth is nonnegative at all points with the aid of the CFL condition (3.14)-(3.16),all the coefficients are nonnegative, then we obtain ℎ +1 + 1  2 ,+ 1 2 ≥ 0. We completed the proof.Remark 3.5.It is worth recalling that, the positivity-preserving property is still valid when using the Unstaggered-Runge-Kutta method to discretize the ODEs in time integration, since its solutions can be seen as a convex combination of several first-order forward Euler solvers.

Adaptive moving mesh strategies
This section proposes an adaptive moving method on quadrangles.We consider the structured irregular quadrangles.We first define the logical domain denoted by Ω  and the physical domain denoted by Ω.As showed in Figure 4, we use   to denote the -th side of the dual cell  + 1 2 ,+ 1 2 containing the vertex ( + 1 2 ,  + 1 2 ).The mapping between the logical domain and the physical domain is denoted as x : Ω  → Ω.
As discussed in [19,23], the mapping x(, ) is the solution of the following mesh PDEs, which is defined as The equation (4.1) is known as the Euler−Lagrange equation for minimizing the energy functional equation discussed in [19], here ∇ := (  ,   ) and  := (W) is a monitor function which will be given later.We assume that the logical domain Ω  is fixed.Notice that, the vertex x of the new mesh is obtained by solving the mesh PDEs (4.1) Integrating the equation (4.1) and using the divergence theorem, we obtain here, n  is the outward normal vector of the side   .We discrete the equation (4.2) as here We first denote by  the iterative times.Since the equation (4.3) is highly nonlinear, we use a iterative method to move the vertex here, ̂︀ and We choose the parameter  to be 0.5 and the iterative times  to be 5 in our all numerical examples when there is no ambiguity.It is mentioning that, the vertex of the new mesh x +1 + 1  2 ,+ 1 2 is contained in the convex domain which consists of the four points as showed in Figure 5, since the equations (4.4) and (4.5) are a convex combination of the old mesh vertexes.Remark 4.1.We use the method proposed in [19,23] and the reference therein to obtain a smoother mesh.The monitor function is given by Remark 4.2.The monitor function in (4.1) can be obtained in [19,21,22,36].In all numerical examples, we use the monitor function proposed in [22,36] which reads as here and |∇ (,) W| denotes the model of the gradient of the numerical solutions, and )︃ , and ∆ and ∆ are the mesh size of the logical domain.The parameters  and  are taken as 50 and 0.25, respectively.

Geometrical conservative interpolation
In this section, we introduce the method discussed in [19,23] for projecting the solutions from old meshes to new meshes.When obtaining the new meshes, we need to compute the cell average of the solutions on new meshes using the information of old meshes.For the shallow water equation with a source term, the designed technique recovering the solutions on new meshes need to satisfy well-balanced and positivity-preserving properties.Especially, it should satisfy the conservation of the water mass.
We denote  the old mesh and  + 1 the new mesh.Using the conservative method in [23], the solutions on the new mesh can be obtained by here, the following values are computed on the old mesh when there is no ambiguity, and We denote by | ̂︀  , | the area of the new mesh ̂︀  , , and   denotes the scanned area of the -th side.For an example, as showed in Figure 6, we compute )︁]︁ .
(4.10)We use the reconstructed polynomial (3.9) to recover the water depth on the new meshes for positivity-preserving property.The method (4.9) can preserve the still-water steady state.Indeed, when the water is at rest and the water surface is assumed to be a constant , we have

It is easily to verify
Then, we can obtain  +1 =   .Next, we discuss the positivity-preserving property of the adaptive moving mesh strategies.We observe that the current strategy admits the following proposition.Proposition 4.3.We consider two-dimensional shallow water equations (1.1) and the corresponding numerical scheme discussed in Section 2. Assume that ℎ  , , ℎ ,  , ℎ ,  are nonnegative for all ,  ∈  .Then, we obtain  Proof.Using the first component of the equation (4.9), we obtain As showed in Figure 7, we can rewrite the cell average of the water depth here ̂︀ ℎ  is taken from (3.11).After elementary calculations give It is mentioning that, in this case, we have Substituting the equations (4.13) and (4.14) into (4.12),we obtain We denote by |  | the size of the scanned area.As showed in Figure 7, we obtain Remark 4.4.The positivity-preserving property of the computed water depth using the adaptive moving mesh methods is still an open problem.It is not easy to prove this property in the theoretical sense.When the adaptive current scheme cannot guarantee the water depth to be nonnegative, the points of the associated with the mesh are fixed.The computed results confirm that, the current adaptive moving mesh strategies are positivity-preserving.
Remark 4.5.As discussed in [23], the method (4.9) is conservative in the sense, Remark 4.6.We would like to point out that the new quadrangle will be not convex.When the new quadrangle is not convex, we fix the related points to overcome this issue.

Numerical examples
In this section, we show several numerical experiments to test the present unstaggered central scheme for shallow water equations with a bed slope source term on irregular quadrangles.We use the CFL number to be 0.3 in all numerical examples when there is no ambiguity."OLD" denotes the numerical solutions obtained by using the discretized source term (2.12); "NEW" denotes the numerical solutions obtained by using the discretized source term (2.19).The gravity is taken as  = 9.8 when there is no ambiguity.
We use transmissive boundary conditions in all directions.The gravitational constant is taken to be  = 9.8.The computational domain is taken as [0, 2] × [0, 1] and the final time is taken to be  = 0.07.A reference solution is computed by the present scheme using 800 × 800 cells.We compute the  1 errors of the water depth and discharges.We use the following formula to test the order here,   is the  1 error obtained on the meshes with size ∆  .As showed in Table 1.We can expect a second-order accuracy in space on regular quadrilateral meshes.This case also confirms that the SSP URK method (2.25) is robust and can improve the accuracy in time.

Dam break over a flat bottom
In this example, we will show the numerical results obtained by the current unstaggered central scheme using the discretized source terms (2.19) and (2.12).The bottom topography is set to be zero.The initial water depth and velocities are given by ℎ(, , 0) = {︃ 0.5, if  < 0.5 0.3, otherwise , (, , 0) = (, , 0) = 0, respectively.We consider a rectangular computational domain, which has the size [0, 1]×[0, 1].There has 200×10 irregular quadrilateral meshes to divide the computational domain.We run this simulation until the final time to be  = 0.1.It is known that the structure of solutions of this initial-value problem is a left-going rarefaction and a right-going shock.The computed water depth and the velocity are shown in Figure 8.One can observe that the "NEW" scheme produced reasonable results, which is more agreement with the exact solutions than the "OLD" scheme.The reason produced this distinction is that, the "New" scheme has the consistent discretized source term with the exact one.However, the "OLD" scheme does not necessarily share this property.

Dam break over a nonflat bottom
In this example, we consider a dam break over a nonflat bottom topography.The bottom topography is given by )︂)︂ .
The computed water surface and the velocity  are shown in Figure 9.One can observe that, the current scheme can correctly reflect these wave patterns.This example also confirmed the well-balanced property of the current scheme.
Finally, we change the initial data to be ℎ(, , 0) = {︃ 2 − (, ), if  < 0.5 0, otherwise , (, , 0) = (, , 0) = 0, for testing the robust property of the current scheme when simulating the water flows over a dry area.The initial mesh is set to be 100 × 10 irregular quadrangles.The final time is taken to be  = 0.05.We showed the numerical results in Figure 10.One can see that the adaptive current scheme can adaptively capture the wave structures, in which the solution has large gradient.This case verified the robust property of the current adaptive scheme.

Circular dam break
In this example, we consider a circular dam break problem.The bottom topography is zero.The computational domain is taken to be [0, 1] × [0, 1].We consider the initial data as ℎ(, , 0) = {︃ 0.5, if  < 0.5 0.1, otherwise , (, , 0) = (, , 0) = 0, here,  = √︀ ( − 0.5) 2 + ( − 0.5) 2 .We compute numerical results using the current unstaggered central scheme with the discretized source term (2.19) and (2.12) on 100 × 100 irregular quadrilateral meshes until the final time  = 0.1.One can see that, a symmetrical wave emerged using both discretized source term (2.19) and (2.12), as showed in the Figure 11.However, the unstaggered central scheme equipped with the discretized source term (2.19) produced higher resolutions.Further, we showed the water depth along the line  = 0.5 in Figure 12.The reference solution was obtained by refining meshes to be 800 × 800 irregular quadrangles.Since the discretized source term (2.12) is not consistent with the homogenous shallow water equation, which produced lower resolutions.
We consider the computational domain to be [0, 2] × [0, 1].The computational mesh is shown in Figure 13.We first test the ability of preserving stationary solutions.Then, the small perturbation is first taken to be  = 0.The computed results are obtained by the current scheme using 200 × 100 irregular quadrangles until the final time  = 0.1.One can see that, the computed discharges are almost zero but within the machine accuracy in Figure 14.This confirms that the current scheme is capable of preserving the stationary solutions.Next, we investigate the ability of the current scheme in capturing small perturbations by taking the parameter to be  = 10 −3 .We run this simulation using the current adaptive unstaggered central scheme equipped with the discretized source term (2.19) on 100 × 50 and 200 × 100 irregular quadrangles until the output time  = 0.6, 1.2, 1.5, 1.8.The numerical solutions are shown in Figures 15-18.Notice that, no spurious oscillations can be observed and the small perturbations can be captured sharply.The adaptive unstaggered central scheme produced more resolutions than the unstaggered central scheme on 100×50 irregular quadrangles.The computed results also confirm that the adaptive unstaggered central scheme is well-balanced and positivity-preserving.
Finally, we test the well-balanced property of the current scheme when the computational domain contains wet-dry fronts.The initial data is taken to be ℎ(, , 0) = max(0.7 − (, ), 0.0), (, , 0) = (, , 0) = 0.The numerical solutions are obtained by the scheme on 100 × 50 irregular quadrangles until the final time  = 0.01.One can see that in Figure 19, the computed velocities are large greater than the machine accuracy.This confirms that the current scheme cannot preserve the "lake at rest" steady states when the computational domain contains wet-dry fronts.

Water drop phenomena
We consider water drop phenomena taken from [8].The computational domain is taken as [0, 1] × [0, 1].We first take the bottom topography to be (, ) = 0.The initial data reads as (, , 0) = 1 + 1 10 e −100[(−0.5) 2 +(−0.5) 2 ] , (, , 0) = 0, (, , 0) = 0. (5.1) We adopt the current adaptive unstaggered central scheme using 100 × 100 irregular quadrangles to run this simulation.The final time is taken to be  = 0.15, 0.3, 0.5, 0.7.We apply solid walls in all directions.The computed results are shown in Figure 20.One can see that, the obtained water surface agrees with that of in [8].We showed the final mesh in Figure 21.The adaptive moving mesh strategies can correctly capture the larger gradient regimes of the solutions.In order to further verify the robust property of the adaptive unstaggered central scheme, we run this simulation on 200 × 200 irregular quadrangles.The computed results are shown in Figure 22.One can observe that, the complex wave structures can be correctly reflect.Next, we take the bottom topography to be The other conditions are not changed.The computed results are shown in Figure 23.We can see more complex waves due to the effect of the bottom topography.It is remarkable that the computed water surface by the current unstaggered central scheme agrees with that of [8].The final meshes are shown in Figure 24.We also recompute the solution using the adaptive unstaggered central scheme on 200 × 200 irregular quadrangles.The  numerical results are shown in Figure 25.The complex wave patterns are captured by the scheme.Finally, we test the conservation errors of the water depth and the bottom topography when using the adaptive unstaggered central schemes.Numerical results are shown in Figure 26.One can see that the adaptive unstaggered central schemes can preserve the conservation of the water depth for the water drop without the bottom topography and conservation errors can be seen for the water drop with the bottom topography.These conservation errors of the water depth is caused by the conservation errors of the bottom topography.

Flows in converging-diverging channels
This test case was taken from [8,11], which models the shallow water flow over an open converging-diverging channel that has a symmetric construction with length 1 at the center.The length of the channel is taken to be 3 and the width is defined by the function   (), which reads as and  is the minimum channel breadth.We first take the parameter  to be 0.9.The sketch of the channel is shown in Figure 27.
We first consider a flat bottom topography that is taken to be (, ) = 0.An inflow boundary with  = 2 is taken at the left end side of the channel and an outflow boundary locates at the right end side of the channel.The solid walls are taken in -directions.Notice that, since the discretized source term (2.19) is consistent with the homogenous shallow water equation, then it is to be zero.However, the discretized source term (2.12) does not satisfy this property.
As discussed in [8,11], the computational quadrilateral meshes are obtained from the structured one using the mapping otherwise.
The numerical solutions are obtained by using two different meshes: 100 × 100 grids and 200 × 200 grids.We run this simulation until  = 2.We also run this simulation using the adaptive unstaggered central scheme on  100 × 100 grids for verifying the robust property.We showed computed results in Figure 28.One can see that computed results match well with that of [8,11].The current adaptive unstaggered central scheme produced more resolution than that without using the adaptive strategies.
Next, we consider a nonflat bottom topography to be (, ) = e −10(−1.9) 2 −50(−0.2) 2 + e −20(−2.2) 2 −50(+0.2) 2 . (5. 2) The initial conditions and boundary conditions are the same as the flat bottom topography.We showed the numerical results obtained by the current unstaggered central scheme on 100 × 100 and 200 × 200 grids until the final time  = 0.7 in Figure 29.We again run this simulation using the adaptive unstaggered central scheme on 100 × 100 grids for verifying the robust property.We can see more complex wave structures due to the effect of    the nonflat bottom topography.Even though the adaptive unstaggered central scheme produced more resolution than that of the unstaggered central scheme, the computed results are in agreement with that of [8,11].

Mach reflection problems
This final example is taken from [38].We take the bottom topography to be (, ) = 0.The initial data is given as (5.3) In this case, we take the gravitational constant to be  = 9.8.We consider solid walls in -directions and free boundaries are considered in -directions.The numerical results are obtained by the current unstaggered central scheme on 100 × 100 and 200 × 200 quadrangles.We consider the computational time to be  = 1.5.One can observe the computed results in Figure 30, in which a Mach stem emerged at the low boundary.Notice that, a lower resolution can be seen when using 100 × 100 quadrangles, it can be improved when using 200 × 200 quadrangles.In order to verify the robust and efficient properties of the adaptive unstaggered central scheme, we also run this simulation on 100 × 100 quadrangles.As showed in Figure 30, one can see that the adaptive unstaggered central scheme improved the resolution than that without using the adaptive strategies on 100×100 quadrangles.

Conclusions
We proposed a well-balanced and positivity-preserving adaptive unstaggered central scheme for the twodimensional shallow water equation with nonflat bottom topography on irregular quadrilateral meshes.We constructed piecewise bilinear functions for obtaining the second-order accuracy in the space and constructed the SSP Unstaggered-Runge-Kutta method to obtain the second-order accuracy in the time integration.We proved that the current adaptive unstaggered central scheme could preserve the stationary solution within the machine accuracy.The computed water depth was proved nonnegative when applying the current scheme.We

Appendix C. Calculations of the cell average of the solution on unstaggered cells
We compute the cell average of the solution on the unstaggered cell as

Figure 2 .
Figure 2. Schematic: stencils of constructing slopes on the unstaggered cell.

Figure 4 .
Figure 4. Moving mesh.Schematic of the mapping between the logical domain and physical domain.

Figure 5 .
Figure 5.The convex domain consists of the four points.

Figure 6 .
Figure 6.A sketch of a moving typical quadrangle.

Figure 7 .
Figure 7.A sketch of a typical quadrangle.

∑︁
, since two neighboring quadrangles share the opposite signs flux but with same absolute values.

Figure 8 .
Figure 8. Dam break over a flat bottom.Left column: computed water depth; Right column: computed velocity.

Figure 9 .
Figure 9. Dam break over a nonflat bottom.Left column: computed water surface; Right column: computed velocity .

Figure 10 .
Figure 10.Dam break over a nonflat bottom.Left column: computed water surface; Right column: final meshes.

Figure 20 .
Figure 20.Water drop phenomena.The bottom topography is taken as zero.The computed water surfaces are shown at different final times.

Figure 21 .
Figure 21.Water drop phenomena.The bottom topography is taken as zero.The total cells are 100 × 100.The final mesh is shown at different final times.

Figure 22 .
Figure 22.Water drop phenomena.The final time is  = 0.3 and the total cells are 200 × 200.The bottom topography is taken as zero.Left: computed water surface; Right: final mesh.

Figure 23 .
Figure 23.Water drop phenomena.The bottom topography is not flat.The computed water surfaces are shown at different final times.

Figure 24 .
Figure 24.Water drop phenomena.The bottom topography is not flat.The total cells are 100 × 100.The final mesh is shown at different final times.

Figure 25 .
Figure 25.Water drop phenomena.The final time is  = 0.3 and the total cells are 200 × 200.The bottom topography is taken as zero.Left: computed water surface; Right: final mesh.

Figure 26 .
Figure 26.Water drop phenomena.The final time is  = 0.3.Left: conservation errors of the water depth of the water drop without the bottom topography; Right: conservation errors of the water depth and the bottom topography of the water drop.

Figure 27 .
Figure 27.Flows in Converging-Diverging Channels.The final mesh is shown without adaptive strategies.

Table 1 .
Order checking:  1 -norm and numerical orders of accuracy.