Averaging Planet Earths or Averaging Planet Xs
Update January 23, 2011. I fixed the broken links.
There are several discussions floating around on the subject of comparing GCM-calculated numbers with experimental data. Climate Audit and Lucia both have several threads. There are too many threads to give links; let me know if you need a specific thread. One focus of these discussions is the ensemble-average approach that is considered to be necessary for the comparisons. The ensemble is made up of the results from the various different versions of GCMs that calculate the results. Not only are GCMs different from each other, but the suggested approach is to make perturbations in the initial conditions between the runs.
I have mentioned that the results do not reflect the effects of perturbations in only the initial conditions. Everything about the approach is different between the runs. The GCMs are based on different continuous-equation models, numerical solution methods, application procedures, run-time options, and users, and maybe other aspects. How can we be sure that the results are all appropriate for Earth?
How do we know we are averaging different results for Earth as opposed to averaging different planets. I have suggested that some kind of filtering mechanism is needed in order to eliminate the calculated numbers that might not be Earth. Why, for example, can’t a model that gives a constant value for the Global Mean Surface Temperature be used in the ensemble.
It is my impression that there are no filtering mechanisms in place. The only requirement for being a part of the ensemble is that your code has calculated the numbers. All the available numbers get tossed into the ensemble. A very rough approach to validation of models, in my opinion.
The ensemble average approach is used, as I understand it, because of the non-deterministic nature of the models / methods / codes / applications / users. Additionally, this heuristic, ad hoc, osmosis-like hypothesis is generally based on the chaotic response observed with the original Lorenz equation system of 1963. Those equations exhibit chaotic response for some ranges of the values of the parameters in the system. These parameters are roughly the Rayleigh and Prandtl Numbers. And while most analyses of the system focus on perturbations in the initial conditions, it is well known that changes in the parameters and the use of different step sizes in the numerical integration also exhibit chaotic response. The Lorenz systems, in general, are ill-posed initial value problems.
For this discussion let’s take chaotic response to mean that given small perturbations, the solution trajectories do not approach fixed equilibrium state points and do not exhibit periodic response; aperiodic response is obtained. It is equally well known that for certain ranges of the parameters, the original Lorenz equation system exhibits well-behaved response and the dependent variables attain fixed equilibrium state points. And, small changes in the coefficients of some of the terms in the original system gives systems that rapidly approach fixed equilibrium states. Finally, we note that while the parameters have constant values in the Lorenz system, the corresponding quantities in the GCMs will vary during the course of the calculations.
So, now back to the original question. Given the wide range of observed responses how do we know that we’ve included only different versions of Earth into the ensemble. How do we know we have not introduced different planets into the ensemble?
I have discussed the original Lorenz system and associated numerical solution methods in several posts on this blog. The equation system is summarized here, the numerical solution methods are summarized here, and some results are summarized here. For this post I have used the explicit Euler method although I’m certain that any other of the methods would produce the same general results. I used the same step size for the independent variable for all the runs discussed here.I have made several runs as follows:
(1) A base-case with selected values of the initial condition and parameters.
(2) A run with a small perturbation in the initial values for the first dependent variable and all other held at the base-case values.
(3) A run with a small perturbation in the initial values for the second dependent variable and all other held at the base-case values.
(4) A run with a small perturbation in the initial values for the third dependent variable and all other held at the base-case values.
At this point I made an ensemble and calculated the average values of the dependent variables. I will get to some details of the results below.
I then made additional runs so as to include a few more versions of Earth into the ensemble.
(5) A run with the base-case values of the initial conditions and the Prandtl Number, but with a change in the Rayleigh Number.
(6) A run with the base-case values of the initial conditions and the Rayleigh Number but with a change in the Prandtl Number.
(7) A run with the base-case values of the initial conditions and changes in both the Prandtl and Rayleigh Numbers.
At this point I made another ensemble including the seven cases. Note that there are a very large number of combinations of changes that could be made, and if we include the magnitude of the changes as an additional parameter, the total number of possibilities grows without bound.
I have up-loaded a plot of the results here. If you successfully get the plot, here’s what you’ll see.
The bottom of the plot shows the results of the 4-run ensemble for the Y(3) dependent variable. All four runs are shown plus the average value, which is the thick black line running through the cloud of results. The top part shows the average value for Y(3) as given by the 7-run ensemble. I didn’t include the component results because the cloud would be even more cloudy.
Well, what has happened? It looks like I’ve picked up a Planet X somewhere in my calculations for Planet Earth. And that Planet dominates the calculated results. Full disclosure; I have not yet looked at any of the other dependent variables in this detail. Looking at the results for Y(3) I see that the run for which the value of both parameters were changed has made Planet X.
This was my zeroth-order cut at trying to find an answer to the question, How do we know we are averaging different results for Earth as opposed to averaging different planets. I think it provides some validity to the issue.
I have suggested that some kind of filtering mechanism is needed in order to eliminate the calculated numbers that might not be Earth.
Models Methods Software 
Dan I don’t see clearly why you talk about Earth as opposed to a planet X .
Everything you did are Earths .
Why ?
Because the only thing that didn’t change in the 7 scenarios are the equations .
Initial conditions and constant parameters may have changed but the equations didn’t .
As those equations represent the invariant natural laws governing the Earth (in this case a convection case) , they describe correctly the Earth in every case .
Of course as it is well known , the Lorenz system (and its physical equivalent a convecting fluid) exhibits deterministic chaos for certains values of the constant parameters also called control parameters.
So once you have selected the value of those parameters so that the system is in a chaotic regime , you will observe the necessary and obvious consequences as follows .
1) The trajectories of Y(3) will diverge exponentially . That means that for a small difference in the initial conditions (= Y1(3,0) – Y2(3,0)) the distance between the 2 curves Y1(3,t) – Y2(3,t) will evolve like exp(L.t)where L is the relevant Lyapounov coefficient for Y(3) .
2) As Y(3) is bounded , this exponential divergence will not go forever but will be bounded by the greatest size of the attractor along the Y(3) dimension . The speed at which it diverges depends on L which is only dependent on the value of the control parameters .
3) In your scenarios where only the initial conditions vary you have Earths with same and constant Lyapounov coefficients . You have also the same attractor . In the scenarios where you changed the control parameters , the Lyapounov coefficients changed and the attractor did so too . These Earths are still Earths because the natural laws didn’t change but they are Earths with a different amount of chaos in them .
Allow me an analogy .
The first scenarios correspond to a river whose eddies you observe and then you put rocks at different places in the flow . The eddies will change according to the size and placement of the rocks but it is the same river that vehiculates approximately the same energy in all rock placing scenarios .
The second scenarios correspond to the same river bed but you significantly increase (or decrease) the water speed . Here the impact on the form , number and location of the eddies will be important even if you put no rocks in the flow . It is not really the same river as the first one and it vehiculates a different amount of energy but it is still very “riverish” all the same .
4) And now … you make … err … an average of all that …. 🙂
I will oversimplify but it will give a gist of what happens .
Let’s choose an arbitrary reference scenario Y0(3,t) .
Then let’s choose N scenarios where initial conditions vary and give Yi(3,t) and M scenarios where the control parameters vary and give Yj(3,t) .
The … umm … average is Yav(3,t) = [Sum over N (Yi) + Sum over M (Yj)]/(N+M) .
But we know that Yi – Y0 = Ai.exp(L.t) and Yj – Y0 = Bi.exp(Li.t) .
So Yav = Y0 + [Sum over N (Ai.exp(L.t)) + Sum over M (Bi.exp(Li.t))]/(N+M)
That’s for the first time steps .
Then the term Ymax with the largest Lmax attains the bound which is the maximum size of the attractor (it means physically that Ymax and Y0 are in completely different places of the attractor and have long forgotten anything about their initial similitudes) and Ymax – Y0 stops growing exponentially and begins to do whatever it does – decreasing , increasing , plateauing .
It is soon joined by the second largest Li term a bit later .
Etc .
Yav has lost any meaning and any “relationship” with Y0 but as it had none in the first place , we are not that surprised .
And now imagine that in reality the control parameters are not constant but vary with time thus defining a function L(t) and no more a finite , discrete spectrum of Li .
Arrrgh !
Question for you Dan. If a climate model uses “noise” to reconcile projected temperature trends to observed temperature trends, how much much of this noise can or should be traced back to uncertainties in forcings and/or feedbacks? Would uncertainties on the “input” side of the equation place constraints or limits on the amount of noise which can be attributed on the output side?
Layman Lurker, I have not yet even started looking at these issues. I’m still stuck back at the point that the codes and calculations require Independent Formal Verification before the calculated numbers can be assumed to accurately and correctly represent the model equations and numerical solution methods.
Until it is determined that the ‘noise’ is indeed a proper representation of the model equations, the large potential that it is in fact actually noise and unrelated to any physical phenomena and processes cannot be dismissed.
Generally it has been my experience that the wiggles produced by codes are trying to tell us something and cannot be simply assigned to be ‘noise’ associated with physical phenomena and processes.