Three-Dimensional Unstructured Grid Finite-Volume Model for Coastal and Estuarine Circulation and Its Application

: We developed a three-dimensional unstructured grid coastal and estuarine circulation model, named the General Ocean Model (GOM). Combining the finite volume and finite difference methods, GOM achieved both the exact conservation and computational efficiency. The propagation term was implemented by a semi-implicit numerical scheme, the so-called 𝜃 scheme, and the time-explicit Eulerian – Lagrangian method was used to discretize the nonlinear advection term to remove the major limitation of the time step, which appears when solving shallow water equations, by the Courant – Friedrichs – Lewy stability condition. Because the GOM uses orthogonal unstructured computational grids, allowing both triangular and quadrilateral grids, considerable flexibility to resolve complex coastal boundaries is allowed without any transformation of governing equations. The GOM was successfully verified with five analytical solutions, and it was also validated when applied to the Texas coast, showing an overall skill value of 0.951. The verification results showed that the algorithm used in GOM was correctly coded, and it is efficient and robust.


Introduction
This study focuses on the development of a new three-dimensional coastal and estuarine circulation model, and we named this model the General Ocean Model (GOM). GOM was developed by combining finite difference and finite volume numerical schemes, taking advantage of the computational efficiency of the finite difference method (FDM), the exact conservation of finite volume method (FVM), and the flexibility of representing complex geometry with an orthogonal unstructured grid system. The advantage of the unstructured grid system over the structured grid system is obvious, but it requires more simulation effort; i.e., the unstructured mesh system well resolves complex boundaries, however, the structured grid system struggles to resolve complex geometries but has a regular structured algebraic equation system and thus has an efficient solution technique. The refining model grid to better represent a complex coastal geometry enforces modelers to use a small simulation time step to ensure numerical stability. Even though it has been common to use a high-performance computing system, implementing either distributed memory message passing interface (MPI) or shared memory open multi-processing (OpenMP), it is also true that it requires a higher simulation cost and longer simulation time as it demands more dramatic grid refinement. Thus, more efficient numerical schemes and algorithms are required to adapt the unstructured grid system. Either significant time constraint with traditional explicit schemes or wave damping with implicit schemes arises when solving shallow water equations. However, the limitation is now well overcome with the semi-implicit approach, the so-called method, which was successfully adapted in several ocean models (e.g., the Unstructured nonlinear Tidal Residual Inter-tidal Mudflat model (UnTRIM) by Casulli and Walters [1]; Stanford Unstructured Nonhydrostatic Terrain following Adaptive Navier-Stokes Simulator (SUNTANS) by Fringer et al. [2]; the Semi-Implicit Cross-Scale Hydroscience Integrated System Model (SCHISM) by Zhang et al. [3,4]); and the Finite Volume Coastal Ocean Model (FVCOM) by Chen et al. [5,6]. Then, another significant bottleneck appears in the nonlinear advection term, and the bottleneck can be successfully removed by using the timeexplicit Eulerian-Lagrangian method (ELM), which is also known as the semi-Lagrangian (SL) method in the field of the atmospheric modeling. The ELM, which is an unconditionally stable scheme-even though it is an explicit method [7]-has been gaining more attention in the ocean modeling community since Casulli and Walters [1] adapted the method in their model, UnTRIM and successfully applied it in the San Francisco Bay area (e.g., [8][9][10]).
Distinct features of well-recognized ocean circulation models, which are actively used in the United States of America, are well summarized by Fringer et al. [11]. Each model introduced by Fringer et al. [11] has different approaches in horizontal/vertical coordinates systems, numerical schemes, and algorithms when solving governing equations based on the model development purpose. Among those approaches, we greatly benchmarked UnTRIM of Casulli and Walters [1], and we developed a new three-dimensional (3D) estuarine circulation model including the following features to apply our model in general coastal water bodies: (1) an unstructured orthogonal triangular and/or quadrilateral horizontal mesh system; (2) a vertical z-grid system; (3) inclusion of wind stress, atmospheric pressure, Coriolis, horizontal/vertical diffusion, and bottom friction; (4) FVM/FDM for equation discretization; (5) ELM for the nonlinear advection equation; (6) the semi-implicit method for tidal propagation; and (7) wetting and drying. The model we developed here, GOM, is based on well-proven numerical techniques; thus, it is robust, accurate, and fast.

