Abstract

Use of horizontal diffusion of temperature and salinity in numerical ocean models causes spurious diapycnal transfers—the “Veronis effect”—leading to erosion of the thermocline and reduced poleward heat transports. The authors derive a relation between these diapycnal transfers and the dissipation of vorticity gradients. An increase in model resolution does not significantly reduce the diapycnal transfers since vorticity gradients cascade to smaller scales and must ultimately be dissipated to maintain numerical stability. This is confirmed in an idealized primitive equation ocean model at a variety of resolutions between 1° and 1/8°.

Thus, the authors conclude that adiabatic dissipation schemes are required, even in eddy-resolving ocean models. The authors propose and implement a new biharmonic form of the Gent and McWilliams scheme, which adiabatically dissipates at the grid scale while preserving larger-scale features.

1. Introduction

Coarse-resolution ocean models used as part of coupled climate model simulations poorly maintain the sharp thermocline found within the upper kilometer of the ocean (e.g., Johns et al. 1997; Manabe et al. 1991). Spurious erosion of the modeled thermocline makes validation with observations very difficult—for example, the present trend of increasing temperatures at middepths, as observed by Parrilla et al. (1994) and Bryden et al. (1996), will be overwhelmed.

A recent analysis of North Atlantic Ocean simulations by Böning et al. (1995) suggests that erosion of the thermocline is associated with the “Veronis effect” (Veronis 1975), illustrated schematically in Fig. 1. Here, spurious horizontal diffusion of temperature across a boundary current leads to warming and cooling on its flanks, driving vertical motion as shown. The upwelling of deep waters into near-surface boundary currents reduces the transport of cold, deep water equatorward and consequently reduces poleward heat transports. The results of Böning et al. suggest that a simple increase in resolution, with a corresponding decrease in horizontal diffusion coefficients, only partially resolves this problem.

Fig. 1.

Schematic diagram illustrating the Veronis effect. Spurious horizontal diffusion of density across the front leads to warming onshore and cooling offshore, thereby inducing upwelling and downwelling, respectively. Here the dashed lines are isopycnals.

Fig. 1.

Schematic diagram illustrating the Veronis effect. Spurious horizontal diffusion of density across the front leads to warming onshore and cooling offshore, thereby inducing upwelling and downwelling, respectively. Here the dashed lines are isopycnals.

The aims of this paper are to show that

  • spurious diapycnal transfers associated with the Veronis effect are related to the dissipation of vorticity gradients in a model;

  • an increase in model resolution does not significantly reduce the diapycnal transfers since vorticity gradients cascade to smaller scales and must ultimately be dissipated to maintain numerical stability; and

  • consequently, adiabatic subgrid dissipation schemes are required, even in eddy-resolving models.

In section 2 a theoretical analysis of the Veronis effect is presented, explicitly relating the spurious diapycnal transfers to the dissipation of vorticity at the grid scale. Scale analysis suggests that the diapycnal transfers within the boundary currents are reduced as the grid spacing is reduced; however, at these higher resolutions additional diapycnal transfers occur in the interior of the basin, associated with the cascade of vorticity gradients to the grid scale. In section 3 the diapycnal transfers in an idealized ocean model are diagnosed, and it is shown that these are remarkably insensitive to model resolution between 1° and 1/8°. In section 4 various strategies for eliminating the spurious diapycnal transfers are explored; in particular a new scale-selective version of the Gent and McWilliams (1990) scheme is proposed that adiabatically dissipates grid-scale structures in an eddy-permitting model while leaving larger-scale features intact. We discuss the results from this work in section 5.

2. Theoretical background

a. Diapycnal transfer in numerical models

Consider the potential temperature equation from a typical eddy-resolving numerical ocean model,1

 
formula

which incorporates a biharmonic horizontal diffusion, a vertical diffusion, and external buoyancy forcing 𝒬⁠. Here the total derivative D/Dt = ∂/∂t + u· + w∂/∂z, where is the horizontal gradient operator, u is the horizontal velocity vector, w is the vertical velocity, θ is the potential temperature, Aθ is the biharmonic diffusion coefficient for temperature, and Kθ is the coefficient for vertical diffusion of temperature. The boundary conditions on the biharmonic diffusion are n· (Aθ2θ) = n·θ = 0 and ensure a global sink of tracer variance.

Following the approach of Walin (1982), we define the diapycnal velocity

 
formula

which involves three separate contributions,

 
formula

The first term, wH, is the spurious diapycnal transfer due to horizontal diffusion of temperature—the “Veronis” effect.

The second term, wV, parameterizes the diapycnal transfer associated with the breaking of internal waves and other diapycnal mixing processes. Results from tracer release experiments (e.g., Ledwell et al. 1993) suggest Kθ ∼ 10−5 m2 s−1, implying a diapycnal velocity wVKθ/H ∼ 10−8 m s−1, where H ∼ 1 km is a typical vertical scale.

The final term, w𝒬, represents the diapycnal transfer driven by external buoyancy forcing and is negligible beneath the surface mixed layer.

