## Abstract

This paper exploits a new observational atlas for the near-global ocean for the best-observed 3-yr period from December 2003 through November 2006. The atlas consists of mapped observations and derived quantities. Together they form a full representation of the ocean state and its seasonal cycle. The mapped observations are primarily altimeter data, satellite SST, and Argo profiles. GCM interpolation is used to synthesize these datasets, and the resulting atlas is a fairly close fit to each one of them. For observed quantities especially, the atlas is a practical means to evaluate free-running GCM simulations and to put field experiments into a broader context. The atlas-derived quantities include the middepth dynamic topography, as well as ocean fluxes of heat and salt–freshwater. The atlas is publicly available online (www.ecco-group.org). This paper provides insight into two oceanographic problems that are the subject of vigorous ongoing research. First, regarding ocean circulation estimates, it can be inferred that the RMS uncertainty in modern surface dynamic topography (SDT) estimates is only on the order of 3.5 cm at scales beyond 300 km. In that context, it is found that assumptions of “reference-level” dynamic topography may yield significant errors (of order 2.2 cm or more) in SDT estimates using in situ data. Second, in the perspective of mode water investigations, it is estimated that ocean fluxes (advection plus mixing) largely contribute to the seasonal fluctuation in heat content and freshwater/salt content. Hence, representing the seasonal cycle as a simple interplay of air–sea flux and ocean storage would not yield a meaningful approximation. For the salt–freshwater seasonal cycle especially, contributions from ocean fluxes usually exceed direct air–sea flux contributions.

## 1. Introduction

Estimates of the oceanic state, whether instantaneous or averages over varying time intervals, have a variety of uses. In particular, “climatologies” are meant to represent the long-term mean ocean state. Such climatologies (e.g., Stephens et al. 2002; Gouretski and Koltermann 2004; Curry 2001) are widely used to initialize and test model solutions and to describe flows and transports. Almost without exception, such climatologies are based upon some form of objective mapping involving weighted averages of the data in both time and space. Despite their evident utility, a limitation of these climatologies is that the mapped temperatures, salinities, and/or derived flows do not follow from known equations of motion or kinematics. An early attempt to produce time-averaged fields satisfying known equations was given by Wunsch (1994), but the methodology was not practical on a global scale.

Another difficulty with climatologies is the extreme inhomogeneity in both space and time of the underlying oceanic sampling [e.g., Forget and Wunsch 2007; Wunsch et al. (2007) for a discussion of sampling issues]. The recent increase in in situ data coverage due to the advent of the Argo array is such that, over large fractions of the ocean, more samples have been collected over the last few years than over the previous 100 yr (e.g., Forget and Wunsch 2007). It seems obvious that this dense database yields relatively robust estimates of the contemporary ocean state, but it is unclear whether it may yield significant improvements in representing the long-term mean ocean state (i.e., the ocean climatology). Unless the temporal bias in data coverage is handled properly, new climatology estimates may indeed only be remotely representative of the long-term mean ocean state and boil down to estimates for the few recent years instead. This issue is clearly illustrated by Fig. 1 showing the mean date of observation for two hypothetical climatologies—one computed before the advent of Argo (top panels) and one computed afterward (bottom panels). The complexity of Fig. 1’s maps is challenging, even when the recent years’ data are not included (top panels). It suggests a need for caution in interpreting in situ climatology estimates as representations of the long-term mean ocean state in general. In that respect, the likely temporal bias toward present-time climatology estimates that include the vast recent Argo data is of most concern (bottom panels).

Here, the intention is to provide a mean monthly atlas over the comparatively intensely sampled recent period 2004–06. The result is a simple time average of solutions to known numerical equations of motion based upon the Massachusetts Institute of Technology (MIT) GCM (Marshall et al. 1997; Adcroft et al. 2004). Observational maps are created by solving for a nonlinear least squares fit of GCM trajectories to a variety of observations (listed in Table 1). Spatial interpolation is provided by the GCM itself, and the underlying time evolutions are consistent with the model equations. Availability of Argo data, in particular, motivates the choice of time interval. The greatly increased data coverage means that this interval should be a suitable reference time against which to measure *future* shifts in the oceanic state.

Before presenting the results, we first summarize the observations employed (section 2) and briefly review the estimation approach (section 3). Further methodological details are provided in appendixes at the end of the paper. Section 4 establishes the reliability of the atlas for observed quantities. Sections 5 and 6 analyze oceanographic implications of the results. Section 7 concludes the paper.

## 2. Datasets

The near-global ocean datasets included as constraints in the GCM interpolation problem are (i) Argo subsurface profiles of temperature (hereafter *T*) and salinity (hereafter *S*); (ii) monthly mapped sea surface temperature (SST) satellite measurements from infrared sensors (Reynolds and Smith 1994) and microwave sensors (maps produced by Remote Sensing Systems); (iii) along-track sea level anomaly measurements by satellite altimeters (see Table 1); and (iv) mean surface dynamic topography based on the CLS01 mean sea surface and the EIGEN-GRACE-03S geoid model (but no in situ data) provided by M. H. Rio (Rio et al. 2005). Regional in situ datasets are also included as observational constraints (see list in Table 1). The following simplifications are used: 1) geoid and altimeter data are combined to form a single (along track) absolute sea level constraint; 2) mapped SST data are used as constraints rather than raw pixel SST data; and 3) no sea-ice data is used as constraint.

An adjusted version of the *World Ocean Atlas* 2005 (*WOA05*; Locarnini et al. 2006) climatology for *T* and *S* is also used as part of the GCM interpolation problem, to constrain regions where Argo does not provide observations (typically below 2000-m depth and over high-latitude ice-covered regions). To avoid encroaching on the fit to Argo observations elsewhere, *WOA05* maps are adjusted to Argo observations using a univariate interpolation method (see appendix A).^{1} The resulting *WOA05*–Argo blend is then simply applied as a background constraint everywhere in the GCM interpolation problem. As Argo data are introduced twice (directly as profiles and via the *WOA05*–Argo blend) in the GCM interpolation problem, other observational constraints are upweighted by a factor of 2 to maintain a balance between data types. This least squares upweighting is omitted in section 4 for simplicity.