Governing Equations in GOM
Using the right-handed Cartesian coordinate system ( , , ) with the vertical origin at the water surface and -axis points upward, the primitive equation of motion, for the incompressible fluid, with the parameterization of viscous stress terms, with the Boussinesq approximation, and with the hydrostatic assumption, can be written as Momentum equations: Continuity equation: where ( , , , ), ( , , , ), and ( , , , ) are the velocity components in the horizontal , , and vertical direction; is time; is the free surface elevation measured from the vertical origin; is pressure; is atmospheric pressure; is the density of the reference fluid; is the gravitational acceleration; is the Coriolis parameter; and ℎ and are the horizontal and vertical eddy coefficients, respectively. Now, we need another form of the continuity equation for the free surface water, and it can be obtained by integrating the original continuity Equation (4) over the depth, applying the kinematic surface (i.e., motion at the surface) and bottom boundary conditions (i.e., motion at the bottom). Then, the original continuity Equation (4) can be successfully rewritten as a new continuity form for free surface flows (i.e., no rigid boundary at the surface), and this form is called as the "free surface equation"

Unstructured Orthogonal Mesh and Index
We used an unstructured orthogonal mesh, both triangular and quadrilateral, in this model GOM. Once we construct the orthogonal meshes for the study site, the grid system has polygons (i.e., the total number of cells or elements), each having three (if it is triangular) or four sides (if it is quadrilateral), i.e., = 3 4, where = 1, 2, … , . The total number of sides in a grid system is defined as , and the length of each side is defined as . Each side of the th polygon is identified by ( , ), = 1, 2, … , , and the corresponding side lengths are defined as ( , ) . Finally, the distance between the centers of neighboring cells, which share the th side, is denoted by ( Figure 1); note that a definition of an orthogonal unstructured grid is well explained in several previous studies (e.g., [2,12]).
As shown in Figure 2, along the vertical direction, a finite difference discretization, which is not necessarily uniform, is adopted; the th vertical layer has a height of Δ , i.e., the distance between levels − 1 and . The horizontal velocities and water surface elevation are defined at staggered locations as follows: The water surface elevation , the density , the salinity , and the temperature are located at the center of the th polygon; note that the salinity and temperature are not yet included in the current version, but we show this in Figure 2 for the clarity and to explain the baroclinic gradient term in the later section. The velocity component normal to each face of a prism is defined at the point of intersection between the face and the segment joining the centers of the two prisms that share the face (i.e., the face center).

Figure 1.
Orthogonal unstructured mesh and its notation used in the General Ocean Model (GOM), where is the th polygon; ( ,1~) is the length of the th side when calling from the th polygon; and is the distance between the centers of neighboring cells that share the th side.

Figure 2.
Location of computational variables, where , , , and denotes the water surface elevation, water density, salinity, and water temperature, respectively, which are defined at the center of the th element; ℎ is the water depth which is defined at the th side; , * and , * are the face normal and tangential velocities which are defined at the th side; and is the vertical layer index.

Momentum Equations in a New Coordinate
Equations (1) through (5) are invariant under solid rotation of the -and -axis on the horizontal plane. If we introduce a new coordinate system, * and * , regarding the cell face, and by using the invariant property of the equations, the horizontal momentum Equations (1) and (2) where * and * are the horizontal velocity components in a new coordinate system, * and * ; note that we omit asterisk * at the * -and * -coordinates. The relations of velocity components between true horizontal velocities, and , and coordinate transformed along each cell face, horizontal velocities, * and * , are = * − * , and = * + * , where is the angle between the new -axis, * , and the traditional -axis, , (Figure 3). For simplicity, * is omitted for the equations from here onwards. is the angle between these two axes.
Each term in Equation (6) can be discretized using different methods for the sake of accuracy and efficiency. When solving these equations, the two most significant bottlenecks arise from the barotropic gradient term and the vertical mixing term. Thus, we used the semi-implicit scheme for the barotropic gradient term and the implicit scheme for the vertical diffusion term. On the other hand, the advection, horizontal diffusion, Coriolis, and baroclinic terms are discretized by explicit methods; note that the semi-implicit scheme is also used in the air pressure term. More details in each term are explained in the following sections.

