Pat Frank and Error Propagation in GCMs

This post reviews recent controversy over some strange theories of Pat Frank. It reviews his WUWT posts, blog discussion, some obvious errors, and more substantively, how error propagation really works in solution of numerical partial differential equations (pde, as GCM climate models are), why it is important, and why it is well understood. It's a bit long so I'll start with a TOC.

Contents:

• The story
• Responses
• ad Absurdum
• How error really propagates in solving differential equations
• Eigenvalues
• Summary

The story

There has been a long running soap opera about attempts by Pat Frank to publish a paper which claims that climate modellers are ignoring something he calls "propagation of errors", which he claims would yield extraordinarily large error bars, and invalidate GCM results. These attempts have been unsuccessful, about which he has been becoming increasingly strident. He says the six journals and many referees that have rejected his papers, often with scathing reviews, are motivated only by conflict of interest, desire to preserve their funding etc. Publication of Pat Frank's papers would bring that all down. The most recent outburst at WUWT was here, where a seventh journal promptly rejected it. James Annan and Julie H were involved.

What he means by "propagation of errors" is just the process of forming an expression combining variables with associated error, and expressing this in terms of a total derivative, with the error magnitudes then combining in quadrature (they are assumed independent). That is certainly not something scientists are ignorant of, but it doesn't express what is happening here. There is a lot more involved in differential equation solution.

Responses

Far too much time has been taked up with this nuttiness. First of course by the journals, who have to deal with not only normal reviewing, but a stream of corrosive correspondence. ATTP has had, I think, several threads, the latest here. James Annan expands on his cameo here. And Patrick Brown even put together a video. Further back, similar stuff made the rounds of skeptic sites, eg here. The more mathematical folks there thought it was pretty nutty too.

As reference materials, Pat Frank posted his paper and also a listing of the journal reviews and correspondence here at WUWT. His link to the paper is here. It's not a short paper - 13.4Mb. His link to the zipfile of reviews (44 Mb) is here.

ad Absurdum

As you might expect, I have been joining in at WUWT too. My view, which I'll expand on below, is that propagation of error in numerical partial differential equations is very important, necessarily gets much attention, and is nothing like what he describes. But when folks like Eric Worrall at WUWT are solemnly chiding me for not understanding that propagation is a random walk, I can see expounding the correct theory there would be a waste of time. So I instead try reductio ad absurdum on a few rather obvious points. But the audience at WUWT is not sensitive to absurdum.

In the first recent post I picked up on a criticism of a reviewer that Pat Frank had insisted on a change of units when he averaged quantities over 20 years. He said that the units became quantity/year. Others had raised this, so in the text he had doubled down. Surely this qualifies as absurd?

I again tried a somewhat peripheral approach, saying, why years? Why not regard it as a period of 240 months, and say quantity/month. He said, well the quantity would change. 4 W/m2/year would become 1/3 W/m2/month. But the problem there is that he was using data from someone else's paper (Lauer and Hamilton, 2013) who simply said the average was 4 W/m2. No time interval specified.

It actually relates to the objection many referees raised, including James Annan. If you're going to say the error adds like a random walk, what is the time interval? Why per year? It obviously affects the result how long between steps. He gave some strange replies. Here he expounds the difference between a "measurement average" and a "statistical average":

Several measurements of the height of one person: meters. A measurement average.
Average height of people in a room: meters per person. A statistical average.
Middle school math, and a numerical methods PhD [moi] stumbles over it.

So, what if you sample in Eurpoe to see if Dutchmen are taller than Greeks. You get 1.8 m/Dutch and 1.7 m/Greek. Can you compare if they are in different units? (h/t JA)

On the second thread, I picked up on a claim from PF's review of JA's review:

He wrote, “… ~4W/m^2 error in cloud forcing…” except it is ±4 W/m^2 not Dr. Annan’s positive sign +4 W/m^2. Apparently for Dr. Annan, ± = +.
...
How does it happen that a PhD in mathematics does not understand rms (root-mean-square) and cannot distinguish a “±” from a “+”?

Thius was a less central point, but it throws down the gauntlet. Either JA or PF doesn't have a clue. So I thought, surely there are some there who are familiar with RMS, or the particular case of standard deviation? Surely they know they are always written as positive numbers, with good reason. But to reinforce it, I challenged them to find somewhere, anywhere, where an RMS quantity was written with a ±. But no-one could. PF did come up with a number of locations, which wrote the number as positive, but used it in an expression like a±σ. I couldn't get across the distinction. Yes, 3±2 is an interval, but 2 is a positive number. It isn't JA who is clueless.

The depressing thing to me was that there were people there who seemed to have some acquaintance with statistics or math, and would never before have thought of adding a ± to express a RMS. But now it was all so obvious to them that you should. The same conversion happened with the average units.