The in situ data coverage of the upper ocean (above 2000-m depth) for the period from 2004 to 2006 is illustrated in Fig. 2. Most of the global upper ocean has been observed during that period, which is largely due to Argo floats. Hence an atlas for that period may improve upon (long-term mean) climatologies over most of the upper ocean. In Fig. 2, the median number of observations per grid point is 3, while 15% of the grid points show more than 10 observations, and 30% of the grid points show no observations. Of concern are some relatively wide regions that lack observations altogether, because they present a particular challenge to ocean observation by drifting floats [such as sea ice coverage (e.g., in the Southern Ocean) or topographic obstacles (e.g., in the Indonesian seas, or near coasts in general)]. Unlike for the majority of the global upper ocean, a 2004–06 atlas may not improve much upon climatologies in those regions that lack contemporary in situ observations (Fig. 2). The same is true for the entire deep ocean (below 2000-m depth) that Argo floats do not reach.

## 3. GCM interpolation problem

The constrained least squares problem (formulated in appendix A) is similar to that solved by Wunsch and Heimbach (2007). The reader is referred to Wunsch and Heimbach (2007) for a discussion of the method of Lagrange multipliers used to solve the constrained least squares problem. As for the latest Estimating the Circulation and Climate of the Ocean (ECCO) Global Ocean Data Assimilation Experiment (GODAE) estimates (Wunsch and Heimbach 2008), a sea ice model and a surface atmospheric layer model are included. The GCM configuration is summarized in appendix B. The single most important difference compared to Wunsch and Heimbach (2007, 2008) is a reduction of the estimation interval(s), from 14–16 yr to 16 months here. Each calculation estimates the oceanic state over one calendar year plus the previous 4 months, hence covering slightly more than one seasonal cycle. The results of three such (overlapping) calculations are time averaged (see appendix A for details) first to form a 3-yr daily time series, then to form a mean monthly atlas representing an average seasonal cycle for the period from 2004 to 2006. The resulting atlas is hereafter called the Ocean Comprehensible Atlas (OCCA). The term “comprehensible” denotes that the underlying equations of motion are known. Unsurprisingly, the shorter estimation interval(s) produce a closer fit to observations than does the 14–16-yr ones used by Wunsch and Heimbach (2007, 2008) (Fig. 3), as initial conditions are more efficient controls over the shorter period. The reduction of the estimation interval(s) to 16 months can be thought of as a practical means to confine model error accumulation [as explained by Forget et al. (2008b); see also appendix A]. OCCA consists of a set of fields for each month of the year. It does not resolve interannual fluctuations. The underlying 3-yr daily time series may be used to study those, but this is beyond the scope of this study.

## 4. The fit to observations

The reliability of OCCA maps for observed quantities must first be established. OCCA is meant to synthesize mainly three datasets: Argo profiles, SST data, and altimeter data. For each one, the accuracy of alternative univariate atlases (due to other research groups) will be taken as a reference. Formal error estimates are generally lacking in ocean state estimates, whatever their origin. The present study is no exception. Although this major issue must eventually be resolved, it lies beyond the scope of the present study. Yet in practice, normalized estimate-to-observation distances can be used to crudely compare the accuracy of various large-scale estimates. This framework shall first be laid out (next paragraph), then applied to OCCA and alternative univariate atlases (Figs. 3 –8).

Assume that an ensemble of *N* observations, {*x _{i}*,

*i*= 1,

*N*}, of a quantity

*X*are available with normally distributed noise. Let

*X̃*be an estimate of

*X*, the true value of

_{T}*X*. Let

*σ̃*be an estimate of

*σ*, the true noise standard deviation. We can define a normalized estimate-to-observation distance as

_{T}It is clear that the expected values for and are 0 and *σ _{T}*

^{2}, respectively. It follows that the expected value of

*J*(

*X̃*) is

*σ*

_{T}^{2}/

*σ̃*

^{2}+ (

*X̃*−

*X*)

_{T}^{2}/

*σ̃*

^{2}, where

*X̃*−

*X*is the actual estimate error.

_{T}^{2}This relation implies that (i)

*J*(

*X̃*) increases with the actual estimate error amplitude; (ii) for a given set of estimates, {

*X̃*,

_{j}*j*= 1,

*K*}, the range of

*J*(

*X̃*) decreases as the estimated noise level (

_{j}*σ̃*) increases; and (iii)

*J*(

*X̃*) does not simply depend on the number of observations. In using pointwise measurements to derive large-scale state estimates, representation error includes eddy signals. The amplitude of eddy signals shows strong spatial variations so that

*σ̃*must vary in space (e.g., Forget and Wunsch 2007). Domain averages of

*J*(

*X̃*) would otherwise only represent the most energetic regions. Property 1 makes the comparison of

*J*(

*X̃*) values a practical method to compare the accuracy of large-scale estimates, even though it should not be mistaken for an analysis of formal errors (that would not behave according to properties 2 and 3). We note that the following analysis could be extended by examining distributions of (

*X̃*−

*x*)/

_{i}*σ̃*in detail, which are largely Gaussian (not shown).

For the subsurface hydrography, OCCA is evaluated relative to the *WOA01* climatology (Stephens et al. 2002), taken as a best-reference estimate of the long-term mean ocean state from pre-Argo observations. We postpone comparisons to the latest climatology estimates (such as *WOA05*) that are a blend of recent Argo data and historical in situ data. It is indeed unclear how these climatology estimates should be interpreted, since they may show a large temporal bias toward present time (see Fig. 1).

