Interannual variability in upper ocean heat content, temperature, and thermosteric expansion on global scales

2.1. Profile Data

[8] The set of in situ data was compiled from several different archives and includes temperature profiles from XBTs, CTDs, profiling floats, moored buoys (primarily from the TAO array), and autonomous pinniped bathythermographs. Approximately 1,000,000 unique profiles were compiled from the historical archives for the period from 1993 through 2003. Profiles were retrieved from the World Ocean Database 2001 (WOD01) [Conkright et al., 2002] as well as the Global Temperature and Salinity Profile Program (GTSPP), both of which are maintained by the National Oceanographic Data Center. In addition, float profiles were obtained from both the World Ocean Circulation Experiment and the Argo project. Although these data sets contain substantial overlap, each set contained enough unique profiles to warrant retrieval and processing. It was therefore necessary to detect and remove duplicate profiles. Profiles were considered duplicates if they were colocated within about 1 km in space and 2.5 hours in time. Figure 1 shows the profile distribution for a typical year. Note the paucity of profiles in the Southern Hemisphere. Also shown in Figure 1 is a plot of profile availability for each year.

[9] Quality control flags provided by the data centers were used to remove spurious profiles whenever possible. Some data, however, did not undergo quality control, and enough obvious errors remained that further quality control was deemed necessary. Profiles were grouped into modified, 10° × 10° World Meteorological Organization (WMO) squares. Squares that contained small amounts of data in geographically similar regions were consolidated. In addition, some squares were split in order to group profiles with similar properties. The result was about 300 geographically distinct regions containing anywhere from a few hundred to tens of thousands of profiles. All profiles in each square were visually inspected en masse in order to remove gross outliers. A mean profile was then computed, and any profiles with points more than 6 standard deviations away from the mean were removed. These procedures resulted in the removal of about 3% of all the profiles.

[10] In addition to removing spurious data, profiles with no data below 350 m were also discarded. A significant number of the remaining profiles ended at depths shallower than 750 m, the depth to which the analysis was performed. About one fourth of the profiles are complete to 750 m, and slightly more than half have data to 500 m. A regression technique described by Smeed and Alderson [1997] was used to “extend” these short profiles to 750 m. This was done separately for each modified WMO square in order to allow the regression coefficients to vary from region to region. Enough complete profiles were available in almost all regions to reconstruct the short temperature profiles as well as characterize the error due to the reconstruction. For depth-integrated quantities such as heat content, the errors due to the reconstruction represent about 8% of the variability in the regional boxes. This translates to approximately 8 × 10⁸ J/m² in heat content or 0.9 cm in thermosteric expansion. In temperature, errors below the cutoff depth are approximately 0.2°C or 30% of the variability.

[11] Finally, it was necessary to correct a large number of profiles for an error that was discovered in the fall rate equations of a widely used set of XBT probes. The linear depth correction suggested by Hanawa et al. [1995] was applied to all profiles that were determined to need correction. Prior to mid-1996, however, it was not common for information regarding XBT probe type to be recorded, making it difficult to determine exactly which profiles should be corrected. Both the GTSPP and WOD01 data sets provide information to assist users in determining whether or not to correct a given profile. In the GTSPP database, each XBT profile was assigned a status of either “needs correction,” “no need to correct,” or “unknown probe type.” The correction was applied to all profiles determined by GTSPP to need it. In the WOD01 data set, only information regarding probe type was given. In this data set, all XBT probes of type T-4, T-6, T-7, and Deep Blue were corrected. Unknown profile types in both data sets were treated in the same way. Following the criterion used to produce the WOA01 gridded products, all unknown XBT profiles with maximum depths less than 840 m were assumed to need correction [Conkright et al., 2002]. In all, XBTs account for about 400,000 profiles, and just over half of these required correction. About 100,000 XBT profiles were of unknown probe type, with almost all of these occurring prior to 1996. Since equipment that automatically applied the corrected depth equations was not in widespread use prior to 1996, it is likely that almost all of the unknown profiles require correction, and as expected, the maximum depth criterion selects all but about 10,000 of the unknown profiles to receive the correction.

