GU Juan, PENG Xindong, DAI Yongjiu, and CHE Yuzhang

Fourth-Order Conservative Transport on Overset Grids Using Multi-Moment Constrained Finite Volume Scheme for Oceanic Modeling

GU Juan1), 3), PENG Xindong2), *, DAI Yongjiu4), and CHE Yuzhang5), 6)

1) College of Global Change and Earth System Science, Beijing Normal University, Beijing 100875, China 2) State Key Laboratory of Severe Weather, Chinese Academy of Meteorological Sciences, Beijing 100081, China 3) College of Air Traffic Management, Civil Aviation Flight University of China, Guanghan 618307, China 4) School of Atmospheric Sciences, Sun Yat-Sen University, Guangzhou 510275, China 5) Plateau Atmosphere and Environment Key Laboratory of Sichuan Province, Chengdu University of Information Technology, Chengdu 610103, China 6) Department of Mechanical Engineering, Tokyo Institute of Technology, Tokyo 152-8550, Japan

With an increase in model resolution, compact high-order numerical advection scheme can improve its effectiveness and competitiveness in oceanic modeling due to its high accuracy and scalability on massive-processor computers. To provide high-quality numerical ocean simulation on overset grids, we tried a novel formulation of the fourth-order multi-moment constrained finite volume scheme to simulate continuous and discontinuous problems in the Cartesian coordinate. Utilizing some degrees of freedom over each cell and derivatives at the cell center, we obtained a two-dimensional (2D) cubic polynomial from which point values on the extended overlap can achieve fourth-order accuracy. However, this interpolation causes a lack of conservation because the flux between the regions are no longer equal; thus, a flux correction is implemented to ensure conservation. A couple of numerical experiments are presented to evaluate the numerical scheme, which confirms its approximately fourth-order accuracy in conservative transportation on overset grid. The test cases reveal that the scheme is effective to suppress numerical oscillation in discontinuous problems, which may be powerful for salinity advection computing with a sharp gradient.

multi-moment constrained finite volume scheme; oceanic modeling; overset grid; conservation; numerical transport

1 Introduction

Advection is the most essential process in oceanic dynamics. Ocean models, such as the modular, Princeton, Hybrid coordinate, and regional ocean models, are widely used for high-resolution oceanic simulation. The models address the difficulty of sharp gradient description and detailed structure simulation of mesoscale eddies with popularly used second-order schemes. Highly accurate transport schemes are attractive in numerical model development. Over the last three decades, a common practice in computational fluid dynamics is to construct high-order schemes by using locally increased degrees of freedom (DOFs), such as finite element methods (including continuous and discontinuous Galerkin methods (Bastian, 2014), spectral element method (Taneja and Higdon, 2018)), finite difference method (Sun, 2007), and finite volume method (Wang, 2004; Ii and Xiao, 2007, 2009; Chen and Xiao, 2008; Xie, 2014). The development of higher-order schemes requires additional spatial information to establish a proper high-order polynomial. Increasing the order of spatial computation accuracy involves two ways. One is to employ a wider stencil (Liu, 1994; Jiang and Shu, 1996), which requires more halo grids in the parallel computation. A scheme beyond the second order usually results in a significant increase in algorithmic complexity and computational cost. Another way is to increase the amount of local information (Tanaka, 2000; Xiao and Yabe, 2001; Yabe, 2001; Xiao, 2002; Peng, 2004; Imai, 2008; Onodera, 2011), for example, the multi-moment constrained finite volume (MCV) method used in this study (Ii and Xiao, 2009), which fully utilizes the DOFs of a cell. This method achieves higher-order reconstruction on a more compact stencil that describes the discontinuity better with less spatial numerical diffusion, and therefore benefits the detailed structure captured by the model. The MCV method, which is developed under finite volume definition, has been applied to computational fluids based on structured and unstructured grids (Li, 2008, 2015b; Ii and Xiao, 2010; Chen, 2014; Xie and Xiao, 2014; Deng, 2017a, 2017b; Jin, 2018). The fourth-order three-point MCV formulationwithboundarygradientswitching(MCV3-BGS)(Deng, 2017a) set the local DOFs in each cell as the constraint conditions to ensure oscillation-free computation (see Ii and Xiao (2009) for further details). The multi-moment quantities include point value (PV), volume-integrated average (VIA), and derivatives of physical quantity. The BGS method is an essentially non-oscillatory (ENO)- like limiter that redefines boundary fluxes. With the help of BGS, numerical oscillations are successfully suppressed, and the MCV3-BGS scheme has become a promising tool for oscillation-free high-order accuracy computing, which is attractive for high-resolution simulation of ocean salinity and sharp temperature gradient; high-resolution global climate prediction is necessary to improve China’s extended range prediction (Zhang, 2018).

