## Abstract

A simplified version of the near-boundary eddy flux parameterization developed recently by Ferrari et al. has been implemented in the NCAR Community Climate System Model (CCSM3) ocean component for the surface boundary only. This scheme includes the effects of diabatic mesoscale fluxes within the surface layer. The experiments with the new parameterization show significant improvements compared to a control integration that tapers the effects of the eddies as the surface is approached. Such surface tapering is typical of present implementations of eddy transport in some current ocean models. The comparison is also promising versus available observations and results from an eddy-resolving model. These improvements include the elimination of strong, near-surface, eddy-induced circulations and a better heat transport profile in the upper ocean. The experiments with the new scheme also show reduced abyssal cooling and diminished trends in the potential temperature drifts. Furthermore, the need for any ad hoc, near-surface taper functions is eliminated. The impact of the new parameterization is mostly associated with the modified eddy-induced velocity treatment near the surface. The new parameterization acts in the depth range exposed to enhanced turbulent mixing at the ocean surface. This depth range includes the actively turbulent boundary layer and a transition layer underneath, composed of waters intermittently exposed to mixing. The mixed layer, that is, the regions of weak stratification at the ocean surface, is found to be a good proxy for the sum of the boundary layer depth and transition layer thickness.

## 1. Introduction

Ocean general circulation models (OGCMs) routinely used in long climate simulations cannot afford to explicitly resolve the geostrophic eddies at the mesoscale. These eddies contain most of the kinetic energy of the ocean and need to be parameterized to include their effects on the ocean circulation and climate. There is also a demonstrated need for adiabatic eddy parameterizations in eddy-permitting (Roberts and Marshall 1998) and even in eddy-resolving (Smith and Gent 2004a) models because increased model resolution does not necessarily resolve the full mesoscale spectrum. An eddy parameterization that has been extensively used during the last decade is the Gent and McWilliams (1990, hereafter GM90) isopycnal transport parameterization where the tracers are advected by an additional, eddy-induced velocity and diffused along isopycnals.

In the upper ocean within the surface mixed layer, the mesoscale eddy fluxes remain important and may even dominate the tracer budgets, as documented in various observational, theoretical, and modeling studies (e.g., Schmitz 1996; Treguier et al. 1997; Phillips and Rintoul 2000; Weller et al. 2004; Marshall 2005; Kuo et al. 2005). Here, the eddy velocities are constrained to follow the boundary, and the isopycnals often outcrop, thus producing diabatic fluxes. The GM90 eddy parameterization is developed for the quasi-adiabatic interior, and it is not valid near the ocean boundaries where the diabatic processes are significant.

In GM90, the eddy-induced velocity is directly proportional to the eddy diffusivity and isopycnal slope and inversely proportional to the isopycnal layer thickness. Therefore, as the ocean surface is approached and isopycnal slopes tend to get large and layer thicknesses tend to get small, the GM90 scheme yields unrealistically large eddy-induced velocities. To our knowledge, there is no observational evidence for such strong eddy-induced circulations at the boundaries: eddy fluxes have low-mode vertical structure and do not show signs of surface trapped circulations. So, both to reduce these spurious circulations and yet remain faithful to the quasi-adiabatic constraint, various tapering functions have been introduced to attenuate the strength of the eddy flux, particularly near the ocean surface. For example, Large et al. (1997) applies a taper based on preferred eddy length-scale arguments. For the same purpose, Hirst and McDougall (1996) and Griffies et al. (2005) use a maximum slope of 0.002, either globally or locally only in the upper ocean. These and any other tapering approaches (e.g., Gerdes et al. 1991) very effectively and quite strongly reduce both isopycnal diffusion and eddy-induced velocities near the ocean surface. As discussed in Griffies et al. (2005), Gnanadesikan et al. (2007), and Ferrari et al. 2008 hereafter FMCD), the OGCM results are sensitive to the choices of the maximum slope value and the tapering functions. In addition, the zero eddy flux boundary conditions at the ocean surface and bottom that produce the vanishing normal component of the eddy-induced velocity are generally enforced by setting the diffusivities to zero at these boundaries. In a widely used discrete implementation of the scheme (Griffies et al. 1998), a tracer grid box is subdivided into a top and bottom half in the vertical, with each half possibly carrying a different value for the diffusivity coefficient. So to impose, for example, the surface boundary condition, the diffusivities for the top halves of all surface tracer grid boxes are set to zero, thus completely turning off any mixing there. This approach forces all eddy fluxes to vanish—not just the normal component—and results in no eddy mixing within the half cells. These reductions in eddy-induced velocity and eddy flux within the surface mixed layer are clearly contrary to observational evidence of large surface eddy fluxes (e.g., Robbins et al. 2000; Price 2001).

The treatment of boundary conditions, particularly near the ocean surface, has been the focus of many recent studies (e.g., Tandon and Garrett 1996; Treguier et al. 1997; McDougall and McIntosh 2001; Killworth 2005). Among the suggested modifications is the addition of horizontal diffusion within the surface layer to represent diapycnal mixing, complementing reduced isopycnal diffusivity (Treguier et al. 1997). Smith and Gent (2004a) follow this suggestion and report much reduced numerical noise in the upper ocean in their eddy-resolving model. Griffies et al. (2005) further extend this approach and add both horizontal diffusion and an eddy-induced velocity with zero shear wherever isopycnal diffusivities are reduced due to tapering.

