A Two-Way Nesting Unstructured Quadrilateral Grid, Finite-Differencing, Estuarine and Coastal Ocean Model with High-Order Interpolation Schemes

: The balance between the need of improving horizontal resolution in simulating local small-scale ocean processes and computational costs makes it desirable to reﬁne model mesh locally. A three-dimensional, two-way nesting unstructured quadrilateral grid, primitive equations, ﬁnite-differencing, estuarine and coastal ocean model is developed for multi-scale modeling. Because the model grid is capable of multi-area nesting and multi-level reﬁnement at each subdomain, the model is highly compatible with simulations involved in complex topography and studies of local small-scale ocean processes. The two-way information exchange is achieved by a virtual grid method, and its basic idea is to implement numerical integrations of variables at nesting interfaces with the support of virtual grid variables, which are interpolated or updated from actual grid variables. The model is novel for two interpolation schemes: the high-order spatial interpolation at the middle temporal level (HSIMT) parabolic interpolation scheme and HSIMT advection-equivalent interpolation scheme, and they have high-order accuracy and good consistency with the advection scheme applied to solving the tracer equations. The conservation of both volume and tracer contents is ensured via a ﬂux correction algorithm. The two original interpolation schemes are examined in an ideal salinity advection experiment in the peak preservation skill, stability, and conservation properties. A realistic application to the Deep Waterway Project area in the Changjiang Estuary showed that the nested grid model can reproduce the hydrodynamic processes at the observed sites successfully while it failed to maintain the performance with the structured grid model in simulating the variance of salinity, for which the enforced conservation had primary responsibility. The HSIMT parabolic interpolation scheme was distinguished from other schemes for its outstanding performances in simulating salinity.


Introduction
The balance between the need of improving horizontal resolution in simulating local small-scale ocean processes and computational costs makes it desirable to refine model mesh locally. Distinguished by horizontal discretization, there are basically two approaches to realize this idea. One natural and trendy way is to use unstructured mesh models. The first work that applied finite elements and unstructured meshes to ocean modeling probably dates back to 1975 [1]. A lot of unstructured mesh ocean models have been developed over the last several decades. There are some older names like ADCIRC [2,3], QUODDY [4], UnTRIM [5], and TELEMAC [6,7]. Plenty of models, such as FVCOM [8], FE-SOM [9,10], ICOM/Fluidity [11][12][13], SLIM [14][15][16][17], MIKE 21 & MIKE 3 [18], SUNTANS [19], SCHISM [20], MPAS [21,22], and GOM [23], featuring improvements in many different of cells), which also means that the neighborhood information of each cell is required to be stored for the high-order finite-differencing numerical scheme. This numbering system is totally different from that used by common structured-grid models where all the cells are numbered with two indices (i and j) representing two horizontal directions, while similar with that applied by some unstructured-grid models instead. Moreover, that is why, we call it unstructured quadrilateral grid. The nesting is in two-way interaction implemented by the virtual-grid method, which will be detailed in the following section of this paper. The model utilized an advanced and robust HSIMT-TVD (high-order spatial interpolation at the middle temporal level coupled with a TVD limiter) scheme to solve the advection terms in the tracer equations [41]. Two novel interpolation schemes were tested in UFDECOM: HSIMT parabolic interpolation (HPI) scheme and HSIMT advectionequivalent interpolation (AEI) scheme. They are developed both based on the special advection scheme. The HPI scheme directly apply the upwind parabola function which is used to describe the scalar distribution in the scheme, while the HSIMT AEI scheme is constructed from the HSIMT-TVD formula via the idea of AEI scheme. They both have high-order accuracy and good consistency with the advection scheme, which provides an important advantage for UFDECOM compared to other nested grid models. Overall, there is no precedent for an ocean model featuring all these techniques.
There are basically two approaches for the treatment of the fast mode associated with external gravity waves. The first and commonly used one is the split-explicit time integration method of Blumberg and Mellor [42]. In this case, the barotropic mode and the baroclinic mode are integrated separately with different time steps. The other way is an implicit time-stepping method which breaks the limits of the Courant-Friedrich-Lewy (CFL) stability criterion for the external mode, at the cost of solving an elliptic system [43,44]. Due to the grid topography and numbering approach of UFDECOM, the implicit time-stepping method is no longer available, and there lie the consistency issues in the mode-splitting method [43], so a fully explicit time integration scheme is adopted except for the vertical eddy viscosity and diffusion terms, with the optional OpenMP parallelization to relieve the computational pressure. However, the mode-splitting option can be added to UFDECOM to face the need for faster simulation in the near future.
The rest of this paper is organized as follows: in Section 2, we focus on the technical issues of UFDECOM, such as fundamental model characteristics, algorithms of nesting procedure, and embedding techniques. An ideal salinity advection experiment is performed to evaluate different interpolation schemes in Section 3. For validation and assessment of novel interpolation schemes (Section 4), the new model is applied to a real case in the Changjiang Estuary and compare the model results with that of a structured-grid model and observational data. Conclusions are listed in Section 5.

