Sea level data from Church appear be integrated as I(1).

d Root ADF Padf
[1,] 0 0.9713052 -0.8354583 0.9561317
[2,] 1 -0.2771277 -5.8808801 0.0100000
[3,] 2 -1.1410606 -8.1287823 0.0100000

As does Jevrejeva’s data set from 1700.

d Root ADF Padf
[1,] 0 0.7552908 -2.106932 0.5312376
[2,] 1 -0.4415736 -9.329505 0.0100000
[3,] 2 -1.3634252 -12.083777 0.0100000

And while the correlation is high when sea level is added into the linear model, the sea level almost blocks out all the other variables:

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.2567157 0.0258570 -9.928 < 2e-16 ***
Solar 0.3310376 0.2199565 1.505 0.1350
StratAer 0.0236418 0.0232146 1.018 0.3106
deltaAGW 1.0072159 1.1876318 0.848 0.3981
Residuals 1.0680420 0.4129881 2.586 0.0109 *
SeaLevel 0.0024081 0.0004251 5.665 1.07e-07 ***
—
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Multiple R-squared: 0.7067, p-value: < 2.2e-16

And according to the ADF test they don't cointegrate (ie the residuals are not I(0).

d Root ADF Padf
[1,] 0 0.5935631 -2.530107 0.3563188
[2,] 1 -0.3420658 -6.627341 0.0100000
[3,] 2 -1.1384317 -9.886875 0.0100000

However, the pp.test disagrees:

Phillips-Perron Unit Root Test on l$residuals
Dickey-Fuller Z(alpha) = -31.0654, Truncation lag parameter = 4, p-value = 0.01
alternative hypothesis: stationary

Here is a plot of the residuals of a linear model of temperature and sea level. It does also look stationary, which it would be if it were a 60 year cycle. That is the requirement for variables to cointegrate, for the residuals to be stationary, though its based on the AR(1) coefficient and not the trend.

A difference between tests was also noted in Beenstock’s paper on testing the order of integration of CO2.

Test 1 shows that according to all three test statistics rfCO2 is not trend stationary. Test 2 shows that according to the PP statistic rfCO2 is marginally difference stationary, but the KPSS and ADF statistics clearly reject this hypothesis.

The upshot is that sea level is well correlated with global temperature and essentially moves in lock-step, at least over this time span. Temperature and sea level convey effectively the same information, apart for a possible 60 year cycle. It is not informative to include them in a linear model with other variables. Which underlines the need to understand what linear regression and the tests are doing, and not apply them blindly. Global temperature plus sea level (blue below).

The cointegration of sea level an temperature also means that there is no evidence in the current data to support the view that sea level will ‘detach’ from temperature and rise non-linearly to high levels. A similar rise in temperature this century, as occurred last century, will result in a rise in sea level of the same amount, around 20cm.

While I think a proxy for ocean heat should have been included in the Beenstock analysis, as the ocean is such an important part of the climate system, as it could have caused an upset, it doesn’t appear to contradict their analysis.

The abstract of the Phillips Perron test for a unit root is below.

This paper proposes new tests for detecting the presence of a unit root in quite general time series models. Our approach is nonparametric with respect to nuisance parameters and thereby allows for a very wide class of weakly dependent and possibly heterogeneously distributed data. The tests accommodate models with a fitted drift and a time trend so that they may be used to discriminate between unit root nonstationarity and stationarity about a deterministic trend. The limiting distributions of the statistics are obtained under both the unit root null and a sequence of local alternatives. The latter noncentral distribution theory yields local asymptotic power functions for the tests and facilitates comparisons with alternative procedures due to Dickey & Fuller. Simulations are reported on the performance of the new tests in finite samples.