2.2. Altimeter Data

[12] A merged, gridded product produced by AVISO and containing TOPEX/Poseidon, Jason 1, and ERS 1 and 2 data [Ducet et al., 2000] was used for altimetric height anomalies. This product provides sea surface height anomalies relative to a 7-year mean from 1993 through 1999. It consists of maps produced every 7 days on a 1/3° × 1/3° Mercator grid. Because the maps contain data from multiple altimeters, they are able to resolve variability on scales as small as 150–200 km with an accuracy of a few centimeters over most of the globe [Ducet et al., 2000]. The high grid resolution of the maps substantially oversamples this variability so that linear interpolation (as opposed to a more complicated spline fit or other technique) was found to be adequate for estimating altimetric height anomalies at the time and location of individual profiles.

2.3. Interannual Variability and the Difference Estimate

[13] The time-mean and seasonal cycle were removed from both the altimeter and in situ data prior to analysis. Removal of the time-mean was necessary in order to reduce error in the altimeter data associated with uncertainty in the geoid. The seasonal cycle was removed in order to emphasize interannual variability and the decadal trend. Although altimeter data were provided as anomalies relative to a 7-year mean, a 10-year mean over the anomalies was also calculated and subtracted. This was done to ensure that altimetric height anomalies and anomalies calculated from the in situ data were referenced to the same temporal baseline (January 1993 through December 2002). The seasonal cycle was estimated at each grid point from 3-month bin averages over all years of altimeter data. The seasonal cycle was calculated in this way to ensure agreement with in situ anomalies that were computed relative to the WOA01 seasonal cycle, as described below.

[14] To remove the time-mean from the in situ data, objective maps were made using profile data from January 1993 through December 2002. The maps were computed using the World Ocean Atlas 2001 (WOA01) climatology as a first guess in order to minimize the variance lost by the objective interpolation. Mapping was done for heat content as well as temperature on level depths every 10 m from the surface to 750 m. As in work by Willis et al. [2003], a two-scale covariance function was used in the mapping procedure,

and was chosen to be consistent with the global wavenumber–frequency spectrum for altimetric height published by Zang and Wunsch [2001]. The maps were estimated at each grid point using the nearest 300 data points along with an additional 50 points selected randomly from a 15° × 15° region to help resolve variability at larger scales. These procedures and mapping parameters were used for all objective maps throughout the work. To remove seasonal variability from the in situ data, the WOA01 seasonal anomalies were used. Seasonal cycles in heat content and thermosteric expansion were calculated from the WOA01 temperature anomalies. The seasonal anomalies are provided as 3-month bin averages, consistent with the processing for the altimeter data.

[15] Maps depicting interannual variability in heat content, temperature, and thermosteric expansion were calculated using the “difference estimate” developed by Willis et al. [2003],

Here XBT represents the quantity to be estimated (heat content, temperature, or thermosteric expansion), AH is satellite altimetric height, α is the time-averaged regression coefficient of AH onto XBT, and brackets denote objective mapping. The two bracketed terms on the right-hand side of (2) are referred to as the “difference field” and the “synthetic estimate,” respectively. The estimate amounts to an objective map of the quantity XBT using the synthetic estimate (α AH) as the initial guess for its variability. Note that the synthetic estimate has the same temporal and spatial resolution as the altimeter data, whereas the difference field can only be evaluated where in situ data are present. In regions where profiles are abundant, the objective map pulls the initial guess into agreement with the in situ data in a smooth and statistically consistent way. In regions with no in situ data, the map tends to zero and the difference estimate reverts to the synthetic estimate. As discussed below, error bars on the difference estimate were calculated using a twin experiment in which altimeter data were treated as truth.