Basic Model Characteristics
Under the hydrostatic assumption and the Boussinesq approximation, we introduce a horizontal non-orthogonal curvilinear coordinate transformation (ξ = ξ(x, y), η = η(x, y)) and a vertical σ coordinate transformation (σ = z−ζ H+ζ , σ varies from −1 at z = −H to 0 at z = ζ) to the primitive equations to form the basis of this model. The governing equations are first derived by Chen et al. [45]. The momentum equations are given as and the continuity equation is ∂ζ ∂t where The conservation equation for salinity and temperature can be written as

∂JDS ∂t
∂JDθ ∂t The density equation is where In the above equations, H is the water depth, and D = H + ζ with ζ, the surface elevation. f is the Coriolis parameter, and g is the gravitational acceleration. J is the Jacobian function, which can be calculated by J = x ξ y η − x η y ξ , where the subscripts indicate partial derivatives. h 1 , h 2 , and h 3 are defined as u 1 and v 1 are the ξ and η components of the velocity, respectively, and the relationship between them and the x and y components (u and v) is U andV are defined asÛ In Equations (5) and (6), S is the salinity and θ is the potential temperature; ρ and ρ 0 are the perturbation and reference densities; the vertical eddy viscosity coefficient and the thermal vertical eddy diffusion coefficient are denoted as K m and K h , and the modified Mellor-Yamada level 2.5 turbulence closure scheme [46,47] is applied to calculate them; F represents the horizontal diffusion terms, which can be expressed as where A m and A h are the horizontal counterparts of K m and K h , respectively, and are obtained by Smagorinsky's formula [48].
The initial elevation and water current are usually set to zero, and the initial salinity and temperature filed are interpolated from climatological or observational data.
The surface and bottom boundary conditions for the momentum, salinity, and temperature equations are given by ω σ=0 or−1 = 0; In Equation (17), τ 0ξ , τ 0η are the ξ and η components of surface wind stresses, which are calculated with the neutral steady-state drag coefficient [49]. τ bξ , τ bη are the ξ and η components of bottom stresses, and they can be obtained with the following formula: where C d is the bottom drag coefficient which can be determined by the logarithmic law [50]. P net and Q net are the surface net heat and salinity flux, respectively. The lateral boundary conditions for closed basins are specified as where t and n represent the tangential and normal direction to the boundary, respectively. The open boundary conditions are yielded by measured data or radiation boundary conditions.
The governing equations are discretized using the finite difference method on an Arakawa C-grid [51] (Figure 1). The fractional-step method [52] is applied to the time integration of the discretized momentum equations. In the first step, all the terms except for the vertical diffusion terms are treated explicitly with the simple Euler scheme. In the second step, the implicit treatment of the vertical diffusion terms leads to a tri-diagonal matrix, which can be solved by the tri-diagonal matrix algorithm [53], and after that, the horizontal velocity field is obtained. The water elevation is calculated by substituting the n+1 step u and v into the continuity equation. An HSIMT-TVD (high-order spatial interpolation at the middle temporal level coupled with a TVD limiter) scheme with highorder accuracy and non-oscillation is utilized to solve the advection terms in the tracer equations [41].

Figure 1.
A coarse grid cell divided into nine fine grid cells with the refinement ratio of 3 in both horizontal directions and variable distributions on a C-grid. The coarse-resolution grid (CRG) variables and high-resolution grid (HRG) variables are represented by the hollow and solid symbols, respectively.

Implementation of Nesting
For most of the nested structure-grid models, the nesting procedure is implemented by running multiple models [40], which means the HRG model is driven by the boundary conditions interpolated from the CRG solution and meanwhile, the CRG model is affected by the HRG solution through updating at the nesting interfaces. Associated with mesh refinement, time refinement is usually adopted. The time integration of different-resolution-mesh models will be coupled like the treatment of split-explicit models to the barotropic mode and baroclinic mode. Hence, the nesting algorithm can be simplified as follows: UFDECOM does not adopt any time-splitting method, and our aim is to run one multi-scale model for better consistency. The grid can be as concise as those used by unstructured-grid models ( Figure 2). However, for the numerical integration on one certain grid cell at the nesting interfaces, values of neighboring cells cannot be directly used due to different spatial resolutions. Thus, the virtual-grid method is devised to solve this problem. Figure 2 is a schematic of the virtual-grid method. To put it simply, we set some extra grids next to the nesting interfaces in both horizontal directions. The number of circles in one virtual grid depends on the numerical schemes (it is 4 for the HSIMT-TVD advection scheme). Variables at these grids are obtained not by calculating but by interpolating or updating, so we name them virtual grids. When a virtual grid is used by a finer actual grid, it is an external virtual grid. Otherwise, it is internal. In the numerical integration, if some variable is stepped to the nest time step, the values of virtual grids are interpolated/updated immediately by an information-exchanging subroutine. Therefore, values of actual and virtual grids are renewed almost synchronously and solutions can be obtained in the whole model domain.

