GENERALIZED VERTICAL COORDINATES FOR EDDY-RESOLVING GLOBAL AND COASTAL OCEAN FORECASTS

Abstract : Numerical modeling studies over the past several decades have demonstrated progress both in model architecture and in the use of rapidly advancing computational resources. Perhaps the most notable aspect of this progression has been the evolution from simulations on coarse-resolution horizontal and vertical grids that outline basins of simplified geometry and bathymetry and that are forced by idealized surface fluxes, to fine-resolution simulations that incorporate realistic coastline definition and bottom data on relatively short time scales.


N U M E R I C A L M O D E L I N G S T U D I E S over the past
several decades have demonstrated progress both in model architecture and in the use of rapidly advancing computational resources.Perhaps the most notable aspect of this progression has been the evolution from simulations on coarse-resolution horizontal and vertical grids that outline basins of simplifi ed geometry and bathymetry and that are forced by idealized surface fl uxes, to fi ne-resolution simulations that incorporate realistic coastline defi nition and bottom topography and that are forced by observational data on relatively short time scales (Hurlburt and Hogan, 2000;Smith et al., 2000;Chassignet and Garraffo, 2001;Maltrud and McClean, 2005).
The Global Ocean Data Assimilation Experiment (GODAE) is a coordinated international effort envisioning a global system of observations, communications, modeling, and assimilation that will deliver regular, comprehensive information on the state of the oceans in a way that will promote and engender wide utility and availability of this resource for maximum benefi t to the community.Specifi c objectives of GODAE are to apply state-of-the-art ocean models and data assimilation methods to produce shortrange open ocean forecasts, initial and boundary conditions to extend the predictability of coastal and regional subsystems, and to provide initial conditions for climate forecast models (International GODAE Steering Team, 2000).In the United States, a broad partnership of institutions 1 is addressing these objectives by building global and basin-scale ocean prediction systems using the versatile HYbrid Coordinate Ocean Model (HYCOM) (more information available at http://www.hycom.org).
The choice of the vertical coordinate system in an ocean model remains one of the most important aspects of its design.
In practice, the representation and parameterization of the processes not resolved by the model grid are often directly linked to the vertical coordinate choice (Griffi es et al., 2000).Oceanic general circulation models traditionally represent the vertical in a series of discrete intervals in either a depth, density, or terrain-following unit.Recent model comparison exercises performed in Europe (Willebrand et al., 2001) and in the United States (Chassignet et al., 2000) have, however, shown that the use of only one vertical coordinate representation cannot be optimal everywhere in the ocean.These and earlier comparison studies (Chassignet et al., 1996;Marsh et al., 1996, Roberts et al., 1996) show that all the models considered were able to simulate large-scale characteristics of the oceanic circulation reasonably well, but the interior water-mass distribution and associated thermohaline circulation are strongly infl uenced by localized processes that are not represented equally by each model's vertical discretization.

HYBRID VERTICAL COORDINATES
Hybrid vertical coordinates can mean different things to different people: they can be a linear combination of two or more conventional coordinates (Song and Haidvogel, 1994;Ezer and Mellor, 2004;Barron et al., 2006) or they can be truly generalized (i.e., aiming to mimic different types of coordinates in different regions of a model domain) (Bleck, 2002;Burchard and Beckers, 2004;Adcroft and Hallberg, 2006;Song and Hou, 2006).
The generalized vertical coordinates in HYCOM deviate from isopycnals (constant density surfaces) wherever the latter may fold, outcrop, or generally provide inadequate vertical resolution in portions of the model domain.HYCOM is at its core a Lagrangian layer model, except for the remapping of the vertical coor-dinate by the hybrid coordinate generator after all equations are solved (Bleck, 2002;Chassignet et al., 2003;Halliwell, 2004) (Adcroft and Hallberg, 2006).The ability to adjust the vertical spacing of the coordinate surfaces in HYCOM simplifi es the numeriz σ ρ Schematic of an ocean basin illustrating the three regimes of the ocean germane to the considerations of an appropriate vertical coordinate.Th e surface mixed layer is naturally represented using fi xeddepth z (or pressure p) coordinates, the interior is naturally represented using isopycnic ρ (density tracking) coordinates; and the bottom boundary is naturally represented using terrain-following σ coordinates (after Griffi es et al., 2000).

BOX 1: OCEAN REGIMES AND VERTICAL COORDINATES
stiel School of Marine and Atmospheric Science, Division of Meteorology and Oceanography, University of Miami, Miami, FL, USA.cal implementation of several physical processes (e.g., mixed-layer detrainment, convective adjustment, sea-ice modeling) without harming the model of the basic and numerically effi cient resolution of the vertical that is characteristic of isopycnic models throughout most of the ocean's volume (Bleck and Chassignet, 1994;Chassignet et al., 1996).(Halliwell, 2004).The choice of the vertical mixing parameterization is also of importance in areas of strong entrainment, such as overfl ows (Xu et al., submitted).

COMPUTATIONAL ASPECTS
Before one can numerically solve (i.e., discretize) the model's equations, one must decide on a global projection and on how to treat the singularity associated with the North Pole (the South Pole, being over land, is not an issue).
There are several projections that al-low the Arctic to be included in a global ocean model.The HYCOM global confi guration uses an Arctic dipole patch matched to a standard Mercator grid at 47°N (Figure 2).Because our target horizontal resolution is 1/12° at the equator (i.e., 9 km), the location of the dipoles at 47°N gives us a good resolution at mid latitude (i.e., 7 km) as well as in the Arctic Ocean (i.e., 3.5 km at the North Pole) where the Rossby radius of deformation is smaller.An advantage of this pole-shifting projection, as opposed to most others, is that all the grid points below 47°N remain on regular grid spacing (i.e., x-y grid ratio of 1).
The corresponding numerical array size is 4500 by 3298 with 32 hybrid layers in the vertical.The complete system will include the Los Alamos Sea Ice Model (CICE) (Hunke and Lipscomb, 2004) on  (Halliwell et al., in preparation).Th is fi gure illustrates the transition among pressure, terrain-following, and isopycnic coordinates in 1/25° West Florida Shelf simulations nested within a 1/12° North Atlantic confi guration, and it demonstrates the fl exibility with which vertical coordinates can be chosen and the capability of adding additional vertical resolution.
the same grid.The ocean and ice models will run simultaneously, but on separate sets of processors (a much smaller number for CICE), communicating via an Earth System Modeling Framework (ESMF)-based coupler (Hill et al., 2004).
HYCOM's basic parallelization strategy is two-dimensional domain decomposition (i.e., the entire model domain is divided up into smaller sub-domains, or tiles, and each processor "owns" one tile).A halo region is added around each tile to allow communication operations (e.g., updating the halo) to be completely separated from computational kernels, greatly increasing the maintainability and expandability of the code base.Rather than the conventional oneor two-grid points-wide halo, HYCOM has a six-grid-point-wide halo, which is "consumed" over several operations to reduce halo communication overhead.
For global and basin-scale applications, it is important to avoid calculations over land.HYCOM fully "shrink wraps" calculations on each tile and discards tiles that are completely over land (Bleck et al., 1995).HYCOM goes farther than most structured grid ocean models in land avoidance by allowing more than one neighboring tile to the north and south.Figure 2    Even together, these data sets are insufficient to determine the state of the ocean completely, so it is necessary to exploit prior knowledge in the form of statis-two-dimensional Modular Ocean Data Assimilation System (MODAS) SSH analysis (Fox et al., 2002).The MODAS analysis is an optimal interpolation technique that uses complex covariance functions, including spatially varying length and time scales as well as propagation terms derived from many years of altimetry (Jacobs et al., 2001).The model SST is relaxed to the daily MODAS SST analysis that uses daily Multi-Channel Sea Surface Temperature (MCSST) data derived from the 5-channel Advanced Very High Resolution Radiometers (AVHRR)-globally at 8.8 km resolution and at 2 km in selected regions.
To properly assimilate the SSH anomalies determined from satellite altimeter data, the oceanic mean SSH over the altimeter observation period must be determined.In this mean, it is essential that the mean current systems and associated SSH fronts be accurately represented (position, amplitude, and sharpness).Unfortunately, Earth's geoid is not presently known with suffi cient accuracy for this purpose, and coarse hydrographic climatologies (~1° horizontal resolution) cannot provide the spatial resolution necessary.HYCOM, therefore, uses a mean SSH from a previous fully eddy-resolving ocean model simulation that was found to have fronts in the correct position and is consistent with hydrographic climatologies (Chassignet and Garraffo, 2001).Several satellite missions are either underway or planned to determine a more accurate geoid, but until the accuracy reaches a few centimeters on horizontal scales of about 30 km, the present approach will be necessary.The North Atlantic system runs weekly every Wednesday and produces a 10-day hindcast and a 14-day forecast.The Navy Operational Global Atmospheric Prediction System (NOGAPS) (Rosmond et al., 2002) provides the atmospheric forcing, but in the 14-day forecasts, the forcing linearly reverts toward climatology after fi ve days.During the forecast period, the SST is relaxed toward climatologically corrected persistence of the nowcast SST with a relaxation time scale of onefourth the forecast length (i.e., one day for a four-day forecast).The impact of these choices is discussed by Smedstad et al. (2003) and Shriver et al. (in press).For an evaluation of the North Atlantic system, the reader is referred to Chassignet et al. (2005, in press).
The near real-time North Atlantic basin model outputs are made available to the ocean science community within 24 hours via the HYCOM Consortium data server (more information available at http://www.hycom.org/dataserver)using a familiar set of tools such as OPeNDAP, Live Access Server (LAS), and fi le transfer protocol (FTP).These tools have been modifi ed to perform with hybrid vertical coordinates to provide HYCOM subsets to coastal or regional nowcast/ forecast partners as initial and boundary conditions.The LAS has been implemented with an intuitive user interface to enhance the usability of ocean prediction system outputs and to perform diagnostics.

BOUNDARY CONDITIONS FOR REGIONAL MODELS
The chosen horizontal and vertical resolution for the above HYCOM prediction system only marginally resolves the coastal ocean (7 km at mid latitudes, with up to 15 terrain-following σ coordinates over the shelf), but provides an excellent starting point for even higher-resolution coastal ocean prediction systems.
The model resolution should increase to 1/25° (3-4 km at mid-latitudes) by the end of the decade.An important attribute of the data assimilative HYCOM simulations is, therefore, its capability to provide boundary conditions to regional and coastal models.The nesting of higher-resolution models within the basin-scale HYCOM therefore allows limited-area regional forcings (terrestrial buoyancy inputs, tides), physics (wetting and drying), and coastal geometry (tidal inlets, estuaries) to add value to the larger-scale HYCOM oceanstate estimates.
The quasi-operational regional-scale modeling system developed at the Uni-  (Blanton, 2003) uses the fi nite element coastal ocean model QUODDY (Lynch et al., 1996).
Terrestrial buoyancy inputs to the continental shelf, strong tides, and a vigorous western boundary current contribute to the complexity of this region.To simulate the density-dependent dynamics, the UNC-SAB modeling system is nested within the HYCOM GODAE near-realtime system.Boundary and initialization data for the SAB regional-scale model are obtained from the HYCOM GODAE Live Access Server and are mapped to the fi nite element regional model domain.
For each forecast, the system spins up the  .Th e three-dimensional solution for this date is used to initialize the University of North Carolina-South Atlantic Bight (UNC-SAB) fi nite element modeling system (middle, every 5 th model vector).Th e surface temperature and velocity are shown, after addition of the regional tides and atmospheric fl uxes.Th e limited-area fi nite element implementation (right, every 3 rd model vector) includes the estuary and tidal inlets along the Georgia/South Carolina coast and extends to the shelf-break.Th is snapshot also shows surface temperature and velocity and is for three days after initialization.Th e color scales between the HYCOM and nested plots are not the same.
to the tides.
Evaluation of the near real-time HY-COM outputs relative to available observations in the SAB consists of comparisons to National Ocean Service water levels along the coast and mid-shelf temperature from an established continental shelf observational network (SABSOON) (Seim, 2000).Variability in observed coastal water levels is due to tides, windstress-driven and other lower-frequency fl uctuations (deep-ocean contributions to shelf-wide sea level).In the SAB, tides account for at least 90% of the total water level variability.Figure 5a

1
University of Miami, Naval Research Laboratory (NRL), National Aeronautics and Space Administration/Goddard Institute for Space Studies (NASA/GISS), National Oceanic and Atmospheric Administration/National Centers for Environmental Prediction (NOAA/NCEP), National Oceanic and Atmospheric Administration/Atlantic Oceanographic Meteorological Laboratory (NOAA/AOML), National Oceanic and Atmospheric Administration/Pacifi c Marine Environmental Laboratory (NOAA/PMEL), Planning Systems, Inc. (PSI), Fleet Numerical Meteorology and Oceanography Center (FNMOC), Naval Oceanographic Offi ce (NAVOCEANO), Service Hydrographique et Océanographique de la Marine (SHOM), Laboratoire des Ecoulements Géophysiques et Industriels (LEGI), Open Source Project for a Network Data Access Protocol (OPeNDAP), University of North Carolina, Rutgers, University of South Florida, Fugro GEOS, Roff er's Ocean Fishing and Forecasting Service (ROFFS), Orbimage, Shell, ExxonMobil.Oceanography Vol. 19, No. 1, Mar. 2006 119 Because none of the three main vertical coordinates (depth, density, and terrain-following) (see Box 1 for details) provide universal optimality, it is natural to envision a hybrid approach that combines the best features of each vertical coordinate.Isopycnic (density-tracking) layers work best for modeling the deep stratifi ed ocean, levels at constant fi xed depth or pressure are best to use to provide high vertical resolution near the surface within the mixed layer, and terrain-following levels are often the best choice for modeling shallow coastal regions.In HYCOM, the optimal vertical coordinate distribution of the three vertical coordinate types is chosen at every time step.The hybrid vertical coordinate generator makes a dynamically smooth transition among the coordinate types using the continuity equation.
The default confi guration of HYCOM is isopycnic in the open stratifi ed ocean, but makes a dynamically smooth transition to terrain-following coordinates in shallow coastal regions and to fi xed pressure-level coordinates in the surface mixed layer and/or unstratifi ed seas (Figure 1).In doing so, the model takes advantage of the different coordinate types in optimally simulating coastal and openocean circulation features.A user-chosen option allows specifi cation of the vertical coordinate separation that controls the transition among the three coordinate systems.The assignment of additional coordinate surfaces to the oceanic mixed layer also allows the straightforward implementation of multiple vertical mixing turbulence closure schemes

Figure
Figure 1 illustrates the transition among pressure, terrain-following, and isopycnic coordinates in 1/25° West Florida Shelf simulations nested within a 1/12° North Atlantic confi guration, and it demonstrates the fl exibility with which vertical coordinates can be chosen and the capability of adding additional vertical resolution.The original vertical discretization used in the 1/12° North Atlantic confi guration is compared to

Figure 1 .
Figure 1.Cross sections of layer density and model interfaces across the West Florida Shelf in a 1/25° West Florida Shelf subdomain covering the Gulf of Mexico east of 87°W and north of 23°N and embedded in a 1/12° Atlantic basin HYCOM simulation(Halliwell et al.,  in preparation).Th is fi gure illustrates the transition among pressure, terrain-following, and isopycnic coordinates in 1/25° West Florida Shelf simulations nested within a 1/12° North Atlantic confi guration, and it demonstrates the fl exibility with which vertical coordinates can be chosen and the capability of adding additional vertical resolution.
shows the current tiling for the 1/12° global domain with equal sized tiles that (a) allows rows to be offset from each other if this gives fewer tiles over the ocean and (b) allows two tiles to be merged into one larger tile if less than 50% of their combined area is ocean.The Mediterranean region of Figure 2 illustrates both these optimizations.In the example presented in Figure 3, the global model is initialized from an oceanic climatology of temperature and salinity.The model is forced by prescribed climatological atmospheric wind, thermal, and precipitation fi elds while evaporation is computed using the modeled sea surface temperature (SST).There is also a weak relaxation to climatological surface salinity.These simulations use a simple thermodynamic sea-ice model (CICE is running standalone on the global Arctic patch grid, but we are awaiting ESMF-based coupling between HYCOM and CICE before running a coupled case).Even without ice dynamics, the seasonal cycle of ice coverage is good overall (not illustrated).One reason for the good agreement is that the atmospheric forcing is based on an accurate ice extent, which provides a strong tendency for the ocean/sea-ice system to form ice appropriately.Figure 3 compares the sea surface height (SSH) variability from the climatologically forced 1/12° global HYCOM to the Oct. 1992-Nov.1998 SSH variability based on Topex-Poseidon, ERS-1, and ERS-2 altimeter data (derived by Collecte Localisation Satellite (CLS), France).Overall, the modeled regions of high variability are in reasonably good agreement with the observations, especially in the Antarctic Circumpolar Current region.The equatorial Pacifi c is an exception because the altimeter data include interannual variability not present in the model (e.g., the large variability associated with the 1997-1998 El Niño).OCEAN PREDICTION Although HYCOM is a relatively sophisticated model that includes a large suite of physical processes and incorporates numerical techniques that are optimal for dynamically different regions of the ocean, data assimilation is still essential for ocean prediction because (a) many ocean phenomena are due to nonlinear processes (i.e., fl ow instabilities) and thus are not a deterministic response to atmospheric forcing, (b) errors exist in the atmospheric forcing, and (c) ocean models are imperfect, including limitations in numerical algorithms and in resolution.Substantial information about the ocean surface's space-time variability is obtained remotely from instruments

Figure 2 .
Figure 2. Current tiling for the 1/12° global domain with equal-sized tiles that (a) allows rows to be off set from each other if this gives fewer tiles over the ocean and (b) allows two tiles to be merged into one larger tile if less than 50 percent of their combined area is ocean.Th e Mediterranean region especially illustrates both these optimizations; out of the original 1152 (36 by 32) approximately equalsized tiles, 371 tiles are entirely land and are discarded, leaving 781.Th e Arctic "wraps" between the left and right halves of the top edge.

Figure 3 .
Figure 3.Comparison between the (observed) Oct. 1992-May 2005 sea surface height (SSH) variability based on Topex-Poseidon, ERS-1, and ERS-2 altimeter data (top) (derived by Collecte Localisation Satellite [CLS], France) and three years of (modeled) sea surface height (SSH) variability from the climatologically forced 1/12° global HYCOM (bottom).Overall, the modeled regions of high variability are in good agreement with the observations, especially in the Antarctic Circumpolar Current region.Th e equatorial Pacifi c is an exception because the altimeter data include interannual variability not present in the model, for example, the large variability associated with the 1997-1998 El Niño.
To increase the predictability of coastal regimes, several partners within the HYCOM consortium are developing and evaluating boundary conditions for coastal prediction models based on the HYCOM data assimilative system outputs.The inner nested models may or may not be HYCOM, so the coupling of the global and coastal models must be able to handle dissimilar vertical grids.Coupling HYCOM to HYCOM is now routine via one-way nesting (Zamudio et al., in preparation).Outer model fi elds are periodically interpolated to the horizontal mesh of the nested model (specifi ed by the user, but typically daily) and stored in an archive fi le.The number of coordinates can be increased to augment the vertical resolution of the nested model and to ensure that there is suffi cient vertical resolution to resolve the bottom boundary layer.The nested model is initialized from the fi rst archive fi le; the entire set of archives provides boundary conditions during the nested run, ensuring consistency between initial and boundary conditions.Coupling HY-COM to other fi nite difference models, such as the Navy Coastal Ocean Model (NCOM) or the Regional Ocean Model System (ROMS), has already been demonstrated, and coupling of HYCOM to unstructured grid/fi nite element models is in progress.We now describe the use of near-realtime HYCOM nowcasts and forecasts as boundary and initial-condition providers to a nested coastal simulation in the South Atlantic Bight (SAB) region of the eastern U.S. coast.The 1/12° North Atlantic HYCOM does not necessarily include all forcing and physics relevant at the coastal region scale (e.g., lack of tidal forcing in the present simulation).
regional tides for fi ve model days, with the density fi eld held fi xed.The atmospheric (buoyancy and momentum) and river fl uxes are turned on.The timing is synchronized such that the fl uxes are active for one day at the time the HYCOM fi elds are valid.An example of this nested system is shown in Figure 4.The HYCOM nowcast for November 28, 2005 is used to initialize the regional model domain, onto which the tides and high-resolution atmospheric fl uxes are applied.The resulting solution is used to initialize and drive the limited-area, estuary-resolving fi nite element implementation.The effects of both the Gulf Stream, as provided by the HYCOM initial and boundary conditions, and the local tides are seen in the limited-area model (Figure 4, right).Strong along-shelf and poleward fl ow is seen at the shelfbreak.The poleward, offshelf-directed fl ow on the shelf is due

Figure 4 .
Figure 4. (Left) 1/12° North Atlantic HYCOM-GODAE sea surface temperature (SST) nowcast for November 28, 2005.Th e three-dimensional solution for this date is used to initialize the University of North Carolina-South Atlantic Bight (UNC-SAB) fi nite element modeling system (middle, every 5 th model vector).Th e surface temperature and velocity are shown, after addition of the regional tides and atmospheric fl uxes.Th e limited-area fi nite element implementation (right, every 3 rd model vector) includes the estuary and tidal inlets along the Georgia/South Carolina coast and extends to the shelf-break.Th is snapshot also shows surface temperature and velocity and is for three days after initialization.Th e color scales between the HYCOM and nested plots are not the same.
Figure 5. (a) Water level comparison for two stations in the SAB.Daily HYCOM best-estimate water levels (blue) and observed NOS water levels (red) are shown for Fort Pulaski, Georgia and Virginia Key, Florida.Th e two stations are arbitrarily off set for clarity.(b) Observed, near-surface temperature (blue) and HYCOM mixed-layer temperature (red) for the SABSOON mid-shelf station R2.Th e root mean square (RMS) error at this location is 1.4°C.Jun03 Sep03 Dec03 Mar04 Jun04 Sep04 Dec04 Mar05 Jun05 -0.2 and for the fact that there is a non-