Overset grids are suitable for fluid simulations on a complex surface for easy grid generation and high-performance computation, such as two-way nesting simulation of oceanic properties. The yin-yang grid (Kageyama and Sato, 2004), a quasi-uniform overset grid on sphere, has been used in several global atmospheric and oceanic models (Baba, 2010; Qaddouri and Lee, 2011; Zerroukat and Allen, 2012; Li, 2015a; Li and Peng, 2018). Using the overset yin-yang grid helps improve the uniformity of grid spacing and computational efficiency of the models. Interpolation is required in the overset grid to ensure a global computation, and interpolation in the overlapped zone inevitably introduces numerical errors and spurious mass sources or sinks. Therefore, an additional mass corrector or flux constraint on the boundary is necessary for mass conservation in long-term integration. To date, various interpolation schemes, such as spline and Lagrange interpolation (Li, 2006, 2015b), have been used in the numerical models, while they often span several cells for high-order applications. Through the use of multi-moment quantities in a cell, the MCV facilitates the construction of high-order polynomial interpolation functions within a single cell (Li, 2008). Ii(2005) constructed a cubic polynomial that has fourth-order accuracy over a single unstructured mesh by using the surface-integrated average and PV.

Even though high-order local reconstruction reduces numerical errors related to mass conservation in the overlapped grid (Li, 2008), the total mass variation damages the long-term integration of an ocean model. Non-conservative computation results in numerical instability and physically meaningless results in long-term ocean model integrations. Thus, the conservation of physical quantities is critical in computational fluid dynamics, including the atmospheric and oceanic models. The flux-form advection is naturally powerful to guarantee local conservation when discretized under the finite volume definition. Additional processing is required to achieve total mass conservation on an overset grid. Wang (1995) proposed a simple conservation constraint to transform two overlapped zones into a pair of patched zones. The global mass conservation is achieved easily by enforcing one of the two intersecting boundaries to be conserved. Peng(2006) proposed a conservative constraint scheme for the yin-yang grid for the first time. In this scheme, the fluxes on the boundaries of two grid components are constrained identical to ensure global and local mass conservation. This constraint provides a simple solution to the conservative computation of global atmospheric and oceanic models on the yin-yang grid with advantages of economic and efficient flux remapping. Two-way nesting is popularly used in the global and limited-area ocean model for partly high-resolution simulation. Conservative computation is essential to avoid circulation mode shifting in long-term integration of oceanic models.

In this study, we tested the fourth-order MCV-BGS scheme on overset grids in the Cartesian coordinate to simulate idealized 1D and 2D continuous and discontinuous problems, respectively. To ensure global and local conservative transport, we used the boundary flux constraint proposed by Peng(2006) to enforce the identity of mass across boundaries. A fully cubic polynomial is defined to describe the subgrid-scale mass distribution. With the help of a high-order interpolation polynomial, a new interface flux remapping ensures that MCV-BGS is conservative and fourth-order on the overset grid, which is potentially used in nesting oceanic simulation.

The rest of this paper is organized as follows. Section 2 presents a brief introduction to the MCV scheme, BGS limiter, and implementation in addition to the interpolation algorithm in 1D and 2D cases. Then, a 2D cubic polynomial is developed with the help of DOFs of different quantities. A flux corrector is then implemented into the scheme to ensure mass conservation for the conservation law, which modifies the total boundary fluxes to be identical between the overset grids. In Section 3, numerical experiments are discussed to evaluate the accuracy, stability, and convergence speed of the scheme. Finally, Section 4 discusses the results obtained and concludes the paper.

2 Implementation of MCV3-BGS to Overset Grids

2.1 1D MCV3 Scheme

2.1.1 MCV3 algorithm

Finite volume method is popularly used in the oceanic modeling community. We consider the 1D flux-form advection equation

whereis transported quantity, flux,is time, andis spatial position. Initial and boundary conditions are necessary to solve the advection problem.

The three-point MCV scheme (Ii and Xiao, 2009) is adopted in this study. A brief notation of the MCV3 is also given here. Three local DOFs are defined in each control volume, as shown in Fig.1a. In the 1D problem, PVs at the boundaries and cell center are regarded as computational variables and can be used to develop the cell-wise MCV-Lagrangian polynomial

The PV and VIA moments on cellare defined as

The PVs at the cell boundary are updated through a differential form of the governing equations as follows:

wheref?1/2is the spatial derivative of the numerical flux at=x?1/2, and can be obtained by solving the Riemann problem

The VIA is rigorously conservative because of its flux-form computation with the finite volume definition

where ?x=x+1/2?x?1/2and fluxesf+1/2,f?1/2at the cell boundary are computed by the PVs and velocity,,f+1/2=u, 3q, 3,f?1/2=u, 1q, 1. Here, the center point q, 2replaces the VIA and is set as the new calculation variable, and the VIA remains the forecast quantity and is regarded as the constraint condition to updateq, 2.

The following equations can be used to update PVs:

a more detailed numerical procedure can be found in Ii and Xiao (2009).

where A?1/2is the Jacobian matrix at=x?1/2, and[·]is given by

and

An oscillatory-free TVD-type polynomial is defined as

wheredis an adjustable parameter with desired numeri-cal properties, and its specific equation can be found in the work of Deng(2017a). This TVD polynomial is used in the BGS limiter. The dmin and absmin functions are defined as follows:

2.1.2 Temporal discretization