Implementation of Nesting
For most of the nested structure-grid models, the nesting procedure is implemented by running multiple models [40], which means the HRG model is driven by the boundary conditions interpolated from the CRG solution and meanwhile, the CRG model is affected by the HRG solution through updating at the nesting interfaces. Associated with mesh refinement, time refinement is usually adopted. The time integration of different-resolutionmesh models will be coupled like the treatment of split-explicit models to the barotropic mode and baroclinic mode. Hence, the nesting algorithm can be simplified as follows:

1.
Step the CRG model forward to the next time step.

2.
Interpolate the HRG boundary conditions from the CRG model values both spatially and temporally. The time refinement and spatial refinement usually share the same ratio.

3.
Step the HRG model forward to the same physical time as the CRG model.

4.
Update the values of the overlapped CRG points with the new HRG model values.
UFDECOM does not adopt any time-splitting method, and our aim is to run one multi-scale model for better consistency. The grid can be as concise as those used by unstructured-grid models ( Figure 2). However, for the numerical integration on one certain grid cell at the nesting interfaces, values of neighboring cells cannot be directly used due to different spatial resolutions. Thus, the virtual-grid method is devised to solve this problem. Figure 2 is a schematic of the virtual-grid method. To put it simply, we set some extra grids next to the nesting interfaces in both horizontal directions. The number of circles in one virtual grid depends on the numerical schemes (it is 4 for the HSIMT-TVD advection scheme). Variables at these grids are obtained not by calculating but by interpolating or updating, so we name them virtual grids. When a virtual grid is used by a finer actual grid, it is an external virtual grid. Otherwise, it is internal. In the numerical integration, if some variable is stepped to the nest time step, the values of virtual grids are interpolated/updated immediately by an information-exchanging subroutine. Therefore, values of actual and virtual grids are renewed almost synchronously and solutions can be obtained in the whole model domain. Combining the numerical scheme and the virtual-grid method, the general algorithm of UFDECOM is illustrated in a flow diagram ( Figure 3). The information-exchanging subroutine is inserted between any two independent processes. Combining the numerical scheme and the virtual-grid method, the general algorithm of UFDECOM is illustrated in a flow diagram ( Figure 3). The information-exchanging subroutine is inserted between any two independent processes. J. Mar. Sci. Eng. 2021, 9, x FOR PEER REVIEW 8 of 25

Embedding Techniques
Obviously, update and interpolation schemes are the most crucial techniques for nested grid models. However, the embedding procedure of updating and interpolating may cause some other problems, which can affect the accuracy, stability, and conservation of solution, and requires careful treatment.

Embedding Techniques
Obviously, update and interpolation schemes are the most crucial techniques for nested grid models. However, the embedding procedure of updating and interpolating may cause some other problems, which can affect the accuracy, stability, and conservation of solution, and requires careful treatment.

Update Schemes
Since update is feeding the HRG solutions back to the CRG variables, a natural principle will be maximizing the information well revolved on the CRG and damp the small scales as strongly as possible at the same time.
We tested five update schemes in UFDECOM. Here is a brief description.

1.
Directly-replacing (DR) scheme In this scheme, all internal-virtual-grid (IVG) variables are directly replaced by the values of their respective nearest actual-grid (AG) points. This scheme is extremely simple and holds no smoothing or conservative characteristics, which shows great instability in our ideal experiments.

2.
Inverse distance weighting interpolation (IDWI) scheme A certain IVG point is taken as the origin, and we search for its closest AG point in four quadrants. Then, the value is calculated with a weighted average of these four points. This scheme can be expressed as where φ IVG is the value at the IVG point and φ AG i is the value at the closest AG point in quadrant i. m i is the free surface mask, which is updated at each time step. The weight is denoted as w i , which is determined by the distance d i between φ AG i and φ IVG .

1.
Area-averaging (AVE) scheme The coarse IVG value in the cell is replaced by the area-weighted sum of the AG values in the same cell where A is the area of the fine AG cell and n is the total number of AG cells in the coarse cell.

9-point Shapiro filtering (SF) scheme
The 9-point Shapiro filtering scheme is derived by Shapiro [54]. The transfer function in 2D is presented here: where S is the smoothing element, which is 0.5 in UFDECOM.

Full-weighting operator (FWO) scheme
Debreu, Marchesiello, Penven, and Cambon [36] discovered that using a larger stencil can construct higher-order restriction operators, so they expanded the stencil to 25 points and derive the full-weighting operator scheme. The 2D transfer function is given as where S = 2/3. Among all the update schemes above, the directly-replacing scheme and inverse distance weighting interpolation scheme have low-order accuracy and are coded only for possible tests at the early developing stage, while the area-averaging scheme and 9-point Shapiro filtering scheme have been widely used and tested by nested grid models [33][34][35][36]55,56]. Debreu, Marchesiello, Penven, and Cambon [36] used a larger stencil to develop the high-order full-weighting operator scheme, which were designed to properly damp subgrid scale features. On the contrary, the area-averaging scheme and the Shapiro-filtering scheme require artificial increasing diffusion near the nesting interfaces (sponge layer) to damp the small-scale modes. Ideal experiments indicated that the fullweighting operator could yield a significant improvement over the solution based on the area-averaging scheme [36]. So further experiments in this work adopted the full-weighting operator scheme as the update scheme considering previous research.

