# Numerical Modeling Approaches

All hydrodynamic models are based upon the Newtonian law that acceleration is the result of imposed forces. However, the full Navier-Stokes equations of continuum fluid motion cannot be solved without some simplifications. Numerical models discretize, or solve simplified versions of these equations at defined points in space and time and simulate unresolved features of the fluid flow, such as turbulence, with parametrized approximations. Although many different approaches have been formulated to perform the discretization and simplification of the continuum equations, all models simulate the same fundamental processes.

### Hydrodynamic Approximations

The Navier-Stokes equations include processes that occur at time and space scales that are much smaller than most real-world coastal and oceanographic problems. The equations are therefore simplified before they are discretized. These simplifications broadly rely on characteristics of ocean and estuarine systems such as that the horizontal scale of problems is much larger than the vertical scale, that is, the oceans are much wider than they are deep. They also focus on spatial scales that are large or slow enough to be influenced by the Earth's rotation. These include simplifications associated with the use of an average or mean density except where the vertical response of the fluid to gravity is computed (the Boussinesq approximation). Some side effects of this approximation are that almost all ocean models are volume rather than mass conserving, and also incompressible, which removes the potential for sound waves. Many models also use the hydrostatic approximation, which assumes that vertical accelerations are not forced by vertical pressure gradients. This assumption is valid when the horizontal scales of ocean flows are much greater than vertical scales. However, new 'nonhydrostatic' models are being developed to simulate problems with small spatial scales, such as convection and outflows and even at larger scales (e.g., MIT general circulation model). The rigid lid assumption removes sea-surface perturbations including gravity waves, but most modern hydrodynamic models no longer use this assumption, which would exclude most coastal phenomena where tides are important.

The basic equations that most ocean models use, including the hydrostatic and Boussinesq approximations in Cartesian coordinates, are shown below with an explanation of the symbols and physical meaning of each term:

Here, eqns — are the momentum equations for describing changes to the momentum of the ocean at a fixed Cartesian point. a/x, v/y, w/z are the zonal or east-west, meridional or north-south, and vertical or up-down velocity/distance components, respectively. The terms on the left side of the momentum equations describe the evolution or rate of change over time of the velocity field at a fixed point. Note that this evolution over time is expressed as the material derivative, that is, it includes the effects of gradients in velocity that flow past a fixed point.

Equations  and  also include terms to describe the Coriolis acceleration: fa or fv. The Coriolis terms represent the effect of the Earth's rotation on the ocean, which appears like a 'force' in the equation, but is only an effect of viewing the oceans from a reference frame that is rotating. The value of f (eqn ), depends upon the rate of rotation of the Earth, and the distance from the equator, e, the latitude. Many models simplify or linearize this equation in a further simplification called the /3-plane approximation.

On the right-hand side of the momentum, eqns - are the forces that drive changes in the velocity fields, p is the local pressure, and p is the density of the water. Where p0 is used, an average density is chosen for simplicity. Only in the vertical momentum equation is the actual density, p, used. These terms that involve the horizontal change in pressure are called the pressure gradient terms, and with the Coriolis terms represent the first-order balance (geostrophic balance) at large spatial scales. \$t represents a body force due to the gravitational pull of lunar and solar motions that drives the ocean tides. The tides are generally used in coastal applications, but not in open ocean or climate models. KM, KT, and KS represent the vertical kinematic viscosity or eddy diffusion for momentum, temperature, and salt, respectively. These terms parametrize the effects of turbulence on vertical mixing of the water column. AM, AT, and AS represent horizontal eddy diffusion of momentum, heat, and salt, respectively. These terms parametrize the effects of friction at the sea surface or sea floor, as well as include the effects of horizontal unresolved turbulent stirring. Simplifications to the equations (discussed above) cause the vertical momentum equation  to appear quite different from eqns  and . Here, the local vertical pressure gradient is balanced by the gravitational acceleration, g. The significance of this difference is that the horizontal equations can often be numerically solved as a two-dimensional problem for the sea-surface height, and then the vertical variations of U, V, S, and T can be solved as a one-dimensional profile problem at every location in the horizontal.

Equations  and  express the conservation of heat and salt, analogous to the momentum equations. T and S are temperature and salinity, respectively. Terms on the left side of the equation express changes in the local temperature or salinity, while those on the right are due to unresolved turbulence stirring together thermal and salinity gradients. St and Ss are source and sink terms that include, for example, radiative inputs of heat at the sea surface for example, or changes to the salt balance due to evaporation and precipitation of freshwater that are represented in many models as a salt flux rather than a change to total mass of fluid in the system. Finally, eqn  expresses the nonlinear and complex dependence of seawater density on temperature, pressure, and salinity.

