A Wet/Dry Point Treatment Method of FVCOM, Part I: Stability Experiments

: A 3-dimensional wet/dry point treatment method was developed for the unstructured-grid Finite-Volume Community Ocean Model (FVCOM). Analytical equations were derived to examine discretized errors that occurred during the ﬂooding/drying process by the wet/dry point treatment for the ﬂooding/drying process. Numerical experiments were carried out for an idealized estuary, including the inter-tidal zone. The model results show that if the ratio of internal to external mode time steps ( I split ) is appropriately selected, FVCOM was capable of simulating the ﬂooding/drying process with sufﬁcient accuracy to ensure the mass conservation. The up-bound limit of I split was restricted by the bathymetric slope of the inter-tidal zone, external mode time step, horizontal/vertical resolution, and amplitude of tidal forcing at the open boundary, as well as the thickness of the viscous layer speciﬁed in the model. Criteria for time steps via these parameters were derived from these experiments, which provide a helpful guide in selecting I split for applying FVCOM to realistic geometric estuaries.


Introduction
One of the most challenging problems in estuarine modeling is to provide an accurate simulation of the water transport flooding into and draining out of the intertidal zone. In many estuaries, this transport is comparable to or larger than the transport carried in the main river channel, so its effect on the overall tidal dynamics and salt or nutrient fluxes can be substantial. In the past, two types of approaches have been developed to solve this problem. One is the moving boundary method, and the other is the wet/dry point treatment method. In the first method, the computational domain is bounded by an interface line between land and water, where the total depth and transport normal to the boundary are equal to zero [1]. Because the boundary moves over flooding and drying cycles, the model grid must be re-generated at every time step. This method works for idealized estuaries with simple geometry [2,3], but it is not practical for realistic estuaries with complex geometries, including tidal creeks, islands, barriers, and inlets. In the second method, the computational domain covers the maximum flooding area, with each grid point determined to be either wet or dry as the flow field evolves in time. The wet and dry points are distinguished by the local total water depth D(x, y) = H(x, y) + ζ(x, y, t) (H is the reference water depth, and ζ is the surface elevation). A grid point with D > 0 is wet. Otherwise, it is dry. At a dry point, water velocities are automatically set to zero, but the salinity remains the same from the previous time step. Since this method is relatively simple, it has been widely used to simulate the water transport over inter-tidal areas in estuarine models [4][5][6][7][8].
The wet/dry point treatment technique works only in the case where a finite-value solution of the governing equations exists as D approaches zero. In a z-coordinate system, to ensure numerical stability of a 3-dimensional (3D) model, the thickness of the layer closest to the surface must be greater than the amplitude of the tidal elevation [9]. This makes it challenging to apply this method to a shallow estuary in which the amplitude of the tidal elevation is the same order of magnitude as the local water depth. Because irregular bathymetry is difficult to resolve in a z-coordinate model, it is hard to believe that a z-coordinate model can accurately simulate water transport in an intertidal area with realistic, complex bathymetry during both flooding and drying periods [10][11][12][13][14]. In a terrain-following coordinate system, the wet/dry point treatment is no longer valid as D becomes zero. A straightforward way to avoid numerical singularity is to specify zero velocity at all dry points. This approach, however, does not ensure volume conservation in the numerical solution due to discontinuous water removal from elements that turn from wet to dry over one time step. An alternative way is to add a viscous boundary layer (D min ) at the bottom and redefine wet/dry points using a sum of D and D min . The grid is treated as a wet point for D > D min . Otherwise, it is a dry point. In terms of the vertical structure of turbulent mixing, a viscous sublayer always exists below the log boundary layer near a solid wall [15]. However, to avoid adding additional water transport into a dynamic system, the viscous sublayer should be sufficiently thin to satisfy a motionless condition. Good examples of the application of this method can be found in [7,8].
No matter which methods are used to simulate the flooding/drying process over the intertidal zone in an estuary, they should be validated with respect to volume conservation. In our previous experiences in implementing a wet/dry point treatment to the semiimplicit Estuarine-Coastal Ocean Model (ECOM-si) [16], a semi-implicit version of the Princeton Ocean Model (POM), we found it is difficult to ensure the salt conservation over the intertidal area [8]. Therefore, checking the volume conservation should be the first step to validating the dry/wet point treatment method. The easiest way is to specify the temperature and salinity constant everywhere as the initial condition and run the model to see if the model can maintain these constant values over the time integration. In any of these methods, the dry and wet points are determined using some empirical criterion so that the computed water transport in the dry-wet transition zone depends on (1) the criterion to define the wet/dry points, (2) the time step used for numerical integration, (3) the horizontal and vertical grid resolutions, (4) the surface elevation amplitude, and (5) the bathymetry. In a terrain-following coordinate system, the transport may also depend on the thickness of the viscous bottom sublayer (D min ). To our knowledge, these issues have not been discussed in detail in previous modeling studies of estuarine intertidal flooding and draining. This paper describes a wet/dry point treatment method introduced into the unstructuredgrid, 3-dimensional (3D) free-surface primitive equation Finite-Volume Community Ocean Model (FVCOM). FVCOM was initially developed by [17] and has been upgraded by the FV-COM development team at the University of Massachusetts and the user community [18][19][20][21][22][23]. We have successfully applied this method to simulate the flooding/drying process in different estuarine systems [24][25][26] and storm-and tsunami-induced coastal inundations [27][28][29][30][31] in the last decades, but the wet/dry treatment methodology implemented in FVCOM was never published. This paper is focused on the stability analysis for an idealized semi-enclosed estuary with an intertidal zone. The relationships of the time step with grid resolution, the amplitude of external forcing, the slope of the inter-tidal area, and the thickness of the viscous bottom layer are discussed, and the criterion for the selection of the time step is derived. The rule used in the model validation is mass conservation, which, we believe, is a prerequisite condition for an objective evaluation of the wet/dry point treatment technique in estuarine models.