Barotropic Gradient and Vertical Mixing: (Semi-) Implicit Treatment
For the two implicit scheme applied terms, barotropic gradient and the vertical diffusion terms, the finite difference discretization for the velocity component normal to each vertical face of a prism with a staggered grid system can be derived from Equation (6) and takes the following forms: The barotropic gradient term with a semi-implicit scheme: The vertical diffusion term with the implicit scheme: Then, Equation (6) can be written as where = , + 1, … , , thus and denote the lower and upper limit for the vertical index at -th side and time step . The , is the horizontal velocity component normal to the -th side of the -th mesh (we omit the subscript in equations for convenience), at vertical level and time step ; is an explicit finite difference operator, which includes the remaining terms in Equation (6): nonlinear advection, baroclinic gradient, air pressure, horizontal diffusion, and Coriolis terms. Calculations of the term , will be discussed in the following sections.

Nonlinear Advection: ELM
The Eulerian-Lagrangian Method (ELM) is used for solving the nonlinear advection terms in Equations (6) and (7) to take advantage of its simplicity, the enhanced stability, and accuracy. Most of the fundamental equations in fluid dynamics can be derived from principles in either Eulerian or Lagrangian forms. Consider the diffusion-free non-conservative advection equation in 1D where the ( , ) represents any scalar quantities or velocity vectors, is the velocity field, and ∇ is a gradient operator. Then, the one-dimensional Eulerian advection Equation (11) can be expressed in the Lagrangian form as The mathematical equivalence of Equations (11) and (12) follows from the definition of the total derivative, = + (13) and the definition of the velocity, where = 1, 2, and 3, and * represents a linear interpolant between time step and + 1. When a Lagrangian numerical treatment is applied to Equation (12), the computational grid will be continuously deformed in the general case when is non-constant. For operational advantage, however, we will discretize Equation (12) on a fixed Eulerian grid system. Let us consider the onedimensional time-dependent grid in the Eulerian grid shown in Figure 4. A finite difference scheme for Equation (11) is simply where − = [( − )∆ , ∆ ] and is time step, is any mesh or grid point, is the index at a grid point, and is the Courant-Friedrichs-Lewy (CFL) number. In general, the CFL number is not an integer, therefore, ( − ) is not the index of a grid point, and a proper interpolation formula should be used to define − . The stability, numerical diffusion, and unphysical oscillations of Equation (15) depend on the interpolation formula chosen. If a linear interpolation between ( − 1) and ( ) is used to estimate − , one obtains the first-order upwind scheme. If a quadratic polynomial fit is used to interpolate between ( − 1) , ( ), and ( + 1) , one obtains Leith's [13] method. The ELM is well introduced by other researchers (e.g., [14][15][16]), thus, we do not re-introduce this method in detail here. The final form of the nonlinear term ( ) in the governing equation that we used in GOM is approximated as where and are the horizontal and vertical CFL numbers, respectively.

Coriolis: ELM
The Coriolis term in the equation is treated with the explicit ELM to take advantage of no time limitation for numerical stability. Consider the Coriolis term in Equation (6) = By the explicit ELM, Equation (17) can be discretized as where * is the tangential velocity component in a right-hand coordinate system obtained by the Eulerian-Lagrangian method (see Figure 2). Finally, the Coriolis term ( ) in Equation (6) is discretized by

Horizontal Diffusion: The Finite Volume Cell-Centered Method
The conservation laws of fluid motion may be expressed mathematically in either differential or integral forms. When the integral form of the equation is utilized, the discretization of the equations is the finite volume method. To generalize the method, consider a two-dimensional heat conduction equation where ∇ 2 is the Laplace operator. Then, Equation (20) can be written as a conservative form Let = ( ) and = ( ), then Equation (21) is written as Equation (22) is integrated over an element's area such as quadrilateral or triangular mesh, then where denotes the -th polygon. Subsequently, Green's theorem is applied to the right-hand side of Equation (23). Recall that Green's theorem converts area integrals to line integrals. Thus, Equation (23) is written as Finally, Equation (24) can be approximated as Thus, the horizontal diffusion term ( ) in Equation (6) is approximated as:

Atmospheric Pressure: Semi-Implicit Method
The atmospheric pressure term in Equation (6) is approximated by the semi-implicit, method in time and forward difference in space. Consider the atmospheric pressure term in Equation (6) = − 1 The semi-implicit in time and forward differencing in spatial derivative produces the following finite difference equation of the form , +1 − , and the atmospheric pressure term ( ) in Equation (6) is approximated as The stability analysis indicates that this method is stable when 0.5 ≤ ≤ 1. The values of + are not calculated but provided explicitly as data at each time step.

Baroclinic Gradient: Explicit
The baroclinic gradient term (BG) in Equation (6) can be discretized with the explicit scheme , Then, the final numerical discretization of explicit term in Equation (10) is written as

Surface and Bottom Boundary Treatment for the Vertical Diffusion Term
While the vertical diffusion term in Equation (10) applies for a middle of a water column, new relations are required at both surface and bottom boundary layers to accelerate by the wind stress and to retard by the bottom friction. Using the shear stress, the vertical boundary conditions at the surface and the bottom yield the following formulae: where is the hydraulic radius, and is Manning's roughness coefficient; in shallow estuaries, the hydraulic radius can be approximated by the total water depth .
Similarly, the boundary conditions at the free surface are specified by the prescribed wind stresses where and are the wind stresses in and directions at the free surface; and are the density of the water and the air, respectively; and are the components of wind speed measured at some distance above the free surface; and is the surface drag coefficient. Then, the surface drag coefficient is now calculated by the empirical relationships such developed by either Garratt [17] or Smith [18], which are adopted in this model.
where is the wind speed measured at 10 m above the water surface.

Treatment of Open Boundary Conditions
The tidal open boundary condition for the surface elevation is prescribed by where , , and are the amplitude, period, and phase angle of each tidal constituent, respectively. When open boundary conditions are provided in terms of the surface elevation, the normal velocity component is assumed to be of a zero slope while the tangential velocity component may be either (1) zero, (2) zero slope, or (3) computed from the momentum equations. In this model, it is assumed that the velocity gradients are zero at the open boundary.

Free Surface Equation
To obtain complete local and global mass conservation and stability, the free surface Equation (5) is discretized by the finite volume method. Consider a uniform rectangular mesh in the Cartesian coordinate, as shown in Figure 5. To obtain a semi-implicit finite volume equation, integrate Equation (5) over an area of an element , then Now, a semi-implicit finite volume discretization for Equation (40) gives where Δ = ( ,2) = ( ,4) and Δy = ( ,1) = ( ,3) . If Equation (41) is generalized for any shape of polygons in the unstructured grid system, a semiimplicit finite volume discretization for the free surface equation at the center of each polygon is taken to be [1]: where denotes the area of the -th polygon, i.e., = Δ × ∆ in Figure 5, and , is a sign function of flows at each side, which is defined as [1]  are the length of each cell face.

Solution of Governing Equations
Including wind stress at the surface layer and bottom stress at the bottom layer, and re-arranging the original x-momentum Equation (10), then we have the following tridiagonal matrix form for the -th water column from the surface to the bottom layer at the -th side