Afive-stagefourth-orderstrongstability-preserving(SSP) Runge-Kutta (RK) method (Spiteri and Ruuth, 2002) with a TVD property is used to achieve high-order accuracy and prevent oscillations during temporal integration processes. The following formulations show the computational procedure of the SSP RK scheme. The coefficients are reported in Table 1.

Table 1 Coefficients of optimal SSPRK (5, 4) scheme

where()=?/?represents the spatial discretization in Eq. (8). At each time step, the flux correction is in the final stage, and any oscillations that are generated in the correction process are suppressed in the later steps.

2.1.3 Configuration of 1D overlapped grids

Fig.2 Schematic of 1D overset grid in a fixed region.

2.2 2D MCV3 Scheme

2.2.1 Brief description of 2D MCV3 algorithm

In a limited-area ocean model, map projection, such as Lambert and Mercator projection, is typically used to project a spherical area to a Cartesian system. The regularity of the Cartesian coordinate simplifies the extension of the numerical scheme from one dimension to multiple dimensions. The 2D transportation is expressed as a hyperbolic system

whereis transported quantity, and flux components,are in theanddirections, respectively.

The distribution of the moments in a 2D problem is shown in Fig.3b. In the 2D case, a dimensional splitting procedure is used to deal with the advection by employing the algorithm developed in the 1D case. Three kinds of nodes are found, namely, points on the cell vertexes, center points of cell interfaces, and center points of cells, which are marked with N1, N2, and N3, respectively. The quantities on N1 nodes (q?1/2,?1/2,q+1/2,?1/2,q?1/2,+1/2,q+1/2,+1/2) can be updated with the following equations:

The symbols withandare the same as those in the 1D case, butindirection.

The quantities on nodes N2,q?1/2,j,q+1/2,jandq,j?1/2,q,j+1/2are updated with different equations. The former two can be updated with

where

Similarly,q,?1/2andq,+1/2are updated with

and

The evolution equation for VIA can be expressed as

where ?=??and1,2,3and4are fluxes at the cell interfaces, which were obtained using the three-point Simpson formula. For details, please refer to Ii and Xiao (2009). Fluxes on the overlapped boundaries are further fixed to guarantee conservation.

Combining Eqs. (18), (20), and (22) with (24), we find that the quantity on N3 node is

2.2.2 Design of 2D overset grid

The idealized 2D transport problem is proposed on a horizontal domain [?1, 1]×[?1, 1] in Cartesian coordinate. To implement an overset grid system into this problem, we first adopt the major uniform grid A to cover the domain with a hole H left (the box with the bold line in Fig.3a) in the center. In another uniform grid B, the minor grid is designed to cover hole H and overlaps with grid A. Grid B is not required to be parallel with grid A, which is a popular nesting grid in the mother grid A. The grid system can be imagined as a simple nesting ocean model with only advection process. The two grids may overlap with an arbitrary cut angle. With respect to the stencils needed in the MCV3 scheme, at least two cells are expended as halo regions (shown using dashed lines in Fig.3a) in the boundaries. The combined grid covers the entire computational domain [?1, 1]×[?1, 1]. In Fig.3a, the transport-computation domain is plotted with solid mesh, while the halo region is shown using dashed cells. All the PVs on the cell center, corner, and midpoint of the cell interface (indicated by dotted lines in Fig.3b) are obtained using the interpolating procedure from another component grid. The transformation of spatial coordinates between grids A and B reads as follows:

whereA,AandB,Bare the coordinates of a fixed point in grids A and B, respectively.=π/6is specified here.

2.2.3 2D cubic interpolation

To achieve fourth-order accuracy interpolation in the 2D overlapped region, we construct a cubic polynomial by fully utilizing the multi-moment quantities in related cells. A cubic interpolation polynomial in cell (,) is first proposed as

wherea(=1, ···, 10) are 10 coefficients to be defined.

As illustrated in Fig.3b, nine quantities can be used to define polynomial (28) in a 2D cell. The nine independent quantitiesare insufficient to determine the 2D cubic polynomial. In addition, only three PVs out of the four in vertex points can be used in the determination process owing to the geometrically isotropic property of the cubic function. Two more constraint conditions are required to construct the perfect cubic polynomial. The spatial derivativesF, jandF, jare naturally suitable candidates. To ensure fourth-order accuracy of the derivative, additional information from the adjacent cells should be used. On the stencils {x?1,x?1/2,x,x+1/2,x+1} and {y?1,y?1/2,y,y+1/2,y+1} in the x and y directions, the fourth-order differencing representation of the derivativesf, jandf, jon the cell center can be written as

Integration of the two derivatives in Eq. (29) and the eight PVs with one at the upper right vertex point rejected in each element enables the determination of the coefficients in Eq. (28). By substituting the 10 moments into Eq. (28), we obtain

whereq, jshows the transported quantities on cell (,), andF, j(x,y) andF, j(x,y) are the derivatives on the cell center. The coefficients are then determined uniquely as follows by solving Eqs. (30)–(39):

and the cubic polynomial (28) is determined with the known quantities on the cell.

2.3 Conservation Constraint for Overset Grid