How error really propagates in solving differential equations

It's actually a very important practical issue for anyone who solved large systems of differential equations (de's), as I do. Because you always have error, if only from rounding. And if it can grow without limit, the solution is unstable.

A system of partial de's, as with GCMs, can be first discretised in space, to give a large but finite set of nodal values (and perhaps derivatives, or other quantities) with relations between them. There is a simplification that usually only nearby locations are related, which means the resulting set of ordinary de's (in time) are sparse. If the quantities are A (big vector) and the drivers F, the de could be written linearised as
dA/dt = -S*A + F
where S is the big sparse matrix. The time step will need to be taken short enough so that S varies slowly between steps, and it is enough to take it as constant. It can be expressed in terms of eigenvectors. So the driver, which could include an error term, is then also expressed in terms of the eigenvectors, and you can think of the equation in each direction as just being a scalar equation
da/dt = -s*a + f
The solution of this is a = exp(-s*t) ∫ exp(s*τ)*f dτ
where the range of integration is from arbitrary t0 to t.
As you go backward in time from t, the multiplier of f dies away exponentially, if s>0. So a "propagates" f to only a limited extent. It's actually an exponentially smoothed version of f. S may not be symmetric, so the eigenvalues could be complex, and the criterion then is real part Re(s)>0.

The same thing persists after you discretise in time, so integration is replaced by summation in discrete steps. This is actually where care is needed to ensure that a stable de converts to a stable recurrence. But that can be done.

So to prevent instability, it must certainly be true that eigenvalues s≥0. The interesting cases are when Re(s) is small. I'll look at the possible spectrum.

Eigenvalues

In the recurrences that result from a GCM, there are millions of eigenvalues. The salvation is that most of them have relatively large real part, and correspond to rapidly decaying processes (typically, hours for GCM). They are all the local flows that could be generated on the element scale, which lose kinetic energy rapidly through shearing. That is why the solution process can make some sense, because one can focus on a much smaller eigenspace where decay is either absent or slow on the scale, perhaps decades, of the solution being sought. I'll look at three classes:

• Waves. Here s has zero (or near) real part, but large imaginary part, so decay is slow and period is comparable to the time step. In classical CFD, these would be sound waves. They are important, because they create action at a distance, propagating the pressure field so the medium behaves like a not very compressible fluid (elasticity). In GCM's, they merge into gravity waves, because a pressure peak can move the air vertically. That actually reduces the propagation velocity. The waves do have some dissipation (eg surface shear) so the ideal s has small positive real part.
  
  These waves are important, because to the extent they aren't well resolved in time after time discretisation, s can have negative real part, and if it outweighs dissipation, the wave will grow and the solution fail. That is why GCMs have a maximum time step of about half an hour, which demands powerful computers.
• Secular motion and slow decay. The secular part is the desired solution responding to the driver F; I've bundled in the slow decay processes because they are hard to separate. The usual way is to wind back, so that only the slowest remain from the initial perturbation. But some does remain, and that is why GCMs tend to have apparently similar solutions which may differ by a degree or two in temperature, say. There are two remedies for this
  • use of anomalies, which subtracts out the very slow decay, leaving he secular variation
  • Increased care with initial conditions. The smaller the component of these modes that is present initially, the smaller the later problem.
• Conserved quantities - energy and mass (of each component). These totalled correspond to zero eigenvalues; they shouldn't change. But in this case small changes do indeed accumulate, and that is a problem. GCMs include energy and mass "fixers" (see CAM 3 for examples). The totals are monitored, discrepancies noted, and then corrected by adding in distributed amounts as required.
  
  I remember a post at WUWT where Willis Eschenbach discovered this and cried malpractice. But it is perfectly legitimate. GCMs are models which incorporate all the math that the solution should conform to, and conservation is one such. It is true that as an extra constraint it appears to make more equations than variables. But large systems of equations are always ill-conditioned in some way; there aren't as many effective equations as it seems. That is why the conservations were able to erode, and the fixers just restore that. "Fixing" creates errors in the first class of eigenvalues, because while care should be taken to do it without introducing perturbations in their space, this can't be done perfectly. But again the salvation is that such errors rapidly decay.

Summary

So what do we have with error propagation, and why is Pat Frank so wrong?

• Inhomogeneous equation systems are driven by a mix of intended drivers and error.
• The effect of each is generally subject to exponential decay. The resulting solution is not a growing function of past errors, but in effect, like an exponential smooth of them. This may increase by a steady factor, but will also attenuate noisy error.
• Propagation is not something that is ignored, as Pat Frank claims, but is a central and much studied part of numerical pde.