Treatment of Flooding and Drying
Once the free surface elevation has been computed throughout the computational domain, before proceeding to the next time step, some of the vertical grid spacings Δ , +1 must be updated to account for the new surface location. As shown in Figure 6, at each time step, the new total water depth at the polygon's sides, +1 , may be defined as According to the new total water depth, the vertical grid spacing Δ +1 should be updated. Thus, an occurrence of zero value for the total water depth +1 implies that all the vertical faces separating prisms between the water column ( , 1) and ( , 2) are dry and may become wet at a later time when +1 becomes positive [1]. When the total water depth is equal to zero, the friction factor at that point will be assumed to be infinity, and, accordingly, the corresponding velocity or across the side of the cell is forced to vanish. The occurrence of a zero value for the total water depth in one side of a cell implies zero velocity or zero mass flux until the total water depth becomes positive i.e., the boundary shorelines, which are varying with time, are defined by the condition of no mass flux. This guarantees mass conservation over the computational domain. An element is considered a dry cell only if the total water depths at all sides are zero.
To reduce computational noise or oscillation due to a very small wet element, a minimum critical dry depth is defined except using "0" in Equation (48). Therefore, elements are considered a wet element when the total water depth is greater than the critical depth and dry elements when the total water depth is less than or equal to the critical dry depth. The drying process can take place not only along the coast but also in interior regions such as shoals. For a shallow estuary, even under moderate wind stresses, some interior points over shoals can be dried completely, whereas the surrounding elements are still wet. Similarly, when the sea surface elevation is at a previously flooded location, this element returns to a dry cell. Then, an isolated dry element will not be turned into a wet cell until at least one of the neighboring elements also turns into a wet element.

Model Verification with Analytical Solutions
In previous sections, Sections 2 to 5, we showed governing equations, grid systems, numerical schemes, and solution algorithms used in GOM. These processes may be similar or identical to other models, such listed by Fringer et al. [11]. However, conveying this information into a computer machine language is all different, and we used a standard FORTRAN 90 for this. Even though these processes are well-proven and robust approaches, any mistake in the final coding process leads to failure. To check the final success of this project, we compared a newly developed model, GOM, to well-known analytical solutions.

Wind Setup
The steady-state wind-induced setup due to constant wind stress in a rectangular basin can be written as where is the setup of the water surface, is the applied wind stress, and are the depth and length of the basin, respectively, and the is the distance from the origin. The horizontal grid used in the wind setup test is 21 × 5 cells with a length of 21 and a width of 5 . The depth of the computational domain is a constant 5 , and the horizontal grid spacing is 1 in both directions. Constant wind stress of 0.1 / 2 (i.e., 1 / 2 ) is applied in the positive -direction, and a time step of 30 s is used in the simulation. Both 2-dimensional (with one vertical layer) and 3dimensional (with 5 vertical layers) setups are used to verify the model. Both analytical and numerical results are calculated/extracted from 0.5, 10.5, and 20.5 from the origin and compared, and both 2D and 3D simulation results reached a steady-state condition in one day and are the same as the analytical results (Table 1 and Figure 7).

Tidal Propagation with Constant Depth
Tide is the most significant phenomenon in the ocean, and thus tidal simulation is one of the most important applications of a coastal and estuarine hydrodynamic model. Lynch and Gray [19] derived the analytical solutions for tidally forced estuaries of various geometries and depths. If we consider only tidal propagation term, i.e., the barotropic gradient term, in the original momentum Equation (1), and the continuity Equation (5) where is a tidal amplitude, is an angular frequency of the given tide, and 0 is the initial water surface elevation measured from the undisturbed surface. With these initial and boundary conditions, the analytical solutions of the one-dimensional shallow wave equations for water surface elevation and velocity are [20]  , is the basin's length, is the number of nodes, is the interested location measured from the origin, and is the time.
To compare the numerical solution with a linear analytical solution, a rectangular basin with a constant water depth of 10 m is considered. The model domain is discretized with 929 mixed elements (quadrilateral elements at the left and triangular elements at the right, ∆ = ∆ =~2 ) and one vertical layer, as shown in Figure 8; note that the only reason we constructed a mixed grid in this problem was to prove our model's capability to handle both triangular and rectangular grids. Input tidal forcing along the open boundary was set with 0.5 m of amplitude and a period of 12.42 h. The simulation time step was set to 30 min, and the test simulations were run with = 1 and = 0.5, respectively. It is noted that when using the implicit numerical scheme, i.e., = 1 , numerical diffusion is introduced. Thus, the numerical solution should correspond to the first mode of the analytical solution, which corresponds to the first terms of Equations (53) and (54). If a semi-implicit scheme, i.e., = 0.5, is used, the numerical solution should be compared to the complete Equations (53) and (54) which include the higher mode solution [20,21].
The model results were extracted from three locations and compared with analytical solutions: ). Model results with = 1, which should be compared with the first term in the analytical solution, agree well with analytical solutions (Figure 9). When = 0.5 is used, which should be compared with the full equation, the model results had a little discrepancy with analytical solutions but mostly agreed well ( Figure 10).