For a numerical model to maintain a realistic thermocline, a prerequisite for climate simulations, it is essential that the spurious first term is small relative to the remaining terms; that is, we require

 
wHwV, w𝒬.
(3)

b. Relation of diapycnal transfer to vorticity dissipation

Neglecting variations in the Coriolis parameter at the grid scale, the diapycnal transfer due to horizontal diffusion of temperature can be rewritten:2

 
formula

where we have used the hydrostatic equation and assumed a thermal wind balance. Here ζg = (∂xυg − ∂yug) is the geostrophic relative vorticity, ρ0 is a typical density, f is the Coriolis parameter, N2 = ∂b/∂z is the buoyancy frequency, and b is buoyancy. Thus, diapycnal transfer induced by the Veronis effect is intimately connected with the dissipation of vorticity gradients in a numerical ocean model. This result is significant for two reasons:

  • The sharpest vorticity gradients in numerical ocean models are maintained adjacent to coastlines by the no-slip (or free-slip) boundary condition. For example, in a northward flowing boundary current, low values of vorticity are advected northward within the inertial part of the boundary current, but large positive values of vorticity are injected at the coastline to satisfy the lateral boundary condition, as illustrated schematically in Fig. 2a. Hence, one expects the diapycnal transfers to be large within the no-slip sublayer adjacent to a western boundary. This is the classical Veronis effect identified in previous studies (e.g., Veronis 1975; Böning et al. 1995).

  • In models that resolve eddies, there is a turbulent cascade of vorticity gradients from large to small scales (Rhines 1977), as shown schematically in Fig. 2b. Thus, in an eddy-resolving model, we suggest that this turbulent cascade will lead to additional diapycnal transfers in the interior of an ocean basin, associated with the dissipation of vorticity filaments at the grid scale.

Fig. 2.

The Veronis effect as it relates to the dissipation of vorticity gradients in an ocean model. (a) In a western boundary current, the presence of large relative vorticity (shaded) within the no-slip sublayer leads to intense upwelling adjacent to the coastline. (b) The turbulent cascade leads to additional diapycnal transfers in the interior of an ocean basin, associated with the dissipation of vorticity filaments (shaded) at the grid scale.

Fig. 2.

The Veronis effect as it relates to the dissipation of vorticity gradients in an ocean model. (a) In a western boundary current, the presence of large relative vorticity (shaded) within the no-slip sublayer leads to intense upwelling adjacent to the coastline. (b) The turbulent cascade leads to additional diapycnal transfers in the interior of an ocean basin, associated with the dissipation of vorticity filaments (shaded) at the grid scale.

Hereafter in this paper we interpret the “Veronis effect” to include both boundary and interior contributions.

c. Scaling of diapycnal transfers with grid spacing

1) Boundary contribution

The largest vorticity gradients will be found within the no-slip sublayer. Classical homogeneous ocean circulation theory predicts the width of this layer to be δNS ∼ (Au/β)1/5, where Au is the coefficient of momentum dissipation (assumed here to be biharmonic) and β = df/dy is the gradient of the Coriolis parameter. In practice, ocean modelers choose Au such that δNS is indistinguishable from the grid spacing Δx, and thus here we assume δNS ∼ Δx.

The net diapycnal transfer within the no-slip layer can be obtained by integrating (6) between the coastline and the edge of the no-slip layer, giving

 
formula

where we have applied Aθ = 0 at the coastline (x = 0).

We assume that the biharmonic temperature diffusion is scaled such that

 
formula

where τθ is the dissipative timescale for thermal structures at the scale of the grid. Therefore, we find

 
formula

where H is a typical vertical scale height, and we have taken ζgυBC/(Δx) where υBC is the maximum boundary current velocity.

In an inertial boundary current, there is no reason to expect υBC to scale with the dissipation coefficients, and thus

 
formula

Conversely in a frictional boundary current, as would be obtained in a coarse-resolution model, one expects υBC ∝ (Δx)−1, giving

 
formula

2) Interior contribution

The diapycnal transfers induced by the vorticity cascade within the interior of the basin are more difficult to quantify. One can obtain an order of magnitude estimate of |wH| by assuming a mixing length scale l, giving a typical relative vorticity anomaly ζ′ ∼ βl. Thus,

 
formula

However the scaling is particularly sensitive to the mixing length scale, which is a highly nonlinear function of the grid spacing and the two dissipation coefficients. In particular, one would expect l to increase sharply as the grid spacing is decreased and the eddy field becomes more energetically active. In fully resolved turbulence, one expects l to tend toward the Rhines scale, (U/β)1/2 (Rhines 1975), where U is a typical background velocity.

It should be noted that these scalings only give an estimate for the root-mean-square diapycnal transfer, not its global integral. However, they do suggest three ways of reducing the bogus diapycnal motion in numerical ocean models:

  • increasing the resolution, which allows both Aθ and Au to be reduced;

  • increasing the horizontal viscosity to broaden the no-slip sublayer and reduce the eddy mixing length-scale;or

  • removing the horizontal diffusion of temperature altogether.

Each of these suggestions is now examined in a series of idealized numerical experiments.

3. Numerical experiments

