A new sea surface temperature (SST) analysis on a centennial time scale is presented. In this analysis, a daily SST field is constructed as a sum of a trend, interannual variations, and daily changes, using in situ SST and sea ice concentration observations. All SST values are accompanied with theory-based analysis errors as a measure of reliability. An improved equation is introduced to represent the ice–SST relationship, which is used to produce SST data from observed sea ice concentrations. Prior to the analysis, biases of individual SST measurement types are estimated for a homogenized long-term time series of global mean SST. Because metadata necessary for the bias correction are unavailable for many historical observational reports, the biases are determined so as to ensure consistency among existing SST and nighttime air temperature observations. The global mean SSTs with bias-corrected observations are in agreement with those of a previously published study, which adopted a different approach. Satellite observations are newly introduced for the purpose of reconstruction of SST variability over data-sparse regions. Moreover, uncertainty in areal means of the present and previous SST analyses is investigated using the theoretical analysis errors and estimated sampling errors. The result confirms the advantages of the present analysis, and it is helpful in understanding the reliability of SST for a specific area and time period.
Sea surface temperature (SST) is a basic variable in climate studies. At present, more than 150-yr-long SST time series are available thanks to international cooperation to form a marine meteorological data archive: the International Comprehensive Ocean–Atmosphere Data Set (ICOADS; Woodruff et al. 2011). This database contributes to the understanding of past climate changes over the globe. In addition, climate modeling studies benefit from this activity because it provides a boundary condition for atmospheric models (Gates 1992), atmospheric reanalyses (Kalnay et al. 1996; Uppala et al. 2005; Onogi et al. 2007), and validation of climate model outputs on historical time scales (e.g., Nozawa et al. 2005).
In our previous study, a system was developed that produces centennial-scale objective analyses of marine meteorological variables [Centennial In Situ Observation-Based Estimates of the Variability of SST and Marine Meteorological Variables (COBE); Ishii et al. 2005, hereafter ISSM05]. This system is helpful in understanding not only long-term changes in the variables but also the quality of historical observations. Simultaneously, we evaluated the impact of the Kobe collection (Manabe 1999) on the reduction of estimation errors in SST and sea level pressure.
Removing biases from historical observations is a delicate issue. All marine meteorological variables may suffer from some biases, attributable to changes in instrumentation and observation methods. The most critical aspect is that they can yield unrealistic variations in time series for the global oceans and in some particular areas, and they are sometimes comparable to climatic signals. For SST, cold biases in uninsulated-bucket observations must be eliminated in order to obtain global mean SST time series of homogeneous quality (Folland and Parker 1995, hereafter FP95; Smith and Reynolds 2002). A similar example occurs in oceanographic subsurface temperature observations by expendable bathythermograph (Gouretski and Koltermann 2007; Wijffels et al. 2008; Ishii and Kimoto 2009; Levitus et al. 2009). The biases in this case cause a large change in the time series of global ocean heat content from the 1970s to the 1980s.
In this report, a new historical SST analysis is presented and compared with the previous version and other widely used SST analyses. The analysis method has been upgraded. Also, we attempt to correct historical SST biases. The quality of the analysis depends closely on the spatiotemporal distribution of the observations used. Sampling errors in areal mean SSTs are discussed by conducting a virtual analysis with observations that have daily distributions resembling those of the actual observations in the past.
In the new analysis, the latest version of ICOADS (release 2.5; Woodruff et al. 2011) is used. As in the previous analysis, a set of extra SST observations for the area around Japan before World War II is merged with ICOADS. This dataset is archived by the Japanese Fishery Agency (Tomosada 1982). An operational database compiled by the Japan Meteorological Agency is also used for the period after 2008. Figure 1 shows time series of the total number of SST observations in each year. The number of SST observations increases with time, and 100–1000 times more observations are available after 1950 than in the 1850s. Bucket observations compose the majority of observations before World War II and are gradually replaced by engine room intake (ERI) observations during the war. Bucket measurements appear to have been widely used even for the postwar period, although the method of measurement is not recorded for more than half of all SST observations. Recent SST measurements are conducted primarily by drifting buoy, ERI, and hull-contact sensor with the number of bucket and unknown-type observations decreasing since the 1980s.
Figure 2 depicts the annual coverage of observational SST data over the global oceans used in ISSM05 and the present study. The coverage based on all observations has increased over time from 20% in the 1850s to 90% in the 2000s. These time series are based on the metadata recorded in the datasets. Additional metadata from 1955 to 1994 compiled by Kent et al. (2007) are used to compensate for a lack of metadata in ICOADS. The source of these metadata is the World Meteorological Organization (WMO) publication 47. This improves the coverage of ERI during this period in the new database, although data of an unknown type are still dominant. The observational database offers better temporal coverage than was available previously for the 1880s for unknown-type measurements and the 1960s for bucket measurements. The data for these periods are from the U.S. Marine Meteorological Journals and the U.K. Marine Data Bank, respectively (Woodruff et al. 2011).
To construct unobserved variability in data-sparse regions, satellite observations are incorporated into the present objective analysis scheme. However, the satellite observations are used only for constructing empirical orthogonal functions (EOFs) that represent interannual-to-decadal SST variations; they are not used in the final product. This is because the use of satellite data for the whole analysis makes the gridwise variability of the SST analysis larger by 10%–20%, compared with that without the data. The EOFs are used in a new analysis scheme as described in section 6. For the computation of EOFs, we conducted an SST analysis using satellite and in situ observations as described below. The satellite observations used here are derived from the Advanced Very High Resolution Radiometer (AVHRR) Pathfinder version 5 of the Physical Oceanography Distributed Active Archive Center (PO.DAAC) of the National Aeronautics and Space Administration (NASA). First, large-scale satellite cold biases (Reynolds et al. 2002) are estimated on a daily global 2° × 2° grid from observed differences between matched-up satellite and in situ data from 1981 to 2010. Daytime and nighttime biases of satellite are computed separately. As done in the daily analysis of ISSM05, daily changes in the biases are computed by an optimal interpolation (OI) scheme. The spatial and temporal decorrelation scales are set to 1500 km and 60 days, respectively, based on statistics of the matched-up differences. The same background errors and observational errors as used in ISSM05 are prescribed for the OI scheme except weights for observational errors of individual observational types, which are redefined in this study as described in section 3. In the satellite bias analysis, an error variance for the observed satellite minus in situ SST difference is the sum of these two observational errors. Using the analyzed biases, the satellite observations are adjusted to the in situ data. Second, we conduct a daily SST analysis on a global 1° × 1° grid with the same OI scheme as in ISSM05, using daily differences between satellite and in situ data. Here, both of the data are bias corrected. This SST analysis is used for the EOF computation mentioned above. The satellite biases are excluded successfully, because the global mean SSTs based on satellite analysis agree with an in situ–only SST analysis within the analysis error. Biases in in situ observations are also removed before the above analysis. The details of the bias corrections of in situ data are presented in section 4.
The SST analysis is constrained by sea ice so that SSTs in the north and south polar regions are represented by a proxy of SST estimated from sea ice concentration (SIC). This ensures that SST fields are represented continuously between regions with and without sea ice. Satellite observations used here are from the Nimbus-5 Scanning Multichannel Microwave Radiometer (SMMR) from October 1978 to July 1987 and the Special Sensor Microwave Imager (SSM/I) and the Special Sensor Microwave Imager/Sounder (SSMIS; Armstrong et al. 2012) from August 1987 to December 2010. Satellite sea ice analyses were conducted using the NASA team algorithm (Cavalieri et al. 1984, 1991) and a bootstrap method (Comiso et al. 1997) on a daily and 1° × 1° grid. The latter is not applicable to SMMR, because it requires an additional channel that is unavailable in the SMMR data. The NASA team algorithm slightly underestimates SIC particularly in summer; thus, we opted for the bootstrap method. For the period covered by SMMR, sea ice fields calculated according to the NASA team algorithm are modified by adding climatological averages of SIC differences between the NASA team and bootstrap methods with SSM/I and SSMIS observations from 1988 to 2000. Sea ice over the Arctic Ocean before the satellite era from January 1870 to November 1978 is provided by Walsh and Chapman (2001). The monthly climatology of the Walsh sea ice for 1968 to 1977 is adjusted to that of a detrended satellite analysis for 1979 to 1988. In this dataset, SIC is unavailable for the most part of the Sea of Okhotsk. Therefore, we patched this region using the climatology of the detrended satellite sea ice mentioned above. Similarly, the Antarctic sea ice prior to November 1978 is filled using the climatology. Time series of the Arctic and Antarctic sea ice extent are shown in Fig. 3. Decreasing and weak increasing trends appear in the Arctic and Antarctic time series, respectively. The time series of the present study agrees well with Hadley Centre Sea Ice and Sea Surface Temperature dataset (HadISST) during the satellite era, although the amplitude is slightly smaller than that of HadISST particularly in recent years. Although the Antarctic sea ice before 1978 is the mean from 1979 to 1988, the sea ice extension of the Antarctic climatology is larger than the mean of the extension for the same period as that of the mean. This is due to wide-spreading low SIC of the climatology. The HadISST time series shows a large discontinuity of the Antarctic sea ice extent before 1979, as a result of the use of additional various and heterogeneous sea ice observations (Rayner et al. 2003). Melting of a large portion of the Arctic sea ice is observed in summer of 2007, and variations in sea ice extent have increased in subsequent years. The maximum sea ice extent is bounded by the geographical shape of the coasts surrounding the Arctic Ocean. This limitation is observed in cases when the summer sea ice extends to the coast, as results from a dynamical model (Goosse et al. 2009) suggests. After 2007, the sea ice edges are away from land and the limitation may become ineffective.
3. Quality control
The quality-control procedure used in this study is the same as that described by Ishii et al. (2003) and ISSM05. Observations with plausible dates and locations are compared with climatology and nearby observations to assess whether their values are realistic. Of observations passing the above checks, those in each 1° × 1° box and in the same calendar day are merged into one observation with their mean value at their mean location and time. Here, considering correlation in errors between observations reported by the same ship, such observations are first averaged ship by ship before merging. The observational error variance of the merged observation, which is prescribed for the objective analysis scheme, is reduced by a factor proportional to the number of ships and data (Ishii et al. 2003). Observational standard deviations vary in space and they are identical to those of interannual SST variation in this study, as introduced later together with the climatology. Weights according to instrument types are also considered, as listed in Table 1. The observational errors shown in the table are expressed as ratios to the background error. A blacklist of ship and buoy is introduced in order to remove biased observations that cause analyzed anomalies being abnormally enhanced along the tracks of ship and buoy, and a moving-speed check is also incorporated. Details of these two checks are described in ISSM05.
4. Bias corrections
a. Bucket biases
Correcting biases in historical SST observations is crucial for the reproduction of secular SST changes. However, a lack of metadata complicates this problem. Uninsulated-bucket observations have sizable negative biases to be removed, compared with ERI and hull-contact sensor observations (FP95). Such bucket observations are dominant for the period before World War II. The number of bucket data dramatically reduces during the war, but increased again for the postwar period (Fig. 1). Uninsulated buckets seem to have been used widely even after the war, decreasing gradually to 5% of all SST observations by around 1970 according to the literature and metadata, although the exact time is uncertain (Kennedy et al. 2011). FP95 proposed seasonally varying bucket biases from 1856 to 1941, which are attributable to changes in ship speed and the fraction of measurements originating from uninsulated buckets. In this analysis, the FP95 correction adjusted to a new climatology of this study is applied to all SST observations except ERI observations prior to 1941. The adjusted correction is also referred to as FP95, and how to adjust it is described later (section 4c).
The FP95 correction of 1941 is used repeatedly after 1942 to remove SST biases in uninsulated-bucket observations. This is necessary to avoid a large discontinuity in global mean SST in the mid-twentieth century (Thompson et al. 2008). In fact, there are few metadata in existence (Figs. 1, 2). Kennedy et al. (2011) estimated the fractions of different measurement methods by using information about the country that equipped the ship. By contrast, in this study, time-varying ratios of uninsulated-bucket observations to all unknown-type observations are inferred by comparing SSTs between unknown and known types, as presented below.
During the 2000s, observations by measurement methods except ERI, hull-contact sensor, and moored and drifting buoys seem to be biased cold, although these do not constitute a major fraction of the observational database. The mean bias relative to drifting buoy increases gradually toward the end of the decade, and the maximum is estimated at −0.2°C. Therefore, the linearly increasing biases from 0° to −0.2°C are subtracted from the questionable data for the 2000s, and the constant bias of −0.2°C is assumed for the data after the 2000s.
b. Mean biases
Recent SST observations are conducted primarily by drifting buoys deployed in the global oceans (Figs. 1, 2). The buoys measure SST directly without moving seawater onto deck or to the inside of a ship. Therefore, buoy observations are thought to be more accurate than either bucket or ERI data. In this study, moored buoy observations are grouped with drifting buoy data as they exhibit similar characteristics, and hull-contact sensor observations are regarded as ERI (Kennedy et al. 2011). In ISSM05, uninsulated-bucket observations before 1942 were adjusted to ERI, but the differences between ERI and drifting buoy data were not considered. In the present study, we regard this difference as a bias in the ERI measurements, and no biases in drifting buoy observations are assumed. A mean ERI bias is computed by comparing global means of ERI and drifting buoy SSTs averaged in the same monthly 5° × 5° boxes for a period from 1996 to 2010. The box averages are computed by arithmetic averaging of observed anomalies from a daily and spatially interpolated climatology. The climatology of ISSM05 is used here, and it will be updated by using bias-corrected in situ observations later (section 6) before the SST analysis.
The mean ERI bias of +0.13°C is obtained and is within the range for the global region listed in Table 5 of Kennedy et al. (2011). The biases appear to vary regionally and seasonally with large biases in the central North Pacific and the southern oceans. However, sampling data are insufficient to attribute these features to the ERI bias. Thus, only the global mean bias is used. A mean bias, +0.07°C, is estimated for insulated-bucket biases. This value is derived from the mean ERI bias and a global mean ERI–bucket difference, +0.06°C. The latter is computed with the same method as used in the ERI bias for the period from 1986 to 2005 in which most of bucket observations are expected to be insulated (section 4a). It is also expected here that the metadata for recent decades are reliable and that measurement random errors reduce substantially by arithmetic averaging. In the subsequent section, we introduce time-varying biases of unknown-type measurements, that takes into consideration the fact that strains appear in global mean SSTs owing to a lack of metadata.
c. Bias correction scheme
As shown in Fig. 1, observations of unknown type occupy a substantial part of the database after World War II. It is expected that the unknown type is one of insulated and uninsulated buckets and ERI. Moreover, bucket-identified observations are a mix of uninsulated and insulated buckets. According to these circumstances, two independent parameters are introduced for the SST bias corrections. That is, when bucket or ERI or unknown observation are given, bias is defined as
where α is the fraction of all unknown observations that are bucket observation; β is the fraction of all bucket observations that are uninsulated; and δiB, δuB, and δERI are biases in insulated bucket, uninsulated buckets, and ERI, respectively. In the case of ERI, α is zero, while a bucket bias is computed with α = 1 and time-varying β. When the type of observation is unknown, the bias is given using time-varying α and β.
As defined in the previous subsection, δiB and δERI are set to 0.07° and 0.13°C, respectively. We used the FP95 correction as the uninsulated-bucket bias δuB in the present study. Actually, the correction is constructed as a function of ship speed and the mixture of insulated and uninsulated buckets. However, the spatial pattern of the correction is reasonable: for example, there are large corrections in the Kuroshio and the Gulf Stream regions in winter, and it appears dominant even though the ship speed and the mixture change in time. In applying this to bucket observations from 1942 onward, we assume that the FP95 correction for 1941 is typical for uninsulated-bucket observations after 1941, although the exact ship speed and the mixture is unknown over the period. Note that all bucket data are subject to the correction with δiB in order to adjust them to drifting buoys.
Two-step computations are performed to obtain values of α and β. First, observations identified as bucket should be corrected with δuB or not. Values of β for only bucket-identified observations are estimated as a function of year from 1942 to 2010. Comparing observed anomalies between bucket and the other known types in the same monthly 5° × 5° boxes, β for each year is uniquely defined such that the global mean of bucket anomalies is equivalent to that of the other known types. We do not rigidly follow the evidences in literature listed in Kennedy et al. (2011) but aim for consistency among existing observations. To reduce observational uncertainty in estimated β, the data are averaged over 5 yr, centered at the year of the computation. Furthermore we assume that the ratio, β, is preserved for bucket observations reported as unknown type, expecting that the estimated β values reflect reality. Second, values of α are estimated by comparing the box-averaged data of unknown type with those of known types as did in computing values of β. Here, the observations identified as bucket are also used with the values of β obtained above.
For the period after the 1970s, more than 80% of the bucket observations are used in the computation of β. The estimation is probably unreliable particularly for the period before 1960 because of the sparseness of metadata (cf. Fig. 2). The situation in the computation of α is similar to that of β. Hence, we introduced two additional constraints for the computation of α and β. One is that the global annual mean SST time series is close to that of box-averaged global mean nighttime air temperature (NAT) for this period. The other is that values of α and β are prescribed for specific periods at first; all unknown data before 1939 are required to apply the FP95 correction, 95% of unknown is ERI during World War II from 1942 to 1945, and bucket subject to the FP95 correction occupies less than 5% of the total after 1975. Moreover, during the period from 1939 to 1941, the FP95 correction needs modification (Rayner et al. 2006) because of ERI observations became dominant gradually.
For years before the 1970s, β is modified by comparing global mean time series of SST anomaly of all type with that of global mean NAT anomaly, assuming that the most part of the global mean SST–NAT difference is originated from the uninsulated-bucket observations. Next, α is computed by the aforementioned procedure with the modified β. If necessary, β are adjusted by every 5% so that α takes a value in the reasonable range. Furthermore, it is assured that α and β are continuous with the prescribed values. Large biases in NAT, as reported by Bottomley et al. (1990), are subtracted before comparison by applying the correction proposed by Rayner et al. (2003), but the global mean NAT was scarcely used over the period of World War II. As a result, the coefficients were determined by the prescribed values in wartime that connect with the postwar estimates continuously in time.
Since the FP95 correction is relative to SST climatology for the period from 1951 to 1980, it has to be adjusted to a new climatology defined by bias-corrected SST observations. This adjustment is achieved by subtracting a global mean difference between the old and new climatology of SST analysis from the original FP95 correction. The value of the difference is estimated at +0.02°C from the mean biases of ERI, and insulated and uninsulated buckets, taking into account spatiotemporal distributions of individual observational types (Fig. 4a) and the weights (Table 1). The FP95 correction adjusted to the new climatology is used for the above-mentioned computation of α and β.
Alternatively, another approach using “local” SST–NAT differences is applicable to determine α and β for the entire period, as Smith and Reynolds (2002) did. However, air temperature observations are not always reported simultaneously with SST observations, and instantaneous SST–NAT differences are affected by large air–sea interaction sometimes. Therefore, we preferably relied on consistency among instrument types of SST in determining α and β.
d. Time series of bias
Using values of α and β estimated above, annual mean weights of individual observational types used for constructing global mean SST are computed as shown in Fig. 4a. The weights are a function of the number of data, spatial data distribution, and the instrument type (see section 3). There are large changes in the contribution of bucket observations to global mean SST until the late 1970s, and they contribute around one-third by 1990. Bucket observations are dominant in the global SST before 1940, in the latter half of the 1940s, and around 1970, primarily because bucket observations were particularly common during these periods. The ERI observations determine around a half of the global means except before the 1940s and during the 1960s and 1970s. The contribution of drifting buoy has grown since the 1980s and it reaches approximately 70% of global mean SST in recent years. These features agree with Fig. 2 of Kennedy et al. (2011). The bias correction scheme is not intended to attribute each observational report to a plausible instrument type, but it appears to be consistent with earlier estimates.
Figure 4b shows time series of annual global mean SST bias correction. The total corrections are within ±0.1°C with the exception of the period before 1940, and negative corrections are dominant from the mid-1960s to 2010. Biases in uninsulated buckets dominate the period before 1960, except for the period of World War II. The global mean corrections to bucket observations are mostly zero after the 1980s, with the ERI corrections becoming dominant since the mid-1960s. The estimated confidence intervals of the annual mean correction at the 95% level are shown in the figure, assuming that the majority of the uncertainty originated from uncertainty in the global mean SSTs used for the computation of α and β.
5. Ice–SST relationship
In historical SST analyses conducted so far (Reynolds et al. 2002; Rayner et al. 2003; ISSM05), SIC has been used as a proxy for SST observation at high latitudes where SST observations are scare. In ISSM05, the proxy SST is estimated using quadratic functions of SIC that reflect the empirical relationship between SIC and observed SST at the same positions on the same day. Despite the wide use of the ice–SST relationship in the historical SST analysis, little attention has been paid to the dependency of freezing points on sea surface salinity. In this study, such salinity-induced contrasts of freezing point in space and seasons are introduced to the empirical functions used in the previous analysis. The freezing point of seawater (FPS) is high at low salinity. In the previous analysis, a constant FPS of −1.85°C was adopted.
The ice–SST relationship is given by the following quadratic function of SIC I,
where T is an estimated SST and a, b, and c are coefficients of the function. The coefficients vary on a monthly basis on a 5° × 5° grid. The second term on the right-hand side of Eq. (2) is not considered in ISSM05. This term can represent colder SSTs than those found in regions of higher SIC. As mentioned above, an additional constraint is introduced when determining the coefficients
where a freezing point Tf at salinity S is given by Millero (1978) and Im is maximum SIC, which is 0.95 in this study. Here, S is the climatological monthly salinity from the World Ocean Atlas 2005 (WOA05; Boyer et al. 2006). As salinity observations are much scarcer than those of SST, the interannual variation of salinity is not considered here. Compared with the constant of −1.85°C, the salinity-dependent FPS is higher by 0.2°–0.5°C off the Siberian coast and east of the Hudson Bay, by 0.2°C in the Beaufort Sea, and by more than 0.5°C in the Gulf of Bothnia, while it is almost the same as the constant in the oceans around the Antarctic continent. The coefficients are determined in a least squares sense by comparing in situ SST observations with SIC observations on a daily and global 1° × 1° grid.
Figure 5 shows the observation frequency of in situ SSTs in highly ice-covered regions, where SIC is greater than 0.7, as a function of FPS of the climatological salinity. Seawater freezes at a certain temperature below 0°C. Many pairs of observed SST and climatological FPS data are distributed along a line that indicates that the two temperatures are identical. Because the observed SSTs are marginally higher than the climatological FPS, the peaks of sampling in most 0.5°C-divided bins appear to be located above the line. This is because the samples were collected in regions where the sea surface is not fully covered by sea ice. The figure suggests that the globally constant FPS of −1.85°C results in cold biases in regions of low salinity and that the empirical ice–SST conversion used in this study provides better spatial SST distributions at high latitudes than that of ISSM05. The ice–SST samples are rather scarce at freezing point of −1.0°C, because the salinity of this range appears only in the Laptev Sea.
In the objective analysis, the SST proxy translated by Eq. (2) is treated as an SST observation. The observational error of the proxy is larger by 50% than that of ship observation when observed SIC is used and 3 times larger when only climatological sea ice is available (Table 1). These factors are determined by considering the errors in SSTs derived from the ice–SST relationship [Eq. (2)]. Because SIC data are given at 1° × 1° grid points and hence correlated highly among the data, a data-thinning procedure is employed. Here, local peaks and peaks along coasts and sea ice margins are kept with the highest priority, and no more than one observation is selected within 60 nautical miles. This also benefits a reduction of the computational cost. The SST around sea ice is optimally computed by the objective analysis so that it is close to the climatological FPS. This does not necessarily yield the same SST as the freezing point at the maximum SIC when in situ observations are available there.
6. Objective analysis
The new SST analysis is given daily as the sum of secular trend, interannual variations, and daily changes on a global 1° × 1° grid from 1845 onward. We refer to this analysis as “multi-time-scale analysis” (MTA). The climatology of SST is redefined as the mean of an SST analysis by OI of ISSM05 with bias-corrected in situ observations for 1961–2005. The new climatology is slightly cooler than the old one. The standard deviations of interannual SST variations are the same as those used in ISSM05, which provide background and observational errors for the objective analysis. A temporally variable climatology (i.e., the climatology plus the secular trend) is used in the analysis scheme and in the quality-control schemes (section 3). This may help us to realize an unbiased analysis, in particular, by avoiding this problem that the objective analysis scheme reduces amplitudes of climate signals of interest too much. In some cases, observations deviating far from a temporally fixed climatology are salvaged in the quality-control procedure. The secular trend is defined in the subsequent section. Interannual variations of SST are reconstructed with empirical orthogonal functions defined from an SST analysis with in situ and satellite observations over the period from 1961 to 2005. The reconstruction generates SST variations in areas devoid of observations. The daily analysis provides not only day-to-day changes but also variations that are represented by higher EOF modes than those used in the reconstruction. Further details of each analysis are provided below. As in ISSM05, the analysis errors are also computed daily on each grid point.
The satellite observations are not used in the final products in order to preserve the homogeneity of the SST analysis for more than 100 yr. If satellite data are included in the MTA, our experience is that there is an inhomogeneity in the gridwise variances of SST analysis, which increases by between 10% and 20% at the point at which the satellite data start.
An OI scheme is also used for comparison with MTA. This scheme is the same as that of ISSM05, except that the secular trend is prescribed as in MTA. The SST analysis with the satellite and in situ data is conducted according to the OI scheme.
When interpolating the monthly climatology to the day of the objective analysis, we introduce a procedure adopted in the Atmospheric Model Intercomparison Project (Taylor et al. 2000). To do this, another monthly climatology is introduced; it is constructed from the original climatology such that the monthly mean of daily SSTs interpolated from the new climatology is identical to the original. Thereby, the amplitude of the annual cycle of the new climatology is larger than that of the original particularly in the North Pacific and the North Atlantic Oceans. If we do not use a climatology with a boosted seasonal cycle, sizable biases exceeding 0.5°C can appear in these regions of the analyzed SST fields.
As in ISSM05, the other marine meteorological variables, sea level pressure (SLP), surface wind, air temperature, cloudiness, and dewpoint temperature, are also analyzed. However, only an SLP analysis (COBE-SLP2) is presented below.
a. Trend analysis
The trend is the leading EOF of annual mean SST anomalies of bias-corrected in situ observations averaged in 5° × 5° boxes from 1845 to 2010. The missing values are replaced by estimates from a reconstruction scheme. Here, differently from the reconstruction of MTA (section 6b), the EOFs are computed from 5° box averages of satellite and in situ SST data from which the global means are subtracted. Data over sea ice regions are not used for the computation. It is confirmed that the global means of the reconstructed 5° data are kept unchanged within the uncertainty. From the reconstructed box averages, a trend EOF was computed. By doing this, the trend EOF is expected to be representable during the whole period from the mid-nineteenth century onward; otherwise, the spatial pattern of EOF is enhanced by gridwise trends for later data-rich decades. Figure 6 shows the spatial pattern of the leading EOF mode and time series of the principal components. Similar to linear trends of HadISST (Rayner et al. 2003) and extended reconstructed sea surface temperature, version 3b (ERSST.v3b; Smith et al. 2008), areas of positive trends (e.g., in the Indian Ocean, northern Pacific, and the tropical and southern Atlantic) are commonly seen in the leading EOF. The time series does not change monotonically over the entire period but an increasing trend is dominant from 1970 onward.
Values of the secular trend at grid points on the date of objective analysis are given simply by multiplying the EOF by a principal component temporally interpolated to the analysis date. From 2011 onward, the values of the principal component of 2010 are used.
From the detrended SST analysis by OI with satellite and in situ observations from 1961 to 2005, we computed EOFs explaining 90% of the total SST variance on interannual-to-interdecadal time scales. The purpose of reconstructing SST on a historical base is to prepare the SST analysis useful for various climate studies with dynamical models. Therefore, we use higher EOF modes so that the reconstructed SSTs contain as much variances as those for recent decades. The reconstruction with the higher modes produces partly signals and partly uncertainty, depending on the data distribution on the analysis date. Before the computation of the EOFs, the variance of the satellite SST analysis is adjusted to that of in situ SSTs by multiplying the former by a globally unique and seasonally changing factor of about 0.8. Thanks to this procedure, jumps of gridwise SST variances between before and after 1981 were minimized less than 1% of standard deviation of interannual SST variation on average; otherwise, they were 15%. With these EOFs, SST anomalies on a monthly time scale are estimated daily using in situ observations within 15 days either side of the analysis day. To explain 90% of variance of SST variation requires 133 EOF modes. An optimum SST anomaly field x minimizes the following cost function J,
where is a background error covariance matrix, is a diagonal matrix of observational error variances, is a bilinear interpolation operator from the grid point to an observation location, and y is a vector containing observed anomalies. In the present analysis, is given by the set of EOFs : that is,
where is a diagonal matrix containing SST variances associated with the individual EOF modes. The observational errors are defined the same as in ISSM05, which are standard deviations of the interannual SST variation multiplied by weights that depend on the observation methods employed. When merging nearby data, correlations between observations are considered as described in section 3. The solution x is computed by using a conjugate gradient method (Derber and Rosati 1989; Ishii et al. 2003). The computational cost of searching for the solution is low because it is realized in a phase space of 133 degrees of freedom spanned by the EOFs.
c. Day-to-day analysis
Finally, a daily changing component is added to the trend plus the interannual SST fields computed above. Although the reconstruction is conducted on a daily basis, it does not produce substantial daily SST changes. Day-to-day changes relative to 31-day mean SSTs are computed using an OI method as departures from those on the previous day. Background variances of the daily changes are represented by variances of deviations from 31-day moving averages of the daily SST analysis with satellite and in situ observations. Here, the deviations are multiplied by the aforementioned factor of 0.8. The standard deviations of daily change are 30%–100% of those of monthly SST, and the daily variability is large in summertime (Fig. 7) when ocean surface mixed layers are thin. The spatial and temporal decorrelation scales are set to 600 km and 5 days, respectively.
The daily analysis is incorporated for inclusion of high-frequency SST variations in the analysis particularly for recent data-rich decades. Thus, the quality of the daily analysis is not homogeneous in both time and space in the product. For data-sparse periods, an additional reconstruction for reasonable day-to-day changes must be needed for a homogeneous quality, but this will be a theme of our future study. The day-to-day analysis compensates some amount of the missing part (10%) of monthly SST variance in the EOF reconstruction. The background errors prescribed for the OI scheme include this amount. Accordingly, the sum of the daily analysis over the month is not equal to zero. In addition, we confirmed that the daily analysis compensates the lack of monthly SST variability reconstructed by the EOFs.
d. Analysis error
This section summarizes the analysis error estimation in MTA. As done in ISSM05, the analysis errors are calculated for the SST analysis at all grid points every day. When the OI schemes are used, the analysis error is computed following the OI theory (e.g., Ghil and Malanotte-Rizzoli 1991). Daily analysis errors of MTA are computed as a sum of errors in the EOF reconstruction and the day-to-day analysis. In the case of the reconstruction, analysis error matrix is given by
Similar to the case of Eq. (4), diagonal components of , are computed directly in the phase space spanned by the 133 EOFs. The analysis error of the day-to-day analysis is computed by the OI scheme. When the errors of the monthly mean SST analysis are computed, the correlation between the daily analyzed SSTs is taken into account. Here, a temporal decorrelation scale is set to 15 days as used in ISSM5. The uncertainty in the secular trend is not considered in the present definition of analysis error. An expected global mean error is 0.02°–0.03°C at a 95% confidence level, as shown in Fig. 6b. In addition to this error, further uncertainty is expected in local trends, but they are not specified in this study.
In this section, we describe the basic performance of the new analysis, focusing on major climate variables such as global mean SST. The SST analyses of the Hadley Center [Hadley Centre Sea Surface Temperature dataset, version 3 (HadSST3); Kennedy et al. 2011] and the NOAA/National Climatic Data Center [ERSST.v3b (Smith and Reynolds 2004) and National Climatic Data Center (NCDC) OI.v2 (Reynolds et al. 2002)] are used for comparison.
a. Global mean SST
The global mean SST time series from 1850 to 2011 is compared with those of HadSST3 and ERSST.v3b in Fig. 8. Since HadSST3 includes missing values in cases where the observations are unavailable, the monthly global means of the present analysis (COBE-SST2) and ERSST.v3b are computed from SSTs only at grid points where HadSST3 is given. Throughout the period shown in the figure, except the 1950s, the global means agree with each other to within 0.1°C. An SST increase of about 0.6°C during the past 150 yr is common to all three analyses. This is also apparent in the global mean of the trend (section 6a; Fig. 6), implying that the bias corrections and trend estimation are crucial in reproducing the global warming trend. Furthermore, HadSST3 and COBE-SST2 agree well despite the different approaches adopted. These two global mean SSTs decrease over the period from 1940 to 1980. Conversely, the previous analysis (COBE-SST) exhibits an increasing trend in this period, similar to ERSST.v3b. In addition, COBE-SST is generally colder than the other time series before 1940 and its SST warming in recent years is rather weak. This is due to the bias corrections after 1941, which are not adopted in ISSM05.
b. ENSO indices
Figure 9 shows SST anomalies averaged over the Niño-3 region (150°–90°W, 5°S–5°N; referred to as Niño-3 index here) of HadISST and COBE-SST2, alongside the Southern Oscillation index (SOI) provided by the University Corporation for Atmospheric Research (UCAR; Trenberth 1984) and that computed from COBE-SLP2. SOI of COBE-SLP2 is computed as the difference in SLP at the two grid points nearest to Darwin and Tahiti. The time series after 1990 are not shown because it is in good agreement with each other. The time series of the two indices correspond favorably to those of the other centers from the late nineteenth century onward. However, the amplitudes of Niño-3 SST anomalies are small and SOI becomes noisy for the period before the mid-1870s. This is probably the effect of data sparseness during the period and will be discussed in section 8. There are some gaps in SOI between the UCAR and COBE-SLP2 data for the 1880s and around 1930. The SOI of COBE-SLP2 appears to vary in phase with the Niño-3 SST index during these periods.
c. Comparison with satellite SST
Our present and former SST analyses are compared with a satellite SST analysis of NCDC OI.v2 on a monthly basis. Figure 10 illustrates better agreement between the present analysis and NCDC OI.v2 in comparison with the previous analysis. Especially in the southern oceans at middle and high latitudes, the root-mean-square differences (RMSDs) reduce significantly and the poor anomaly correlation coefficients (ACCs) seen in the previous analyses disappear in COBE-SST2 mostly. It is difficult to discern what has contributed to the reduction of RMSDs and the increase of ACC in COBE-SST2 from the figure because many parts of the analysis scheme as well as the observational data have been changed. According to additional experiments in which MTA is replaced by OI, the reconstruction scheme for interannual SST variations is one reason for the improvement and produces realistic SST anomalies in the southern oceans. Although the satellite observations are not used in COBE-SST2, large SST variability in the Kuroshio–Oyashio extension regions and along the Gulf Stream is represented better in the present analysis than the previous analysis. We also observed that the new ice–SST relationship did not affect the SST variations in the polar regions adversely. The impact of the use of the salinity-dependent FPS on the SST analysis is barely perceptible in the figure because the mean fields are subtracted from the analyzed SSTs. In general, the agreement of the new SST analysis with NCDC OI.v2 becomes better over the polar regions. Note that the period of the statistics of ACC and RMSD overlaps with that of the EOFs used for the reconstruction.
d. 1-month lag
The persistence of SST anomalies is an important aspect for the improvement of understanding oceanographic conditions and for the use in various climate studies with dynamic models (Hurrell and Trenberth 1999). The lag 1-month autocorrelation of the SST anomaly of COBE-SST2 is presented in Fig. 11 together with that of NCDC OI.v2. The lag correlations appear to reflect oceanographic thermal conditions. For the period from 1982 to 2011, large lag correlations are observed in the tropics, the North Pacific, and the North Atlantic, as well as in areas around 65°S and 120°W, where SST is affected by ENSO through the atmospheric bridge (Alexander et al. 2002). In contrast, the lag correlations are relatively small along the western boundary currents and in areas where ocean eddy activity is large.
In the present analysis during recent decades, the high lag correlations are reproduced successfully by the reconstruction of interannual SST variations as seen in NCDC OI.v2 (Fig. 11, right), although the satellite observations are not used in MTA. Meanwhile, in our previous analysis, low correlations are dominant south of 30°S because the analysis is consisted only of ship observations by OI. The lag correlations tend to be larger globally when the daily changes are not incorporated. The day-to-day analysis increases the variance on spatially small and short-term scales.
Data sparseness is rather severe for the early twentieth century, particularly in COBE-SST, but both SST analyses still capture spatial structures of the lag correlations characterized by the oceanographic conditions (Fig. 11, left). The lag correlations are still large in the tropical oceans and the region around 120°W and 50°S during this period. In summary, the present analysis seems to exhibit a considerable advantage in its reproduction of the autocorrelation, although the lag correlations for this period cannot be verified.
e. Differences between SST analyses
The performance of the SST analysis and divergences among SST analyses of the different centers are examined on a monthly basis for the period from 1901 to 1930 with the same measures as used in Fig. 10 (Fig. 12). Here, COBE-SST2, ERSST.v3b, and HadISST are compared with each other. In the tropical and northern Atlantic and the northern part of the Indian Ocean, the divergence between the three analyses is small. In contrast, the divergence appears large in the western North Pacific, and furthermore the SST analyses are poorly correlated with each other at high latitudes. These features are common in the statistics for the period before 1900. The data sparseness is expected to be severe during this period, causing large divergence in SST analysis, particularly in areas of low SST variability such as the western tropical Pacific and in eddy-active regions.
8. Uncertainty and reliability
Uncertainty in the analyzed SST on a monthly basis are discussed in this section, focusing on the advantages of the present analysis compared with our previous analysis, the reliability of areal mean SSTs, and uncertainty in the analyzed SSTs. To this end, we conducted a test in which additional SST analyses are performed for years from 1840 to 1990 with MTA and OI, accessing uncertainty due to spatiotemporal data distribution in the two analyses.
Before the test, we constructed another observational dataset; this contains observations in recent years whose daily distributions are similar to those in days from 1840 to 1990. A period of 5 yr from 2006 to 2010 containing abundant observations is chosen as the recent years for the test. With the subsampled observations, analyses by OI and MTA are conducted in the same way as that described in section 6. These analyses are compared with the SST analysis of the recent years based on satellite and in situ data. The subsampled data and their analysis are referred to below as pseudo-observations and pseudoanalysis, respectively. Quality-controlled and merged data used in the pseudoanalysis are subsampled on a daily basis, which are nearest to the date (month and day only) and position of the past data. Duplicate sampling within 15 days is allowed. The number of subsampled observations is almost the same as that of the original datasets for the past years. The 5-yr period from 2006 to 2010 covers almost one ENSO cycle and its related variations, although the ENSO during this period is probably modulated by decadal changes at low latitudes or by the North Pacific Gyre Oscillation (Lee and McPhaden 2008; Lorenzo et al. 2010). Note that the SST analysis of this 5-yr period is statistically independent of the EOFs used for reconstruction of the interannual-to-interdecadal variations relative to the trend. It is also noteworthy that uncertainty in the past SST analysis based on the EOFs representing recent SST variations cannot be evaluated in this test. Here, we focus only on sampling errors.
Root-mean-square errors in both global monthly means and Niño-3 SST estimated from the pseudoanalyses (thick lines) decrease gradually over time (Fig. 13) for both OI and MTA and are substantially smaller than the detrended interannual standard deviations (dotted lines) for the whole period except the 1840s and around the 1860s. Accordingly, it is expected that the global mean SST and Niño-3 SST time series are reliable from the 1880s onward. The Niño-3 SST suffers from data sparseness before 1880, although the Niño-3 SST of COBE-SST2 agrees well with HadISST and SOI in Fig. 9. After 1920, for the global means and after 1960 for Niño-3 SST, the errors are comparable with or smaller than mean root-mean-square (RMS) difference among HadISST, ERSST.v3b, and COBE-SST2 (short and long dashed lines). These RMS differences include systematic errors derived from differences in quality-control and optimization schemes and bias corrections among the centers.
Although the OI method is not a bad choice for estimating global mean SST and Niño-3 SST, the figure confirms several advantages of MTA. The OI scheme suffers from data sparseness during World Wars I and II and in the nineteenth century, while MTA produces more accurate SST fields than OI. The EOF reconstruction is successful in reproducing SST fields under unfavorable conditions.
The figure also shows verification of the analysis error computed according to the objective analysis theory. The theory-based analysis errors (thin lines) given at each grid point are averaged, taking spatial correlation of the errors into account (Ishii et al. 2003). The spatial decorrelation scale is 600 km, which is the same as that used in the OI scheme. The decreasing trends and temporal changes in errors are similar to those of the pseudoanalysis and are much pronounced in MTA than in OI. The uncertainty in Niño-3 SST by MTA is generally equivalent to the RMSD of the pseudoanalysis except for the period before the 1870s. Therefore, the analysis errors generated by MTA are appropriate as a measure of the reliability of SST analysis, although the analysis errors cannot represent large uncertainty in the analysis over the 1860s.
Next, sampling errors in the Kuroshio–Oyashio extension are discussed. This region is a key area in climate change owing to strong air–sea interactions (e.g., Nonaka and Xie 2003; Kwon et al. 2010) and the central location of Pacific decadal oscillation, which are worth monitoring and predicting (Mantua et al. 1998; Mochizuki et al. 2010). A particular target region is located from 120°E to 150°W between 30° and 50°N, where large divergence is observed between the three SST analyses for 1901–30 (Fig. 12). Figure 14 shows RMSD and ACC for MTA and OI against the SST analysis with satellite and in situ–based psudo-observations for the historical period. From the mid-1910s to 1990, excluding the period around World War II, the two pseudoanalyses are less uncertain than in an earlier period before 1910 and are highly correlated with the satellite SST analysis. Long-term climate changes are probably detectable with satisfactory accuracy during the twentieth century. An advantage of MTA over OI is seen also in this case, but RMSEs before 1910 exceed the standard deviation of the climate signal. Over the whole period, there appears a systematic difference between the pseudoanalyses and the satellite analysis because of the large SST variability of the satellite analysis. Such variability is not seen in either ERSST.v3b or HadISST. Thus, a standard deviation of the excess of variability due to the satellite data is estimated at about 0.2°C.
Our previous historical sea surface temperature (SST) analysis (ISSM05) is replaced by a new analysis scheme, MTA, consisting of the three analyses on secular, interannual-to-decadal, and daily time scales. The secular trends are estimated first, and interannual-to-interdecadal SST variations are reconstructed on a daily basis using empirical orthogonal functions defined from a 45-yr OI analysis with satellite and in situ observations. Moreover, day-to-day changes are analyzed with observations available during the day of the analysis. The analysis errors are provided by the theory of objective analysis. An empirical equation representing the ice–SST relationship is also updated, taking climatological salinity distribution into account. The new analysis has an advantage over the previous method in climate studies with dynamical models, because additional variances of SST are presented even in regions with data-sparse coverage. In a test evaluating sampling errors, the present analysis scheme reduces the uncertainty in analyzed Niño-3 SST anomalies during World Wars I and II. Concurrently, the analysis errors should be taken into account when using the analyzed SSTs because the SST analysis is not always reliable.
Bias corrections of bucket and engine room intake (ERI) are crucial in determining secular trends. However, a lack of metadata prevents us from applying an appropriate bias correction directly to the observations in question. A method of the bias corrections proposed in this study is to remove expected mean biases in individual measurement types from historical observations. When the measurement type is not recorded in the observational database, the bias is defined as that for a mixture of ERI and insulated- and uninsulated-bucket observations. The ratios of these three types vary in time. Consequently, the biases are estimated from all available SST and nighttime temperature observations so that all global means of each type or element are consistent.
Our future study will focus primarily on uncertainties in historical SST changes and applications of the analysis errors to climate studies. There remains room for improvement in the handling of historical SST bias correction. In particular, the accuracy of absolute SST values should be improved upon further.
Observational data used in this study are obtained from the following centers: ICOADS from CDC/NOAA; AVHRR data from PO.DAAC/NASA; SSM/I, SSMIS, and SSMR brightness temperature data from NSIDC of the University of Colorado; HadSST3 and HadISST from the Met Office; ERSST.v3b and satellite SST analysis from NCDC/NOAA; SOI index from UCAR; and Arctic sea ice concentration data from the University of Illinois. The authors are grateful to Dr. Kent (NOCS/UK) for the latest information regarding WMO publication 47. They also thank Dr. Komuro (JAMSTEC) for valuable information on recent sea ice changes. Their appreciation is extended to two anonymous reviewers for many insightful comments. This work is supported by the Innovative Program of Climate Change Projection for the 21st Century and the Program for Risk Information on Climate Change of the Ministry of Education, Culture, Sports, Science, and Technology, Japan.