Tidal Propagation with Nonlinear Advection
If nonlinear advection terms are included in the two-dimensional shallow water equations, it is not feasible to obtain an analytical solution. However, if the equations are reduced to one-dimension, we can obtain the harmonic series solution. Including only nonlinear advection and propagation terms, the original momentum Equation (1) and the continuity Equation (5) can be reduced as

Quarter-Annular Harbor Test with a Sloping Bottom
Lynch and Gray [19] derived the analytical solutions for a tidally forced estuary with a flat bottom and a sloping bottom. The quarter-annular test, which contains spatially varying geometry and bathymetry, is one of the well-known test cases for testing the integrated numerical schemes of a developed ocean circulation numerical model. Numerical results will be poor with spurious oscillations or with excessive numerical dissipation if poor numerical schemes are used.
The analytical solutions for the Quarter-annular harbor test are well described by Lynch and Gray [19] and Lynch and Officer [22]. Neglecting nonlinear, Coriolis, and horizontal diffusion terms and assuming linear bottom friction, we can obtain the following vertically averaged equations where is water surface elevation, and are vertically averaged velocities, is gravitational acceleration, and is the linear bottom friction coefficient. With the following boundary conditions: no flow boundaries at the closed end; = 1 , which is the inner radius for the quarter-annular harbor at the open boundary; = 0 , which is the outer radius for the quarter-annular harbor; and periodic tidal forcing = 0 cos( ), where 0 is amplitude, the analytic solutions for surface elevation and horizontal velocity are [19,22]

Wetting and Drying Test on Tidal Flats
Some shallow parts of coastal water bodies become wet and dry periodically in response to the tide, and, thus, the correct reproduction of the wetting and drying of the tidal flats is one of the desirable features of numerical tidal flow models based on shallow water equations. There are several different approaches to solve the wetting and drying process, and the major difference between all the models is the way to determine the drying and wetting cells and depths. The method to determine drying cells in the present model is based on Casulli and Walters [1], as explained in Equation (48).
To validate the wetting and drying scheme implemented in the developed model, we compared our model results with the analytical solution developed by Carrier and Greenspan [24]. The derivation of the analytical solution is well-reviewed in several papers (e.g., [20,25,26]). To compare with the analytical solution, a rectangular basin with a linearly sloping bottom is considered. The length and the width of the basin are set to 55 km and 100 m, respectively. The water depth at the origin, = 0 m, is set to 50 m, and the bottom slope, , is 1:1000; thus, the initial land and water interface is at = 50 km ( Figure 14). The computational domain is divided by 550 × 1 quadrilateral elements (∆ , ∆ = 100 ). A periodic tide with an amplitude 0.2 and a period of 1 h is applied at the open boundary, at = 0 m. Analytical results, which are originally given with non-dimensional variables, are converted to dimensional variables and compared with model results. Both analytical and numerical results are calculated and extracted every 1 6 tidal period and compared. Numerical simulation results show that the analytical and the numerical solutions are in good agreement, showing that the wetting and drying process implemented in GOM works well (Figure 15 and Figure  16).

Model Application to the Texas Coast
Even though a model is verified with some analytical problems, it is sometimes not enough since each analytical problem is focused only on a specific term. For this reason, it is also important to verify the model with a realistic field problem, which includes all the physical aspects. The developed model, GOM, is applied to highly complex geometry, the northern Texas coast, to assess the performance of the developed model. Then, the model results are compared with the results of one of the well-known and widely used ocean circulation models, the Semi-Implicit Cross-Scale Hydroscience Integrated System Model (SCHISM). One of the main reasons why we chose the SCHISM here is its similarity in the horizontal grid system and its numerical scheme, the ELM, which is applied in the nonlinear advection term.