Figure 4 shows that OCCA is more consistent with contemporary Argo profiles than is the *WOA01* long-term mean ocean state estimate (Stephens et al. 2002). This result is encouraging because it suggests that OCCA reflects some of the ocean state characteristics specific to the period 2004–06. Unsurprisingly, the biggest differences in *J*(*X̃*) occur in the Southern Ocean where the climatology is based upon a very thin database (Forget and Wunsch 2007). The full model time series (from which OCCA is a time average) is marginally closer to the Argo profiles than is OCCA, which suggests that some of the model variability that is averaged out upon moving from the time series to the atlas may be realistic. A somewhat opposite conclusion arises when Fourier truncation is applied to smooth the climatology seasonal cycle (see appendix A for details), which marginally reduces misfits to Argo profiles. This behavior is thought to reflect seasonal cycle noise in the original monthly climatology (see section 6).

Misfits between observations and truly optimal large-scale estimates ideally should consist of randomly distributed noise (e.g., with no large-scale or long-term patches). Relative to a long-term mean climatology such as *WOA01*, misfits between OCCA and observations do usually consist of small-scale noise to a good extent. This is illustrated for salinity at an intermediate depth in the Southern Ocean by Fig. 5. At first glance, the differences in hydrographic structures between the two atlases are relatively subtle (Fig. 6). One noticeable feature is that the salinity minimum along the Antarctic Circumpolar Current (ACC; between 40° and 65°S) shows almost the same value at all longitudes in OCCA, but not in the climatology. An important question is whether such differences can be attributed to climate fluctuations, which is left for further investigation.

For sea surface temperature (SST; Fig. 7) and sea level anomaly (SLA; Fig. 8), comparison of the distances to observations yields similar conclusions about the respective accuracies of OCCA and *WOA01* estimates. Note that Fig. 8 does not address the annual mean sea level that will be discussed in section 5. Also, for *WOA01*, the seasonal cycle in dynamic height is used as an approximation of the seasonal cycle in SLA. This approximation was tested for OCCA and makes almost no difference for *J*(*X̃*) (not shown). The range of *J*(*X̃*) is the narrowest for bin-averaged along-track SLA observations (Fig. 8) and the widest for mapped SST data (Fig. 7). This simply reflects (see property 2 above) that the signal-to-noise ratio (defined as for the *X̃ _{j}*,

*j*= 1,

*K*estimates) differs from one dataset to another.

Slightly more insight can be obtained by comparing *J*(*X̃*) for OCCA and for *satellite atlases*. First, for SST, a mean monthly satellite atlas is formed by time averaging the monthly data provided by Remote Sensing Systems and Reynolds and Smith (1994). Accuracy of OCCA is similar to that of this satellite atlas (Fig. 7), reflecting that OCCA captures most of the SST seasonal cycle. Second, for SLA, a mean monthly satellite atlas is formed from Archiving, Validation, and Interpretation of Satellite Oceanographic data (AVISO)/CLS maps. These maps synthesize altimeter observations down to the mesoscale (on a ¼° × ¼° mesh). OCCA, however, may only effectively resolve signals that span a few grid points on a 1° × 1° mesh. Before comparison with OCCA, AVISO/CLS maps are smoothed to filter out spatial scales smaller than 3°. The mean monthly satellite atlas is formed by time averaging these smoothed maps. Consistent with our large-scale objectives, the accuracy of OCCA is similar to that of this satellite atlas (Fig. 8).

In summary, OCCA includes accurate maps for observed quantities, synthesizing Argo, SST, and altimeter data. In resolving the large scales, the accuracy of OCCA maps is similar to that of univariate maps. The key additional value of the new atlas follows from the inclusion of dynamical principles, which relate the various fields to one another. OCCA includes fields not only for observed quantities but also for many other quantities (e.g., volume, heat, and salt–freshwater fluxes). As part of the GCM interpolation procedure, all quantities (observed or otherwise) are made consistent with the full suite of observations. OCCA fields for quantities that are not directly observed may be thought of as *derived estimates*. These fields allow for an extended interpretation of observations. The following sections provide examples of derived estimates and discuss oceanographic implications.

## 5. OCCA annual mean circulation

The classical problem of determining the “reference-level” circulation is embedded in GCM dynamics. Previous studies (Forget et al. 2008a,b) have provided evidence that sets of Argo profiles allow successful inversions for the barotropic and baroclinic circulations in a GCM interpolation framework. The focus in this section is on a few important aspects of the annual mean circulation, while the following section assesses seasonal fluctuations in the circulation.

### a. Surface dynamic topography

In the absence of sufficiently accurate geoid estimates (Vossepoel 2007), improving estimates of time-mean surface dynamic topography (SDT) remains an outstanding issue. By definition, SDT is the difference between sea surface height and geoid. Many SDT estimates consist *simply* of a difference between a time-mean sea surface height estimate and a geoid estimate (hereafter MSS-geoid estimates). Some recent SDT estimates (Stammer et al. 2002; LeGrand et al. 2003; Maximenko and Niiler 2004; Rio et al. 2005; among others), however, employ additional observational and/or dynamical constraints.

The OCCA SDT estimate is thus an adjusted version of a MSS-geoid estimate provided by M. H. Rio. The adjustment is obtained, as part of the GCM interpolation procedure, by fitting altimetric and in situ data, along with the MSS-geoid SDT. Rio et al. (2005) have also estimated an SDT adjustment by adding altimetric and in situ data constraints, but using a very different method. We shall here compare the two estimated adjustments. While computational details are reported in appendix C, two aspects are worth emphasizing here: (i) Rio et al. (2005) and this study both start from the same MSS-geoid, except for changes of reference time periods; (ii) scales smaller than 300 km are filtered out of both adjustments for comparison purposes.

