Modeling Ocean Circulation

Albert J. Semtner*

Abstract: Ocean numerical models have become quite realistic over the past several years as a result of improved methods, faster computers, and global data sets. Models now treat basin-scale to global domains while retaining the fine spatial scales that are important for modeling the transport of heat, salt, and other properties over vast distances. Simulations are reproducing observed satellite results on the energetics of strong currents and are properly showing diverse aspects of thermodynamic and dynamic ocean responses ranging from deep-water production to El Niño. Now models can represent not only currents but also the consequences for climate, biology, and geochemistry over time spans of months to decades. However, much remains to be understood from models about ocean circulation on longer time scales, including the evolution of the dominant water masses, the predictability of climate, and the ocean's influence on global change.

---

The ocean is now well known to play a dominant role in the climate system because it can initiate and amplify climate change on many different time scales. The best known examples are the interannual variability of El Niño (1) and the potential modification of the major patterns for oceanic heat transport as a result of increasing greenhouse gases (2). Yet the ocean has been very much undermeasured for most of the history of ocean science. Even though systematic observations began in the 1880s with pioneering observations by Nansen and others (3), the seagoing and theoretical efforts were mainly oriented toward describing large-scale circulation (4), which was often regarded as steady for lack of more detailed information. It was not until the 1960s, when long-distance tracking of drifting buoys at mid-depth showed currents to be highly variable on quite small spatial scales (5), that oceanographers became aware of the immensity of their task.

In dealing with an undersampled natural fluid system, it is extremely helpful to invoke the simulation capability of supercomputers to improve understanding of basic processes and their interconnectedness, as well as to help interpret sparse observations. This article describes the development of ocean modeling as it has been influenced by improved methods, more powerful computers, and data from the latest measuring systems. Recent successes that pave the way to better physical understanding and more reliable ocean prediction are emphasized, although issues of running coupled atmosphere-ocean models are not discussed. Some remaining problems are highlighted briefly.

The Evolution of Models

An early paradigm. A first significant modeling effort extended over the decade of the 1960s at the Geophysical Fluid Dynamics Laboratory (GFDL) (now located at Princeton University and operated by the National Oceanic and Atmospheric Administration). Bryan built an ocean model that was to be used with atmospheric models to study climate. In rapid succession, he and co-worker Cox developed a 2D model, a 3D box model, and then a model of full circulation, with variable density as well, for the world ocean with its complex coastline and bottom topography (6).

The models used space-centered explicit finite differencing, and the last version used a variable number of vertical levels to represent the complex geometry. All of the models used an assumption known as the "rigid lid" approximation to eliminate high-speed external gravity waves and allow a longer time step: the upper surface was held fixed but could support pressure changes related to waves of lower speed and currents of interest. As a result, ocean tides and other waves having the speed of tsunamis were filtered out. Use of the rigid lid created an ancillary 2D equation at every time step to obtain the vertically averaged component of the currents. Bryan's 1969 model used partial differential equations, formulated in rotating spherical coordinates, for horizontal velocity, temperature, and salinity, plus hydrostatic balance, incompressibility, and an empirical equation for density. Processes occurring on scales smaller than the grid size were represented by horizontal and vertical viscosity and diffusion, using four different coefficients.

The first application with specified global geometry was done in the early 1970s (7). Cox designed a 2° latitude-longitude grid with up to 12 vertical levels at each point. The model was run for hundreds of hours on the fastest available computer but had only reached a few simulated years after many months. Because the model had been initialized with smoothly gridded fields of the observed temperature and salinity, it successfully diagnosed currents that were consistent with those fields, but then it began to evolve toward a new state. Realizing the impossibility of reaching thermodynamic equilibrium at this grid resolution, Bryan, Cox, and others turned the model formulation to coarser grid (5°) global studies and to regional studies. However, currents in the global studies were far too broad and sluggish, and pathways throughout the global ocean were omitted from regional studies. However, because of the generality of this "GFDL model," it was sought out and became a standard tool for ocean simulations at many other institutions.

