This paper examines nine analyses of global ocean 0-/700-m temperature and heat content during the 43-yr period of warming, 1960–2002. Among the analyses are two that are independent of any numerical model, six that rely on sequential data assimilation, including an ocean general circulation model, and one that uses four-dimensional variational data assimilation (4DVAR), including an ocean general circulation model and its adjoint. Most analyses show gradual warming of the global ocean with an ensemble trend of 0.77 × 108 J m−2 (10 yr)−1 (=0.24 W m−2) as the result of rapid warming in the early 1970s and again beginning around 1990. One proposed explanation for these variations is the effect of volcanic eruptions in 1963 and 1982. Examination of this hypothesis suggests that while there is an oceanic signal, it is insufficient to explain the observed heat content variations.
A second potential cause of decadal variations in global heat content is the uncorrelated contribution of heat content variations in individual ocean basins. The subtropical North Atlantic is warming at twice the global average, with accelerated warming in the 1960s and again beginning in the late 1980s and extending through the end of the record. The Barents Sea region of the Arctic Ocean and the Gulf of Mexico have also warmed, while the western subpolar North Atlantic has cooled. Heat content variability in the North Pacific differs significantly from the North Atlantic. There the spatial and temporal patterns are consistent with the decadal variability previously identified through observational and modeling studies examining SST and surface winds. In the Southern Hemisphere large heat content anomalies are evident, and while there is substantial disagreement among analyses on average the band of latitudes at 30°–60°S contribute significantly to the global warming trend. Thus, the uncorrelated contributions of heat content variations in the individual basins are a major contributor to global heat content variations.
A third potential contributor to global heat content variations is the effect of time-dependent bias in the set of historical observations. This last possibility is examined by comparing the analyses to the unbiased salinity–temperature–depth dataset and finding a very substantial warm bias in all analyses in the 1970s relative to the latter decades. This warm bias may well explain the rapid increase in analysis heat content in the early 1970s, but not the more recent increase, which began in the early 1990s.
Finally, this study provides information about the similarities and differences between analyses that are independent of a model and those that use sequential assimilation and 4DVAR. The comparisons provide considerable encouragement for the use of the sequential analyses for climate research despite the presence of erroneous variability (also present in the no-model analyses) resulting from instrument bias. The strengths and weaknesses of each analysis need to be considered for a given application.
This is an examination of the temperature in the upper 700 m of the global ocean during the 43-yr period of 1960–2002 as represented in nine gridded historical analyses. Despite the limited and nonstationary character of historical temperature sampling and problems with instrument bias, a number of studies now show a gradual warming of the global oceans (e.g., Levitus et al. 2000, 2005; Carton et al. 2005; Ishii et al. 2006). The analyses used in these studies also show that the multidecadal trend is modified by rapid fluctuations in the warming rate with increases occurring in the early 1970s and again beginning in the 1990s (Fig. 1). These fluctuations have been interpreted as being either an artifact—the result of inadequacies of the data or analysis method (Gregory et al. 2004; AchutaRao et al. 2006; Gouretski and Koltermann 2007)—or real, for example, the result of volcanic eruptions (Ramaswamy et al. 2001; Church et al. 2005; Delworth et al. 2005), variations in solar emissions (White et al. 2003), or natural variability of the climate system (Miller and Schneider 2000; Lozier et al. 2008). Here we explore decadal heat content variability as it appears in a suite of nine historical analyses.
In individual ocean basins the character of decadal variability differs significantly from the global average. Observational and modeling studies in the North Pacific show decadal variability dominated by a cooling of subtropical SST in the central basin and in the Kuroshio–Oyashio Extension region, and warming in the east from the tropics to midlatitudes in the winter months beginning in 1976/77 (Mantua et al. 1997; Miller and Schneider 2000). These changes are followed by a reduction of heat in the waters of the western subtropical gyre north of 30°N associated with changes in the wind patterns, advection, and surface heat flux (Miller et al. 1998; Miller and Schneider 2000). Along the equator in the Pacific, Zhang et al. (1997) and Luo and Yamagata (2001) document a complex pattern of warm anomalies west of 160°E between 1965 and 1980 followed by a period of generally cool anomalies, while farther east warm anomalies are present during the late 1970s and again through the 1990s. These and subsequent papers link the appearance of subtropical and tropical heat anomalies through subduction and advective processes.
The North Atlantic, in contrast, shows rapid steady warming in the subtropics and midlatitudes and some cooling at subpolar latitudes, but with a basin-average rise of heat content according to Levitus et al. (2005) of 2 × 108 J m−2 (10 yr)−1. This basin-average warming must reflect a net downward surface heat flux of 0.7 W m−2. However, the decadal redistribution of heat within the basin has been ascribed to other processes. Dickson et al. (2000) examine cooling in the southern Labrador Sea and conclude that it results partly from decadal increases in convection and changes in freshwater input associated with winter storms. Studies farther east suggest that these changes in the Labrador Sea have acted to cool the broader subpolar gyre (Read and Gould 1992). Lozier et al. (2008) relate both the subpolar cooling and subtropical warming to changes in surface transports associated with changes of the North Atlantic Oscillation and associated shifts in the winter storm tracks. Farther north observational studies have reported a gradual warming of the Eurasian portion of the Arctic Ocean. Grotefendt et al. (1998) show the warming trend in the Barents Sea through the mid-1990s to be the result of a combination of local warming and heat advection resulting from a poleward extension of warm Atlantic Water. This warming trend apparently reversed in the late 1990s (Gunn and Muench 2001).
In the Southern Hemisphere decadal variations in heat content are horrendously undersampled. For example, of the 1.1 million temperature profiles collected in the 1960s, only 78 000 were collected south of 30°S [based on data described in Boyer et al. (2006)]. Alory et al. (2007) avoid this limitation by examining the much better-sampled SST dataset together with coupled model output and suggest that some of the most rapid increase of upper-ocean heat content during the past five decades has occurred in the South Indian Ocean. This suggestion is consistent with Willis et al. (2004) who detect substantial warming in the southern midlatitudes during recent years (1993–2002).
In addition to undersampling, two additional limitations of the historical observational dataset complicate the determination of decadal heat content variability. The first is the presence of changing observation bias resulting from the evolution of the observing system. The primary temperature dataset consists of 7.9 million vertical profiles. Prior to the early 1970s the majority were collected with a device called a mechanical bathythermograph (MBT), while between the 1970s and the early 2000s (the end of our period of interest) the majority were increasingly collected with the expendable bathythermograph (XBT; Boyer et al. 2006). The XBTs obtain their depth observations from the timing of the instrument fall rate. The equation converting fall time to depth that was used prior to the 1990s is known to underestimate depth by a few percent (Hanawa et al. 1995; Ingleby and Huddleston 2007). However, there are ongoing debates regarding the efficacy of the fall-rate corrections in current use, thus leaving open the possibility that the increase in the rate of warming in the early 1970s noted in Fig. 1 may result from observation bias (AchutaRao et al. 2006).
A second limitation of the historical temperature dataset is its changing vertical sampling. MBTs measure temperature above 280 m, while XBTs generally extend to either 400 or 700 m. Thus, of the 1.1 million profiles collected during the 1960s when MBTs were the primary instrument, only 100 000 extended to 500 m (Boyer et al. 2006). The impact of this shallow observation sampling on global heat content estimates has been considered recently by AchutaRao et al. (2006) and Gouretski and Koltermann (2007). Gouretski and Koltermann, for example, use spatial variability to estimate an uncertainty of 2.4 × 108 J m−2 in global heat content.
Another explanation for the presence of decadal variations in global heat content is that they are real and the result of volcanic aerosols injected into the stratosphere (a time series of global aerosols is included in Fig. 1). Modeling studies suggest a reduction of global ocean heat content of ∼1 × 108 J m−2 for recent major eruptions beginning in the year following the eruption (Delworth et al. 2005; Church et al. 2005). However, in addition to their effect on global heat content we may anticipate distinct latitudinal responses to different eruptions. Mount Agung, which erupted in 1963, distributed much of its aerosols into the Southern Hemisphere (Hansen et al. 2005), and El Chichón, which erupted in 1982, distributed its stratospheric aerosols into the Northern Hemisphere, while Mount Pinatubo, which is located nearly as far north as El Chichón, distributed the aerosols from its 1991 eruption roughly symmetrically about the equator.
In this study we compare 0-/700-m heat content computed from nine separate analyses of the historical temperature record. The methodologies used in the nine analyses vary, but can roughly be divided into three approaches. The first is the “no model” analyses, which use temperature observations to update a first guess provided by climatological monthly estimates of temperature. The second is the sequential data assimilation analyses, which march forward in time from a previous analysis using a numerical simulation of the evolving temperature and other variables produced by an ocean general circulation model. The simulation provides the first guess of the state of the ocean (temperature, salinity, etc.) at the next analysis time, while corrections are made to this first guess based on observations of variables such as temperature, salinity, or sea level. The third approach is four-dimensional variational data assimilation (4DVAR), which in this implementation uses the initial conditions and surface forcing as control variables to be modified in order to be consistent with the observations as well as a numerical representation of the equations of motion through iterative solution of a giant optimization problem [see Kalnay (2003) for a more detailed description of these three methodologies]. Comparison among these analyses allows us to address the question raised by Gregory et al. (2004) and AchutaRao et al. (2006) regarding the impact of the analysis methodology on the determination of decadal heat content in the global ocean.
2. Methods and data
Here we briefly review the datasets and methodologies used in constructing the nine analyses (summarized in Table 1). The sources of the subsurface temperature and salinity datasets are one or another version of the National Oceanic and Atmospheric Administration’s (NOAA’s) World Ocean Database (WOD). The release in 1998, WOD98, contains a total of 5.3 million profiles, including 2.1 million MBTs (Boyer et al. 2006, see their Table 1.2). WOD98, together with recent updates from the Global Temperature–Salinity Profile Program (GTSPP; Wilson 1998), forms the subsurface temperature dataset used in the Global Ocean Data Assimilation System (GODAS). By the WOD01 release (Conkright et al. 2002), which is used in Centre Européen de Recherche et de Formation Avancée en Calcul Scientifique [CERFACS; i.e., the European Centre for Research and Advanced Training in Scientific Computation; see Davey (2005)], Istituto Nazionale di Geofisica e Vulcanologia (INGV; see Davey 2005), Ishii et al.’s (2006) analysis (hereafter ISHII), Levitus et al.’s (2005) analysis (hereafter LEVITUS), and the U.K. Forecasting Ocean Assimilation Model (U.K. FOAM; see Bell 2000; Bell et al. 2004), the number of profile observations had grown by one-third to 7.0 million, with the number of MBTs increasing by 14%. Most of these analyses have also included recent data from GTSPP. The most recent release, WOD05, used in the Simple Ocean Data Assimilation (SODA; Carton and Giese 2008), has seen an additional 16% increase to 7.9 million profiles, while the number of MBTs increases by an additional few percent along with improvements in quality control.
As mentioned in the introduction, XBT data are known to contain a warm bias resulting from inaccurate modeling of the instrument fall rate. ISHII, Geophysical Fluid Dynamics Laboratory (GFDL; Sun et al. 2007), GODAS (Behringer 2005), INGV, LEVITUS, SODA, and U.K. FOAM all include a correction to this fall-rate equation provided by Hanawa et al. (1995). Three analyses—CERFACS, INGV, and U.K. FOAM—are part of the Enhanced Ocean Data Assimilation and Climate Prediction (ENACT) program (Davey 2005) and use a common dataset that includes some additional bias corrections (Thadathil et al. 2002). For XBT profiles collected after 1996 there is some confusion in the observational record regarding which XBTs have had a fall-rate correction applied and which have not (see Levitus et al. 2005 for a discussion). Analyses, such as GFDL, that show rapid warming after 1996 may possibly be suffering from a lack of universal fall-rate correction during this recent period (A. Rosati 2007, personal communication).
In addition to the profile data two remotely sensed datasets enter many of the analyses. The first is satellite altimeter sea level, which is available continuously since 1991 from a succession of satellites. CERFACS, Global Estimation of Circulation and Climate Experiment (GECCO; Köhl and Stammer 2008), GODAS, INGV, and U.K. FOAM all include these data. The second is satellite SST, which is available nearly continuously since 1981. For CERFACS, GECCO, GFDL, GODAS, and INGV these data enter indirectly through the use of a gridded SST product. SODA assimilates the nighttime SST observations directly. In addition to these datasets the analyses based on sequential data assimilation and 4DVAR use surface meteorological forcing provided by the 40-yr European Centre for Medium-Range Weather Forecasts (ECMWF) Re-Analysis (ERA-40; Uppala et al. 2005), the National Centers for Environmental Prediction–National Center for Atmospheric Research (NCEP–NCAR) Global Reanalysis 1 (NCEP-1) or NCEP/Department of Energy (DOE) Global Reanalysis 2 (NCEP-2; see Kalnay et al. 1996; Kanamitsu et al. 2002), or the Met Office operational product. Each meteorological product contains biases that must be considered when interpreting the analyses (e.g., Hines et al. 2000).
Next we briefly describe the methodologies used to construct the analyses. The two no-model analyses, ISHII and LEVITUS, begin with a first guess of the climatological monthly upper-ocean temperature based on climatologies produced by the NOAA/National Oceanographic Data Center (see Levitus et al. 2005; Ishii et al. 2006). For each analysis the differences between the observations and the first guess are mapped onto the analysis levels. LEVITUS uses the technique of Cressman (1959) and Barnes (1964) with a homogeneous decorrelation scale of 555 km to objectively map the differences between the temperature observations and the climatology onto a uniform grid. ISHII uses an alternative three-dimensional variational data assimilation (3DVAR) scheme to carry out this objective mapping with a somewhat smaller decorrelation scale in midlatitudes (300 km) that elongates in the zonal direction by a factor of 3 at equatorial latitudes.
We examine the following six analyses: CERFACS, GFDL, GODAS, INGV, SODA, and U.K. FOAM, which use sequential data assimilation. In sequential data assimilation an ocean general circulation model provides a first guess, while a version of optimal interpolation or 3DVAR is used to provide corrections to that first guess based on a set of empirically determined error covariances (see Kalnay 2003).
The three ENACT assimilation analyses are CERFACS, INGV, and U.K. FOAM. CERFACS and INGV use an ocean model based on Océan Parallélisé (OPA) version 8.2 (Madec et al. 1998) numerics with 2° × 2° horizontal resolution in the midlatitudes, decreasing to 2° × 0.5° in the deep tropics and 3DVAR assimilation. ERA-40 surface winds and climatological fluxes are used by both. CERFACS uses an assimilation algorithm based on 3DVAR, while INGV uses the mathematically similar optimal interpolation (Davey 2005). U.K. FOAM uses a different ocean model evolved from Cox (1984), with 1° × 1° × 33 level resolution and 16 levels in the upper 700 m. Surface forcing is provided by the Met Office numerical weather prediction winds and fluxes, along with climatological river discharge. The data assimilation is the timely optimal interpolation scheme of Bell (2000). In addition to the ENACT subsurface data, the analysis assimilates satellite altimetry using a modified version of the algorithm of Cooper and Haines (1996).
The GODAS analysis of Behringer (2005) uses an ocean model based on GFDL Modular Ocean Model (MOM) version 3 numerics with 1° × 1° horizontal resolution in midlatitudes, decreasing to 1° × 1/3° resolution in the deep tropics. The Arctic Ocean is not included in this grid. Surface forcing is provided by the NCEP-2 winds and fluxes. This analysis is the shortest considered here, spanning the period of 1979–2005. A related analysis produced by GFDL uses MOM version 4 numerics (Griffies et al. 2003) with a similar grid, but includes a complete Arctic Ocean. Winds, heat, and freshwater flux are provided by the NCEP–NCAR reanalysis (Kalnay et al. 1996). The assimilation is also a version of 3DVAR, but one that differs in detail from that used for GODAS.
The last of the sequential analyses is the Simple Ocean Data Assimilation, version 2.1, which uses a model based on the Parallel Ocean Program version 2.1 numerics with global average 0.25° × 0.25° resolution (0.25° × 0.4° resolution in the tropics) and a complete Arctic Ocean. This analysis is driven by ERA-40 winds through 2001, with an extension beyond 2001 using satellite scatterometer winds.
Finally, the GECCO 4DVAR analysis is based on the forward and adjoint versions of the Massachusetts Institute of Technology (MIT) ocean circulation model (Marshall et al. 1997) with 2° × 2° horizontal resolution, with initial flux estimates provided by the NCEP–NCAR reanalysis. Climatological temperature and salinity were among the constraints. Among the control variables are the initial conditions and atmospheric fluxes (Stammer et al. 2004).
Each analysis output is reduced to monthly averaged temperature (except LEVITUS, which is only available as an annual average) for whatever portion of our period of interest (1960–2002) is covered by the analysis (see Table 1). Heat content is then computed by integrating the temperature vertically 0/700 m using linear interpolation, while anomalies are computed relative to the 25-yr period of 1966–95. For the GODAS analysis, which does not span the full base period, anomalies are computed relative to the 1979–95 average.
To explore the potential contribution of volcanic aerosols to decadal variability of heat content we compare the difference in analysis heat content for 4-yr periods after and before the three major eruptions in 1963, 1982, and 1991. Because we have a limited number of eruptions to examine one concern is that the apparent response of the ocean will be affected by unrelated climate variability, particularly ENSO. In this comparison we attempt to mitigate the impact of ENSO by first carrying out a regression of the heat content analyses against the Southern Oscillation index. The part of the heat content variability that is correlated with this index at zero time lag is removed prior to computing the 4-yr differences. The impact of a linear trend in heat content is also removed prior to computing the 4-yr differences.
To aid in identification of heat content changes directly resulting from volcanic aerosols we compare the 4-yr differences computed from the analyses with similar difference calculations performed on output from GFDL’s Climate Model version 2.1(CM2.1) global coupled model (Delworth et al. 2005). To filter out the effects of random weather and climate variations from these coupled model differences we exploit the availability of multiple simulations and examine an average of five ensemble members.
Finally, as discussed in the introduction, there is serious concern that components of the observing system are biased and that bias in the observing system may lead to time-dependent bias within the analyses. To identify such bias within the analyses we compare the analyses to what we believe is an unbiased dataset, the WOD05 historical salinity–temperature–depth and conductivity–temperature–depth (CTD) dataset for the following three decades: 1970–99. This dataset contains 38 000 profiles in the North Atlantic and North Pacific. For each profile we vertically integrate to produce an estimate of heat content and then linearly interpolate the monthly heat content estimates from the analyses to the same geographic location to produce a set of analysis-minus-CTD observation differences. These differences are then composited into 10-yr averages. The interpretation of these statistics is discussed by Hollingsworth and Lonnberg (1989). No attempt is made to compute analysis-minus-CTD observation differences for LEVITUS because this analysis is only available as an annual average.
In most analyses global average heat content has two periods of rapid growth—in the late-1960s through early 1970s and again beginning in the early 1990s, separated by quiescent or cooling periods in the early to middle 1960s and in the 1980s (Fig. 1). The scatter among the individual analyses generally lies within ±1 × 108 J m−2. Because of the general similarity of the no-model and sequential analyses, we define an ensemble average based on these eight analyses. The amplitude of the decadal anomalies is 1–2 × 108 J m−2 (a 1 × 108 J m−2 change in 5 yr implies a net heating of 0.63 W m−2), while the multidecadal trend of the ensemble is 0.76 × 108 J m−2 (10 yr)−1 or 0.24 W m−2. The major exception is the GECCO 4DVAR analysis, which cools until the mid-1970s and then begins a multidecadal period of warming. Two of the analyses—U.K. FOAM and GFDL—show significant additional warming beginning in the mid-1990s.
The spatial distribution of the linear trend for each analysis is presented in Fig. 2. Several regions show rapid warming. These include the western half of the subtropical and midlatitude North Atlantic, the South Atlantic, the subtropical North and South Pacific, and the northern parts of the Indian Ocean, as well as the Indian Ocean south of 30°S. Interestingly, there is a cooling trend in the Kuroshio–Oyashio Extension region of the North Pacific in most analyses except GECCO. Several analyses have cooling trends in the central and western equatorial Pacific, including INGV, ISHII, LEVITUS, and notably GECCO. Weak cooling is also evident in most analyses just south of the equator in the Indian Ocean.
We next consider the impact of the three major eruptions: Mount Agung in 1963, El Chichón in 1982, and Mount Pinatubo in 1991. For each eruption the 4-yr heat content average following the eruption minus the 4-yr average up to and including the year of the eruption (Fig. 3), reveals a complex and varied pattern of change. Only five of the nine analyses extend back to 1960, prior to Mount Agung. Of these, the no-model analyses show only weak anomalies in its aftermath, while the data assimilation analyses show cooling north and south of the equator in the equatorial Pacific, with a suggestion of warming in the North Pacific and North Atlantic. The coupled model results interestingly show an antisymmetric response in the tropical Pacific with cooling north of the equator and warming south of the equator.
For El Chichón all nine analyses are available (Fig. 3) and the analyses as well as the coupled model are consistent in showing cooling on and/or just south of the equator in the tropical Pacific. Some cooling is also evident in many of the analyses in the southwestern tropical Indian Ocean. In the North Pacific many of the analyses show cooling in the Kuroshio–Oyashio Extension region. This cooling, which is not evident in the coupled model, may be an aspect of North Pacific climate variability unrelated to aerosol forcing. For Mount Pinatubo most analyses show general warming except in the western equatorial Pacific (Fig. 3).
We next consider the decade-by-decade heat content anomalies in each analysis (Fig. 4). To develop an understanding of the potential biases within the analyses we also examine corresponding decadal estimates of heat-content-analysis-minus-CTD observation differences (Fig. 5). These analysis-minus-CTD observation comparisons are limited to the Northern Hemisphere and to the three decades 1970–99 in order to ensure sufficient data coverage.
Examination of the heat content anomalies in Fig. 4 shows the North Atlantic to be warming rapidly in nearly all analyses similar to the behavior described in Levitus et al. (2000). In the 1960s most analyses are anomalously cool in the western tropics, subtropics, and midlatitudes by up to 6 × 108 J m−2 with anomalously warm temperatures in the eastern subtropics and subpolar regions. The main exception is GECCO, which is warm in the central and eastern subtropics, and along with INGV is anomalously cool throughout the tropics. In the 1970s and 1980s heat content anomalies are generally weak again, with the exception of GECCO, where the warm subtropics and cool midlatitude and tropical anomalies persist throughout the 1970s. However, the GECCO anomalies are superimposed on a mean state that is warmer than the CTD observations throughout much of the Atlantic so the analysis-minus-CTD observation differences are positive (Fig. 5).
By the 1990s the western subtropics to midlatitudes in most analyses have warmed (Fig. 4). In the tropics INGV has rather a strong warm anomaly (warmer than the CTDs, see Fig. 5), while GODAS has cool anomalies in the west. At subpolar latitudes many analyses show cool anomalies in the 1990s concentrated in the west and extending into the Labrador Sea, consistent with independent observations, as discussed in the introduction. Interestingly, half of the analyses show the cooling in the 1990s to be concentrated not at the latitude of the Labrador Sea, but farther south off Newfoundland. In most analyses the warming in the 1990s is concentrated in the upper 350 m, as noted by Levitus et al. (2000). Interestingly three analyses, CERFACS, INGV, and U.K. FOAM, also show significant warming extending to thermocline depths of 700 m, while two analyses, GECCO and GFDL, show the warming to be concentrated in the upper 200 m (Fig. 6).
Farther north in the Arctic Ocean all seven analyses that include a full Arctic Ocean (see Table 1) have cool heat content anomalies in the Eurasian sector in the 1960s and warm anomalies in the 1990s (not shown). The strongest warm anomalies, exceeding 6 × 108 J m−2, appear in GFDL and SODA. In other analyses the warming in the 1990s is primarily confined to the Greenland and Barents Seas where the warming has been previously documented, as discussed in the introduction.
The Caribbean Sea and Gulf of Mexico, which is coarsely resolved in most analyses, is cool in half the analyses in the 1960s (CERFACS, INGV, ISHII, and LEVITUS) and warm in the rest (Fig. 4). By the 1970s most analyses show cool anomalies, although examination of the analysis-minus-CTD observation differences show that these cool anomalies are still too warm compared to the observations (Fig. 5). By the 1990s the cool anomalies have been replaced by warm anomalies in all analyses except GECCO. However, in all analyses the warm anomalies in the Gulf of Mexico in the 1990s are too warm relative to the CTDs.
We next consider the decadal heat content of the other relatively well-sampled basin, the North Pacific (Fig. 5). While the North Atlantic is dominated by long-term warming, the North Pacific shows a modest warming trend, but with strong decadal heat content variability. Many analyses (with the exceptions being GFDL and GECCO) have cool heat content anomalies in the 1960s, particularly in the northern tropics, while all show neutral or warm heat content anomalies along the equator. Basin-average heat content in most analyses increases to a peak in 1980, and then cools in the Kuroshio–Oyashio Extension region until the late 1980s. In the 1990s the midlatitude ocean begins to warm again in most analyses while the western tropics cool and the eastern tropics warm.
Two analyses, GECCO and GFDL, behave somewhat differently in the North Pacific. GECCO is warm in the 1960s with the heat concentrated in the tropics, cool in the 1970s, and finally warm in the late 1980s and 1990s in the midlatitudes. GFDL is also very warm in the 1960s in the tropics, but by 1970 it resembles the other analyses. INGV and the no-model analyses have weaker anomalies, but patterns similar to CERFACS, SODA, and U.K. FOAM. Examination of the analysis-minus-CTD observation differences shows all analyses to be biased warm throughout the 1970s–90s, with the most pronounced bias evident in GECCO (Fig. 5).
The changes in the North Pacific heat content bear resemblance to the observed decadal variability of winter SST in the North Pacific (Pacific decadal oscillation; Mantua et al. 1997). Until the late 1970s SSTs were cool in the northwest basin and warm in the eastern tropics and along the west coast of North America. In the late 1970s until 1999 this pattern of SST reversed. Examination of the leading empirical orthogonal eigenfunction of heat content anomalies (computed for 20°–60°N, analogous to the definition of the PDO index) shows positive values along the eastern side of the basin and northwestern tropics in many of the analyses and negative values throughout the interior subtropics in most analyses (Fig. 7), which are similar in structure to the PDO. The explained variance of the leading mode ranges from a low of 7% for the highest-resolution analysis (SODA) to 17% and 19% for U.K. FOAM and INGV. The principal component time series of most analyses also show considerable similarity to the PDO index (Fig. 7). The exception is GECCO for which neither the first nor second empirical eigenfunction resemble the Pacific decadal oscillation pattern.
We next consider the poorly sampled South Atlantic and South Pacific Ocean. The no-model analyses, LEVITUS and ISHII, have only weak 3 × 108 J m−2 warm anomalies in the southwestern Pacific and Atlantic in the late 1960s and early 1970s, and then begin to cool (Fig. 8). Many of the sequential analyses (CERFACS, INGV, SODA, and U.K. FOAM) show patterns with similar spatial structures to the no-model analyses, but with larger 3–6 × 108 J m−2 amplitudes and with some evidence of warming in the southern midlatitude oceans by the 1990s. GECCO also shows anomalies that resemble the no-model analyses, but with even larger amplitudes, in excess of 15 × 108 J m−2, and with some other differences such as the presence of warm anomalies in the subtropical South Atlantic in the 1960s.
In the Indian Ocean the no-model analyses also show weak anomalies with cool anomalies in the southern subtropical west and in the Arabian Sea in the 1960s and warm anomalies south of 30°S in the 1990s (coinciding with the observed strong warming of SSTs discussed in the introduction; see Fig. 9). The sequential analyses show much stronger heat content anomalies, but with patterns generally resembling the no-model analyses. The most pronounced differences among the sequential analyses occur south of 30°S where CERFACS, SODA, and U.K. FOAM show strong warming, while GFDL shows no anomalies and GODAS shows warm and cool anomalies. Along the equator most of the sequential analyses show weak equatorial cooling, but SODA shows a very pronounced equatorial cool anomaly in the eastern half of the basin in the 1990s as a result of intense upwelling along the eastern equator in 1994 and 1997 (Grodsky et al. 2001). GECCO, alternatively, shows dramatic basin-scale warming north of 30°S except in the Arabian Sea, and cooling between 30° and 60°S.
Examination of the nine analyses shows that the global ocean 0/700 m has been warming at a rate of 0.76 × 108 J m−2 (10 yr)−1 (=0.24 W m−2) during 1960–2002, confirming the results of studies with individual analyses (e.g., Levitus et al. 2000, 2005; Carton et al. 2005; Ishii et al. 2006). Many uncertainties surround these individual estimates because of inadequacies in the historical observation network and the instruments, the techniques used to construct gridded analyses, and even the applicability of estimating a linear trend from the global data. By comparing the analyses to each other and to the historical observation set this paper we intended to shed light on these uncertainties and the presence of natural climate variability and to improve understanding of the utility of the analyses for decadal climate research.
The analyses we consider include two no-model analyses—ISHII and LEVITUS—in which a first guess of the climatological monthly cycle of upper-ocean temperature is updated based on historical ocean profile observations. We consider six analyses based on sequential data assimilation—CERFACS, GFDL, GODAS, INGV, SODA, and U.K. FOAM—in which an ocean general circulation model provides a first guess that is updated based on the ocean observations. Finally, we consider one analysis based on a 4DVAR algorithm—GECCO—in which a forward model and its adjoint are used to adjust initial conditions and surface forcing based on misfits to the ocean observations.
The analyses differ in choice of initial conditions, surface forcing, model physics, and resolution, and to a lesser extent which observations are used. Also, the analyses are affected by the different purposes for which they were created. Two of the analyses, GFDL and GODAS, were created to provide initial conditions for seasonal forecasts and thus have been designed to provide their best results in the upper levels of the tropics. Both of these have spurious anomalies in the first few years of the analysis period, and GODAS has very large anomalies in the Southern Ocean. Many of the analyses have 1° × 1° horizontal resolution and thus do not resolve the ocean mesoscale. However, several have enhanced tropical resolution to resolve the intense tropical current systems and one, SODA, has global eddy-permitting resolution in order to resolve features of the midlatitude and subpolar circulation. Seven of nine analyses (CERFACS, SODA, GFDL, INGV, ISHII, LEVITUS, and U.K. FOAM) include the Arctic Ocean. Five of nine analyses (CERFACS, GECCO, GODAS, INGV, and U.K. FOAM) assimilate satellite altimeter sea level observations as well as the historical profile dataset. Of these, GODAS and GECCO show large changes in the late 1980s to early 1990s, which may be related to the introduction of altimeter sea level data.
Despite these differences, the analyses do have many common features. Most analyses show the following two periods of rapid heating: first in the early 1970s and again beginning in the late 1980s. One possible explanation for these periods of rapid warming is that they occur following cool periods associated with the major tropical volcanic eruptions of Mount Agung in 1963 and El Chichón in 1982 (Ramaswamy et al. 2001; Church et al. 2005). Indeed, coupled modeling results reported in Delworth et al. (2005) show reductions in ocean heat content of 1–2 × 108 J m−2 lasting for a few years resulting from the impact of volcanic aerosols (less in response to Mount Agung). The spatial pattern of the heat content anomalies produced by their coupled model does indicate the possibility of a cooling of ocean heat content in the equatorial Pacific. However, the aerosol effect is unlikely to explain the nearly decadal pattern of heat content anomalies that appear in the global time series. We note, for example, that the analysis heat content estimates do not seem to reflect an impact from the Mount Pinatubo eruption of 1991.
A second possible explanation for these periods of rapid warming is that it is the result of combining the “natural” variability of the climate system operating fairly independently in different basins. Basin-average heat content time series for the North Atlantic (Fig. 10a), for example, shows a rapid, fairly constant increase at a rate of 1.7 × 108 J m−2 (10 yr)−1 [Levitus et al. (2005) finds a similar nearly linear trend of 2.2 × 108 J m−2 (10 yr)−1 for the full Atlantic basin]. The horizontal and vertical structure of the anomalies are similar for the no-model and sequential analyses, with most of the heating occurring in the subtropics and midlatitudes and weak cooling in and south of the Labrador Sea. Much of this additional heat is confined to the upper 400 m. Two regions of particular interest, the Eurasian sector of the Arctic Ocean and the Gulf of Mexico, likewise show significant warming in the 1990s. In the North Pacific most of the basin-average time series (Fig. 10b) show strong decadal variability, that is, warm in the mid-1970s to early 1980s and again beginning in the early 1990s. This North Pacific variability is very similar to that examined in previous studies, such as Mantua et al. (1997) and Zhang et al. (1997).
Heat content anomalies in the Southern Hemisphere in many of the analyses are large and thus are very important contributors to the global average, but are not well documented previously. Although the analyses show many differences, they also show some common features. In contrast to their Northern Hemisphere counterparts, Southern Hemisphere anomalies extend zonally across ocean basins. The 1960s are characterized by warm anomalies in the 0°–10°S and 30°–40°S latitude bands, with cool anomalies in between. By the 1990s much of the Southern Hemisphere has warmed in most analyses, with the most pronounced warming occurring in the Atlantic and Indian sectors. Among the analyses the no-model analyses (ISHII and LEVITUS) and INGV show only weak anomalies. GECCO and GODAS have more vigorous anomalies, with GECCO showing strong warming and GODAS showing strong cooling in recent decades.
A third possible explanation for the periods of rapid global warming is that they are artificial, resulting from changes in geographic coverage and the introduction of an uncorrected warm bias resulting from the introduction of the XBTs in the 1970s. The warm bias, which varies by batch and type of XBT, would affect all analyses, but not the coupled models, which do not show such periods of rapid warming (Gregory et al. 2004; Delworth et al. 2005). We explore this explanation by comparing the analyses to the unbiased CTD observation dataset during three decades of 1970–99 in the relatively well-sampled North Atlantic and North Pacific. This analysis-minus-CTD observation comparison, shown geographically in Fig. 5 and with time in Fig. 11, reveals a warm bias that is most pronounced in the 1970s and appears to diminish by the 1990s. The change in bias from the 1970s through the 1990s is more pronounced in the North Pacific and likely contributes to the warmth of the 1970s relative to the multidecadal trend, but does not explain the rapid rise in global heat content beginning in the 1990s and continuing to the present (see Fig. 11). These results confirm the results of Gouretski (2008, personal communication) and Wijffels et al. (2008), who directly compare XBT and CTD observations.
A comparison of the results of the current study with a preliminary attempt to compare ocean analyses by Chepurin and Carton (1999) highlights the tremendous progress that has been made in reconstructing historical variations in ocean heat content. These results indicate that the sequential analyses have scientifically interesting subseasonal variability in both hemispheres despite problems resulting from measurement bias (a problem that afflicts all current ocean analyses). Still, the current study does show continuing differences in representation of climate anomalies, particularly in regions of poor historical observation coverage like the Southern Hemisphere. Addressing the causes of these differences ultimately may require reconsideration of both surface meteorology and assimilation strategies in this region. Further improvements to the analysis of ocean temperature likely will also require consideration of salinity because of its contribution to the stability of the water column and thus future studies will need to consider both variables.
The authors want to express their appreciation to the lead authors of the analyses used in this study: David Behringer, Mike Bell, Alessio Bellucci, Michael Davey, Armin Köhl, Masayoshi Ishii, Tony Lee, Sydney Levitus, Chaojiao Sun, and their colleagues. This work was motivated by a series of workshops organized by Professor Detlef Stammer as part of the CLIVAR Global Synthesis and Observations Panel. Many thanks to Dr. Gennady Chepurin for processing the CTD data and to Dr. Semyon Grodsky for illuminating discussions. Financial support for this research has been provided by the National Science Foundation (OCE0351319). JAC gratefully acknowledges University of Maryland and the NOAA/Geophysical Fluid Dynamics Laboratory for providing sabbatical time to pursue this project.
Corresponding author address: Dr. James A. Carton, Department of Atmospheric and Oceanic Science, University of Maryland, College Park, College Park, MD 20742. Email: email@example.com