The two estimated adjustments show RMS values of 3.5 cm (OCCA) and 3.3 cm (Rio et al. 2005). They both may be regarded as estimates of errors in the MSS-geoid SDT. Hence we attribute an RMS uncertainty between 3.3 and 3.5 cm to the MSS-geoid SDT at scales beyond 300 km. This inference is consistent with the RMS uncertainty range suggested by Vossepoel (2007), which is from ≈3.5 to ≈4.2 cm for MSS-geoid SDT at scales beyond 334 km (see their Fig. 4). The two adjustments show fairly contrasting patterns though (Fig. 9), so that their cross correlation is only 0.3, and this raises suspicion about the details of both adjustments. The RMS difference between the OCCA SDT and Rio et al. (2005) adjusted SDT is 4.0 cm, while that between the OCCA SDT and MSS-geoid SDT is 3.5 cm. We thus attribute an RMS uncertainty between 3.5 and 4.0 cm to the OCCA SDT, which is basically the state of the art.

### b. Middepth dynamic topography

Under the hydrostatic approximation, dynamic height (DH) for the layer between the surface and 1600-m depth (0–1600) equals the difference between surface dynamic topography and 1600-m dynamic topography (DT)(1600). These three quantities [i.e., SDT, DT(1600), and DH(0–1600)] are estimated jointly in OCCA according to geoid and altimeter constraints (of SDT), Argo constraints (of dynamic height), volume conservation, and other constraints. As part of the three-term hydrostatic balance, DT(1600) is the term that is less directly related to observations and may therefore be regarded as the *inverted quantity*. Note that inversion using GCM interpolation involves no ad hoc assumption of reference-level circulation.

The estimated DT(1600) reflects expected middepth circulation pathways (Fig. 10), including North Atlantic subpolar gyre recirculations (Lavender et al. 2005) and Southern Hemisphere subtropical gyres (Davis 2005). For the Antarctic Circumpolar Current, the estimated difference in DT(1600) across the current is 47 cm, which is a third of the surface value (150 cm).^{3} Other noteworthy features include widespread midlatitude contrasts in DT(1600) between oceanic basins: 5 cm between the Pacific and Atlantic, and 3 cm between the Pacific and Indian.^{4} Geostrophic streamlines tend to align along continental margins in the DT(1600) map, except in regions of significant vertical motions. In the subpolar North Atlantic, streamlines appear near margins due to (well known) high-latitude deep downwelling. In the subtropical western North Atlantic, streamlines vanish as the deep southward flow partly upwells toward intermediate depths.

Over the global ocean, DT(1600) shows a standard deviation of 24 cm, mostly due to ACC and isolated basin values. For the region between 40°S and 60°N and excluding the Mediterranean and the Japan Sea, DT(1600) shows a standard deviation of 2.2 cm. The magnitude of DT(1600) contrasts (Fig. 10) is not small compared with SDT uncertainties and estimated SDT adjustments (Fig. 9). Hence we conclude that erroneous assumptions of reference-level circulation, corresponding to the middepth dynamic topography, may seriously undermine efforts to improve SDT estimates using in situ data.

### c. Example of transport estimates

Provided that the large-scale flow is near geostrophic, observational constraints of the hydrography and surface dynamic topography can provide efficient constraints of transports. Lacking direct transport observations, the new atlas transports can hardly be tested though—only compared with independent estimates. The emphasis here is on volume transports for two of the best-documented sheared flows.

The estimated 3-yr mean ACC transport through Drake Passage is 152 Sv.^{5} The estimated shear consists of a 42-Sv cell reversing about 1400-m depth (Fig. 11). As in previous studies, the flow is now described using two layers separated by 2500-m depth. The flow in both layers (140 Sv above 2500 m and 12 Sv below) is slightly larger than the reference estimate (125 ± 10 Sv above 2500 m and 9 ± 5 Sv below) of Whitworth and Peterson (1985) and Whitworth (1983) by about 15%. The estimated ACC variability at Drake Passage (5.7 Sv) is mostly unsheared, while the reversing cell strength barely varies (0.9 Sv).^{6} An eddying Southern Ocean state estimate (so-called SOSE; M. Mazloff 2008, personal communication) shows a fairly similar mean flow structure (154 Sv top to bottom, and 55 Sv for the reversing cell) and variability structure (4.8 Sv top to bottom, and 1.6 Sv for the reversing cell).

At 27°N the estimated North Atlantic meridional overturning rate is 16 Sv.^{7} It includes a southward contribution of 7 Sv below 3000 m, which is large compared with unconstrained *z*-coordinate model solutions (see, e.g., Willebrand et al. 2001; Forget et al. 2008b). To facilitate comparisons with other observational estimates, the flow is now decomposed into an ocean interior flow (between the Bahamas and the Africa shelf) and a Florida Strait flow. The Florida Strait northward transport (26.5 ± 3 Sv) is 20% smaller than cable measurements indicate (31.5 ± 3 Sv). The estimated total ocean interior flow compensates the estimated Florida Strait transport by construction and includes a 3.7-Sv northward Ekman transport. The ocean interior shear consists of a 9-Sv cell reversing at 500-m depth, with a northward deep branch (Fig. 12, left). Figure 12 shows that the geostrophic ocean interior shear in OCCA (12 ± 2 Sv) is 20% smaller than the value (15 ± 3 Sv) computed from the Rapid Climate Change (RAPID) mooring data (Fig. 12, right). A difference in deep geostrophic shear is evident between the RAPID and OCCA estimates (Fig. 12, right). A test computation for OCCA (gray curves) shows that this difference is likely a consequence of the mooring subsampling. Additional RAPID array data (e.g., over the mid-Atlantic ridge) might help resolve this issue.

For volume transports (OCCA and others), the above analysis suggests 20% uncertainties. That figure is merely an educated guess, however, based upon only two sheared flows.

### d. Concluding remarks