Retooling for computational limitations. In the mid-1970s, there was growing awareness that most ocean currents have cross-stream dimensions that are a fraction of a degree, that is, less than 100 km. Most oceanic high-frequency motions are channeled into slowly varying ocean currents over a distance equal to the speed of the internal gravity wave divided by the local Coriolis frequency. This length scale, called the internal radius of deformation, varies from about 100 km in the tropics to 10 km or less at very high latitudes. The currents then meander, and some produce closed or cutoff circulations on scales of tens of kilometers. These features, called mesoscale eddies, are the oceanic analogs of atmospheric storms. They are known from instability theory to form at a scale near the radius of deformation, but as they develop, they scatter and mix momentum and properties in ways that elude theoretical description. Subsurface drifters show these disturbances to be common throughout the world ocean; in fact, they account for the bulk of the kinetic energy in the ocean (8).

If strong currents and eddies are to be treated adequately in numerical models, then the horizontal grid size has to be small enough to represent or resolve them, meaning that grid spacing should be approximately 20 km in middle latitudes. This could only be achieved effectively in the mid-1970s by shrinking the domain sizes to be much smaller than almost all ocean basins, reducing the number of levels, and further filtering the equations in advance to remove internal gravity waves and allow only the major currents and low-frequency eddies. The resulting two- and three-layer quasi-geostrophic models of Holland, McWilliams, and others became workhorses to understand how eddies affect the circulation in box domains (9) and in channels (10). These models were adiabatic ones with fixed-density layers of varying thickness and were driven only by specified winds. Usually they had a scale-selective viscosity and diffusion to reduce damping.

By means of these adiabatic models, the origin of mesoscale eddies in the strong currents, the complicated mixing of properties produced by eddies, and the generation of deep time-mean currents by vertical eddy transfers were able to be understood, and to a far better extent than by observations alone. In addition, one adiabatic layered model by O'Brien and his students, which did retain internal gravity waves so that equatorial and coastal problems involving these waves could be treated, led to an initial understanding of El Niño in terms of those waves (11).

Toward realism with more choices. By the mid-1980s, progress had been made in simulating more aspects of real-world ocean circulation, much of it resulting from increases in computing power. Quasi-geostrophic models were applied to larger domains and run longer and for more cases, thereby giving a better idea of how physical parameters affect model dynamics. Models for the Indian Ocean and progressively larger fractions of the Pacific were constructed following O'Brien's approach and were run with the observed time-varying wind patterns. Equatorial currents related to the Indian monsoon and the Southern Oscillation were better understood as a consequence. The GFDL-type models were able to have eddies marginally included in small, idealized domains, which revealed how eddies both transport and homogenize certain properties (12). Global GFDL-type models could be run for centuries of simulation with grids of about 3° spacing, which demonstrated that the complete structure of temperature and salinity could be crudely simulated but with a systematic overestimate of the depth of the temperature structure (13).

During this time, a number of alternate model formulations arose, such as ones using curvilinear coordinates to follow the local coastline (14) and with different methods for treating the varying depths (15). In a particularly promising alternative, Bleck used coordinate surfaces of constant density that are able to move freely within the water column (16); this method has the potential to represent the preferred spreading of water masses along such isopycnal surfaces and to reduce excessive horizontal diffusion that can occur in the GFDL model.

In the late 1980s, simulations could finally be undertaken using the GFDL formulation with eddies marginally resolved over extensive domains and with observed winds and some atmospheric influence on density. These studies of the Southern Ocean south of latitude 25° (17), the North Atlantic (18), and the World Ocean without the Arctic (19) were hailed as groundbreaking because they were the first ones of planetary scale that were realistic enough to invite side-by-side comparison with data. The grid sizes of these models were 1/4°, 1/3°, and 1/2°, respectively, and the models were run for decades of simulated time. Output from each of the simulations was scrutinized by a number of investigators who were independent of the model development; the collective efforts became known as community modeling. Each of the three community efforts has produced more than 20 scientific papers dealing with model evaluation and interpretation of physical results. Planetary-scale studies were also conducted with adiabatic models for the Pacific and Atlantic, mainly for the purpose of developing shorter term forecast models but also with applications to wind-forced events like El Niño (20). However, nonadiabatic changes in density such as in the GFDL or isopycnal models are needed to reproduce vertical circulations and associated heat transports that govern the long-term climate.

A free surface and less smoothing. Early in the 1990s, the fraction of computer time devoted to the 2D ancillary problem associated with the rigid lid approximation was becoming excessive in cases of large domains and resolved eddies. For this reason, and to allow options to include tidal effects or height data from satellites, methods were developed to predict the height and pressure of the ocean surface directly. The goal was to retain the free surface without incurring a tremendous reduction in allowable time step due to fast-moving external gravity waves. A straightforward method was first developed to treat the free surface and the vertically averaged velocity using many small steps in time for each single step of the full 3D model (21). Another method developed at Los Alamos National Laboratory solves the same 2D equations using an implicit method for the free surface (22). Both methods are quite efficient.