a. Model parameters and forcing

We solve a Bryan–Cox primitive equation ocean model (Bryan 1969; Cox 1984) in a rectangular domain of size 1000 km by 2000 km, centered at 24°N, y = 0 (see Fig. 3). There is a flat bottom at depth 5500 m and 20 unevenly spaced levels in the vertical (see Table 1), with higher resolution near the surface to resolve the thermocline. The wind stress varies sinusoidally with a maximum magnitude of 0.1 N m−2 to give a double-gyre circulation. We consider horizontal resolutions of 1°, 1/2°, 1/4°, 1/6°, and 1/8°.

Fig. 3.

A schematic diagram of the model domain, which consists of a rectangular basin centred at 24°N (y = 0). A sinusoidal wind stress spins up subtropical and subpolar gyres as illustrated.

Fig. 3.

A schematic diagram of the model domain, which consists of a rectangular basin centred at 24°N (y = 0). A sinusoidal wind stress spins up subtropical and subpolar gyres as illustrated.

Table 1.

Table of the depths of each model level and the thickness of each of the vertical grid cells.

Table of the depths of each model level and the thickness of each of the vertical grid cells.
Table of the depths of each model level and the thickness of each of the vertical grid cells.

Our strategy is to integrate for 10 years from a state of rest and an initial exponential temperature profile. The e-folding scale is 1000 m, and the temperature ranges from 22°C at the surface to 2.2°C at the bottom. We quantify the extent to which the initial temperature profile is retained over the 10 years.

The model equations are

 
formula

Here p is pressure, ρ is density, g is the acceleration due to gravity, and Ku is the coefficient for vertical diffusion of momentum. We adopt the UNESCO equation of state (UNESCO 1981), with salinity held constant at 35.0 psu. There is no thermal forcing. The following values are employed:

 
g = 9.81 m s−2,  Ku = 1.0 × 10−3 m2 s−1,ρ0 = 103 kg m−3.

The lateral mixing of momentum by subgrid turbulence Fu is represented in the model by biharmonic horizontal diffusion:

 
Fu = −∇2(Au2u).
(17)

The lateral turbulent mixing of temperature is initially parameterized by a biharmonic diffusion:

 
Fθ = −∇2(Aθ2θ);
(18)

this form ensures that temperature variance is dissipated. Alternative formulations to horizontal diffusion are investigated later. The diffusion coefficients for momentum and temperature, Au, Aθ respectively, are varied: details are given with each experiment.

b. Sensitivity of Veronis effect to horizontal resolution

First, we consider a series of experiments in which both biharmonic diffusion coefficients are scaled roughly with the fourth power of the grid spacing,

 
formula

where τθ = 24 days and τu = 95 days. The precise values are given in Table 2. In these integrations, the vertical diffusion remains constant at 3 × 10−5 m2 s−1, broadly consistent with observations of diapycnal mixing.

Table 2.

Biharmonic diffusion coefficients, Aθ and Au, for experiments detailed in section 3b. All experiments have Kθ = 3 × 10−5 m2 s−1.

Biharmonic diffusion coefficients, Aθ and Au, for experiments detailed in section 3b. All experiments have Kθ = 3 × 10−5 m2 s−1.
Biharmonic diffusion coefficients, Aθ and Au, for experiments detailed in section 3b. All experiments have Kθ = 3 × 10−5 m2 s−1.

The time-mean (years 6–10) streamfunction for the depth-integrated circulation from the various resolution models is shown in Fig. 4. As the resolution increases, inertial effects become more important and alter the basic Sverdrup gyres; tight recirculating subgyres appear to the north and south of the separated inertial jet. The emergence of a vigorous geostrophic eddy field at the highest resolutions curtails the penetration of the time-mean separated jet.

Fig. 4.

Barotropic streamfunction (Sv) averaged between years 6 and 10 for resolutions (a) 1°, (b) 1/2°, (c) 1/4°, (d) 1/6°, and (e) 1/8°.

Fig. 4.

Barotropic streamfunction (Sv) averaged between years 6 and 10 for resolutions (a) 1°, (b) 1/2°, (c) 1/4°, (d) 1/6°, and (e) 1/8°.

The water mass budget, that is, the volume of water contained between any two isopycnals, can only be changed through horizontal diffusion, vertical diffusion, and implicit numerical diffusion since there is no surface buoyancy forcing in the model. The integrated effect of these processes is quantified in Fig. 5, which plots the mean depth of each isopycnal, z(ρ), at years 0 and 10 for the various model resolutions. The lightest water masses are lost in each case, their heat is diffused downward and warms the main thermocline (defined to be approximately 300–600 m) by approximately 1°. In contrast to the findings of Böning et al. (1995), there is remarkably little sensitivity to an increase in resolution. The contribution of vertical diffusion to these water mass budgets is quantified in Fig. 6, which shows z(ρ) after 10 years in a 1° integration with no vertical diffusion, but with all other parameters unchanged from the previous experiment. The loss of light water masses is only slightly reduced, ruling out vertical diffusion as a major contributor to the diapycnal transfers. One might expect implicit numerical diffusion to become important at the higher resolutions as vertical velocities increase and grid Peclet numbers become larger. To exclude this, we have diagnosed the water mass budgets from a 1/6° integration with an enhanced vertical resolution of 30 levels;the differences are found to be very small. (Further evidence that horizontal diffusion is the major contribution to the diapycnal transfers is obtained from the adiabatic integrations described later: these are subject to the same truncation errors, but show little evidence of implicit numerical diffusion.)