[16] To produce difference estimates of the various quantities, heat content anomaly, thermosteric expansion, and temperature anomaly were first computed from each in situ profile. Heat content was calculated as the integral from 0 to 750 m of ρC_pT(z)dz. Thermosteric expansion was calculated as the integral over depth of ΔTγ(S_c,T_c + ΔT/2,P)dz. Here γ is the thermal expansion coefficient of seawater, S_c is climatological salinity (S(x, y, z)) from WOA01, T_c is temperature from the 10-year mean, and ΔT is the in situ temperature anomaly relative to the 10-year mean. The mapping procedure discussed above was used to create maps depicting 1-year average anomalies of each quantity on a 1° × 1° grid. A -year resolution in time was used in order to oversample the interannual variability.

[17] The coefficient of regression α in (2) was calculated separately for each of the modified WMO squares. This allowed the coefficient to vary slowly with latitude and longitude in order to reflect the varying vertical structure of temperature anomalies occurring in different regions. Figure 2 shows a contour plot of the regression coefficient for heat content. Also shown is the correlation coefficient (r), the square of which represents the fraction of variance explained by the regression. The regression coefficient is fairly uniform throughout the tropics, increasing somewhat through the midlatitudes in most basins. A notable exception, however, is the Labrador Sea, where the regression and correlation coefficients drop to zero. This reflects an interannual signal in this region described by Antonov et al. [2002] marked by a haline contraction that largely compensates for the thermal expansion present there. The resulting signal in altimetric height is small and poorly correlated with heat content. Regression coefficients for temperature and thermosteric expansion have similar geographic distributions and are therefore not shown.

3.1. Heat Content

[18] The “difference estimate” combines in situ data with satellite data to produce improved estimates of heat content, temperature variability, and thermosteric expansion. The estimate relies on a time-averaged linear regression of the variable of interest onto altimetric height to provide an initial guess for its variability. The objective map pulls the initial estimate into agreement with in situ data (where it is available) in a smooth and statistically consistent way. The correction to the synthetic estimate, or difference field, represents upper ocean variability that cannot be adequately described using only the altimeter data and regression coefficients.

[19] Figure 3 shows the interannual variability of globally averaged, upper ocean heat content computed using the difference estimate. A considerable warming trend is visible in the record. The 10-year heat increase from mid-1993 to mid-2003 implied an average warming rate of 0.86 ± 0.12 watts per square meter of ocean, or 0.29 ± 0.04 pW for this period. Also shown in Figure 3 are the heat content signals of the difference field and the synthetic estimate, the two terms on the right-hand side of (2). A significant amount of interannual variability is present in the time series of the difference field as well as a negative trend of −0.49 W/m².

[20] In order to calculate error bars for the global heat content time series, the two bracketed terms on the right-hand side of (2) are considered separately. First, we consider the final term of (2), the synthetic estimate. Recall that the synthetic estimate is given by the linear regression of altimetric height onto heat content. Regression error will be accounted for in the difference field term, so the error due to the synthetic estimate is given by the random error in the 1-year mean of globally averaged altimetric height. Since the satellites provide near-global coverage and resolve variability on scales greater than a few hundred kilometers, this error is expected to be small. Note that the effect of potential systematic drifts in the altimeter data is discussed in 8. With the annual cycle removed, the RMS variability in globally averaged sea level about the 1-year mean is 2.5 mm. The TOPEX/Poseidon altimeter samples the global mean once every 10 days, giving a standard error on the mean of 0.4 mm for the 1-year mean of globally averaged sea level. Multiplying by α, this gives a random error of 0.6 × 10⁷ J/m² for the second term in (2).

