Bart: Cloud feedback is negative – ocean response is around 4.88 years

I’m not sufficiently up on stats to really understand this, but I’m collating the relevant comments from Climate Audit and Roy Spencers site here, because it looks important, and moderation at CA seems to be impeding the flow of the conversation so it has become disjointed.

Is there anyone with a voice in these matters who could just do a cross spectral estimation of the two series and read off the phase relationship directly?

You cannot diagnose feedback by performing a linear regression on a phase plane plot when the driving input is all over the place and the phase response is nonlinear. As I explained on the last thread, and have written up here.

What I have found is significantly greater evidence that the feedback relationship is, indeed, negative, contrary to Dessler. The phase shift at low frequencies, which determines the sign of the feedback, is very clearly near 180 degrees. Saying an input is 180 degrees from the output is a long winded way of saying that the one is the negative of the other…I’m sitting back and watching people argue about techniques which are entirely unsuited to the problem in the first place. It’s insane.

 I need other knowledgeable people to become aware of the result and determine what the effect of a -9.5 W/m^2/degC feedback would be. I think, though, that the IPCC models require the overall feedback to be positive to get any significant warming, so it appears to me that this result could, potentially, kill CAGW.

There could be a valid criticism that the span of data is too short to have high confidence in the result. But, I would argue that the onus is on them to prove that the overall feedback is positive, because I think this establishes that the running assumption should be that it is likely negative.

.

If you have an input x(t) and an output y(t), the cross spectrum is the expected value of X’(f)*Y(f), where X(f) and Y(f) are the Fourier transforms and “‘” denotes the complex conjugate. Since Y(f) = X(f)*H(f), where H(f) is the transfer function, the cross spectrum is C(f) = P(f)*H(f), where P(f) is the power spectral density of x(t). Since P(f) is a zero phase, positive real function, the phase response of C(f) is identical to that of H(f). If it is near 180 degrees at low frequency, then you’ve got negative feedback. If near zero, it is positive.

The dispersion in the phase plane makes it very difficult to diagnose the sign of the feedback – the nonlinear phase characteristic of H(f) ensures that you will get a messy cloud unless you focus on a particular frequency range and filter everything else out. For sure, the cross spectrum will be polluted with other processes and noise, too, but I expect it should be a pretty good indicator of the dominant process occurring, if there is one.

Caveat: estimation of the cross spectrum is not as easy as multiplying the FFTs of the two data streams. The outcome of that will be highly erratic and of questionable value. A good method is to do that, then inverse transform to get the impulse response. Window that impulse response with a good window (e.g., Hamming or other) at a truncated limit before it visibly starts breaking up, then transform the windowed impulse response function back to get the estimate of C(f). There are, of course, other reasonable methods. A mathematical analysis tool like MATLAB has its own built-in functions for doing it.

Well, here are some preliminary results. I hope I used the right quantities. I pulled this data from the Spencer link. I used column 9 for the HADCRUT3 temperature anomaly, and I used column 5 minus column 8 for the cloud response. The relationship between these variables most assuredly has a negative dc gain.

Here is a plot of the estimated frequency response. If this holds up, I think it’s going to shock some people. The response at high frequency is a jumble, and probably due to independent processes going on. But, the low frequency region is dominated by a fairly well defined 2nd order response with natural frequency of about 0.0725 year^-1 and a damping ratio of about 0.45, which indicates a time constant of about 4.88 years. Yes, years.

The time constant is 4.88 years, meaning the step response is 63% of the way to the final value in that time. Convergence is exponential in the standard, usual way.

The impulse response is:

I am talking about the time it takes for clouds to react to temperature changes. That is, if you increase global temperature by 1 degC (and could hold it there), within 4.88 years, you will be 1-exp(-1) = 63% of the way to creating an opposing 9.5 to 10 W/m^2 reduction in insolation. The unit step response is plotted below

I have posted the transfer function and impulse response estimates here and here, respectively, created using artificially generated data according to the prescription here.

These show it is, indeed, possible to pull out the long term correlation using only 124 monthly data points. It doesn’t always work, but sometimes, it does.

It would be very nice to come up with a more robust deconvolution scheme which works the majority of the time, and reanalyze the data to prove the veracity of the estimated functions beyond doubt, and I encourage any interested parties to look into this.

Here is the MATLAB code used to compute all this:

% Import monthly data
temp = data(:,9);
dR = data(:,5)-data(:,8);
N = length(dR);

% Sample period
T = 1/12; % years

% Pad time series with zeroes to prevent time aliasing of impulse response
Nsamp = 8192;
Npad = Nsamp-N;
X = fft([temp;zeros(Npad,1)]);
Y = fft([dR;zeros(Npad,1)]);