This study adopts the basic idea of Wang (1995) to transform the overlapped grid into a patched grid. The overset grid in Fig.3a becomes a patched form in Fig.4a, where the boundary of H is marked with Γ. Clearly, H is the common region of A and B; thus, the key point of the conservative constraint is to enforce the fluxes on Γ to be identical in grids A and B. The mass flux is then balanced in the transport process, which guarantees the total mass conservation. As opposed to the conservative constraint proposed by Peng(2006), our study proposes a 2D cubic-function structure of mass and a 1D third-order function of flux on the interfaces instead of the cell-wise constant mass and piece-wise constant flux to improve the constraint accuracy.

Fig.4b shows the intersection of boundaries of grid A with cell interfaces of grid B. One of the cells in grid B cut by Γ is referred to as ABCD, the intersection points are labeled as E and F, and GHQP is one of the virtue hole-boundary cells of grid A. The fluxes through the interfaces of ABCD and GHQP are denoted asAB,BC,CD,ADandGH,HQ,PQ,PG, respectively. The reconstruction of the boundary fluxes to be identical is shown as follows.

and it transforms to the following form by multiplying both sides withAEFD

where(,) is the constructed cubic interpolation polynomial at time step. As shown in Fig.4b, the integrated mass on AEFD is then

whereA,E,F,AandDare the coordinates of A, E, F, and D in grid B.

The fluxEFis then obtained. The cubic polynomials(,) and three-point Simpson formula describes the subgrid structure of the transported quantity and fluxes on the cell interfaces of each of the cut cells. A direct integration of the polynomial, which is the same as the Simp-son formula, increases the computational accuracy compared with the uniform distribution of mass (Peng, 2006) at the cost of a slight increase in the computational time.

4) All the fluxes are integrated through Γ within the cut cells in grid B. The summation is noted as ∑ΓB.

5) The sum of original fluxes across Γ is marked in grid A with ∑ΓA, and substitutes ∑ΓAwith ∑ΓBto make the total flux across Γ identical between the overset grid components. The flux remapping is used to re-compute the boundary-cell VIA quantities. A subtraction of the total flux residual ∑Γ(A?B)on Γ helps to correct the conservative error. The subtraction of flux is in proportion to the total cell incomes (GH?PQ+HQ?GP) on Γ. Thus, the correct flux is written as

where denotes the original flux on Γ. Then, on Γ in Eq. (25) is substituted by in the final stage of Runge-Kutta. Any oscillation that is generated is suppressed in the next stage.

3 Numerical Results

Kuroshio, known as the strongest warm current in the western North Pacific Ocean, forms a thermal vortex and island-induced ocean vortex train behind the Green Island (Hsu, 2015). Zhang(2016) first captured the 3D structures of one cyclonic and two anticyclonic eddies in the northeastern South China Sea west of the Luzon Strait. Due to the existence of general circulation flow (for example, Kuroshio) and meso-scale eddy in the ocean, local sharp variation of temperature and salinity enhance the difficulty of accurate advection computing. Idealized smooth and extremely sharp distribution of the transported quantities is considered in this study. Several widely used numerical experiments are tested in the 1D and 2D overset grids. The results are compared with those obtained on the non-overset grid system to evaluate the con- servative constraint by using the computational accuracy, convergence rate, and total mass variation with time.

3.1 1D Flux-form Transport

3.1.1 Sine wave transport

In the following 1D test cases, the computational domain is set at [?1, 1], and a periodic boundary condition is adopted. The MCV3-BGS is first tested in a 1D case of a smooth sine wave transported by a constant flow on the overset grid. The initial profile is given as

The numerical integration is conducted on gradually refined grids (grid point numbers=20, 40, 80, 160, and 320) with a constant advection velocity=1.0 and CFL=0.1 to check the individual numerical accuracy and convergence rate. The exact solution of this transport problem is known as

The numerical accuracy of the scheme on the overset grid can be estimated by comparing the numerical results with the exact ones. Normalized errors benefit the evaluation of the numerical scheme.1and∞,,

with the resolution increasing, wherecan be either1,2, or∞; and Δ is the grid spacing in the tests. The numerical errors and convergence rates in the two grids are similar, and clearly reach fourth-order accuracy. The∞displays a smaller convergence rate in low-resolution cases due to its special interest on local extremes. The1error is shown as 1.97E?06 and 1.01E?06, the2error is 3.31E?06 and 1.60E?06, and the∞error is 1.10E?05 and 5.21E?06 on the non-overlapped and overlapped grids, respectively, with 80 grid points. The convergence rate on the overlapped grid is slightly lower than that of the non-overlapped one, and this is mainly due to the additional boundary interpolation. The elapse time is also shown in Table 2. The program in the overlapped grid runs slightly slower than that in the non-overlapped grid and the elapsed time is imaginably increased with additional computations.

3.1.2 Advection of a square wave

Besides the smoothly distributed tracers in the transport problem, the discontinuous density issue is popular in fluid dynamics and geophysical issues. The robustness of the advection scheme is a major concern in the numerical computation. The numerical transport of a discontinuous quantity is then tested. The initial profile of the transported square wave is given by

Table 2 Numerical errors and convergence rate in 1D sine wave advection