FMCD shows that the diabatic nature of the eddy fluxes can be retained with some simple modifications of the extant schemes near the boundaries and proposes a new eddy parameterization for these near-boundary regions. In the turbulent boundary layer (BL) the eddy-induced velocity is set parallel to the boundary and has no vertical shear, as expected in a well mixed layer. In addition there is an eddy diffusion of buoyancy along the boundary as well as along isopycnals. In the interior the parameterization satisfies the adiabatic constraint as, for example, in GM90. The two forms are matched through a transition layer that separates the quasi-adiabatic interior with isopycnally oriented eddy fluxes from the near-boundary regions. Treguier et al. (1997) have argued for a similar approach but they did not derive an explicit parameterization. Griffies et al. (2005) made a first attempt at implementing the ideas of Treguier et al. (1997). FMCD revisited the problem and provided a careful theoretical derivation of the parameterization forms. A novel result of their work is the appearance of the transition layer that separates interior and BL physics.

The goal of the present study is to assess any potential impacts of the FMCD near-boundary eddy flux parameterization in our OGCM solutions. We propose to start with a preliminary, simplified version of the FMCD scheme that requires only minor modifications of our existing model. Furthermore, our current implementation is only for the surface boundary because we expect the model response to be larger there. In section 2 we briefly describe the OGCM. The implementation details are presented in section 3 and the appendix. Section 3 also includes a discussion of our simplifications and the related assumptions in comparison to FMCD. The results are given in section 4. Finally, a summary and concluding remarks are given in section 5.

## 2. Ocean model

The model is the ocean component of the National Center for Atmospheric Research (NCAR) Community Climate System Model version 3 (CCSM3). It is a level-coordinate model based on the Parallel Ocean Program (POP 1.4) of the Los Alamos National Laboratory (Smith et al. 1992; Dukowicz and Smith 1994; Smith et al. 1995). The model solves the primitive equations in general orthogonal coordinates in the horizontal with the hydrostatic and Bousinesq approximations. A linearized, implicit free-surface formulation that allows variations of the surface layer thickness is used for the barotropic equations. However, the global integral of the ocean volume remains constant because the freshwater fluxes are treated as virtual salt fluxes using a constant reference salinity. Here, we present a summary of the ocean model setup, physics, and parameters. Further model details are given in Smith and Gent (2004b) and Danabasoglu et al. (2006).

We use the nominal 3° horizontal resolution version of the ocean model detailed in Yeager et al. (2006). This coarse-resolution version has 100 (zonal) × 116 (meridional) × 25 (vertical) grid points. As in other versions of the ocean model, the grid North pole is displaced into Greenland. The finest meridional resolution occurs at the equator with 0.6°. The vertical resolution monotonically increases from 8 m near the surface to about 500 m in the abyssal ocean. The minimum and maximum ocean depths are 26 and 5000 m, respectively.

The model tracer equations use the GM90 isopycnal transport parameterization. Some relevant aspects of this scheme and the modifications due to the implementation of the FMCD parameterization near the ocean surface are presented in section 3. The anisotropic horizontal viscosity formulation for the momentum equations follows that of Large et al. (2001), as generalized by Smith and McWilliams (2003). Unless lower values are required for numerical stability, the minimum horizontal viscosity is 1000 m^{2} s^{−1}. The K-profile parameterization (KPP) of Large et al. (1994), as modified by Danabasoglu et al. (2006), is used as the vertical mixing scheme, including its double-diffusive mixing part. In the oceanic interior, similar to Bryan and Lewis (1979), the background internal wave mixing diffusivity varies in the vertical from 0.1 × 10^{−4} m^{2} s^{−1} near the surface to 1.0 × 10^{−4} m^{2} s^{−1} in the deep ocean with the increase occurring at about 1000-m depth to crudely represent enhanced mixing over rough topography (Ledwell et al. 2000). The vertical viscosity has the same shape but is 10 times larger.

The bulk forcing scheme described in Large et al. (1997) and Large and Yeager (2004) is used to compute the surface fluxes of heat, salt, and momentum. The present simulations utilize the *normal year forcing* developed by Large and Yeager (2004). This consists of single annual cycles of all the needed datasets and it can be used repeatedly without initiating any spurious transients. When observed sea surface temperatures [a blending of Levitus et al. (1998) and Steele et al. (2001)] are used with these *normal year forcing* datasets, the net heat flux shows a negative bias of about 5 W m^{−2}. In the present work this bias is crudely eliminated by adding uniformly 5 W m^{−2} to the longwave downward heat flux component.

The absorption of solar radiation is based on spatially varying, monthly mean climatologies of ocean surface chlorophyll concentrations inferred from limited satellite ocean color measurements (Ohlmann 2003). A climatological river runoff distribution described in Large and Yeager (2004) contributes to the surface salt fluxes. We do not use a sea ice model. Instead, the daily mean sea ice concentrations from the Special Sensor Microwave Imager (SSM/I) dataset (Comiso 2006) are used to determine the ice extent. To prevent any unbounded local salinity trends due to the lack of any appreciable feedbacks between the salt fluxes and model surface salinities, we apply a weak salinity restoring (corresponding to an 8-month time scale over 8 m) to observed monthly mean climatology [a blending of Levitus et al. (1998) and Steele et al. (2001)] with zero global mean. Further details of the forcing formulation, including the treatment under ice-covered regions, are given in Large and Yeager (2004) and Danabasoglu (2004).

## 3. Implementation of a simplified version of the FMCD parameterization

The representation of eddy fluxes in POP presently consists of an isopycnal diffusion (Redi 1982) and an eddy-induced velocity (GM90) represented as a skew flux (Griffies 1998) proportional to the isopycnal slope. In all experiments we use 800 m^{2} s^{−1} for both the isopycnal, *A _{I}*, and thickness,

*A*

_{ITD}, diffusivities (the thickness diffusivity multiplied by the isopycnal slope gives the skew flux). Following Large et al. (1997), the diffusivities are tapered to zero when either the isopycnal slope, |

**S**|, is too steep or the ocean surface is approached using