Fig. 5.

Area-averaged depth of density surfaces z(ρ) after 10 yr of integration at various model resolutions, as detailed in section 3.b/Table 2. Only the upper 900 m is shown, and the initial conditions are shown by the solid line.

Fig. 5.

Area-averaged depth of density surfaces z(ρ) after 10 yr of integration at various model resolutions, as detailed in section 3.b/Table 2. Only the upper 900 m is shown, and the initial conditions are shown by the solid line.

Fig. 6.

Area-averaged depth of density surfaces for 1° model integrations with different vertical diffusivities after 10 yr.

Fig. 6.

Area-averaged depth of density surfaces for 1° model integrations with different vertical diffusivities after 10 yr.

To estimate the magnitude of the water mass transformation induced by the spurious horizontal diffusion we examine the change in the mean depth of an isopycnal over the 10 years in Fig. 5. For example, the mean depth of the 26.0 isopycnal (located within the main thermocline) changes by approximately 50 m over the 10 years. This implies a loss of 1014 m3 of water lighter than 26.0 over the integration, equivalent to a water mass transformation of approximately 0.3 Sv (Sv ≡ 106 m3 s−1). Scaled up to the global ocean this gives a transformation of 50 Sv, which far exceeds estimates of diapycnal upwelling inferred from observations.

To explain these results a series of instantaneous relative vorticity plots is shown in Fig. 7 at the various model resolutions. As shown in section 2b the diapycnal transfers associated with the Veronis effect are related to relative vorticity gradients through (6). At 1° and 1/2°, there are no eddies, and hence the largest vorticity gradients are within the frictional boundary current. As the resolution is increased and eddies form, we see the typical signature of geostrophic turbulence, which is a cascade of vorticity gradients to the smallest scales. In particular, in the 1/8° integration the gyre interior is filled with vorticity filaments; the vorticity gradients within the no-slip sublayer are also markedly increased.

Fig. 7.

Instantaneous relative vorticity fields at 300 m after 10 yr for different model resolutions of (a) 1°, (b) 1/2°, (c) 1/4°, (d) 1/6°, and (e) 1/8°. At higher resolutions there is a cascade of vorticity filaments to the grid scale. Shading indicates positive values.

Fig. 7.

Instantaneous relative vorticity fields at 300 m after 10 yr for different model resolutions of (a) 1°, (b) 1/2°, (c) 1/4°, (d) 1/6°, and (e) 1/8°. At higher resolutions there is a cascade of vorticity filaments to the grid scale. Shading indicates positive values.

These results are at odds with the scaling arguments presented in section 2c, which suggest that the diapycnal transfers should be reduced at higher resolution, and with the findings of Böning et al. (1995). To explain this paradox, in Fig. 8 we plot the local upwelling within the western boundary current of the subtropical gyre3. As in Böning et al., we find a decrease in boundary layer upwelling at higher resolutions. However, this does not imply a decrease in the global diapycnal transfer for three reasons: (i) the local upwelling is not equivalent to the diapycnal transfer; (ii) the diapycnal transfers within the boundary currents of the subtropical and subpolar gyres are of opposite sign, and it is unclear how their residual should vary with the grid spacing; and (iii) the formation of vorticity filaments at the grid scale leads to enhanced interior diapycnal transfers at higher resolutions.

Fig. 8.

Net upwelling (Sv) through the first grid point in the western boundary layer of the subtropical gyre at 360 m, as a function of model resolution. The upwelling is reduced at high resolutions consistent with the findings of Böning et al. (1995). The upwelling is also reduced in the lowest 1° resolution, due to the transition between frictional and inertial boundary current regimes, and associated changes in the boundary current transport and extension.

Fig. 8.

Net upwelling (Sv) through the first grid point in the western boundary layer of the subtropical gyre at 360 m, as a function of model resolution. The upwelling is reduced at high resolutions consistent with the findings of Böning et al. (1995). The upwelling is also reduced in the lowest 1° resolution, due to the transition between frictional and inertial boundary current regimes, and associated changes in the boundary current transport and extension.

c. Sensitivity to diffusion ratios

One method for reducing the diapycnal transfer might be to increase the magnitude of the momentum diffusion relative to the temperature diffusion (see, e.g., Böning et al. 1995). This will both broaden the no-slip sublayer and reduce the interior vorticity gradients. Both will act to decrease the diapycnal transfer given by (6). Figure 9 shows three 1/4° integrations with different ratios of Aθ/Au as shown in Table 3. Only a very slight improvement in the retention of light water masses is achieved by increasing the coefficient of momentum diffusion Au. Curiously, the impact of reducing the coefficient of temperature diffusion is also only slight. It would appear that the only way to remove the spurious diapycnal transfers is to remove the horizontal diffusion of temperature altogether.