[21] Next, we consider the error in the first term of (2), the difference field. Recall that the difference field is equal to the misfit of the linear regression of altimetric height onto in situ data. The difference field can only be evaluated, however, at the times and locations where profiles are available. Undersampling of the global mean by the in situ data is therefore the largest source of error in the difference field term. A twin experiment that treated the altimeter data as truth was used to estimate the effect of this undersampling. Although the altimeter data do not contain variability on scales shorter than about 200 km, they do provide a very accurate measure of globally averaged sea level and are therefore an excellent tool for estimating the effect of subsampling on globally averaged quantities. First, altimetric height was interpolated to the time and location of each in situ profile. The subsampled altimeter data were then objectively mapped using the same technique as in the difference estimate. The resulting 1-year average maps were then compared to 1-year average maps from the complete AVISO altimeter data set. The differences represent the errors due to undersampling by the in situ data set. Averaging over all years gives an RMS error of about 1.7 mm for the global mean, or about 24% of the RMS signal strength of the time series of globally averaged altimetric height. This suggests that the relative error in predicting the global average of a quantity from in situ data alone is about 24%. Thus the relative error in the difference field is expected to be 24% of its variance. Estimating the variance of the difference field, however, is somewhat difficult. The objective maps tend to underestimate the variance of the difference field because they tend toward zero far away from the in situ data. To account for the reduced variance, we turn again to the altimeter twin experiment. In the twin experiment, the subsampled estimate of the global mean suffered a 14% reduction of variance relative to the AVISO maps. To account for this in the error bar, the estimated variance of the difference field was increased by an appropriate factor before the 24% error bars were estimated. This resulted in a 1.2 × 10⁷ J/m² error bar for the difference field term of the globally averaged heat content time series. Combining this in quadrature with the 0.6 × 10⁷ J/m² error in the last term of (2) gives a random error of 1.3 × 10⁷ J/m² for the difference estimate. Random error in the estimates of the other quantities was calculated in the same way.

[22] To calculate the error in the average warming rate, the random error in the global mean heat content was used as follows:

where Δ(HC) is the error in an individual point of the heat content time series. This implies a 0.06 W/m² error bar for the warming rate. As discussed at the end of 8, however, we must account for the possibility of a small systematic error, or bias, of approximately 0.06 W/m² in the warming rate. Because this error is systematic, it adds directly to the 0.06 W/m² random error bar, giving the 0.12 W/m² error bar quoted above. Finally, the error bar for each point of the heat content time series is also increased using (3) to account for the potential systematic error. This gives the roughly 2.4 × 10⁷ J/m² error bars shown in Figure 3.

[23] Although the altimeter data provide an excellent tool for testing the effects of undersampling on the estimation of long-term trends, the 0.12 W/m² error bar is small, and another check was performed to ensure that the suggested error was reasonable. To do so, a global map of the 10-year change in heat content was first calculated using the difference estimate. For comparison, a similar map was produced from an estimate of heat content made using only in situ data. For the “in situ only” estimate, heat content was objectively mapped directly from the profile data for each year in the time series, again using the mapping technique discussed above. The 10-year change in heat content was then calculated from the resulting time series of in situ only maps. Figure 4 shows the two maps of 10-year heat content change along with the data availability over the time series. Also shown are the zonal integrals of the maps in watts per meter of latitude. The zonal integral was used rather than the zonal average in order to reflect the contribution of each latitude band to the globally integrated (or averaged) rate of heat storage. It is clear that most of the signal in the 10-year change can be resolved from in situ data alone. Only in the Southern Ocean, where data are sparse, does the in situ only estimate fail to recover much of the signal. The global rate of heat storage from the in situ only estimate is 0.81 W/m². This agrees within error bars with the 0.86 ± 0.12 W/m² warming rate calculated using the difference estimate. Furthermore, the smaller in situ only warming rate is consistent with the variance reduction implied by the altimeter twin experiment when in situ data are used to estimate a globally averaged quantity. In addition to supporting the 0.12 W/m² error bar quoted above, this calculation illustrates how the use of altimeter data to form an initial guess reduces the problem of variance lost in the objective interpolation of undersampled data.

[24] Figure 5 shows the interannual variability in globally averaged heat storage (the time derivative of heat content). The curves shown in Figure 5 were calculated from 1-year, centered differences of the heat content curves in Figure 3. Note that the heat storage signal contains significant interannual variability. In particular, the time series shows a rapid warming during the onset of the 1997–1998 El Niño, followed by a slight cooling in the latter part of 1998. Error bars for the difference estimate of heat storage can also be calculated using (3) by changing the time step to 1 year. This gives an error bar of 0.6 W/m² for the time series of the difference estimate shown in Figure 5.