the asterisk denotes the untapered values of these diffusivities, 800 m^{2} s^{−1}. The first tapering, *f*_{1}, is necessary to prevent diffusive numerical instabilities (Cox 1987) and its form is

where *S*_{max} is an allowable maximum slope: *f*_{1} is set to zero for |**S**| ≥ *S*_{max}. We note that (2) is a computationally less expensive alternative to the *f*_{1} function given in Danabasoglu and McWilliams (1995) and creates similar distributions. In our experience, we find that the commonly used *S*_{max} value of 0.01 can occasionally produce *A _{I}* =

*A*

_{ITD}≈ 0, particularly in isolated regions of the abyssal ocean. This behavior can lead to spurious numerical oscillations and instabilities. To suppress these, we prefer not to introduce any horizontal diffusivity in the ocean interior that results in undesirable diapycnal diffusion. Instead, we allow the model to use isopycnally oriented mixing to the extent possible with the small-slope approximation of the mixing tensor, that is, |

**S**|

^{2}≪ 1. In our view, this restriction is sensibly obeyed when |

**S**|

^{2}≈ 0.1, thus producing

*S*

_{max}≈ 0.3 (Danabasoglu et al. 2006). With this value of

*S*

_{max},

*f*

_{1}is zero by |

**S**| = 0.18 (see also Fig. B1 in Large et al. 1997).

The second tapering, *f*_{2}, is constructed to attenuate both isopycnal diffusion and eddy-induced velocities near the ocean surface where the KPP mixing scheme and the isopycnal transport parameterization have competing effects: the former reduces stratification, while the latter increases stratification by slumping isopycnals. The rationale for this approach is that eddy transport is associated with correlations between eddy velocities and isopycnal thickness displacements. Thickness variations are inhibited when isopycnals outcrop. Thus, a tapering is applied at all depths, *d*, within a distance *D* from the surface, given by

