## 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.

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}

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

which involves three separate contributions,

The first term, *w*_{H}, is the spurious diapycnal transfer due to horizontal diffusion of temperature—the “Veronis” effect.

The second term, *w*_{V}, 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} m^{2} s^{−1}, implying a diapycnal velocity *w*_{V} ∼ *K*_{θ}/*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

### 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}

where we have used the hydrostatic equation and assumed a thermal wind balance. Here *ζ*_{g} = (∂_{x}*υ*_{g} − ∂_{y}*u*_{g}) is the geostrophic relative vorticity, *ρ*_{0} is a typical density, *f* is the Coriolis parameter, *N*^{2} = ∂*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.

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} ∼ (*A*_{u}/*β*)^{1/5}, where *A*_{u} 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 *A*_{u} 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

where we have applied *A*_{θ} = 0 at the coastline (*x* = 0).

We assume that the biharmonic temperature diffusion is scaled such that

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

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

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

#### 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 |*w*_{H}| by assuming a mixing length scale *l,* giving a typical relative vorticity anomaly *ζ*′ ∼ *βl.* Thus,

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*A*_{u}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}°.

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

Here *p* is pressure, *ρ* is density, *g* is the acceleration due to gravity, and *K*_{u} 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:

The lateral mixing of momentum by subgrid turbulence **F**^{u} is represented in the model by biharmonic horizontal diffusion:

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

this form ensures that temperature variance is dissipated. Alternative formulations to horizontal diffusion are investigated later. The diffusion coefficients for momentum and temperature, *A*_{u}, *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,

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} m^{2} s^{−1}, broadly consistent with observations of diapycnal mixing.

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.

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.)

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 10^{14} m^{3} of water lighter than 26.0 over the integration, equivalent to a water mass transformation of approximately 0.3 Sv (Sv ≡ 10^{6} m^{3} 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.

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 gyre^{3}. 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.

### 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*_{θ}/*A*_{u} 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 *A*_{u}. 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.

## 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) by^{4}

where

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 m^{2} 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.

### 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:

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 m^{2} s^{−1} and *γ* = 2 × 10^{8} m^{4} 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.

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 × 10^{11} m^{4} 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 m^{2} 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.

## 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

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):

After two integrations by parts, (23) becomes

where we have assumed the boundary conditions **∇**[*γ*(∂*b*/∂*z*)^{−1}∇^{2}*b*]·**n** = **∇***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

where **∇**_{b} is the gradient operator evaluated along an isopycnal, *z*_{b} = ∂*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

_{2}. Part I: Annual mean response. J. Climate, 4, 785–818.

### 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,

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,

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 *S*_{max} = 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:

where the eddy transfer coefficient *γ* is again tapered to zero as the magnitude of the slope approaches *S*_{max}, 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.

Email: mjroberts@meto.gov.uk

^{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, *w*_{H}, 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, *w*_{H}, 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.