Cette annexe illutre en 5 papiers et une *Newsletter* certaines de mes activités dé-
crites dans les pages précédentes.

[P14] *Influence of numerical schemes on current-topography interactions in ¼° global ocean simulations* (2007, 73 citations) concerne la calibration de la configuration DRAKKAR ORCA025 (global ¼°) de NEMO, en l'occurrence l'impact bénéfique des "partial steps" et de la formulation de l'opérateur d'advection de quantité de mouvement sur la circulation générale.

[P22] *Impact of global ocean model resolution on sea-level variability with emphasis on interannual time scales* (2010, 69 citations ; avec M. Juza) se concentre sur la sensibilité de la variabilité interannuelle de SLA à la résolution horizontale de simulations globales. Il montre notamment qu'à chaque doublement de résolution (2°-1°-½°-¼°), cette variabilité interannuelle tend vers la distribution et les niveaux observés par altimétrie, mais se décorrèle temporellement des séries temporelles observées lors du passage au régime turbulent (¼°).

[P25] *Sea-Level Expression of Intrinsic and Forced Ocean Variabilities at Interannual Time Scales* (2011, 29 citations ; avec M. Juza) examine et confirme l'hypothèse émise en conclusion du papier précédent : une variabilité intrinsèque interannuelle de nature chaotique est spontanément émise par l'océan global en régime turbulent (simulation climatologique au ¼°) ; ce phénomène explique dans certaines régions la quasi-totalité de la variance de SLA obtenue dans le *hindcast* correspondant (avec forçage interannuel et synoptique), et pourrait expliquer la décorrélation temporelle mentionnée plus haut.

[P33] *Intrinsic variability of sea-level from global 1/12° ocean simulations: spatio-temporal scales* (2015 ; avec G. Sérazin et S. Grégorio) prolonge l'étude précédente. Sur la base de paires de simulations globales (1 *hindcast*, 1 run climatologique) au ¼° et au 1/12°, nous y montrons notamment que la variabilité de SLA est très largement intrinsèque et chaotique à mésoéchelle, et qu'elle le demeure (malgré une contribution moindre à la variance totale) aux échelles décennales et de milliers de kilomètres dans les zones de forte activité mésoéchelle. Cette correspondance spatiale suggère que l'interaction entre tourbillons pourrait localement produire cette variabilité chaotique basse fréquence grande échelle.

[P34] *Intrinsic variability of the Atlantic meridional overturning circulation at interannual-to-multidecadal timescales* (2015 ; avec S. Grégorio et G. Sérazin) explore l'empreinte de la variabilité intrinsèque basse fréquence sur la circulation d'overturning Atlantique (AMOC) et son transport de chaleur associé, notamment le long des sections où ces quantités sont ou seront mesurées (RAPID, OVIDE, O-SNAP, SAMOC). Nous y montrons notamment le caractère chaotique et l'intensité substantielle de la variabilité intrinsèque d'AMOC, sur des périodes atteignant plusieurs décennies et l'intégralité du bassin.

Enfin, [A27] *Ensembles of eddying ocean simulations for climate* (2014 ; avec L. Bessières, G. Sérazin, et S. Grégorio) présente les fondements, les attendus et la stratégie du projet OCCIPUT. Cette *Newsletter* a fait suite à ma conférence invitée lors du High-Resolution Ocean Climate Modeling workshop du WGOMD (CLIVAR) à Kiel en 2014.
<figure><img src="image_27.png" /></figure>

# Influence of numerical schemes on current-topography interactions in 1/4° global ocean simulations

T. Penduff<fnref n="1" />, J. Le Sommer<fnref n="1" />, B. Barnier<fnref n="1" />, A.-M. Treguier<fnref n="2" />, J.-M. Molines<fnref n="1" />, and G. Madec<fnref n="3" />

Received: 11 June 2007 – Published in Ocean Sci. Discuss.: 25 June 2007

Revised: 25 October 2007 – Accepted: 1 November 2007 – Published: 20 December 2007

**Abstract.** The combined use of partial steps and of an energy-enstrophy conserving momentum advection scheme was shown by Barnier et al. (2006) to yield substantial improvements in the surface solution of the DRAKKAR 1/4° global sea-ice/ocean model. The present study extends this investigation below the surface with a special focus on the Atlantic and reveals many improvements there as well: e.g. more realistic path, structure and transports of major currents (Gulf Stream, North Atlantic Current, Confluence region, Zapiola anticyclone), behavior of shedded rings, narrower subsurface boundary currents, stronger mean and eddy flows (MKE and EKE) at depth, beneficial enhancement of cyclonic (anticyclonic) flows around topographic depressions (mountains). Interestingly, adding a no-slip boundary condition to this improved model setup cancels most of these improvements, bringing back the biases diagnosed without the improved momentum advection scheme and partial steps (these biases are typical of other models at comparable or higher resolutions). This shows that current-topography interactions and full-depth eddy-admitting model solutions can be seriously deteriorated by near-bottom sidewall friction, either explicit or inherent to inadequate numerical schemes.

## 1 Introduction

The three main types of Ocean General Circulation Models essentially differ by their vertical coordinate system (Griffies et al., 2000): geopotential coordinates (in so-called z-level models) and terrain-following coordinates (in so-called sigma models) remain fixed in time, while isopycnic coordinates follow density surfaces in time. More recent hybrid coordinates aim at taking advantage of these various approaches in specific regions. Intercomparison experiments

(e.g. Roberts et al., 1996; DYNAMO Group, 1997; Chassignet et al., 2000; Treguier et al., 2005) have contributed to improve basin-scale ocean models by a careful evaluation of their respective strengths and weaknesses. More specifically, these studies revealed that current-topography interactions are essential for the simulation of realistic ocean states, that these processes are handled differently by various classes of models, and that the increase of resolution tends to lessen the discrepancy between various solutions. The DYNAMO experiment (Willebrand et al., 2001) was focused on 1/3° (so-called eddy-admitting) ocean models able to partly resolve the mesoscale activity. Along with the ongoing (and beneficial) increase of basin- and global-scale ocean model resolution (Paiva et al., 1999; Smith et al., 2000; Eden and Böning, 2002; Masumoto et al., 2004;nut et al., 2007), eddy-admitting ocean models require further investigation since the 1/4° resolution is the target of future ocean climate models. In this context, important questions remain open concerning the sensitivities of such models, and the benefits of (partially) resolving oceanic eddies in future climate assessments.

The DYNAMO program revealed specific tendencies of widely-used z-level models in the eddy-admitting regime. Further investigations made during the French CLIPPER experiment (presented in Tréguier et al., 1999) highlighted additional biases in 1/3° and 1/6° eddy-admitting z-level solutions with respect to observations and sigma-coordinate solutions. It was conjectured (Penduff et al., 2002, 2005) that the robust lack of kinetic energy found at depth, possibly due to spurious topographic friction, may adversely affect current-topography interactions and explain the circulation biases found at the surface. These circulation biases (and their conjectured origin) are likely not specific to CLIPPER; they are met in other setups, other eddy-admitting and eddy-resolving models, including the first version of the 1/4° global DRAKKAR model based on CLIPPER-like numerics (Barnier et al., 2006). This latter study describes how two numerical changes (the use of partial steps and of an

<footnote n="1">¹Laboratoire des Ecoulements Géophysiques et Industriels, CNRS, UJF, INPG, Grenoble, France</footnote>

<footnote n="2">²Laboratoire de Physique des Océans, CNRS, UBO, IFREMER, Brest, France</footnote>

<footnote n="3">³Laboratoire d'Océanographie et de Climat par Expérimentation et Approche Numérique, CNRS, UPMC, IRD, MNHN, Paris, France</footnote>

<footnote>*Correspondence to:* T. Penduff
(thierry.penduff@hmg.inpg.fr)</footnote>
<table><caption>**Table 1.** Differences between the four model experiments (see text). The same simulations were investigated by Barnier et al. (2006), but are renamed here for clarity.</caption><thead><tr><th>Simulation name in Barnier et al. (2006)</th><th>Simulation name in the present paper</th><th>Momentum advection scheme</th><th>Topography</th><th>Sidewall boundary conditions</th></tr></thead><tbody><tr><td>ORCA025-G04</td><td>ENSf</td><td>ENS</td><td>Full steps</td><td>Free-slip</td></tr><tr><td>ORCA025-G03</td><td>EENf</td><td>EEN</td><td>Full steps</td><td>Free-slip</td></tr><tr><td>ORCA025-G22</td><td>EENp</td><td>EEN</td><td>Partial steps</td><td>Free-slip</td></tr><tr><td>/</td><td>EENp-ns</td><td>EEN</td><td>Partial steps</td><td>No-slip</td></tr></tbody></table>

enstrophy-energy-conserving momentum advection scheme) largely improved this initial DRAKKAR solution, removing a substantial part of the biases found usually in this class of model solutions and making it comparable to higher resolution standards. These improvements were especially noticed in regions of strong eddy activity and topographic influence, suggesting a beneficial modification of topographic constraints on the eddy and the mean flow.

Barnier et al. (2006) described the successive impact of both numerical changes on DRAKKAR solutions near the surface with respect to observed mean sea-surface height (Niiler et al., 2003) and altimeter-derived Eddy Kinetic Energy (EKE, Ducet et al., 2000). The present paper extends the former assessment toward the subsurface and makes use of an additional simulation to investigate the impact of momentum advection schemes, partial steps, and sidewall boundary conditions on the simulated dynamics and on current-topography interactions. We aim at complementing Barnier et al. (2006)'s study, illustrating the major role of discrete topographic constraints on the simulated dynamics, interpreting the origin of the circulation biases mentioned in CLIPPER-like models and of the improvements seen in DRAKKAR setups. A detailed investigation of certain improvements reported here and in Barnier et al. (2006) is presented from a numerical point of view in Le Sommer et al. (2007). The numerical setup and the four sensitivity experiments are presented in Sect. 2. The sensitivity of the model to changes in the momentum advection scheme and to the use of partial steps in then described, with a focus on circulation features (Sect. 3), on the vertical structure of the flow and topographic constraints (Sect. 4). Section 5 discusses the impact of sidewall boundary conditions and sheds light on the former results. Our results are summarized and discussed in Sect. 6.

## 2 Model, configurations, experiments

The reader is referred to Barnier et al. (2006) for a complete description of the global 1/4° model setup: only the main features are recalled here. This configuration (ORCA025) is part of the DRAKKAR hierarchy of models (DRAKKAR

Group, 2007), developed by a consortium of several research groups in Europe and with the Mercator-Ocean Operational Oceanography Reasearch and Development Team. This global sea-ice ocean model, implemented here on a 46-level 1/4°-resolution tripolar ORCA grid, is based on the NEMO modeling system<fnref n="1" /> which includes the z-coordinate free-surface ocean code OPA9 and the multi-layered sea-ice code LIM2. The simulations presented here use a Flux Corrected Transport (FCT) scheme (Lévy et al., 2001) for tracer advection. Climatological initial conditions for temperature and salinity are taken in January from PHC2.1 (Steele et al., 2001) at high latitudes, MEDATLAS (Jourdan et al., 1998) in the Mediterranean, and Levitus et al. (1998) elsewhere. The 10-year simulations presented here were forced through bulk formulae by a climatological seasonal forcing without any relaxation of sea-surface temperature (SST) or salinity (SSS). These climatological atmospheric fields, close to those used by Timmermann et al. (2005) and labeled “DRAKKAR forcing set 1”, combine wind stresses from ERS and NCEP, daily air temperature from NCEP/NCAR reanalysis, monthly precipitations from CMAP (Xie and Arkin, 1997), air humidity from Trenberth et al. (1989), cloud cover from Berliand and Strokina (1980), and river runoff derived from the UNESCO database (E. Rémy, personal communication).

This paper investigates the model sensitivity to various numerical parameters from four simulations (see Table 1) which share the following parameterizations:

- Laplacian isopycnal tracer diffusion is applied with a coefficient that decreases poleward proportionally to the grid size (300 m² s⁻<fnref n="1" /> at the equator).

- Momentum is dissipated horizontally at small scales by a biharmonic viscosity operator: −1.5×10¹<fnref n="1" /> m⁴ s⁻<fnref n="1" /> at the equator and decreasing poleward as the cube of the grid size. Unresolved turbulent processes exert a lateral viscosity along the equator and are parameterized over the upper 100 m by an additional laplacian operator (500 m² s⁻<fnref n="1" />, see Arhan et al., 2006).

<footnote n="1">¹http://www.lodyc.jussieu.fr/NEMO/</footnote>
- Surface boundary layer mixing and interior vertical mixing are parameterized by a 1.5-order turbulent closure model (Blanke and Delecluse, 1993). Background vertical diffusion is set to $10^{-5} \text{ m}^2 \text{ s}^{-1}$. Background vertical viscosity $K^v$ is set to $10^{-4} \text{ m}^2 \text{ s}^{-1}$. In case of static instability, vertical viscosity and diffusivity are raised to $10 \text{ m}^2 \text{ s}^{-1}$.

- Quadratic bottom friction is introduced as a boundary stress with a simple parameterization of the effect of residual tidal currents (Willebrand et al., 2001): $K^v \frac{\partial u}{\partial z} = C_d \cdot u \sqrt{(u^2 + u_0^2)}$ at the bottom, with $C_d = 10^{-3}$ as in Treguier (1992), and $u_0 = 5 \text{ cm s}^{-1}$.

Various momentum advection schemes (MASs) are available in NEMO-OPA9. In an idealized non-divergent shallow water framework, the horizontal components of the schemes referred to as ENS (Sadourny, 1975) and EEN (Arakawa and Lamb, 1981) conserve potential enstrophy and both energy and potential enstrophy, respectively. The conservation of potential enstrophy is not ensured in the present realistic case. These latter studies show that imposing constraints on the momentum advection scheme truncation error significantly improves the solutions obtained in marginally resolved regimes, e.g. the present eddy-admitting case. In particular, Arakawa and Lamb (1981) show that energy and potential enstrophy conservation suppresses a spurious energy cascade found in the presence of steep topography when only enstrophy is conserved. Both ENS and EEN schemes are written in the vector-invariant form; an extensive description of their formulations and of their dynamical behavior in a realistic setting similar to ours is given in Le Sommer et al. (2007). Simulations ENSf and EENf (see Table 1) will be compared to each other to evaluate the dynamical impact of momentum advection schemes.

Topography is represented as staircases in all runs. It is distorted in the cases when local depths are approximated to the closest model level (i.e. full steps, simulations ENSf and EENf), and much less with partial steps (Pacanowski and Gnanadesikan, 1998) where the thickness of bottom grid cells is adjusted to the real depth (simulations EENp and EENp-ns). Topography representation is the only difference between EENf and EENp that will thus be compared to study the global impact of partial steps. Note that the CLIPPER simulations were performed with the ENS scheme and full step topographies (like the present ENSf run). The present study, along with those by Barnier et al. (2006) and Le Sommer et al. (2007), describe the impact of these new choices in DRAKKAR models, and contributes to explain the improvements between CLIPPER and DRAKKAR solutions.

The first three runs listed in Table 1 were performed with free-slip boundary conditions. The impact of sidewall friction in the eddy-admitting regime will be evaluated by comparing EENp to EENp-ns, this latter run differing only from the former by a no-slip boundary condition (see Table 1).

<figure><img src="image_28.png" /><figcaption>**Fig. 1.** (a) Maximum intensity of the meridional overturning streamfunctions (Sv) computed in z-coordinates over years 8–10 of the four simulations in the northern hemisphere; shown for the upper and lower cells (upper and lower plots) of the global (left column) and Atlantic (right column) basins. (b) Barotropic transports at 7 straits (see map) in the four simulations averaged over years 8–10.</figcaption></figure>

less stated otherwise, the four simulations will be compared in the following in terms of averages and variances of various quantities computed from the last three years (8 to 10) of 10-year integrations.

## 3 Momentum advection and topography: circulation features

The path and structure of simulated horizontal currents respond quickly to modified numerics, i.e. within about 1 year. Modified advection might then make the three simulations diverge in terms of T-S structure; indeed this becomes
<table><caption>**Table 2.** Summary of the model sensitivity in the Atlantic (see text and Figs. 2, 3, 4, 6).</caption><thead><tr><th></th><th>ENSf → EENf<br>ENS → EEN</th><th>EENf → EENp<br>F. steps → P. steps</th><th>EENp → EENp-ns<br>Free-slip → No-slip</th></tr></thead><tbody><tr><td>Transport of<br>Labrador Current</td><td>~40 Sv → ~50 Sv</td><td>~50 Sv → ~60 Sv</td><td>~60 Sv → ~40 Sv</td></tr><tr><td>Northwest Corner</td><td>Appears</td><td>Remains present</td><td>Disappears</td></tr><tr><td>GS: Northern<br>Recirculation<br>Gyre</td><td>~10 Sv → ~20 Sv</td><td>~20 Sv → ~30 Sv</td><td>~30 Sv → ~0 Sv</td></tr><tr><td>GS: Southern<br>Recirculation<br>Gyre</td><td>Very local, ~10 Sv<br>→<br>Local, ~20 Sv</td><td>Local, ~20 Sv<br>→<br>Elongated, ~30 Sv</td><td>Elongated, ~30 Sv<br>→<br>Very local, ~10 Sv</td></tr><tr><td>Deep Western<br>Boundary Current<br>off Brazil</td><td>Weak eddies appear</td><td>Stronger eddies</td><td>Weaker eddies</td></tr><tr><td>Transport of<br>Zapiola Anticyclone</td><td>~0 Sv → ~0 Sv</td><td>~0 Sv → ~200 Sv</td><td>~200 Sv → ~100 Sv</td></tr><tr><td>Northernmost<br>latitude of<br>Malvinas Current</td><td>39° S → 34° S</td><td>34° S → 34° S</td><td>34° S → 39° S</td></tr><tr><td>Behavior of<br>Agulhas Rings</td><td>More dispersion<br>More zonal drift</td><td>Less dispersion</td><td>EENp-ns similar<br>to ENSf</td></tr></tbody></table>

noticeable after about 5 years. On even longer timescales, SST distributions might have diverged enough to affect the surface buoyancy forcing (despite unchanged atmospheric variables) through bulk formulae, and eventually the meridional overturning circulation. Our 10-year integrations are too short to study the ultimate impact of these numerical changes, but well suited to study how various numerical choices establish different (and persistent) circulation patterns. This section is focused on the sensitivity of the flow at different depths to the change of MAS and topographic representation.

## 3.1 Meridional overturning circulation and horizontal transports

Figure 1a shows for the Global and Atlantic Oceans the maximum intensity of the upper and lower cells of the MOC. The magnitude of the upper (lower) cell quantifies the northward transport of warm surface (cold bottom) waters; in the Atlantic for all simulations, these extrema are located around 37° N/1000 m and 15° N/4000 m respectively (Barnier et al., 2006, their Fig. 3). As noted in this latter paper, the upper

mean MOC averaged over the last three years of the simulation is relatively strong in our simulations (23–24 Sv), and the sensitivity of some regional circulation patterns to the model numerics might be related to this regime. Neither the MAS nor the topographic representation substantially modify the maximum intensity of the Atlantic and global upper MOCs: the extrema remain close to [35° N–1000 m] and change by less than 2%. No change of the upper meridional circulation is noticeable in the Indian and Pacific Oceans either. This supports the fact that numerically-induced circulation changes did not significantly feed back to the surface buoyancy fluxes via a large-scale change of SST advection over 10 years. On the contrary, the deep limb of the MOC (associated with the northward spreading and upwelling of Antarctic Bottom Water, AABW) is strongly affected by our numerical changes over this 10-year period. Its maxima, found around [30° S–4000 m] globally and [15/20° N–4000 m] in the Atlantic, substantially increase in the ENSf-EENf-EENp sequence: respectively +8/+28% (global/Atlantic) when ENS is replaced by EEN, and additional +11/+18% when partial steps are used. An increase of the deep MOC is observed in the Indo-Pacific Ocean as well.
Both numerical changes thus induce a direct enhancement of the deep circulation in the Global Ocean.

Figure 1b shows the depth-integrated time-averaged transports at choke points in ENSf, EENf and EENp simulations. These transports increase when the ENS MAS is replaced by EEN (ENSf vs. EENf), and, to a lesser extent, when the full step topography is replaced by partial steps (EENf vs. EENp). The currents across Florida-Bahamas, Florida-Cuba, and Mozambique Straits (labeled a, b, c in Fig. 1b) mainly carry upper-layer waters poleward, part of which contribute to the upper MOC; their enhancement is not seen in the upper MOCs, suggesting that these increments recirculate horizontally.

In summary, highly-confined and deep flows get stronger with the EEN MAS, and even more so with partial steps. We may thus conjecture at this point that both numerical changes reduce the dissipation of horizontal momentum along topographies. This hypothesis is further examined below.

## 3.2 Atlantic circulation features

Since Rhines (1977)'s pioneering paper, many studies have confirmed that mesoscale eddies are produced by and interact with the general circulation, and that these non-linear eddy-mean flow interactions are highly sensitive to sloping topographies. Eddy-active regions such as western boundary currents and their eastward extensions are locations where eddy-mean flow and current-topography interactions are strong, and thus are expected to be particularly sensitive to the choice of the momentum advection scheme and the representation of topography. The impact of the EEN MAS and partial step topography, shown by Barnier et al. (2006) to be beneficial at the surface, is now investigated in more detail at depth in the Atlantic. Table 2 summarizes how both changes affect the model solution in this basin (also see Barnier et al., 2006's Figs. 5, 7 and 8).

### 3.2.1 North Atlantic

Our numerical changes do not modify the transport of Greenland-Scotland overflows, the magnitude, structure and depth of the Atlantic upper MOC, nor the transport of the Deep Western Boundary Current (DWBC, shown at 1500 m in Fig. 2 for all simulations). The structure of the DWBC is, however, modified. Between the Reykjanes Ridge and the tip of the Grand Banks, the DWBC becomes faster and narrower from ENSf to EENf, and from EENf to EENp. Poleward countercurrents are better-defined in EENf and EENp offshore the DWBC in the southwest Labrador Sea, and east of the Grand Banks. In EENf and EENp, most of the DWBC remains trapped along the topographic slope between the Grand Banks and Cape Hatteras (10–15 cm s⁻¹), whereas ENSf simulates a broad vein that remains away from the slope.

<figure><img src="image_29.png" /><figcaption>**Fig. 2.** Current vectors and velocities (colors, m s⁻¹) at 1500 m along the northwestern boundary of the Atlantic, averaged over years 8–10 in the four simulations. Only one vector out of 16 is represented.</figcaption></figure>