The computation domain is defined as [?1, 1], and a fixed wind u=1.0. The grid spacing is dx=0.01, and the time step is selected so that CFL=0.3. The numerical re- sults in one round (t=2s) on the overset grid are shown in Fig.5 compared with the true solution. No obvious oscillation appears around the discontinuous points in this case. This condition is highly similar to the results of the non-overlapped grid in Deng et al. (2017a). The normalized errors L1, L2, and L∞ are 1.60E?02, 5.49E?02, and 0.29, respectively, at the end of the one-round advection.

3.1.3 Advection of 1D complex wave

The complex wave is a comprehensive configuration that is employed to test the capability of the MCV3-BGS scheme to simulate continuous and discontinuous problems. This configuration reflects a more difficult condition to deal with the advection on the overset grid. The initial profile is given as

whereandare defined by(,,)=exp(?(?)2) and

respectively, with the coefficients=0.5,=0.7,=0.005,=10.0, and=log2/(362). The computational domain, velocity, resolution, and CFL number are set as the same as those in the square wave experiment.

The numerical solution after one round (=2s) on an overset grid is displayed in Fig.5. The thin square wave is simulated in a manner that is similar to that in the last test, and a slightly diffusive and monotonic result appears at the jumping point even though the wavelength is much shorter. The undershoot on the peak of the sharp wave is small, and the well-simulated smooth wave is impressed on the overset grid. An obvious numerical error is shown in which a sharp gradient exists, and a highly diffusive result is produced by the MCV3-BGS scheme to suppress the numerical oscillations on the 1D overset grid. The error norms1,2, and∞are calculated as 4.34E?02, 7.82E?02, and 0.29, respectively, which confirms the flexibility and robustness of the scheme to describe the complex wave with strong discontinuity on the overset grid.

3.2 2D Advection

3.2.1 2D sine wave transportation

The MCV3-BGS scheme was also tested in the 2D advection problem on the overlapped grid. The computational domain is [?1, 1]×[?1, 1] with periodic boundary conditions in both directions on grid A. The most important factor is the numerical accuracy in the real computation. To evaluate the performance of the 2D cubic interpolation function and the impact of the conservative constraint scheme that was implemented, we conducted the smooth 2D sine wave advection to simulate passive advection of continuous distribution in the ocean at different resolutions using grid refinement. The1,2, and∞errors are estimated on the overlapped grid with and without the conservative constraint. A careful comparison of the results with those on the non-overlapped grid helps to show the suitability of the developed scheme on the overlapped grid.

In the test case, the initial condition is defined as

and the exact solution is known as

A steady advection velocity (,)=(1, 1) is appointed throughout the grid A. The velocity on grid B is given using the coordinate transformation Eq. (27). The entire computational domain is defined on [?1, 1]×[?1, 1] with a hole [?0.2, 0.2]×[?0.2, 0.2] in grid A. Grid B covers a region of [?0.3, 0.3]×[?0.3, 0.3]. A series of numerical advections are conducted with gradual refinement of grids. The grid point numbers areA=10×10, 20×20, 40×40, 80×80, 160×160 in grid A andB=3×3, 6×6, 12×12, 24×24, 48×48 in grid B. The CFL numbers are 0.1 in the tests. After being transported for one round, the normalized errors and corresponding convergence rates with/without a conservative constraint on the overlapped grid compared with that on the non-overlapping grid at=2s are listed in Table 3. Compared with the results on the non-overlapped grid, the1,2and∞errors on the overlapping grid increase slightly due to additional computation. The details can be observed in Table 3. The convergence rates in all three cases reach approximately the fourth-order. In the absence of conservative constraint, more errors appear in the overlapped grid case. The results imply that high-quality advection is possibly achieved on the overset grid with the help of fourth-order boundary interpolation. The compact high-order interpolation scheme also benefits the advection of a sharp gradient quantity on the overlapped grid. Similar to the 1D cases, the computing time of the overlapped grid is 30% longer than that of the non-overlapped grid in a 160×160 grid case. However, in the overlapped grid cases, the elapse time with and without conservation constraintare quite close. In a high-resolution computation, the portion of constraint on the overset grid boundary becomes extremely small with respect to the transport. The increase of CPU time is not a serious issue for high-resolution simulations on the overset grid.

The flux-form advection guarantees mass conservation in the finite volume definition. However, the boundary interpolation destroys the conservation in overset grids because of the change in the cell-integrated value. The application of the conservative constraint in the last section aims to patch the hole in grid A with the corresponding mass on grid B. In addition, the 2D cubic interpola- tion polynomial that is used to determine the mass on cut cells naturally simulates the density function within a cell, and it therefore ensures the accurate mass integration of the cut cell. In Fig.6, we illustrate the temporal variation of the total mass during one-round transportation in the case where the grid number is 80×80 in grid A and 24×24 in grid B. The exact conservation of the total mass is clearly shown in the computation, where the total mass is a summation of the cell integration on the patched grid (Fig.4a).

Table 3 Numerical errors and convergence rate in 2D advection problem.

Fig.6 Evolution of total mass in a one-round transport. Dot line shows the error relative to the initial total mass.