[25] Although they contain significant differences, the maps of heat content used to produce these time series are qualitatively similar to 1-year average maps of altimetric height. Figure 6 shows maps of heat storage calculated by differencing the heat content maps.

[26] A large amount of interannual variability is present in the maps with the largest signals due to the onset and decay of the 1997–1998 El Niño. It is evident from the maps that interannual variability in the global average is strongly influenced by variability in the tropics, particularly in the Pacific and Indian Oceans. Indeed, much of the interannual variability in the time series of the globally averaged heat storage is related to the El Niño Southern Oscillation (ENSO). For instance, the 1995 panel shows a substantial warming in the western tropical Pacific that marks the end of persistent El Niño conditions that lasted through the early 1990s. A similar pattern of warming is evident in 1999 and marks the onset of a persistent La Niña that remained through the latter part of the time series.

[27] In order to better illustrate the patterns of variability depicted in Figure 6, it is helpful to compute zonal integrals over the maps of heat content. Figure 7 shows a time-latitude plot of zonally integrated heat content variability in Joules per meter of latitude. Again, the most prominent feature is the 1997–1998 El Niño. Subsequent to this event, it is also possible to see the poleward propagation of a positive heat content anomaly, beginning in mid-1997 at about 20°S and continuing to 30°S by the beginning of 2000. A similar, although weaker, signal can be seen in the Northern Hemisphere at approximately the same time, suggesting that heat from the tropics propagated poleward into midlatitudes subsequent to the large 1997–1998 El Niño event. The importance of the tropical Pacific to interannual variability on global scales is further illustrated by Figure 8, which shows heat content in Joules integrated over the region from 20°N to 20°S. For reference, global heat content is also shown. The heat content in the tropics remained fairly stable through the peak of the El Niño event from late 1997 through early 1998. It then decreased rapidly through the end of 1998 and the first half of 1999. The global integral also shows heat lost during this time, but somewhat less so, suggesting that some of the heat may have been exported to the midlatitudes. The cooling, both globally and in the tropics, gives way to rapid warming in mid-1999 with the onset of La Niña. Without similar estimates of oceanic heat transport and/or surface heat flux, it is not possible to determine whether the propagation features seen in Figure 7 are the result of oceanic heat transport or are simply forced responses to signals in the atmosphere. It is clear, however, that ENSO-related variability can be large enough to affect the global oceanic heat budget on interannual timescales.

[28] In addition to ENSO, Figure 7 also shows steady warming trends in several regions. In particular, a strong, fairly linear warming trend is visible in the Southern Hemisphere, centered on 40°S. This region accounts for a large portion of the warming in the global average. Warming is also apparent at high latitudes in the Northern Hemisphere due to warming in the North Atlantic as previously documented by Levitus et al. [2000a].

3.2. Subsurface Temperature

[29] Maps of interannual temperature variability were also produced using the difference estimate at depths from the surface to 750 m with 10 m resolution in the vertical. Although the temperature maps are too cumbersome to publish here in their entirety, they do provide some insight into the vertical and horizontal distribution of the warming signal shown in Figure 3. To visualize the warming, the 10-year trend in temperature was calculated for each point in the maps.

[30] Figure 9 shows the slope of the trend zonally averaged as a function of latitude and depth. The most rapid warming occurs in the Northern Hemisphere and is due to warming in the North Atlantic. Although this warming is quite rapid, it contributes only marginally to the global trend because the latitude band centered at 60°N represents a small fraction of the global ocean. Strong warming also occurs between the equator and 20°N at a depth of about 200 m. This warming is colocated with the interannual variability associated with ENSO and illustrated in Figure 6. Most of this warming trend is caused by the particular phases of ENSO sampled over the 10 years of the time series and is not related to decadal or long-term variability. As mentioned above, the early part of the time series contains several years of persistent El Niño conditions. The more recent years saw persistent La Niña conditions, and the large 1997–1998 El Niño occurred near the middle. For the 10-year trend, this results in a La Niña-like warming pattern in the western tropical Pacific, coupled with cooling along the equatorial eastern Pacific. Although this pattern of warming is not associated with a long-term trend, it does provide a convenient means of distinguishing the ENSO related variability in the temperature trend maps.