In (3), *R* is the Rossby deformation radius, representing the preferred horizontal length scale of the baroclinic eddies. For simplicity, it is determined from *R* = *c*/*f*, subject to an additional restriction of 15 km ≤ *R* ≤ 100 km. Here, *c* = 2 m s^{−1} is a typical value for the first baroclinic wave speed, and *f* is the Coriolis parameter. With this prescription, *R* is constant at 100 km equatorward of 8° latitude and no other equation is used for the equatorial deformation radius. The *f*_{2} taper is constructed using

Without this taper, the near-surface eddy-induced circulation can become unrealistically large, that is, *O*(100 Sv; Sv ≡ 10^{6} m^{3} s^{−1}), especially when the vertical resolution is *O*(10 m) near the surface. As discussed in appendix B of Large et al. (1997), this function substantially reduces, but does not completely eliminate, the near-surface eddy fluxes. We use both tapering functions in our CONTROL case.

To implement the FMCD parameterization, we must estimate two vertical length scales: the boundary layer depth (BLD) and the transition layer thickness (TLT). We define their sum as the diabatic layer depth (DLD) over which the upper-ocean eddy fluxes depart from their interior formulas. Here, BLD is determined by the KPP vertical mixing scheme as the shallowest depth at which a bulk Richardson number exceeds a specified critical Richardson number for the first time (Large et al. 1994). Five passes of a grid-scale, five-point spatial filter are applied to this depth to eliminate any small-scale features. For simplicity, we choose to apply no filtering in time. FMCD defines the transition layer as the layer containing all isopycnals within an averaging area and time interval that are intermittently exposed to strong turbulent mixing. Thus, TLT is defined by the range of isopycnals that can be lifted into the BL by subgrid-scale eddy heaving and/or the subgrid-scale (subsynoptic) variations of BLD. FMCD show that the expression given by (3) is a reasonable proxy for both of these effects. We compute *D* at each grid point in a vertical column below the BL and search for the shallowest depth where *D* does not reach the BL. This is equivalent to finding the depth *d* where

is satisfied for the first time. TLT is simply obtained from

Because we do not impose any maximum slope restrictions in these computations, TLT can become large or even reach the ocean bottom when the column is still very weakly stratified below the BL. Conversely, TLT can vanish if |**S**| = 0.

We use our model to obtain the boundary layer and transition layer distributions presented in Figs. 1a and 1c, respectively. We label this first experiment as DB. The BL is deepest in the northern North Atlantic, along the Southern Ocean, and in the tropical central Pacific. It is relatively shallow in low latitudes and in the Arctic Ocean. The global-mean BLD is 48 m. The time-mean TLT has a global-mean thickness of 15 m. It exceeds 40 m in the northern North Atlantic, western subtropical North Atlantic and North Pacific, east of Australia, and the Indian Ocean sector of the Southern Ocean. The central subtropical regions of all basins show extensive areas of TLT between 15 and 30 m.

We find that the mixed layer depth (MLD) is usually close to the sum of BLD and TLT (Fig. 1). Following Large et al. (1997), we define MLD as the shallowest depth where the local, interpolated buoyancy gradient matches the maximum buoyancy gradient between the surface and any discrete depth within that water column. In contrast to the BLD, which can exhibit rapid fluctuations during which mesoscale eddies are unlikely to change, the MLD represents a lower-frequency, that is, time-filtered, envelope of the surface boundary layer region. In particular, the mixed layer in our simulations records the maximum depth of the BL after sustained deep mixing events, and it can be considered a measure of the deepest penetration of turbulent mixing in the recent past (see section 2c of FMCD for further discussions).

The finding that MLD ≈ BLD + TLT indicates that MLD includes most of the transition layer. Furthermore, it implies that the tracers are well mixed in the vertical and that horizontal density gradients are not likely to veer too much with depth in the transition layer. It is quite likely that the lack of vertical structure in the transition layer is a result of the coarse horizontal and vertical model resolution. However, for present purposes we can safely neglect any vertical variation in tracer gradients within MLD (≈ BLD + TLT). In this limit the FMCD parameterization for the eddy-induced streamfunction reduces to

where **S*** _{I}* is the isopycnal slope at the base of the transition layer and

**z**is the vertical coordinate, positive upward. The vertical structure function

*G*(

*z*) is defined by

with

where *λ* is the vertical length scale for the eddy fluxes below the transition layer, *ρ* is the local potential density, and **∇*** _{h}* is the horizontal gradient operator. We note that the approximation in (8c) is valid if ∂

*|*

_{z}**∇**

*| is small. The parameterization for tracer diffusion reduces to*

_{h}ρwhere *τ* is a generic tracer.

With these parameterization forms, the eddy-induced velocity has no vertical shear within the BL and it develops a linear vertical shear within the TLT to match the interior values at its base. In the appendix we give the full parameterization expressions as implemented in the model. We note that with the new scheme we do not need the *f*_{2} taper any longer. We thus eliminate *f*_{2} from (1)

in the integrations with the FMCD parameterization. In addition, *f*_{1} is only applied in the interior below the transition layer. We note that further sensitivity experiments that also use *A*_{ITD} = *A**_{ITD} in the ocean interior, that is, no *f*_{1}, lead to numerical instabilities with our present parameter choices.

To test the sensitivity of the parameterized eddy transport to the structure function *G*(*z*), we run an alternative experiment DM, where we set BLD = MLD. In this way we eliminate any eddy-induced shear through the entire MLD. We still allow shear in a transition layer, identified by (3), below the mixed layer to avoid any sharp discontinuities in eddy-induced velocities between the base of the mixed layer and the interior. In DM, we apply the same spatial filter to MLD as for BLD to remove lateral variability at the grid scale. In both DB and DM, the TLT values (Figs. 1c and 1d) are in the range 2–100 m. As in other ocean models, the vertical resolution (15–28 m in the 49–130-m depth range) is not sufficient to resolve (i.e., have at least half a vertical grid cell) the lower end of this range. This is not the same as having 0 m for TLT because these nonzero TLT values are used in all of the equations. Nevertheless, it implies a somewhat abrupt transition in eddy fluxes between the diabatic surface layer and the adiabatic interior. Finally, we run a simulation, denoted as TLT0, where BLD = MLD and TLT is set to zero so that there is no shear in both the boundary and the transition layers. Simulations DB, DM, and TLT0 give very similar results, suggesting that at the vertical resolution used in this paper vertical shears are weak in the transition layer and, hence, they do not affect the simulations much.

The distributions of the starting depth of the interior region (i.e., DLD) presented in Fig. 1e for DB and for the DM − DB difference in Fig. 1f clearly verify MLD ≈ BLD + TLT. One exception to this is the northern North Atlantic where the deeper MLD in DM results in a 100-m-deeper DLD than in DB. The time-mean TLT from DM (Fig. 1d) is uniformly shallower than in DB. The patterns of MLD (Fig. 1b) are very similar to the ones shown in Fig. 1a for the BLD, but the MLD with a global mean of 71 m is deeper than the BLD. Because the BLD is shallower than the MLD, TLT computations in DB are exposed to larger isopycnal slopes in this depth range than in DM, resulting in thicker transition layers. In DM, this large slope range is already within the mixed layer.

### a. Summary of numerical experiments

Table 1 lists the coarse-resolution numerical experiments, including their integration lengths. CONTROL is the only case with the original formulation where *A _{I}* and

*A*

_{ITD}are given by (1). Experiments DB and DM use BLD and MLD, respectively, to represent the well-mixed surface layer. Because their solutions are very similar, we somewhat arbitrarily choose DM as our primary experiment to document the sensitivity of the model solutions to the near-surface eddy flux parameterization. We note that because TLT is rather small, DM is consistent with our assumption that vertical variations of momentum and tracers are weak in the transition layer. We also consider DM the base case for a sensitivity experiment, labeled TLT0, in which the transition layer has no thickness. All of the cases are initialized with January-mean climatological potential temperature and salinity (Levitus et al. 1998; Steele et al. 2001 in the Arctic Ocean) and no velocity. Most of our analysis is based on the time-mean data obtained using the last 20 years of integrations.

We also briefly compare our present results with the new scheme to the solutions from a 1/10° eddy-resolving global simulation (labeled as ER) with the POP model (Maltrud and McClean 2005). This integration was recently extended to calendar year 2000, using daily National Centers for Environmental Prediction–National Center for Atmospheric Research (NCEP–NCAR) reanalysis data (Kalnay et al. 1996) and various monthly observational data products (see Maltrud and McClean (2005) for forcing details). We use the 3-yr mean fields for years 1998–2000 in our analysis.

## 4. Results

We first show the resulting distribution of time- and zonal-mean, global *A _{I}* from DM along with its differences from CONTROL and DB in the upper ocean in Fig. 2. As expected, the largest changes from CONTROL occur within DLD, where

*A*in DM is smaller than in CONTROL by as much as 400 m

_{I}^{2}s

^{−1}(Fig. 2b). Another region of reduction by 200–300 m

^{2}s

^{−1}is at about 60°N, resulting from deeper seasonal DLD due to enhanced wintertime mixing in the northern North Atlantic. We note that these locations of reduced

*A*essentially indicate where the downgradient horizontal diffusion contributes to, or even dominates, tracer mixing. Figure 2c shows that

_{I}*A*in DM is also less than in DB, particularly in the tropical regions as a consequence of deeper DLD in DM (see Fig. 1). Below about 1000 m, the

_{I}*A*distributions for CONTROL and DB are basically the same as the distribution from DM given in Fig. 2a. The smaller

_{I}*A*magnitudes, particularly evident just above the bottom topography, are due to the discrete implementation of the no-flux boundary conditions on tracers at the ocean bottom where both

_{I}*A*and

_{I}*A*

_{ITD}are set to zero in the bottom halves of the bottom grid cells just above the topography as discussed in section 1.

The most significant and notable effects of the new parameterization occur in the near-surface eddy-induced meridional overturning circulation (MOC). Figures 3a and 3b present the zonally integrated eddy-induced MOC from DM and CONTROL, respectively, showing the elimination of the shallow near-surface circulations with FMCD. In CONTROL, these circulations exceed 20 Sv in the latitude band of the Southern Ocean and reach 10 Sv elsewhere. Unfortunately, observational (inferred) estimates of an eddy-induced circulation are very limited. In particular, details of its near-surface structure, specifically the possibility of strong, shallow recirculation cells, are further limited by the vertical resolution of the measurements and their geographical sparseness. An observational estimate, with a fairly good near-surface resolution, is provided by Roemmich and Gilson (2001) along a repeated hydrographic line in the North Pacific. Their work does not support the presence of strong and shallow near-surface eddy-induced circulations. Furthermore, such surface-trapped circulations are also absent in the total (i.e., the sum of eddy-induced and Eulerian-mean) MOC distributions from ER given in Fig. 3g. Indeed, the total MOCs from ER and DM (Fig. 3e) are very similar. In contrast, the CONTROL circulation (Fig. 3f) clearly shows that the eddy-induced circulation is strong enough to change the sense of the total MOC near the surface. Therefore, we believe that FMCD produces a better representation of the ocean circulation. Below about 200-m depth, both DM and CONTROL have very similar eddy-induced circulations (Figs. 3a and 3b). The accompanying Eulerian-mean MOC from DM and its difference from CONTROL are given in Figs. 3c and 3d, respectively. The distribution of differences shows that the DM and CONTROL MOCs are alike with differences rarely exceeding 1 Sv. The circulation associated with the North Atlantic Deep Water has identical maxima (18.1 Sv) in both, but its circulation is slightly lower in DM by about 1 Sv farther north. There is a similar weakening of the cell associated with the Antarctic Bottom Water (AABW) in DM where its maximum is about 0.8 Sv smaller than in CONTROL (10.6 and 11.4 Sv, respectively).

The time-mean global northward heat transport from DM and its difference from CONTROL are shown in Fig. 4, including the contributions of the transport components. Also given for comparison in Fig. 4a is an estimate for the total transport from Trenberth and Caron (2001) computed using residually derived surface fluxes based on the NCEP–NCAR data. In the figure, the diffusive component includes both isopycnal and horizontal diffusion contributions in DM. North of about 40°S, the total transport is dominated by the Eulerian-mean advection component and the differences between DM and CONTROL are much too small to produce any significant changes in these transports. The maximum global transport is 1.32 PW at 19°N. The principal contribution to this transport is from the Atlantic basin (not shown) where the maximum northward transport is 0.85 PW at 15°N. In this basin the direct estimates for this maximum transport are 1.2–1.3 PW between 14° and 24°N with estimated errors of 0.3 PW (Bryden and Imawaki 2001). Consequently, both the global and Atlantic transports are underestimated in all model experiments. This is typical of our coarse-resolution ocean model and can be partially attributed to both the *normal year* datasets used to force the model and a third-order upwind tracer advection scheme. In the Southern Ocean the southward transport is relatively small with a maximum of 0.36 PW at 56°S, and contributions from all components are important. In DM the local minimum at 42°S is 0.1 PW, representing an increase of 0.05 PW from CONTROL, mostly due to the larger diffusive transport (Fig. 4b). As in Danabasoglu and McWilliams (1995), the Eulerian-mean and eddy-induced advection components partially compensate each other, particularly at high southern latitudes. This is also present in the difference plot (Fig. 4b) where large fractional changes in the Eulerian-mean and eddy-induced advection act largely to cancel each other poleward of about 35°S. Therefore, the diffusion component dominates the total difference here. Further partitioning of the diffusive heat transport component to its isopycnal and horizontal diffusion parts is discussed at the end of this section. We note that the new parameterization tends to reduce the southward transport by the eddy-induced advection, likely due to the elimination of the strong near-surface eddy-induced circulations.

Arguably, the most dramatic effects of the elimination of the strong near-surface eddy-induced circulation with FMCD are evident in the vertical structure of the heat transports (or more accurately the zonally integrated total temperature fluxes in °C). Figure 5 shows the vertical profiles of the zonally integrated, time-mean total, that is, Eulerian-mean plus eddy-induced, advective heat transport at 49.4°S for the upper ocean from CONTROL and DM in comparison with the profile from ER. These profiles are typical of the latitude band between 30° and 60°S where the eddy-induced transport is particularly strong in CONTROL (see Fig. 3b). They are given in TW m^{−1} so that their vertical integrals produce the transports given in Fig. 4. Although both CONTROL and DM have similar integrated total advective transports (−0.145 and −0.138 PW, respectively), Fig. 5 shows that their vertical structures differ substantially. In particular, the DM profile is in remarkably good agreement with the profile from ER, both in magnitude and shape. In contrast, the CONTROL profile has alternating northward and southward transports in the upper 200 m, reflecting the dominance of the near-surface eddy-induced circulation. Supported with the ER data, we believe that the DM profile is more realistic and will improve the upper-ocean structures of all tracers.

We present the annual-mean MLD DM − CONTROL difference distribution in Fig. 6. With FMCD, the MLDs increase significantly in the Southern Ocean and northern North Atlantic, including the Labrador Sea, by as much as 60 and 250 m, respectively. There are only a few isolated regions with noticeable reductions in MLD in DM. At low latitudes the figure shows no substantial differences. We compute the global-mean MLDs as 63.7 and 71 m, respectively, for CONTROL and DM. The boundary layer differences (not shown) are very similar to the mixed layer ones, both in magnitude and spatial distributions. An analysis of the seasonal cycle of MLDs (not shown) reveals that the deepening in DM occurs in the nonsummer months in both hemispheres, with the winter months showing the largest deepening. These deepened mixed layer regions of DM compare more favorably than in CONTROL with observational estimates based on the monthly mean Levitus et al. (1998) and Steele et al. (2001) data. However, given that the data coverage is extremely poor at high latitudes, particularly during the winter months, we believe that a more detailed comparison of MLDs with observations is not justified.

Another, unexpected change due to the new parameterization occurs in the long-term evolution of the model’s global-mean potential temperature, 〈*θ*〉, presented in Fig. 7. In CONTROL, after a short initial warming period in which 〈*θ*〉 warms by about 0.05°C compared to the initial value (also representing the observed global mean), the model continuously loses heat. By year 2000, 〈*θ*〉 is more than 0.4°C colder than observed. We compute a heat flux value of −0.06 W m^{−2} based on the 〈*θ*〉 trend during the last 20 yr of this case. In contrast, all of the integrations with the new scheme (Table 1) show a longer and larger initial warming period followed by a much more gradual cooling. For example, both DB and DM get warmer than observed by 0.14° and 0.17°C, respectively, during the first 200 yr. Their 〈*θ*〉 values differ very little from the observed at year 2000. Again, based on the last 20 yr, we obtain a heat flux of only −0.02 W m^{−2} for both, a factor of 3 reduction compared to the CONTROL tendency.

The time- and zonal-mean, global *θ* from CONTROL and DM are differenced from the observations [a blending of Levitus et al. (1998) and Steele et al. (2001)] in Fig. 8. In CONTROL the difference plot shows that the bulk of the cooling (discussed above) occurs below 1000-m depth where the model solution is colder than observed, by more than 1.2°C. In contrast, this cold bias is significantly reduced both in spatial extent and magnitude in DM, resulting in better agreement with observations. South of about 45°S, the DM zonal-mean *θ* is now a little warmer than observed. In the upper ocean, above 1000 m, there are no appreciable differences between CONTROL and DM, both showing almost identical warm biases compared to observed, similar to some previous studies (e.g., Large et al. 1997). Due to the inadequate representation of the ocean bottom topography and, in particular, the Icelandic Ridge in this coarse-resolution model, both CONTROL and DM display relatively large differences from observations at about 65°N, with the DM difference being slightly larger and displaced somewhat deeper than in CONTROL.

To gain further insight into the abyssal warming of the cases with the FMCD parameterization, we next examine changes in ocean ventilation, considering the ideal age tracer (Thiele and Sarmiento 1990). This tracer is governed by the same conservation equations as potential temperature and salinity. It is set to zero at the ocean surface and has a source term of one unit per year of model integration. So, the younger waters indicate recent contact with the surface, while regions of little ventilation contain the oldest waters. In all cases the model is initialized with zero ideal age. We present the time- and zonal-mean ideal age distribution from DM and its difference from CONTROL in Fig. 9. After 2000 years of integration, the ideal age has equilibrated only in the upper few hundred meters, and there are still some adjustments in the abyssal ocean. The oldest waters have a zonal-mean age of more than 1500 yr, occurring in the Pacific basin (not shown). The difference plot reveals that ideal age in DM is uniformly older than in CONTROL below 2000 m. This figure also suggests that the aging is associated with a slightly weaker AABW circulation in DM (Fig. 3d), where ideal age is already older by more than 75 yr by 500-m depth south of 65°S. In contrast, the near-surface mixing appears to be more vigorous in DM than in CONTROL, as implied by younger waters in the upper 250 m. These waters then appear to spill over to the lower latitudes, but they do not penetrate below 2000 m.

Figure 10 shows DM − CONTROL differences for surface density, temperature, and salinity in the latitude band of the Southern Ocean that includes the AABW formation regions. The surface waters in DM are broadly lighter than in CONTROL south of 60°S. Indeed, these lower densities extend farther north to 50°S in the Atlantic and Indian sectors of the Southern Ocean. There are, however, a few small regions with higher surface densities just east of the Antarctic Peninsula and off the Antarctic continent. Because the salinity contraction coefficient is about an order of magnitude larger than the thermal expansion coefficient at these temperatures, the changes in the surface density are largely dictated by the salinity changes (Fig. 10c). In surface temperature (Fig. 10b), there are no well-defined latitudinal difference patterns, as in salinity. Furthermore, the DM and CONTROL zonal-mean surface temperature (not shown) differ by only 0.01°C south of 60°S. Nevertheless, the reduced surface density due to fresher surface waters in DM may imply a weaker equator to pole buoyancy torque for the overturning circulation, thus providing a possible rationalization for the weaker AABW circulation.

The associated surface heat flux DM − CONTROL zonal-mean difference (not shown) reveals slightly more negative heat flux, that is, cooling, in DM than in CONTROL at these southern high latitudes. This cooling is compensated by an increased southward heat transport in DM (Fig. 4b) due to the diffusion contribution. In Fig. 11a, we further partition this component to its horizontal and isopycnal parts in DM. The figure shows that DM and CONTROL have similar isopycnal transports and that the horizontal diffusion component is the reason for increased southward transport in DM. The vertical profiles of the zonally integrated, time-mean diffusive transports at 49.4°S (typical of the latitude band between 20° and 65°S) are plotted in Fig. 11b. In CONTROL the isopycnal diffusion is the only component, approaching zero near the surface. In DM it has a similar behavior. In contrast, there is a significant near-surface southward transport due to the horizontal diffusion component that dominates the total near-surface diffusive transport in DM. We believe that this increased southward transport more than compensates the increased surface cooling and plays a crucial role in the abyssal ocean warming seen in cases with the FMCD parameterization. However, we do not have a detailed mechanistic understanding of the relationship between the near-surface properties and the abyss.

All of the improvements resulting from the FMCD parameterization are also present in DB and TLT0, including the changes in the long-term evolutions (Fig. 7) and distributions of *θ*. These three cases have in common similar distributions and magnitudes of the upper depth of the interior region, that is, DLD. In DB and DM, this is a result of the partial compensation of the surface layer depths and TLT as noted earlier in section 3. Between TLT0 and DM, however, it is essentially due to the overall thinness of TLT in DM, which can be approximated as zero thickness, as in TLT0.

## 5. Summary and concluding remarks

We have implemented a simplified version of the near-boundary eddy flux parameterization of FMCD for the surface boundary in our OGCM. The extant eddy transport schemes commonly used in OGCMs (e.g., Gent and McWilliams 1990, GM90) are generally developed for the quasi-adiabatic interior and are not valid near the boundaries. Therefore, the usual practice is to taper the effects of parameterized eddy fluxes as the surface (or any other boundary) is approached. This attenuation of the eddy effects is not physically consistent with observational evidence, particularly near the surface where diabatic mesoscale fluxes may dominate mixing. FMCD includes the effects of these upper-ocean diabatic fluxes and eliminates the ad hoc, near-surface taper functions.

The experiments conducted with FMCD show significant improvements compared to available observations and to solutions from an eddy-resolving model and from a CONTROL case that uses the usual tapering approach. The most significant of these improvements is the elimination of the unphysical, shallow, and strong near-surface eddy-induced circulations that are not supported by any observational evidence. This subsequently leads to the disappearance of the alternating northward and southward heat transports present in the CONTROL upper-ocean heat transport profiles. The profile with FMCD matches that of the eddy-resolving model quite remarkably. We believe that this may lead to further improvements in the upper-ocean distributions of passive tracers, and possibly produce some changes in the associated surface fluxes.

Another substantial change occurs in the abyssal distributions of potential temperature. In particular, the deep cold bias of CONTROL is significantly reduced with FMCD, both in spatial extent and magnitude, resulting in better agreement with observations. This change is also reflected in the global-mean potential temperature time series that show very reduced trends with FMCD compared to CONTROL. The meridional overturning circulations show a slightly weaker AABW circulation with FMCD, leading to older abyssal waters. We believe that this weaker circulation allows more time for diapycnal diffusion to heat the abyssal ocean. An examination of the heat transport components reveals that the increased southward heat transport due to the horizontal diffusion component more than compensates the increased surface cooling in the AABW formation regions. We speculate that this plays an important role in the abyssal ocean warming with FMCD.

Our simplifications of the FMCD parameterization are based on the findings that the mixed layer depth (MLD) is usually a good approximation for the sum of the boundary layer depth (BLD) and the transition layer thickness (TLT) and that the shears and vertical tracer gradients in the transition layer are weak. These indicate that MLD includes most of the transition layer. Furthermore, the implication is that the tracers are well mixed in the vertical and that horizontal density gradients are not likely to veer too much with depth in the transition layer. It is quite likely that the lack of vertical structure in the transition layer is a result of the coarse horizontal and vertical model resolution. Nevertheless, we assume that 1) the horizontal gradients do not veer much within TLT, 2) they remain constant across the BLD and TLT, 3) the vertical profile of fluxes is linear in the vertical in both BLD and TLT, and 4) the vertical component of the diffusive fluxes can be ignored in the BLD and TLT. The last assumption is based on our analysis that shows that the vertical diffusive fluxes are usually rather small in comparison with the horizontal diffusive fluxes within the boundary layer. The primary consequence of these assumptions is that the eddy-induced velocity is not guaranteed to reduce the mean potential energy, as in the FMCD parameterization, because the direction of the eddy-induced velocity is set by the buoyancy gradients at the transition layer base in our simplified implementation instead of being set by the local buoyancy gradients (see FMCD). Another discrepancy is the missing continuity of the diffusive fluxes across the transition layer base in the present implementation compared to their full continuity in FMCD.

Obviously, the largest uncertainty in FMCD is in the determination of the TLT. Our choice here is based on the local isopycnal slopes and the preferred horizontal length scale of the baroclinic eddies. The isopycnal slopes get steeper as the surface is approached. Therefore, our formulation produces thicker transition layers when the surface layer is too shallow to include these steep slope areas. Conversely, the deeper surface layers lead to thinner transition layers. We think that this is a simple and reasonable formulation, but there is a clear need for better observational guidance. One such observational estimate has recently been provided by Johnston and Rudnick (2007, manuscript submitted to *J. Phys. Oceanogr.*), using a SeaSoar dataset and ADCP measurements in the North Pacific. They find that the TLT is typically on the order of 10%–20% of a specified BLD. Their result is not very sensitive to how they choose to define their boundary layer depths. This estimate gives us some confidence that the TLT distributions and magnitudes given in Figs. 1c and 1d are not completely unrealistic. Using 0 m for TLT is numerically simple and attractive, but not supported by observations.

Our simplified version of the FMCD parameterization reduces to a similar form described in Griffies et al. (2005) in the limit of no transition layer. Although our sensitivity experiment with 0-m TLT produces similar results to those from finite thickness transition layer cases, without a transition layer, the eddy-induced velocities become discontinuous at the base of the boundary layer. In addition, Griffies et al.’s (2005) implementation appears to require a convective adjustment scheme within the surface diabatic layer after the application of the KPP vertical mixing scheme. We note that no such unphysical convective adjustment is necessary in the present study, even in the vanishing TLT case. We do not have a clear understanding of why such a procedure is necessary in Griffies et al. (2005), because such convective instabilities (presumably generated by mesoscale eddies) are not physically justified, and they will adversely modify the exchange of fluxes between the boundary layer and the interior. We speculate that the details of how the surface diabatic layer depth is determined are important as there are known differences between the present and Griffies et al. (2005) approaches. For example, the diabatic layer depth computations in Griffies et al. (2005) involve an artificial critical slope, which can significantly change this depth.

Finally we address a numerical issue associated with the finite thickness of the transition layer. Given its temporal and spatial variability, an OGCM cannot be expected to resolve the transition layer at all time steps and at all locations. Here, we use *resolve* loosely to mean that the transition layer has at least half a vertical grid cell. In fact, our TLT formulation can often give thicknesses that are less than half a vertical grid cell, but this nonzero value is used in all of our equations and is not the same as having 0-m TLT. One obvious remedy is to use a minimum TLT that can be resolved by the model’s vertical grid regardless of depth, but this implies an undesirable grid dependency. We believe that resolving TLT is not very crucial to benefit from the improvements due to FMCD.

## Acknowledgments

The authors thank W. G. Large, P. R. Gent, and J. Marshall for many helpful discussions and suggestions; F. O. Bryan for kindly providing the 1/10° eddy-resolving model data; and S. Griffies and an anonymous reviewer for helpful comments. This research was partially supported by NSF Grant OCE-0336827 for the Climate Process Team on Eddy Mixed-Layer Interactions (CPT-EMILIE). The computational resources were provided by the Scientific Computing Division of the National Center for Atmospheric Research (NCAR). NCAR is sponsored by the National Science Foundation.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

### APPENDIX

#### Implementation of a Simplified Version of the FMCD Parameterization

We split the eddy-induced vector streamfunction given by Eq. (7) into its boundary layer, **Ψ**_{BL}, and transition layer, **Ψ**_{TL}, expressions as follows:

and

The eddy-induced velocity is obtained using **∇** × **Ψ**. In the above equations *z* is the vertical coordinate, positive upward, and *η* is the displacement of the free surface relative to *z* = 0. Also, BLD and TLT are the boundary layer depth and transition layer thickness, respectively, and the upper (starting) depth of the interior region, DLD, is given by DLD = BLD + TLT. The two functions **Ψ**_{o} and **Φ** are chosen such that **Ψ** and its vertical derivative are continuous across the base of BLD and the base of TLT. These constraints then yield

and

where subscript *z* denotes partial differentiation. In (A2a) and (A2b), **Ψ*** _{I}* is the interior eddy-induced streamfunction at the base of the transition layer given by the GM90 parameterization; that is,

**Ψ**

*=*

_{I}**Ψ**

_{GM}(

*z*= −DLD), where

In (A3) *ρ* is the local potential density, **∇*** _{h}* is the horizontal gradient operator, and

*A*

_{ITD}is the thickness diffusivity. We note that (A3) is equivalent to (7) of section 3 with ∂

_{z }**Ψ**

_{I}≈ −

*A*

_{ITD}(1/

*ρ*)

_{z}

_{z }**z**×

**∇**

*.*

_{h}ρIn (A1a) **Ψ**_{o} represents the vector streamfunction value at the base of BLD. We note that **Ψ**_{BL} is linear within BLD, going to zero at the ocean surface. This implies constant eddy-induced velocities with no vertical shear within BLD. These eddy-induced velocities, however, must develop a shear in the transition layer to match the interior values. For this purpose, (A1b) represents the simplest choice for **Ψ**_{TL}, which is parabolic, or equivalently linear in *z* for the eddy-induced velocities.

In our present implementation we neglect any variations of the free surface height in the above equations by setting *η* = 0 because all of the horizontal fluxes between *z* = 0 and *z* = *η* are already neglected in POP due to the linearization of the barotropic continuity equation. The existing discrete implementation of the isopycnal transport parameterization readily subdivides a vertical grid cell into a top and a bottom half. We naturally take advantage of this *doubled* vertical grid resolution and use the depths of the midpoints of the top and bottom halves as reference depths to determine if a particular half is within BLD, TLT, or the interior. We then evaluate separate **Ψ**_{I} for the appropriate top and bottom halves of the vertical grid cells that straddle the base of the transition layer. These **Ψ**_{I} are also used to compute ∂_{z}**Ψ**_{I} along the transition layer–interior interface.

With **Ψ**_{I} and ∂_{z}**Ψ**_{I} now available at *z* = −DLD, **Ψ**_{BL} and **Ψ**_{TL} can be evaluated using (A1) and (A2). The interior eddy-induced streamfunction values below the base of the transition layer are then determined from (A3). We finally note that in the experiments that use the mixed layer depth (MLD) to represent the well-mixed surface layer, we replace BLD with MLD in the above equations.

The mesoscale eddy fluxes still mix tracers along isopycnal surfaces in the ocean interior, as represented by the isopycnal diffusion tensor, **R**(*A _{I}*,

*τ*) =

**∇**· [

*A*

_{I}**K**·

**∇**

*τ*]. Here

*τ*is a generic tracer and

**K**is a second-rank mixing tensor with the small-slope approximation (see Redi 1982; Cox 1987, GM90). Within the boundary and transition layers, the parameterization for diffusive tracer flux is given by

In (A4), *A**_{H} represents the untapered horizontal eddy diffusivity. Of course, the rate of mixing within both the interior and BL remains the same, so *A**_{H} has the same value as the interior *A _{I}*. The vertical function

*c*(

*z*) is defined by

Note that we drop the vertical component of **F**(*τ*) within BL because it is found to generate negligible flux divergences at the resolution used in this paper.

## Footnotes

*Corresponding author address:* Dr. Gokhan Danabasoglu, National Center for Atmospheric Research, P.O. Box 3000, Boulder, CO 80307. Email: gokhan@ucar.edu