Historical bathythermograph datasets are known to be biased, and there have been several efforts to model this bias. Three different correction models of temperature bias in the historical bathythermograph dataset are compared here: the steady model of Hanawa et al. and the time-dependent models of Levitus et al. and Wijffels et al. The impact of these different models is examined in the context of global analysis experiments using the Simple Ocean Data Assimilation system. The results show that the two time-dependent bias models significantly reduce warm bias in global heat content, notably in the 10 years starting in the early 1970s and again in the early 1990s. Overall, the Levitus et al. model has its greatest impact near the surface and the Wijffels et al. model has its greatest impact at subtropical thermocline depths. Examination of the vertical structure of temperature error shows that at thermocline depths the Wijffels et al. model overcompensates, leading to a slight cool bias, while at shallow levels the same model causes a slight warm bias in the central and eastern subtropics and at thermocline depths on the equator in the Pacific Ocean as a result of reduced vertical entrainment. The results also show that the bias-correction models may alter the representation of interannual variability. During the 1997/98 El Niño and the subsequent La Niña the Levitus et al. model, which has its main impact at shallow depths, reduces the 50-m temperature anomalies in the eastern equatorial Pacific by 10%–20% and strengthens the zonal currents by up to 50%. The Wijffels et al. correction, which has its main impact at deeper levels, has much less effect on the oceanic expression of ENSO.
Studies of the 50-yr record of global ocean temperature observations show decadal variations, with pronounced spatial patterns, superimposed on a gradual warming trend (Levitus et al. 2000; Carton and Santorelli 2008). These decadal variations, including the apparent sudden warming and subsequent cooling of the global ocean during the 15-yr period of 1971–85, have been explained as a time-dependent bias in the historical bathythermograph records (Gouretski and Koltermann 2007; Carton and Santorelli 2008; Wijffels et al. 2008, hereafter W08; Domingues et al. 2008; Ishii and Kimoto 2009; Levitus et al. 2009, hereafter L09). This bias, which has been noted in many historical instrument comparisons, can exceed 0.3°C in the upper 150 m (Gouretski and Koltermann 2007; L09).
The historical subsurface ocean temperature dataset begins in the early twentieth century and consists of vertical profiles from a variety of instruments (Boyer et al. 2006). From the 1950s to the early 1970s the majority of the observations were collected with a device called the mechanical bathythermograph (MBT), whereas between the 1970s and the early 2000s (our period of interest for this study is 1967–2002) the majority of the profiles were collected with versions of the expendable bathythermograph (XBT). Neither instrument was designed to monitor climate, and their a posteriori application to climate research has revealed a pattern of bias.
One of the first quantifications of an XBT bias in the refereed literature dates from a comparison study by Flierl and Robinson (1977), who used observations made with the Sippican Corporation T-7 and contemporaneous unbiased conductivity–temperature–depth (CTD) measurements in the Sargasso Sea. Flierl and Robinson found that XBT measurements overestimate the depth of isotherms in the upper 450 m while measurements below 450 m underestimate the depth of isotherms. Both differences appear to be a result of an error in the equation used to relate time measurements with depth using the instrument fall rate. Several subsequent studies confirm the findings of Flierl and Robinson, including the widely cited study of Hanawa et al. (1995) that is based on observations collected during a limited 8-yr period. However, Hanawa et al. and the other bias studies provide corrections only for limited batches of instruments and so the corrections are not representative of the global XBT dataset. In addition, the method by Hanawa et al. was designed for the estimation of a depth bias only, with no attempt to correct for the pure thermal bias. Recent studies by Gouretski and Koltermann (2007), W08, Ishii and Kimoto (2009), and L09 attempt to examine patterns of bias in the whole bathythermograph dataset and to produce models of that bias.
The first application of these new bias models has been to correct the historical time series of global ocean heat content, where they have the effect of reducing decadal variability (Domingues et al. 2008; W08; L09; Ishii and Kimoto 2009). Although the various models all produce a multidecadal increase in global averaged heat content, the vertical and latitudinal structure of the temperature correction is strikingly different. The zonally and temporally (1967–2002) averaged difference in observed temperature for L09 minus the World Ocean Database 2005 (WOD05; Boyer et al. 2006) and for W08 minus WOD05 are shown in Fig. 1. Much of the correction proposed by L09 [as well as by the earlier study of Gouretski and Koltermann (2007)], based on comparisons of collocated instrument records, occurs near the surface. This correction cools the XBTs by as much as 0.1°C in the time and spatial average, with much of the cooling (>0.2°C) occurring in the 1970s. In contrast, W08 and Ishii and Kimoto (not shown) include a theoretical model of the vertical distribution of the correction based on the assumption that the cause of the bias is unrepresented errors in the instrument fall-rate formulas. As a result of this assumption, much of the correction as proposed by either study (a cooling that exceeds 0.15°C on average) occurs at thermocline depths in the tropics and subtropics.
These differences in the spatial structure of the bias models can have a very different impact on the geographic and time evolution of estimates of temperature that will spread geographically and to different variables when included in a data assimilation system. To explore these secondary effects and to allow a more complete comparison of instrument records, we examine identical twin experiments. The control experiment uses the Simple Ocean Data Assimilation (SODA) analysis procedure of Carton and Giese (2008) and the time-independent depth bias model of Hanawa et al. (1995), whereas the two twin experiments use the same analysis procedure, forcing, and observation sets but employ two alternative time-dependent instrument bias models (those of W08 and L09) applied to the same dataset. The method is described in section 2, with results from the experiments, including an error analysis, in section 3 and discussion in section 4.
2. Data and methods
The configuration of the SODA system used in this study is similar to that described by Carton and Giese (2008) except that no SST data are used. The ocean model is based on the Los Alamos National Laboratory Parallel Ocean Program, version 2.0.1, numerics (see online at http://climate.lanl.gov/Models/POP/UsersGuide.pdf). The horizontal resolution, averaged globally, is 0.40° in longitude × 0.25° in latitude with 40 vertical levels (with 10-m resolution near the surface) and a fully resolved Arctic Ocean. Vertical mixing is based on the Large et al. (1994) K-profile parameterization, and horizontal mixing is biharmonic. Rivers are included with climatological seasonal discharge. There is no explicit sea ice model, although surface heat flux is modified to simulate its climatological seasonal impact. The model includes an explicit free surface.
Initial conditions and surface forcing fields are identical for all experiments. The 40-yr European Centre for Medium-Range Weather Forecasts re-analysis (ERA-40) daily surface winds (Uppala et al. 2005) are used from 1 January 1958 through 31 December 2001, and Quick Scatterometer (QuikSCAT) satellite winds are used for 2002. For the period 1979–2002 the surface freshwater flux is provided by the Global Precipitation Climatology Project monthly satellite–gauge merged product (Adler et al. 2003) combined with evaporation obtained from the same bulk formula used to calculate latent heat loss. For years prior to 1979, climatological seasonal precipitation is used. The other terms in the surface heat flux boundary conditions are determined using bulk formulas with parameters obtained from the National Centers for Environmental Prediction–National Center for Atmospheric Research analysis but are strongly influenced by temperature updating within the mixed layer. Sea surface salinity includes a relaxation to the climatological seasonal sea surface salinity of Boyer et al. (2002).
The assimilation is carried out sequentially using a 10-day update cycle with model error covariances determined from a simulation that does not include assimilation. The error covariances evolve in time as a function of the local velocity field and mixed layer depth. Updating is done incrementally following Bloom et al. (1996) to suppress excitation of spurious variability. Output variables are averaged monthly and then mapped onto a uniform global 0.5° × 0.5° horizontal grid using the horizontal grid spherical coordinate remapping and interpolation package of Jones (1999). All experiments begin 1 January 1967 after a 9-yr spinup of the model beginning with climatological initial conditions.
Updating in all three experiments described in Table 1 relies on the set of temperature and temperature–salinity profiles from MBTs, XBTs, ocean station data (OSD), and CTD data contained in WOD05 (Boyer et al. 2006). The spatial coverage of this dataset, shown binned into 1° × 1° × 5-day bins and limited to our analysis time period, is concentrated in the Northern Hemisphere, with the heaviest coverage along routes by commercial ships of opportunity (Fig. 2, top panel). Most of the profiles along these routes were collected using bathythermograph instruments. The distribution of the combined CTD plus OSD observation set is more spatially homogeneous but comprises fewer than 30% of the total number of binned observations (Fig. 2, middle panel). Most of these are OSD observations, particularly for the period before 1990 (Fig. 3). After 1995 the situation reverses and CTDs dominate the combined dataset (our period of interest is limited to 2002 to exclude the influence of the more recent Argo data).
The depths of the XBT profiles used in the control experiment, SODA2.1.0, have been corrected by application of the steady Hanawa et al. (1995) bias model. No correction was made to the MBT data in the control run. The second experiment, SODA2.1.2, uses the bias-correction model of L09 for the MBTs and XBTs. The third, SODA2.1.4, uses the bias-correction model of W08 for the XBTs as shown in their Table 1. Since W08 has no MBT correction, we retain the L09 MBT correction. We have chosen to use the W08 correction as applied to the same WOD05 dataset used in our other experiments rather than the dataset described in W08 because a preliminary comparison of the two datasets shows that they differ somewhat (<3%) and we want to exclude the effects of differing data selection from our comparison.
Two types of comparisons are presented. In the first we present differences between the twin experiments SODA2.1.2 or SODA2.1.4 and the control experiment (SODA2.1.0). These differences show the impact of using time-dependent bathythermograph bias models. In the second comparison we present differences between the gridded temperature fields from the experiments and the set of CTD/OSD station temperature casts. Although the CTD/OSD data are used in the reanalysis, the data are the same for both experiments, so that the only variable is the correction to the XBT/MBT data. This difference procedure is similar to the approach used by L09, who compare XBT profiles with collocated CTD data to develop their bias-correction algorithm. Using the reanalysis output allows us to directly compare our results with all available CTD/OSD observations, which doubles the number of comparison pairs to ∼2000 per month relative to L09. These differences are then mapped onto a uniform 1° × 1° monthly grid at each model level using a kriging procedure to reduce the influence of the inhomogeneous spatial distribution of the observations and their differences. This monthly mapped error is the basic dataset we use for our error analysis.
The control experiment SODA2.1.0 global heat content (0–700 m) shows decadal variability similar to that of previous studies such as Carton and Santorelli (2008), including the dramatic warming by 1023 J in the early 1970s and the early 1980s (Fig. 4a). Here heat content is calculated by integrating temperature times density and heat capacity in the upper 700 m. The anomaly is found by removing the time mean. However, the fact that the corresponding estimate of heat content error shows the same multidecadal variations as are seen in heat content (Fig. 4b) suggests that the decadal variability is erroneous. Since decadal variability appears in studies using many different mapping techniques, it is indicative that the error is in the hydrographic datasets rather than in the mapping. We confirm this result by repeating the experiment using the L09 and W08 time-dependent bias models to correct the hydrographic datasets (SODA2.1.2 and SODA2.1.4).
Global heat content from both of the twin experiments SODA2.1.2 and SODA2.1.4 shows much-reduced decadal variability and reduced error, confirming the results reported by L09 (the latter is included in Fig. 4 for comparison). Comparison of the results using the alternative bias models shows that SODA2.1.4, which uses W08, has somewhat lower error than that in SODA2.1.2, which uses L09. However, the question of which bias model is superior is more complicated, as discussed below. Of interest, the heat content error in SODA2.1.4 becomes negative after the mid-1990s, indicating that the W08 procedure slightly overcorrects the warm bias in the bathythermographs after the mid-1990s, perhaps because of the increasing numbers of deep XBTs. Although the Domingues et al. (2008) data (labeled D08 in Fig. 4a) show a smaller rate of increase than do the L09 data during the 1990s, the experiment using W08 (SODA 2.1.2) shows a much stronger cooling during the 1990s.
Next we consider the heat content trends in individual basins: the North Pacific and North Atlantic Ocean basins, where the data coverage is much greater than in other basins (Figs. 5a,b), and the tropical Pacific, because of its link to ENSO. In the North Pacific the control experiment SODA2.1.0 heat content has strong decadal variability with essentially no trend [also evident in previous studies such as Levitus et al. (2005) and Carton and Santorelli (2008)]. When either time-dependent bias model is applied, however, decadal heat content variability gets weaker (notably an elimination of an apparent warm anomaly in the mid-1970s to early 1980s and a cool anomaly in the late 1980s) and the basin-average time series develops a pronounced warming trend (0.08 × 1022 J yr−1 for SODA2.1.2 and 0.05 × 1022 J yr−1 for SODA2.1.4).
In the North Atlantic (Fig. 5b) all three experiments show similar gradual warming trends. The slight warm anomaly in the mid-1970s to early 1980s, previously noted in the North Pacific, is also evident in this ocean and is also reduced when using either of the two bias models.
The impact of the bias models on the vertical structure of the analyses in the North Pacific is shown in Fig. 6. The W08 correction used in SODA2.1.4 has a bigger impact than the L09 correction in SODA2.1.2 and has its greatest impact in the thermocline and below it. Therefore, vertical structure of the time-mean SODA2.1.4–SODA2.1.0 difference averaged over 20°–30°N shows that the W08 model cools the thermocline by up to 0.3°C in the central and eastern parts of the basin. The W08 model also results in a counterintuitive warming (<0.15°C) of the upper 100 m due to reduced entrainment cooling from below the mixed layer, an effect that is even more pronounced in the southern subtropics (Fig. 7). The impact of the W08 temperature bias correction on the mean velocities in the subtropics appears to be small, reducing the strength of the Kuroshio by less than 1 cm s−1 as a result of a reduction of the zonal pressure gradient. The L09 correction has only a weak impact on time-mean temperature and an even weaker impact on velocity in the subtropics.
In the equatorial Pacific, the W08 correction in SODA2.1.4 acts, on average, to cool temperatures over much of the upper 400 m but again acts to warm the thermocline temperature by 0.3°C (Fig. 8). The counterintuitive warming is indirectly due to a strengthening of eastward currents caused by an increase in eastward temperature advection relative to the control experiment. In contrast, the L09 correction of SODA2.1.2, whose main impact is near the surface, has a negligible impact on mean temperature, except close to the eastern boundary where the thermocline shallows.
The error in the time-average temperature fields, evaluated by comparison with the CTD/OSD dataset, is shown in Fig. 9. The control experiment thermocline is too warm by 0.05°–0.1°C. The L09 model of SODA2.1.2 reduces the temperature error in the upper ocean by a factor of 2 without increasing the error below the thermocline. In contrast, the W08 model of SODA2.1.4 leaves the thermocline water somewhat too warm and the subthermocline water slightly too cold.
Although the impact of the L09 bias model in SODA2.1.2 on time-averaged currents is modest, the instantaneous impact can be considerably larger. For example, during the 1997/98 El Niño event 50-m temperature changed by 1.0°C or more and currents changed by 20 cm s−1 (Fig. 10). These changes act to reduce temperature anomalies associated with both warm El Niños and cool La Niñas by 10%–20%. The impact on zonal velocity is even larger, accelerating currents in the eastern basin in the late stages of the 1997/98 El Niño and extending into the 1998 La Niña by up to 50%. A comparison of 5-day-averaged 50-m zonal currents at 140°W during the 1997/98 El Niño shows that the rms difference between the observed zonal velocity from the ADCP and SODA 2.1.2 is 30.6 cm s−1, whereas the rms difference for SODA 2.1.0 is 40.7 cm s−1 and the rms difference for SODA 2.1.4 is 48.4 cm s−1.
Studies of the changing climate of the ocean are limited by the fact that the ocean observing system was not designed for climate monitoring. One consequence is that until recently relatively little attention was paid to the presence of a time-dependent warm bias in the bathythermograph dataset. However, motivated by the need for more accurate estimates of global ocean heat content, Gouretski and Koltermann (2007), W08, Ishii and Kimoto (2009), and L09 introduced models to correct the bias in the bathythermograph dataset a posteriori. These bias studies rely on a combination of theoretical understanding of the causes of the bias (e.g., assuming that it results from changes in the fall-rate equation) and direct comparison of the bathythermograph profiles with collocated unbiased instruments such as CTDs.
The resulting models have very different vertical structures, particularly below 400 m, even though their vertical integrals are similar. For example, L09 introduce most of their correction near the surface based on comparisons of bathythermographs and independent CTD/OSD observations. Although the L09 adjustment may be correcting mostly for a fall-rate error, the underlying assumption is that this error is too complex to represent as a change in the fall-rate equation. In contrast, W08 and Ishii and Kimoto (2009) assume that the error is the result of instrument fall rate and use instrument comparisons to derive an altered fall-rate equation. In this study we explore the consequences of these two alternative types of bias models in comparison with the widely used time-independent correction of Hanawa et al. (1995) to examine the impact of correcting time-dependent bias on the reconstructed ocean state produced by a full data assimilation system.
Comparison of data assimilation experiments using the L09 and W08 time-dependent bias models with a control experiment using the older time-independent bias model of Hanawa et al. (1995) shows that the new models reduce decadal variability of both the global ocean heat content (and in particular the warm anomaly from the early 1970s to early 1980s) and the decadal variability in individual basins such as the North Pacific. Both models also reduce the temperature error of the analysis in comparison with independent OSD/CTD observations. The W08 bias model tends to undercompensate for temperature error above the thermocline and overcompensate below the thermocline, however, resulting in a slight negative temperature bias in the thermocline. The L09 bias model also undercompensates for temperature bias in the near surface but does not overcompensate in the thermocline. The W08 bias model cools the observations yet results in warming near the surface in the subtropics because of reduced entrainment and at thermocline depths along the equator because of enhanced warm advection. The L09 bias model has little impact on the time-averaged conditions in the tropical Pacific but has a significant impact on the ENSO response in the eastern basin, reducing 50-m temperature anomalies by 10%–20% and causing changes in 50-m currents by up to 20 cm s−1 (Fig. 10). The latter changes improve the comparison of analysis currents with those directly observed from the Tropical Ocean and Global Atmosphere/Tropical Atmosphere–Ocean current data at 140°W in the Pacific Ocean.
Based on the results of this study, we conclude that both types of time-dependent bias models represented by L09 and W08 reduce the estimated heat content error relative to the widely used steady-state bias model of Hanawa et al. (1995). However, experiments with the SODA data assimilation system show that the bias in the XBT records is better represented by L09, which uses a statistical temperature correction that makes no assumptions about the cause of the bias, rather than by W08, which uses a statistical depth correction based solely on the assumption of fall-rate error as the cause of the XBT bias. As new correction schemes become available, such as the one proposed by Gouretski and Reseghetti (2010), we will continue to use SODA as an evaluation framework. The issue of bias correction in the observational records continues to be a topic of great importance for climate studies and requires continued study with independent observations to reduce the biases.
Financial support for this research has been provided to BSG and HFS by NOAA (Grant NA06OAR4310146) and NSF (Grant OCE-0351804). Support for JAC and GAC is provided by the National Science Foundation (Grant OCE-0351319). The authors acknowledge the Texas A&M Supercomputing Facility (http://sc.tamu.edu/) and, in particular, Ping Luo for providing computing resources and advice that were essential in the research reported in this paper.
Corresponding author address: Benjamin S. Giese, Dept. of Oceanography, Texas A&M University, College Station, TX 77843–3146. Email: firstname.lastname@example.org