% Compute impulse Response
h = real(ifft(Y./X))/T;

% Window Impulse Response
Nc = Nsamp/2^2;
w = [ones(Nc/2,1);(1 + cos(pi*(0:(Nc/2-1))’/(Nc/2-1)))/2];
w = [w;zeros(Nsamp-Nc,1)];
hw = h.*w;

% Plot smoothed impulse response
c = [1:15 15:-1:1]/(15*16);
figure(1)
hs = flipud(filter(c,1,flipud(h)));
t = (0:(length(hs)-1))*T;
plot(t,hs)
grid
xlim([0 50])
title(‘Cloud-Temperature System Smoothed Impulse Response’)
xlabel(‘time (years)’)
ylabel(‘W/m^2/^oC/year’)

% Compute Frequency Response and plot
H = T*fft(hw);
f = (0:(length(H)-1))’/Nsamp/T;

% Create 2nd order model
s = sqrt(-1)*2*pi*f;
w0 = 2*pi*0.0725;
zeta=0.45;
Hmod = -9.5./((s/w0+2*zeta).*(s/w0)+1);

figure(2)
subplot(211)
loglog(f,abs([H Hmod]),’LineWidth’,2)
grid
xlim([1e-3 5e-1])
title(‘Cloud-Temperature System Magnitude Response’)
ylabel(‘W/m^2/^oC’)
xlabel(‘frequency (years^-^1)’)
legend(‘From Data’,’Model’,’Location’,’SouthWest’)
subplot(212)
semilogx(f,faz([H Hmod]),’LineWidth’,2)
grid
xlim([1e-3 5e-1])
title(‘Cloud-Temperature System Phase Response’)
ylabel(‘deg’)
xlabel(‘frequency (years^-^1)’)

• Nick Stokes Posted Sep 10, 2011 at 5:32 PM | Permalink

  Bart,
  I translated your code to R, up to the impulse response plot. I get the same result.

  But if I apply a Hanning taper to temp and dR (down to zero at each end of the data window) it is very different. I think that suggests that the finite data length is having a big effect.

  I’ve put up a html page here which has some explanation, the impulse response with and without Hanning, and the R code.

• Bart
  Posted Sep 10, 2011 at 8:23 PM | Permalink

  Nick – yes, it can change things a little through loss of resolution, but it isn’t really all that significant and, importantly, the feedback is still negative (180 deg phase shift at low frequencies).

  I have other responses which may be of interest to you and which are directed to your queries but, unfortunately, the moderator seems to be taking the weekend off. Your insights are cogent, but I think I feel relatively confident now in asserting that this response is real, and the feedback is, clearly, negative. But, please, continue to investigate.

• Bart
  Posted Sep 10, 2011 at 8:32 PM | Permalink

  “You have to accept low frequency problems with finite data…”

  Actually, Nick, abrupt transitions are a high frequency phenomenon. Multiplication in the time domain is the same as convolution in the frequency domain, so the tapering tends to lower the resolution of the frequency response estimate. Which is good for the higher frequency portion of the spectrum, but not so good for the low. Yes, a truncated data record is bad for resolution of low frequency stuff. But, providing a taper only really helps the high frequency portion. And, so, your impulse response so estimated has the slight ringing smoothed out.

• Bart
  Posted Sep 10, 2011 at 8:58 PM | Permalink

  I’m going to try one more time to get through the blocker on this. Here is a discrete time model for generating data with the desired correlation which you can play with. I assume the input is white noise. You should find that, occasionally, it is possible to generate a short time series for which the analysis works, other times, not. Which is why I said: “My preliminary conclusion right now is that the Spencer data just happens to be lucky.”

  a = [1.000000000000000 -1.967462610776618 0.968691947164695];
  b = -[0.617926899846966 0.611409488230977]*1e-2;
  temp=randn(10000,1);
  dR = filter(b,a,temp);
  temp = temp((10000-123):10000);
  dR = dR((10000-123):10000);

  I generate a lot of data before truncating to be sure of eliminating transient start up effects.

  Bart says:
  September 11, 2011 at 2:40 am (Edit)

  neill says:
  September 10, 2011 at 9:33 pm

  “Bart, if you would be so kind as to translate these findings for the less mathy/scientific among us, it would be much appreciated.”

  Basically, clouds exert a negative feedback of about -9.5 W/m^2 per degree C of increased temperature (and, vice versa, +9.5 W/m^2 per degree C of decreased temperature). Which means they contribute at least partly to a climate thermostat mechanism such as Willis Eschenbach has proposed here at WUWT from time to time.