Interpolation Schemes
In contrast to update schemes, plenty of interpolation schemes have been used in nested grid models [33]. Modelers have been seeking a high-order interpolation scheme, which is meanwhile conservative. High-order interpolation schemes may produce spurious oscillations especially in areas of sharp gradients [57]. UFDECOM has seven available options.

1.
Zeroth-order uniform interpolation (UNI) scheme All external virtual grid (EVG) points within an AG cell are assigned the same value of that cell.

2.
Inverse distance weighting interpolation (IDWI) scheme It is the same scheme as the IDWI update scheme.

Inverse bilinear interpolation (IBI) scheme
In this scheme, an EVG point value is interpolated in a quadrilateral defined by four surrounding AG points. We choose the inverse bilinear interpolation scheme, which applies to a random quadrilateral, in case that the nested grid has an even refinement ratio.

4.
Quadratic interpolation (QI) scheme Clark and Farley [55] developed this scheme based on the conservation relation suggested by Kurihara et al. [58]. One is referred to Clark and Farley [55] for detailed expressions.

HSIMT parabolic interpolation (HPI) scheme
A key issue in the construction of a flux-based advection scheme is to determine the flux at the cell interface, while this problem reduces to estimating the concentration of some specific substance (φ) for a C-grid based model. In the HSIMT-TVD scheme, an upwind parabola function is used to describe the distribution of φ, which inspires us to apply the same method to interpolating the values of external virtual points. This idea can maintain the consistency of variable distribution for both actual grids and virtual grids. Following the expressions of Wu and Zhu [41], the local distribution is in the form of where ξ = ε + 1/2 and ε is the normalized local coordinate, which varies from −1 on φ AG the three coefficients will be as follows Under the opposite condition, the expressions will be The basic idea of the HPI scheme is actually originated from the piecewise parabolic assumption in the piecewise parabolic method (PPM) [59]. The PPM is a higher-order extension of Godunov's method [60,61], which is of a type first introduced by van Leer in his MUSCL algorithm [62]. The MUSCL algorithm was devised to construct a numerical scheme for solving the advection terms in the trace equations, while the HSIMT parabolic scheme is developed for information exchange in the nested grid model. Thus, the MUSCL algorithm and the HSIMT parabolic scheme serve for totally different purposes. Generally speaking, the HSIMT parabolic scheme can be seen as a transplantation of a high-order version the MUSCL approach to the nested grid model. However, nested grid model with such upwind piecewise parabolic interpolation scheme has been unseen in previous research.

6.
Upwind advection-equivalent interpolation (AEI) scheme The basic idea of the AEI scheme is to design a shape-preserving interpolator from a variety of monotone advection algorithms [63]. Thus, a formal equivalence of the advection scheme and interpolator can be obtained on discrete meshes. Starting from an advection scheme formula, the local Courant number vector is replaced by a normalized local coordinate whose origin coincides with the center of a coarse AG cell (for scalar variables). The accuracy of the derived interpolation scheme is maintained the same as that of the advection scheme utilized.
The upwind AEI scheme here is a little different from that described by Alapaty, Mathur, and Odman [57]. The original AEI scheme based on the upwind advection scheme is actually the same as the linear interpolation. We make it upwind physically by retaining the velocity sign in the formula. The scheme in 1D can be expressed as

HSIMT advection-equivalent interpolation scheme
Since the advection scheme adopted by UFDECOM is the HSIMT-TVD scheme, a natural idea is to develop a desirable AEI scheme from this robust advection scheme. One can refer to the paper of Wu and Zhu [41] for the derivation and details of this flux-based advection scheme. We omit the derivation process and directly give the final formulae. The HSIMT AEI scheme can be obtained by replacing −u∆t/∆X by ε to yield the following where otherwise. In Equations (35) and (36),

Conservation
If numerical schemes are written in flux form, conservation is maintained in a structured grid model. However, this property vanishes for two-way interactive nested models because of the unavoidable mismatch of fluxes between different-resolution grids at nested interfaces. Although enforced conservation can cause problems on computational accuracy and stability, it is required for long-term integrations. A flux correction algorithm [36,64,65] is applied for volume and tracer conservation. Its basic idea is to integrate the coarse grid variables using the summations of the HRG fluxes at the nested interfaces.

Ideal Salinity Advection Experiment
The HPI scheme and HSIMT AEI scheme described above are evaluated here using an idealized salinity advection test case. The experiment design imitates the test problem used to compare different spatial interpolation schemes in nested grid models in Alapaty, Mathur, and Odman [57].

Experiment Design
Three grids are involved in the ideal experiment. A coarse grid has 60 × 60 cells with a resolution of 30 km, and a fine grid consists of 180 × 180 cells with a resolution of 10 km. The nested grid is produced by refining the central square area of the coarse grid with a refinement ratio of 3 (Figure 4a). The initial salinity field consists of a cylindrical distribution with a radius of four coarse-grid cells near the west boundary of the whole numerical domain. The salinity decreases linearly from 4 psu in the center to 0 psu on the edge (Figure 4b). The uniform eastward 1 m/s flow is specified as the background hydrodynamic condition, and the friction as well as the horizontal/vertical mixing are excluded. Thus, numerical integrations are only determined by the advection terms in the salinity equation. Four simulations which applied the Upwind AEI scheme, QI scheme, HPI scheme, and HSIMT AEI scheme as the interpolation scheme are conducted on the nested grid. Meanwhile, two more numerical cases run on the two structured grids are used to be compared with the nested grid solutions.