The OCCA estimate of the annual mean ocean circulation is in reasonable agreement with alternative estimates of dynamic topography and transports. For transports especially, this result reflects that global observing systems provide efficient constraints of the density field (Argo, etc.) and of the surface dynamic topography (altimeters, etc.).

The atlas middepth dynamic topography, DT(1600), is proposed as a practical means to convert observed transects of *T* and *S* (e.g., along ship tracks) into absolute geostrophic transport transects. DT(1600) lies well within the Argo-sampled layer and the flow is already relatively weak at this depth. By construction, DT(1600) is a solution to the global ocean circulation inverse problem including grid-scale to global conservation principles. Direct application of DT(1600) as an estimate of geostrophic reference-level streamfunction is thus a convenient alternative to challenging synoptic inversion procedures. DT(1600) (section 5b) implies realistic surface dynamic topography (section 5a) and transports (section 5c) when it is combined with hydrography observations (section 4).

## 6. OCCA seasonal cycle

The real ocean seasonal cycle of heat content (freshwater content) has to balance fluxes of heat (freshwater), and this dynamical principle is built into OCCA. One of the interpretations of the GCM interpolation problem is an attempt to infer air–sea fluxes and ocean fluxes that are responsible for observed fluctuations (e.g., Stammer et al. 2004). Let us first provide an overall assessment of OCCA mean monthly surface fluxes of heat and freshwater. The RMS mean monthly differences between OCCA and National Centers for Environmental Prediction (NCEP) surface fluxes are shown in the bottom panels of Fig. 13. These maps are regarded as crude estimates of uncertainties that may equally apply to both flux estimates (OCCA and NCEP fields). While OCCA fluxes are adjusted to best fit ocean observations, these uncertainties (bottom panels) cover a relatively large fraction of the very fluxes (top panels). Keeping this in mind, let us then focus on the seasonal fluctuation (i.e., deviations from the annual mean) and extend the discussion to inferences of fluxes within the ocean.

### a. Implied fluxes

For any monthly hydrographic atlas, the varying flux divergence that is required to explain the seasonal cycle of heat and freshwater content can be computed as appropriately scaled differences between successive months. We shall refer to the result of this diagnostic computation as “implied fluxes,” and focus on their magnitude. This diagnostic computation does not address annual mean fluxes, but only deviations from the annual mean. We also note that implied fluxes may equally be attributed to fluxes through the ocean surface (air–sea fluxes, ice–sea fluxes, and runoff) or fluxes within the ocean (advection and mixing). In section 6b, we will distinguish between those two contributions.

The implied fluxes for *WOA01* and OCCA are strikingly different: *WOA01* usually shows much larger implied fluxes than does OCCA (Fig. 14). It is suggested hereafter that *WOA01* implied fluxes might be too large because of a lack of observational and dynamical constraints. It shall not be excluded that OCCA implied fluxes are too small, and further investigation is suggested on those grounds. Artificially large implied fluxes might reflect noise in the seasonal cycle, that is, month-to-month fluctuations inherited from eddy or instrumental data noise, especially if the database is thin. Such noise is partly precluded in OCCA because seasonal fluctuations are imposed by the GCM to occur gradually as a response to air–sea interactions. Temporal smoothness constraints that are statistical rather than dynamical can alternatively be added. Fourier series truncation to the first two annual harmonics is here chosen and applied to *WOA01* (see appendix A for computational details). This operation largely reduces implied fluxes (Fig. 14, top and middle panels). A significant degradation of the seasonal cycle estimate would imply an increase in estimate-to-observation misfits. However, Fourier truncation tends to reduce *WOA01*-to-observation misfits (see section 4). This behavior argues in favor of the smaller implied fluxes shown by OCCA. For *T* north of 30°S, there is good agreement in implied fluxes between Fourier truncated *WOA01* (middle) and OCCA (bottom), where they both rely on a relatively dense database. This result again argues in favor of the smaller implied fluxes shown by OCCA. For *T* south of 30°S and for *S* almost everywhere, Fourier truncated *WOA01* still shows much larger implied fluxes than does OCCA. *WOA01* relies on a thin pre-Argo in situ database there, while OCCA relies on good coverage by Argo floats, satellite SST, and altimeters. Except for *T* north of 30°S, Fig. 14 suggests that the pre-Argo in situ database is not sufficient to estimate accurately the seasonal cycle.

### b. Surface versus ocean fluxes

OCCA not only yields implied flux diagnostics, but also provides detailed estimates of surface fluxes (air–sea fluxes, ice–sea fluxes, and runoff) and ocean fluxes (advection and mixing) contributing to the seasonal cycle of heat and freshwater content. For simplicity, all advective fluxes (i.e., Eulerian and parameterized eddy effects) and all diffusive fluxes (i.e., interior and boundary layer effects) are hereafter presented jointly as ocean fluxes. Our focus is on splitting the seasonal cycle of heat and freshwater content into surface and ocean fluxes.

First, mixed layer budgets are considered, taking the Southern Ocean as an example. The study of Dong et al. (2007) provides a reference for heat. There is a satisfactory agreement between OCCA and Argo profiles for mixed layer depth diagnostics (Table 2). The study of Dong et al. (2007) and OCCA show remarkably similar mixed layer heat budgets (cf. their Fig. 6a; our Fig. 15 top), despite major differences in parameterizations of surface and ocean fluxes. Both results show that heat content tendency and surface fluxes best match each other in austral winter and spring, whereas ocean fluxes cause a sizable cooling in austral summer (opposing warming from the atmosphere). Furthermore, there is a good quantitative agreement between the two results. While the previous *consistency check* for heat is rather encouraging, OCCA readily allows an extension of the analysis to salt–freshwater (Fig. 15, right). Unlike for temperature, surface and ocean fluxes contribute about equally to the seasonal cycle in mixed layer salinity. Typically, the mixed layer salinity increases in austral spring because of ocean fluxes, and decreases in austral fall because of surface fluxes.