### Horizontal Grids

Generally, the horizontal model grid can be classed as structured (all cells have the same number of sides and neighbors, Figure 1a) or unstructured (using geometrical shapes optimally pieced together to represent complex geometries, Figure 1b). Each method has advantages and disadvantages. Structured grids are -76.4 -76.3 -76.2 -76.1 -76 -75.9 -75.8 -76.4 -76.3 -76.2 -76.1 -76 -75.9 -75.8

Longitude Longitude

Figure 1 (a) Structured vs. (b) unstructured horizontal grids for the lower Chesapeake Bay.

-76.4 -76.3 -76.2 -76.1 -76 -75.9 -75.8 -76.4 -76.3 -76.2 -76.1 -76 -75.9 -75.8

Longitude Longitude

Figure 1 (a) Structured vs. (b) unstructured horizontal grids for the lower Chesapeake Bay.

employed in most large-scale and global hydrodynamic model applications because these models do not generally need to resolve intricacies of the coastal topography and because structured grids are numerically and conceptually simpler. Spatial variations in grid resolution are commonly employed in structured grids, primarily targeted at better resolving western boundary current regions or the equatorial domain; however, these higher-resolution regions must extend throughout the domain in a regular fashion, which may lead to computational waste.

Unstructured grids have been more widely applied in the coastal ocean where they can resolve complexities of the topography and provide fine control over the spatial variation of grid resolution. Unstructured grids allow resolution of small coastal features, down to the scale of docks and jetties while also containing large elements that span tens of kilometers allowing simultaneous resolution of the tidal dynamics of coastal oceans. Two problems associated with the unstructured grids are inaccurate representation of the geostrophic balance (i.e., the pressure gradient and Coriolis terms), and the potential for unphysical wave scattering associated with highly variable grid spacing.

### Vertical Grids

Much recent effort has focused on the vertical representation of the model grid. The problem of accurately modeling the horizontal variation of the vertically variably pressure field has been solved, with varying success, using a variety of vertical discretization schemes. The major groupings fall into: (1) fixed vertical level models, where the thickness of each vertical layer is uniform and does not vary in time (Figure 2a); (2) sigma-coordinate  P = P2 P - P3 P - P4 J P - P5 /

Figure 2 Schematic vertical sections of a (a) fixed vertical level grid structure, (b) a sigma coordinate vertical grid structure, and (c) a density coordinate grid structure where the latter has temporally and spatially varying thickness.

or terrain-following models, which resolve the vertical structure based on a modal decomposition over depth (Figure 2b); and (3) density-coordinate models, which are composed of layers of uniform density with temporally and spatially varying thickness (Figure 2c). Models that merge these different coordinate systems into a hybrid model and target them to regions where they function most effectively have also been developed (e.g., the Hybrid Coordinate Ocean Model, HYCOM). Table 1 includes models from each of these classes, but is not an exhaustive listing.

Generally, fixed vertical coordinate models (e.g., the MOM) are computationally simple, and they provide a straightforward representation of the equation of state and the pressure gradient. The surface mixed layer can also be well resolved with a specified number of fixed vertical levels. Disadvantages include difficulty calculating diffusion and advection along inclined density surfaces (which has been solved through complex algorithms), and the representation of topography and boundary layers associated with topography. Sigma-coordinate (or s-coordinate) models (e.g., the Princeton Ocean Model, POM, and the Regional Ocean Modeling System, ROMS) have been widely used in coastal applications due to their natural representation of bottom topography and ability to resolve the large shear of the turbulent bottom boundary layer. They have a similar but exacerbated problem representing advection and diffusion along inclined density surfaces, and they have problems representing the pressure gradient term on sloping isopycnal surfaces due to the potentially steep angles of intersection of the coordinate surfaces with density or geopotential surfaces. Density coordinate models (e.g., the MICOM) have variable vertical grid spacing in both time and space, providing a natural representation of fronts and of advection of tracer substances which occur primarily along isopycnal surfaces. The horizontal pressure gradient is well represented, as is topography. Disadvantages of density coordinate models include challenges of representing the equation of state so that the vertical coordinate is not multivalued, and of representing the surface mixed layer which tends to be poorly resolved in the vertical because it is inherently unstratified or lacks vertical changes in density. There are also challenges associated with the computational aspects of allowing density layers to vanish, which is why this method is generally not used in coastal seas and estuaries, where density is highly variable in space and time. HYCOM uses fixed coordinates in the upper ocean mixed layer and potentially in the deep ocean where weak stratification can result in very coarse resolution, density coordinates in the ocean interior, and sigma coordinates below the mixed layer in the coastal ocean. 