Figure 3 shows the vertical structure of the mean flow across the Gulf Stream (GS) and the DWBC at 69°N as measured by Joyce et al. (2005) and diagnosed from the DRAKKAR simulations. The ENSf solution is clearly unrealistic there: the GS is located too far south (around 37°N), the equatorward DWBC is wide and slow, it reaches the surface at latitudes (38°N) where the GS is actually observed, and an unrealistic shallow vein of warm water flows poleward along the American shelf. The use of the EEN scheme (EENf) does not suppress this shallow warm current, but lessens the model-observation mismatch: it confines the shallow current horizontally, the DWBC becomes a more realistic vein trapped along the topographic slope, and an eastward barotropic jet appears at the observed GS location. The addition of partial steps (EENp) does not totally suppress the shallow vein, does not help the DWBC match its observed depth range (still 700 m too shallow) and maximum velocities reach their observed magnitudes in the near-surface GS (50 cm/s instead of 80) and in the DWBC (about 10 cm/s instead of 30). However, the EENp solution highlights major
<figure><img src="image_30.png" /><figcaption>**Fig. 3.** Mean normal velocity across the 69° N section shown on top, derived from current meter measurements (Joyce et al., 2005) and averaged over years 8–10 from the four simulations. Contour intervals: 5 cm s⁻¹ below −10 cm s⁻¹, like color limits between −10 and +20 cm s⁻¹, 10 cm s⁻¹ above +20 cm s⁻¹. The zero contour is highlighted in white and black in the observations and simulations, respectively.</figcaption></figure>

improvements from surface to bottom: the GS becomes a surface-intensified front with realistic structure and location, the intensity of the shallow vein is substantially reduced, and the DWBC remains at a correct location.

Most of the poleward flow within the NAC steers around the quasi-stationary Northwest Corner before heading toward the eastern North Atlantic. This anticyclonic circulation is weak in the ENSf solution and appears at the observed location in EENf and EENp (44° W, 51° N). The Northwest Corner transport reaches only half of Lazier (1994)'s 50 Sv estimate (Figs. 4, Table 2), probably because the resolution remains too modest to favour strong recirculations. The deep signatures of the separated GS (Fig. 2), the shape and intensity of its Northern and Southern Recirculation Gyres (NRG and SRG) also get increasingly realistic in the sequence ENSf-EENf-EENp (Fig. 4, Table 2). Both recirculation gyres in EENp are remarkably close to those described in Dengg et al. (1996)'s review, in terms of location, shape, and internal structure (sub-gyres). These recirculation gyres have realistic transports of about 30 Sv, but two local asymmetries (NRG stronger than SRG, and southward shift of the GS velocity maximum with increasing depth) are not simulated there. Frolov et al. (2004) showed from idealized numerical experiments that in presence of a sloping topography, the observed tilt of the GS velocity maximum in the meridional-vertical plane induces an asymmetric mixing of sub-thermocline potential vorticity mixing that enhances the NRG. The simultaneous absence of both asymmetries in EENp is consistent with such a dynamical link.

A basin-scale section of EKE was done along 48° N by

Colin de Verdière et al. (1989) in the North Atlantic and is superimposed on its model counterparts in Fig. 5. The western part of the section intersects the poleward NAC where the EKE increases in the sequence ENSf-EENf-EENp. The EKE minimum observed right above the MAR, clearly absent from ENSf, is correctly reproduced in EENp. The surface-trapped EKE associated with the NAC extension east the Mid-Atlantic Ridge (MAR) is slightly better reproduced (more intense) in EENp. To sum up, the use of EEN and partial steps largely improve the mean circulation and associated eddy field throughout the northwestern Atlantic up to the surface, where their patterns tend toward those obtained at much higher resolution (see Barnier et al., 2006).

### 3.2.2 South Atlantic

ENSf, EENf and EENp simulate quite comparable DWBC's south of the GS region until about 8° N where the observed DWBC breaks up into coherent eddies (Dengler et al., 2004). Successful simulations of this transition have required horizontal resolutions of at least 1/6° so far, like in CLIPPER, up to 1/12° in the FLAME model (Dengler et al., 2004). Fig. 6 shows that the DWBC remains laminar with the ENS scheme, but generates eddies with the EEN scheme. Despite the modest resolution, the addition of partial steps improves the solution and brings EENp's eddies close to those simulated by the 1/12° FLAME model in terms of coherence, temperature anomaly (0.2°C), quasi-circular shape, mutual spacing, and size.
Mean barotropic streamfunction in the Atlantic

<figure><img src="image_31.png" /><figcaption>**Fig. 4.** Mean barotropic streamfunction in the Gulf Stream region (left) and in the Confluence region (right) from simulations ENSf, EENf, EENp, EENp-ns (from top to bottom). Anticlockwise circulations are shaded. Contour interval is 10 Sv (20 Sv) in the left (right) columns.</figcaption></figure>

The right panels in Fig. 4 (also see Table 2 here, and Fig. 9 in Barnier et al., 2006) highlight in the Confluence region the two-step improvement of the model solution in the sequence ENSf-EENf-EENp. The EEN advection scheme enhances by 5° the northward excursion of the Malvinas Current (which reaches about 34° S, consistently with Olson et al., 1988) and intensifies stationary recirculations between 50 and 55° W. The use of partial steps further improves the solution with the appearance around [45° W, 45° S] of the eddy-driven topographically-trapped Zapiola anticyclone (De Miranda et al., 1999). Its simulated shape, mean transport (~200 Sv), surface signature and impact on the EKE field (see Barnier et al., 2006) get in very good agreement with observations. Eddy-mean flow interactions with topography play a major role in other regions such as the Cape basin (impact of the Walvis Ridge on the Agulhas Rings,

see Kamenkovich et al., 1996; Beismann et al., 1999; Penduff et al., 2002). The combined use of the EEN advection scheme and partial steps led to substantial improvements there as well (Table 2; Barnier et al., 2006).

# 4 Momentum advection and topography: vertical structure of the flow

Our changes in numerical schemes modify the vertical structure of the mean and eddy flow. These changes are described hereafter, first in terms of orientation of the mean circulation with respect to topography, then in terms of intensity of horizontal and vertical (mean and eddy) velocities.
<figure><img src="image_32.png" /><figcaption>Fig. 5. Eddy kinetic energy section in the Atlantic across 48° (cm² s⁻²) from the four DRAKKAR simulations (color) and from the current-meter measurements by Colin de Verdière et al. (1989) (contours).</figcaption></figure>

## 4.1 Direction of the mean flow: topostrophy

Eddy flows interact with the ambient barotropic potential vorticity field $\frac{f}{H}$ whose gradient is influenced by sloping topographies: cross-$\frac{f}{H}$ eddy flows basically induce a downgradient, turbulent flux of barotropic potential vorticity. This nonlinear process is able to generate mean barotropic flows, such as the Zapiola Anticyclone (Dewar, 1998), which keep deep waters to the left (right) in the northern (southern) hemisphere, i.e. in the direction of $\mathbf{u}^{\text{topo}}=-\mathbf{f} \times \nabla H$. Statistical theories also predict the existence of such mean barotropic flows driven by eddy-topography interactions (Holloway, 1992). Topostrophy (Merryfield and Scott, 2007) is a statistical measure of the alignment at every cell $m$ of mean normalized currents vectors $\mathbf{u}_m$ ($|\mathbf{u}_m|=1, \forall m$) with the barotropic circulation $\mathbf{u}_m^{\text{topo}}$. Topostrophy $T_s(x, y)$ was computed globally from years 8–10 of the four DRAKKAR simulations and binned in classes of total depth $H_j$ and level depth $d_k$ to yield $T_s(H_j, d_k)$ as done by Merryfield and Scott (2007) (unsmoothed velocity and topography, weighting $=|\nabla H| \cdot |f|$ according to their notation):

$$T_s(H_j, d_k) = \frac{\sum_m (\mathbf{u}_m \cdot \mathbf{u}_m^{\text{topo}}) \cdot \delta_{jk,m} \cdot dV_m}{\sum_m |f \nabla H|_m \cdot \delta_{jk,m} \cdot dV_m}, \quad (1)$$

where $dV_m$ is the grid cell volume, $\delta_{jk,m}$ equals 1 where the bottom depth is $H_j$ and the level depth is $d_k$, 0 otherwise.

For each $(H, d)$ class, $T_s=1, T_s=-1$, and $T_s=0$ respectively correspond to simulated flows everywhere in the direction of $\mathbf{u}^{\text{topo}}$, opposite to $\mathbf{u}^{\text{topo}}$, and uncorrelated with $\mathbf{u}^{\text{topo}}$. The closer $T_s$ to unity, the more the simulated currents circulate in the direction mentioned above, and the better the model agrees with topographic rectification theories. Merryfield and Scott (2007) examined topostrophies computed from global simulations performed at various resolutions (1°, 2°, 0.4°, and 0.1°). The large difference found between coarse and fine resolution topostrophies supports the aforementioned theory, i.e. mesoscale turbulence tends to align mean currents in the direction of $\mathbf{u}^{\text{topo}}$, with values of $T_s$ increasing toward the bottom in all total depth classes.

Global topostrophies $T_s(H, d)$ are shown in Fig. 7 for the four DRAKKAR simulations. Within most depth ranges in ENSf, EENf and EENp, $T_s$ is positive and monotonically increases toward the bottom: simulated currents increasingly tend to align in the direction of topographic Rossby waves with increasing depth. $T_s$ in EENf is about twice its typical value in ENSf, over a 1000m-thick bottom layer (especially within the 0–2000 m bottom depth range). The EEN advection scheme substantially enhances $T_s$ near the surface as well. The effect of partial steps is qualitatively comparable but weaker: this further enhances deep values of $T_s$ in the 2000–4000 m bottom depth range. As a result, $T_s$ in EENp becomes rather uniform over the bottommost 1000 m
<figure><img src="image_33.png" /><figcaption>**Fig. 6.** Deep Western Boundary Current temperature and velocity snapshots in the four DRAKKAR simulations at 2500 m (22 July of year 10) and in a 1/12° FLAME simulation at 1900 m (Dengler et al., 2004).</figcaption></figure>

within the 0–3500 m bottom depth range. As compared to Merryfield and Scott (2007)'s Fig. 7, the $T_3$ pattern in ENSf shares some similarities (weak values becoming negative in deep regions) with those derived from non-eddyking models; our numerical changes, especially the EEN scheme, make the topostrophy structure more comparable to those derived from 0.1°models (stronger, positive almost everywhere, homogeneous within a bottommost 1000 m-thick layer). As discussed by Merryfield and Scott (2007), results of this kind do not demonstrate that eddy-topography interactions in ENSf and EENp are responsible for the changes in $T_3$, but they are strongly suggestive of such an effect.

## 4.2 Vertical structure of horizontal mean and eddy velocities

The temporal average $u$ and standard deviation $u'$ of horizontal velocity magnitudes were computed at every gridpoint over years 8–10 of each run and averaged quasi-zonally on the model grid ($j$=constant lines) to build meridional sections; these averages are noted $|u|$ and $|u'|$ respectively. Their relative differences between ENSf and EENf and between EENf and EENp are expressed in % in Fig. 8. Replacing the ENS scheme by the EEN scheme leads to a clear enhancement of mean deep velocities, especially at mid-latitudes (roughly 30–60°) in both hemispheres. The change of momentum advection scheme increases sub-thermocline mean

velocities by about 10 to 100% at the bottom in these latitude bands (+20% up to 45% on global average). This sensitivity, investigated in detail by Le Sommer et al. (2007), is particularly clear in eddy-active regions influenced by bottom topography. A similar enhancement occurs on subthermocline eddy velocities $|u'|$ as well, but to a lesser extent (about +10% in slightly narrower latitude bands). The right panels in Fig. 8 show that the addition of partial steps reinforces the deep mean and eddy flows, in a similar but stronger and more homogeneous way: everywhere north of 60° S including low latitudes, subthermocline mean (eddy) velocities increase by 10–50 (10–30)% up to 100% near the bottom.

Most circulation biases of the ENSf solution (mislocated currents, no Zapiola anticyclone, overshoot and/or spreading of western boundary currents, straight path of Agulhas eddies as mentioned in Table 2, etc.) were found in the 1/3° and 1/6° CLIPPER model solutions; they are in fact typical of many eddy-admitting and eddy-resolving ocean models (see Barnier et al., 2006). Both CLIPPER models cited above share the ENSf numerics (ENS, full steps); mean and eddy velocities were proven much too weak at depth in CLIPPER with respect to WOCE current meter data (Penduff et al., 2002, 2005). The sub-thermocline velocity increases induced by the EEN scheme and further by the partial steps are thus beneficial.
<figure><img src="image_34.png" /><figcaption>**Fig. 7.** Left: global topostrophy *T* (see text) from simulations ENSf, EENf, EENp and EENp-ns. Right: differences $T_{EENf} - T_{ENSf}$ (impact of momentum advection scheme) and $T_{EENp} - T_{EENf}$ (impact of partial steps).</figcaption></figure>

With the exception of the bottom left panel, Fig. 8 shows that the enhancement of $|u|$ and $|u'|$ due to EEN and to partial steps strongly increases with depth. It is interesting to note also that at the surface, where Barnier et al. (2006) compared these results with observations, the modifications in mean and eddy velocities induced by both numerical changes are actually aligned with those occuring below the thermocline, and much weaker than their near-bottom expressions. Eddy velocities also tend to change in the depth and latitude ranges where mean velocities do so. All these remarks strongly suggest, again, that the realism of surface mean and eddy flows in EENp are a consequence of more intense currents (and more consistent dynamics, see previous section) at depth. This explanation is consistent with a numerically-induced modification of mutual interactions between the mean flow, the eddy flow, and topography.

## 4.3 Vertical structure of vertical mean and eddy velocities