[31] In contrast, the warming around 40°S appears to be much steadier over the course of the time series, as seen in Figure 7. In addition, this warming extends deeper and is more uniform over the water column than the signal in the tropics. This is also illustrated in Figure 10, which shows maps of temperature trend at depths of 200 m and 750 m. Although the ENSO related signal in the tropics is very strong at 200 m, it is completely absent at 750 m. At depth, the strongest warming occurs around Australia. In addition, there is warming throughout most of the southern Indian and Pacific Oceans, and below the southern tip of Africa. All of these regions contribute to the Southern Hemisphere warming signal in the zonal average. Also at depth, a strong warming signal is present in the far North Atlantic where deep water formation occurs, although as mentioned above, this has fairly small implications for the global average. Finally, although it is worth noting that the estimate shows little warming in the far south near Antarctica, scarcity of subsurface data means that the estimate relies heavily on the altimeter data in that region. In addition, the altimetric correlation coefficients in that region are small (see Figure 2b), making it difficult to reach any strong conclusions about the possibility of warming there. Nevertheless, some warming is apparent south of 45°S, and is consistent with the findings of Gille [2002], which showed rapid warming at depths around 1000 m for these latitudes over the past 40 years.

3.3. Thermosteric Expansion

[32] As with heat content and temperature, the difference estimate was used to compute maps of global, thermosteric sea level anomaly in the upper 750 m of the water column. The maps were then area-averaged to produce the time series shown in Figure 11a. This curve represents the change in global mean sea level due solely to thermal expansion. The 10-year difference from mid-1993 to mid-2003 gives an average rate of 1.6 ± 0.2 mm/yr for thermosteric sea level rise. Figure 11a also shows a recent estimate of total change in global mean sea level as calculated by Leuliette et al. [2004]. This curve has been smoothed with a 1-year boxcar filter for comparison with the thermosteric estimate. For clarity, all three curves were shifted vertically so as to cross zero in mid-1993, the point that represents the average over the first full year of TOPEX/Poseidon data.

[33] Also shown in Figure 11a is the difference between the total and thermosteric sea level estimates. This residual curve potentially contains signals due to changes in ocean salinity (halosteric), changes in ocean mass, or changes in ocean temperature below 750 m. As suggested by Munk [2003], the halosteric component is expected to be small compared to the changes in ocean mass. Deep ocean temperature changes over the past 10 years, however, are very difficult to estimate given the present data. Nevertheless, the work of Levitus et al. [2000a] suggests that most of the warming signal on decadal timescales occurs in the upper 1000 m of the ocean. With these caveats, we may consider the residual curve in Figure 11a to be an inferred estimate of the contribution to sea level rise due to increase in ocean mass, which we refer to as the eustatic contribution. From the Leuliette et al. [2004] curve, the 10-year average rate of total mean sea level rise is 2.4 mm/yr. This implies an average rate of 0.8 mm/yr eustatic sea level rise during the 10 years of the record. In addition, the eustatic curve shows a large peak coincident with the 1997–1998 ENSO event.