Model Domain and Grid Generation
We developed a model grid from the lower Galveston Bay to Sabine Lake ( Figure 17). We used National Oceanic Atmospheric Administration (NOAA) 'Medium Shoreline' data for the land/sea boundary and 'Coastal Bathymetry 2013' from the Texas Natural Resources Information System (TNRIS) to correctly identify the locations of ship channels and the Intracoastal Waterway. The horizontal grid resolution varies from about 10 m at ship channels to about 3000 m at offshore boundaries. Most of the ship channels are resolved with quadrilateral elements, and otherwise, most of the model domains are covered by triangular elements. The offshore boundary for the model domain was set at about 15 km offshore from the coastline, and the final horizontal grids consist of 210,510 elements and 126,474 nodes.
After the horizontal grid was generated, the digital elevation models (DEMs) were interpolated onto the computational grid using the inverse distance weight (IDW) interpolation method. We used two different sources of DEMs: (1) the ETOPO1 Global Relief Model and (2) the Coastal Relief Model (CRM). The ETOPO1 is a 1 arc-minute global relief model of Earth's surface that integrates land topography and ocean bathymetry, and it covers the entire Gulf of Mexico. The CRM has finer spatial resolution than ETOPO1; it provides one of the 3 arc-second, 1/3 arc-second, or 1/9 arc-second resolution DEMs depending on the area. Some portions of the study area are covered by either 1/3 or 1/9 arc-second DEMs, thus, both data are interpolated onto the grid; note that the minimum water depth was set to 1.0 m throughout the model domain. Figure 17. The model domain and the horizontal grid. The ten blue dots denote included river boundaries, and the twelve red squares show National Oceanic Atmospheric Administration (NOAA) tide stations; the numbers are the corresponding station numbers in Table 2. The upper panel is the zoom-in of the selected coastal area at the lower Galveston Bay.

Forcing Conditions and Model Setup
The developed model was validated for the one-month period in July 2017 and was forced by the tides, river discharge, and atmospheric wind stress. Three NOAA tide stations data, which are relatively close to the open boundary grids, were interpolated to the model boundary grids: Freeport (NOAA 8772447), Galveston Bay Entrance (NOAA 8771341), and Sabine Pass (NOAA 8770822) ( Figure 17 and Table 2). Daily freshwater inflows obtained from the United States Geological Survey (USGS) gaging stations were specified at 10 river boundaries ( Figure 17 and Table 2). These USGS stations are mostly located several kilometers upstream from the model river boundary nodes, but we assumed that there were no additional sources and sinks between the river gages and model river boundary nodes. For the wind stress, the North American Mesoscale forecast system (NAM) 6-hourly reanalysis data were interpolated onto the entire horizontal model grids. Spatially uniform bottom friction with Manning's coefficient of 0.016 was used.
The still water condition was considered at the beginning, i.e., no water surface elevation disturbance and no water movement. The implicitness parameter was set to 0.6, and a twodimensional simulation was considered. Finally, the model simulation time step was set to 600 s. Note that we also ran the SCHISM model with identical model setups to compare the simulation results.

Model Validation
Model simulation results of the developed model, GOM, and the SCHISM were compared with the hourly measured water surface elevation data at twelve NOAA tide stations, which are shown in Figure 17 and Table 2. Figure 18 and Figure 19 show that the results of both models are almost identical and give a good agreement with observed water surface elevations. Note that there are some differences at the beginning of the simulation, but these were caused by applying a tide ramping up option only for the SCHISM simulation.
To better evaluate the model simulation results, we used two types of quantitative error analysis methods: mean absolute error (MAE) and predictive skill (it is also called skill or index of agreement), which was introduced by Willmott [27]. They are defined as follows [28,29]: where and are th modeled and observed data, respectively; is the total data number compared, and ̅ is the mean observed value. The mean absolute error explains how much the modeled data deviate from the observed data, the value of skill explains how much the model can reproduce the observed data, and the value Skill = 1 indicates the perfect agreement [29].
When performing the error analysis, the first 5-day simulation results were excluded since the model required times to be stabilized. Overall, the developed model reproduces water surface elevations well throughout the entire model domain, and the results are quite similar to the SCHISM ( Table 3). The overall MAE is 5.6 cm for the GOM simulation, varying from 2.1 cm to 10.5 cm. The model skill value varies from 0.887 to 0.995, and the overall value is 0.951, indicating that the developed model performs well in a complex geometric study area with realistic forcing conditions. Note that there are some possibilities that we can improve the model simulation results by modulating bathymetry, bottom friction, and tidal boundary conditions. Applying spatially varying bottom friction is a common remedy to increase the model simulation results, and it is better to apply the correctness factors for the tidal amplitude and phase when the reference tide stations are off from the model's boundary nodes. However, we ignored them since the goal of this simulation was not to improve the model simulation results but to validate the developed model's performance.