Conservation
If numerical schemes are written in flux form, conservation is maintained in a structured grid model. However, this property vanishes for two-way interactive nested models because of the unavoidable mismatch of fluxes between different-resolution grids at nested interfaces. Although enforced conservation can cause problems on computational accuracy and stability, it is required for long-term integrations. A flux correction algorithm [36,64,65] is applied for volume and tracer conservation. Its basic idea is to integrate the coarse grid variables using the summations of the HRG fluxes at the nested interfaces.

Ideal Salinity Advection Experiment
The HPI scheme and HSIMT AEI scheme described above are evaluated here using an idealized salinity advection test case. The experiment design imitates the test problem used to compare different spatial interpolation schemes in nested grid models in Alapaty, Mathur, and Odman [57].

Experiment Design
Three grids are involved in the ideal experiment. A coarse grid has 60 × 60 cells with a resolution of 30 km, and a fine grid consists of 180 × 180 cells with a resolution of 10 km. The nested grid is produced by refining the central square area of the coarse grid with a refinement ratio of 3 (Figure 4a). The initial salinity field consists of a cylindrical distribution with a radius of four coarse-grid cells near the west boundary of the whole numerical domain. The salinity decreases linearly from 4 psu in the center to 0 psu on the edge (Figure 4b). The uniform eastward 1 m/s flow is specified as the background hydrodynamic condition, and the friction as well as the horizontal/vertical mixing are excluded. Thus, numerical integrations are only determined by the advection terms in the salinity equation. Four simulations which applied the Upwind AEI scheme, QI scheme, HPI scheme, and HSIMT AEI scheme as the interpolation scheme are conducted on the nested grid. Meanwhile, two more numerical cases run on the two structured grids are used to be compared with the nested grid solutions.

Results and Discussion
The maximum, the minimum, and the sum of the salinity field are used to evaluate the performance of each interpolation scheme. The maximum indicates the peak preservation skill, and the minimum value shows the oscillations introduced by the interpolation process. Total salt of the whole numerical domain can measure the conservation properties of each scheme. Figure 5 illustrates time series of the three indices for all six cases. The saline water mass reaches the western nesting interface at about hour 94, while completely leaves the refined domain at hour 327 approximately. The maximum value declines due to numerical diffusion. However, the CRG case shows a faster decreasing rate than the HRG one because the numerical diffusion is smaller in the HRG than in the CRG. As seen in Figure 5a, the nested grid solutions exhibit strong ability to retain the salinity peak due to the increased resolution after the scalar distribution is advected to the fine-grid subdomain. The two AEI schemes result in an overshoot of the maximum, which is indictive of poor performance of interpolation scheme, whereas the QI scheme and HPI scheme present better continuity during transition from the coarse-grid subdomain to the fine-grid subdomain.
the performance of each interpolation scheme. The maximum indicates the peak preservation skill, and the minimum value shows the oscillations introduced by the interpolation process. Total salt of the whole numerical domain can measure the conservation properties of each scheme. Figure 5 illustrates time series of the three indices for all six cases. The saline water mass reaches the western nesting interface at about hour 94, while completely leaves the refined domain at hour 327 approximately. The maximum value declines due to numerical diffusion. However, the CRG case shows a faster decreasing rate than the HRG one because the numerical diffusion is smaller in the HRG than in the CRG. As seen in Figure 5a, the nested grid solutions exhibit strong ability to retain the salinity peak due to the increased resolution after the scalar distribution is advected to the finegrid subdomain. The two AEI schemes result in an overshoot of the maximum, which is indictive of poor performance of interpolation scheme, whereas the QI scheme and HPI scheme present better continuity during transition from the coarse-grid subdomain to the fine-grid subdomain.  A significant decrease in minimum occurs several hours before the saline water entering the fine-grid subdomain, and the minimum value recovers afterward (Figure 5b). The magnitude of decrease in the HPI case is the smallest and the minimum quickly resumes to zero at hour 110, which suggests that the HPI scheme is the most stable scheme among the four schemes. However, the low-order upwind AEI scheme delivers the worst performance with the minimum falling to nearly −0.07 psu. The domain minimum behaves similarly for the QI scheme and the HSIMT AEI scheme. Figure 5c presents a comparison of mass conservation properties for each scheme via the temporal variance of relative total salt. The relative total salt undergoes a concave slope because of increased resolution when the saline water mass is entering the fine grid from the coarse grid, and then, it almost remains the same while the salinity distribution is resolved completely in the high-resolution subdomain. A convex variance appears during transition from the fine-grid subdomain to the coarse-grid subdomain. The upwind AEI scheme overpredicted the total salt with a relative value of 100.4%, while on the contrary, the high-order HPI scheme results in an underprediction with a 0.2% deficit. The QI scheme and the HSIMT AEI scheme give more reasonable performances in conserving the total mass compared with the above two schemes.
The differences between the HRG solution and any nested grid solution in the refined subdomain are taken as the absolute errors. Figure 6 illustrates the distribution of absolute errors for each interpolation scheme at hour 212 of simulation (when the saline water mass is at the center of the refined subdomain). For all the interpolation schemes, underpredictions are observed at the center of the distribution, however, overpredictions appear toward both the upwind and downwind edges of the distribution, which can be owed to the numerical diffusion of the advection scheme. The temporal variance of normalized root-mean-square error (NRMSE) is depicted in Figure 7. In fact, four nested-grid simulations show very slight differences between each other in the NRMSE, which is also evidenced by Figure 6.
haves similarly for the QI scheme and the HSIMT AEI scheme. Figure 5c presents a comparison of mass conservation properties for each schem the temporal variance of relative total salt. The relative total salt undergoes a con slope because of increased resolution when the saline water mass is entering the fine from the coarse grid, and then, it almost remains the same while the salinity distribu is resolved completely in the high-resolution subdomain. A convex variance appears ing transition from the fine-grid subdomain to the coarse-grid subdomain. The upw AEI scheme overpredicted the total salt with a relative value of 100.4%, while on the trary, the high-order HPI scheme results in an underprediction with a 0.2% deficit. Th scheme and the HSIMT AEI scheme give more reasonable performances in conservin total mass compared with the above two schemes.
The differences between the HRG solution and any nested grid solution in the ref subdomain are taken as the absolute errors. Figure 6 illustrates the distribution of abs errors for each interpolation scheme at hour 212 of simulation (when the saline water is at the center of the refined subdomain). For all the interpolation schemes, underpr tions are observed at the center of the distribution, however, overpredictions appea ward both the upwind and downwind edges of the distribution, which can be owe the numerical diffusion of the advection scheme. The temporal variance of norma root-mean-square error (NRMSE) is depicted in Figure 7. In fact, four nested-grid sim tions show very slight differences between each other in the NRMSE, which is also denced by Figure 6.