3.2.2 Slotted-cylinder problem

The slotted-cylinder problem (Zhang and Juang, 2013) with a large spatial gradient is widely used to verify numerical advection schemes to capture the jump discontinuities. This approach is used to test the scheme in case of extremely sharp variation in a small scale of oceanic quantities, such as horizontal variation of temperature or salinity. Schemes with obvious dissipative errors fail to simulate the spatial structure, which shows diffusive or oscillatory results. The initial condition is

where the constant0is the height of the cylinder,is the cylinder radius, ands,sare the width and length of the slot, respectively.,are the coordinates relative to the center (x?cos(),y?sin()), anddenotes the distance between the centers of the computational domain and cylinder. The parameters,, andare defined as

Here, dimensionless quantities are used. We set the coefficients as=1,0=1,=20,s=10,s=25, and=25. The coordinates of the flow center are (50, 50). The computational domain is prescribed on [?1, 1]×[?1, 1] with a hole [?0.35, 0.35]×[?0.35, 0.35] in grid A. Grid B spans a domain of [?0.5, 0.5]×[?0.5, 0.5]. The rotating velocity is (,)=(?2π, 2π). CFL=0.3 is selected for the computation. The slotted cylinder is rotated for one round after=1s integration on the uniform overset grid, with 100×100grid points in grid A and 50×50grid points in grid B. The numerical results in one-round rotation of the cylinder are shown in Fig.7. A clear cylinder is displayed in the numerical test, and spurious oscillations are suppressed effectively in the computation. Sharp edges are mainly captured by the numerical scheme, which benefits from the local high-order interpolation and remapping procedure. The high-order interpolation scheme with short stencil may capture the sudden jump in the transported quantity. In addition, the mass redistribution related to the cell-wise fourth-order function in the conservative constraint contributes largely to the high-order transportation. In the test cases, the normalized errors1,2, and∞are 2.38E?02, 3.29E?02, and 0.12, respectively, which generally reflect the diffusive errors at the jumping points.

Fig.7 (a) Initial condition of dimensionless slotted-cylinder problem and numerical solutions (b) without conservative constraint and (c) under the conservative constraint after one-round advection at t=1s on the overset grid. A total of 100×100grid points are in grid A and 50×50 are in grid B.

3.2.3 Advection of 2D complex wave

The MCV3-BGS scheme on the overset grid is further verified using the numerical example of 2D non-dimensional complex waves that extend from the 1D problem in Jiang and Shu (1996) to test the performance in the case of very strong discontinuities. The initial condition is given as follows:

where

The coefficients are set as=0.01,=5, and=log2/(362). The grid spacing, rotation velocity, and CFL number are the same as in the previous test. After one-round rotation, the results at=1s are shown in Fig.8. Compared with the initial condition, all of the wave structures are effectively simulated without any visible spurious undershoot and overshoot around the discontinuities. The waves that are simulated on the overlapped and non-overlapped grids are generally the same quality. The cubic interpolation based on the multi-moment quantities produces similar results in the overlapped region as the non-overlapped mesh. The error norms show a similar magnitude as in the slotted-cylinder test. The1,2, and∞norms are 2.24E?02, 3.26E?02, and 0.11, respectively, in the overlapped and non-overlapped grids.

Fig.8 Same as Fig.7 but for dimensionless complex wave.

4 Discussion and Conclusions

The fourth-order non-oscillatory MCV3-BGS advection scheme is used on overset grids for transport computation in the 1D and 2D cases, thereby showing the potential for use in ocean models for nesting simulation. In the overset grids, boundary interpolation is required for the smooth integration of the transport equation. Boundary interpolation shares the cubic polynomial with the MCV3 scheme in the 1D case, which was constructed to determine the cell-interface flux. However, in 2D cases, a newly built interpolation polynomial is required. With the help of multi-moment quantities, including existing DOFs and newly defined derivatives on one cell, a cubic polynomial was proposed for 2D interpolation. The use of the same order of interpolation accuracy as that of the advection scheme guarantees the accurate transport computation on an overset grid system. The present fourth-order conservative constraint can be expended to the yin-yang grid for global simulation of oceanic circulation.

Both 1D and 2D test cases show approximately fourth-order transportation with the MCV3-BGS scheme on the overset grids, which is confirmed using the smooth sine wave advection problem. A direct cubic interpolation on boundaries ensures similar numerical results to those on non-overlapped grids, where the1norm indicates the identical general error1in both grids. The application of boundary interpolation and conservative constraints introduce local errors into the computation, and these errors significantly affect the∞norm. Although the conservative constraint makes the∞convergence rate of the advection scheme decrease slightly, an approximate fourth-order accuracy is maintained, with the total mass being conservative. The BGS limiter effectively suppresses numerical oscillations where strong discontinuity exists; for example, the diffusion at the jump points in the square wave, slotted-cylinder problem, as well as complex wave cases. This monotonic property is attractive for the positive tracer advection in atmospheric and oceanic simulations. Conservative positive-defined advection is achieved with the MCV3-BGS scheme on the overset grids in 1D and 2D Cartesian coordinates, thereby showing a way for direct high-order computation of transport in limited-area ocean models.