Fig. 9.

Area-averaged depth of density surfaces for model resolutions with different ratios of Aθ/Au after 10 yr (see Table 3).

Fig. 9.

Area-averaged depth of density surfaces for model resolutions with different ratios of Aθ/Au after 10 yr (see Table 3).

Table 3.

Biharmonic diffusion coefficients for experimentsdetailed in section 3c. All experiments have Kθ = 3 × 10−5 m2 s−1.

Biharmonic diffusion coefficients for experimentsdetailed in section 3c. All experiments have Kθ = 3 × 10−5 m2 s−1.
Biharmonic diffusion coefficients for experimentsdetailed in section 3c. All experiments have Kθ = 3 × 10−5 m2 s−1.

4. Adiabatic subgrid parameterization

The incorporation of the new adiabatic eddy parameterization scheme devised by Gent and McWilliams (1990, hereafter GM90) into coarse-resolution ocean models has lead to dramatic improvements in model behavior, for example, a sharper thermocline, improved heat transports, and a restriction of deep convection to places where it is known to occur (Danabasoglu et al. 1994). Most of these improvements can be attributed to the complete removal of horizontal diffusion of temperature and salinity, which eliminates the Veronis effect. Our results suggest that an adiabatic subgrid parameterization will be necessary, even in models that resolve or partially resolve mesoscale eddies.

One possibility might be to use a shape-preserving advection scheme (e.g., the flux-corrected transport algorithm; Boris and Book 1973), which would enable explicit temperature diffusion to be excluded. However, Thuburn (1995) has shown that the implicit numerical diffusion acting on a grid-scale structure as it is advected across a grid cell by a shape-preserving algorithm is finite and comparable to that obtained using standard biharmonic diffusion. Thus, it is not apparent that such schemes will remove the Veronis effect to a satisfactory degree.

In this section, we investigate the extent to which the retention of water masses can be improved by implementing the adiabatic GM90 scheme in our idealized numerical model. We then propose a new scale-selective variation of the scheme for use in eddy-permitting models.

a. Gent and McWilliams scheme

Gent and McWilliams (1990) parameterize the lateral turbulent mixing term Fθ in (16) by4

 
formula

where

 
formula

and b is buoyancy.

The effectiveness of the scheme in preserving the model’s water masses is shown in Fig. 10. Here we plot the mean depth of the isopycnals after 10 years in runs employing GM90 at 1° and 1/4° resolution. For comparison the equivalent results from integrations employing biharmonic horizontal diffusion are also shown. The transfer coefficients κ used in the GM90 runs are 400 and 100 m2 s−1 at 1° and 1/4° respectively; the numerical discretization is given in the appendix. There is no vertical diffusion in each of these experiments, and thus the initial water mass profile should be exactly preserved in a perfect model. We find that the use of GM90 greatly reduces the loss of light water masses from the model. The temperature of the main thermocline is also much better preserved using GM90 at 1° resolution, but the thermocline cools at 1/4° resolution due to vertical motions that develop at the grid scale around the lateral boundaries.

Fig. 10.

Area-averaged depth of density surfaces for model integrations with biharmonic horizontal diffusion or GM90 schemes after 10 yr. These integrations have no vertical diffusion.

Fig. 10.

Area-averaged depth of density surfaces for model integrations with biharmonic horizontal diffusion or GM90 schemes after 10 yr. These integrations have no vertical diffusion.

b. Scale-selective parameterization

The GM90 scheme has been developed primarily as a parameterization of mesoscale eddies for use in coarse-resolution ocean models. While GM90 has been used in eddy-resolving models (Haines and Wu 1998, submitted to J. Mar. Sys.), it is not ideally suited to this purpose: GM90 is less scale selective than biharmonic diffusion and strongly damps eddies and fronts in addition to grid-scale structures.

The GM90 scheme is based on two key assumptions:(i) eddies flux isopycnal layer thickness downgradient and (ii) the eddies are dissipated adiabatically, that is, the eddy energy is dissipated by friction rather than by diapycnal wave breaking (Tandon and Garrett 1996). In an eddy-resolving model, one hopes to partially resolve downgradient transfer; however, the cascade of vorticity, potential enstrophy, and other tracerlike quantities to smaller scales means that one must always parameterize the dissipation at the grid scale.

Hence, there is a need for a scale-selective form of GM90 that will remove enstrophy adiabatically at the grid scale. The simplest form for such a scheme, as a development of GM90, proposes that the velocities of (21) are redefined:

 
formula

Just as biharmonic diffusion is more scale selective than Laplacian diffusion, so in this case the ∇2 operator acting on the slope of the isopycnals ensures that our scheme acts preferentially at the grid scale while leaving larger-scale features intact.