The FVCOM
FVCOM encompasses seven governing equations: three for momentum, one for incompressible continuity, two for temperature and salinity, and one for density. It is closed mathematically using the Generalized Ocean Turbulent Model (GOTM) for vertical mixing [32] and a Smagorinsky eddy parameterization for horizontal diffusion [33]. The model uses the generalized terrain-following transformation in the vertical to provide a smooth and accurate representation of irregular bottom topography [34].
FVCOM subdivides the horizontal numerical computational domain into a set of non-overlapping unstructured triangular cells. Each triangle is comprised of three nodes, a centroid, and three sides (Figure 1), on which the x and y components of the horizontal velocity (u and v) are placed at centroids, and all scalar variables, such as ζ (water elevation), H (reference water depth), D (the total water depth), ω (vertical velocity at the σ-surface), S (salinity), T (temperature), ρ (density), K m (vertical eddy viscosity), K h (vertical thermal diffusion coefficient), A m (horizontal eddy diffusion coefficient), and A h (horizontal thermal diffusion coefficient) are placed at nodes. Scalar variables at each node are determined by a net flux through the sections linked to centroids and the middle point of the edges in the surrounding triangles (called the tracer control element: TCE), while u and v at centroids are calculated based on the net flux through three sides of that triangle (called the momentum control element: MCE).