With the help of the cubic polynomial, the conservative constraint shows good results without significantly decreasing the accuracy of advection. The integration of the cubic polynomial exhibits the subgrid-scale features of the transported quantities, thereby producing highly realistic advection results. The presence of the multi-moment quantities also helps to refine the interface flux computation accurately. Finally, the high-order remapping of the mass and flux leads to high-quality achievements on over- lapped grids, thereby improving the nesting simulation of oceanic properties.

Acknowledgements

The authors thank Dr. X. Deng at the Tokyo Institute of Technology for the useful discussions on the MCV scheme. We also appreciate the helpful input by Dr. X. L. Li at the China Meteorological Administration. This study was supported by grants from the National Natural Science Foundation of China (Nos. 41575103 and 91637210).

Baba, Y., Takahashi, K., and Sugimura, T., 2010. Dynamical core of an atmospheric general circulation model on a yin-yang grid., 138 (10): 3988-4005, DOI: 10.1175/2010MWR3375.1.

Bastian, P., 2014. A fully-coupled discontinuous Galerkin method for two-phase flow in porous media with discontinuous capillary pressure., 18 (5): 779-796, DOI: 10.1007/s10596-014-9426-y.

Chen, C. G., and Xiao, F., 2008. Shallow water model on cubed-sphere by multi-moment finite volume method., 227 (10): 5019-5044, DOI: 10.1016/j.jcp.2008.01.033.

Chen, C. G., Li, X. L., Shen, X. S., and Xiao, F., 2014. Global shallow water models based on multi-moment constrained finite volume method and three quasi-uniform spherical grids., 271: 191-223, DOI: 10.1016/j.jcp.2013.10.026.

Deng, X., Sun, Z. Y., Xie, B., Yokoi, K., Chen, C. G., and Xiao, F.,2017a. A non-oscillatory multi-moment finite volume schemewith boundary gradient switching.,72: 1-23, DOI: 10.1007/s10915-017-0392-0.

Deng, X., Xie, B., and Xiao, F., 2017b. A finite volume multi-moment method with boundary variation diminishing principle for Euler equation on three-dimensional hybrid unstructured grids., 153: 85-101, DOI: 10.1016/j.compfluid.2017.05.007.

Hsu, T. W., Doong, D. J., Hsieh, K. J., and Liang, S. J., 2015. Numerical study of monsoon effect on green island wake., 31 (5): 1141-1150, DOI: 10.2112/JCOASTRES-D-14-00206.1.

Ii, S., and Xiao, F., 2007. CIP/multi-moment finite volume method for Euler equations, a semi-Lagrangian characteristic formulation., 222 (2): 849-871, DOI: 10.1016/j.jcp.2006.08.015.

Ii, S., and Xiao, F., 2009. High order multi-moment constrained finite volume method. Part I: Basic formulation., 228 (10): 3669-3707, DOI: 10.1016/j.jcp.2009.02.009.

Ii, S., and Xiao, F., 2010. A global shallow water model using high order multi-moment constrained finite volume method and icosahedral grid., 229 (5): 1774-1796, DOI: 10.1016/j.jcp.2009.11.008.

Ii, S., Shimuta, M., and Xiao, F., 2005. A 4th-order and single-cell-based advection scheme on unstructured grids using multi-moments., 173 (1-2): 17-33, DOI: 10.1016/j.cpc.2005.07.003.

Imai, Y., Aoki, T., and Takizawa, K., 2008. Conservative form of interpolated differential operator scheme for compressible and incompressible fluid dynamics., 227 (4): 2263-2285, DOI: 10.1016/j.jcp.2007.11.031.

Jiang, G. S., and Shu, C. W., 1996. Efficient implementation of weighted ENO schemes., 126: 202-228, DOI: 10.1006/jcph.1996.0130.

Jin, P., Deng, X., and Xiao, F., 2018. A direct ALE multi-moment finite volume scheme for the compressible Euler equations., 24 (5): 1300-1325, DOI: 10.4208/cicp.OA-2017-0189.

Kageyama, A., and Sato, T., 2004. The ‘yin-yang grid’: An overset grid in spherical geometry., 5 (9): Q09005, DOI: 10.1029/2004GC000734.

Li, X. H., and Peng, X. D., 2018. Long-term integration of a global non-hydrostatic atmospheric model on an aqua planet., 32 (4): 517-533, DOI: CNKI:SUN:QXXW.0.2018-04-001.

Li, X. H., Peng, X. D., and Li, X. L., 2015a. An improved dynamic core for a non-hydrostatic model system on the yin-yang grid., 32 (5): 648-658, DOI: 10.1007/s00376-014-4120-5.

Li, X. L., Chen, D. H., Peng, X. D., Takahashi, K., and Xiao, F., 2008. A multimoment finite-volume shallow-water model on the yin yang overset spherical grid., 136 (8): 3066, DOI: 10.1175/2007mwr2206.1.

Li, X. L., Chen, D. H., Peng, X. D., Xiao, F., and Shen, X. S., 2006.Implementationofthesemi-Lagrangianadvectionscheme on a quasi-uniform overset grid on a sphere., 23 (5): 792-801, DOI: 10.1007/s00376-006-0792-9.