The effectiveness of the biharmonic scheme can be demonstrated in a simple two-dimensional model. Starting with a simple front upon which grid-scale noise is superimposed, we integrate for 3600 days with the GM90 and biharmonic GM schemes (Fig. 11). Here the domain width is 200 km, and the grid spacing is 5 km;we choose the values κ = 8 m2 s−1 and γ = 2 × 108 m4 s−1, both of which give dissipative timescales of about a month for thermal structures at the scale of the grid. After 5 days, the biharmonic GM scheme has removed the grid-scale noise, while it remains in the GM90 case. After 360 days, the grid-scale noise has been removed from both runs; however, the standard GM90 scheme also slumps the large-scale front, while this is relatively unchanged with the biharmonic GM scheme. This is emphasized after 3600 days, in which the front has virtually disappeared in the GM90 case, but remains largely unaltered in the biharmonic GM case.

Fig. 11.

Results from a two-dimensional calculation initialized with a density front upon which grid-scale noise has been superimposed. The front is shown after 5, 360, and 3600 days of integration using the standard GM90 scheme (left column) and the new biharmonic GM scheme (right column). The streamfunctions for the (u*, w*) velocity are superimposed on the day 360 front (note different contour intervals are used for each: the biharmonic cell is an order of magnitude weaker).

Fig. 11.

Results from a two-dimensional calculation initialized with a density front upon which grid-scale noise has been superimposed. The front is shown after 5, 360, and 3600 days of integration using the standard GM90 scheme (left column) and the new biharmonic GM scheme (right column). The streamfunctions for the (u*, w*) velocity are superimposed on the day 360 front (note different contour intervals are used for each: the biharmonic cell is an order of magnitude weaker).

We have implemented the new biharmonic adiabatic scheme in our idealized primitive equation ocean model at 1/6° resolution with no vertical diffusion (γ = 1.6 × 1011 m4 s−1); the numerical discretization is again given in the appendix. In Fig. 12c, we show the instantaneous relative vorticity after 10 years of integration, together with similar fields from 1/6° integrations using biharmonic horizontal diffusion and GM90 (κ = 40 m2 s−1) (Figs. 12a,b). The biharmonic scheme removes the grid-scale noise and also up/downwelling cells found around the boundary which are apparent in the standard GM90 scheme. At the same time, the biharmonic scheme allows the vorticity field to cascade and form filaments to at least the same extent as the horizontal diffusion case. In contrast, the cascade is prematurely curtailed in the GM90 integration, even with such a small value of κ; a larger value would reduce the grid-scale noise but, as shown in the two-dimensional model, it would also greatly damp resolved features including boundary currents. Figure 13 shows how effective the new scheme is in preserving the model’s water masses at 1/6° resolution; the 1/4° GM90 integration is also plotted for comparison. The new scheme greatly improves water mass retention compared to biharmonic horizontal diffusion (compare with Fig. 5) and also preserves the temperature of the main thermocline compared to the standard GM90 scheme (due mainly to the removal of the spurious grid-scale boundary cells). However, there is a greater loss of light near-surface water mass in the biharmonic GM scheme than in GM90.

Fig. 12.

Instantaneous relative vorticity at 300 m after 10 yr at 1/6° resolution for (a) biharmonic diffusion, (b) standard GM90 scheme, and (c) biharmonic GM scheme.

Fig. 12.

Instantaneous relative vorticity at 300 m after 10 yr at 1/6° resolution for (a) biharmonic diffusion, (b) standard GM90 scheme, and (c) biharmonic GM scheme.

Fig. 13.

Area-averaged depth of density surfaces after 10 yr for a 1/6° model integration with biharmonic GM. To ease comparison with Fig. 10, we also include the curve from the 1/4° integration with GM90. These integrations have no vertical diffusion.

Fig. 13.

Area-averaged depth of density surfaces after 10 yr for a 1/6° model integration with biharmonic GM. To ease comparison with Fig. 10, we also include the curve from the 1/4° integration with GM90. These integrations have no vertical diffusion.

5. Discussion

The spurious diapycnal transfers inherent in most coarse-resolution ocean models cause gross errors in the representation of water masses, thermocline gradients, and northward heat transports. Most of the spurious diapycnal transfers occur adjacent to boundary currents as a consequence of the horizontal diffusion of temperature and salinity (the Veronis effect). Previously, it has been assumed that as resolution is increased and diffusion coefficients tuned downward, the Veronis effect becomes less significant.

In this paper we have shown that the diapycnal transfers associated with the Veronis effect are related to the dissipation of vorticity gradients. Increasing the model resolution allows these gradients to cascade to smaller scales, where additional diapycnal transfers occur across vorticity filaments. Our idealized experiments suggest that the magnitude of the Veronis effect is surprisingly insensitive to the model resolution when the dissipative coefficients are scaled in the conventional manner. In particular, buoyant surface water masses are lost from the model, and there is a compensating decrease in density within the main thermocline. The former may be masked by air–sea transfer and/or surface mixing in a“realistic” model, but the fluxes of heat across the air–sea interface and across the base of the mixed layer will nevertheless be in error. Thus, we conclude that it is necessary to use an adiabatic dissipation scheme, even in an eddy-resolving model.