Investigators from Los Alamos (Smith and Dukowicz) noted that a major benefit derives from either type of free-surface approach. The bottom topography can be specified independently at each location and without regard for how severely it may vary from point to point. In contrast, the rigid lid approach is ill behaved without one or more passes of a smoothing operator over the initial depth field. The investigators at Los Alamos also did a thorough reevaluation of existing methods and made several changes to increase the allowable time step. Moreover, their methods were developed to be used on emerging massively parallel computers. Finally, they initiated a treatment of horizontal pressure forces that doubled the time step and still fit smoothly into the free-surface method. The net speedup of the GFDL formulation as a result of their efforts was a factor of 4 in typical applications.

Two recent simulations have capitalized on the new GFDL free-surface formulations applied to the global ocean (23, 24). Both use Mercator grids that progressively decrease grid size at high latitudes; both models have a maximum of 20 levels for the unsmoothed topography; and they are driven by high-frequency observed winds of 1985 to 1995 and by monthly atmospheric heating. Averaged over their full latitude ranges, the models have 1/4° and 1/6° grid spacing, based on equatorial grid sizes of 2/5° and 0.28°, respectively. Another recent free-surface simulation, this time for the North Atlantic and using Bleck's model with isopycnal coordinates, has a Mercator grid based on an equatorial grid of 1/12° (25). Thus, models are approaching the radius of deformation in grid spacing for global modeling and are surpassing it in one basin-size grid.

Recent related modeling studies that compare simulated eddy energies with observed satellite estimates (26) and examine the grid size needed in local models of the unstable Gulf Stream (27) suggest that small grid sizes such as those above are required to reproduce accurately the instabilities and many other aspects of currents throughout the world ocean. This is supported by snapshots of the surface features in two of the latest models. The warm Gulf Stream (Fig. 1A) is seen to separate at the proper location near Cape Hatteras and to follow its known path with pronounced eddy activity; lower resolution models had failed at these tasks. The simulated salinity near the southern tip of Africa (Fig. 1B) is filled with swirls and filaments and has a pronounced pair of counterrotating features in the Agulhas Current, which connects the Indian and Atlantic Oceans; these are further consequences of high resolution. Both panels in Fig. 1 look similar to the highest quality satellite images of ocean surface features.

Including unresolved effects. There are some processes that can never be included in basin-scale and global models: breaking surface waves are an example. Sometimes, net effects can be specified in terms of known quantities, in which case the method of doing so is called a parameterization. For breaking surface waves, the mixing they produce is probably quite amenable to a parameterization, as are other vertical mixing effects related to convection and small-scale internal waves.

Fig. 1. (A) Instantaneous ocean surface temperature from an Atlantic isopycnal model with a 1/12° grid. [Courtesy of R. Bleck] (B) Instantaneous ocean surface salinity [in parts per thousand (ppt)] from a global model with an average grid spacing of 1/6°. Click on the image for an expanded view

Whether mesoscale eddies that make up the bulk of the ocean's kinetic energy can be parameterized is much less certain. It will be difficult to quantify the small-scale correlations of important variables that provide for much of the north-south transport of properties. Most meteorologists have abandoned efforts to parameterize atmospheric storms and the effects that are described as negative viscosity; instead, atmospheric modelers simply resolve the storms explicitly. Parameterizing ocean eddies would be an attractive way to conduct long-term simulations at a relatively coarse ocean grid size of about 1°. A recently developed method of this type is being used in some simulations (28). Earlier model deficiencies in portraying the vertical extent of upper ocean structures have been remedied by this method.

Harnessing the Power of Supercomputers

With the typical explicit time stepping of almost all ocean models, one gains a factor of 8 in computational efficiency for each doubling of the horizontal grid size. (Fully implicit methods allow longer time steps, but only at the expense of iterative solution techniques that tend to squander any potential gains.) If there are no eddies, as in coarse-grid models, additional savings can come from a technique to accelerate convergence by artificially reducing the heat capacity of the deep levels. Without this acceleration technique, not even a 5° model could have run to the equilibrium of 1000 simulated years until the 1980s.