Validation and Assessment
The newborn UFDECOM was applied to a realistic case to test its reliability and robustness in this section. The observational data were used to validate the nested model, and the nested solutions were compared with a reference solution obtained on a structured grid to further evaluate different interpolation schemes, in which the HPI scheme and HSIMT AEI scheme are totally original.

Observational Data
A full description of the measured data can be found in Li et al. [66] and only the necessary information is summarized here. Four synchronous surveys were conducted at three sites (A, B, and C) (Figure 8) in the North Passage of the Changjiang Estuary during four tidal phases, which were the spring tide (16-18 November 2013), the moderate tide after spring tide (MTAS) (20-22 November 2013), the neap tide (24-26 November 2013), and the moderate tide after neap tide (MTAN) (29 November-1 December 2013), respectively. The current and salinity data were adopted for validation.

Validation and Assessment
The newborn UFDECOM was applied to a realistic case to test its reliability and robustness in this section. The observational data were used to validate the nested model, and the nested solutions were compared with a reference solution obtained on a structured grid to further evaluate different interpolation schemes, in which the HPI scheme and HSIMT AEI scheme are totally original.

Observational Data
A full description of the measured data can be found in Li et al. [66] and only the necessary information is summarized here. Four synchronous surveys were conducted at three sites (A, B, and C) (

Validation and Assessment
The newborn UFDECOM was applied to a realistic case to test its reliability and robustness in this section. The observational data were used to validate the nested model, and the nested solutions were compared with a reference solution obtained on a structured grid to further evaluate different interpolation schemes, in which the HPI scheme and HSIMT AEI scheme are totally original.

Observational Data
A full description of the measured data can be found in Li et al. [66] and only the necessary information is summarized here. Four synchronous surveys were conducted at three sites (A, B, and C) (Figure 8) in the North Passage of the Changjiang Estuary during four tidal phases, which were the spring tide (16-18 November 2013), the moderate tide after spring tide (MTAS) (20-22 November 2013), the neap tide (24-26 November 2013), and the moderate tide after neap tide (MTAN) (29 November-1 December 2013), respectively. The current and salinity data were adopted for validation.

Model Configuration & Experiments
The structured grid covers the Changjiang Estuary, and its adjacent seas (Figure 9a), with the minimal cell size reaching nearly 800 × 300 m 2 in the Deep Waterway Project area (Figure 9b). A part of the mesh covering this area was refined with the ratio of 3 (Figure 9c), and the cells were renumbered to make a new grid file for the calculation of UFDECOM. The models have 10 uniform σ layers in the vertical direction.

Model Configuration & Experiments
The structured grid covers the Changjiang Estuary, and its adjacent seas (Figure 9a), with the minimal cell size reaching nearly 800 × 300 m 2 in the Deep Waterway Project area (Figure 9b). A part of the mesh covering this area was refined with the ratio of 3 ( Figure  9c), and the cells were renumbered to make a new grid file for the calculation of UF-DECOM. The models have 10 uniform layers in the vertical direction.  The tidal constants of 16 astronomical constituents (M 2 , S 2 , N 2 , K 2 , K 1 , O 1 , P 1 , Q 1 , MU 2 , NU 2 , T 2 , L 2 , 2N 2 , J 1 , M 1 , and OO 1 ), which are derived from the NaoTide dataset (https://www.miz.nao.ac.jp/staffs/nao99/, accessed on 25 September 2019), were utilized to calculate the tidal level as the open boundary elevation. The open boundary condition for salinity was obtained from the monthly mean salinity of a larger domain model [67].
The initial salinity distribution was interpolated from the Ocean Atlas in Huanghai Sea and East China Sea (Hydrology) [68] outside the river mouth and from recently collected field data in the estuary. The daily Changjiang River discharge was measured at the Datong Hydrological Station. The 6-hourly sea surface wind data from the European Center for Medium-Range Weather Forecasts (ECMWF) reanalysis data products with the resolution of 0.125 • × 0.125 • were used to force the models.
Five experiments were conducted to validate the nested model and assess the performance of four interpolation schemes. A simulation was run on the structured grid in Experiment 0 (Exp. 0) to produce a reference solution, while Experiment 1-4 (Exp. [1][2][3][4] were conducted on the nested grid with the FWO scheme as update scheme and four different interpolation schemes, respectively ( Table 1). The flux correction method was applied to the nested grid model. All the simulations started from 1 October 2013 and ran for 70 days.

Validation and Interpolation Scheme Assessment
Performance of the models was quantified by the correlation coefficient (CC), the root-mean-square error (RMSE), and the skill score (SS) [69]. The formulae are given as where X mod and X obs are the variables of model results and observational data, respectively, and X denotes the time average of the variable. The SS was introduced by [70] as a statistical metric to describe the degree to which the observed deviations from the observed mean correspond to the simulated derivations from the observed mean. Its value ranges from 0 to 1, representing total disagreement to perfect agreement between model results and observational data. The validation statistics of flow velocity are summarized in Table 2. In Exp. 0, the value of CC at the three sites ranged from 0.79 to 0.85 on the surface layer and from 0.76 to 0.88 on the bottom layer; the RMSE ranged from 0.35 to 0.41 on the surface layer and from 0.29 to 0.32 on the bottom layer. The reference experiment showed excellent performances at both sites A and B with the same SS values (0.91/0.84), while the SS at site C yielded slightly lower values, which were 0.88 on the surface layer and 0.73 on the bottom layer. Among the four experiments conducted on a nested grid, Exp. 3 was a notch above the other three experiments with the HSIMT parabolic interpolation scheme especially at site C, where the simulation obtained an SS of 0.72 on the bottom layer, which was the closest to that of Exp. 0. The performances of the other three experiments differed very little from each other. Thus, a deeper comparison was made among the observed data, the reference solution, and the results of Exp. 3 (Figures 10 and 11). Overall, the temporal variance of the surface and bottom velocity in Exp. 3 agreed well with the results of Exp. 0 and the measured data except for that Exp. 3 presented higher peaks of speed on the surface layer at some moments, e.g., at the 4th hour of the MTAN (site B) (Figure 11f), at the 18th hour of the spring tide (site C) (Figure 10i).      Table 3 shows the validation statistics of salinity. The CC of Exp. 0 was between 0.73 and 0.82 on the surface layer, while between 0.77 and 0.82 on the bottom layer. The surface RMSE ranged from 2.44 to 3.21, whereas the bottom layer obtained more significant RMSE ranging from 4.94 to 6.03. Since the minimum SS value was 0.77 on the surface layer (site A) and 0.76 on the bottom layer (site C), we can make a judgment that the reference experiment simulated the salinity transport process quite successfully. However, the performances of the nested models degraded greatly compared with that of Exp. 0. In Exp. 1, 2, and 4, the SS values are about 0.3 lower than that of the reference solution on the surface layer, the gap even increased to 0.45 on the bottom layer of site C (between Exp. 1 and Exp. 0). Even so, Exp. 3 still stood out from the four nested grid simulations. The surface and bottom salinity at three sites in Exp. 3 showed quite consistent temporal variances with the results of Exp. 0 and observational data during the spring tide and MTAS ( Figure  12a,b,e,f,i,j). The nested grid solution was absolutely higher than the reference solution on  Table 3 shows the validation statistics of salinity. The CC of Exp. 0 was between 0.73 and 0.82 on the surface layer, while between 0.77 and 0.82 on the bottom layer. The surface RMSE ranged from 2.44 to 3.21, whereas the bottom layer obtained more significant RMSE ranging from 4.94 to 6.03. Since the minimum SS value was 0.77 on the surface layer (site A) and 0.76 on the bottom layer (site C), we can make a judgment that the reference experiment simulated the salinity transport process quite successfully. However, the performances of the nested models degraded greatly compared with that of Exp. 0. In Exp. 1, 2, and 4, the SS values are about 0.3 lower than that of the reference solution on the surface layer, the gap even increased to 0.45 on the bottom layer of site C (between Exp. 1 and Exp. 0). Even so, Exp. 3 still stood out from the four nested grid simulations.
The surface and bottom salinity at three sites in Exp. 3 showed quite consistent temporal variances with the results of Exp. 0 and observational data during the spring tide and MTAS (Figure 12a,b,e,f,i,j). The nested grid solution was absolutely higher than the reference solution on both layers during the neap tide (Figure 12c,g,k), while the salinity variance in Exp. 3 gradually fitted in with that in Exp. 0 during the MTAN.
both layers during the neap tide (Figure 12c,g,k), while the salinity variance in Exp. 3 gradually fitted in with that in Exp. 0 during the MTAN.    In summary, the five experiments had little difference in reproducing the hydrodynamic process at the Deep Waterway Project area. As for simulating the temporal variance of salinity, Exp. 3 provided the best performance among the nested grid simulations, which suggested that the HSIMT parabolic interpolation scheme could be the best interpolation scheme for UFDECOM in the test case. However, the performances of nested grid models were not even comparable to that of the reference model.

Discussion
The unignorable deterioration of the nested grid model compared to the structured grid model in simulating the salinity transport process undermines the robustness of UFDECOM to some extent. Considering that the flux correction method adjusts the flux of the CRG at the nesting interface artificially, this conservation constraint is suspected to cause the accuracy loss. An extra experiment (Exp. 5) was conducted to investigate the influence of the flux correction method on salinity simulation. The experiment settings were inherited from Exp. 3 except that the conservation restriction was removed.
The validation statistics indicated that Exp. 5 achieved significant improvement compared with Exp. 3 ( Table 4). The modeled salinity yielded an average CC of 0.77 on the surface layer and 0.78 on the bottom layer. The RMSE values reduced to a normal level and the surface/bottom maxima both occurred at site C (3.59/5.11). Most importantly, the SS values had notable increases especially at sites B and C, and the bottom SS values were even a little higher than that of Exp. 0.  Figure 13 illustrates the salinity variance of measured data, Exp. 0 and Exp. 5 at three observational sites during the survey periods. Compared with salinity variance in Exp. 3 (Figure 12), the surface salinity of the new experiment increased, while the bottom salinity decreased, which indicated that the nested grid solution was much more consistent with the reference solution. The blue dashed line went under the solid line in Figure 13c,d,g,h,k,l, which meant that the results of Exp. 5 were closer to the observed salinity at the bottom layer during the neap tide and MTAN. Basically, the model reproduced the salinity variances at three sites much better without enforced conservation.

Conclusions
A three-dimensional, two-way nesting unstructured quadrilateral grid, primitive equations, finite-differencing, estuarine and coastal ocean model was described. The model grid is capable of multi-area nesting and multi-level refinement, and a virtual-grid method was derived to implement two-way interaction between coarse and fine grids. Virtual grids, of which independent variables are obtained through updating and interpolating from actual grid variables, are created next to nesting interfaces for the integration of adjacent actual grid variables. Plenty of update and interpolation schemes are available in UFDECOM, and among the interpolation schemes, the HSIMT parabolic in-

Conclusions
A three-dimensional, two-way nesting unstructured quadrilateral grid, primitive equations, finite-differencing, estuarine and coastal ocean model was described. The model grid is capable of multi-area nesting and multi-level refinement, and a virtual-grid method was derived to implement two-way interaction between coarse and fine grids. Virtual grids, of which independent variables are obtained through updating and interpolating from actual grid variables, are created next to nesting interfaces for the integration of adjacent actual grid variables. Plenty of update and interpolation schemes are available in UFDECOM, and among the interpolation schemes, the HSIMT parabolic interpolation scheme and HSIMT AEI scheme are newly developed with high-order accuracy and great compatibility with the HSIMT-TVD advection scheme. The model adopts a flux correction algorithm for enforced volume and tracer conservation.
The two novel schemes are tested in an ideal salinity advection experiment and they are compared with the upwind AEI scheme and the QI scheme in the peak preservation skill, stability, and the mass conservation characteristic. The QI scheme and the HPI scheme show minimal discontinuity of the maximum during transition, while the peak was overpredicted in the two AEI scheme cases. The upwind AEI scheme generates the largest oscillation when the saline water is reaching the nesting interface, whereas the HPI scheme provides the best performance in stability. The best conservation properties are observed in the cases with the QI scheme and the HSIMT AEI scheme. All schemes show similar field distribution and temporal variance of the normalized RMSE.
The model was applied to a realistic simulation for validation and interpolation scheme assessment. The nested grid model basically remained at the same level as the structured grid model in simulating current, while it fell behind when modeling the tracer transport process, during which the experiment with the HSIMT parabolic interpolation scheme delivered the best performance among the four experiments conducted on a nested grid. A comparison between the results of Exp. 3 and that of Exp. 5 with no enforced conservation revealed that the flux correction method accounted for the degradation of salinity simulation to a large extent, which also warned us that great care should be taken in applying this technique despite its necessity in long-term integrations.

Data Availability Statement:
No new data were created or analyzed in this study. Data sharing is not applicable to this article.