The FVCOM
FVCOM encompasses seven governing equations: three for momentum, one for incompressible continuity, two for temperature and salinity, and one for density. It is closed mathematically using the Generalized Ocean Turbulent Model (GOTM) for vertical mixing [32] and a Smagorinsky eddy parameterization for horizontal diffusion [33]. The model uses the generalized terrain-following transformation in the vertical to provide a smooth and accurate representation of irregular bottom topography [34].
FVCOM subdivides the horizontal numerical computational domain into a set of non-overlapping unstructured triangular cells. Each triangle is comprised of three nodes, a centroid, and three sides (Figure 1), on which the x and y components of the horizontal velocity (u and v) are placed at centroids, and all scalar variables, such as ζ (water elevation), H (reference water depth), D (the total water depth), ω (vertical velocity at the σsurface), S (salinity), T (temperature), ρ (density), Km (vertical eddy viscosity), Kh (vertical thermal diffusion coefficient), Am (horizontal eddy diffusion coefficient), and Ah (horizontal thermal diffusion coefficient) are placed at nodes. Scalar variables at each node are determined by a net flux through the sections linked to centroids and the middle point of the edges in the surrounding triangles (called the tracer control element: TCE), while u and v at centroids are calculated based on the net flux through three sides of that triangle (called the momentum control element: MCE). FVCOM is solved using the mode-split (external and internal modes) scheme with the updating of the semi-implicit scheme as an alternative choice [21]. The discussion in this paper is aimed at the mode-split method with a σ-coordinate transformation: an opensource version that has been successfully used to simulate the flooding/drying process in various estuaries (http://fvcom.smast. umassd.edu/medm-publications/ (accessed on 30 May 2022)). In this method, the external mode is composed of the vertically integrated transport form of the governing equations. The transport equations are solved for the sea surface elevation by an explicit difference scheme in which the time step ∆ is constrained by min ( , , , where , , and are the three side lengths of the smallest size triangle, g is gravity, D is the total local water depth, and is the local shallow water wave speed. The internal mode is comprised of the fully 3-dimensional governing FVCOM is solved using the mode-split (external and internal modes) scheme with the updating of the semi-implicit scheme as an alternative choice [21]. The discussion in this paper is aimed at the mode-split method with a σ-coordinate transformation: an opensource version that has been successfully used to simulate the flooding/drying process in various estuaries (http://fvcom.smast.umassd.edu/medm-publications/ (accessed on 30 May 2022)). In this method, the external mode is composed of the vertically integrated transport form of the governing equations. The transport equations are solved for the sea surface elevation ζ by an explicit difference scheme in which the time step ∆t E is constrained by min (l 1 , l 2 , l 3 ) gD, where l 1 , l 2 , and l 3 are the three side lengths of the smallest size triangle, g is gravity, D is the total local water depth, and gD is the local shallow water wave speed. The internal mode is comprised of the fully 3-dimensional governing equations and solved numerically using a longer time step ∆t I . The internal mode time step is constrained by the phase speed of the lowest interval wave mode. The linkage between external and internal modes is through ζ. The surface pressure gradients at each ∆t I time step in the internal mode are calculated directly by the output of ζ from the external mode. A second-order accuracy upwind finite-difference scheme is used for flux calculation in the integral form of the advective terms [35,36], and the modified fourth-order Runge-Kutta time-stepping scheme is used for time integration.
In a mode-split model, an adjustment must be made in every internal time step to ensure consistency in the vertically integrated water transport produced by the external and internal modes. In FVCOM, the vertical velocity at the σ-surface (ω) is calculated by where i, k, and n are the indices of the ith node point, kth σ level and nth time step, respectively, and Ω n i is the area of the ith TCE at the nth time step. v n N,k indicates the velocity component normal to the boundary of a TCE with a length l. To conserve the volume in the ith TCE, the vertically-integrated form of Equation (1) must satisfy where kb is the total number of the σ-levels. Since ζ n+1 i is calculated through the vertically integrated continuity equation over I split external mode time steps (where Because the numerical accuracy depends on the time step, the left and right sides of (3) are not precisely equal unless I split = 1. Therefore, to ensure the conservation of the volume transport throughout the water column of the ith TCE, the internal velocity in each σ layer must be calibrated using the difference of inequality of vertically integrated external and internal water transport before ω is calculated. Since the vertically-integrated transport features the barotropic motion, the most straightforward calibration is to adjust the internal velocity by distributing the difference or "error" transport uniformly throughout the water column. ω, which is calculated with adjusted internal velocities through (1), thus guarantees that volume transport is conserved in both the whole water column and the individual TCE volume in each σ layer. Because the water temperature and salinity or other tracers are calculated in the same TCE volume as that used for ω, this transport adjustment also conserves mass in an individual TEC volume in each σ layer.

An Unstructured-Grid Wet/Dry Point Treatment Method
The following wet/dry point treatment method has been developed for FVCOM. Similar to our previous practice with ECOM-si [8], we introduce a viscous sublayer with a thickness of D min into FVCOM to avoid the occurrence of singularity when the local water depth approaches zero. The wet or dry criterion for node points is and for triangular cells is where h B is the bathymetric height related to the edge of a river where the mean water depth is zero ( Figure 2) andî,ĵ, andk are the integer numbers to identify the three-node points of a triangular cell. When a triangular cell is treated as dry, the velocity at the centroid of this triangle is specified to be zero, and no flux is allowed through the three side boundaries of this triangle. This triangular cell is removed from the flux calculation in the TCE. For example, the integral form of the continuity equation in FVCOM is written as where u and v are the x and y components of the vertically-averaged velocity. In a dry/wet point system, only wet triangles are taken into account in the flux calculation in a TCE since the flux on the boundaries of the dry triangle is zero (see Figure 3). This approach always ensures volume conservation in a TCE that contains the moving boundary between the dry and wet triangles. The same procedure is used to calculate the tracer flux (temperature, salinity, and other scalar tracers) and momentum flux in an MCE. and for triangular cells is where ℎ is the bathymetric height related to the edge of a river where the mean water depth is zero ( Figure 2) and , , and are the integer numbers to identify the three-node points of a triangular cell. When a triangular cell is treated as dry, the velocity at the centroid of this triangle is specified to be zero, and no flux is allowed through the three side boundaries of this triangle. This triangular cell is removed from the flux calculation in the TCE. For example, the integral form of the continuity equation in FVCOM is written as where and are the x and y components of the vertically-averaged velocity. In a dry/wet point system, only wet triangles are taken into account in the flux calculation in a TCE since the flux on the boundaries of the dry triangle is zero (see Figure 3). This approach always ensures volume conservation in a TCE that contains the moving boundary between the dry and wet triangles. The same procedure is used to calculate the tracer flux (temperature, salinity, and other scalar tracers) and momentum flux in an MCE.  One critical issue in applying the wet/dry point treatment technique in a split-mode model is to ensure mass conservation of individual TCEs which change "state (wet or dry)" one or more times during the I split external mode integrations made during each internal mode time step; Ωn i may change within I split external time integrations due to the occurrence of dry triangles and is treated as zero when D is less than D min , In this case, the external and internal mode adjustment through (3) cannot guarantee that ω reaches zero at σ = −1 for the internal mode. To ensure volume conservation, an additional adjustment for ζ n+1 i must be made in the TCE when ω is calculated by (1). The additional sea-level adjustment generally works, but fails in the case where ζ is very close to D min (for example, ∆ζ = ζ − D min < 10 −1 (cm)). When this happens, minor errors in calculating the volume flux can rapidly accumulate through nonlinear feedback of tracer advection terms and eventually lead to significant errors in mass conservation. One critical issue in applying the wet/dry point treatment technique in a split-mode model is to ensure mass conservation of individual TCEs which change "state (wet or dry)" one or more times during the Isplit external mode integrations made during each internal mode time step; Ω may change within external time integrations due to the occurrence of dry triangles and is treated as zero when D is less than , In this case, the external and internal mode adjustment through (3) cannot guarantee that ω reaches zero at = −1 for the internal mode. To ensure volume conservation, an additional adjustment for must be made in the TCE when ω is calculated by (1). The additional sea-level adjustment generally works, but fails in the case where is very close to Dmin (for example, ∆ = − 10 (cm)). When this happens, minor errors in calculating the volume flux can rapidly accumulate through nonlinear feedback of tracer advection terms and eventually lead to significant errors in mass conservation. When = 1 (external and internal time steps are identical), FVCOM guarantees mass conservation in the wet-dry transition zone. When computational efficiency requires that > 1, several criteria are required to ensure mass conservation in FVCOM. In addition to the criterion for numerical stability, the choice of depends on flooding/drying situations: the surface elevation, bathymetry, thickness of the viscous bottom layer, and horizontal/vertical resolutions. A discussion on the relationship of with these factors is given next through numerical experiments for several idealized cases.

Design of Numerical Experiments
The numerical experiments were conducted for an idealized semi-enclosed channel with a width of 3 km at the bottom, a length of 30 km, a constant depth of 10 m, and a lateral slope of about 0.033 ( Figure 4). This channel is oriented east to west, connected to a relatively wide and flat-bottom shelf to the east and intertidal zones along the northern and southern edges. The inter-tidal zone is distributed symmetrically to the channel, with a constant slope of α and a width of 2 km. Our focus was to examine numerical instability due to the nonconservative volume issue resulting from the 2D and 3D mode adjustment in the wet-dry transition zone under the case of 1. As aforementioned, the finite When I split = 1 (external and internal time steps are identical), FVCOM guarantees mass conservation in the wet-dry transition zone. When computational efficiency requires that I split > 1, several criteria are required to ensure mass conservation in FVCOM. In addition to the criterion for numerical stability, the choice of I split depends on flooding/drying situations: the surface elevation, bathymetry, thickness of the viscous bottom layer, and horizontal/vertical resolutions. A discussion on the relationship of I split with these factors is given next through numerical experiments for several idealized cases.

Design of Numerical Experiments
The numerical experiments were conducted for an idealized semi-enclosed channel with a width of 3 km at the bottom, a length of 30 km, a constant depth of 10 m, and a lateral slope of about 0.033 ( Figure 4). This channel is oriented east to west, connected to a relatively wide and flat-bottom shelf to the east and intertidal zones along the northern and southern edges. The inter-tidal zone is distributed symmetrically to the channel, with a constant slope of α and a width of 2 km. Our focus was to examine numerical instability due to the nonconservative volume issue resulting from the 2D and 3D mode adjustment in the wet-dry transition zone under the case of I split = 1. As aforementioned, the finite volume discrete algorithms used in FVCOM ensure volume conservation. However, this conservation law was not valid in the wet-dry transition zone. Choosing a simple idealized geometric problem helped us distinguish the I split -related issues from other types of nonlinear instability due to irregular bathymetry changes and complex current interaction over the estuarine-tidal creek-wetland complex. The results obtained from this idealized benchmark testing problem can provide a guide to select I split in real estuarine applications. It was demonstrated in the FVCOM simulation of the Satilla River Estuary, Georgia, US [24].
The computational domain was configured with unstructured triangular grids in the horizontal and σ-levels in the vertical. Numerical experiments were conducted for cases with different horizontal and vertical resolutions. The comparison between these cases was made based on differences from a standard run with a horizontal resolution of about 500 m in the channel, 600 to 1000 m over the shelf, and 600 to 900 m over the inter-tidal zone (Figure 4). In the rest of the text, the standard run refers to the case with the horizontal grid shown in Figure 4 and the following parameters: ∆t E = 4.14 s, I split = 6, D min = 5 cm, kb (the number of σ-levels) = 6, and α =7 × 10 −4 (a change of the height of the inter-tidal zone up to 0.8 m over a distance of 2 km in the cross-channel direction). ∆t E was chosen to satisfy the linear instability criterion constrained by min(l 1 , l 2 , l 3 )/ gD. In this case, the minimum value of the horizontal resolution over the inter-tidal zone (∆L) is 600 m.
was made based on differences from a standard run with a horizontal resolution of about 500 m in the channel, 600 to 1000 m over the shelf, and 600 to 900 m over the inter-tidal zone (Figure 4). In the rest of the text, the standard run refers to the case with the horizontal grid shown in Figure 4 and the following parameters: ∆ = 4.14 s, Isplit = 6, Dmin = 5 cm, kb (the number of σ-levels) = 6, and α =7 × 10 −4 (a change of the height of the inter-tidal zone up to 0.8 m over a distance of 2 km in the cross-channel direction). ∆ was chosen to satisfy the linear instability criterion constrained by ( , , ) / . In this case, the minimum value of the horizontal resolution over the inter-tidal zone (∆ ) is 600 m. Specifying an initial distribution of water temperature T = 20 • C and water salinity S = 35 PSU everywhere, we ran the model prognostically with T and S calculated at every internal mode time step. This numerical exercise was done under a pre-condition of numerical stability for the cases with different parameters of ∆t E , I split , ζ o , D min , α, ∆L, and ∆σ. In each case, the model was run continuously for 10 model days. According to the law of mass conservation, T and S should remain unchanged over the entire numerical integration for all these cases since there are no external heat and salt sources or sinks in the model problem. Therefore, any change in T and S in the dry-wet transition zone must be due to mass conservation errors during the integration. Since such errors depend primarily on I split , the upper bound of I split for given geometric and physical forcing conditions can be defined as the value for which changes in T and S remain small (∆T > 10 −3 • C or ∆S > 10 −3 PSU) over the entire (long-term) model integration. Based on this criterion, we can use this idealized case exercise to derive a quantitative relationship of the upper-bound value of I split with ζ o , D min , α, ∆L, and ∆σ.

The Standard Case Flooding/Drying Process
The numerical experiments conducted in this study clearly show that the wet/dry point technique introduced into FVCOM is capable of simulating tide-induced flooding and drying. An example of the model-predicted variation of the water elevation over an M 2 tidal cycle is shown in Figure 5 for the standard case with ζ o = 1.5 m. During flood tide, the water elevation gradually rises as the inflow of the tidal current increases ( Figure 5: the upper-left panel). The water floods onto the intertidal zone when the water level is over the local height ( Figure 5: the upper-right panel). The entire inter-tidal zone is covered by water at the transition time from flood to ebb ( Figure 5: the lower-left panel0, and then the water drains back into the channel during the ebb tide ( Figure 5: the lower-right panel). This flooding/drying process repeats every M 2 tidal cycle, and the model runs stably with ∆T or ∆S significantly smaller than 10 −3 • C or PSU for all time.  Selecting sites A (at the channel center), B (at the edge of the main channel and intertidal zone), and C (within the intertidal zone), we examined the changes in water elevation ( ), along-channel current (U), and cross-channel current (V) over tidal cycles ( Figure 6). The water elevation at sites A and B were almost identical. They changed stably with an M2 tidal period after a two-tidal cycle spin-up. At site A, the along-channel current was predominant, changing nearly symmetrically over the flood and ebb-tidal periods. At site B, although the water elevation varied symmetrically over the flood and ebb-tidal periods, the currents over the ebb-tidal period varied slowly and lasted much longer. Site C became wet during the flood-tidal period. At that site, the water elevation reached its maximum synchronously with the water change within the channel. The water was drained rapidly during the ebb-tidal period. The change in the cross-channel velocity varied with M2 tidal Selecting sites A (at the channel center), B (at the edge of the main channel and intertidal zone), and C (within the intertidal zone), we examined the changes in water elevation (ζ), along-channel current (U), and cross-channel current (V) over tidal cycles ( Figure 6). The water elevation at sites A and B were almost identical. They changed stably with an M2 tidal period after a two-tidal cycle spin-up. At site A, the along-channel current was predominant, changing nearly symmetrically over the flood and ebb-tidal periods. At site B, although the water elevation varied symmetrically over the flood and ebb-tidal periods, the currents over the ebb-tidal period varied slowly and lasted much longer. Site C became wet during the flood-tidal period. At that site, the water elevation reached its maximum synchronously with the water change within the channel. The water was drained rapidly during the ebb-tidal period. The change in the cross-channel velocity varied with M 2 tidal cycles due to temporospatial variation of vertical viscosity and horizontal diffusion.

Relationship with α and
The upper bound of Isplit varies with the slope of the intertidal zone (α) and amplitude of tidal forcing ( ) (Figure 7). Considering the standard case but with a smaller slope of α = 4.0 × 10 , the upper bound of Isplit gradually becomes smaller as increases. It is This result suggests that if I split is appropriately selected, FVCOM can simulate the flooding/drying process with sufficient accuracy that mass conservation errors are insignificant. Clearly, by this criterion, the upper bound for I split exceeds 6 for the standard case conditions. Next, we consider how the upper bound of I split depends on the other parameter values.

Relationship with α and ζ o
The upper bound of I split varies with the slope of the intertidal zone (α) and amplitude of tidal forcing (ζ o ) (Figure 7). Considering the standard case but with a smaller slope of α The change in the upper bound of Isplit relative to α and is consistent with the physics of the flooding/drying process. The flooding speed to cover the intertidal zone during flood increases as α decreases. For a given time step of the external mode, any inconsistency between the external and internal modes becomes more critical to mass conservation as α is decreased, thus the upper bound of Isplit must become smaller to conserve mass during the flooding and drying cycle. Since the flooding speed also increases with increasing tidal forcing amplitude , the inequality shown in Equation (7) becomes more significant for larger due to a numerical cutoff of ζ when D < Dmin. Thus, the upper bound of Isplit generally becomes smaller as becomes larger. In these experiments, ∆ = 4.14 s, kb = 6, and Dmin = 5 cm.

Relationship with Dmin
In general, under given tidal forcing, vertical/horizontal resolution, and external mode time step, the upper bound of Isplit increases as the thickness of the viscous layer Dmin becomes larger (Figure 8). In the standard case with = 1.5 m and α = 4.0 × 10 −4 , for example, Isplit must be smaller or equal to 9 when Dmin = 5 cm but 22 when Dmin = 20 cm. Isplit could be much larger for a greater bottom slope α of the intertidal zone. In the cases with α = 7.0 × 10 −4 and 9.0 × 10 −4 , the upper bound of Isplit could be up to 10 and 13, respectively, for Dmin = 5 cm and up to 28 and 33, respectively, for Dmin = 20 cm.
The relationship of Isplit with Dmin is not a linear function. The upper bound of Isplit tends to converge towards a constant value for Dmin > 20 cm in the three cases with In these experiments, ∆t E = 4.14 s, kb = 6, and D min = 5 cm.
The change in the upper bound of I split relative to α and ζ o is consistent with the physics of the flooding/drying process. The flooding speed to cover the intertidal zone during flood increases as α decreases. For a given time step of the external mode, any inconsistency between the external and internal modes becomes more critical to mass conservation as α is decreased, thus the upper bound of I split must become smaller to conserve mass during the flooding and drying cycle. Since the flooding speed also increases with increasing tidal forcing amplitude ζ o , the inequality shown in Equation (7) becomes more significant for larger ζ o due to a numerical cutoff of ζ when D < D min . Thus, the upper bound of I split generally becomes smaller as ζ o becomes larger.

Relationship with D min
In general, under given tidal forcing, vertical/horizontal resolution, and external mode time step, the upper bound of I split increases as the thickness of the viscous layer D min becomes larger (Figure 8). In the standard case with ζ o = 1.5 m and α = 4.0 × 10 −4 , for example, I split must be smaller or equal to 9 when D min = 5 cm but 22 when D min = 20 cm. I split could be much larger for a greater bottom slope α of the intertidal zone. In the cases with α = 7.0 × 10 −4 and 9.0 × 10 −4 , the upper bound of I split could be up to 10 and 13, respectively, for D min = 5 cm and up to 28 and 33, respectively, for D min = 20 cm. J. Mar. Sci. Eng. 2022, 10, x FOR PEER REVIEW 12 of 18 different slopes. Increasing Dmin directly reduces the occurrence of a negative D in a TCE at the wet-dry edge and thus enhances the numerical consistency between external and internal modes. However, a blind increase of Dmin would destroy the physics of flooding/drying by adding too much additional water transport in the wet area, particularly in the shallower inter-tidal zone where the maximum water level might be only 1 to 2 m or less. A value of Dmin = 5 cm was found to work well in recent FVCOM simulations of several tidal estuaries in Georgia and South Carolina [24,25], where the maximum flushing area is defined approximately by the 2-m elevation line. The ratio of Dmin/max(hB) equaled 40, sufficiently small to ensure adequate mass conservation in these realistic model applications.

Relationship with kb
The upper bound of Isplit with respect to vertical resolution is sensitive to Dmin ( Figure  9). For the standard case with kb = 6, α= 7.0 × 10 −4 , and = 1.5, for example, the upper bound of Isplit is 10 for Dmin = 5 cm and about 28 for Dmin = 20 cm. Keeping the same forcing and geometry but increasing the number of sigma levels kb to 11, the upper bound of Isplit remains almost the same for Dmin = 5 cm but drops as Dmin increases.
In the case with α = 7.0 × 10 −4 , defining H = 0 at the edge of the channel, the height of the land at the boundary of the inter-tidal zone is 1.4 m. Considering an amplitude of = 1.5 m, when the intertidal area is entirely covered by the water during the flood tide, the water elevation at the edge of the intertidal zone should be 10 cm. In the actual numerical computation, however, the water elevation is adjusted by an additional Dmin in the region after becoming wet. For a case with Dmin = 10 cm, for example, adding Dmin is equivalent to double the water elevation at the outer edge boundary of the intertidal zone. In this case, The relationship of I split with D min is not a linear function. The upper bound of I split tends to converge towards a constant value for D min > 20 cm in the three cases with different slopes. Increasing D min directly reduces the occurrence of a negative D in a TCE at the wetdry edge and thus enhances the numerical consistency between external and internal modes. However, a blind increase of D min would destroy the physics of flooding/drying by adding too much additional water transport in the wet area, particularly in the shallower inter-tidal zone where the maximum water level might be only 1 to 2 m or less. A value of D min = 5 cm was found to work well in recent FVCOM simulations of several tidal estuaries in Georgia and South Carolina [24,25], where the maximum flushing area is defined approximately by the 2-m elevation line. The ratio of D min /max(h B ) equaled 40, sufficiently small to ensure adequate mass conservation in these realistic model applications.

Relationship with kb
The upper bound of I split with respect to vertical resolution is sensitive to D min (Figure 9). For the standard case with kb = 6, α = 7.0 × 10 −4 , and ζ o = 1.5, for example, the upper bound of I split is 10 for D min = 5 cm and about 28 for D min = 20 cm. Keeping the same forcing and geometry but increasing the number of sigma levels kb to 11, the upper bound of I split remains almost the same for D min = 5 cm but drops as D min increases.
for kb = 6 and 11, the thickness of each σ layer varies from 4 cm and 2 cm (at the outer edge of the intertidal zone) to 32 cm and 16 cm (at the edge of the channel), respectively. When the water is drying in the ebb tide, there are at least 5 layers or 2 layers in a viscous bottom layer of Dmin at the wet-dry transition zone for the cases with kb = 11 or 6. Like other coastal ocean models such as POM (Princeton Ocean Model) [37] and ROMs (Regional Ocean Model system) [38], a log-layer dynamics is incorporated into FVCOM, and also this layer could be numerically resolved when the vertical resolution is sufficiently high. Because of the existence of the log layer at the bottom, the model predicts a stronger vertical shear of the current in an additional layer of Dmin for kb = 11 than for kb = 6. This shear directly contributes to vertical mixing and destroys the mass conservation in the TCE in the wetdry transition zone. This is one of the reasons that we found that the up-bound value of Isplit reduces as kb increases in a relatively thicker layer of Dmin. Figure 9. The model-derived relationship of ∆ /∆ with Dmin for the two cases with kb = 6 and 8, respectively. In these two cases, ∆ = 4.14 s, α= 4.0 × 10 −4 , and = 1.5 m.

Relationship with ∆ and ∆
For a given ∆ , the upper bound of Isplit decreases as horizontal resolution increases ( Figure 10). For the standard cases with = 1.5 m, for example, the upper bound of Isplit is cut off at 10 for Dmin = 5 cm and at 28 for Dmin = 20 cm. However, these values drop to 7 for Dmin = 5 cm and to 21 for Dmin = 20 cm when ∆ is decreased from 600 m to 300 m ( Figure  9: upper panel). Similarly, for a given ∆ = 300 m, the upper bound of Isplit increases as ∆ decreases, e.g., Isplit jumps from 7 to 15 for Dmin = 5 cm and from 21 to 45 for Dmin = 20 as ∆ decreases from 4.14 s to 2.07 s (Figure 9: lower panel). In the range of Dmin shown in Figure 9, for the two cases with (∆ )1 and (∆ )2, (Isplit)2 is approximately estimated by . (8) Figure 9. The model-derived relationship of ∆t I /∆t E with D min for the two cases with kb = 6 and 8, respectively. In these two cases, ∆t E = 4.14 s, α = 4.0 × 10 −4 , and ζ o = 1.5 m.
In the case with α = 7.0 × 10 −4 , defining H = 0 at the edge of the channel, the height of the land at the boundary of the inter-tidal zone is 1.4 m. Considering an amplitude of ζ o = 1.5 m, when the intertidal area is entirely covered by the water during the flood tide, the water elevation at the edge of the intertidal zone should be 10 cm. In the actual numerical computation, however, the water elevation is adjusted by an additional D min in the region after becoming wet. For a case with D min = 10 cm, for example, adding D min is equivalent to double the water elevation at the outer edge boundary of the intertidal zone. In this case, for kb = 6 and 11, the thickness of each σ layer varies from 4 cm and 2 cm (at the outer edge of the intertidal zone) to 32 cm and 16 cm (at the edge of the channel), respectively. When the water is drying in the ebb tide, there are at least 5 layers or 2 layers in a viscous bottom layer of D min at the wet-dry transition zone for the cases with kb = 11 or 6. Like other coastal ocean models such as POM (Princeton Ocean Model) [37] and ROMs (Regional Ocean Model system) [38], a log-layer dynamics is incorporated into FVCOM, and also this layer could be numerically resolved when the vertical resolution is sufficiently high. Because of the existence of the log layer at the bottom, the model predicts a stronger vertical shear of the current in an additional layer of D min for kb = 11 than for kb = 6. This shear directly contributes to vertical mixing and destroys the mass conservation in the TCE in the wet-dry transition zone. This is one of the reasons that we found that the up-bound value of I split reduces as kb increases in a relatively thicker layer of D min .

Relationship with ∆L and ∆t E
For a given ∆t E , the upper bound of I split decreases as horizontal resolution increases ( Figure 10). For the standard cases with ζ o = 1.5 m, for example, the upper bound of I split is cut off at 10 for D min = 5 cm and at 28 for D min = 20 cm. However, these values drop to 7 for D min = 5 cm and to 21 for D min = 20 cm when ∆L is decreased from 600 m to 300 m (Figure 9: upper panel). Similarly, for a given ∆L = 300 m, the upper bound of I split increases as ∆t E decreases, e.g., I split jumps from 7 to 15 for D min = 5 cm and from 21 to 45 for D min = 20 as ∆t E decreases from 4.14 s to 2.07 s (Figure 9: lower panel). In the range of D min shown in Figure 9, for the two cases with (∆t E ) 1 and (∆t E ) 2 , (I split ) 2 is approximately estimated by The fact that the upper bound of Isplit drops as ∆ decreases can be easily understood in view of the geometric change of the tracer control element. For example, when ∆ is cut in half, a triangle will be divided into 4 triangles. For a given ∆ , the possibility for individual triangles to change state (switch from wet to dry or vice versa) is much higher in a high horizontal resolution case than in a low horizontal resolution case. Therefore, Isplit should be more restrictive as ∆ becomes smaller. Similarly, for a given ∆ , the possibility for individual triangles to change state should be less as ∆ becomes smaller, so that Isplit should be less restrictive in the case with a smaller ∆ .
In addition to satisfying the criterion of numerical stability, how to select ∆ and how to determine ∆ over the intertidal zone will directly affect the accuracy of the simulated flooding time over the intertidal area. Theoretically speaking, ∆ and ∆ should be chosen in such a way that the numerical advection velocity (defined as ∆ /∆ ) is equal to or smaller than the actual water advection velocity. The failure to satisfy this prerequisite can cause an unrealistic faster flooding speed and exaggerate the horizontal diffusion, thus making the modeled time required for flooding smaller than its actual value.  The fact that the upper bound of I split drops as ∆L decreases can be easily understood in view of the geometric change of the tracer control element. For example, when ∆L is cut in half, a triangle will be divided into 4 triangles. For a given ∆t E , the possibility for individual triangles to change state (switch from wet to dry or vice versa) is much higher in a high horizontal resolution case than in a low horizontal resolution case. Therefore, I split should be more restrictive as ∆L becomes smaller. Similarly, for a given ∆L, the possibility for individual triangles to change state should be less as ∆t E becomes smaller, so that I split should be less restrictive in the case with a smaller ∆t E .
In addition to satisfying the criterion of numerical stability, how to select ∆t E and how to determine ∆L over the intertidal zone will directly affect the accuracy of the simulated flooding time over the intertidal area. Theoretically speaking, ∆L and ∆t E should be chosen in such a way that the numerical advection velocity (defined as ∆L/∆t E ) is equal to or smaller than the actual water advection velocity. The failure to satisfy this prerequisite can cause an unrealistic faster flooding speed and exaggerate the horizontal diffusion, thus making the modeled time required for flooding smaller than its actual value.

Discussion
The flooding and drying numerical experiments conducted in this study indicate that the upper bound of I split =∆t I /∆t E depends on factors related to bathymetry, horizontal/vertical resolution, the external mode time step ∆t E , boundary forcing, and the thickness of the viscous bottom layer. Although the relationship of I split with these factors is based on simple experiments with an idealized tidal channel with an intertidal zone and a simple criterion of the small net change in conservative water properties like heat and salt, the general patterns drawn from these results should provide a helpful guide for modelers who are interested in applying FVCOM to simulate flooding/drying in estuaries and/or the coastal ocean. We highly recommend that any modeler starting such a study first conduct test experiments (like those described here) with spatially uniform temperature and salinity as a check on their model setup and choice of I split . If the model temperature and salinity remain constant over time within small limits, then the mass is conserved in the intertidal zone, and FVCOM will produce accurate temperature and salinity fields in experiments with realistic non-uniform, initial temperature and salinity conditions. If the model temperature and salinity change significantly over time from their uniform initial distribution in the test experiment, then the mass is not conserved, and a smaller value of I split should be considered.
It has been recognized for many years that the Mellor and Yamada [39] level 2.5 turbulence model (called MY-2.5) abruptly shuts down mixing in stratified water when the Richardson number (R i ) rises above 0.25 [40]. When applying MY-2.5 to this idealized estuary (with R i = 0), we found that it produces a substantial vertical shear of vertical eddy viscosity (K m ) and thermal diffusion coefficient (K h ) in the wet-dry transition zone. Without the inclusion of wind forcing, the model-predicted vertical structure of K m and K h reduces to a parabolic shape with a maximum value in the water column no matter how shallow the water is. Thus, the vertical gradient of K m and K h can be very large in the wet-dry transition zone, where the value of D-D min could be just an order of a few centimeters. When this happens, the implicit calculation of the vertical diffusion terms in either momentum or tracer equations can lead to unrealistic velocity or tracer concentration. The unrealistic current speed produced by the MY-2.5 turbulence closure model could be higher than a few meters per second in local wet-dry transition cells. As a result, it can cause the model to blow up numerically within a few time steps.
To avoid this problem, we adjust K m and K h in the wet-dry transition region. A simple way is to use the vertically-averaged values of the model-predicted K m and K h in the wet-dry transition region. This approach works well in the idealized case discussed in this study and an application to natural estuaries in Georgia and South Carolina. In these and many other realistic applications (e.g., [26]), the height of the inter-tidal zone with respect to the height of the river edge is less than 2 m. A more straightforward alternative approach is to use the vertically-averaged values of the model-predicted K m and K h everywhere in the inter-tidal area where it is flushed. Test results show that this approach also works well to avoid the occurrence of unrealistic values of velocity and tracer concentration.
From a physical point of view, the parabolic shape of K m and K h is incorrect in a shallow region with vigorous mixing where the vertical depth is an order of a few centimeters. This parabolic distribution is produced by a zero value condition of q 2 l (where q 2 is the turbulent kinetic energy and l is the turbulent macroscale) at the surface and bottom. Clearly, this boundary condition is not valid in a very shallow water region. This raises the dynamic question of whether MY 2.5 or other turbulence closure models should be applied to a shallow estuary with an extensive intertidal zone. This issue must be addressed in future estuarine model development.

Summary
A 3-dimensional wet/dry point treatment method has been developed for the unstructured grid finite-volume coastal ocean model FVCOM and tested according to the prerequisite of volume conservation. The model equations were used to examine the reasons for numerical errors produced by mass and volume inconsistencies between the external and internal modes during the simulation of the flooding/drying process. These ideas were tested using numerical experiments for an idealized estuary with an intertidal zone.
The model results indicate that in order to ensure volume and mass conservation during flooding/draining cycles, the ratio of internal to external mode time steps (I split ) is constrained by the amplitude of water elevation, bottom slope, and horizontal/vertical resolution as well as the thickness of the viscous bottom layer specified in the model. The upper bound of I split increases with an increase in the bottom slope of the intertidal zone (α) or a decrease in the amplitude of tidal forcing ( ζ o ). For a given ∆t E , the upper bound of I split decreases as horizontal resolution (∆L) increases. Similarly, for a given ∆L, the upper bound of I split increases as ∆t E decreases. Although the upper bound of I split becomes more flexible as the thickness of the viscous bottom layer (D min ) increases, larger values of D min can lead to an unrealistic simulation of the water transport in the intertidal area and make I split sensitive to vertical resolution.
These numerical experiments indicate that MY 2.5 leads to unrealistic velocity shears during the flooding/drying process over the intertidal area, particularly in the wet-dry transition zone where the total water depth can be just a few centimeters. Two simple approaches have been introduced to avoid this problem in the inter-tidal area. The question of a more appropriate turbulent mixing scheme to simulate the flooding/drying process in a shallow estuary or coastal area with an extensive flooding area needs further investigation.