Vector machines. The advent of computers that rapidly perform the same operation on strings of numbers, so-called "vector computers" like the original CRAY-1A, was crucial to ocean modeling. By 1980, well-programmed ocean models were running a hundred times faster than they had a decade earlier. The memory size was beginning to be a barrier at that point. By the mid-1980s, the memory problem was ameliorated by the addition of secondary memory that had to be managed like disk storage but was much faster. Then, multiple processors were introduced, starting with the four-processor CRAY X-MP in 1986 and continuing with the eight-processor Y-MP in 1989. These were the first computers that were both fast enough and large enough in memory to allow global ocean modeling at a grid size that is approaching the radius of deformation. Ocean models account for some of the biggest applications on such machines, where they can run at speeds of a billion numerical operations per second.

Trillions of operations per second. A great hope among computational scientists is that a machine will be produced by the year 2000 that will run complex problems at a trillion operations per second. Such a teraflop (10^12 floating point operations per second) computer would almost certainly have hundreds or even thousands of processors, making it a massively parallel machine. Already, there are massively parallel machines that rival the vector machines with their multiple processors in attaining a present performance level of about 10 billion operations per second. Those are Thinking Machine's 1024-processor Connection Machine 5 and the CRAY T3D. Both the isopycnal and the GFDL models already mentioned have versions for these machines. Thus, adaptations have already been made to exploit teraflop machines for ocean simulations when those machines arrive.

Right now, a calculation at the radius of deformation, using a global grid of about 1/10° on average, would require about 3 days of the best machine's time per simulated year. Therefore, a teraflop computer could perform a 1000-year simulation to full oceanic equilibrium in a month. It is likely that such a calculation can and will be done early in the next century. When that happens, it will be the first time the turbulent ocean equilibrium, along with the associated 3D structure of its properties, is completely determined from simple physical principles. However, computer power already exists to portray the global ocean quite realistically over decades of simulated time. This means that a variety of important problems dealing with climatic, geochemical, and biological processes in the ocean are already tractable.

Global Data Sets

Forcing data. A major forcing field for all models, including adiabatic forecast models, is wind stress. Early modeling relied on ship anemometer records to build a climatological seasonal cycle of wind data. Present modeling relies on higher quality analyses of global meteorological data that are done twice daily at major weather forecasting centers. These fields give information about surface conditions over the last 15 years. Because weather analyses are made from sparse data over some mid-ocean regions, a better alternative that is slowly becoming available is satellite vector wind stress derived from radar backscatter.

Surface density changes are needed to drive properly the large overturning circulations that influence climate. This requires that surface exchanges of both heat and fresh water be known. Net heat exchange is available from crude estimates based on seasonally varying climatological data, but net water exchange resulting from evaporation minus precipitation is poorly known. Both of these fields will have to be specified from atmospheric analyses, with limited hope for improvements from satellite measurements in the near future.

In situ data. Historical observations over the past 100 years consist mainly of vertical profiles of temperature and salinity collected by oceanographic ships and maps of near-surface temperature from ships and weather satellites. They give an estimate of the long-term means and standard deviations of temperature, salinity, and density, but their quality is variable in different regions. Large international observing programs have been under way to provide much better data with which to test models. Among these are the Tropical Ocean and Global Atmosphere (TOGA) and the World Ocean Circulation Experiment (WOCE). In particular, WOCE is giving a modern description of the global ocean by extensive vertical profiling, deep moored current meters, surface drifting buoys, intermediate depth floats, and mapping from ships. TOGA provides information concentrated near the equator in the Pacific.

Surface height: A key variable. The ocean surface height is exceedingly important because its slopes yield the near-surface velocity and its time variations are measures of internal instabilities and wind-forced current fluctuations. It is observable by radar from space with little atmospheric obstruction. Surface height variations are now being measured from satellites such as TOPEX/POSEIDON and ERS-1 and 2. These are used as primary validation fields for ocean models. In addition, satellite altimetry can map the surface currents on a frequent basis, such as every 10 days for TOPEX and 35 days for the ERS satellites. The information can be used as a forcing function, relying on models to assimilate the surface height and progressively build reliable knowledge, even for features below the surface. This technique is still under development but promises to provide much better information about instantaneous ocean states. In that case, other observations could be more easily interpreted, and ocean forecast models could be properly initialized. Furthermore, a time sequence of such ocean states would constitute a useful data set of its own for understanding ocean physics.

Recent Highlights