Implementation of the recent Gent and McWilliams scheme improves the retention of water masses in our model. However, for sensible choices of transfer coefficient, the scheme appears to generate excessive grid-scale noise while prematurely curtailing the turbulent cascade. We therefore propose a new biharmonic version of GM90 for use in eddy-permitting models; it acts strongly at the grid scale, arresting the turbulent cascade, but weakly at larger, resolved scales. This scheme dramatically improves the retention of water masses in an eddy-permitting model, while providing all the benefits of a traditional biharmonic diffusion.

There are a number of remaining issues that we briefly mention here.

a. Potential energy considerations

An attractive property of the GM90 scheme is that it provides a domain-averaged sink of potential energy. The biharmonic scheme proposed in (22) is convenient in computational terms; however, it does not ensure a domain-averaged sink of potential energy.

The equation for the rate of change of potential energy ρgz due to the biharmonic GM scheme is

 
formula

where we have integrated by parts, and u*3 and 3 denote three-dimensional vectors. To ensure that the new scheme provides a sink of potential energy requires that the rhs of (23) be negative. This can be achieved with a modified formulation (P. Gent 1997, personal communication):

 
formula

After two integrations by parts, (23) becomes

 
formula

where we have assumed the boundary conditions [γ(∂b/∂z)−12bn = b·n = 0 on solid boundaries. However, we have yet to find a numerical discretization that preserves this integral property and which is computationally stable. In contrast, the scheme employed in section 4 appears remarkably robust to date.

b. General implementation for tracers

In this paper we have focused on the idealized limit in which salinity is uniform, and the adiabatic scheme can therefore be expressed entirely in terms of an advective flux. However, more generally, it is necessary to incorporate both an advective and diffusive flux of tracer along isopycnals (Gent and McWilliams 1990). In the full implementation of the new biharmonic scheme, (20) should be replaced by

 
formula

where b is the gradient operator evaluated along an isopycnal, zb = ∂z/∂b, and C is an arbitrary tracer. The form of the biharmonic diffusion term ensures a global sink of tracer variance, although other forms are possible. However, in practice, implementation of any biharmonic isopycnal diffusion operator is likely to be problematic and will require further investigation.

c. Flux-corrected transport schemes

There may be alternative approaches for eliminating the Veronis effect. For example, one can reasonably argue for a shape-preserving advection algorithm that requires no explicit horizontal diffusion to maintain numerical stability. However, it remains to be shown that the implicit numerical diffusion acting on grid-scale structures is sufficiently small to damp the Veronis effect to a satisfactory level (see, e.g., Thuburn 1995). Moreover, the computational overhead of shape-preserving schemes can be prohibitive for use in long-range climate integrations. In principle, it should be possible to formulate our biharmonic scheme as part of a skew-symmetric diffusion (Griffies 1998), in which case its overhead should be minimal.

d. Intermediate resolutions

We regard our biharmonic scheme primarily as a numerical device to adiabatically truncate the turbulent cascade and thereby avoid the spurious diapycnal transfers associated with the Veronis effect. Our results suggest that such schemes are required even in models that resolve the mesoscale eddy field.

However at intermediate resolutions, in which the mesoscale is marginally resolved, there may be a role for using the biharmonic scheme in conjunction with the standard Gent and McWilliams scheme [with variable transfer coefficients (Visbeck et al. 1997)]; the former to dissipate at the grid scale, and the latter to parameterize subgrid-scale eddies. For example, in the Hadley Centre model run at 1.25° resolution, an artificially large eddy transfer coefficient is required to damp grid-scale noise using a GM90-type scheme; this has an adverse impact on the structure of boundary currents. Preliminary experiments suggest that the grid-scale noise is efficiently damped if the biharmonic scheme is introduced in parallel.

Acknowledgments

We are grateful to Peter Gent and an anonymous reviewer for their constructive suggestions; these have lead to a significantly improved manuscript. We also appreciate useful communications with John Thuburn, Stephen Griffies, and Brian Hoskins. This work has been carried out under DOE Contract PECD/7/12/37.

REFERENCES

REFERENCES
Böning, C. W., W. R. Holland, F. O. Bryan, G. Danabasoglu, and J. C. McWilliams, 1995: An overlooked problem in model simulations of the thermohaline circulation and heat transport in the Atlantic Ocean. J. Climate, 8, 515–523.
Boris, J. P., and D. L. Book, 1973: Flux-corrected transport, I. SHASTA: A fluid transport algorithm that works. J. Comput. Phys., 11, 38–69.
Bryan, K., 1969: A numerical method for the study of the circulation of the world ocean. J. Comput. Phys., 4, 347–376.
Bryden, H. L., M. J. Griffiths, A. L. Lavín, R. C. Millard, G. Parrilla, and W. M. Smethie, 1996: Decadal changes in water mass characteristics at 24°N in the subtropical North Atlantic Ocean. J. Climate, 9, 3162–3186.
Cox, M. D., 1984: A primitive equation, three dimensional model of the ocean. Ocean Group Tech. Rep. 1, Geophysical Fluid Dynamics Laboratory, Princeton, NJ, 143 pp.
Danabasoglu, G., and J. C. McWilliams, 1995: Sensitivity of the global ocean circulation to parameterizations of mesoscale tracer transports. J. Climate, 8, 2967–2987.
——, ——, and P. R. Gent, 1994: The role of mesoscale tracer transports in the global ocean circulation. Science, 264, 1123–1126.
Gent, P. R., and J. C. McWilliams, 1990: Isopycnal mixing in ocean circulation models. J. Phys. Oceanogr., 20, 150–155.
Griffies, S. M., 1998: The Gent–McWilliams skew-flux. J. Phys. Oceanogr., 28, 831–841.
——, A. Gnanadesikan, R. C. Pacanowski, V. D. Larichev, J. K. Dukowicz, and R. D. Smith, 1998: Isoneutral diffusion in a z-coordinate ocean model. J. Phys. Oceanogr., 28, 805–830.
Haines, K., and P. Wu, 1998: GCM studies of intermediate and deep waters in the Mediterranean. J. Mar. Sys., in press.
Johns, T. C., R. E. Carnell, J. F. Crossley, J. M. Gregory, J. F. B. Mitchell, C. A. Senior, S. F. B. Tett, and R. A. Wood, 1997: The second Hadley Centre coupled ocean–atmosphere GCM: Model description, spinup and validation. Climate Dyn., 13, 103–134.
Ledwell, J. R., A. J. Watson, and C. S. Law, 1993: Evidence for slow mixing across the pycnocline from an open-ocean tracer-release experiment. Nature, 361, 701–703.
Manabe, S., R. J. Stouffer, M. J. Spelman, and K. Bryan, 1991: Transient responses of a coupled ocean–atmosphere model to gradual changes of atmospheric CO2. Part I: Annual mean response. J. Climate, 4, 785–818.
Parrilla, G., A. Lavín, H. Bryden, M. García, and R. Millard, 1994:Rising temperatures in the subtropical North Atlantic Ocean over the past 35 years. Nature, 369, 48–51.
Rhines, P. B., 1975: Waves and turbulence on a β-plane. J. Fluid Mech., 69, 417–443.
——, 1977: The dynamics of unsteady currents. The Sea, Vol. VI. Wiley, 189–318.
Tandon, A., and C. Garrett, 1996: On a recent parameterization of mesoscale eddies. J. Phys. Oceanogr., 26, 406–411.
Thuburn, J., 1995: Dissipation and cascades to small scales in numerical models using a shape-preserving advection scheme. Mon. Wea. Rev., 123, 1888–1903.
UNESCO, 1981: Tenth report of the joint panel on oceanographic tables and standards. UNESCO Tech. Papers in Marine Sci. No. 36, UNESCO, Paris.
Veronis, G., 1975: The role of models in tracer studies. Numerical Models of the Ocean Circulation, Natl. Acad. of Sci., 133–146.
Visbeck, M., J. Marshall, T. Haine, and M. Spall, 1997: On the specification of eddy transfer coefficients in coarse-resolution ocean circulation models. J. Phys. Oceanogr., 27, 381–402.
Walin, G., 1982: On the relation between sea-surface heat flow and thermal circulation in the ocean. Tellus, 34, 187–195.

APPENDIX

Numerical Implementation of Adiabatic Schemes

The dissipative properties of the biharmonic Gent and McWilliams (1990) scheme rely on a robust numerical implementation. Here we provide details of the numerical discretizations employed in section 4.

GM90 scheme

The GM90 scheme is discretized on a C-grid,

 
formula

where superscripts indicate an average over adjacent gridpoints and subscripts denote centered differentiation. This discretization avoids simultaneous applications of averaging and differentiation operators in the same spatial direction and thereby prevents computational modes from developing [see Griffies et al. (1998) for a related discussion]. The vertical component is obtained from integration of the continuity equation on a C-grid,

 
formula

and the eddy transfer coefficient κ is set to zero at the ocean top and bottom so that w* is zero here; κ is also tapered to zero as the magnitude of the slope approaches a value Smax = 0.01; this is accomplished using a hyperbolic tangent function, as described in Danabasoglu and McWilliams (1995).

Biharmonic scheme

The discretization of the new biharmonic adiabatic scheme is a straightforward extension of the above:

 
formula

where the eddy transfer coefficient γ is again tapered to zero as the magnitude of the slope approaches Smax, and w* is evaluated using the continuity equation.

Footnotes

Corresponding author address: Dr. Malcolm J. Roberts, Hadley Centre for Climate Prediction and Research, Meteorological Office, Room H022, London Road, Bracknell, Berkshire RG12 2SY, United Kingdom.

1

In this paper, we consider an ocean of uniform salinity and define the diapycnal transfers in terms of potential temperature θ. In general, however, the diapycnal transfers should be defined in terms of buoyancy b.

2

The relation between the diapycnal transfer, wH, and dissipation of vorticity is a general result, although the details of this relation are dependent on the form of horizontal temperature diffusion.

3

This is a proxy for the diapycnal transfer, wH, which is difficult to diagnose accurately.

4

In general, one should also include a diffusion of temperature, salinity, and other water mass properties along density surfaces. In our model, however, this term is identically zero due to the use of a constant salinity.