Second, top to bottom integrated budgets are considered (Figs. 16 and 17), following on Fig. 14. For heat content, surface fluxes (Fig. 16, top left) are only found to explain most of the seasonal cycle (Fig. 14, bottom left) in open-ocean subtropical and subpolar regions. Ocean fluxes (Fig. 16, bottom left) provide major contributions in the tropics and in the vicinity of intergyre jets. Strikingly for salt–freshwater, ocean fluxes (bottom right) are found to dominate over surface fluxes (top right) almost everywhere, except over high-latitude regions that experience sea ice coverage fluctuations. These results clearly stress the importance of fluctuating ocean transports in efforts to understand the seasonal cycle, for most regions.

It should be noted that the maps of Fig. 16 address the local behavior of the ocean over cells of 1° × 1°. Both for heat and salt–freshwater, surface fluxes usually account for a larger part of the seasonal cycle budget over wide regions (Fig. 17, bottom panels) than they do over 1° × 1° cells (top panels). This contrast is observed because fluctuations in air–sea fluxes correlate over wide regions more than do fluctuations in ocean fluxes. However, even for wide regions, the seasonal cycle can only be roughly approximated as a two-term balance (between storage and air–sea fluxes) for *T* poleward of the subtropics (bottom left panel). Even for wide regions, contributions of ocean fluxes to the seasonal cycle are relatively large for *T* in the tropics and subtropics, and for *S* everywhere (bottom panels).

### c. Concluding remarks

The seasonal cycle of heat and salt–freshwater content in OCCA is consistent with ocean in situ observations (Argo, etc.) and satellite observations (SST and altimetry). Additional dynamical constraints prevent overfits of data noise, so that the new atlas does not show large unphysical seasonal fluctuations in *T* and *S* (see section 6a). Unlike univariate maps based on a purely statistical model, OCCA includes a full budget for the seasonal cycle of heat and salt–freshwater content. Hence it provides a physical explanation for the observed seasonality in heat and salt–freshwater content, which is a combination of surface fluxes (adjusted from NCEP data) and ocean fluxes (resolved and parameterized physics).

The present data analysis (section 6b) stresses the importance of fluctuating ocean fluxes, which largely contribute to the seasonal cycle. Tentative applications of the new atlas include studies of mode water seasonality and its dynamics. For example, the North Atlantic Eighteen Degree Water is investigated in separate papers using OCCA (Maze et al. 2009; Forget et al. 2009, manuscript submitted to *J. Phys. Oceanogr.*). The OCCA ocean flux fields might also provide a practical way to define side boundary conditions in regional models. They are thought of as a reasonable first guess inferred from global sets of observations to improve upon thanks to regional observational efforts.

## 7. Summary and discussion

This paper exploits a near-global mean monthly ocean atlas (OCCA) synthesizing a variety of ocean observations for the 3-yr period from December 2003 through November 2006. OCCA is shown to be acceptably close to satellite altimeter observations, satellite sea surface temperature data, and in situ Argo profiles. OCCA thus adequately synthesizes all three datasets. Additional value of the new atlas follows from the use of a GCM interpolation framework. First for any individual quantity, observed or otherwise, the atlas takes account of the whole suite of observations as well as dynamical constraints. A key point is that the additional constraints do not preclude a close fit to each individual dataset. Second, the various mapped variables are all consistent with one another according to known dynamical principles, namely a discretized version of the primitive equations. Third, maps of observed quantities are augmented with derived estimates for quantities that are not directly observed. These estimates allow for an extended interpretation of observations and additional applications of the atlas. The atlas is publicly available online (www.ecco-group.org) and will later be extended beyond 2006.

This paper sheds new light on two oceanographic problems that are the subject of vigorous ongoing research. First, improving combined estimates of the ocean circulation and the geoid (e.g., Wunsch and Stammer 2003) remains a challenge. The present analysis suggests that the RMS uncertainty range for modern geoid estimates ought not to exceed 3.5 cm at scales of 300 km and more. The results imply that erroneous assumptions about the deep dynamic topography (or the reference-level circulation) can yield significant errors of a few centimeters in combined circulation and geoid estimates. Second, current field experiments, such as the Climate Variability and Predictability (CLIVAR) Mode Water Dynamic Experiment (CLIMODE), aim to refine the understanding of seasonal cycle budgets. The present observational synthesis suggests that seasonal fluctuations cannot in general be well approximated as a one-dimensional balance between air–sea fluxes and storage into the ocean. Ocean fluxes (i.e., advection plus mixing) are found to be key contributors to seasonal cycle budgets, especially for salt–freshwater and for local budgets.

One decisive difference between OCCA and previous atlases (such as *WOA01*) is that known dynamics are accounted for in OCCA. The primitive equations are indeed encoded in the GCM, which is used as an extensive box model including volume, momentum, and tracer conservation, and geostrophy has a limit case, etc. In contrast, ordinary objective mapping relies on a purely statistical model of the ocean behavior—a time–space decorrelation model typically. This omission of dynamics all together is obviously a very crude approximation. GCMs, while being imperfect, are regarded as refined approximations compared with decorrelation models.

A concern with GCMs is that model errors tend to grow as time evolves, and may grow in convoluted nonlinear ways. For example, in the North Atlantic, high-latitude overflow errors are expected to result in broad midlatitude circulation errors, 3 yr later or so (see Doscher et al. 1994; Gerdes and Köberle 1995). For that matter, our pragmatic approach (as in Forget et al. 2008b) is to confine model error growth by stopping integration of the GCM after each seasonal cycle. This approach readily allows a close fit to observations, which is the defining quality of an atlas, while imposing dynamics in the seasonal cycle. Model error is thus confined, even though it is not precluded. A notorious, long-lasting issue is the Gulf Stream path: many GCM solutions misplace it, or show a distorted momentum balance there (see Stammer et al. 2004). The OCCA dynamical balances may be a reasonable approximation in general, but there likely are such cases of distorted balances that will require further refinement.

