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

. 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 ◦ , 12 ◦ , and 14 ◦ resolutions. Simulated mean sea-surface heights (MSSH) and sea-level anomalies (SLA) are collocated over the period 1993–2004 onto the AVISO dataset. MSSH ﬁelds are compared with an inverse estimate. SLA datasets are ﬁl-tered and compared over various time and space scales with AVISO regarding three characteristics: SLA standard deviations

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 bandpassed 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 dependance on frequency and latitude.Benefits of enhanced resolution are greatest in the 1 • -1 2 • and 1 2 • -1 4 • transitions, in the 14-180 day range, and within eddy-active mid-and highlatitude regions.In the real ocean, most eddy-active areas are

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: subgridscale 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 1 • 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 Published by Copernicus Publications on behalf of the European Geosciences Union.
T. Penduff et al.: Sensitivity of sea-surface height features to resolution 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 1 2 • -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 Budich, 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/seaice simulations representative of laminar and eddy-admitting classes (2 • and 1 • , and 1 2 • and 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-2004hindcasts (DRAKKAR Group, 2007) were performed by the DRAKKAR consortium 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-1 http://www.ifremer.fr/lpo/drakkar/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 (modelobservation 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.

Model configurations
Our model setups are based on the NEMO code (Madec, 2008) and differ by their horizontal resolutions ( 1 4 1 shows that in the 2 • and 1 • configurations, the meridional resolution dy is increased toward the equator where it reaches 1 2 • and 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. (2006Brodeau et al. ( , 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   Griffies et al. (2009).In all runs, the bottom topography is discretized as partial steps for an accurate represention of topographic slopes and f 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 topographicallyrectified 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 2 s −1 in the 1 4 • and 1 2 • models respectively, and de-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 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 11 and 12 × 10 11 m 4 s −1 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 freeslip 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 47year (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)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004) when altimeter observations are available; this ensures a relatively long (35-year) and identical period of spinup in the four integrations.
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 timedependent 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 Ocean Sci., 6, 269-284, 2010 www.ocean-sci.net/6/269/2010/"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 observed 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 oscillations, in one and two dimensions.Figure 3 summarizes the space-time scales considered in the following.

Model-observation comparison statistics
Collocated SSH timeseries at location (i,j ) are noted η A (t) for AVISO and η m (t) for the m th model4 .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: (1) The temporal standard deviations of AVISO SLAs are denoted as σ A (i,j ); the ratio σ m /σ A will be called resolved 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: 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 approach requires the computation of local autocorrelation functions at lag τ of all η m (t) and η A (t) timeseries; these autocorrelation functions are noted ρ m (i,j,τ ) and ρ A (i,j,τ ), respectively.Individual temporal correlations C m t (i,j ) are then considered significant at the 95% confidence level where 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 significant temporal correlations, i.e. satisfying Eq. 3, are kept in the following.Resolution-induced changes in significant temporal correlations are shown in Fig. 5 for each frequency band (columns).
Mean and variable surface topography features are compared to their observed counterparts within 28 periodic latitude bands labeled λ ∈ (1;28) spanning the (70 • S-70 • N) range by 5 • intervals.
-Temporal standard deviation σ m (i,j ) maps are split into each latitude band and rearranged as a long onedimensional vector of length l containing elements noted σ m (λ,l).This vector is then averaged along dimension l to yield σ m (λ) = σ m (λ,l) l , i.e. the average of σ m in each latitude band.These quantities and their observed counterparts (noted σ A (λ)) are shown in various colors in the first row of Fig. 6 for each frequency band (columns).
www.ocean-sci.net/6/269/2010/Ocean Sci., 6, 269-284, 2010 -The same procedure is applied to (significant) temporal correlation maps C m t (i,j ).This provides C m t (λ), 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 m t (λ) 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 (σ A (i,j ) and σ 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 σ A and σ m maps in this latitude band, is then computed as follows: Here, φ l denotes the spatial average of any variable φ over a latitude band.The spatial standard deviation of σ m within each latitude band is computed as ; its AVISO counterpart is noted α A (λ).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 σ A fields in distinct frequency ranges will be noted • 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 fieds that were collocated at the same spatio-temporal resolution (

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.Spatiallysmoothed 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. 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.
(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 ° to 1 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.

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.

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 globallyaveraged observed magnitude are simulated at 2 • , 1 • , 1 2 • and 1 4 • 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.

Distribution of variability
Spatial correlations C m s 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 • to 1 4 • transition (green) improves the distribution of high-frequency SLA variability as much as the 1 2 • -1 4 • (red) transition in the Southern Ocean and the northern tropical oceans.

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 Ocean Sci., 6, 269-284, 2010 www.ocean-sci.net/6/269/2010/oceans.Significant temporal correlations C m t 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 1 2 • -1 4 • 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 1 4 • .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 • , 1 2 • , and 1 4 • resolutions seem to have comparable skill in resolving these relatively fast, forced nearcoastal processes.Whether model resolutions finer than 1 4 • would further increase (through a more accurate representation of forced signals) or decrease (by admitting more chaotic mesoscale turbulence) near-coastal C m t values remains to be determined.

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 quasiannual and mesoscale bands is large both in the AVISO dataset (0.76), and in the 1 4 • model simulation (0.72).Accordingly, the RMS difference between global, observed σ A maps of quasi-annual and mesoscale variability is relatively small (2.3 cm), and is the same in the 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 1 4 • datasets.

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 (σ A ∼6cm) and 10 • S (σ 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 1 4 • simulation, but falls to around 50% at 1 2 • resolution, and below 40% in both laminar runs.In this quasi-annual band, the 1 • 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).

Distribution of variability
Figure 6 (second column, third row) shows that the agreement between the observed and 1 4 • maps of quasi-annual variability is greatest between 20 • S and 20 • N (C 1/4 • s ∼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) 1 2 s , and (iii) the largest improvements in quasi-annual SLA variability maps are mostly obtained in the 1 • resolution transitions poleward of 30 • (second column, third row in Fig. 7).

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 everytime 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 • -1 2 • 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 m t values increase with resolution between 2 • and 1 2 • , but drop by about 0.06 when resolution reaches 1 4 • , i.e. back to their 2 • values.However, these changes occur over a cross-shore scale (about 1 3 • ) that seems small compared to those (about 1 • ) of coastallytrapped processes that may prevail at quasi-annual timescales (e.g.upwellings).The interpretation of these differences is thus unclear and would require dedicated investigations.The global spatial correlation coefficient (0.72) between observed variability maps σ A in the interannual and quasiannual 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 σ distributions in the 3 frequency bands is also clear in the 1 4 • 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 eddyactive regions (KS, GS-NAC, Agulhas Retroflection, ACC, Brazil-Malvinas Confluence), both in the real ocean and in the 1

4
• 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.

Variability levels
Interannual timescales are those at which the 1 4 • 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 1 2 • and 60% in both laminar models.A similar decrease with resolution is seen at quasi-annual timescales (76, 68, 62, 61%, respectively, from 1 4 • to 2 • ).Resolved interannual variabilities undergo comparable increases in the 1 • -1 2 • and 1 2 • -1 4 • resolution transitions: about +10% north of 30 • N and more than 20% in the Southern Ocean.

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 m s and C m t 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 • -1 2 • (and also 1 2 • ) 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).

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 1 4 • 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 1 4 • SLA datasets.In other words, a substantial part of the observed and 1

4
• 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 1 4 • 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 1 4 • 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 To summarize, extending the spectrum of resolved scales toward the first internal deformation radius progressively improves the representation of interannual variability over a broad range of spatial scales, including those much larger than the grid size.

Conclusions
The main aim of this study was to assess and better understand the effects of resolution on global ocean model solutions, as considered through sea surface height.The outputs from four 47-year global ocean/sea-ice simulations, performed at 2 • , 1 • , 1 2 • , and 1 4 • resolutions, and driven by identical surface forcing, have been collocated onto AVISO sea level anomaly (SLA) observations.This yielded a consistant 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), geographical distribution of sea-level variability (using spatial correlations between observed and simulated maps of SLA variability), and phase agreement (significant temporal correlations) 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 benchmark regarding the comparative skills of two climateoriented laminar models (2 • and 1 • , i.e. resolutions typical 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 • 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 levels in global average).The transition between 1 2 • and 1 4 • resolution yields the strongest enhancement in eddyactive regions, as expected.Resolved SLA variability 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 observed 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 interannual 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 interannual SLA timeseries from their AVISO counterparts.Since the atmospheric forcing is unchanged, these sensitivities probably come from the more non-linear dynamics simulated at 1 4 • .To clearly identify the origin 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 hypothesize that both the strong enhancement and (realistic) confinement in eddy-active areas of SLA variabilities at all timescales, and the aforementioned temporal • 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 climateoriented 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.

Fig. 2 .
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.

Fig. 3 .
Fig. 3. Scheme of the 4 bands in which model-observation comparisons are performed.The Large-scale Interannual band is shown in light gray.
Fig. 4. SLA standard deviation maps (cm) in AVISO observations (σ A , top) and the 4 model simulations (σ 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.

Fig. 6 .
Fig.6.Observed and simulated SLA standard deviations (σ (λ) in cm, upper panels), significant temporal correlations and density of gridpoints with significant temporal correlations (C m t (λ) in plain lines and p m t (λ) in dashed lines, middle panels), and spatial correlations between observed and simulated maps of SLA (C m s (λ), lower panels).Results are shown in average over latitude bands.See text for details.

Fig. 7 .Fig. 8 .
Fig. 7. Resolution-induced changes in resolved SLA variabilities ( σ m2 σ A − σ m1 σ A , first row), SLA temporal correlations (C m2 t − C m1 t , second row), and SLA spatial correlations (C m2 s − C m1 s , third row) as a function of latitude (abcissae).Results are shown for four resolution transitions m1 → m2 (colors, see legend).Lines in the second row are dashed at latitudes where p m1 t or p m2 t (densities of points with significant temporal correlations) is smaller than 20%.