[34] Figure 11b again shows the estimate of globally averaged thermosteric sea level rise from the present study, along with another estimate of the same based on yearly maps of global temperature anomaly from Levitus et al. [2000b]. These were the maps used by Cabanes et al. [2001] to produce their estimate of thermosteric sea level rise, and the curve shown in Figure 11b is equivalent to the one published there. The Cabanes et al. [2001] curve shows a much more rapid increase in thermosteric sea level at the onset of the 1997–1998 El Niño event. By comparing the maps of thermosteric sea level rise, it was found that most of the discrepancy between the estimate presented in this study and the Cabanes et al. [2001] estimate occurred in the Indian Ocean. Specifically, the El Niño-related cooling pattern in the eastern tropical Indian Ocean is poorly sampled by the data set used to create the Levitus et al. [2000b] maps (World Ocean Database 1998). This causes the Cabanes et al. [2001] curve to appear artificially warm during these years. In the present study, floats deployed as part of WOCE program provided over 1000 profiles per year in the Indian Ocean, beginning in 1995. These data were not included in the Levitus et al. [2000b] maps. By including them here, along with the altimeter data, sampling in the tropical Indian Ocean is substantially improved. The resulting estimate shows warming and cooling patterns similar to those of the tropical Pacific, but opposite in sign.

[35] Because the warming rate and the rate of long-term thermosteric sea level rise are numbers of particular importance, it is wise to check for the presence of systematic errors in the difference estimates of these quantities. We first consider the in situ data. It is possible for a systematic error to be introduced by the hydrographic data through a misapplication of the XBT fall rate correction. Although the fall rate correction was applied to the best of our knowledge, it is nevertheless important to quantify the magnitude of its effect on globally averaged quantities. The correction recommended by Hanawa et al. [1995] is a linear increase in the depth of a profile by 3.36%. Since temperature generally decreases with depth in the ocean, application of the correction usually results in an overall warming of the profile. Furthermore, most profiles of unknown XBT probe type occur in the first few years of the record. To quantify the effect of correcting the profiles of unknown XBT probe type, the thermosteric sea level curve was recalculated without applying the correction to the unknown probes. The resulting curve was lower by about 7 mm in the first 3 years of the record. This increases the average rate of thermosteric sea level rise to 2.3 mm/yr for the record, suggesting a potential systematic error of 0.6 mm/yr or about 35% of the 10-year rate. This is a very large systematic error, and it is unlikely that the depth correction has been improperly applied to such a large number of profiles. Nevertheless, this calculation underscores the sensitivity of the globally averaged time series to the XBT fall rate correction. Since some care was taken in the application of the fall rate correction in the present analysis, this systematic error was not included in the sea level rise error budget. The calculation, however, serves to illustrate the importance of applying the correction carefully and of disclosing the assumptions used in selecting the profiles that received the correction.

[36] We would also like to test whether the rate of thermosteric sea level rise, as estimated by the difference method, is sensitive to systematic drifts in the rate of total sea level rise as measured by the altimeter. To do this, the slope of global mean sea level rise in the altimeter data was artificially adjusted by adding a spatially uniform trend to the raw AVISO altimeter data. Using this “adjusted” altimeter data set, the difference estimate of thermosteric sea level was then recalculated and compared with the curve shown in Figure 11. The rate of thermosteric sea level rise was found to be insensitive to large changes in the trend of the mean sea level in the altimeter data. The rate of global mean sea level rise in the AVISO altimeter data was 2.4 mm/yr. Doubling this rate to 4.8 mm/yr or removing the slope entirely resulted in rates of thermosteric sea level rise of 1.9 mm/yr or 1.5 mm/yr, respectively. This suggests that enough in situ data exist to largely determine the 10-year rate of thermosteric sea level rise and that the difference estimate is relatively insensitive to any systematic errors that may be present in the altimeter data. Furthermore, this supports the results of Figure 4, which show that in situ data alone resolve most of the 10-year trend in the global average. Using a conservative estimate of 1.5 mm/yr error in the altimetric-based measurement of global mean sea level rise implies a systematic error of 0.1 mm/yr in the difference estimate of thermosteric sea level rise. Adding this to the estimate of random error based on the altimeter data gives the 0.2 mm/yr error in the rate quoted above. Including systematic errors, individual points of the thermosteric sea level time series have errors of about 1.2 mm. As mentioned above, a similar error was included in the difference estimate of heat content. Scaling the error appropriately suggests the systematic error in the warming rate of 0.06 W/m² used in 6.