A common limitation of global ocean atlases is their time/space resolution. The OCCA atlas consists of mean monthly fields on a 1° × 1° mesh, implying that sharp and eddying ocean features can only be represented in a smooth form. The atlas omits synoptic and small-scale signals, and the fields of Forget and Wunsch (2007) provide an upper bound estimate for this omission error for *T* and *S*. For that matter and beyond, the most outstanding issue may be our limited quantitative understanding of errors. First, estimates of omission–representation errors need to be refined to better account for synoptic and small-scale signals that global observing systems do not fully resolve. Second, it is likely that observing system deficiencies have not yet all been uncovered, or accounted for in quality control procedures. Third, errors in ocean, surface layer, sea ice, etc. models remain poorly understood and accounted for. Fourth, the extent to which the various ocean state variables may be retrieved from available observations largely remains to be elucidated. Along those lines, the error budget of ocean state estimates, and OCCA specifically, is the subject of an ongoing investigation, which will be reported elsewhere.

## Acknowledgments

Carl Wunsch, Patrick Heimbach, colleagues at MIT, and the two anonymous reviewers greatly helped improve this manuscript. This work was supported by the ECCO GODAE group (www.ecco-group.org), with funding from the National Ocean Partnership Program (NOPP) and the National Aeronautics and Space Administration (NASA). The Argo profile data were collected and made freely available by the International Argo Project (http://www.coriolis.eu.org). The SEaOS profile data were collected by the Southern Elephant Seals as Oceanographic Samplers program (http://biology.st-andrews.ac.uk/seaos/). Microwave OI SST maps are produced by Remote Sensing Systems (www.remss.com). The CLS/AVISO sea level products, which were used for comparison purposes, were produced by SSALTO/DUACS and distributed by AVISO with support from CNES (http://www.aviso.oceanobs.com).

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**(

**,**

**,**

**,**

**,**(

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

### APPENDIX A

#### Interpolation and Averaging Procedures

All interpolation procedures are done within the MIT GCM ECCO framework (Stammer et al. 2002; Wunsch and Heimbach 2007). They consist of the minimization of a least squares distance to observations:

where **x** are the adjustable parameters (hereafter control vector), **x*** _{b}* are a prior estimate for

**x**,

**y**

*is the observation vector,*

_{o}**y**(

**x**) are estimate counterparts to

**y**

*, and 𝗕*

_{o}^{−1}and 𝗥

^{−1}formally are the inverses of error covariance matrices. In practice,

*B*is implemented as a covariance operator using a diffusion operator (after Weaver and Courtier 2001). This covariance operator is set to favor large-scale adjustments in

*x*(see below for details). Here 𝗥 is a diagonal matrix involving only model–data error variance estimates (

*σ̃*

^{2}using notation in section 4). The compact

**y**(

**x**) notation may include dynamics (e.g., GCM time stepping) and interpolation to observed points. We shall distinguish two levels of complexity.

In the GCM interpolation procedure, **x** consists of initial conditions (for *T* and *S*) and weekly mean atmospheric state variables (used in bulk formula forcing; see appendix B), **y*** _{o}* includes all observations (from Argo, altimeters, etc.) and

**y**(

**x**) is the GCM counterpart to

**y**

*. The Weaver and Courtier (2001) operator is set to favor scales of 400 × cos(*

_{o}*ϕ*) km and 300 × cos(

*ϕ*) km (where

*ϕ*is latitude) for initial condition adjustments and atmospheric variable adjustments, respectively. The GCM adjoint is used to force a least squares fit of the GCM trajectory to the observation vector, through iterative adjustments (≈30 iterations) of the control vector.

In a simplified interpolation procedure (hereafter OI), the very GCM dynamical core is omitted and there is no time dependency. In this case, **y**(**x**) boils down to an interpolation at observed points and information is only *propagated* in space by the 𝗕 covariance operator. In a slightly more complex procedure (hereafter OIF), time dependency is represented by a Fourier series expansion, and *x* consists of fields for the various Fourier coefficients. When using a small number of annual harmonics, OIF yields a smooth seasonal cycle. In both OI and OIF, each quantity of interest (e.g., the temperature field) is treated individually. This is in contrast to the GCM interpolation procedure, where various physical quantities are coupled through equations of motion and thermodynamic relations. In OI and OIF procedures, the Weaver and Courtier (2001) operator is set to favor scales of 400 × cos(*ϕ*) km (*ϕ* is latitude). The least squares fit is achieved (iteratively using an adjoint) as for GCM interpolation.

The OCCA atlas results from applications of the full GCM interpolation procedure. The initial result consists of three optimized GCM solutions. Each solution estimates the oceanic state over one calendar year plus the previous 4 months, hence covering slightly more than one seasonal cycle. Solutions for two consecutive years have a 4-month overlap period. First, a 3-yr time series is compiled from the three solutions, with a gradual transition during the overlap period. The gradual transition procedure simply consists of a weighted average of overlapping solutions with time variable weights. The weights are *α* = (*t* − *t*_{0})/(*t*_{1} − *t*_{0}) and *β* = 1 − *α*, where *t* is time, *t*_{0} is 1 September, and *t*_{1} is 31 December. Second, the time series is further averaged to produce mean monthly atlas fields. All atlas fields (temperature, fluxes, mixed layer depth, etc.) are time averaged in the same way.

The choice to perform three 16-month optimizations rather than one 3-yr optimization aims at confining model error accumulation. Given necessarily imperfect GCM dynamics, one expects spurious GCM drifts. In the simplest case when model error is linear, the integrated effect would scale as the experiment duration. For the purpose of producing a mean monthly atlas, continuous resolution of the seasonal cycle is important, but continuous resolution of interannual fluctuations is inessential. Confinement of model error is not to be mistaken for reduction of model error though, which is an important long-term goal beyond the scope of this paper.

The simplified interpolation procedures (OI and OIF) are applied for various purposes, as mentioned in the text body. OI is first applied to compute sample mean *date of observation* maps shown in Fig. 1. In this case, the observation vector **y*** _{o}* consists of pointwise sample mean dates computed from the Boyer et al. (2006) dataset and recent datasets (Argo, XBT, etc.) within 1° × 1° cells. Second, OIF is applied with two annual harmonics to smooth the seasonal cycle of the

*WOA01*and

*WOA05*atlases. In this case, the observation vector

**y**

*consists of monthly atlas fields (with no other data involved). The result for*

_{o}*WOA01*is referred to as

*Fourier truncated WOA01*in the article body. Third, and for

*WOA05*only, OI is applied to adjust the annual mean fields to best fit in situ observations over the period from 2004 to 2006. The same is done for the two annual harmonics using OIF. In this case, the observation vector

**y**

*consists of in situ observations (Argo, XBT, etc.) of temperature or salinity. The result is referred to as*

_{o}*WOA05*–Argo blend in section 2 and Table 1.

### APPENDIX B

#### General Circulation Model Setup

The general circulation model is the MIT GCM (Marshall et al. 1997; Adcroft et al. 2004) and its adjoint is obtained by automatic differentiation (Giering and Kaminski 1998; Heimbach et al. 2005). Aside from the core ocean dynamics, the present experiments include a bulk formula atmospheric surface layer scheme (Large and Yeager 2004) and a sea ice thermodynamic model (Hibler 1980) as implemented in the MIT GCM.

The near-global model domain extends from 80°S to 80°N, has a resolution of 1° in the horizontal, 50 levels in the vertical, and a 1-h time step. The model setup uses a third order upwind advection scheme, an Adam–Bashforth explicit time stepping, along with an implicit scheme for vertical diffusion. The model setup uses hydrostatic and Boussinesq approximations, and a linear implicit free surface. For subgrid-scale process parameterizations, the model setup includes KPP vertical mixing (Large et al. 1994), GM eddy thickness diffusivity (Gent and Mcwilliams 1990), vertical and isopycnic tracer diffusivity, vertical and horizontal friction, and a quadratic bottom drag (coefficients listed in Table B1). The model uses free-slip side boundary conditions and no-slip bottom boundary conditions. The bulk formulas’ scheme handles freshwater as a virtual salt flux. The bulk formula atmospheric state variable inputs are air temperature, specific humidity, wind velocity, shortwave downward flux, and precipitation. Six-hourly atmospheric state variables from the NCEP reanalysis are used as a starting point.

### APPENDIX C

#### Dynamic Topography Computations

Two estimates of time-mean surface dynamic topography from the Rio et al. (2005) study are used here: (i) an estimate including in situ information (hereafter SDT_{Rio}) for the 1993–99 period (hereafter dT1 period) that is distributed by CLS (http://www.cls.fr); (ii) an estimate *not* including in situ information (hereafter ) for the 1993–2004 period (hereafter dT2 period) that was provided to the ECCO GODAE group by M. H. Rio. SDT_{Rio} is only used for discussion (in section 5a), while is also part of the GCM interpolation problem. Both SDT_{Rio} and are based on the CLS01 mean sea surface and the EIGEN-GRACE-03S geoid. SDT_{Rio} is an adjusted version of to account for in situ data (see Rio et al. 2005). Formation of the OCCA atlas involves an analogous adjustment estimated in a different framework.

To compare the two adjustments (in section 5), complication arises from different choices of the reference period (dT1 or dT2) for the SDT itself and for sea level anomaly (SLA) data. Essentially, dT1 is the CLS standard and dT2 is the ECCO GODAE standard, and a necessary correction is computed as

where SLA_{CLS} is the weekly mapped multimission SLA obtained from AVISO/CLS (DT Maps of Sea Level Anomaly (MSLA) reference product), and SLA^{dT} denotes the time average of SLA over dT. The adequate Rio et al. (2005) SDT adjustment is then computed as

and can equally be associated with dT1 or with dT2.

For the new atlas, the analogous adjustment is computed as

where SDT_{OCCA} is the OCCA time-mean dynamic topography estimate for the (dT3) period from 2004 to 2006. The third term on the right-hand side is the CLS estimate of the time-mean anomaly over dT3 referenced to dT1. Then is the same anomaly but referenced to dT2, and finally is an estimate of the SDT over dT2 that we compare to .

Aside from the reference period issue, both ΔSDT_{Rio} and SLA_{CLS} are smoothed in space to allow comparisons with the new atlas (sections 4b and 5a). Using a diffusion operator, scales smaller than 300 × cos(*ϕ*) km (i.e., three GCM grid points) are filtered out, since the GCM cannot properly resolve these scales. While this smoothing simply is convenient for comparison purposes, we have no reason to doubt the MDT_{Rio} or AVISO maps’ skill at smaller scales.

## Footnotes

*Corresponding author address:* Gaël Forget, Dept. of Earth, Atmospheric and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA 02139. Email: gforget@mit.edu

^{1}

In this paper, the term “univariate” simply denotes that a physical quantity is individually estimated, unlike in GCM interpolation.

^{2}

The “actual estimate error” is not to be confused with a “formal error estimate.”

^{3}

Average values between the 152-Sv (1 Sv ≡ 10^{6} m^{3} s^{−1}) barotropic streamline and the Antarctica shelf.

^{4}

Median values between 40°S and 40°N.

^{5}

Top to bottom integrated value.

^{6}

Values are standard deviation of daily transports.

^{7}

Maximum value of the time-mean overturning streamfunction at 27°N.