Ocean heat content (HC) is one of the key indicators of climate variability and also provides ocean memory critical for seasonal and decadal predictions. The availability of multiple operational ocean analyses (ORAs) now routinely produced around the world is an opportunity for estimation of uncertainties in HC analysis and development of ensemble-based operational HC climate indices. In this context, the spread across the ORAs is used to quantify uncertainties in HC analysis and the ensemble mean of ORAs to identify, and to monitor, climate signals. Toward this goal, this study analyzed 10 ORAs, two objective analyses based on in situ data only, and eight model analyses based on ocean data assimilation systems. The mean, annual cycle, interannual variability, and long-term trend of HC in the upper 300 m (HC300) from 1980 to 2009 are compared.
The spread across HC300 analyses generally decreased with time and reached a minimum in the early 2000s when the Argo data became available. There was a good correspondence between the increase of data counts and reduction of the spread. The agreement of HC300 anomalies among different ORAs, measured by the signal-to-noise ratio (S/N), is generally high in the tropical Pacific, tropical Indian Ocean, North Pacific, and North Atlantic but low in the tropical Atlantic and extratropical southern oceans where observations are very sparse. A set of climate indices was derived as HC300 anomalies averaged over the areas where the covariability between SST and HC300 represents the major climate modes such as ENSO, Indian Ocean dipole, Atlantic Niño, Pacific decadal oscillation, and Atlantic multidecadal oscillation.
Ocean reanalyses (ORAs) are now routinely produced at operational centers around the world particularly for initializing the ocean as part of dynamical seasonal forecast systems (Balmaseda et al. 2009). The current generation of operational ORAs has improved in quality because of advances in data assimilation schemes, improvements in ocean models, increased model resolutions, and dramatic improvements in the global ocean observing system (Behringer and Xue 2004; Usui et al. 2006; Zhang et al. 2007; Balmaseda et al. 2008; Drévillon et al. 2008; Yin et al. 2011; Vernieres et al. 2012).
Operational ORAs are now routinely used at national climate centers for El Niño–Southern Oscillation (ENSO) monitoring, and prediction, efforts. It is also becoming evident that modes of ocean climate variability other than ENSO also have significant influence over different areas of the globe, and therefore, there is an increasing requirement to monitor, predict, and understand such modes (Xue et al. 2010).
Increasing use of ORAs in the areas of societal relevance, however, raises questions about their fidelity and suitability: what are the similarities and differences among ORAs; could better monitoring products, together with errors bars, be developed based on an ensemble of ORAs? Indeed, the OceanObs’09 Conference called for an intercomparison of operational ORAs, and development of ocean monitoring products based on an ensemble of operational ORAs (Xue et al. 2010). As a follow up to that recommendation, the focus of this paper is to analyze the depiction of the upper ocean heat content variability from an ensemble of ORAs.
The upper-ocean heat content (HC) provides the ocean memory for potential predictability for climate variability on interannual to decadal time scales. For example, HC in the tropical Pacific is the ocean memory for the long-lead predictability of ENSO (Ji et al. 1998; Balmaseda and Anderson 2009; Xue et al. 2000; Clarke and Van Gorder 2003). HC variability is also associated with the Indian Ocean dipole (IOD) (Saji et al. 1999; Rao et al. 2002) and SST variability in the southwestern Indian Ocean (Huang and Kinter 2002; Xie et al. 2002), as well as the Atlantic Niño (Zebiak 1993; Ruiz-Barradas et al. 2000; Keenlyside and Latif 2007). An improved monitoring of HC variability can help us better understand, and forecast, those SST modes.
On decadal to multidecadal time scales, HC variability in the North Pacific is found to be associated with decadal variability of SST, for example, the Pacific decadal oscillation (PDO; Mantua et al. 1997), and ecosystem indicators along the west coast of North America (Mantua et al. 1997; Di Lorenzo et al. 2008). HC variability in the North Atlantic has been linked to the Atlantic meridional overturning circulation (AMOC) (Hakkinen and Rhines 2004; Böning et al. 2006; Huang et al. 2011) and Atlantic multidecadal oscillation (AMO) (Zhang 2008). Therefore, an accurate ocean initialization of HC in coupled climate models could be an important factor in enhancing skill in decadal climate prediction efforts (Meehl et al. 2009; van Oldenborgh et al. 2012).
A common practice for model-based operational ORAs is to assimilate various ocean observations into an ocean model that is generally forced by fluxes at the air–sea interface from an independent atmospheric reanalysis (Balmaseda et al. 2009). For ocean data assimilation a hierarchy of methods has been adopted that range from optimal interpolation (OI) (Balmaseda et al. 2008), three-dimensional variational (3DVAR) methods (Behringer et al. 1998; Usui et al. 2006), and Kalman filter (Drévillon et al. 2008; Keppenne et al. 2008; Yin et al. 2011). Several groups have also been developing coupled ocean and atmosphere data assimilation systems (Zhang et al. 2007; Sugiura et al. 2008; Saha et al. 2010).
Accuracy of HC analyses based on ocean data assimilation systems can be affected by many factors that include uncertainties in surface fluxes, biases in ocean models, deficiencies in data assimilation schemes, and changes in ocean observing systems. Several studies have been performed to compare HC variability from multiple ocean data assimilation systems. For example, Carton and Santorelli (2008) examined the consistency of the decadal variations of the ocean heat content in the upper 700 m (HC700 hereafter) in 1960–2002 in nine ocean reanalyses. They found some significant disagreements in the representation of climate anomalies, particularly in regions of poor historical observation coverage like the Southern Hemisphere but also in other areas. Lee et al. (2009) and Stammer et al. (2010) provided a summary of results of a comparison of ocean reanalysis carried out by the Climate Variability and Predictability (CLIVAR)–Global Synthesis and Observations Panel (GSOP) community. Stammer et al. (2010) examined the consistency of the decadal variations in HC700 in the past 50 years in 11 ocean reanalyses and noted the surprisingly large spread across products toward the end of the time series despite the fact that this period is best observed. Recently, Zhu et al. (2012) conducted an ensemble estimation of HC variability in the tropical Atlantic using six ocean reanalysis products. HC700 from multiple ocean reanalysis products have also been compared in the context of studying the impacts of different XBT bias correction schemes (Levitus et al. 2009; Lyman et al. 2010; Giese et al. 2011). Those studies suggest that it is very difficult to remove all XBT biases associated with the assumptions used for the XBT fall rate.
In this study, we compared HC in the upper 300 m (HC300 hereafter) variability from 10 ORAs, two objective analyses based on in situ data only, and eight analyses based on ocean data assimilation systems. The mean, annual cycle, interannual variability, and long-term trend of HC300 from 1980 to 2009 from each of the products are compared. Our study differs from earlier studies in that we only included ocean reanalysis products that are updated operationally on a real-time basis. Our goal is not only to conduct a comprehensive comparison of HC300 analyses but eventually to extend the comparison in real time in support of climate monitoring, climate assessment and climate prediction. A goal of this study is to lay down the foundation and implementation strategy for the development of operational ensemble ocean heat content monitoring products.
In our study, we did not attempt to quantify errors in each analysis by validating against in situ observations. Instead, we focused on documenting and understanding the differences and similarities among HC300 analyses, with the premise that errors in the analysis of HC300 can thus be estimated. We selected one analysis as the reference and calculated various statistics relative to the reference analysis. We also used the ensemble method of Lee et al. (2009), in which the signal is estimated as the temporal standard deviation of the ensemble mean, and noise as the spread across different ocean analyses. Zhu et al. (2012) also demonstrated that the ensemble mean of multiple analyses is superior to any individual analysis.
We investigated the influence of evolving ocean observing systems on the spread of HC300 analyses. In principal, when more ocean observations are available, as different ORAs are better anchored by increasing observations, ensemble spread should decrease. Therefore, ensemble spread can be used to measure influences of ocean observing systems on the convergence of multiple ocean reanalysis products. The ensemble spread also helps us identify outliers, which would encourage analysis centers to understand the sources for large departures in subsequent studies.
Section 2 describes the 10 ORAs, and the methodologies used in the comparison. Section 3 discusses the ensemble spread among HC300 analyses, and how it varies with changing ocean observing systems. Section 4 discusses the capability of the current generation of operational ocean reanalyses in describing HC300 variability associated with the major modes of climate variability on interannual, decadal and long-term time scales. The summary and discussion are given in section 5.
2. Operational ORAs and analysis methods
The main features of the 10 operational ORAs used in the study are summarized in Table 1, and a brief description of each ORA is provided in the appendix. Information on operational ORAs can also be found in the review paper by Balmaseda et al. (2009) and in individual references for each ORA.
In the study, HC300 is defined as the average temperature in the upper 300 m. To facilitate comparison, the monthly mean HC300 for 1980–2009 was first calculated on the native model grid of each ORA and then linearly interpolated to a common 1° by 1° grid. We defined a common ocean mask that is shared by all ORAs, which roughly covers the domain from 70°S to 70°N.
The monthly climatology for each ORA was obtained by taking the first three harmonics of the monthly averages over the period 1985–2009. The first few years (1980–84) were not included in the climatology since a large initial drift was found in a few ORAs, for example, in the Global Ocean Data Assimilation System (GODAS), Climate Forecast System Reanalysis (CFSR), and Geophysical Fluid Dynamics Laboratory (GFDL). The anomalous HC300 (HC300a) was derived by removing the 1985–2009 climatology for each ORA separately. For the comparison of HC300 it is useful to choose a reference analysis. Since the Met Office analysis (EN3) is based on an objective analysis of the in situ data only, and has a monthly resolution, it was used as the baseline for various comparison statistics described below. It is important to emphasize that this does not mean that EN3 is viewed as “the truth.”
In the comparison of HC300 across ORAs, we frequently looked at the spread. For spatially and temporally varying HC300, the spread among ORAs was calculated as the standard deviation of the departure of ORAs from the ensemble mean at each grid point and each month in 1980–2009. For temporally fixed maps such as the climatological mean, linear trend, and anomaly correlation of HC300, the spread was simply the standard deviation from the ensemble mean for the respective quantity.
In the analysis, we also frequently used a signal to noise ratio (S/N hereafter), where signal is the temporal standard deviation of the ensemble mean of HC300 over a period, and noise the average of standard deviation of the spread of ORAs relative to the ensemble mean over the same period. When S/N is larger than 1, the signal in the quantity analyzed is larger than the noise and more confidence can be placed in the ensemble mean as an estimate of the climate signal.
To measure how well HC300 variability in each ORA resembles that of EN3, the anomaly correlation (AC) with EN3 HC300 was calculated. To identify the regions where HC300 variability is highly correlated with SST variability, HC300 from each ORA was correlated with the monthly mean observed SST derived from the weekly Optimal Interpolation SST version 2 (OISST) (Reynolds et al. 2002). Because HC300 has slower variability than SST does, it represents a potential mechanism for predictability of SST variations on time scales of a season and longer. As expected, those regions with high HC300–SST covariability are associated with ENSO, PDO, AMO, IOD and Atlantic Niño.
3. Uncertainties in HC300 analysis and the influences of ocean observing systems
Estimation of HC300 can be influenced by many factors, of which the distribution of in situ observations plays a critical role. The temperature observations assimilated in different ORAs generally include profiles from XBTs, CTDs, mooring arrays and Argo profiling floats (see Table 1). The number of temperature profiles per month increased significantly around 1985, and increased dramatically around 2002 with the introduction of Argo data (see Fig. 11 from Saha et al. 2010). The altimetry sea surface height (SSH) data, available since 1993, was included in some of the ORAs, for example, in, European Centre for Medium-Range Weather Forecasts (ECMWF), GFDL, Global Modeling and Assimilation Office (GMAO), Mercator-Ocean (MERCATOR), and Japan Meteorological Agency (JMA) (Table 1). The analysis in this section aims to estimate uncertainties in the mean of HC300 analysis and how they vary with changing ocean observing systems.
a. Analysis of mean HC300
The HC300 from EN3 was used as the baseline measure for comparison of the mean HC300 over 1985–2009. Figure 1 shows the mean HC300 for EN3, the offset in the mean for the other nine ORAs relative to EN3, the ensemble mean of the offsets for the model-based ORAs, and the ensemble spread of model-based ORAs.
The mean HC300 from EN3 shows that the average temperature in the upper 300 m of the tropical oceans is generally above 17°C except in the eastern tropical Pacific and Atlantic, and in the southwest tropical Indian Ocean. Those regions with low mean HC300 are also where covariability between HC300 and SST anomalies are highest and will be discussed in next section.
The mean offsets of HC300 from EN3 are generally small in the tropical Pacific except in the GFDL analysis, which is colder along the equatorial waveguide (Fig. 1g). The mean offsets are also small in the extratropical northern oceans except near the western boundary currents and Labrador Sea. In the tropical Indian Ocean, the mean offsets are generally small except in the Bay of Bengal and near the west coast of Australia (Figs. 1f,g,i). However, the mean offsets are relatively large in the extratropical southern oceans and the tropical Atlantic where observations are sparse. The ensemble mean offsets of model-based ORAs are generally smaller than those of any individual ORA, but remain large near the western boundary currents, Gulf of Mexico and extratropical southern oceans (Fig. 1k). The model-based ORAs also show large spread in HC300 climatology in some of these poorly observed areas, like the Antarctic Circumpolar Current (Fig. 1l).
To quantify the temporal variations in the total HC300 in each ORA, 2-yr running-means of HC300 averaged in different ocean basins are shown in Fig. 2. The spread among analyses generally decreases with time, particularly after the early 2000s when Argo data became available. Temporal variations are in best agreement in the northern oceans where the S/N ratio is 2.3, but they are barely significant in the tropical Pacific (Fig. 2a) and tropical Atlantic (Fig. 2c) where the S/N ratio is 1.3 and 1.1, respectively. The S/N ratio in the tropical Indian Ocean and southern oceans is very poor.
b. Spread of HC300 and influences of ocean observing systems
The spread of HC300 across ORAs varies with time due to changes in the coverage of in situ observations. Figure 3 shows the spread averaged over each ocean basin calculated, where spread is computed either using all 10 ORAs (thick black line) or with only nine ORAs with one ORA withheld (thin lines). The latter is an attempt to identify the ORAs that are outliers and have a large contribution to the spread.
Removing one ORA generally reduces the spread. For the tropical Atlantic removing CFSR reduced the spread significantly through the whole period, indicating CFSR as an outlier. CFSR is also an outlier in the tropical Indian Ocean during 1980–85. GODAS and GFDL were identified as outliers in the southern oceans, probably because both analyses have the largest variability in those regions (see discussion related to Fig. 9). In the tropical Pacific, the spread decreased substantially around 1985 and 2003 when the TAO and Argo data began to be assimilated into ORAs. The peaks around 82–84 and 97–99 are associated with the 1982/83 and 1997/98 El Niño events. The advent of the Argo data is clearly seen in the tropical Indian Ocean (Fig. 3b) and extratropical southern oceans (Fig. 3e), but is not as clear in the tropical Atlantic (Fig. 3c) and extratropical northern oceans (Fig. 3d).
The temporal variations of the spread can be partially attributed to the temporal variations of data counts in each ocean basin (Fig. 4). The data counts in Fig. 4 include the instrument casts that contain temperature measurements for at least one depth in the upper 300 m, regardless of the depth(s) at which the measurement(s) was (were) made. Temperature profiles are from bottle/CTD, XBT, MBT, moored buoy and Argo; those from gliders and drifting buoys (other than Argo) were excluded as they may skew the counts due to their very large numbers of profiles in a very small area. Each temperature profile is either from a completely separate instrument, or separate casts of the same instrument.
In the tropical Pacific, the data counts increased substantially in the early 1990s because of the completion of the Tropical Atmosphere Ocean (TAO) array, and then increased further in the early 2000s when the Argo data became available (Fig. 4a). The Argo data also dramatically increased the data counts in the tropical Indian, the tropical Atlantic, and the southern oceans (Figs. 4b,c,e). However, the data counts decreased gradually from 1986 to 2010 in the northern oceans (Fig. 4d) and from 1986 to 2002 in the global ocean (Fig. 4f) because of the decrease of counts from bottles/CTDs and XBTs (not shown).
To see how observations influence the spatial distribution of uncertainties in the mean HC300, Fig. 5 shows the ensemble spread of HC300 analyses in 1985, and the corresponding data counts, the changes of the ensemble spread from 1985 to 1995, and from 1995 to 2006, and corresponding changes in the data count. The data counts are the number of profiles with temperature data in the upper 300 m that passed quality control checks in EN3 averaged into one degree grid boxes.
In 1985, the ensemble spread is above 0.4°C in parts of the tropical Pacific, tropical Indian Ocean and tropical Atlantic, and above 1°C near the western boundary currents and extratropical southern oceans (Fig. 5a), while the data coverage is mostly confined to north of 10°S (Fig. 5d). In 1995, the full implementation of the TAO mooring array dramatically increased the data density in the tropical Pacific (Fig. 5e), which led to a reduction in the ensemble spread by as much as 40% (Fig. 5b). The data increase in the Bay of Bengal, subpolar North Atlantic, and subtropical southern oceans also led to a reduction in the ensemble spread in those regions (Figs. 5b,e). The spread increased across the tropical Indian Ocean and the North Pacific, and corresponds to a decrease in data counts over those regions.
The introduction of Argo in the early 2000s significantly improved the global data coverage (Fig. 5f). From 1995 to 2006, the data coverage increased over most of the global ocean except in the South China Sea, the western tropical Atlantic and midlatitude North Atlantic. There is a good correspondence between the increase of data counts and reduction in the spread, particularly in the tropical Indian Ocean, equatorial eastern Atlantic, and high-latitude southern oceans (Fig. 5c). Despite the reduction in the ensemble spread, the spread in 2006 is still above 0.7°C near the western boundary currents, Gulf of Mexico, western tropical Atlantic, and extratropical southern oceans, indicating a need for further improvement of HC300 analysis in those regions.
c. Uncertainties in HC300 analysis near the equator
Because of the global influence of oceanic variability in the equatorial oceans, time variations of uncertainties in HC300 analysis in the equatorial ocean zones are investigated further. Figure 6 shows the time evolution of the spread of HC300, and the offsets of HC300 from EN3 for each ORA in the equatorial Pacific from 1980–2009. The spread of HC300, that is, the standard deviation around the ensemble mean of the 10 ORAs, reduces significantly around 1993 when the TAO mooring array was almost fully implemented (Fig. 6a). The spread is largest during the 1982/83 and 1997/98 events, during which the differences between EN3 and National Oceanographic Data Center (NODC) are particularly large (Fig. 6b). A further examination of the NODC and EN3 analysis suggests that the differences between EN3 and NODC are probably caused by differences in the TAO buy data retained in each analysis. A future study is needed to understand differences in the quality control procedures in the two analyses.
The GMAO and Bureau of Meteorology (BOM) analyses agree very well with the EN3, particularly after 1990. The GFDL analysis has a significant cold offset before 1993 that gets smaller after 1993 and is mostly gone after 2003 (Fig. 6g). The cold offset identified here was attributed to a too large observation error used in the analysis, which has been fixed in a later version of the analysis. The CFSR analysis has a cold offset before 1999 and a warm offset afterward, which does not diminish in the recent period (Fig. 6f). The switch from cold to warm offset in the CFSR is related to a sudden change in the quality of the trade winds, which were too strong before 1998 but became close to observations once the Advanced Television and Infrared Observation Satellite Operational Vertical Sounder (ATOVS) satellite observations were assimilated in the atmospheric analysis (Xue et al. 2011; Zhang et al. 2012). The GODAS has some warm offsets that do not diminish in the recent period (Fig. 6c). Other analyses (ECMWF, JMA, MERCATOR) have some warm offsets relative to EN3 early in the period.
The spread in the analyses of HC300 in the equatorial Indian Ocean is generally small in the central basin, but relatively large in the western and eastern Indian Ocean (Fig. 7a), which is consistent with the standard deviation of HC300 (Fig. 9a). The spread is reduced significantly around 2003 in the eastern Indian Ocean when the two RAMA mooring lines were installed at 80.5° and 90°E. The NODC analysis agrees well with EN3 except during 1997–2001 when the NODC analysis was colder (Fig. 7b). It is interesting that the model-based ORA offsets from EN3 are mostly negative (positive) in the eastern (western) tropical Indian Ocean before 1997, which suggests that the east-west thermocline slope (Fig. 1a) is smaller in all models than in the EN3 analysis. A reduction in offsets around 1997 might be related to an improvement in atmospheric reanalysis winds due to assimilation of satellite scatterometer winds (Chelton and Freilich 2005). In fact, a similar reduction in the offset around 1997 is found in the control simulation (without data assimilation) of the BOM analysis (Yonghong Yin, personal communication). Soon after, differences between the analyses diminish significantly because of the implementation of Argo and RAMA.
The spread in the analyses of HC300 in the equatorial Atlantic Ocean is unevenly distributed across the basin, with the largest spread confined east of 10°W (Fig. 8a). The spread in the eastern Atlantic was significantly less after 2006, while the spread in the west-central Atlantic decreased considerably after 1995 (Fig. 8a). Note that the Pilot Research moored Array in the Tropical Atlantic (PIRATA) moorings implemented at 35°, 23°, and 10°W on the equator around 1998 have contributed to reduce the spread in the western equatorial Atlantic but had very little effect in the eastern Atlantic. In this region it is only after 2005/06 with the advent of Argo that the spread is substantially reduced.
The offsets relative to EN3 are generally small except in JMA, BOM, GFDL, and CFSR. The JMA and BOM analyses have significant cold offsets in the eastern equatorial Atlantic throughout the whole period, while the GFDL analysis is colder in the early period. The CFSR has alternating warm and cold offsets and is related to model drifts, and discontinuities across the six streams in which the CFSR was executed (Xue et al. 2011). In summary, the HC300 analysis in the equatorial Atlantic has the largest spread (uncertainty) among the three equatorial oceans, but spread reduced significantly in the last few years of the analysis period that may be due to enhanced observations from the Argo floats and PIRATA moorings.
4. Interannual, decadal, and long term variability in HC300
a. Analysis of variability of monthly mean HC300 anomaly
The magnitude of HC300 variability is quantified based on the standard deviation (STD) of monthly mean HC300 anomalies. Further, the STD of ORAs is compared with that from EN3 (Fig. 9). The STD from EN3 shows that HC300 variability is relatively large in the eastern equatorial Pacific and western tropical Pacific. The STD in the tropical Indian Ocean is weaker than that in the tropical Pacific with maximum amplitude in the eastern and western tropical Indian Ocean and south of the equator near 10°S. The STD in the tropical Atlantic is the weakest among the tropical Oceans. In the northern oceans, the STD is relatively high near the Kuroshio-Oyashio Extension (KOE), the Gulf Stream (GS), and subpolar North Atlantic. Large HC300 variability is also observed in the midlatitude South Atlantic near the coast of South America and South Africa.
The STD of NODC is about 15%–30% weaker than that of EN3 (Fig. 9b), which may be partially related to the fact that the NODC (EN3) is a seasonal (monthly) analysis. The STD of the model-based ORAs is generally larger than EN3 near the western boundary currents and in the southern oceans. A few exceptions are as follows: 1) the STD in GMAO, MERCATOR, and BOM is generally weaker than EN3 in the southern oceans, which might be related to a strong relaxation to climatology (Table 1); 2) the STD in MERCATOR is weaker than that of EN3 in the North Pacific; and 3) the STD in CFSR is more than 100% stronger than EN3 in the tropical Atlantic, which is related to the six analysis streams and model drift within each analysis stream (Xue et al. 2011). Near the Antarctic Circumpolar Current (ACC), the HC300 variability is particularly strong in GODAS and GFDL, probably because they do not include a relaxation to climatology (Table 1). In data sparse regions such as near the ACC, the analysis is essentially the model run forced by reanalysis fluxes (GODAS) or the fully coupled model run (GFDL).
Consistency of HC300 variability in different ORAs with EN3 is measured by anomaly correlation (AC). Figure 10 shows the AC of each analysis with EN3, the ensemble mean and the spread in AC calculated over 1985–2009. The AC is generally low near the GS and KOE. This is understandable since the resolution of the climate models cannot represent the ocean dynamics necessary for analyzing the correct location or strength of the GS and KOE (Kwon et al. 2010). On the other hand, the data-only analyses, that is, EN3 and NODC, may not represent the location and strength of the GS and KOE either due to their coarse resolution.
The ACs averaged in different ocean basins are summarized in Table 2. In terms of average ACs, the consistency in the North Pacific (mean = 0.67, spread/mean = 9%, i.e., spread is 9% of the mean) is slightly higher than that in the North Atlantic (mean = 0.62, spread/mean = 11%). In the tropical Pacific, the average AC ranges from 0.74 in CFSR to 0.86 in BOM, and the spread/mean is 4%, indicating a high consistency among ORAs. In the tropical Indian (Atlantic) Ocean, the mean AC is 0.58 (0.45) and the spread/mean is 8% (20%). Therefore, consistency in the tropical Atlantic is the lowest of the three tropical oceans. The low consistency in the tropical Atlantic is largely attributed to the poor AC (0.24) in CFSR, which is much lower than the ACs in other ORAs (>0.42). The consistency of AC in the three extratropical southern oceans is generally low, indicating large uncertainties in HC300 variability over those regions.
We next quantify the strength of the interannual variability in the ensemble mean of monthly mean HC300 versus the amplitude of the spread across the HC300 analyses over 1985–2009. In this calculation, the CFSR was excluded since doing so significantly reduces the ensemble spread in the tropical Atlantic (Fig. 3c). Further, since the S/N in the tropical Indian was improved if NODC and EN3 were also removed, the signal, the noise and S/N shown in Fig. 11 were calculated from the ensemble of 7 ORAs.
The STD of the ensemble mean HC300a (Fig. 11a) is comparable to that of EN3 HC300a (Fig. 9a) in the tropical Pacific, southwestern tropical Indian Ocean and near the western boundary currents, but is much smaller in the western/eastern tropical Indian Ocean, the tropical Atlantic and midlatitude southern oceans, indicating ensemble averaging damps variability in those regions. On the other hand, the noise is the largest near the western boundary currents and midlatitude southern oceans (Fig. 11b). If the S/N is larger than 1, a climate signal represented by the ensemble mean is regarded as significant, and so can be used to monitor an aspect of climate variability.
The S/N in the tropical Pacific is above 4, which is lower than the estimation by Zhu et al. (2012). The S/N is more than 1.5 in the eastern and southwest tropical Indian Ocean, suggesting that the common variability in ORAs is larger than the noise. In contrast, the S/N is near 1 in the tropical Atlantic where there is hardly a consistent signal in HC300 across the ORAs. We also noticed that the S/N is larger than 2 in the eastern North Pacific, where both the signal and noise are relatively weak, and subpolar North Atlantic, where the signal is strong and the noise is weak. We will show next that the regions where the S/N is large coincide with the areas where HC300 variability is also significantly correlated with SST variability.
b. Interannual variability of HC300 associated with major climate modes
The level of anomaly correlation between HC300 and SST helps identify the regions where HC300 variability may have significant influences on SST variability, therefore providing a source of potential predictability for climate variability on seasonal and longer time scales. Figure 12 shows that monthly mean HC300a is significantly correlated with SST anomaly (SSTa) in many areas of the global ocean. The AC is generally above 0.6 in the central and eastern tropical Pacific, eastern North Pacific and subpolar North Atlantic. The high AC in the eastern tropical Pacific is attributable to strong air-sea couplings and ocean dynamics associated with ENSO, while the high AC in the North Pacific and subpolar North Atlantic is probably associated with atmosphere-driven strong vertical mixing and convection resulting in variations in mixed layer, SST and HC300. The high AC between HC300 and SST in the subpolar North Atlantic provide a source of long-lead SST prediction skill (van Oldenborgh et al. 2012). Other regions where the AC is moderately high (>0.4) are the western tropical Pacific, southwest tropical Indian Ocean, and eastern tropical Atlantic.
The spread in AC (Fig. 12i) is relatively small in the tropical Pacific, North Pacific, southwest tropical Indian Ocean, and subpolar North Atlantic, suggesting that uncertainties in HC300-SST covariability are relatively small and the relationship is robust. However, the spread is relatively large near the GS and KOE, in the southeast tropical Indian Ocean, tropical Atlantic, and extratropical southern oceans, indicating that the correlations between HC300a and SSTa have large uncertainties in those regions.
The boxes in Fig. 12 highlight regions of maximum covariability between HC300a and SSTa that are associated with ENSO, IOD, and Atlantic Niño. The ACs averaged over the boxes as well as the box definitions are listed in Table 3. In the Niño-3 region (5°S–5°N, 150°–90°W), the mean AC between HC300a and SSTa is 0.77 and a spread/mean of 5% (Table 3). The AC in the western tropical Pacific is much lower, with a mean of 0.39 and a spread/mean of 11%. In the southeastern tropical Indian Ocean, the AC ranges from 0.16 in EN3 and NODC to 0.39 in MERCATOR. We will show later that the EN3 and NODC analysis are problematic in this region. The AC in the southwest tropical Indian Ocean has a mean of 0.44 and a spread/mean of 11%. The moderately high correlation is expected in the southwest tropical Indian Ocean where the mean thermocline is shallow, and strong upwelling and vertical mixing tend to bring subsurface temperature anomalies near to the surface (Xie et al. 2002). In the equatorial eastern Atlantic, the AC ranges from 0.35 in CFSR to 0.60 in GFDL, and has a mean of 0.50 and a spread/mean of 18%.
With the goal of monitoring interannual variability of HC300a associated with ENSO, IOD, and the Atlantic Niño in real time, a set of climate indices were derived as the average HC300a in the boxes shown in Fig. 12. In Fig. 13, the box average HC300a from individual analysis (thin line) and the ensemble mean (thick black line) are shown along with the S/N displayed on the top right corner of each panel. The S/Ns were derived for the period 1985–2009, excluding the early years when some ORAs had considerable initial drift.
The S/N is high in the eastern (5.5) and western (6.5) tropical Pacific, consistent with that shown in Fig. 11c. The tropical Pacific HC300a has a decadal shift (Figs. 13a,b): variability is much weaker, and the equatorial western Pacific is much warmer after 2000. In the southeast tropical Indian Ocean, the large HC300a variability is associated with the IOD events in 1982, 1994, 1997 and 2006 (Fig. 13c). Note that the NODC and EN3 analyses, without the benefit of the model dynamics responding to the observed surface forcings to compensate for sparse observations in the region, missed the large positive anomaly in 1999 (Fig. 13c). The S/N is 2 when all 10 ORAs were included, while it is 2.7 if NODC and EN3 were excluded, suggesting that model-based ORAs were more consistent. Similarly, the S/N increased from 2.4 to 2.7 in the southwest tropical Indian Ocean (Fig. 13d) without NODC and EN3. The HC300a in the southwest tropical Indian Ocean are highly correlated with those in the tropical Pacific, indicating strong influences of ENSO on HC300a variability over this region. In addition, the large positive HC300a in 1994 and 2006 are associated with IOD events (Rao et al. 2002). The S/N in the equatorial eastern Atlantic (Fig. 13e) is 1.3 when all 10 analyses were included, but it is 1.6 when CFSR was excluded, indicating that HC300a variability associated with the Atlantic Niño is marginally represented by the ensemble mean of ORAs. The Atlantic Niño has higher frequency variability than ENSO during 1980–2000, and there is a warming trend over much of the period from 2000–08.
c. Decadal and long-term variability
Large multidecadal variability of the heat content in the upper 300 m was first reported a decade ago (Levitus et al. 2000). Recent studies suggest that earlier estimations of multidecadal heat content variability, and a global warming trend, might be contaminated by large biases in XBT data (e.g., Domingues et al. 2008). It is very difficult to remove all biases associated with the assumption used for the XBT fall rate. In fact, uncertainty due to choice of XBT bias correction dominates variability across different estimates in global heat content estimation during 1993–2008 (Lyman et al. 2010). Among the 10 ORAs used in the study, two ORAs (GFDL and GMAO) include XBT corrections, while others do not. Although, we do not discuss the impacts of XBT bias on decadal variability and long-term trend of HC300, we focus on a discussion of ensemble mean signal and spread for the long-term trend.
The quasi-global (70°S–70°N) means of HC300 from the 10 ORAs are shown in Fig. 14a. There are differences among ORAs, but there is a consensus that the global ocean increased from 1984 to 1992 followed by a short cooling episode in 1992/93, and then increased from 1994 to 2003/04, followed by flattening or a decrease. The GFDL analysis has the strongest warming trend from 1984 to 2003, and the steepest cooling trend from 2003 to 2009. The BOM analysis also has a substantial cooling trend after 2003 (Fig. 14d). It is noted that the GODAS and CFSR are the warmest in 1980/81, which is related to the initial spin up (Xue et al. 2011). The cooling in the early 1980s may not be reliable due to the large spread (Fig. 14c), but the cooling in 1992/93 are most likely related to the dimming affect of aerosols (Church et al. 2005) associated with eruptions of the Mt. Pinatubo volcano in June 1991.
The variations in HC300a are similar to those in total HC300 except that the estimates are closer together (Fig. 14b) due to removal of differences in mean values. It is noticed that the GMAO and CFSR analyses have a third cooling period around 1997/98 when other analyses remain fairly constant. Most of the analyses have maximum values around 2002–04, while CFSR, JMA, and GMAO have peak values in 2006, 2007, and 2009 respectively (Fig. 14d).
The spread in the mean HC300 (solid line in Fig. 14c) decreases from 1980 to 1995, but increases from 1995 to 2003. The spread then decreases after 2003 when the Argo array reaches a global coverage. It is interesting to note that the spread of HC300a is generally smaller than that of total HC300, and the increase of the spread during 1995–2003 is also less pronounced. The increase of the spread during 1995–2003 is significantly reduced if GFDL is not included. The S/N for HC300a in 1985–2009 is 1.8, indicating the global mean HC300 variability is reasonably well represented by the set of ORAs.
The positive trend in the quasi-global mean of HC300 has contributions from different ocean basins. We calculated the linear trends for the 1985–2009 and 1993–2009 periods separately, and found the trend patterns are similar over the two periods while the trend over 1993–2009 has a stronger amplitude. Figure 15a shows the ensemble mean of the linear trends over 1993–2009. The uncertainty of the trend is estimated by the ratio of the ensemble mean of the trends to the ensemble spread of the trends from 10 ORAs (Fig. 15b). The mean trend can be viewed as a robust estimate if the ratio is large. Next we will discuss HC300 variability in the regions where the trend is robust.
For the North Pacific, Xue et al. (2011) calculated the dominant empirical orthogonal functions (EOFs) of HC300a using the NODC, GODAS, and CFSR products for the period 1979–2009. They found that the first two EOFs and their principal components (PCs) are highly consistent among the three ocean reanalysis products, and also, are very similar to those derived from sea surface height in an ocean model simulation of the 1950–2004 period (Di Lorenzo et al. 2008). The first PC, representing the temporal variation of the PDO (Di Lorenzo et al. 2008), is well described by the HC300a differences in the eastern North Pacific and west-central North Pacific shown as boxes in Fig. 15. It is seen that the PDO was in a positive phase from 1980 to 1998, except for a brief period in early 1990s, and switched to a negative phase around 1999, staying negative since then (Fig. 16a). The S/N of the PDO is 4.2, indicating a high consistency among ORAs, except that the GMAO has slightly weaker amplitude than other analyses in the early 2000s.
The decadal variability of HC300 in the subtropical North Pacific and South Pacific can be described by the HC300a differences in the eastern and western subtropical Pacific shown as boxes in Fig. 15. The S/N for the subtropical North Pacific and South Pacific is 4.9 and 5.4 respectively. The two indices are highly correlated with each other (Figs. 16b,c), and are negatively correlated with that in the tropical western Pacific (Fig. 13b). They are dominated by strong interannual variability from 1980 to 1999, switching to a colder state around 1999, followed by a cooling trend overlain with weak interannual variability.
The subtropical southern Indian Ocean has a weak warming trend from 1993 to 2009, which is robust despite considerable spread among analyses (Fig. 16d). The warming trend in HC300 is consistent with that in SSH in a model simulation of Han et al. (2010), who suggest that the warming trend started in early 1960s. The subtropical North Atlantic has a weak warming trend from 1984 to 2003 and a moderate cooling trend from 2003 to 2009 (Fig. 16e). The S/N is 2, indicating a better consistency in the subtropical North Atlantic than in the subtropical South Indian Ocean. The first EOF of altimetry sea surface height shows a weakening of the North Atlantic subpolar gyre since 1995 (Hakkinen and Rhines 2004). Consistent with Hakkinen and Rhines (2004), the subpolar HC300a increased substantially around 1995, reached a peak value around 2006, and decreased since then (Fig. 16f). The S/N is 4, indicating a high consistency among HC300a analyses.
5. Summary and conclusions
The availability of multiple ORAs that are now routinely produced by operational and research centers provides an opportunity to assess uncertainties in HC300 analyses, to help identify gaps in observing systems as they impact the quality of ORAs (and therefore climate model forecasts), and help identify deficiencies in data assimilation schemes. Multiple ORAs also provide the basis for development of real-time multimodel ensemble HC300 (and other) monitoring products.
We analyzed 10 ORAs, two objective analyses based on in situ data only, and eight model analyses based on ocean data assimilation systems. The mean, annual cycle, interannual variability, and long-term trend of HC300 in 1980–2009 from the 10 ORAs was compared. In this comparison, we selected the EN3 analysis as the reference (without regard to its quality) to assess similarities and differences among analyses. We also defined a signal to noise ratio (S/N), where signal is the temporal standard deviation of the ensemble mean of HC300, and noise the temporal standard deviation of the spread of ORAs.
The mean HC300 in 1985–2009 was highly consistent among ORAs north of 30°S except near the Gulf Stream and Kuroshio-Oyashio Extension, in the tropical Atlantic, and some regional seas (Fig. 1). The model-based ORAs show large spread in their mean value in some of these poorly observed areas, like the Antarctic Circumpolar Current. This is a strong indication that both improvements in the first guess (model+fluxes) as well as in observation coverage are needed to represent H300 in those areas.
The spread in the estimate of the mean HC300 in different ocean basins generally decreased with time and reached a minimum in the early 2000s when the Argo data became available (Fig. 2). There is a warming trend in the northern oceans that is well represented by the ensemble where the S/N is 2.3. In the tropical Pacific and tropical Atlantic, the S/N is 1.3 and 1.1, respectively. There is no robust trend before 2000 because of the large spread among analyses, however there is a rapid warming in early 2000s and a cooling trend after 2005 that are consistent among analyses. The S/N in the tropical Indian Ocean and southern oceans is too low (0.6, 0.2 respectively) to detect any decadal variability or long-term trend.
The ensemble spread, calculated as the standard deviation from the ensemble mean at each grid point and each month, was found to be strongly influenced by the distribution of in situ data. There was a good correspondence between the increase of data counts and reduction of the spread (Fig. 5). In the tropical Pacific, the averaged spread decreased substantially around 1985 and 2003 when the TAO and Argo data began to be assimilated into ORAs (Fig. 3a). The impacts of the Argo data was also clearly seen in the averaged spread in the tropical Indian Ocean (Fig. 3b) and extratropical southern oceans (Fig. 3e), but not as clear in the tropical Atlantic (Fig. 3c) and extratropical northern oceans (Fig. 3d).
The standard deviation of monthly mean HC300 anomaly was generally consistent among ORAs except near the western boundary currents, in the tropical Atlantic and southern oceans (Fig. 9). In the southern oceans, the BOM, GMAO, and MERCATOR have much weaker STD than that of EN3, probably related to a strong relaxation to climatology (Table 1). On the other hand, the HC300 variability is particularly strong in GODAS and GFDL near the Antarctic Circumpolar Current (ACC) since they do not include a relaxation to climatology.
Uncertainties in HC300 anomaly were quantified by anomaly correlation (AC) with EN3 (Fig. 10) and the S/N ratio (Fig. 11). The AC with EN3 is generally high in the tropical Pacific, tropical Indian Ocean, North Pacific, and North Atlantic, but it is low in the tropical Atlantic and extratropical southern oceans where observations are very sparse. The S/N was above 4 in the tropical Pacific, while it was more than 1.5 in the eastern and southwest tropical Indian Ocean and near 1 in the tropical Atlantic. The S/N was larger than 2 in the eastern North Pacific and subpolar North Atlantic. The regions where the S/N is large coincide with the areas where the AC with EN3 is high.
As for a larger spread over the tropical Atlantic, we propose two factors contributed to the relatively low S/N in the tropical Atlantic: 1) the amplitude of signal is weaker and is about 1/3 of that in the tropical Pacific (Fig. 11a); 2) the amplitude of noise is the largest in the tropical oceans (Fig. 11b). The large noise in the tropical Atlantic may be attributed to large uncertainties in the surface forcings used in ocean reanalyses. Balmaseda and Mogensen (2010) showed that the improved surface forcings from the ERA-Interim significantly improved sea surface height anomaly correlation with satellite observations in two ocean reanalysis systems. It is expected that the signal will enhance and the noise will reduce in the tropical Atlantic when better, and more consistent, surface forcings are used in the next generation ocean reanalyses.
Monthly mean variability in HC300a was significantly correlated with SST anomaly (SSTa) in many areas of the global ocean (Fig. 12). The AC was generally above 0.6 in the central and eastern tropical Pacific, North Pacific, and subpolar North Atlantic. Other regions where the AC was moderately high (>0.4) were the western tropical Pacific, southwest tropical Indian Ocean, and eastern tropical Atlantic.
With the goal of monitoring interannual variability, a set of climate indices were derived as HC300a averaged over the areas where the covariability between SST and HC300 was associated with ENSO, IOD, and Atlantic Niño (Fig. 13). The S/N was high in the eastern (5.5) and western (6.5) tropical Pacific. In the southeast tropical Indian Ocean, the S/N was 2 when all 10 ORAs were included, while the S/N was 2.7 if NODC and EN3 were excluded, indicating a consistency in the forced ocean models when data is sparse in this region. The S/N in the equatorial eastern Atlantic was 1.3 when all 10 analyses were included, but it was 1.6 when CFSR was excluded. This indicates that HC300 variability associated with the Atlantic Niño is marginally represented by the ensemble mean of the ORAs.
HC300 has large multidecadal variability and long-term trends. The consensus among ORAs suggests that the mean HC300 in 70°S-70°N increased from 1984 to 1992 followed by a short cooling episode in 1992/93, and then increased from 1994 to 2003/2004, followed by flattening or a decrease. A set of climate indices were derived based on HC300a average in those regions where the trend was robust.
Our analysis suggests that the suite of operational ocean reanalyses can be used to monitor climate signals related to variability on interannual to decadal time scales. The ensemble mean of climate indices represent our best knowledge of climate variability, while the ensemble spread provides an estimation of uncertainties. We envision this activity could be extended to routine exchange of climate information from ORAs, and implementation of operational monitoring of climate indices using multimodel analyses in the near future.
We thank support from the NOAA Climate Observation Division of Climate Program Office for this study. The EN3 analysis was supported by the Joint DECC/Defra Met Office Hadley Centre Climate Program (GA01101); The GMAO analysis was supported by the NASA Modeling, Analysis and Prediction (MAP) Program; The BOM analysis was supported by the Australian Managing Climate Variability Program. We also thank Dr. David Behringer, Dr. Oscar Alves, Dr. Zeng-Zhen Hu, and Dr. Caihong Wen for helpful comments on the manuscript.
Summary of Operational Ocean Reanalyses
A brief description of each ORA is provided below.
The operational GODAS was implemented at National Centers for Environmental Prediction (NCEP), National Oceanic and Atmospheric Administration (NOAA)/United States, in 2003 (Behringer and Xue 2004), providing oceanic initial conditions from 1979 to present for the NCEP Climate Forecast System (Saha et al. 2006). GODAS is based on the Geophysical Fluid Dynamics Laboratory’s Modular Ocean Model version 3 (MOM3) (1° with 0.3° equatorial refinement and 40 levels) with forcings from the NCEP–Department of Energy (DOE) Reanalysis and a 3DVAR data assimilation scheme. At first, only observed temperature and synthetic salinity profiles (Behringer et al. 1998) were assimilated; subsequently, beginning in August 2006, satellite altimetry sea surface height (SSH) was added (Behringer 2007). In the 3DVAR system the influence of the satellite SSH is projected onto the subsurface temperature and salinity fields by including the difference between the satellite SSH and the model dynamic height, linearized in temperature (T) and salinity (S), in the cost function. Only the variable part of the satellite SSH is assimilated, which is accomplished by removing the annual averages for the years 1993–99 from the satellite SSH and the model dynamic height. The model SST is nudged strongly to the NOAA weekly OI SST (OISST, hereafter) product (Reynolds et al. 2002) (with a relaxation coefficient equivalent to time scale of 5 days). The GODAS is updated daily in real-time with a 2-day delay, and pentad and monthly averages are used for real-time global ocean monitoring products (http://www.cpc.ncep.noaa.gov/products/GODAS).
A new reanalysis for the atmosphere, ocean, sea ice and land, known as the CFSR (Saha et al. 2010), was implemented in March 2011, and was used to provide atmospheric and oceanic initial conditions from 1979 to present for the NCEP Climate Forecast System version 2. The oceanic component of CFSR includes many advances over that of GODAS: (i) the MOM version 4 (MOM4) ocean model at ½° with ¼° equatorial refinement and 40 levels with an interactive sea ice model, (ii) 6-h coupled model forecast as the first guess, (iii) inclusion of the mean climatological river runoff, (iv) high spatial (0.5°) and temporal (hourly) model outputs, and (v) strong relaxation to the NOAA daily OISST product (Reynolds et al. 2007). See details in Xue et al. (2011). The execution of the CFSR was completed with six parallel streams with one year overlap across different streams (Saha et al. 2010).
c. GFDL ocean reanalysis
The GFDL ocean data assimilation system consists of an Ensemble Kalman Filter applied to GFDL’s second generation fully coupled climate model (CM2.1) (Zhang et al. 2007). The ocean component is the MOM4 configured with 50 vertical levels and 1° horizontal resolution, telescoping to ⅓° meridional spacing near the equator. The atmospheric component has a resolution of 2.5° longitude and 2 o latitude with 25 vertical levels. The first guess is from a fully coupled model where the atmosphere component is constrained by winds, sea level pressure and temperature data from the NCEP–DOE reanalysis. For the ocean component, observed temperature (T) and salinity (S) profiles, and the OISST (Reynolds et al. 2002) were assimilated. Chang et al. (2011) constructed so-called pseudo salinity profiles from the weighted least squares procedure that minimizes the misfits between the predetermined vertical coupled T–S EOF modes and the observed temperature and altimetry SSH data. These pseudo salinity profiles are assimilated into the GFDL analysis. The altimetry SSH data is partially used for the generation of pseudo salinity profiles, but not directly assimilated to the current analysis (see details in http://www.gfdl.noaa.gov/ocean-data-assimilation).
d. GMAO ocean reanalysis
The GMAO reanalysis uses the Goddard Earth Observing System coupled atmosphere–ocean general circulation model (GEOS-5) and is based on MOM4 (0.5° with ¼° equatorial refinement and 40 levels) and the GEOS-5 AGCM (1° × 1.25° with 72 levels) model. The atmosphere is constrained by the atmospheric fields from the Modern Era Retrospective-analysis for Research and Applications (MERRA) (Rienecker et al. 2011) and the first guess for the ocean comes from a coupled forecast (Vernieres et al. 2012). The ocean data assimilation uses a multivariate ensemble optimal interpolation (EnOI) to infer background-error covariances from a static ensemble of 50 model state–vector EOFs. Observed temperature (T) and salinity (S) profiles, along-track sea level anomalies (SLA) from Archiving, Validation, and Interpretation of Satellite Oceanographic data (AVISO) are assimilated daily. The T/S profiles are from the uncorrected EN3 dataset. The XBT temperature profiles have been corrected according to Levitus et al. (2009). Synthetic salinity profiles are generated for all temperature-only observations. Ten percent of observations (selected randomly) from the NOAA daily OISST of Reynolds et al. (2007) are assimilated daily. The climatological sea surface salinity is also assimilated to compensate for errors in freshwater input from precipitation and river runoff. Ten percent of T and S profiles randomly extracted from the World Ocean Atlas 2009 gridded climatology (Antonov et al. 2010; Locarnini et al. 2010) are also assimilated every five days with very small weights. The SLA assimilation corrects only the model surface height field. To assimilate the SLA a constant offset between the AVISO anomalies and the model’s surface height was removed. This offset was calculated from the mean difference between the AVISO anomalies and an ocean analysis which did not assimilate SLA. More details are available at http://gmao.gsfc.nasa.gov/research/oceanassim/.
e. ECMWF ocean reanalysis
The ECMWF ocean reanalysis, referred to as ORA-S3, has been operational since August 2006, providing ocean initial conditions for the ECMWF seasonal and monthly forecasts since March 2007. The ORA-S3 is based on the Hamburg Ocean Primitive Equation (HOPE) model (1° with 0.3° equatorial refinement and 29 levels), and a 3D multivariate OI scheme to assimilate temperature, salinity, altimeter-derived sea level anomalies, and global sea level trends (Balmaseda et al. 2008). It spans the period 1959–present. Forcing fluxes are provided by the 40-yr ECMWF Re-Analysis (ERA-40) until 2002 and ECMWF operational analysis thereafter. ORA-S3 uses the daily interpolated SST used in ERA-40 from 1959 until 1982 and the OISST (Reynolds et al. 2002) thereafter. The model SST is nudged strongly to the externally analyzed SST product (with a relaxation coefficient equivalent to time scale of 3 days). The T/S profiles are from the EN2 dataset (Ingleby and Huddleston 2007) where the XBTs were not corrected, for the historical period (prior to operational running) and directly from the GTS ever since the system became operational. Both altimeter-derived sea level anomalies maps from AVISO as well as the global sea level trends are assimilated. The model sea level anomalies are computed respect a mean dynamic topography from an equivalent experiment assimilating T and S only. Because of the Bousinesque approximation, the model mean sea level is blind to the global thermal expansion. Therefore, the global mean sea level is removed from model and altimeter maps before assimilation. These “zero-global-mean” sea level anomalies are assimilated using the scheme of Cooper and Haines (1996) as described in Balmaseda et al. (2008) and Vidard et al. (2009). The global mean sea level from AVSIO maps is then compared with the model-diagnosed steric height, and the difference is added as a uniform freshwater flux. The assimilation procedure uses an adaptive bias correction to correct systematic model errors. A selection of historical and real-time ocean analysis products can be found online (http://www.ecmwf.int/products/forecasts/d/charts/ocean).
f. MERCATOR reanalysis
The MERCATOR ocean reanalysis, also referred to as PSY2G2, covers the 1979–present time period and is used at Météo-France for coupled seasonal forecasts. The PSY2G2 is based on the Océan Parallélisé version 8.2 (OPA8.2) ocean model in the ORCA2 global configuration at 2° with 0.5° equatorial refinement and 31 levels. The surface forcing used is ERA-40 until 2002 and ECMWF operational analyses afterward. In situ temperature and salinity profiles, OISST maps are assimilated weekly using a fixed basis reduced order Kalman filter (Drévillon et al. 2008). Along-track altimetry SLA data is assimilated together with a mean dynamic topography (or mean sea surface height) derived from a long model simulation without data assimilation. The underlying assumption is that the mean sea level will not be changed by data assimilation but only its variability. The background covariance error is inferred from a static seasonally variable ensemble of ~300 ocean states.
g. JMA ocean reanalysis
The JMA ocean reanalysis, also referred to as the multivariate ocean three-dimensional variational estimation (MOVE)/ Meteorological Research Institute Community Ocean Model (MRI.COM) global system [MOVE/MRI.COM-G (Usui et al. 2006)], was implemented in March 2008. The analysis system covers the quasi-global ocean (75°S–75°N) at 1° grid with 0.3° equatorial refinement and 50 vertical levels. The analysis scheme adopted in the MOVE system is a multivariate 3DVAR analysis scheme with vertically coupled temperature–salinity EOF modal decomposition of a background error covariance matrix. Temperature and salinity profiles and along-track altimetry SSH anomaly are assimilated. The along-track SSH anomaly is assimilated along with the mean surface dynamic height (SDH) calculated from a preliminary analysis that assimilates temperature and salinity observations only (Usui et al. 2006). The Centennial in situ Observation-Based Estimates (COBE) Analysis SST (COBE-SST) product (Ishii et al. 2005) is assimilated as well. The system provides pentad and monthly fields from 1979 to present (http://ds.data.jma.go.jp/tcc/tcc/products/elnino).
h. BOM ocean reanalysis
The BOM ocean reanalysis, also called the Predictive Ocean Atmosphere Model for Australia (PEODAS) Ensemble Ocean Data Assimilation System (POAMA; http://poama.bom.gov.au is for the period from 1980 to present. It is an approximate form of ensemble Kalman filter system (Yin et al. 2011) wherein a single (central) analysis is computed utilizing an ensemble-based covariance and each of the perturbed ensemble members is nudged toward the central analysis to control the ensemble spread and constrain the mean. Both in situ temperature and salinity observations are assimilated every three days, and corrections to the ocean currents are generated based on the ensemble covariances. See details in Yin et al. (2011).
i. EN3 objective temperature analysis
The EN3 of the Met Office is an objective monthly temperature analysis based on in situ observations with a 1° grid and 42 levels (EN3_v2a; Ingleby and Huddleston 2007). A historical reanalysis for the period 1950 to present is available, and the real time updates have an approximately one-month delay (http://www.metoffice.gov.uk/hadobs/en3). This analysis does not include XBT bias corrections.
j. NODC objective temperature analysis
The NODC analysis is an objective seasonal temperature analysis based on in situ observations. The analysis is on a 1° grid with 16 levels ranging from the ocean surface to 700 m in depth from 1955 to 2009 and includes XBT temperature corrections (Levitus et al. 2009). However, for this comparison study, an uncorrected XBT version of the analysis was generated. That is, the NODC product shown in this paper is from a seasonal temperature analysis similar to that of Levitus et al. (2009) except that XBT temperature corrections were not included. Discrepancies between EN3 and NODC, therefore, should be due to factors other than XBT corrections.