Li, X. L., Chen, C. G., Xiao, F., and Shen, X. S., 2015b. A high-order multi-moment constrained finite-volume global shallow-water model on the yin-yang grid., 141 (691): 2090-2102, DOI: 10.1002/qj.2504.

Liu, X. D., Osher, S., and Chan, T., 1994. Weighted essentially non-oscillatory schemes., 115 (1): 200-212, DOI: 10.1006/jcph.1994.1187.

Onodera, N., Aoki, T., and Kobayashi, H., 2011. Large-eddy simulation of turbulent channel flows with conservative IDO scheme., 230: 5787-5805, DOI: 10.1016/j.jcp.2011.04.004.

Peng, X. D., Xiao, F., and Takahashi, K., 2006. Conservative constraint for a quasi-uniform overset grid on the sphere., 132 (616): 979-996, DOI: 10.1256/qj.05.18.

Peng, X. D., Xiao, F., Takahashi, K., and Yabe, T., 2004. Conservative CIP transport in meteorological models., 47 (4): 725-734, DOI: 10.1299/jsmeb.47.725.

Qaddouri, A., and Lee, V., 2011. The Canadian Global Environmental Multiscale model on the Yin-Yang grid system., 137 (660): 1913-1926, DOI: 10.1002/qj.873.

Spiteri, R. J., and Ruuth, S. J., 2002. A new class of optimal high-order strong-stability-preserving time discretization methods., 40 (2): 469-491, DOI: 10.1137/S0036142901389025.

Sun, Y. Z., Wang, Z. J., and Liu, Y., 2007. High-order multidomain spectral difference method for the Navier-Stokes equations on unstructured hexahedral grids., 2 (2): 310-333, DOI: 10.2514/6.2006-301.

Tanaka, R., Nakamura, T., and Yabe, T., 2000. Constructing exactly conservative scheme in a non-conservative form., 126 (3): 232-243, DOI: 10.1016/s0010-4655(99)00473-7.

Taneja, A., and Higdon, J., 2018. A fully-coupled discontinuous Galerkin spectral element method for two-phase flow in petroleum reservoirs., 352: 341-372, DOI: 10.1016/j.jcp.2017.09.059.

Wang, Z. J., 1995. A fully conservative interface algorithm for overlapped grids., 122 (1): 96-106, DOI: 10.1006/jcph.1995.1199.

Wang, Z. J., 2004. Spectral (finite) volume method for conservation laws on unstructured grids: Basic formulation., 194 (2): 716-741, DOI: 10.1006/jcph.2002.7041.

Xiao, F., and Yabe, T., 2001. Completely conservative and oscillationless semi-Lagrangian schemes for advection transportation., 170 (2): 498-522, DOI: 10.1006/jcph.2001.6746.

Xiao, F., Yabe, T., Peng, X. D., and Kobayashi H., 2002. Conservative and oscillation-less atmospheric transport schemes based on rational functions., 107 (D22): 4609, DOI: 10.1029/2001jd001532.

Xie, B., and Xiao, F., 2014. Two and three dimensional multi-moment finite volume solver for incompressible Navier-Stokes equations on unstructured grids with arbitrary quadrilateral and hexahedral elements., 104: 40-54, DOI: 10.1016/j.compfluid.2014.08.002.

Xie, B., Ii, S., Ikebata, A., and Xiao, F., 2014. A multi-moment finite volume method for incompressible Navier-Stokes equations on unstructured grids: Volume-average/point-value formulation., 277: 138-162, DOI: 10.1016/j.jcp.2014.08.011.

Yabe, T., Tanaka, R., Nakamura, T., and Xiao, F., 2001. An exactly conservative semi-Lagrangian scheme (CIP-CSL) in one dimension., 129 (2): 332-344, DOI: 10.1175/1520-0493(2001)1292.0.CO;2.

Zerroukat, M., and Allen, T., 2012. On the solution of elliptic problems on overset/yin-yang grids., 140 (8): 2756-2767, DOI: 10.1175/MWR-D-11-00272.1.

Zhang, S. Q., Liu, Y., Ma, X. H., Wang, H. N., Zhang, X. F., Yu, X. L., and Lu, L., 2018. The ‘two oceans and one sea’ extended range numerical prediction system with an ultra-high resolution atmosphere-ocean-land regional coupled model., 11: 1674-2834.

Zhang, Y. F., and Juang, H. M. H., 2013. A mass-conserving non-iteration-dimensional-split semi-Lagrangian advection scheme for limited-area modelling., 138 (669): 2118-2125, DOI: 10.1002/qj.1938.

Zhang, Z., Tian, J. W., Qiu, B., Zhao, W., Chang, P., Wu, D. X., and Wan, X. Q., 2016. Observed 3D structure, generation, and dissipation of oceanic mesoscale eddies in the South China Sea., 6 (1): 24349, DOI: 10.1038/srep24349.

. Tel: 0086-10-68409552

E-mail: pengxd@cma.gov.cn

August 13, 2019;

September 30, 2019;

October 12, 2019

(Edited by Xie Jun)