Figure 9a shows the temporal average and standard deviation of vertical velocities just above topography in EENf and EENp in the Cape Basin (the following remains true

in other eddy-active regions). The use of partial steps reduces the extrema of both statistical moments, which appear much smoother at the bottom. The mean ($|w|$) and standard deviation ($|w'|$) of vertical velocities were averaged along the model grid ($j$=constant lines) as $|u|$ and $|u'|$ introduced earlier, and are shown in Fig. 9b. Between about 500 and 4000 m at mid latitudes, the use of partial steps enhances $|w|$ and $|w'|$ by about 10–20%; these increases are collocated with their counterparts in $|u|$ and $|u'|$, respectively. Below 4000 m, partial steps enhance horizontal velocities but substantially reduce mean and eddy vertical velocities (–30% and –15% respectively) on global average at all latitudes. Partial steps thus keep the vertical-to-horizontal “aspect ratios” of mean and eddy velocities unchanged above about 4000 m, but strongly reduce them near the bottom (by 50 and 30%, respectively for $|w|/|u|$ and $|w'|/|u'|$).

Topographic slopes, especially weak ones, are much better represented with partial steps: this suppresses localized straircases where bottom $w$ and $w'$ reach their maxima in the full step case (Fig. 9a). This benefit was mentioned by Pacanowski and Gnanadesikan (1998) (their Fig. 7) and is confirmed here at higher resolution in a global realistic setup. By reducing velocity gradients at grid-scale, partial steps are expected to reduce the explicit dissipation of momentum (and vorticity) and may promote stronger velocities at depth as found in EENp. By improving the representation of $f_H$ contours, partial steps are also expected to yield more consistent vorticity dynamics, as discussed in section 6. Note that $|w|$ and $|w'|$ are not significantly modified by the change of momentum advection schemes (not shown).

To sum up at this point, the use of EEN and partial steps improves the model solution at all depths, enhances subsurface mean and eddy flows beneficially, and makes the orientation of the mean flow more physically-consistent with topographic constraints throughout the water column. Eddy-admitting models with reasonable numerics thus seem able to simulate part of the mean circulation driven by eddy-topography interactions, as expected from statistical theories and barotropic potential vorticity arguments. The multi-resolution DRAKKAR model hierarchy would be adequate to further study the impact on these mechanisms of progressive changes of model resolution, as initiated by Merryfield and Scott (2007) from coarse and eddy-resolving solutions.

# 5 Explicit and spurious sidewall friction

A free-slip boundary condition is used in all simulations, except in EENp-ns which differs from EENp only by a no-slip condition. Figures 1, 2, 3, 4, 5, 6, 7 and Table 2 exhibit many changes between EENp and EENp-ns. The barotropic transport at Denmark Straits and the upper global MOC are stronger in EENp than in ENSf and EENf, and further increase by a few percents in EENp-ns. These features are not fully understood yet, but none is due to a stronger Atlantic
<figure><img src="image_35.png" /><figcaption>**Fig. 8.** Relative differences (%) of $|u|$ (top) and $|u'|$ (bottom) between ENSf and EENf (100×$\frac{\text{EENf}^-\text{ENSf}}{\text{ENSf}}$, left) and between EENf and EENp (100×$\frac{\text{EENp}^-\text{EENf}}{\text{EENf}}$, right). $|u|$ and $|u'|$ are the quasi-zonal averages of the time-mean and standard deviation of current intensities, respectively, computed over years 8–10 from each simulation.</figcaption></figure>

upper MOC (which slightly decreases); there is a slight enhancement of the cyclonic circulation around Iceland, and a possible transfert of the overturning circulation at global scale toward the Pacific.

Besides these exceptions, the main impact of the no-slip condition is to bring the EENp solution back towards ENSf or EENf with respect to all remaining features. Indeed, the no-slip condition strongly reduces topostrophy in all local and bottom depth classes (i.e. decorrelates the direction of currents from that of topographic Rossby waves, Fig. 7), weakens the mean and eddy horizontal velocities at depth at all latitudes (not shown), adversely reduces the equatorward excursion of subpolar boundary currents (NRG and Malvinas Current) by about 5°, exagerates the poleward overshoot of subtropical western boundary currents (GS and Brazil Current), widens along-slope currents north of about 50°N as well as the detached GS, makes the path of Agulhas Rings unrealistically straight as in ENSf (see Barnier et al., 2006), strongly reduces the Zapiola Anticyclone transport (see Table 2 and Fig. 4), weakens the Labrador Current, the NAC east of the Grand Banks, most barotropic transports (in Fig. 1b), the DWBC (Figs. 2, 3), and accordingly, the deep overturning circulation (Fig. 1a).

Interestingly, the no-slip condition introduced in EENp-ns cancels out most improvements due to the use of partial steps and/or the EEN scheme. In other words, the sidewall friction brings back the biases diagnosed in ENSf or EENf, most of which were also mentioned in CLIPPER (Penduff et al., 2005) and seen as well in many models at comparable or higher resolutions (Barnier et al., 2006). This strongly support the hypothesis made previously: not using the EEN scheme and partial steps distorts current-

topography interactions and the model solution in free-slip cases (ENSf, EENf, and ENSf-like CLIPPER runs), just as if a numerically-induced spurious sidewall friction was actually present. Adcroft and Marshall (1998) demonstrated the existence of such a spurious sidewall friction within a 1/4°-resolution, C-grid shallow-water idealised model implemented in free-slip with an enstrophy-conserving momentum advection scheme (same numerics as ENSf). Such an effect had been hypothized to explain the absence of the Zapiola anticyclone (De Miranda et al., 1999; Penduff et al., 2001), underestimated sub-thermocline kinetic energy levels and general circulation biases met in CLIPPER with similar numerics (Penduff et al., 2002, 2005) and in other eddy-admitting z-level simulations (Barnier et al., 2006).

Our results thus confirm the existence of spurious sidewall friction in CLIPPER-like models (and presumably other eddy-admitting models), illustrate its adverse impacts on near-bottom flows, along-slope currents at all vertical levels and subsequent realism of numerical solutions up to the surface, and show that these problems are substantially reduced by the use of EEN and partial steps.

The 1/4° OCCAM global model (Coward and de Cuevas, 2005) shares the same resolution and the use of partial steps with EENp. However OCCAM is based on the B-grid MOM3 code (DRAKKAR is based on the C-grid NEMO code), uses an explicit no-slip sidewall friction (vs. free-slip in EENp), a momentum advection scheme in flux form (vs. vector-invariant), and a different pressure gradient ($\nabla P$) scheme. Despite the use of partial steps, OCCAM compares much better with ENSf, EENf or EENp-ns, than with EENp (Barnier et al., 2006). In particular, OCCAM does not simulate the Zapiola Anticyclone. A sensitivity experiment
<figure><img src="image_36.png" /><figcaption>**Fig. 9.** (a) Perspective views of bottom vertical velocity statistics ($10^{-4}$ m s$^{-1}$) in the Cape Basin from years 8–10. Temporal average of $w$ ($|w|$, left) and standard deviation of $w$ ($|w'|$, right) in simulation EENf (top) and EENp (bottom). (b) Relative differences (%) of $|w|$ (top) and $|w'|$ (bottom) between EENf and EENp ($100 \times \frac{\text{EENp} - \text{EENf}}{\text{EENf}}$) computed from years 8–10.</figcaption></figure>

was performed in DRAKKAR, where the EENp model was run with the same $\nabla P$ scheme as OCCAM's. This did not change the solution, nor affected the Zapiola Anticyclone. The explicit no-slip boundary condition used in the OCCAM run (and in most MOM-based simulations) thus likely explains many discrepancies compared to EENp, with a possible

contribution of the momentum advection scheme written in flux-form (shown to degrade the DRAKKAR solution, as shown in Le Sommer et al., 2007).

The consequences of the spurious “no-slip-like” sidewall friction hypothesized in CLIPPER at 1/6° were much less visible at low latitudes, where local Rossby radii were
properly resolved (Penduff et al., 2005). The OFES global model (Masumoto et al., 2004) is based on MOM, uses partial steps with a no-slip sidewall friction like OCCAM, but simulates a reasonable Zapiola Anticyclone at 1/10° resolution. Both facts suggest that increased resolution reduces the adverse effects of no-slip conditions, either introduced explicitly (e.g. B-grid models like OCCAM or OFES), or spuriously (CLIPPER, ENSf, EENf). This strongly supports Adcroft and Marshall (1998)'s conjecture: "in the limit of both high resolution and very low viscosity, [...] the impact of [spurious form stresses along piecewise-constant coastlines] on the structure of the circulation is likely to be reduced".

## 6 Conclusion and discussion

Four 10-year 1/4° global simulations have been analyzed to estimate the impact of numerical choices on DRAKKAR solutions, and to propose dynamical interpretations. Our analysis extends the study by Barnier et al. (2006) beneath the surface and in terms of current-topography interactions. Starting from CLIPPER-like numerics, two successive changes were shown to substantially improve the model solution: replacing the enstrophy-conserving (ENS) by an energy-enstrophy-conserving (EEN) momentum advection scheme, and the full steps by a partial steps topography. The sensitivity of the model solution to no-slip and free-slip boundary conditions helped identify the dynamical origin of these improvements. Interestingly, the use of either ENS instead of EEN, full steps instead of partial steps, or no-slip instead of free-slip conditions (spurious in CLIPPER and ENSp, see Adcroft and Marshall, 1998, and explicit in EENp-ns) has comparable and adverse consequences:

- it weakens deep overturning cells (northward spreading of AABW), horizontal transports at choke points, and boundary currents (Greenland and Labrador Currents, DWBC, GS, GS's northern recirculation gyre, Zapiola anticyclone, etc.);

- it intensifies the bias diagnosed in CLIPPER (Penduff et al., 2002, 2005), i.e. reduces the eddy and mean kinetic energy levels at depth;

- it adversely affects the realism (path, width, structure, unstable character) of the DWBC and mean currents, especially in eddy-active, topographically-influenced regions (GS and both its recirculation gyres, NAC, Northwest Corner, Brazil and Malvinas Currents, Confluence Region, Agulhas Rings);

- it decreases topostrophy, especially at depth. This reveals a weaker tendency of mean currents up to the surface to align in the direction of topographic Rossby waves, consistently with weaker eddy-topography interactions. This suggests that eddy-admitting models with

adequate schemes can represent topographic rectification processes, at least partly.

Our results indicate that, like partial steps and sidewall boundary conditions, the influence of momentum advection schemes on the simulated dynamics originate near topography (confirmed in Le Sommer et al., 2007) and yield significantly different solutions up to the surface. In other words, the upper-layer circulation improvements reported by Barnier et al. (2006) are surface signatures of near-bottom numerically-induced changes. Adding an explicit sidewall friction, i.e. no-slip boundary condition, to our "best" EENp configuration actually cancels most benefits of the EEN scheme and partial steps, bringing the solution (EENp-ns) back towards CLIPPER-like standards (ENSf) in all aspects. This shows that many discrepancies found in CLIPPER and DRAKKAR simulations performed with enstrophy-conserving (ENS-like) momentum advection schemes come from a spurious sidewall friction inherent to free-slip conditions (Adcroft and Marshall, 1998). Le Sommer et al. (2007) conclude that compared to the EEN scheme, ENS generates substantially more grid-scale velocity shears along topographies, likely to yield more dissipation of deep momentum through biharmonic viscosity, thus to explain weaker kinetic energy at depth.

The use of partial steps decreases the near-bottom intensity of $|w|$ and $|w'|$ (Fig. 9b) by regularizing topographic slopes (see Fig. 3 in Merryfield and Scott, 2007). However, unlike all other quantities under investigation, neither $|w|$ nor $|w'|$ levels are substantially affected near the bottom by the use of EEN instead of ENS, or by the use of no-slip instead of free-slip. It is thus likely that the sensitivity of the mean and eddy flow to these numerical choices is governed more by their rotational part than their divergent part, which directly determines vertical velocities. If we assume that (both explicit and spurious) sidewall friction has comparable effects on deep flows as mesoscale topographic roughness (see the next paragraph), we may conjecture that switching from ENS to EEN, no-slip to free-slip, or full steps to partial steps improves the representation of vorticity dynamics near topographic slopes, and that this effect reaches up to the surface. Mesoscale topographic roughness was shown in quasigeostrophic (QG) models to remove eddy kinetic energy (EKE) from the barotropic mode, i.e. at depth (Treguier and Hua, 1988; Barnier and Le Provost, 1993), by modifying vertical mode interactions and inhibiting the inverse cascade. Switching from the EEN to the ENS scheme, from partial to full steps, and from free- to no-slip conditions interestingly yield very similar effects on DRAKKAR EKE fields, especially in eddy-active topographically-constrained regions. The use of the ENS scheme, full steps, and no-slip conditions reduce the topographic imprint on the general circulation up to the surface, weaken eddy-driven mean features like the Zapiola anticyclone, and topostrophy at all depths. These full-depth changes in the mean circulation might be
driven by both Reynolds stresses and topographic rectification. The first process matches the conjecture made in the companion numerical study by Le Sommer et al. (2007): "more vigourous bottom-intensified eddy motions [...] can interact with the surface-intensified eddy field, resulting in a modification of the mean flow forcing by the surface eddy field" when switching from ENS to EEN. The improved representation of the topographic potential vorticity field ($\frac{\partial \phi}{\partial H}$), provided by replacing full steps by partial steps, is likely to promote more consistent topographic eddy rectification, and thus the contribution of the second process mentioned above. More generally, quasi-geostrophic vorticity dynamics and eddy rectification theories seem consistent with most of our findings.

Like mesoscale topographic roughness (Treguier and Hua, 1988; Barnier and Le Provost, 1993), bottom friction was shown to decrease subsurface EKE levels in QG and primitive equation models (Rivière et al., 2004), inhibit the inverse cascade and thus weaken the barotropic flow (Panetta, 1993). The smoothing of mesoscale topographic roughness in an eddy-admitting CLIPPER (i.e. ENSf-like) model, thus subject to spurious sidewall friction, was also shown (Penduff et al., 2002) to enhance the penetration of EKE at deep levels and the topographic imprint on the mean circulation (e.g. the Zapiola anticyclone). In other words, the increase of topographic roughness in CLIPPER had comparable effects on the flow as introducing either (i) spurious sidewall friction by replacing EEN by ENS in DRAKKAR, (ii) explicit sidewall friction via a no-slip condition in DRAKKAR, (iii) or topographic roughness in a QG model. The assumption made above about the comparable effect of sidewall friction and topographic roughness, and our QG (rotational) dynamical hypotheses are thus well supported by previous studies.

A substantial topographic smoothing had been shown by Penduff et al. (2002) to reduce the dynamical biases diagnosed in eddy-admitting z-level models like CLIPPER and the present ENSf setup. The present results show that the spurious near-topography no-slip-like friction highlighted by Adcroft and Marshall (1998) and subsequent circulation biases (also typical of other models) can be largely reduced without removing high-wavenumber bathymetric structures, i.e. by the use of adequate momentum advection schemes and partial steps. EENp numerics are being used in present DRAKKAR simulations and are advised for other eddy-admitting simulations with C-grid z-level models. These spurious effects might not be as problematic at higher resolution; both numerical choices yield realistic results in finer DRAKKAR setups ($\frac{1}{12}°$, see DRAKKAR Group, 2007).

It is interesting to note that this global z-level EENp solution shares substantial similarities with regional sigma-coordinate model solutions (see e.g. Ezer and Mellor, 1997; Penduff et al., 2001; Willebrand et al., 2001), in particular enhanced kinetic energies at depth, a realistic Zapiola anticyclone and an improved Gulf Stream path compared to full-step z-level solutions at similar resolutions. Com-

pared to full steps, sigma coordinates and partial steps do not only provide a more accurate grid-scale representation of gentle topographic slopes, associated $\frac{\partial \phi}{\partial H}$ fields and ambient potential vorticity gradients. They also allow a direct communication (through advection, diffusion, pressure gradient computations, etc) between adjacent cells along sloping topographies. Both partial-step and sigma-coordinate models thus yield a more continuous representation of bottom Ekman transport and pumping, known to be largely invloved in eddy-topography interactions (see e.g. Dewar, 1998) and topographic steering of mean currents (see e.g. Arhan et al., 1989). Pressure gradients are everywhere computed on the horizontal in NEMO, including with partial steps above topographic slopes (computation at the shallowest of adjacent bottom T-points) without any vertical extrapolation. No sigma-like pressure gradient truncation error nor hydrostatic inconsistency may thus contaminate discrete pressure gradient computations in our partial steps simulations.

Current-topography and eddy-topography interactions are crucially involved in the systematic improvement seen up to the surface in the sequence of experiments ENSf-EENf-EENp. This suggests that topography also plays a major role in shaping real ocean circulations, locating major surface currents at their observed position and thus the regions where major air-sea interactions occur. Several studies have shown that the explicit resolution of ocean eddies improves the consistency of processes involved in ocean-atmosphere and physico-biogeochemical interactions: besides their significant impact on the variability of oceanic heat transport (Hall et al., 2004) and ocean-atmosphere coupled variability modes (Hogg et al., 2006), ocean eddies are directly involved in upper-ocean processes, e.g. mixed-layer dynamics (Levy et al., 2005), subduction (Valdivieso da Costa et al., 2005) or post-convective restratification (Chanut et al., 2007). Our results emphasize the importance of physically-consistent representations of deep currents and topographic constraints as well, since they shape the circulation and eddy activity at the surface where air-sea interactions ultimately take place.

**Acknowledgements.** The authors would like to thank L. Debreu and G. Holloway for stimulating discussions, T. Ezer and an anonymous reviewer whose comments helped clarify the present article. This work is a contribution of the DRAKKAR project. Support to DRAKKAR comes from various grants and programs listed hereafter: French national programs GMMC, PATOM, and PNEDC; PICS 2475 from Institut National des Sciences de l'Univers (INSU) and Centre National de la Recherche Scientifique (CNRS); Kiel SFB460 and CLIVAR-marin (03F0377A/B) supported by Deutsche Forschungsgemeinschaft. Computations presented in this study were performed at Institut du Développement et des Ressources en Informatique Scientifique (IDRIS). Partial support from the European Commission under Contract SIP3-CT-2003-502885 (MERSEA project) is gratefully acknowledged.

Edited by: J. Schröter
References

Adcroft, A. and Marshall, D.: How slippery are piecewise-constant coastlines in numerical ocean models?, Tellus Series A, 60, 95–108, doi:10.1034/j.1600-0870.1998.00007, 1998.

Arakawa, A. and Lamb, V.: A potential enstrophy and energy conserving scheme for the shallow water equations, Mon. Wea. Rev., 109, 18–36, 1981.

Arhan, M., Colin de Verdire, A., and Mercier, H.: Direct observations of the mean circulation at 48°N in the Atlantic ocean, J. Phys. Oceanogr., 19, 161–181, 1989.

Arhan, M., Treguier, A.-M., Bourlès, B., and Michel, S.: Diagnosing the annual cycle of the equatorial undercurrent in the Atlantic ocean from a general circulation model, J. Phys. Oceanogr., 36, 1502–1522, 2006.

Barnier, B., Madec, G., Penduff, T., Molines, J.-M., Treguier, A.-M., Le Sommer, J., Beckmann, A., Biastoch, A., Böning, C., Dengg, J., Derval, C., Durand, E., Gulev, S., Remy, E., Talandier, C., Theetten, C., Maltrud, M., McClean, J., and de Cuevas, B.: Impact of partial steps and momentum advection schemes in a global circulation model at eddy permitting resolution, Ocean Dynam., 56(5-6), 543–567, doi:10.1007/s10236-006-0082-1, 2006.

Barnier, B. and Le Provost, C.: Influence of bottom topography roughness on the jet and inertial recirculation of a mid-latitude gyre, Dyn. Atm. Oceans, 18, 29–65, 1993.

Beismann, J.-O., Käse, R. H., and Lutjeharms, J. R. E.: On the influence of submarine ridges on translation and stability of Agulhas rings, J. Geophys. Res., 104, 7897–7906, doi:10.1029/1998JC900127, 1999.

Berliand, M. and Strokina, T.: Global distribution of the total amount of clouds, hydrometeorological, Technical report, Hydrometeorology Publishing House, Leningrad, (in Russian), 1980.

Blanke, B. and Delecluse, P.: Variability of the tropical Atlantic ocean simulated by a general circulation model with two different mixed-layer physics, J. Phys. Oceanogr., 23, 1363–1388, 1993.

Chanut, J., Barnier, B., Large, W. G., Debreu, L., Penduff, T., Molines, J.-M., and Mathiot, P.: Mesoscale eddies in the Labrador Sea and their contribution to convection and re-stratification, J. Phys. Oceanogr., in press, 2007.

Chassignet, E., Arango, H., Dietrich, D., Ezer, T., Ghil, M., Haid-vogel, D., Ma, C., Mehra, A., Paiva, A., and Sirkes, Z.: DAMEE-NAB: The base experiments, Dyn. Atmos. Oceans, 32, 155–183, 2000.

Colin de Verrière, A., Mercier, H., and Arhan, M.: Mesoscale variability from the western to the eastern Atlantic along 48°N, J. Phys. Oceanogr., 19, 1149–1170, 1989.

Coward, A. and de Cuevas, B.: The OCCAM 66-level model: model description, physics, initial conditions, and external forcing, Internal document 99, Southampton Oceanography Center, Southampton, UK, 2005.

De Miranda, A., Barnier, B., and Dewar, W.: On the dynamics of the Zapiola anticyclone, J. Geophys. Res., 104, 21137–21149, 1999.

Dengg, J., Beckmann, A., and Gerdes, R.: The Gulf Stream separation problem, *The Warmwatersphere of the North Atlantic Ocean*, W. Krauss, ed., Gebrüder Borntraeger, Berlin, 253–290, 1996.

Dengler, M., Schott, F., Eden, C., Brandt, P., Fisher, J., and Zantopp, R.: Break-up of the Atlantic Deep Western Boundary Current into eddies at 8°S, Nature, 432, 1018–1020, doi:10.1038/nature03134, 2004.

Dewar, W.: Topographic and barotropic transport control by bottom friction, J. Mar. Res., 56, 295–328, 1998.

Ducet, N., Le Traon, P.-Y., and Reverdin, G.: Global high-resolution mapping of ocean circulation from TOPEX/Poseidon and ERS-1 and -2, J. Geophys. Res., 105(C8), 19477–19498, 2000.

DRAKKAR Group: Eddy-permitting ocean circulation hindcasts of past decades, Clivar Exchanges, No 42 (vol 12 No 3), 8–10 (Unpublished manuscript), 2007.

DYNAMO Group: Dynamo: Dynamics of North Atlantic models. Simulation and assimilation with high resolution models, Technical report, Ber. Inst. f. Meeresk. Kiel, 294, 1997.

Eden, C. and Böning, C.: Sources of eddy kinetic energy in the Labrador Sea, J. Phys. Oceanogr., 32, 3346–3363, 2002.

Ezer, T. and Mellor, G.L.: Simulations of the Atlantic Ocean with a free surface sigma-coordinate ocean model, J. Geophys. Res., 102, C7, 15,647–15,657, 1997.

Frolov, S. A., Sutyrin, G. G., and Ginis, I.: Asymmetry of an equilibrated Gulf Stream-type jet over topographic slope, J. Phys. Oceanogr., 34, 1087–1102, 2004.

Griffies, S. M., Böning, C., Bryan, F.O., Chassignet, E.P., Gerdes, R., Hasumi, H., Hirst, A., Treguier, A.-M., and Webb, D.: Developments in ocean climate modelling, Ocean Modelling, 2, 123–192, 2000.

Hall, N., Barnier, B., Penduff, T., and Molines, J.-M.: Interannual variation of Gulf Stream heat transport in a high resolution model forced by reanalysis data, Clim. Dyn., 23, 341–351, 2004.

Hogg, A., Dewar, W., Killworth, P., and Blundell, J.: Decadal variability of the mid-latitude climate system driven by the ocean circulation, J. Climate, 19, 1149–1166, 2006.

Holloway, G.: Representing topographic stress for large-scale ocean models, J. Phys. Oceanogr., 22, 1033–1046, 1992.

Jourdan, D., Balopoulos, E., Garcia-Fernandez, M., and Maillard, C.: Objective analysis of temperature and salinity historical data set over the mediterranean basin, Technical report, IEEE, 1998.

Joyce, T., Dunworth-Baker, J., Pickart, R., Torres, D., and Waterman, S.: On the Deep Western Boundary Current south of Cape Cod, Deep-Sea Res. II, 52, 615–625, 2005.

Kamenkovich, V. M., Leonov, Y. P., Nechaev, D. A., Byrne, D. A., and Gordon, A. L.: On the influence of bottom topography on the Agulhas eddy, J. Phys. Oceanogr., 26, 892–912, 1996.

Lazier, J. R. N.: Observations in the northwest corner of the North Atlantic Current, J. Phys. Oceanogr., 24, 1449–1463, 1994.

Le Sommer, J., Penduff, T., Theetten, S., Madec, G., and Barnier, B.: How momentum advection schemes influence current-topography interactions at eddy-permitting resolution, Ocean Modelling, in press, 2007.

Levitus, S., Boyer, T., Conkright, M., O'Brian, T., Antonov, J., Stephens, C., Stathopolos, L., Johnson, D., and Gelfeld, R.: World Ocean database 1998, Technical Report NESDID18, NOAA Atlas, 1998.

Lévy, M., Estubier, A., and Madec, G.: Choice of an advection scheme for biogeochemical models, Geophys. Res. Lett., 28, 3725–3728, 2001.

Levy, M., Gavart, M., Memory, L., Caniaux, G., and Paci, A.: A
4D-mesoscale map of the spring bloom in the Northeast Atlantic (POMME experiment): results of a prognostic model, J. Geophys. Res., 110, C07S21, doi:10.1029/2004JC002588, 2005.

Masumoto, Y., Sasaki, H., Kagimoto, T., Komori, N., Ishida, A., Sasai, Y., Miyama, T., Motoi, T., Mitsudera, H., Takahashi, K., Sakuma, H., and Yamagata, T.: A fifty-year eddy-resolving simulation of the World Ocean: Preliminary outcomes of OFES (OGCM for the Earth Simulator). *J. Earth Simulator*, 1, 35–56, 2004.

Merryfield, W. and Scott, R.: Bathymetric influence on mean currents in two high-resolution near-global ocean models, Ocean Modelling, 16, 1-2, 76-94, doi: 10.1016/j.ocemod.2006.07.005, 2007.

Niiler, P., Maximenko, N., and McWilliams, J.: Dynamically balanced absolute sea dynamically balanced absolute sea level of the Global Ocean derived from near-surface velocity observations, Geophys. Res. Lett., 30 (22), 2164, doi:10.1029/2003GL018628, 2003.

Olson, D., Podesta, G., Evans, R., and Brown, O.: Temporal variations in the separation of Brazil and Malvinas currents, Deep-Sea Res., 35, 1971–1990, 1988.

Pacanowski, R. and Gnanadesikan, A.: Transient response in a z-level ocean model that resolves topography with partial cells, Mon. Wea. Rev., 126, 3248–3270, 1998.

Paiva, A., Hargrove, J., Chassignet, E., and Bleck, R.: Turbulent behavior of a fine mesh (1/12°) numerical simulation of the North Atlantic, J. Mar. Res., 21, 307–320, 1999.

Panetta, R.: Zonal jets in wide baroclinically unstable regions: persistence and scale selection, J. Atmos. Sci., 50, 2073–2106, 1993.

Penduff, T., Barnier, B., Béranger, K., and Verron, J.: Comparison of near-surface mean and eddy flows from two numerical models of the South Atlantic Ocean, J. Geophys. Res., 106, 16, 857–16, 867, 2001.

Penduff, T., Barnier, B., Madec, G., and Molines, J.-M.: On the use of current meter data to assess the realism of ocean model simulations, Ocean Modelling, 11, 399–416, 2005.

Penduff, T., Barnier, B., Verron, J., and Kerbiriou, M.-A.: How topographic smoothing contributes to differentiating the eddy flows simulated by sigma- and z-level models, J. Phys. Oceanogr., 32, 122–137, 2002.

Rhines, P. B.: The dynamics of unsteady currents, in The Sea, Wiley Interscience, NY, vol. 6, 189–318, 1977.

Rivière, P., Treguier, A.-M., and Klein, P.: Effects of bottom friction on nonlinear equilibration of an oceanic baroclinic jet, J. Phys. Oceanogr., 34, 416–432, 2004.

Roberts, M., Marsh, R., New, A., and Wood, R.: An intercomparison of a Bryan-Cox type ocean model and an isopycnic model. Part I: The subpolar gyre and high-latitude processes, J. Phys. Oceanogr., 26, 1495–1527, 1996.

Sadourny, R.: The dynamics of finite-difference models of the shallow-water equations, J. Atmos. Sc., 32, 1975.

Smith, R., Maltrud, M., Bryan, F., and Hecht, M.: Numerical simulation of the North Atlantic Ocean at 1/10°, J. Phys. Oceanogr., 30, 1532–1561, 2000.

Steele, M., Morley, R., and Ermold, W.: PHC: A global ocean hydrography with a high quality Arctic Ocean, J. Climate, 14, 2079–2087, 2001.

Timmermann, A., Goosse, H., Madec, G., Fichefet, T., Ethe, C., and Dulière, V.: On the representation of high latitude processes in the ORCA-LIM global coupled sea-ice ocean model, Ocean Modelling, 8, 175–201, 2005.

Treguier, A. M.: Kinetic energy analysis of an eddy resolving, primitive equation model of the North Atlantic, J. Geophys. Res., 97, No C1, 687–701, 1992.

Treguier, A. M. and Hua, B. L.: Influence of bottom topography on stratified quasi-geostrophic turbulence, Geophys. Astrophys. Fluid Dynamics, 43, 265–305, 1988.

Tréguier, A.-M., Reynaud, T., Pichevin, T., Barnier, B., Molines, J.-M., de Miranda, A., Messager, C., Beismann, J.-O., Madec, G., Grima, N., Imbard, M., and Le Provost, C.: The CLIPPER project: High resolution modelling of the Atlantic, Int. WOCE Newslett, 36, 3–5, 1999.

Treguier, A.-M., Theetten, S., Chassignet, E., Penduff, T., Smith, R., and Talley, L.: The North Atlantic subpolar gyre in four high resolution models, J. Phys. Oceanogr., 35, 757–774, 2005.

Trenberth, K., Olson, J., and Large, W. G.: A global ocean wind stress climatology based on ECMWF analyses, Technical Report NCAR/TN-338+STR, NCAR, 1989.

Valdivieso da Costa, M., Mercier, H., and Tréguier, A.-M.: Effects of the mixed layer time variability on kinematic subduction rate diagnostics, J. Phys. Oceanogr., 35, 427–443, 2005.

Willebrand, J., Barnier, B., Böning, C., Dieterich, C., Killworth, P., Le Provost, C., Jia, Y., Molines, J.-M., and New, A.: Circulation characteristics in three eddy-permitting models of the North Atlantic, Progress in Oceanogr., 48, 123–162, 2001.

Xie, P. and Arkin, P.: Global precipitation: A 17-year monthly analysis based on gauge observations, satellite estimates, and numerical model outputs, Bulletin of the American Meteorological Society, 78, 2539–2558, 1997.
<figure><img src="image_37.png" /></figure>

# Impact of global ocean model resolution on sea-level variability with emphasis on interannual time scales

T. Penduff<fnref n="1" />,<fnref n="2" />, M. Juza<fnref n="1" />, L. Brodeau<fnref n="1" />, G. C. Smith<fnref n="3" />,*, B. Barnier<fnref n="1" />, J.-M. Molines<fnref n="1" />, A.-M. Treguier<fnref n="4" />, and G. Madec<fnref n="5" />,<fnref n="3" />

Received: 1 June 2009 – Published in Ocean Sci. Discuss.: 8 July 2009

Revised: 17 February 2010 – Accepted: 22 February 2010 – Published: 25 February 2010

**Abstract.** Four global ocean/sea-ice simulations driven by the same realistic 47-year daily atmospheric forcing were performed by the DRAKKAR group at 2°, 1°, ½°, and ¼° resolutions. Simulated mean sea-surface heights (MSSH) and sea-level anomalies (SLA) are collocated over the period 1993–2004 onto the AVISO dataset. MSSH fields are compared with an inverse estimate. SLA datasets are filtered and compared over various time and space scales with AVISO regarding three characteristics: SLA standard deviations, spatial correlations between SLA variability maps, and temporal correlations between observed and simulated band-passed filtered local SLA timeseries. Beyond the 2°–1° transition whose benefits are moderate, further increases in resolution and associated changes in subgrid scale parameterizations simultaneously induce (i) strong increases in SLA standard deviations, (ii) strong improvements in the spatial distribution of SLA variability, and (iii) slight decreases in temporal correlations between observed and simulation SLA timeseries. These 3 effects are not only clear on mesoscale (14–180 days) and quasi-annual (5–18 months) fluctuations, but also on the slower (interannual), large-scale variability ultimately involved in ocean-atmosphere coupled processes. Most SLA characteristics are monotonically affected by successive resolution increases, but irregularly and with a strong dependence on frequency and latitude. Benefits of enhanced resolution are greatest in the 1°–½° and ½°–¼° transitions, in the 14–180 day range, and within eddy-active mid- and high-latitude regions. In the real ocean, most eddy-active areas are

characterized by a strong SLA variability at all timescales considered here; this localized, broad-banded temporal variability is only captured at ¼° resolution.

## 1 Introduction

The choice of ocean/sea-ice primitive equation model configurations for global climate-oriented (multidecadal or longer) studies generally results from a compromise between the range of time and space scales to be simulated and the available computer resources.

Laminar ocean models (1°-resolution and coarser) do not explicitly resolve mesoscale eddies and fluxes: subgrid-scale diffusive (e.g. laplacian operators) and advective (e.g. Gent and McWilliams, 1990, noted GM90 hereafter) parameterizations are used to mimick certain down-gradient eddy fluxes. Such parameterizations were not designed to mimick up-gradient fluxes though, nor the inverse cascade fed by non-linear interactions at scales close to the first internal deformation radius (e.g. Scott and Arbic, 2007). Resolving wide western boundary currents in coarse-resolution models also requires strong viscosity values, which damp a substantial part of the currents' variability. Because they are computationally-effective, laminar oceans are being used in most global coupled models to address paleoclimatic and prediction issues (i.e. IPCC).

“Eddy-admitting” ocean models (roughly ½° to 10° resolution) resolve mesoscale eddies where the local ratio between the grid scale and dynamically-unstable scales is small enough. Parameterizations are also used at these resolutions

<figure><img src="image_38.png" /></figure>

Correspondence to: T. Penduff
(thierry.penduff@legi.grenoble-inp.fr)

<footnote n="1">¹Laboratoire des Écoulements Géophysiques et Industriels, UMR 5519, CNRS and Université de Grenoble, BP 53, 38041, Grenoble, France</footnote>

<footnote n="2">²Department of Oceanography, The Florida State University, Tallahassee, USA</footnote>

<footnote n="3">³National Oceanography Centre, Southampton, UK</footnote>

<footnote n="4">⁴Laboratoire de Physique des Océans, UMR 6523, CNRS, IFREMER and UBO, 29280 Plouzané, France</footnote>

<footnote n="5">⁵Laboratoire d'Océanographie et du Climat: Expérimentations et Approches Numériques, UMR 7159, CNRS, IPSL, and UPMC, 75252 Paris, France</footnote>

<footnote>*now at: Northwest Atlantic Fisheries Centre, Fisheries and Ocean Canada, St John's, Canada</footnote>
to represent unresolved downgradient eddy fluxes, but through more scale-selective (often bilaplacian) operators designed to preserve the resolved part of the mesoscale spectrum. Non-linear energy transfers occur at the mesoscale in such models, and may then feed back onto larger space and time scales through inverse cascade or rectification processes (e.g. Zhai et al., 2004; Penduff et al., 2007).

Primitive equation models are being implemented and assessed at even higher resolution; recent studies show that they yield further dynamical improvements (e.g. Smith et al., 2000; McClean et al., 2002; Masumoto et al., 2004; Drillet et al., 2005; Treguier et al., 2005; Kelly et al., 2007; Chanut et al., 2008; Hecht and Smith, 2008; Hecht and Hasumi, 2008). Most of these simulations are presently restricted to either individual basins and/or decade-long integrations. Several years of research and computational power increase will probably be needed before these promising models can eventually be used by a wide community in long-term, forced and coupled global simulations.

The superiority of $\frac{1}{2} - \frac{1}{6}$° over laminar ocean models for large-scale ocean simulations has been demonstrated by many authors since the pioneering works by Holland et al. (1983) and Semtner and Chervin (1988a,b) and largely confirmed since then, mostly in terms of mean states, mesoscale features, and their mutual interactions (e.g. Böning and Budych, 1992; Beckmann et al., 1994; Böning and Bryan, 1996; Dengg et al., 1996; Haidvogel et al., 2000; Gulev et al., 2007; DYNAMO Group, 1997). Eddy-admitting models are also expected to be beneficial when coupled to the atmosphere (see e.g. Fanning and Weaver, 1997; Roberts et al., 2004), and are presently being substituted for laminar models in hindcasts or forecasts of the full climate system. Global eddy-admitting models presently appear as an interesting trade-off between available computer resources, the partial resolution of mesoscale effects, and the need to perform several multi-decadal integrations to study climate-related oceanic changes.

The main goal of this study is to quantitatively compare the realism of four state-of-the-art global ocean/sea-ice simulations representative of laminar and eddy-admitting classes (2° and 1°, and $\frac{1}{2}$° and $\frac{1}{4}$° resolutions, respectively) over a large range of timescales. We particularly focus on interannual timescales at which the impact of ocean model resolution has not been quantitatively assessed. These four ocean/sea-ice 1958–2004 hindcasts (DRAKKAR Group, 2007) were performed by the DRAKKAR consortium<fnref n="1" />. The AVISO altimeter SLA dataset is chosen as our main reference for its unique quasi-global character and its wide range of space-time scales, despite its restriction to the ocean surface. The present assessment is performed on atmospherically-forced models, but our complementary investigation of simulated interannual variabilities at scales larger than 6° may also help illustrate how in-

creased resolution could modify the large-scale behavior of numerical oceans in coupled mode.

Model-observation and model-model SLA comparisons will be performed globally and locally within various frequency ranges addressing the following three considerations: magnitude and spatial distribution of the simulated SLA variability, and temporal correlations between local timeseries of observed and simulated SLAs. This study also complements and extends the assessment of Drakkar simulations (e.g. Barnier et al., 2006; DRAKKAR Group, 2007; Treguier et al., 2007; Penduff et al., 2007; Lique et al., 2009; Lombard et al., 2009) following a dedicated approach (model-observation collocation, dedicated metrics) to yield a quantitative benchmark between various classes of models.

The following section presents the four model setups, the methods used to collocate their outputs and the observations, to filter the results, and the metrics used to compare them together. After a brief comparison in section 3 of the four simulated mean sea-surface height (MSSH) fields with respect to Niiler et al. (2003)'s observational estimate, the effect of resolution is assessed in terms of sea-level anomalies using the three criteria introduced above within three ranges of timescales separated by 5 and 18 months cutoff periods (Sects. 4, 5, 6). In Sect. 7, the comparison of interannual variabilities will be restricted to scales larger than 6°, i.e. the range of scales that is resolved in the four simulations, and in present climate-oriented coupled simulations.

## 2 Model configurations and assessment procedure

### 2.1 Model configurations

Our model setups are based on the NEMO code (Madec, 2008) and differ by their horizontal resolutions ($\frac{1}{4}$°, $\frac{1}{2}$°, 1°, 2°). Figure 1 shows that in the 2° and 1° configurations, the meridional resolution *dy* is increased toward the equator where it reaches $\frac{1}{2}$° and $\frac{1}{3}$°, respectively. The four simulations share the same vertical discretization (46 geopotential levels whose spacing progressively increase from 6m at the surface to 250m at the bottom), the same parameterization of unresolved vertical mixing and convection processes (TKE turbulent closure scheme), the same linearized free-surface formulation (Roullet and Madec, 2000), and the same surface forcing. The four runs were driven over this 47-year period by the same hybrid forcing function (referred to as DFS3), which is fully described and compared to the CORE forcing (Large and Yeager, 2004; Griffies et al., 2009) in Brodeau et al. (2006, 2010): precipitation and radiative fluxes come from satellite products; air-sea and air-ice turbulent fluxes are computed through bulk formulae from surface model variables and corrected 10-m atmospheric state variables from ECMWF (ERA40 reanalysed fields before 2002, ECMWF analysed fields afterwards). Uncertainties in precipitation fields, along with the need to limit drifts requires in all runs a

<footnote n="1">¹http://www.ifremer.fr/lpo/drakkar/</footnote>
<figure><img src="image_39.png" /><figcaption>**Fig. 1.** Zonally-averaged zonal (dashed) and meridional (plain) resolution of the four model grids (km, log scale). The meridional resolution is enhanced near the equator in the 2° and 1° models. Dots indicate the zonally-averaged first Rossby radius (Chelton et al., 1998).</figcaption></figure>

moderate (60-day timescale over the upper 10m, i.e. 600-day over 100 m mixed layers) relaxation of sea-surface salinity toward the monthly Levitus et al. (1998) climatology. The technical report by Molines et al. (2006) describes the parameterizations and numerical choices made in the 1/4° simulation<fnref n="2" />. Timmermann et al. (2005) and Cravatte et al. (2007) provide complementary information on the 2° model (despite a different vertical resolution and forcing); the same 2° model is compared to other ocean/sea-ice components of climate models by Griffies et al. (2009).

In all runs, the bottom topography is discretized as partial steps for an accurate representation of topographic slopes and 1/H contours. The momentum advection scheme (Arakawa and Lamb, 1981) conserves both energy and potential enstrophy. These latter two choices were shown to yield a remarkably realistic 1/4° global model solution in preliminary climatological simulations, thanks to improved numerical schemes and subsequent eddy-topography interactions (Penduff et al., 2007; Le Sommer et al., 2009). We do not expect such a benefit at coarser resolutions where weak or parameterized mesoscale turbulence cannot drive realistic topographically-rectified mean flows. The same quadratic bottom friction parameterization is used in all simulations (see Penduff et al., 2007).

Although as many numerical and physical parameters as possible (e.g. initial states, surface forcing, bulk formulae, vertical physics, bottom friction, etc) were kept identical in the four configurations, certain parameters needed to be adjusted according to the numerical and physical requirements of each configuration, in order to yield the most consistent and realistic solution in their specific dynamical regimes. In both eddy-admitting simulations, temperature and salinities are mixed along isopycnals through a laplacian operator; the associated diffusion coefficients at the equator equal 300 and 600 m²s⁻¹ in the 1/4° and 1/2° models respectively, and de-

crease proportionally to the grid size. Horizontal viscosity is achieved by bilaplacian operators at 1/4° and 1/2° resolutions; associated coefficients at the equator equal 1.5 × 10¹¹ and 12 × 10¹¹ m⁴s⁻¹ respectively, and vary as the cube of the gridsize. A free-slip sidewall boundary condition is used in the 1/4° and 1/2° models.

In the 2° and 1° simulations, the isopycnal laplacian tracer mixing is complemented by a GM90 parameterization, and momentum mixing is performed by a horizontal laplacian operator. Details about the spatial distribution of the associated coefficients may be found in Cravatte et al. (2007) for the 2° model; these coefficients were simply divided by two in the 1° simulation. The sidewall boundary condition is free-slip at 1° resolution, but no-slip at 2° resolution, as is usually done in climate-oriented simulations with this configuration.

All simulations are started from rest in 1958 on the first of January, and are archived on their native grids over the 47-year (1958–2004) model runs as successive 5-day averages labeled by the central date (e.g. the 5-day average between January first and fifth is dated January thrid at noon). The present study focuses on the last 12 years (1993–2004) when altimeter observations are available; this ensures a relatively long (35-year) and identical period of spinup in the four integrations.

## 2.2 Model-observation and model-model comparison methodology

### 2.2.1 Collocation and filtering

The series of simulated SSH 5-day averages are first linearly interpolated in time and space at the resolution of AVISO SLA fields (weekly and on a 1/3° × 1/3° Mercator grid)<fnref n="3" /> on the same space-time domain (1993–2004 period between both polar circles). Each pair of collocated maps (simulated

<footnote n="2">²For clear reference, the 1/4°, 1/2°, 1°, and 2° simulations investigated in the present study are referred to as ORCA025-G70, ORCA05-G70.113, ORCA1-R70, and ORCA246-G70 in the DRAKKAR simulation ensemble, respectively</footnote>

<footnote n="3">³To ensure consistency over the period 1993–2004, we used delayed-time SLA maps generated from two simultaneously operating altimeters (Topex/Poseidon or Jason-1+ERS1/2 or Envisat). These maps resolve spatial scales as short as 100 km at midlatitudes, according to AVISO.</footnote>
<figure><img src="image_40.png" /><figcaption>**Fig. 2.** Upper left panel: mean sea surface height (MSSH, contours and colors) from Niiler et al. (2003). Lower 4 panels: MSSH from each model simulation (contours) and spatially-smoothed MSSH differences between simulations and Niiler et al.'s reference (colors); regions where simulated MSSHs are higher (lower) than the reference are shown in red (blue). Units are cm; the five maps share the same contours (from -260 to 160 cm by 20 cm intervals); the lower four panels share the same colorscale. Upper right panel: spatial correlation between each simulated MSSH map and Niiler et al.'s reference in 28 latitude bands. See text for details.</figcaption></figure>

SSH, observed SLA) is then masked where either real or simulated sea-ice is present at that date. The 1993–1999 temporal averages of both timeseries are then subtracted at each grid point, as done routinely in AVISO; the global average of these fields is finally subtracted every 5 days from the observed and simulated fields, as recommended by Great-batch (1994). This processing provides the time-dependent collocated 5-member (observations plus 4 simulations) SLA dataset.

The time averages subtracted during this process are simply the mean SSH (MSSH) fields simulated by each model, available on the AVISO Mercator grid, and computed over the period 1993–2004. These simulated MSSHs are compared to Niiler et al. (2003)'s observation-based, 1/2-resolution, 1992–2002 MSSH, after collocation of the latter field onto the AVISO grid. The small temporal shift between

observed and simulated MSSH time-averaging periods is not expected to adversely affect our comparisons given the large temporal overlap. We finally subtracted from observed and simulated MSSH maps a globally-uniform scalar taken from each individual map at (180° W, 0°). Contours in Fig. 2 show the global distribution of MSSHs derived from observations (top left) and simulations (4 panels below), and will be commented upon in Sect. 3.

To evaluate the skill of our models (and the impact of resolution) at different time scales, the 5-member time-dependent SLA dataset presented above is split into 3 frequency bands with two passes of a low-pass Lanczos filter. The first pass uses an 18-month cutoff period and yields the collocated interannual fluctuations of simulated and observed SLAs. The second pass uses a 5-month cutoff period to split the remaining signals into "quasi-annual" and
<figure><img src="image_41.png" /><figcaption>Fig. 3. Scheme of the 4 bands in which model-observation comparisons are performed. The Large-scale Interannual band is shown in light gray.</figcaption></figure>

”mesoscale” components. In summary, the timescales *T* of the interannual, quasi-annual, and mesoscale bands correspond to *T* > 18 months, 5 months < *T* < 18 months, and 14 days < *T* < 5 months, respectively; these bands are close to those chosen by Berloff and McWilliams (1999) to analyze the nonlinear response of idealized models at increasing resolution.

This paper puts some emphasis on simulated interannual
variabilities; the statistics described below are computed in
this frequency band from collocated SLA timeseries both
(i) directly, and (ii) after an additional low-pass Lanczos
isotropic 2-D filtering in space of each collocated SLA map.
Comparisons between observed and simulated SLA transects
showed that our coarsest-resolution model (2°) isotropically
resolves spatial scales larger than about 6°×cos(latitude).
This cutoff wavelength was chosen to further compare ob-
served and simulated interannual SLA characteristics over
the same range of resolved spatial scales, i.e. 6°-to-global;
this spatially and temporally low-passed filtered dataset will
be denoted as "large-scale interannual", and will be analyzed
in Sect. 7. The Lanczos filtering technique is described in
Duchon (1979); it was chosen for its relative simplicity and
its ability to provide "clean" signals devoid of Gibbs oscilla-
tions, in one and two dimensions. Figure 3 summarizes the
space-time scales considered in the following.

2.2.2 Model-observation comparison statistics

Collocated SSH timeseries at location (i, j) are noted η^A(t) for AVISO and η^m(t) for the m^th model ⁴. Let φ̃^t denote the time average of any variable φ. Simulated and observed MSSH maps, noted η̃^m^t(i, j) and η̃^A^t(i, j) respectively, are shown as contours in Fig. 2. Differences between simulated and observed MSSHs are shown in colors in the same figure, after spatial filtering.

- SLA temporal standard deviations at location (i, j) for model *m* are computed in each frequency band as:

$$
\sigma^m(i, j) = \sqrt{[\eta^m(t) - \overline{\eta^m}^t]^2}. \quad (1)
$$

The temporal standard deviations of AVISO SLAs are
denoted as σ^A(i, j); the ratio σ^m/σ^A will be called re-
solved variability. Fig. 4 shows σ^A(i, j) maps (first row)
and σ^m(i, j) maps for each simulation (other rows) in
all frequency bands (columns).

- Temporal correlations between simulated and observed SLA timeseries are computed at each (i, j) and for each model *m* in each frequency band as:

$$
C_t^m(i, j) = \frac{\overline{(\eta^m(t) - \overline{\eta^m}^t)(\eta^A(t) - \overline{\eta^A}^t)}{\sigma^m(i, j)\sigma^A(i, j)}. \quad (2)
$$

Associated significance level maps are then estimated
as proposed by von Storch and Zwiers (1999, their
equation 12.55) for each simulation, each frequency
band, and each wet point (i, j) independently. This ap-
proach requires the computation of local autocorrela-
tion functions at lag τ of all η<sup>m</sup>(t) and η<sup>A</sup>(t) timeseries;
these autocorrelation functions are noted ρ<sup>m</sup>(i, j, τ) and
ρ<sup>A</sup>(i, j, τ), respectively. Individual temporal correla-
tions C<sub>t</sub><sup>m</sup>(i, j) are then considered significant at the 95%
confidence level where

$$
C_t^m(i, j) > 2 \sqrt{\frac{1}{nt} \left[ 1 + 2 \sum_{\tau=1}^{nt-1} \rho^m(i, j, \tau) \rho^A(i, j, \tau) \right]} \quad (3)
$$

Note that this method does not assume that individual
timeseries are white or have the same autocorrelations,
and provides significant temporal correlation maps for
each simulation and each frequency band. Only signifi-
cant temporal correlations, i.e. satisfying Eq. 3, are kept
in the following. Resolution-induced changes in signif-
icant temporal correlations are shown in Fig. 5 for each
frequency band (columns).

Mean and variable surface topography features are com-
pared to their observed counterparts within 28 periodic lat-
itude bands labeled $\lambda \in (1; 28)$ spanning the (70° S–70° N)
range by 5° intervals.

- Temporal standard deviation $\sigma^m(i, j)$ maps are split into each latitude band and rearranged as a long one-dimensional vector of length $l$ containing elements noted $\sigma^m(\lambda, l)$. This vector is then averaged along dimension $l$ to yield $\sigma^m(\lambda) = \overline{\sigma^m(\lambda, l)}^l$, i.e. the average of $\sigma^m$ in each latitude band. These quantities and their observed counterparts (noted $\sigma^A(\lambda)$) are shown in various colors in the first row of Fig. 6 for each frequency band (columns).

<footnote>$^4$i ∈ [1; 1080], $j \in [1; 915]$, $m \in [1; 4]$, $t \in [1; nt]$, $nt = 625$.</footnote>
<figure><img src="image_42.png" /><figcaption>Fig. 4. SLA standard deviation maps (cm) in AVISO observations ($\sigma^A$, top) and the 4 model simulations ($\sigma^m$, other rows). In this and the following figures, results are shown from left to right in the mesoscale, quasi-annual, interannual, and large-scale interannual bands. See text for details.</figcaption></figure>

- The same procedure is applied to (significant) temporal correlation maps $C_t^m(i, j)$. This provides $C_t^m(\lambda)$, i.e. the average of temporal correlations between observed and simulated local timeseries in each latitude band, for each simulation and frequency band (second row in Fig. 6). These middle panels also show in dashed lines the density $p_t^m(\lambda)$ of significant temporal correlations within each latitude band, i.e. the number of grid points with significant temporal correlation divided by the number of wet points in each band of latitude ($p = 1$ when Eq. 3 is verified at all wet points).

- We also quantify the spatial correlation between observed and simulated (stationary) SLA standard deviation maps ($\sigma^A(i, j)$ and $\sigma^m(i, j)$) in each latitude band. Both maps are split and rearranged as $l$-long vectors within each band; the correlation between these two vectors, i.e. the spatial correlation between $\sigma^A$ and $\sigma^m$

maps in this latitude band, is then computed as follows:

$$C_s^m(\lambda) = \frac{\overline{(\sigma^m(l) - \overline{\sigma^m l})(\sigma^A(l) - \overline{\sigma^A l})^l}}{\alpha^m(\lambda)\alpha^A(\lambda)} \quad (4)$$

Here, $\overline{\phi}^l$ denotes the spatial average of any variable $\phi$ over a latitude band. The spatial standard deviation of $\sigma^m$ within each latitude band is computed as $\alpha^m(\lambda) = \sqrt{[\sigma^m(l) - \overline{\sigma^m l}]^2}$; its AVISO counterpart is noted $\alpha^A(\lambda)$. These quantities are shown for each run (color) and in each frequency band (columns) in the third row of Fig. 6. Spatial correlation coefficients between two satellite-derived $\sigma^A$ fields in distinct frequency ranges will be noted $C_s^A$. Note that global spatial correlation coefficients between global $\sigma^A(i, j)$ maps are also mentioned in the following; they simply
<figure><img src="image_43.png" /><figcaption>**Fig. 5.** Changes in SLA temporal correlations ($C_t^{m2} - C_t^{m1}$) when the model resolution is enhanced from *m*1 to *m*2. Red (blue) denotes a resolution-induced increase (decrease) in $C_t$; regions where $C_t^{m2}$ and/or $C_t^{m1}$ are insignificant are left in white. Results are shown for 4 resolution transitions: from 2° to 1/4° (row 1), 1/2° to 1/4° (row 2), 1° to 1/2° (row 3), and 2° to 1° (row 4).</figcaption></figure>

correspond to $C_s^m$ or $C_s^A$ computed over a wide latitudinal band extending between 66° S and 66° N.

- A similar spatial correlation computation is performed on MSSH maps ($\bar{\eta}^{A'}$ and $\bar{\eta}^{m'}$) to quantify the agreement between simulated MSSHs and Niiler et al. (2003)'s observational estimate as a function of latitude. This quantity is shown for each run in the upper right panel of Fig. 2.

Fig. 7 complements Fig. 6: it shows in the same latitude and frequency bands the successive impacts of resolution increases (from *m*1 to *m*2) on resolved SLA variables ($\bar{\sigma}^{m2}-\bar{\sigma}^{m1}$, first row), on significant temporal correlations ($C_t^{m2} - C_t^{m1}$, second row), and on spatial correlations ($C_s^{m2} - C_s^{m1}$, third row). Results are shown for resolution transitions 2°–1° (blue), 1°–1/2° (green), 1°–1/4° (red). The impact of the 2°–1/4° transition, shown in cyan, is the sum of the three incremental changes.

Finally, the middle row in Fig. 6 is complemented by Fig. 8, which shows in all frequency bands (columns) and in all simulations (colors) how significant temporal correlations vary as the distance from the coast increases (abcissae, in °).

Simulated and observed MSSHs are compared in the next section. In Sects. 4, 5, 6, the observed and simulated SLA variabilities are compared (without spatial filtering) within the frequency bands defined above in terms of intensity, spatial distribution ($C_s$), and temporal phase ($C_t$). The impact of resolution on the large-scale interannual variability is presented in Sect. 7.

We remind the reader that all statistics are based on observed and simulated SLA fields that were collocated at the same spatio-temporal resolution (1/3° × 1/3° × 7 days) and filtered identically, and that the 1° and 2° model grids are refined near the equator, where all four models share comparable meridional resolutions, as shown in Fig. 1.
<figure><img src="image_44.png" /><figcaption>**Fig. 6.** Observed and simulated SLA standard deviations ($\sigma(\lambda)$) in cm, upper panels), significant temporal correlations and density of gridpoints with significant temporal correlations ($C_t^m(\lambda)$ in plain lines and $p_t^m(\lambda)$ in dashed lines, middle panels), and spatial correlations between observed and simulated maps of SLA ($C_s^m(\lambda)$, lower panels). Results are shown in average over latitude bands. See text for details.</figcaption></figure>

### 3 Impact of resolution on mean sea-level maps

The contours in Fig. 2 present mean SSH fields from Niiler et al. (2003)'s reference and the 4 simulations. Spatially-smoothed MSSH differences between simulations and this reference are shown in color. Barnier et al. (2006) analyzed a previous 1/4° DRAKKAR simulation (referred to as G22, performed with the same numerics but a different forcing), and showed that many regional MSSH patterns were as realistic as in higher-resolution models. Penduff et al. (2007) and Le Sommer et al. (2009) attributed this skill to a more accurate representation of topographic rectification processes 5. Our main focus is on variabilities: only a general description of MSSH differences is given here. Note that Niiler et al.

(2003)'s field was derived by inversion and may not represent the "true" MSSH.

Model-observation spatial correlations, shown as a function of latitude in the upper-right panel, remain above 0.98 between 30° S and 30° N: the most realistic MSSH features are found in this tropical band in the four simulations. The maps show that all simulations tend to slightly overestimate the large-scale MSSH in this tropical band (+5–10 cm with respect to 180° W, 0°), especially in the Indonesian Archipelago, the North Atlantic and Indian Oceans. The main mismatches between simulations and our reference are found in both hemispheres poleward of 60° and between 30 and 40°, where MSSH spatial correlations drop slightly below 0.9.

As in most eddy-admitting models, the Gulf Stream (GS) overshoots at Cape Hatteras at 1/4°. West of about 50° W, its path is in good agreement with our reference (as in Barnier et al., 2006), although the stream tends to recirculate to the southeast; between 50 and 30° W, the northward-flowing North Atlantic Current (NAC) is displaced to the east in this simulation. At coarser resolutions, the GS-NAC system gets much broader than observed. These discrepancies are consistent with the dominance of blue (low MSSH) visible along the GS-NAC system in Fig. 2.

<footnote n="5">⁵In this previous and the present 1/4° simulation, most mean currents have the same paths except the North Atlantic Current (NAC), whose location was more realistic in the earlier study. In this earlier simulation, however, the buoyancy losses in the Nordic Seas, the density of overflows, and the Atlantic Meridional Overturning Circulation were all overestimated. The present 1/4° simulation has more realistic air-sea fluxes, water masses and overturning, but a displaced NAC. A connection between these features is therefore likely, but not fully understood yet.</footnote>
<figure><img src="image_45.png" /><figcaption>**Fig. 7.** Resolution-induced changes in resolved SLA variabilities ($\frac{\sigma_{m2}}{\sigma_A} - \frac{\sigma_{m1}}{\sigma_A}$, first row), SLA temporal correlations ($C_l^{m2} - C_l^{m1}$, second row), and SLA spatial correlations ($C_s^{m2} - C_s^{m1}$, third row) as a function of latitude (abcissae). Results are shown for four resolution transitions $m1 \rightarrow m2$ (colors, see legend). Lines in the second row are dashed at latitudes where $p_l^{m1}$ or $p_l^{m2}$ (densities of points with significant temporal correlations) is smaller than 20%.</figcaption></figure>

<figure><img src="image_46.png" /><figcaption>**Fig. 8.** Significant model-observation temporal correlations ($C_l^m$, ordinates) as a function of the distance from the coast (abcissae, °), in the mesoscale (left), quasi-annual (middle) and interannual (right) frequency bands. All significant $C_l^m$ values found at the same distance from the land mask have been averaged.</figcaption></figure>

In the western North Atlantic subpolar gyre, increases in resolution yield decreases in simulated MSSHs; indeed, the time-averaged cyclonic transport around the Labrador and Irminger Seas monotonically increases with increasing resolution from about 20 Sv at 2° to more realistic values (45 Sv) at $\frac{1}{4}$°. This strong impact of resolution has already been mentioned by Treguier et al. (2005).

Positive MSSH biases are found within and north of simulated Kuroshios; this mismatch is associated with an overshoot of this current in eddy-admitting simulations, and to the underestimation of the western North Pacific subpolar gyre intensity in laminar models.
The minima in MSSH spatial correlations, visible at 40° S north of the Kerguelen Plateau in the four simulations, come from a local 4° northward shift of the northernmost branch of simulated Antarctic Circumpolar Currents (ACC). These shifts, denoted by negative MSSH biases in Fig. 2's lower 4 panels, appear northwest of Kerguelen (60° E) in all runs except at $1/4$° resolution where it is found northeast of Kerguelen (80–100° E). This latter bias was found by Barnier et al. (2006) in the previous $1/4$° DRAKKAR simulation mentioned above. Its appearance in the presence of (partially) resolved eddies supports the hypothesis that in this particular region, our $1/4$° model overestimates topographic constraints associated with eddy-topography interactions (also see Penduff et al., 2007).

MSSH spatial correlations also decrease in the Southern Ocean where mean circulation biases appear at regional scale. All simulations, especially at laminar resolution, tend to overestimate the reference MSSH in the Atlantic sector and around 150° W south of the reference ACC: simulated ACCs tend to spread southwards compared to the reference.

Colors in Fig. 2 do not reveal many resolution-induced changes in large-scale MSSHs. The clearest benefits of higher resolution on simulated MSSH appear most clearly on the contours of Fig. 2: in many regions, the scale and location of mean surface fronts become increasingly realistic with increasing resolution.

# 4 Impact of resolution on the high frequency (mesoscale) SLA variability

In the real ocean, mesoscale SLA standard deviations reach high values (above 13 cm, see the upper left panel in Fig. 4) in the main eddy-active regions: GS, NAC, Kuroshio (KS), Brazil-Malvinas Confluence, East Australian Current, Mozambique Channel, and south of Australia and America. The mesoscale variability approaches 5 cm in zonal average between 35–40° N, 40° S and 55–60° S (upper left panel in Fig. 6). Secondary maxima, exceeding 7 cm locally, are found all along the ACC, within the Indo-Pacific subtropics (20–30° S and 20–30° N), and along the equatorial path of Tropical Instability Waves.

## 4.1 Variability levels

As expected, mesoscale variability is very weak in both laminar models, i.e. from around 20% of observed levels in the tropics up to 45% near 60° S (Figs. 4 and 6); their meridional distributions are also very similar. Increasing resolution (and decreasing dissipation) yields a monotonic increase in mesoscale variability: 29, 35, 42 and 56% of its globally-averaged observed magnitude are simulated at $2^\circ, 1^\circ, \frac{1}{2}^\circ$ and $\frac{1}{4}^\circ$ models, respectively.

The $1/4$° mesoscale variability roughly accounts for about 50% of its observed levels within 40° S–40° N; this relatively

small ratio is explained in part by the still moderate resolution, but might also be due to the absence of small spatial scales in the atmospheric forcing (Milliff et al., 1996); this leaves room for improvement through further resolution increases and more realistic forcing. Both the $1/4$° resolved mesoscale variability, and its increase with resolution, are stronger poleward of about 30–40° (ACC, Agulhas Retroflection, Brazil-Malvinas Confluence, subpolar North Atlantic). At these latitudes, none of our model grids can adequately resolve the first Rossby radius (see Fig. 1) and then the most unstable baroclinic waves (baroclinic instability occurs at larger scales). However, the contribution of barotropic, topographically-influenced fluctuations on the SLA variability increases with latitude at these frequencies (Guinehut et al., 2006; Vinogradova et al., 2007). This resolution-induced increase in high-latitude mesoscale resolved variability might, thus, come from stronger (and/or less damped) barotropic motions, rather than enhanced baroclinic instability.

The upper left panel in Fig. 7 confirms the major contribution of the $1/2 - 1/4$° (red) resolution increase on the intensity of resolved mesoscale SLA variability, in particular in the Southern Ocean.

## 4.2 Distribution of variability

Spatial correlations $C_s^m$ between observed and $1/4$° mesoscale variability maps lie between 0.6 and 0.8 nearly everywhere (lower left panel in Fig. 6). Spatial correlations markedly decrease at $1/2$° resolution, in particular south of 20° S and at the KS-GS-NAC latitudes, where the strong mesoscale activity gets damped and/or distorted. Coarsening the resolution to 1° further deteriorates the high-frequency variability distribution in the Southern Ocean (which almost disappears along the ACC, see Fig. 4), and in the 0–15° N band. Despite their two-fold resolution ratio, both laminar models produce similar high-frequency variability maps which are unrealistic in the regions cited above (weak spatial correlation values). It is likely that the non-eddying character of both grids and/or the use of the GM90 parameterization in both simulations are responsible for this similarity.

The lower left panel in Fig. 7 summarizes the contribution of successive resolution increases on the overall improvement from $2^\circ$ to $1/4^\circ$; the $1^\circ-1/2^\circ$ transition (green) improves the distribution of high-frequency SLA variability as much as the $1/2 - 1/4^\circ$ (red) transition in the Southern Ocean and the northern tropical oceans.

## 4.3 Phase of temporal variability

Left panels in Fig. 5 show that resolution affects mesoscale temporal correlations in diverse ways; for instance, resolution increases induce a quasi-systematic decrease in the eastern tropical Pacific, a systematic increase on the Argentinian Plateau, and rather complex changes over most of the deep
oceans. Significant temporal correlations $C_t^m$ are averaged in latitude bands in the middle left panel in Fig. 6. The red line in the middle left panel on Fig. 7 shows that the $\frac{1}{2}^\circ-\frac{1}{4}^\circ$ transition slightly decorrelates (-0.015) local SLA timeseries at most latitudes in this frequency range; this very likely comes from the stronger mesoscale activity simulated at $\frac{1}{4}^\circ$.

A clear impact of increased resolution in the laminar regime is found in the Northern Hemisphere, where temporal correlations increase by about 0.05 in zonal average in the 2°-1° transition (also see the middle left panel in Fig. 7). The lower left panel in Fig. 5 indicates that these increases mostly occur along the coasts: in the Indonesian Archipelago and along the edges of the subpolar gyres. Provided that AVISO data are accurate enough near the coast for our metrics to be sensible, the left panel in Fig. 8 confirms the latter tendency at global scale: significant temporal correlations at mesoscale frequencies increase toward the continents, but substantially less in the 2° model than at higher resolution. This suggests that the 2° resolution, which is used in many climate prediction systems, is too coarse to properly resolve forced boundary waves at these frequencies (e.g. Kelvin waves), i.e. part of the ocean adjustment to the atmosphere. The ocean dynamics, the influences of topography and forcing simulated at 1°, $\frac{1}{2}^\circ$, and $\frac{1}{4}^\circ$ resolutions seem to have comparable skill in resolving these relatively fast, forced near-coastal processes. Whether model resolutions finer than $\frac{1}{4}^\circ$ would further increase (through a more accurate representation of forced signals) or decrease (by admitting more chaotic mesoscale turbulence) near-coastal $C_t^m$ values remains to be determined.

# 5 Impact of resolution on the medium frequency (quasi-annual) SLA variability

The quasi-annual AVISO SLA standard deviations globally exceed their mesoscale counterparts by about 10%, up to 80% around 35° N (second column, first row in Fig. 6). The spatial correlation between global variability maps at quasi-annual and mesoscale bands is large both in the AVISO dataset (0.76), and in the $\frac{1}{4}$ model simulation (0.72). Accordingly, the RMS difference between global, observed $\sigma^A$ maps of quasi-annual and mesoscale variability is relatively small (2.3 cm), and is the same in the $\frac{1}{4}$ dataset. In other words, and as can be seen in Fig. 4, most quasi-annual and mesoscale variability maxima are found in the same areas (KS, GS-NAC, Agulhas Retroflection, ACC, Brazil-Malvinas Confluence), both in the AVISO and $\frac{1}{4}$ datasets.

## 5.1 Variability levels

Zonally-averaged simulated variability levels nicely follow the observations north of 35° S lying around 60–90% of their observed value, with only a small impact of resolution (first row, second column in Fig. 7). Variability maxima, peak-

ing at quasi-annual timescales, are observed along 10° N ($\sigma^A$~6cm) and 10° S ($\sigma^A$~4cm) in the Indian and Eastern North Pacific basins (Fig. 4). These zonally-extended maxima are present in the four simulations; their magnitude reach about 70 to 85% of observed levels, except in the 2° model where it drops by another 10%.

South of about 45° S, this quasi-annual resolved variability remains above 80% in the $\frac{1}{4}^\circ$ simulation, but falls to around 50% at $\frac{1}{2}^\circ$ resolution, and below 40% in both laminar runs. In this quasi-annual band, the $1^\circ-\frac{1}{2}^\circ$ and $\frac{1}{2}^\circ-\frac{1}{4}^\circ$ resolution transitions contribute to comparable enhancements of the resolved SLA variability, from about +5% in the KS and GS/NAC systems up to +20% in the Southern Ocean (first row, second column in Fig. 7).

## 5.2 Distribution of variability

Figure 6 (second column, third row) shows that the agreement between the observed and $\frac{1}{4}^\circ$ maps of quasi-annual variability is greatest between 20° S and 20° N ($C_s^{1/4} \sim 0.9$) and remains substantial (0.6–0.7) at most latitudes. As in the mesoscale range: (i) quasi-annual spatial correlations are weakest in both laminar simulations, (ii) $\frac{1}{2}^\circ$ results are intermediate between $C_s^{1/4}$ and $C_s^{2\circ}$, and (iii) the largest improvements in quasi-annual SLA variability maps are mostly obtained in the $1^\circ-\frac{1}{2}^\circ$ and $\frac{1}{2}^\circ-\frac{1}{4}^\circ$ resolution transitions poleward of 30° (second column, third row in Fig. 7).

## 5.3 Phase of temporal variability

The impact of resolution on quasi-annual temporal correlations is shown in the second columns of Figs. 5 and 6 (second row). Quasi-annual temporal correlations exhibit a marked increase in the eastern tropical Pacific with resolution, indicating an improvement in the phase of quasi-annual equatorial Rossby waves every time resolution is doubled. A marked decrease of quasi-annual temporal correlations appears in the Southern Ocean from laminar to eddy-admitting resolutions (especially in the $1^\circ-\frac{1}{2}^\circ$ transition south of 55° S, see the second row, second column in Fig. 7). This temporal decorrelation, which concerns only parts of the Southern Ocean, resembles its counterpart found at longer timescales (see next sections) and may share the same origin. This issue is commented upon in the conclusion.

Significant temporal correlations in the quasi-annual band increase toward the continents in all simulations (middle panel in Fig. 8). Near-coastal $C_t^m$ values increase with resolution between 2° and $\frac{1}{2}^\circ$, but drop by about 0.06 when resolution reaches $\frac{1}{4}^\circ$, i.e. back to their 2° values. However, these changes occur over a cross-shore scale (about $\frac{1}{3}^\circ$) that seems small compared to those (about 1°) of coastally-trapped processes that may prevail at quasi-annual timescales (e.g. upwellings). The interpretation of these differences is thus unclear and would require dedicated investigations.
# 6 Impact of resolution on the low frequency (interannual) SLA variability

The global spatial correlation coefficient (0.72) between observed variability maps $\sigma^A$ in the interannual and quasi-annual bands (upper two center panels in Fig. 4) is about as large as its counterpart between quasi-annual and mesoscale maps (0.76, see previous section). This correspondence between observed $\sigma$ distributions in the 3 frequency bands is also clear in the $\frac{1}{4}^\circ$ simulation (second row in Fig. 4): in this simulation, interannual and quasi-annual variability maps have a spatial correlation equal to 0.70; quasi-annual and mesoscale variability maps have a spatial correlation equal to 0.72. In other words, both interannual and quasi-annual SLA variabilities happen to be strong in the main eddy-active regions (KS, GS-NAC, Agulhas Retroflection, ACC, Brazil-Malvinas Confluence), both in the real ocean and in the $\frac{1}{4}^\circ$ simulation. This geographical correspondence, which is somewhat less marked at low latitudes, suggests that in several mid- and high-latitude eddy-active areas, the interannual variability and the mesoscale variability are connected. This hypothesis is further discussed in the conclusion.

## 6.1 Variability levels

Interannual timescales are those at which the $\frac{1}{4}^\circ$ resolution yields the largest and most realistic resolved SLA variability: it reaches 81% globally and at most latitudes (first row, third column in Fig. 6). This global fraction reduces to 71% at $\frac{1}{2}^\circ$ and 60% in both laminar models. A similar decrease with resolution is seen at quasi-annual timescales (76, 68, 62, 61%, respectively, from $\frac{1}{4}^\circ$ to $2^\circ$). Resolved interannual variabilities undergo comparable increases in the $1^\circ-\frac{1}{2}^\circ$ and $\frac{1}{2}^\circ-\frac{1}{4}^\circ$ resolution transitions: about +10% north of 30° N and more than 20% in the Southern Ocean.

## 6.2 Distribution and phase of temporal variability

The distribution and temporal evolution of simulated interannual variabilities are most realistic between 20° S and 20° N (i.e. $C_s^m$ and $C_l^m$ lie around 0.9 and 0.8, respectively). Both metrics, and the magnitude of the SLA variability, are almost identical in the four simulations in this latitude-frequency band. This means that given an adequate forcing, the tropical interannual SLA variability is well simulated when the meridional resolution remains in the 30–60 km range, and remains insensitive to both a 8-fold change in zonal resolution and to the choice of parameterizing baroclinic instability through GM90 instead of resolving it.

The third columns in Figs. 6 and 7 show that poleward of ±30°, the $1^\circ-\frac{1}{2}^\circ$ (and also $\frac{1}{2}^\circ-\frac{1}{4}^\circ$) resolution transitions yield two major effects: an improved spatial distribution of interannual SLA variabilities, and a decrease in significant temporal correlations (third column, second row in Fig. 7) south of 30° S. This temporal decorrelation, which is not pre-

cisely located but distributed at circumpolar scale (third column in Fig. 5), will be commented upon in the conclusion. The right panel in Fig. 8 also shows that over a 1.5°-wide coastal band, interannual temporal correlations are smaller in both eddying simulations compared to their laminar counterparts (0.05–0.1 difference).

# 7 Impact of resolution on the large-scale interannual SLA variability

The comparison of observed and simulated interannual variabilities is now restricted to spatial scales larger than 6°, i.e. on the oceanic patterns that are actually resolved in our 4 models, and at present in most coupled models. We also assess in this way how the explicit resolution of scales smaller than 6° can affect scales larger than this threshold. SLA standard deviations, temporal and spatial correlations are shown in the right columns of Figs. 4 and 6.

In both the observed and $\frac{1}{4}^\circ$ datasets (compare the uppermost two panels in the last two columns in Fig. 4), excluding scales smaller than 6° decreases the interannual variability by up to about 50% locally in the GS, in the Brazil-Malvinas Confluence, in the ACC between 20 and 75° E and in the Pacific sector to a lesser extent. This decrease in interannual variability levels amounts to 14–16% on average over the GS and ACC latitude bands, in both AVISO and $\frac{1}{4}^\circ$ SLA datasets. In other words, a substantial part of the observed and $\frac{1}{4}^\circ$ SLA interannual variability in these eddy-active regions is [i] confined along unstable currents, [ii] collocated with higher-frequency variability, and [iii] accounted for by scales smaller than 6°. At coarser resolutions, the interannual SLA variability is largely accounted for by large-scale motions without significant contribution from scales shorter than 6° (which are not resolved). Therefore, the $\frac{1}{4}^\circ$ resolution appears necessary to admit turbulent eddies (small space and time scales), and to realistically represent the small-scale component of the interannual variability. This may explain why the $\frac{1}{4}^\circ$ resolution yields strong improvements (magnitude and spatial distribution of SLA variability) in all frequency bands simultaneously, in these two eddying regions and especially the Southern Ocean (first row in Fig. 6). Away from eddy-active regions, the last two columns in Fig. 4 are much more similar in both the observations and the four simulations: the interannual variability in quiescent regions is mostly accounted for by motions with scales larger than 6°, and may thus be resolved at coarse resolution.

The third and fourth columns in Fig. 6 show that the impacts of resolution increases that were highlighted earlier in the absence of spatial filtering (i.e. monotonic increase in interannual variabilities, marked improvement of the maps of interannual variability patterns, slight decorrelation in SLA timeseries) are essentially the same when motions with scales smaller than 6° are omitted. Note that choosing a 12° cutoff lengthscale instead of 6° yields the same general results, with
moderate quantitative changes: in this very large-scale inter-
annual band, increased resolution from 2° to 1/4 enhances the
standard deviation of SLA fluctuations (+15–20% poleward
of 40°), decreases temporal correlations by up to 0.2 south
of 40° S, and enhances spatial correlations by about 0.2 pole-
ward of 30°.

To summarize, extending the spectrum of resolved scales
toward the first internal deformation radius progressively im-
proves the representation of interannual variability over a
broad range of spatial scales, including those much larger
than the grid size.

8 Conclusions

The main aim of this study was to assess and better under-
stand the effects of resolution on global ocean model solu-
tions, as considered through sea surface height. The out-
puts from four 47-year global ocean/sea-ice simulations, per-
formed at 2°, 1°, 1/2°, and 1/4° resolutions, and driven by iden-
tical surface forcing, have been collocated onto AVISO sea
level anomaly (SLA) observations. This yielded a consis-
tant 5-member SLA dataset, which was then filtered in time
and space. Model-observation comparisons have been made
within three frequency ranges (interannual, quasi-annual,
mesoscale) considering the following features: magnitude of
SLA variability (using their standard deviations), geograph-
ical distribution of sea-level variability (using spatial corre-
lations between observed and simulated maps of SLA vari-
ability), and phase agreement (significant temporal correla-
tions) between observed and simulated local SLA timeseries.
Simulated mean sea-surface height (MSSH) fields were also
compared to Niiler et al. (2003)'s reference.

Our results provide a quantitative and extensive bench-
mark regarding the comparative skills of two climate-
oriented laminar models (2° and 1°, i.e. resolutions typi-
cal of most ocean components in present climate prediction
models), and two eddy-admitting models (1/2° and 1/4°). Our
main results are summarized below:

- MSSH: Barnier et al. (2006) demonstrated the good agreement between regional MSSH patterns from a comparable 1/4° DRAKKAR simulation and Niiler et al. (2003)'s observation-based estimate. Decreasing resolution quickly deteriorates this regional agreement, but does not strongly modify large-scale MSSH fields except in northern subpolar gyres and north of the Kerguelen Plateau. Spatial correlations between this reference map and its simulated counterparts are largest at 1/4° resolution in most latitude bands, not only at the native 1/3° × 1/3° resolution of the observational product (upper right panel in Fig. 2), but also after removing MSSH features smaller than 6° (not shown).

- Variability levels: increasing resolution from 2° to 1/4° monotonically enhances mesoscale SLA standard devi-

ations in all regions (from 29 to 56% of observed lev-
els in global average). The transition between 1/2° and
1/4° resolution yields the strongest enhancement in eddy-
active regions, as expected. Resolved SLA variabil-
ity levels substantially increase with resolution in the
quasi-annual and interannual ranges as well (i.e. from
61 to 76% and from 59 to 81% of globally-averaged ob-
served levels, respectively). Doubling resolution from
1° to 1/2°, and from 1/2° to 1/4° successively enhances the
resolved SLA variability (up to 20–30% each) at inter-
annual timescales, especially in eddy-active regions.

- Distribution of interannual variability: Most increases in SLA variability levels found in the 1°-1/2° and 1/2°-1/4° resolution transitions are associated over the same latitude bands with large improvements in the spatial distribution of SLA variability, not only at short time scales but at quasi-annual and interannual timescales as well. At 1/4°, many mid- and high-latitude interannual variability maxima are collocated with their quasi-annual and mesoscale counterparts, as in AVISO observations. In other words, the 1/4° model simulates most realistically the broad range of timescales found in confined eddy-active regions (Gulf Stream, North Atlantic Current, Kuroshio, ACC, etc); coarse-resolution models restrict the interannual variability to large spatial scales.

- Resolution-induced improvements in the distribution and magnitude of interannual variability persist (despite minor quantitative changes) when spatial scales smaller than 6° (or 12°) are filtered out of collocated SLA fields. This means that the explicit resolution of small scales substantially improves many aspects of the ocean variability at large scales, i.e. those resolved in present ocean-atmosphere coupled models. This benefit of eddy-admitting resolution had not been quantified previously, and supports the use of eddying ocean components in future climate prediction systems.

- We have showed that the 1/2°-1/4° transition yields large improvements in the intensity and distribution of the interannual variability (especially south of 30° S), but also a weak, yet significant, decorrelation of interan- nual SLA timeseries from their AVISO counterparts. Since the atmospheric forcing is unchanged, these sen- sitivities probably come from the more non-linear dy- namics simulated at 1/4°. To clearly identify the ori- gin of these concomitant changes lies beyond the scope of this study. However, process-oriented numerical studies have shown that mesoscale eddies may drive chaotic low-frequency fluctuations in the ocean (e.g. Jiang et al., 1995; Spall, 1996; Berloff and McWilliams, 1999; Hogg and Blundell, 2006). One may then hy- pothesize that both the strong enhancement and (realis- tic) confinement in eddy-active areas of SLA variabili- ties at all timescales, and the aforementioned temporal
decorrelation, denote the emergence at $1/4$° of a partly chaotic, eddy-driven interannual variability where eddies are strong (hence the localized broad-band temporal variability emerging at $1/4$°). Such an intrinsic interannual variability would likely be out of phase with the actual observations (hence the temporal decorrelation). Testing this hypothesis would require additional simulations and diagnostics, and is left for future studies.

- When resolution gets coarser than 1°, temporal correlations in the mesoscale frequency band decrease along the continents. This suggests that 2° models do not properly resolve the spatial scales of Kelvin waves (whose timescales roughly match this frequency band), hence part of the forced response and adjustment of coastal and shelf areas to the atmosphere or to remote oceanic regions. Resolutions much finer than $1/4$° are needed to properly resolve near-coastal dynamics and forcings (Capet et al., 2004), but enhanced nonlinearities that would emerge on grids much finer than $1/4$° may lead to further decreases in temporal correlations.

- The 2°–1° resolution transition yields slightly enhanced SLA variability levels in the subtropics at interannual timescales, and an enhancement in SLA local temporal correlations at timescales shorter than 5 months. Besides these exceptions, both laminar models yield strikingly similar results in most latitude-frequency bands (Fig. 6), despite a twofold resolution ratio. This similarity may be due to the use in the 2° and 1° models of the GM90 parametrization, which might have a similar damping and distorting influence on quasi-annual and mesoscale variabilities worldwide and on the interannual variability south of 30° S, regardless of the resolution difference. This striking resemblance between both laminar solutions might also be due to their non-eddying character away from the equator, regardless of their particular subgrid-scale parameterizations; in other words, a (partial) resolution of mesoscale eddies and of their dynamical impacts might be necessary to yield substantial improvements in the SLA variability. Determining the individual contribution of resolution and parameterization choices in the laminar regime would require additional simulations (e.g. a $1/2$° simulation with GM90); this is left for future studies.

Eddy-admitting resolution largely improves sea-surface variability patterns (magnitude and distribution) over most of the global ocean, throughout the whole range of timescales considered here (15 day–6 year). Enhanced model resolution also improves the simulated interannual variability at spatial scales (6° to global) that are presently resolved in climate-oriented coupled simulations.

Our conclusions support the pioneering work by Roberts et al. (2004) about the benefits of eddy-admitting resolutions for ocean models, which will certainly remain impor-

tant tools for a long time for both paleoclimatic coupled studies and IPCC-like ensemble climate predictions.

*Acknowledgements.* The authors would like to thank the NEMO Team for numerous technical and scientific interactions, and Jean-Michel Brankart for substantial inputs on statistical issues. Interesting discussions with several colleagues have led to an improvement of this analysis. In particular, the constructive comments from two anonymous reviewers are gratefully acknowledged. This work is a contribution of the DRAKKAR project, which is supported by the Centre National d'Etudes Spatiales (CNES) through the Ocean Surface Topography Science Team (OST/ST), by the Centre National de la Recherche Scientifique (CNRS), the Institut National des Sciences de l'Univers (INSU), the Groupe Mission Mercator Coriolis (GMMC). Computations presented in this study were performed at Institut du Développement et des Ressources en Informatique Scientifique (IDRIS). Partial support from the European Commission under Contract SIP3-CT-2003-502885 (MERSEA project) is gratefully acknowledged. The observed sea-level anomaly maps were obtained from AVISO. The 1992–2002 mean ocean dynamic topography data has been obtained from Nikolai Maximenko (IPRC) and Peter Niiler (SIO).

Edited by: L. Bertino

<figure><img src="image_47.png" /></figure>

INSU
Institut national des sciences de l'Univers

The publication of this article is financed by CNRS-INSU.

## References

Arakawa, A. and Lamb, V.: A potential enstrophy and energy conserving scheme for the shallow water equations, Mon. Weather Rev., 109, 18–36, 1981.

Barnier, B., Madec, G., Penduff, T., Molines, J.-M., Treguier, A.-M., le Sommer, J., Beckmann, A., Biastoch, A., Böning, C., Dengg, J., Derval, C., Durand, E., Gulev, S., Remy, E., Talandier, C., Theetten, S., Maltrud, M., McClean, J., and de Cuevas, B.: Impact of partial steps and momentum advection schemes in a global circulation model at eddy permitting resolution, Ocean Dynam., 56(5-6), 543–567. doi:10.1007/s10236-006-0082-1, 2006.

Beckmann, A., Böning, C., Köberle, C., and Willebrand, J.: Effects of increased horizontal resolution in a simulation of the North Atlantic Ocean, J. Phys. Oceanogr., 24, 326–344, 1994.

Berloff, P. S. and McWilliams, J. C.: Large-Scale, Low-Frequency Variability in Wind-Driven Ocean Gyres, J. Phys. Oceanogr., 29, 1925–1949, 1999.

Böning, C. and Bryan, F. O.: The Warmwatersphere of the North Atlantic Ocean, chap. Large-scale transport processes in high-resolution circulation models, 91–128, Gebrüder Borntraeger, 1996.
Böning, C. and Budich, R. G.: Eddy Dynamics in a Primitive Equation Model: Sensitivity to Horizontal Resolution and Friction, J. Phys. Oceanogr., 22, 361–381, 1992.

Brodeau, L., Barnier, B., Treguier, A.-M., and Penduff, T.: Comparing sea surface atmospheric variables from ERA40 and CORE with a focus on global net heat flux, 2006.

Brodeau, L., Barnier, B., Penduff, T., Treguier, A.-M., and Gulev, S.: An ERA40-based atmospheric forcing for global ocean circulation models, Ocean Model., 31(3-4), 88–104, 2009.

Capet, X., Marchesiello, P., and McWilliams, J.: Upwelling response to coastal wind profiles, Geophys. Res. Lett., 31, L13311, doi:10.1029/2004GL020123, 2004.

Chanut, J., Barnier, B., Large, W. G., Debreu, L., Penduff, T., Molines, J.-M., and Mathiot, P.: Mesoscale Eddies in the Labrador Sea and their Contribution to Convection and Restratification, J. Phys. Oceanogr., 38, 1617–1643, doi:10.1175/2008JPO3485.1, 2008.

Chelton, D. B., de Szoeke, R. A., Schlax, M. G., Naggar, K. E., and Siwertz, N.: Geographical variability of the first-baroclinic Rossby radius of deformation, J. Phys. Oceanogr., 28, 433–460, 1998.

Cravatte, S., Madec, G., Izumo, T., Menkes, C., and Bozec, A.: Progress in the 3-D circulation of the eastern equatorial Pacific in a climate ocean model, Ocean Model., 17, 28–48, 2007.

Dengg, J., Beckmann, A., and Gerdes, R.: The Gulf Stream separation problem, in: The Warmwatersphere of the North Atlantic Ocean, edited by: Krauss, W., 253–290, Gebrüder Borntraeger, Berlin, 1996.

Drillet, Y., Bourdallé-Badie, R., Siefridt, L., and Le Provost, C.: Meddies in the Mercator North Atlantic and Mediterranean Sea eddy-resolving model, J. Geophys. Res., 110, C03016, doi:10.1029/2003JC002170, 2005.

Duchon, C.: Lanczos Filtering in One and Two Dimensions, J. Appl. Meteorol., 18, 1016–1022, 1979.

DRAKKAR Group: Eddy-permitting ocean circulation hindcasts of past decades, CLIVAR Exchanges, 12, 8–10, 2007.

DYNAMO Group: DYNAMO: Dynamics of North Atlantic Models. Simulation and assimilation with high resolution models, Tech. rep., Ber. Inst. f. Meeresk. Kiel, 294, 1997.

Fanning, A. and Weaver, A.: A horizontal resolution and parameter sensitivity study of heat transport in an idealized coupled climate model, J. Climate, 10, 2469–2478, 1997.

Gent, P. R. and McWilliams, J. C.: Isopycnal Mixing in Ocean Circulation Models, J. Phys. Oceanogr., 20, 150–155, 1990.

Greatbatch, R. J.: A note on the representation of steric sea level in models that conserve volume rather than mass, J. Geophys. Res., 99(C6), 12767–12771, 1994.

Griffies, S. M., Biastoch, A., Böning, C., Bryan, F. O., Danabasoglu, G., Chassignet, E. P., England, M. H., Gerdes, R., Haak, H., Hallberg, R., Hazeleger, W., Jungclaus, J., Large, W. G., Madec, G., Pirani, A., Samuels, B., Scheinert, M., Gupta, A. S., Severijns, C., Simmons, H. L., Treguier, A.-M., Winton, M., Yeger, S., and Yin, J.: Coordinated ocean-ice reference experiments (COREs), Ocean Model., 26, 1–46, 2009.

Guinehut, S., Traon, P.-Y. L., and Larnicol, G.: What can we learn from global altimetry/hydrography comparisons?, Geophys. Res. Lett., 33, L10604, doi:10.1029/2005GL025551, 2006.

Gulev, S. K., Barnier, B., Molines, J.-M., Penduff, T., and Chanut, J.: Impact of spatial resolution on simulated surface water mass

transformations in the Atlantic, Ocean Model., 19, 138–160, 2007.

Haidvogel, D. B., Arango, H. G., Hedstrom, K., Beckmann, A., Malanotte-Rizzoli, P., and Shchepetkin, A. F.: Model evaluation experiments in the North Atlantic Basin: simulations in nonlinear terrain-following coordinates, Dyn. Atmos. Oceans, 32, 239–281, 2000.

Hecht, M. and Hasumi, H.: Ocean Modeling in an Eddying Regime, 177, AGU Monograph Series, 2008.

Hecht, M. W. and Smith, R. D.: Ocean Modeling in an Eddying Regime, chap. Towards a Physical Understanding of the North Atlantic: A Review of Model Studies in an Eddying Regime, Geophysical Monograph Series, American Geophysical Union, 2008.

Hogg, A. McC. and Blundell, J.: Interdecadal Variability of the Southern Ocean, J. Phys. Oceanogr., 36(8), 1626–1645, doi:10.1175/JPO2934.1, 2006.

Holland, W., Harrison, D., and Semtner, A.: Eddies in Marine Science, chap. Eddy-resolving numerical models of large-scale ocean circulation, 379–403, Springer-Verlag, 1983.

Jiang, N., Jin, F. F., and Ghil, M.: Multiple equilibria, periodic and aperiodic solutions in a wind-driven, double-gyre, shallow-water model, J. Phys. Oceanogr., 25, 764–786, 1995.

Kelly, K. A., Thompson, L., Cheng, W., and Metzger, E. J.: Evaluation of HYCOM in the Kuroshio Extension region using new metrics, J. Geophys. Res., 112, 2007.

Large, W. G. and Yeager, S. G.: Diurnal to decadal global forcing for ocean and sea-ice models, NCAR technical note, 22 pp., 2004.

Le Sommer, J., Penduff, T., Theetten, S., Madec, G., and Barnier, B.: How momentum advection schemes influence current-topography interactions at eddy-permitting resolution, Ocean Modelling, 29, 1–14, doi:10.1016/j.ocemod.2008.11.007, 2009.

Levitus, S., Boyer, T., Conkright, M., O'Brian, T., Antonov, J., Stephens, C., Stathopolos, L., Johnson, D., and Gelfeld, R.: World Ocean Database 1998, Tech. Rep. NESDID18, NOAA Atlas, 1998.

Lique, C., Treguier, A.-M., Scheinert, M., and Penduff, T.: A model-based study of ice and freshwater transport variability along both sides of Greenland, Clim. Dynam., 33, 5, 685–705, doi:10.1007/s00382-008-0510-7, 2009.

Lombard, A., Garric, G., and Penduff, T.: Regional patterns of observed sea level change: insights from a 1/4° global ocean/sea-ice hindcast, Ocean Dynam., 59, 3, 433–449, doi:10.1007/s10236-008-0161-6, 2009.

Madec, G.: NEMO reference manual, ocean dynamics component : NEMO-OPA. Preliminary version, Tech. Rep. 27, Note du Pôle de modélisation, Institut Pierre-Simon Laplace (IPSL), France, ISSN No 1288-1619, 2008.

Masumoto, Y., Sasaki, H., Kagimoto, T., Komori, N., Ishida, A., Sasai, Y., Miyama, T., Motoi, T., Mitsudera, H., Takahashi, K., Sakuma, H., and Yamagata, T.: A fifty-year eddy-resolving simulation of the world ocean: Preliminary outcomes of OFES (OGCM for the Earth Simulator), J. Earth Simulator, 1, 35–56, 2004.

McClean, J., Poulain, P.-M., Pelton, J. W., and Maltrud, M.: Eulerian and Lagrangian Statistics from Surface Drifters and a High-Resolution POP Simulation in the North Atlantic, J. Phys. Oceanogr., 32, 2472–2491, 2002.
Milliff, R. F., Large, W. G., Holland, W. R., and McWilliams, J. C.: The general circulation responses of high-resolution North Atlantic Ocean models to synthetic scatteromter winds, J. Phys. Oceanogr., 26, 1747–1768, 1996.

Molines, J.-M., Barnier, B., Penduff, T., Brodeau, L., Treguier, A.-M., Theetten, S., and Madec, G.: Definition of the interannual experiment ORCA025-G70, 1958–2004, Tech. Rep. LEGI-DRA-2-11-200, LEGI, available at http://www.ifremer.fr/lpo/drakkar, 2006.

Niiler, P. P., Maximenko, N. A., and McWilliams, J. C.: Dynamically balanced absolute sea level of the global ocean derived from near-surface velocity observations, Geophys. Res. Lett., 30, 22, OCE7.1-OCE7.4, doi:10.1029/2003GL018628, 2003.

Penduff, T., Le Sommer, J., Barnier, B., Treguier, A.-M., Molines, J.-M., and Madec, G.: Influence of numerical schemes on current-topography interactions in 1/4° global ocean simulations, Ocean Sci., 3, 509–524, 2007, http://www.ocean-sci.net/3/509/2007/.

Roberts, M., Banks, H., Gedney, N., Gregory, J., Hill, R., Muller-worth, S., Pardaens, A., Rickard, G., Thorpe, R., and Wood, R.: Impact of an Eddy-Permitting Ocean Resolution on Control and Climate Change Simulations with a Global Coupled GCM, J. Climate, 17, 3–20, 2004.

Scott, R. B. and Arbic, B. K.: Spectral Energy Fluxes in Geostrophic Turbulence: Implications for Ocean Energetics, J. Phys. Oceanogr., 37, 673–688, 2007.

Roullet, G. and Madec, G.: Salt conservation, free surface, and varying levels: a new formulation for ocean general circulation models, J. Geophys. Res., 105, C10, 23927–23942, 2000.

Semtner, A. and Chervin, R.: A simulation of the global ocean circulation with resolved eddies, J. Geophys. Res., 93, 15502–15522, 1988a.

Semtner, A. and Chervin, R.: Ocean general circulation from a global eddy-resolving model, J. Geophys. Res., 97, 5493–5550, 1988b.

Smith, R., Maltrud, M., Bryan, F., and Hecht, M.: Numerical simulation of the North Atlantic Ocean at 1/10°, J. Phys. Oceanogr., 30, 1532–1561, 2000.

Spall, M. A.: Dynamics of the Gulf Stream/Deep Western Boundary Current Crossover, Part II: Low-frequency internal oscillations., J. Phys. Oceanogr., 26, 2169–2182, 1996.

Street, J. O., Carroll, R. J., and Ruppert, D.: A Note on Computing Robust Regression Estimates Via Iteratively Reweighted Least Squares, The American Statistician, 42, 152–154, 1988.

Timmermann, A., Goosse, H., Madec, G., Fichefet, T., Ethe, C., and Dulière, V.: On the representation of high latitude processes in the ORCA-LIM global coupled sea-ice ocean model, Ocean Model., 8, 175–201, 2005.

Treguier, A.-M., Theetten, S., Chassignet, E., Penduff, T., Smith, R., and Talley, L.: The North Atlantic subpolar gyre in four high resolution models, J. Phys. Oceanogr., 35, 757–774, 2005.

Treguier, A. M., England, M. H., Rintoul, S. R., Madec, G., Le Sommer, J., and Molines, J.-M.: Southern Ocean overturning across streamlines in an eddying simulation of the Antarctic Circumpolar Current, Ocean Sci., 3, 491–507, 2007, http://www.ocean-sci.net/3/491/2007/.

Vinogradova, N. T., Ponte, R. M., and Stammer, D.: Relation between sea level and bottom pressure and the vertical dependence of oceanic variability, Geophys. Res. Lett., 34, L03608, doi:10.1029/2006GL028588, 2007.

von Storch, H. and Zwiers, F. W.: Statistical Analysis in Climate Research, Cambridge University Press (ISBN 0 521 45071 3), 484 pp., 1999.

Zhai, X., Greatbatch, R. J., and Sheng, J.: Diagnosing the role of eddies in driving the circulation of the northwest Atlantic Ocean, Geophys. Res. Lett., 31, 1–4, 2004.
# Sea Level Expression of Intrinsic and Forced Ocean Variabilities at Interannual Time Scales

THIERRY PENDUFF,*,+ MÉLANIE JUZA,* BERNARD BARNIER,* JAN ZIKA,* WILLIAM K. DEWAR,+ ANNE-MARIE TREGUIER,‡ JEAN-MARC MOLINES,* AND NICOLE AUDIFFREN@

*Laboratoire des Ecoulements Géophysiques et Industriels, UMR CNRS 5519, Université de Grenoble, Grenoble, France

+Earth, Ocean and Atmospheric Science Department, The Florida State University, Tallahassee, Florida

‡Laboratoire de Physique des Océans, UMR CNRS 6523, IFREMER, UBO, Plouzané, France

@Centre Informatique National de l'Enseignement Supérieur, Montpellier, France

(Manuscript received 5 January 2011, in final form 19 May 2011)

## ABSTRACT

This paper evaluates in a realistic context the local contributions of direct atmospheric forcing and intrinsic oceanic processes on interannual sea level anomalies (SLAs). A 1/4° global ocean-sea ice general circulation model, driven over 47 yr by the full range of atmospheric time scales, is quantitatively assessed against altimetry and shown to reproduce most observed features of the interannual SLA variability from 1993 to 2004. Comparing this simulation with a second driven only by the climatological annual cycle reveals that the intrinsic part of the total interannual SLA variance exceeds 40% over half of the open-ocean area and exceeds 80% over one-fifth of it. This intrinsic contribution is particularly strong in eddy-active regions (more than 70%-80% in the Southern Ocean and western boundary current extensions) as predicted by idealized studies, as well as within the 20°-35° latitude bands. The atmosphere directly forces most of the interannual SLA variance at low latitudes and in most midlatitude eastern basins, in particular north of about 40°N in the Pacific. The interannual SLA variance is almost entirely due to intrinsic processes south of the Antarctic Circumpolar Current in the Indian Ocean sector, while half of this variance is forced by the atmosphere north of it. The same simulations were performed and analyzed at 2° resolution as well: switching to this laminar regime yields a comparable forced variability (large-scale distribution and magnitude) but almost suppresses the intrinsic variability. This likely explains why laminar ocean models largely underestimate the interannual SLA variance.

## 1. Introduction

### a. Context

The dynamical response of the nonlinear ocean to the atmospheric variability is complex: dynamical system approaches and diagnostics of potential vorticity balances have demonstrated that in presence of stochastic, seasonal, or constant atmospheric forcing, oceanic nonlinearity can intrinsically generate and sustain variability of the horizontal circulation at interannual time scales (e.g., Jiang et al. 1995; Spall 1996; Berloff and McWilliams 1999; Hazeleger and Drijfhout 2000; Dewar 2001; Cessi and Louazel 2001; Wirth et al. 2002; Dewar 2003; Simonnet et al. 2003, 2005; Hogg and Blundell 2006; Berloff et al. 2007). There are various manifestations of

this variability, including, for example, fluctuations of inertial gyres and associated mode water pools, western boundary current transitions between high- and low-energy states (studied in closed basins), slow vacillations of turbulent jets, and eddy-induced modulations of topographic steering (studied in periodic channels). In this study we focus on the interannual variability of the near-surface horizontal circulation [see section 2 in Dijkstra and Ghil (2005) for a review of idealized results], which differs from the (slower) variability of the thermohaline circulation (their section 3).

The studies mentioned above made use of numerical models with simplified dynamics (e.g., quasigeostrophic or shallow water, often flat bottomed and/or barotropic) and geometries, resolving a limited range of processes. Such idealized approaches are necessary to perform detailed analyses of key interactions and to demonstrate their possible existence in the real ocean. These results provide important dynamical insights, but they cannot

<footnote>Corresponding author address: Thierry Penduff, LEGI-MEOM, BP53, 38041 Grenoble CEDEX 9, France.
E-mail: thierry.penduff@legi.grenoble-inp.fr</footnote>
be directly compared with actual observations regarding
the distribution and magnitude of the intrinsic variabil-
ity in the real ocean.

Regional and basin-scale manifestations of intrinsic
low-frequency (LF) variability have also been reported
in seasonally forced eddying simulations based on ocean
general circulation models (OGCMs) that solve the prim-
itive equations within realistic geometries with realistic
initial states. These studies have revealed significant
interannual-to-decadal modulations of, for example,
basin-averaged mesoscale activity (Penduff et al. 2004),
interygre heat transport (Hall et al. 2004), sea level
variability (Cabanes et al. 2006), expansion-contraction
of western boundary current systems (Taguchi et al.
2007), and meridional overturning circulation (Biastoch
et al. 2008).

Both idealized and more “realistic” model studies
suggest that at high Reynolds number and with a suffi-
ciently strong atmospheric forcing, the oceanic variability
is likely a complex combination of direct responses to
the atmosphere, and of intrinsic fluctuations (Dijkstra
and Ghil 2005). However, the interplay between forced
and intrinsic variabilities and their relative contributions
to the observed interannual variability remain unclear.
For instance, the Kuroshio undergoes realistic inter-
annual variability patterns both in wind-driven laminar
models (Qiu 2003) and in constantly forced eddying
models (Pierini 2006); Taguchi et al. (2007) suggests
from more realistic simulations that the observed vari-
ability arises from the atmospheric forcing of intrinsic
variability modes. In other words, the contributions of
forced and intrinsic interannual variabilities to the ob-
served signals and their possible interactions may largely
depend on location, but they have been mostly investi-
gated regionally and in idealized frameworks so far.

b. *This study*

To our knowledge, no dedicated mapping and de-
scription of intrinsic and forced variabilities have been
yet attempted at global scale with comparisons to ob-
servations. This paper addresses three main issues in the
(still simplified) context of global ocean–sea ice primi-
tive equation dynamics with realistic geometry, and
from the perspective of the (well observed) sea level
variability:

* What are the contributions of the forced and intrinsic components in the interannual sea level anomaly (SLA) variability?

* How are these contributions distributed spatially at global scale?

* What is the sensitivity of intrinsic interannual vari-
ability to mesoscale activity?

Our aim is to qualitatively (and perhaps quantita-
tively) evaluate the relative contributions of forced and
intrinsic interannual variabilities in the global ocean as
measured by altimeters (at the expense of detailed
process analyses). This complements process-oriented
studies where individual dynamical processes are stud-
ied in detail in highly idealized contexts (at the expense
of direct comparisons with actual observations). We first
evaluate against altimeter observations a “pivot” exper-
iment, that is, a state-of-the-art, 47-yr eddy-admitting
(1/4°) ocean–sea ice global simulation forced by a real-
istic interannual forcing. This evaluation follows the
approach proposed by Penduff et al. (2010, hereafter
P10). The same model is then driven by the mean annual
cycle only to isolate the intrinsic interannual compo-
nent; the variance of the forced component is then es-
timated from the latter two simulations. To illustrate the
possible contribution of nonlinearities on these com-
ponents, the latter comparison is also made between
coarse-resolution (2°), laminar versions of the interan-
nually and seasonally forced experiments.

Archiving, Validation, and Interpretation of Satellite
Oceanographic data (AVISO) is the only existing dataset
that allows a relatively long-term, global, and local de-
scription of the range of scales we are focused on. We
shall therefore focus our analyses on SLAs. The detailed
investigation of unobserved variables and underlying
processes in individual regions and simulations would
require, in turn, a number of idealized studies, which are
left for the future.

Section 2 presents the 2° and 1/4° model configurations,
the interannual and seasonal forcing functions, the
construction of our collocated multimodel and observed
SLA dataset, and comparison procedures. Our pivot
experiment (1/4°, interannual forcing) is assessed against
SLA observations in section 3. In section 4, we compare
this simulation with its seasonally forced counterpart to
describe the features of intrinsic and forced interannual
SLA variabilities at global scale. Both 2° simulations are
compared in section 5 to estimate the possible contri-
bution of nonlinearities at finer resolution. Conclusions
are given and discussed in section 6.

**2. Models, datasets, and comparison procedures**

a. *Model simulations*

The six global ocean–sea ice model simulations men-
tioned in this paper were performed in the framework of
the Drakkar project (http://www.ifremer.fr/lpo/drakkar/),
using the Nucleus for European Modelling of the Ocean
(NEMO) code (Madec et al. 2008). Figure 1 summa-
rizes the relationships between these simulations. Our
sensitivity studies to the forcing were performed using
<figure><img src="image_48.png" /><figcaption>FIG. 1. Summary of the model simulations mentioned in this study, including P10 simulations and those in bold form the core of the present paper. Note that all model parameters except atmospheric fields and bulk formulas remain unchanged within each column and everything except resolution and subgrid-scale parameterizations remain unchanged within each row. The three curved arrows, and associated descriptions below the matrix, summarize configurations changes between successive lines or columns.</figcaption></figure>

eddy-admitting (1/4°) and laminar (2°) configurations. Their forcings aside, these configurations are exactly the same as those described in P10; they differ from each other in terms of horizontal resolution (2° and 1/4°) and subgrid-scale parameterizations. Note that the meridional resolution of the 2° grid gets finer toward the equator, where it reaches 1/3°; this explains why the 1/4° and 2° simulations do not differ much in this region. The 1/4° and 2° configurations are driven by two forcing functions, yielding a total of four simulations.

## 1) *T* FORCING

The first forcing function, labeled *T* since it drives the *total* low-frequency variability (forced plus intrinsic), is based on reanalyzed atmospheric fields and satellite data, and covers the period 1958–2004. It is described in detail [and referred to as Drakkar forcing set 4.1 (DFS4.1)] in Brodeau et al. (2010). In short, DFS4 is based on satellite-derived monthly precipitations, daily radiative heat fluxes, and 6-hourly 10-m atmospheric state variables from European Centre for Medium-Range Weather Forecasts (ECMWF) fields [40-yr ECMWF Re-Analysis (ERA-40) before 2002, ECMWF analysis

afterward]. Turbulent air-sea and air-ice fluxes are computed from ECMWF variables and surface model variables through bulk formulas. This forcing function is used to drive the *T* − 2° and *T* − 1/4° simulations.

To keep the paper concise, the assessment of *T* − 2° and *T* − 1/4° against altimetry will involve comparisons with two earlier simulations, labeled *T* − 2°<sub>P10</sub> and *T* − 1/4°<sub>P10</sub>, respectively, in the following. These were extensively described and evaluated in P10. Simulations *T* − 2° and *T* − 1/4° only differ from these earlier simulations by two features of their forcing function. First, the atmospheric fields (DFS3) used to drive the P10 runs have been updated in *T* − 2° and *T* − 1/4°.<fnref n="1" /> Second, the bulk formulas in *T* − 2° and *T* − 1/4° (as well as in *I* − 2° and *I* − 1/4°, described below) do not take into account surface currents in the computation of wind stress,

<footnote n="1">¹ The calibration of atmospheric variables has been improved from DFS3 to DFS4, and some spurious discontinuities have been removed. The rationale, description, and oceanic impact of individual changes made from DFS3 to DFS4 are discussed in detail in Brodeau et al. (2010).</footnote>
unlike in P10 simulations. This current-stress feedback has been discussed since Pacanowski (1987), but there is still no clear consensus yet on whether and how it should be parameterized in numerical models [see in particular Eden and Dietze (2009), and references therein]. To maximize the consistency between this *T* forcing and the seasonal *I* forcing needed for this study, we chose to ignore this latter feedback in the computation of wind stresses (see the section below).

## 2) I FORCING

The second forcing function was built to preserve the daily mean annual cycle and long-term average of the *T* forcing, but it is devoid of interannual variability; it is labeled *I* since only *intrinsic* low-frequency variability may be excited. This forcing was built by computing 365 daily averages (i.e., a daily mean annual climatology) of 6-hourly DFS4 variables during 50 yr. This daily annual climatology includes precipitations, runoff, radiative fluxes, atmospheric state variables (air temperature, air humidity, wind speed vector *u*<sub>air</sub>), and the main quadratic contributions to air-sea fluxes: bulk exchange coefficients (*C*<sub>d</sub> for momentum, *C*<sub>e</sub> for latent fluxes, and *C*<sub>h</sub> for sensible heat fluxes) and pseudowind stresses *C*<sub>d</sub>|*u*<sub>air</sub>|*u*<sub>air</sub>. As mentioned above, ocean currents are not used to calculate wind stresses. The resulting *I* forcing fields were then low-pass filtered in time (three-point Hanning filter) to remove the remaining, although small, high-frequency (HF) noise. This *I* forcing function is applied for 314 yr on the 2° and 1/4° configurations. It is devoid of any interannual atmospheric variability and of intraseasonal synoptic transients. The corresponding simulations will be referred to as *I* − 2° and *I* − 1/4°, respectively.

This method ensures that the linear and the main quadratic air-sea flux terms (stability of the atmospheric boundary layer, exchange coefficients, wind stresses, turbulent heat and freshwater fluxes, etc.) have the same seasonal cycles and long-term means in the *T* and *I* forcing functions. In other words, this method maximizes the consistency between the interannual and seasonal forcing functions, and between our simulations.

To summarize, the four new simulations investigated in this paper<fnref n="2" /> are *T* − 2°, *T* − 1/4°, *I* − 2°, and *I* − 1/4°. When joined with the two P10 simulations, we have a total of six global runs. Model fields, including sea surface height (SSH), were saved as successive 5-day averages throughout both interannually forced 47-yr simulations.

Model fields were saved at the same 5-day rate after year 300 in both 314-yr seasonally forced simulations but as monthly averages before this. Our choice of outputs and associated constraints are described in appendix A.

## b. Sea level datasets

The collocation and filtering process of simulated and observed SLAs is presented in detail in P10 (their section 2.2.1); only a short summary is given here. Our observational reference is the 1993–2004 series of altimeter sea level anomalies, provided by AVISO as a series of global weekly maps on a 1/3° × 1/3° Mercator grid. Simulated SSHs are converted into SLA maps after a time-space collocation onto this AVISO dataset. The resulting dataset consists of weekly maps of observed (AVISO) and simulated SLAs available on the same 1/3° Mercator grid, on the same dates (the case of seasonally forced simulations is presented in appendix A). This seven-member (AVISO + six simulations) “raw” SLA dataset is then filtered in time using a Lanczos filter. As done and described in P10, each member of this raw dataset is split into its “low frequency” (time scales longer than 18 months) and “high frequency” (time scales shorter than 5 months) parts.<fnref n="3" /> In the present study, linear trends are removed at each grid point from each LF dataset (see appendix B).

## c. Comparison statistics

Our purpose is to characterize at high and low frequency (the main focus being on the latter) the distribution and magnitude of SLA variability in our seven collocated sources, that is, AVISO, the four new simulations, and the two P10 simulations. The processing summarized below is described in detail in P10 (their section 2.2.2).

### 1) $\sigma$, CT, AND CS

For each of these 14 cases (7 datasets × 2 frequency bands), we compute 2D fields of SLA temporal standard deviations [$\sigma(i, j)$]. The $\sigma$ maps are shown in log scale as gray shadings (low frequency) and contours (high frequency) for AVISO and the four recent runs in Figs. 2 and 3 after spatial smoothing for clarity. The reader is referred to P10’s Fig. 4 for $\sigma$ maps in the $T - 2_{P10}^\circ$ and $T - 1_{P10}^\circ$ simulations.

For both frequency bands and each simulation (12 cases), we compute maps of correlation coefficients $C_t(i, j)$

<footnote n="2">² For clear reference, these four simulations are referred to as ORCA246-G83b, ORCA025-B83, ORCA246-MJM01, and ORCA025-MJM01 in the Drakkar database.</footnote>

<footnote n="3">³ LF and HF frequency bands will also be referred to as “interannual” and “mesoscale,” respectively, as in P10. Time scales between 5 and 18 months (referred to as *quasi annual* in P10) are barely sensitive to the shutoff of interannual forcing and are not investigated in the present study.</footnote>
<figure><img src="image_49.png" /><figcaption>FIG. 2. Total variability (log scale, cm) at LF (gray shading) and HF (selected contours: 0.2 in blue, 0.5 in red, 0.8 in orange, 1.1 in yellow): (top) AVISO, and fully forced model simulations at (middle) 1/4° and (bottom) 2°. For clarity, maps have been smoothed.</figcaption></figure>

between (observed and simulated) collocated local SLA time series. Local significance levels at 95% are estimated at each grid point for each Ct(i, j) map from lagged autocorrelations functions of observed and simulated SLA time series, as proposed by von Storch and Zwiers (1999). As noted by P10, these significance level estimates are relevant for our time series, which have different (and arbitrary) spectra and autocorrelations. Insignificant Ct(i, j) values are masked in the following: only significant Ct(i, j) values are considered in further computations and discussed below.

Statistics will also be presented within 14 latitude bands labeled $\lambda \in (1; 14)$, spanning the 70°S–70°N latitude range by 10° intervals (this interval was 5° in P10). To this aim, the 14 $\sigma(i, j)$ and 12 Ct(i, j) maps are simply area weighted and averaged over each band to yield $\sigma(\lambda)$ and Ct(λ) in each frequency band. The spatial correlation Cs(λ) between simulated $\sigma(i, j)$ maps and their AVISO counterpart is computed within each latitude band for all 12 cases to evaluate the agreement between both variability maps as a function of latitude. These latitude-dependent skill estimates are presented for the earlier
<figure><img src="image_50.png" /><figcaption>FIG. 3. As in Fig. 2, but for intrinsic variability obtained from seasonally forced simulations (top) $I - 1/4^\circ$ and (bottom) $I - 2^\circ$.</figcaption></figure>

$(T - 1/4^\circ)_{P10}, T - 2^\circ)$ and new $(T - 1/4^\circ, T - 2^\circ)$ interannually forced simulations in Fig. 5, and for the interannually forced $(T - 1/4^\circ$ and $T - 2^\circ)$ and seasonally forced $(I - 1/4^\circ$ and $I - 2^\circ)$ new simulations in Fig. 6.

The agreement between (unsmoothed) LF and HF variability maps within each individual dataset is evaluated from spatial correlations in successive latitude bands, and shown in Fig. 4 for AVISO and for the four new simulations.

## 2) INTRINSIC AND FORCED LF VARIABILITY

We shall consider in this study that at low frequency, the total SLA variability $T(t)$ locally found in interannually forced simulations $T$ is the sum of the intrinsic component $I(t)$ (diagnosed from $I$ simulations) and of a forced variability component $F(t)$. This simple assumption is tacitly made in idealized studies, where the low-frequency variability simulated under steady or seasonal forcing is considered representative of the intrinsic component of the oceanic variability. This yields the following relationship between the total, intrinsic, and forced local LF variances:

$$ (\sigma^T)^2 = (\sigma^I)^2 + (\sigma^F)^2 + 2\sigma^I\sigma^F C_{I,F}, \quad (1) $$

where $C_{I,F}$ denotes the local temporal correlation between intrinsic and forced SLA fluctuations. As mentioned in the introduction, idealized studies have shown that $I(t)$ has a strong random character; in other words, we also assume at this point that the local time series of $I(t)$ and $F(t)$ at low frequency are not significantly correlated in time (i.e., $C_{I,F} \sim 0$; the validity of this convenient, initial assumption<fnref n="4" /> is assessed and discussed in section 4c(1)]. This second assumption allows us to estimate the local variance of forced low-frequency SLA variabilities as the difference between total and intrinsic LF variances. The forced LF standard deviation of SLA

<footnote n="4">⁴ This assumption was explicitly made by Hazeleger and Drijfhout (2000, p. 13 966), who wrote, “By definition, externally generated oceanic variability is generated independently of oceanic processes, for instance, the passive response of the ocean to variability in the atmospheric forcing.”</footnote>
<figure><img src="image_51.png" /><figcaption>FIG. 4. Skills of earlier ($T - 1/4_{P10}$ and $T - 2_{P10}$, dashed lines) and new ($T - 1/4$ and $T - 2$, solid lines) fully forced simulations as a function of latitude (abscissas) in the (left) HF and (right) LF frequency bands. (top row) Total SLA standard deviation [$\sigma^T(\lambda)$, cm]; AVISO reference is shown in green. (middle row) Significant temporal correlations with local AVISO time series [Ct($\lambda$)]. (bottom row) Spatial correlations [Cs($\lambda$)] between simulated and observed $\sigma^T(i, j)$ maps.</figcaption></figure>

in the eddying regime may thus be estimated at every grid point by

$$ \sigma_{LF}^{F-1/4^\circ} = \sqrt{(\sigma_{LF}^{T-1/4^\circ})^2 - (\sigma_{LF}^{I-1/4^\circ})^2} $$

This latter quantity, and its laminar counterpart $\sigma_{LF}^{F-2^\circ}$ (derived from $\sigma_{LF}^{T-2^\circ}$ and $\sigma_{LF}^{I-2^\circ}$), are shown as maps and in zonal average in Fig. 7. Finally, the percentage of the total LF variance that is generated intrinsically is estimated as

$$ R_{LF}^{I} = 100\% \times (\sigma_{LF}^{I} / \sigma_{LF}^{T})^{2}. $$

These ratios are shown as maps and zonal averages in Fig. 8 for the $1/4^\circ$ and $2^\circ$ resolutions.

### 3. Assessment of SLA variability in the pivot experiment $T - 1/4^\circ$

We first assess with respect to AVISO the geographical distribution, magnitude, and phase of the SLA variability in the $1/4^\circ$ interannually forced simulation ($T - 1/4^\circ$) between 1993 and 2004. The (smoothed) spatial distributions of high- and low-frequency variability in both datasets are shown in the top and middle panels of Fig. 2. Besides the recalibration of its forcing (see Fig. 1), the
<figure><img src="image_52.png" /><figcaption>FIG. 5. Spatial correlations by latitude band (abscissas) between unsmoothed HF and LF $\sigma(i, j)$ maps for each dataset (colors); large values mean large correspondence between spatial structures of $\sigma$ in both frequency bands.</figcaption></figure>

$T - 1/4°$ simulation is the same as $T - 1/4_{P10}$ which was extensively compared to the same AVISO referenced by P10. For conciseness, the fidelity of $T - 1/4°$ with respect to AVISO is not discussed as extensively here but is compared to $T - 1/4_{P10}$ in zonal average (plain and dashed black lines in Fig. 5); in addition to this figure, the interested reader is referred to P10 for an in-depth description of the strengths and weaknesses of $T - 1/4_{P10}$.

### a. HF variability at 1/4°

The middle left panel in Fig. 5 shows in zonal average that the temporal correlations between observed and local HF time series are almost identical in $T - 1/4°$ and $T - 1/4_{P10}$. However, the transition from $T - 1/4_{P10}$ to $T - 1/4°$ enhances the magnitude of the SLA HF variability by 10% to almost 100% at most latitudes in zonal average (cf. dashed and solid black lines, top left panel),<fnref n="5" /> bringing it much closer to the AVISO reference. Despite this large increase, the HF variability in $T - 1/4°$ remains weaker than observed at midlatitudes (cf. black and green lines). As noted by P10 indeed, $1/4°$ eddy admitting models like ours cannot resolve more than the first deformation radius away from the lowest latitudes; hence, they may underestimate the midlatitude eddy activity generated by baroclinic instability.

P10 have shown that the HF surface variability in $T - 1/4_{P10}$ was correctly distributed in most regions of the global ocean; at most latitude ranges (black lines, bottom left panel in Fig. 5), this distribution is further improved in $T - 1/4°$. However, certain currents in both

<figure><img src="image_53.png" /><figcaption>FIG. 6. Skills of the simulations driven by the full DFS4 forcing ($T - 1/4°$ in black, $T - 2°$ in red) and by the seasonal forcing ($I - 1/4°$ in gray, $I - 2°$ in pink) in the LF band as a function of latitude (abscissas). Rows correspond to the same quantities as in Fig. 5. AVISO reference is shown in green.</figcaption></figure>

simulations are locally shifted compared to observations (Fig. 2): the North Atlantic Current (NAC) locally follows the mid-Atlantic Ridge, the Kuroshio does not extend far enough to the east, and the Antarctic Circumpolar Current (ACC) undergoes a localized meridional shift north of the Kerguelen Plateau. Details about these local biases, which are often present in higher-resolution simulations as well, may be found in Barnier et al. (2006).

To summarize at this point, the high-frequency SLA variability in $T - 1/4°$ is typical of eddy-admitting solutions with a few local discrepancies and is more realistic than in $T - 1/4_{P10}$ in terms of intensity. The reader is referred to P10 for a more detailed description of these local biases.

<footnote n="5">⁵ This large increase is consistent with Eden and Dietze's (2009) findings about the suppression of the surface current feedback onto the wind stress.</footnote>
<figure><img src="image_54.png" /><figcaption>FIG. 7. Maps: Forced variability at LF $\sigma_{LF}^F$ (log scale, cm) at (top) 1/4° and (middle) 2° res- olutions. Contours in (top) enclose areas where the intrinsic SLA variance exceeds the total SLA variance by 20% at 1/4° resolution (this does not happen at 2°); see text for details. (bottom) Zonally averaged LF forced variability $\sigma_{LF}^F$ as a function of latitude in the 1/4° (black) and 2° (red) cases.</figcaption></figure>

b. LF variability at 1/4°

The dashed and solid black lines in the right panels in Fig. 5 show at low frequencies how the transition from $T - 1/4^\circ_{P10}$ to $T - 1/4^\circ$ modifies the realism of SLA variabilities. More details about this variability in $T - 1/4^\circ_{P10}$ may be found in P10. The low-frequency SLA variability undergoes a 10%–30% enhancement, which brings its zonally averaged meridional profile remarkably close to the AVISO reference (top right panel). Despite a very small decrease in temporal correlations (middle right panel), the geographical distribution of the LF variability at midlatitudes (bottom right panel) is also improved

compared to $T - 1/4^\circ_{P10}$. In other words, the interannual SLA variability is globally improved and compares well with AVISO, in particular in terms of amplitude.

c. Geographical correspondence between HF and LF variabilities at 1/4°

HF and LF variability maps from each dataset are superimposed in Figs. 2 and 3; their mutual correspon- dence (e.g., spatial correlations within each dataset) is shown as a function of latitude in Fig. 4. Previous studies have shown [see, e.g., the introduction in Chang et al. (2001)], and the AVISO data confirm (top panel in Fig. 2,
$R_{LF}^I$ (%): LF VARIANCE EXPLAINED BY INTRINSIC PROCESSES

<figure><img src="image_55.png" /><figcaption>FIG. 8. Maps: Percentage of the LF variance explained by intrinsic processes ($R_{LF}^I$, shading) at (top) 1/4° and (bottom) 2°. Contours show the LF intrinsic variability in log scale (same colors as in Fig. 2). (bottom) Same percentages averaged zonally and shown as a function of latitude in the 1/4° (red) and 2° (black) cases.</figcaption></figure>

green line in Fig. 4), that most extratropical regions with strong mesoscale variability (HF in contours) also exhibit strong variability at longer time scales, that is, in the interannual (LF, shading) and quasi-annual (not shown; see P10) bands. This is particularly clear along the Gulf Stream (GS), NAC, Kuroshio, ACC, in the Brazil–Malvinas Confluence, in the Agulhas retroflection region, and in the East Australian Current area. P10 showed that a significant part of this localized low-frequency variability is found at scales larger than about 6°.

This observed and simulated correspondence denotes low-frequency fluctuations of both mesoscale activity and regional circulation patterns in the same regions, with probable interactions between both. Namely, eddy

energy is drawn from—hence is expected to follow the slow changes of—the regional circulation; in turn, eddies are likely to feed the slow and large-scale variability locally through, for example, Reynolds stress (Berloff et al. 2007) or inverse cascade processes (e.g., Scott and Arbic 2007). Note that the observed maxima of HF and LF variability are much less collocated at low latitudes (spatial correlation falls to 0.15, green line in Fig. 4). This is consistent with the more linear dynamics at work in the equatorial waveguide, where fast and slow waves are more likely to evolve independently than at higher latitudes and to reach maximum amplitudes at different locations.

These observed features are well reproduced in $T - 1/4°$ (middle panel in Fig. 2, black line in Fig. 4). This realistic
correspondence between HF and LF SLA variabilities further confirms the skill of this simulation; its sensitivity to the forcing is further discussed below.

To summarize this section 3, the SLA variability at all time scales in the new $T - 1/4°$ simulation is comparable to its earlier $T - 1/4°_{P10}$ version (extensively described in P10) in terms of phase and spatial distribution, and is much closer to altimetric observations in terms of magnitude. Except in localized areas, our $1/4°$ model reproduces many features of the real ocean's SLA variability under realistic forcing, in particular at interannual time scales. Its interannual variability without direct interannual forcing is now examined.

## 4. Intrinsic and forced low-frequency variability at $1/4°$

We now describe the sea level interannual variability that is generated intrinsically at $1/4°$ under seasonal forcing (in the $I - 1/4°$ simulation). The forced interannual variability and its contribution relative to the intrinsic component are presented afterward. Note that our diagnostics are local but that the interannual variability might be forced (or intrinsically generated) remotely in space and previously in time. The shading in Figs. 2 and 3 show the maps of total $\sigma^T$ and intrinsic $\sigma^I$ variabilities, respectively. Figure 6 shows the zonally averaged profiles of $\sigma$, Ct, and Cs for the total and intrinsic variabilities (in black and gray, respectively).

### a. Intrinsic LF variability at $1/4°$

Zonally averaged, significant temporal correlations Ct between low-frequency SLA time series in AVISO and in $I - 1/4°$ (gray line, middle panel of Fig. 6) are much smaller than in the $T - 1/4°$ case (black). This means that the phase of the interannual variability is significantly different from observed, consistent with its intrinsic generation. Note, however, that the decrease in Ct due to the shutoff of interannual forcing is small in the Southern Ocean (where Ct is already small with interannual forcing).

### 1) SPATIAL DISTRIBUTION

The top panel in Fig. 6 shows that the intrinsic interannual variability in $I - 1/4°$ (gray line) is very small at low latitudes. Without direct forcing indeed, our numerical ocean cannot generate the main modes of interannual variability that largely involve air-sea coupling in these regions: the El Niño-Southern Oscillation (ENSO; see Philander 1990) or the Indian dipole (e.g., Smith 2000). Relatively linear low-latitude ocean dynamics are also unlikely to locally generate low frequencies through inverse cascade processes. The top panel in Fig. 3 suggests

that this weak intrinsic interannual variability may be generated remotely—for example, by the North Brazil Current retroreflection or the Indian Ocean Great Whirl—and reach the equator as coastal Kelvin waves. Further regional investigations are needed to test this hypothesis.

The same panel shows that in the extratropics, the largest interannual variability is mostly found in the same areas in $T - 1/4°$ and $I - 1/4°$—that is, along the ACC, in western boundary current extensions (GS, NAC, Kuroshio, East Australian Current), in the Brazil-Malvinas Confluence, and in the Agulhas region. In other words, many extratropical regions where strong interannual variability is found with full forcing (and in AVISO to a large extent) are also regions with strong intrinsic interannual variability. This geographical agreement poleward of 20° is quantified by comparing the gray and black lines in the bottom panel of Fig. 6.

Comparing the contours and shadings superimposed in Figs. 2 and 3 shows that the spatial correspondence between eddy-active areas and regions of large low-frequency variability is stronger in $I - 1/4°$ than in $T - 1/4°$ and in AVISO (this is confirmed by the black, gray, and green lines in Fig. 4). The enhancement of this observed correspondence as the low-frequency forcing is switched off and its confinement in eddy-active midlatitude currents ($\pm 40°$) demonstrates that the broad-band temporal variability found in SLA along midlatitude eddy-active currents is intrinsically generated in the ocean, rather than forced.

### 2) MAGNITUDE AND REGIMES

The ratio $R_{LF}^{I-1/4°}$ between intrinsic ($I - 1/4°$) and total ($T - 1/4°$) interannual SLA variances at $1/4°$ resolution, defined in section 2c(2), is shown in Fig. 8 (map in the top panel, zonal average as a gray line in the bottom panel). Large and small values of $R$ correspond to predominant contributions of intrinsic and forced interannual variances, respectively. The top panel in Fig. 8 shows that $R$ remains smaller than 20% in the $10°S–10°N$ band, except in the North Brazil Current retroreflection and the area surrounding the Indian Peninsula where ocean-only processes bring the proportion of intrinsic-to-total variability up to 40%–80%. The standard deviation of intrinsic interannual variability is about 2 cm in these areas (top panel in Fig. 3).

Largest values of $R$ are found within the main eddying currents (ACC, Kuroshio, GS, NAC). In zonal average (bottom panel in Fig. 8), 80%–90% of the interannual variability simulated in $T - 1/4°$ is generated intrinsically between the ACC's Sub-Antarctic Front and the sea ice limit (top panel), in particular in its Atlantic and Indian sectors; this feature is further discussed in section 4c below. Relatively large contributions of intrinsic variability
(i.e., 60%-80%) are also found along the path of the GS-NAC and in the North Atlantic subpolar gyre, in the Kuroshio and in the northern Alaskan gyre, and east of Australia. As mentioned above, strong mesoscale activity is also found in these latter regions (contours in the same panel), suggesting the existence of dynamical interactions between both ranges of time scales. Indeed, animations and Hovmöller diagrams from *I* − 1/4° reveal that both HF and LF SLA variabilities are eddy sized in these regions, further supporting the idea that mesoscale eddies fluctuate over a wide range of time scales near the main currents.

In their Fig. 10b, Taguchi et al. (2007) show in the Kuroshio Extension a ratio $R^* = \sigma_{LF}^I/\sigma_{LF}^T$ that relates to our ratio $R$ in a simple way: $R \sim 100\% \times (R^*)^2$. Here $R^*$ was computed from a different model simulation [OGCM for the Earth Simulator (OFES)] at finer resolution (1/10°) over a longer, 42 yr, period with a different method. The eastward penetration of our simulated Kuroshio is shorter than observed (Fig. 2) and Taguchi et al.'s (2007) but comparing $R^*$ and $R$ reveals two clear similarities in this region. First, the LF variability in both simulations is largely intrinsic ($R^*$ and $R$ are close to 1 and 100%, respectively) in a broad region surrounding the Kuroshio Extension. Second, $R^*$ and $R$ exceed 1 and 100%, respectively, just poleward of the Kuroshio's variability maximum with peaks locally exceeding 1.1 and 120%, respectively ($R > 120$ within Fig. 7's cyan contours; see below). In other words, these quite different model setups confirm that intrinsic processes largely contribute to the total LF variability in this eddy-active region.

The low-frequency intrinsic variability is not as strong in the quiescent subtropics of both hemispheres; however, it nevertheless exhibits zonally elongated maxima between about 20° and 35°, both in absolute value (Fig. 3) and compared to the total LF variability (Fig. 8). In these latitude bands, which roughly correspond to the equatorward flanks of subtropical gyres, animations and Hovmöller diagrams reveal a different picture than that described in eddy-active regions: the LF intrinsic SLA variability consists of large-scale interannual Rossby wave-like signals that propagate to the west at a few centimeters per second, that are emitted in eastern basins and intensify toward the west throughout the whole 314-yr integration. The scale, phase speed, period, westward intensification, meridional location, and confinement of these first-mode-like baroclinic features are consistent with those found by Hazeleger and Drijfhout (2000) in their idealized study of intrinsic subtropical variability. This qualitative agreement is noteworthy and may indicate similar generation processes; however, the study of these intrinsically generated midlatitude

Rossby waves would require a complete characterization of the 3D variability, comparisons with various theories, and perhaps a series of sensitivity experiments in idealized or realistic contexts. This lies beyond the scope of this global description and is left for future studies.

In summary, our results show that a large percentage of the (realistic) SLA interannual variability present in the fully forced case is generated intrinsically: this percentage is actually larger than 40% over half of the global ocean area and exceeds 80% over one-fifth of it.

b. Forced LF variability at 1/4°: Main features

Comparing the black and gray lines in the top panel of Fig. 6 shows that switching on the interannual forcing (i.e., the transition between *I* − 1/4° and *T* − 1/4°) yields the largest increase in zonally averaged interannual variability between about 20°S and 20°N. The intertropical interannual variability in the real ocean, largely due to ocean–atmosphere interactions that are not explicitly resolved in *T* − 1/4° [see section 4a(1)], is accurately forced by the DFS4 forcing, in particular at these low latitudes (see section 3b).

Figure 7 shows the forced $\sigma^F$ variability map at 1/4° (top panel) and its zonal average (black line, bottom panel). This map shows that the intertropical interannual variability is markedly influenced by the atmosphere in the Pacific and Indian Oceans, that is, much more than in the Atlantic. This is consistent with the low-latitude distribution of ocean–atmosphere interannual variability modes: ENSO in the Pacific and the Indian dipole (Saji et al. 1999), which are partly connected (Schott and McCreary 2001) and have a larger interannual SST (hence, presumably, SLA) signature than their Atlantic equivalent (Zebiak 1993).

Noting that ENSO and the Pacific decadal oscillation (PDO) have comparable imprints on the North Pacific wind stress curl interannual anomalies (Gershunov and Barnett 1998), and that interannual atmospheric anomalies are stronger in the Pacific than in the Atlantic, the direct influence of the atmosphere on the midlatitude interannual variability [via, e.g., Qiu's (2003) PDO-related mechanism] is expected to be stronger in the Pacific than in the Atlantic. The top panel in Fig. 7 shows that, indeed, the forced midlatitude interannual variability is stronger in the eastern North Pacific ($\sigma_{LF}^{-1/4^\circ} > 2.5 - 3$ cm) than in the Atlantic.

Forced interannual variability is found along eastern boundaries, for example, along the northeastern boundary of the Indian Ocean (2 cm). The meridional extent of this maximum suggests that it is forced by equatorial air–sea coupling and propagates northward. Part of this atmospherically driven near-coast interannual variability is also likely to excite westward-propagating Rossby
waves (e.g., along the large-scale local maxima originating from eastern boundaries in the southern Indian and Pacific Oceans) as those mentioned above at midlatitudes.

### c. Intrinsic and forced LF variability at 1/4°

#### 1) POSSIBLE RELATIONSHIPS

As mentioned above, the cyan contours in the top panel of Fig. 7 enclose regions (mostly in the Southern Ocean and around 25°N) where the LF intrinsic variability found in $I - 1/4°$ exceeds by 20% or more the LF total variability found in $T - 1/4°$. This does not question the validity of our simulations but questions the two simple assumptions made in section 2—that is, 1) the total variability is the sum of the forced and intrinsic components and 2) both components are uncorrelated; at least one of these initial assumptions is invalid there. It is possible that assumption 1 is wrong (e.g., in the highly variable Southern Ocean) because of strong nonlinear interactions between both components. However, it is formally possible that assumption 1 is valid everywhere in $T - 1/4°$ if temporal correlations between the intrinsic and forced components are negative in contoured regions [see Eq. (1)]. In other words, both components might add up everywhere in the fully forced simulation but be partly out of phase in certain regions, hence yielding a total variance that is smaller than the intrinsic variance.

These two simple assumptions may not be formally valid in our complex simulations; however, they proved useful in estimating the forced variability. Further understanding the actual relationships between the forced and intrinsic LF variabilities would certainly require dedicated investigations, preferably in idealized contexts first.

#### 2) CASE OF THE SOUTHERN OCEAN

Figures 2 and 3 have shown that the observed distribution of the interannual SLA variability in the Southern Ocean is correctly simulated by our $1/4°$ model, with and without interannual forcing (in $T - 1/4°$ and $I - 1/4°$, respectively). Strong interannual variabilities in AVISO and in the latter two simulations are found in eddy-active regions (see Fig. 4), that is, along ACC fronts. The top panel in Fig. 9 shows the interannual SLA variability simulated with full forcing ($\sigma_{LF}^{T-1/4°}$) and black dots show the location of its maximum at each longitude. Away from the Kerguelen Plateau (75°–90°E), this maximum follows the Subtropical Front between about 30° and 65°E and the Sub-Antarctic Front between about 65° and 140°E [based on the map by Orsi et al. (1995)]. Figure 9's bottom two panels show the contributions of

the intrinsic and forced interannual variabilities at given meridional distances (ordinates) from the fully forced variability maximum shown as black dots in the top panel.

The right bottom panel in Fig. 9 shows that the intrinsic variability maximum ($I - 1/4°$, red peak) is slightly displaced southward (50 km) compared to the total variability maximum ($T - 1/4°$, black peak). Although set at first order by the topography (Sallée et al. 2008), the position of the ACC fronts may vary on interannual time scales, either because of the forcing (Fyfe and Saenko 2006) or just intrinsically. As the front constitutes a large gradient in SSH, forced and intrinsic shifts in its position should induce large SLA variabilities there. The forced and intrinsic variabilities, indeed, reach their (comparable) maximum amplitudes close to the main front.

Interestingly, most of the interannual SLA variability is forced by the atmosphere north of the total maximum (therefore including the area where the Subantarctic Mode Water is formed and subducted). On the contrary, the interannual SLA variability is almost fully intrinsic (red line) south of the total maximum (black peak); this is consistent with the aforementioned small temporal correlations and with their weak sensitivity to the presence of interannual forcing (middle panel in Fig. 6). The left bottom panel in Fig. 9 as well as Figs. 3 and 7 exhibit this meridional contrast over most of this 130° range of longitude.

The origin of this marked contrast is not clear, but we offer a hypothesis. The maximum in LF variability shown as black dots in the top panel of Fig. 9 approximately follows the southern edge of the so-called Southern Hemisphere supergyre [see Fig. 3 in Ridgway and Dunn (2007)]. This wide circulation pattern is predicted by Sverdrup dynamics (de Ruijter 1982) and is hence presumably sensitive to low-frequency fluctuations of the large-scale wind stress curl. Our results show that, indeed, the interannual SLA variability north of the black dots has large spatial scales (larger than 6°, not shown) and is mostly forced by the atmosphere. In the ACC south of this subtropical cell, the interannual SLA variability is not likely to be governed by the Sverdrup balance (Rintoul et al. 2001; LaCasce and Isachsen 2010), and hence it is presumably not to be structured at the basin scale. Indeed, we find that the interannual SLA variability south of the maximum is largely independent of the forcing and has much smaller scales than north of it.

To summarize, our $1/4°$ simulations exhibit sharp changes in the origin and spatial scales of the SLA interannual variability where subtropical Sverdrup-based dynamics change into circumpolar dynamics. More investigation and additional sensitivity experiments should help verify whether these concomitant transitions are linked physically.
<figure><img src="image_56.png" /></figure>

FIG. 9. The LF SLA variability in the ACC at 1/4° resolution. (top) Total variability ($\sigma_{\text{LF}}^{T-1/4°}$, cm). The sinuous black line follows the maximum of this field in $T - 1/4°$ at each longitude (excluding the Kerguelen Plateau) and corresponds to ordinate 0 in the other two panels. (bottom left) Percentage of variance explained by the intrinsic component ($R_{\text{LF}}^{\text{I, color}}$, color), as a function of longitude (abscissas) and of the meridional distance to the maximum of $\sigma_{\text{LF}}^{T-1/4°}$ (ordinate, km); the Kerguelen region is omitted; contours as in Fig. 7. (bottom right) Total (black), intrinsic (red), and forced (blue) LF standard deviations are averaged between 30° and 163°E (omitting the Kerguelen Plateau) at constant distances from the maximum of $\sigma_{\text{LF}}^{T-1/4°}$, and shown as a function of this distance.

**5. Role of resolved versus parameterized eddies**

We now isolate from the 2° simulations ($T - 2°$ and $I - 2°$) the intrinsic and forced components of the interannual SLA variability, and compare them with their 1/4° equivalents. Our main purpose here is to evaluate the possible sensitivity of these components to the explicit (although partial) resolution of mesoscale eddies at 1/4° versus their parameterization at 2°. The total interannual variability in $T - 2°$ is evaluated against altimetry (using $T - 2°_{P10}$ as a reference for conciseness) in section 5a below. Its intrinsic and forced components are described and compared to their $T - 1/4°$ counterparts afterward.

a. *Assessment of SLA variability in T - 2°*

The red solid lines lie on or above red dashed lines in all panels of Fig. 5. This means that at 2° resolution as well, the forcing update has a beneficial (or neutral) impact on the magnitude and distribution of SLA variability at all time scales. Therefore, the $T - 2°$ interannually forced simulation has globally more skill than

its earlier “P10” version in reproducing the observed SLA variability at all time scales and differs from $T - 1/4°$ in the way described in P10 regardless of forcing update. We now focus on the comparison of $T - 2°$ and $I - 2°$ to isolate the intrinsic and forced variabilities simulated at 2°.

b. *Intrinsic and forced low-frequency variability at 2°*

1) *INTRINSIC LF VARIABILITY AT 2°*

The maps in Fig. 3, the top panel of Fig. 6 (gray and pink lines), and Fig. 8 show that the transition from eddying to laminar regimes almost suppresses the intrinsic variability, bringing its contribution to the total interannual variability to very small percentages (smaller than 3% over half of the global ocean). This confirms that (even partly) resolved nonlinear eddies (and/or reduced dissipation) play a critical role in generating (and/or not dissipating) intrinsic LF variability. Although much smaller, the contribution $R_{\text{LF}}^{\text{I}}$ of intrinsic variability to the total reaches 20%–30% at 2° in a few regions where it is large at 1/4°, suggesting that noneddying
processes (e.g., large-scale Rossby or Kelvin waves) are involved there. However, in areas of strong eddy activity (ACC, western boundary currents), the decrease in the intrinsic contribution due to the coarser resolution is greater, further supporting the large involvement of nonlinear eddy processes in producing intrinsic variability.

## 2) FORCED LF VARIABILITY AT 2°

The shading in the middle panel of Fig. 7 shows the forced LF variability at 2°: the large-scale patterns of this component appear much less sensitive than the intrinsic component to the decrease in model resolution. Differences between $\sigma^{F-1/4°}$ and $\sigma^{F-2°}$ are found in the Southern Ocean; nonetheless, the 2° counterpart of Fig. 9's bottom right panel (not shown) indicates that the forced variability remains stronger north of the ACC front than south of it at 2° as well. This weak sensitivity of the forced variability to resolution strongly suggests that at first order in our forced simulations, and unlike the intrinsic component, the forced LF variability involves relatively linear dynamics; these include, for example, Rossby waves, Kelvin waves, Sverdrup-type responses, or wind-forced topographically trapped Taylor columns that are likely responsible for the southeastern Pacific maximum (Webb and de Cuevas 2003) visible at both resolutions.

Besides this general agreement, the zonally averaged forced LF variability (bottom panel in Fig. 7) is up to 20% larger at 1/4° than at 2° at mid- and high latitudes, this is, where eddies are most energetic. Two arguments may be proposed to explain this difference. First, horizontal viscosity is much reduced at 1/4° and may have less impact on the dissipation of horizontal motions, including forced motions. Second, the main currents—such as the ACC (see above), the Gulf Stream, or the Kuroshio—have a stronger signature in SSH gradients at 1/4°; interannual migrations of these fronts due to the forcing and will therefore yield a larger forced SLA variability at 1/4°.

# 6. Conclusions and discussion

The main purpose of this study was to estimate and describe the sea level expression of intrinsic and forced variabilities in the ocean at interannual time scales in a global, realistic context (using a 1/4° eddying OGCM with realistic geometry, initial condition, forcing, etc.). We have compared the interannual ($T > 18$ months) sea level anomalies (SLAs) simulated with two different forcing functions: a realistic forcing with temporal scales ranging from 12 h to decades (simulation $T - 1/4°$) and

the climatological seasonal cycle derived from the former (simulation $I - 1/4°$).

A quantitative comparison of the $T - 1/4°$ simulation with AVISO observations demonstrated the realism of the total sea surface variability (global distribution, magnitude) between 1993 and 2004, in particular at interannual time scales. The intrinsic interannual variability was diagnosed in $I - 1/4°$, and the forced component was deduced by subtracting the intrinsic from the total variance fields. Our main results concerning the intrinsic and forced components are summarized below.

* Without any interannual forcing, $I - 1/4°$ generates an intrinsic SLA variance that exceeds 40% of its fully forced counterpart over half of the global ocean. This intrinsic contribution is largest in strong eddy-active regions: its zonal average reaches 70%-75% south of the ACC and 50%-55% over the main western boundary currents and their extensions. Since their phase is not constrained by the forcing, correlation coefficients between intrinsic SLA fluctuations and AVISO time series at low frequency are small everywhere. The intrinsic component is weaker in absolute value but remains large in proportion within the 20°-35° latitude band in both hemispheres; there, preliminary diagnostics reveal the presence of interannual Rossby waves continuously emitted despite the absence of direct interannual forcing.

* The LF variability is largely influenced by the atmosphere (i.e., forced by it in the absence of air-sea coupling) in the tropical Pacific and Indian Oceans. The weaker interannual variability forced at low latitudes in the Atlantic (and perhaps at higher latitudes too) is likely due to the imprint of ENSO, which is largest in the Pacific and Indian Oceans. The strongest direct influence of the atmosphere on interannual SLAs far from the equator is found in the northeastern Pacific.

* As predicted by idealized studies and as seen in altimeter observations, the interannual variability is strong where mesoscale eddies are strong in our fully forced simulation. In other words, unstable currents exhibit a strong broadband temporal variability. This spatial correspondence is even clearer without direct interannual forcing ($I - 1/4°$), strongly suggesting the oceanic (eddy driven) origin of this interannual variability. The interannual SLA variability is largely intrinsic and has a broad range of spatial scales south of the ACC in the Indian Ocean Basin; in contrast, this variability is mostly forced by the atmosphere north of it and has typical scales larger than 6°. Explaining this sharp meridional contrast is left for the future.

* The intrinsic interannual variability is almost absent from our 2° simulations; however, the forced interannual