Conclusions
We developed a three-dimensional numerical model with an orthogonal unstructured grid, either triangular or quadrilateral, for coastal and estuarine circulation, and we named it the General Ocean Model (GOM). This model is implemented with a combined finite difference and finite volume method. To eliminate the major simulation time step constraint which arises when solving shallow water equations, we used a semi-implicit method, which is stable when the implicitness factor, , is between 0.5 and 1.0 for the wave propagation term. To moderate another time constraint, we adopted the Eulerian-Lagrangian method for the nonlinear advection term. In addition to the elimination of time constraints, exact numerical conservation is achieved by using the finite volume method.
We benchmarked an algorithm developed by Casulli and Walters [1], and this algorithm is simple, stable, and efficient. The algorithm is easy to expand from 2D to 3D equations, and by implementing ELM and semi-implicit schemes, this model can be used with fine spatial and vertical resolutions of the grid with relatively large time steps. We aimed to develop a model that can be used in coastal and estuarine regions, and, thus, the wetting and drying process is naturally achieved. In addition, this model can be used for a storm surge simulation since the atmospheric pressure term is included. The developed model, GOM, successfully passed five analytical tests we chose, and finally, it was verified with a realistic application to the Texas coast. The model simulation results were also compared with SCHISM to double-check the model performance.
There was another model that we initially considered to compare with our model, FVCOM by Chen et al. [5]. FVCOM is one of the leading and probably the most widely used finite volume ocean circulation models, thus, it might be worth comparing GOM with FVCOM, and they are as follows: First, FVCOM uses so-called "mode splitting", while GOM solves 3D momentum equations directly from the integral forms, thus, no splitting error exists in GOM. Second, adopting the semi-implicit scheme in the recent FVCOM, users can have more freedom to apply this model with a highly refined unstructured model grid having a relatively large simulation time step. However, there is one difference in the horizontal mesh system when compared to our model-FVCOM only allows users to use a triangular mesh grid, while both triangular and quadrilateral grids can be used in our model, GOM. The study area, the Texas coast, is mostly shallow and fully connected with relatively narrow ship channels and the Gulf Intracoastal Waterway. This ship channel system can be better represented by quadrilateral grids than triangular grids to reduce the total grid size. Note that our model has not yet been parallelized, thus, the total number of the elements in the final grid is relatively important to reduce the total simulation time. Another advantage of using quadrilateral grids over triangular grids is their computational efficiency. We believe that adapting the advantages of both grids, i.e., grid refinement flexibility from triangular grids and efficiency from quadrilateral grids, gives users more accurate results and freedom to use this model. Another important difference in numerical schemes used in our model and FVCOM can be found in the treatment of the nonlinear advection term. Our model adopted the Eulerian-Lagrangian Method, while FVCOM uses the second-order explicit upwind method [5,30]. Even though this model is successfully tested, there are some aspects we want to improve when comparing it to well-recognized ocean circulation models, as mentioned by Fringer et al. [11]. With the Z-grid system (our model uses the current version of this), the model saved computational time in 3-dimensional simulations because the number of vertical grid cells is less in the shallow water area. However, the model requires one big surface layer, and this unwillingly reduces 3D to 2D in shallow water or an inundation area, thus, it is difficult to resolve 3D motions in the shallow area; this drawback will be improved by implementing the sigma ( ) grid. Another missing puzzle in this model, which aims to solve general coastal water motions, is the lack of transport equations. Even though this model includes the baroclinic gradient term, the salinity transport equation is not yet included. Implementation of these in this model will be our next task.