The Oil ConunDRUM

Wind speeds of the World

2012-02-10T22:25:00.000-08:00

See if you can follow this argument for the neatest little derivation of a wind speed estimator. I believe that this thing can model the distribution of wind speeds over the entire planet.

In The Oil Conundrum, we modeled the wind speed for a single location over a span of time. This turned into a simple maximum entropy estimator, which assumed a kinetic energy for a wind volume

$$ E \sim v^2 $$

and then a probability distribution function which assumes an average energy over a region i:

$$ P(E \mid E_i) = e^{\frac{-E}{E_i}} $$

but we know that the Ei can vary from region to region, so we leave it as a conditional, and then set that as a maximal entropy estimator as well

$$ p(E_i) = \alpha e^{- \alpha E_i} $$

then we integrate over the conditional's range according to standard practice.

$$ P(E) = \int_0^\infty P(E \mid E_i) p(E_i) \,dE_i $$

This results in a simple lookup in your favorite comprehensive table of cataloged integration formulas, which leads to the following solution:

$$ P(E) = 2 \sqrt{\frac{E}{\bar{E}}} K_1(2 \sqrt{\frac{E}{\bar{E}}})$$

where K1 is the modified BesselK function of the second kind, in this case of order 1, which is found in any spreadsheet program (such as Excel).

First of all, note that this is the same function that we used in TOC to determine the distribution of terrain slopes for the entire planet as well! There, I took a slightly different route and derived it based on a cumulative probability distribution function (CDF) instead of a probability density function (PDF). No matter, in both cases it reduces to the simplest model one can imagine -- a function with only one adjustable parameter, that of the mean energy of a sampled wind measurement. In other words the expression contains a single number that represents an average planetary wind speed. That's all we need.

So let us see how it works.

I used wind data from Bonneville Power Administration, which has over 20 meteorological stations set up around northern Oregon. The download consisted of over 2.5 million data points collected at 5 minute intervals, archived over the span of a little less than 2 years.

Fig. 1: Cumulative distribution function of wind energies from Bonneville with model fit.

I actually didn't have to even vary the parameter, as the average energy was derived directly by computing the mean over the complete set of data separately. This corresponded to a value of 12.4 MPH, and I placed a pair of positive and negative tolerances to give an idea of the sensitivity of the fit.

As this is a single parameter model, the only leeway we have is in shifting the curve horizontally along the energy axis, and since this is locked by an average, the fit becomes essentially automatic with no room for argument. The probabilities are automatically normalized.

I also show the log-log plot next, which reveals a departure at high wind speeds.

Fig. 2: Cumulative distribution function of wind energies on a log-log plot.
At large wind speeds, the fit diverges from the data as the peak 5-minute wind speed observed was 60 MPH.

I am not too worried about this slight divergence as it only shows that excessive gale force winds (greater than 60 MPH) did not occur over the extended region during a span of two years. The next logical step is to average this over longer times and larger spatial regions.

I have to admit that this investigation into the Bonneville dataset was prompted by a provocation and a dare on another blog. In my opinion, anti-renewable-energy advocates seem to crawl out of the woodwork, and tend to spread their misery however they can. This commenter was named P.E. (possibly the professional engineer moniker), and what set me off was his assertion that wind "power output is essentially random":

P.E. | December 30, 2011 at 5:39 pm
I did a little research on the Bonneville system a year or so ago, because they have all their SCADA data available for several years online with excellent resolution. The results were downright shocking. And not in a good way. And this is real live system-wide data, not modeled results. Anyway, be sure to check Bonneville’s website out; it’s low hanging fruit.
Bonneville’s own words: “power output is essentially random”.
David Wojick
Predictability is not the issue. Even if we knew just when the wind would not blow (about 70% of the time) we would still have to have a fueled power plant to produce the needed power, either that or some humongous expensive power storage systems. Renewables reduce emissions by idling existing fueled power plants. That is the most they can do. They are not an alternative to fueled power plants.
P.E.
In the specific case of Bonneville, they have enough hydro in the system to where reserve/backup isn’t an issue. The windmills allow them to use less water when it’s available, making more hydro power available when it’s needed. That’s the ideal situation for wind (and solar). But not that many places in the world can do that.
WebHubTelescope
PE, You are so full of it as to be embarrassing. I have done statistical analysis of wind energy and it follows maximum entropy principles. You can read all about it and more to come, as I seem to have the knack for environmental modeling. Too bad you missed the boat on this one.
P.E.
Web, go do the analysis yourself, and then come back after you’ve looked at some real data. They seem to have changed their web site around, but here’s a place to start: http://www.intellectualtakeout.org/library/chart-graph/load-and-wind-generation-are-essentially-random?library_node=69696
And this: http://pjmedia.com/blog/electric-grid-myths-part-ii-the-effect-of-alternatives/
There used to be some huge excel spreadsheets as evidenced by the two articles above (that my computer would choke on), but they’re either moved to another part of the site, or no longer available to the public.
Maybe I’ll do a FOIA.
P.E.
Correction: they’re right here:
http://transmission.bpa.gov/Business/Operations/Wind/default.aspx
Item #5. Now go to work, or eat your words.
P.E.
Oh, and Web – if you need any help with Excel, I hear Phil Jones is really good at it.
WebHubTelescope
P.E., That’s the problem with your ilk. You are basically ashamed about not being able to do the analysis yourself so lash out at someone like me that knows something about probability and statistics. It is all pathetic, this transparently petty jealousy.
hunter
WHT, who do we believe? Someone who cites a third party that is actually running a windfarm, or a self declared genius like you who has not bothered to challenge the citation, but blames the person who referred to it?
WebHubTelescope
Hunter, The problem is that he shows a spreadsheet with raw data of a few weeks. I have analyzed data from wind farms, compiling the data over the span of a year in terms of a probability distribution. I have it right here on my blog
http://mobjectivist.blogspot.com/2010/05/wind-energy-dispersion-analysis.html
And it is also in my recent book which you can click on my handle to read for free.
You have to keep up with renewable research, otherwise you will be left behind.
hunter
Web, wind is crap , delivering z small fraction of its rated capacity.
P.E.
For the record, that’s not a few weeks, that’s five years. And just for the edification of the ignorant, Bonneville Power Authority is the largest federal hydropower agency in the US, including the Grand Coulee dam at something like 7 gigawatts. This isn’t some small-potatoes outfit.
But keep kvetching.
WebHubTelescope
Kiddies all jealous over the fact that I know how to do the analysis.

The larger point is that the randomness can be modeled, which the garden-variety dogmatists don't seem to understand. This is really no different than understanding black-body Planck's response for radiative emission. Sure, the radiation appears random, but it also must still follow a distribution function with respect to wave-number, which is proportional to energy and thus follows statistical mechanical laws. The same thing holds for wind energy distributions, and we are seeing statistical mechanics, superstatistics, and maximum entropy principles hold sway.

Bottomline is that these environmental skeptics are renowned for all talk and no action when it comes to discussing energy and climate. I admit that I was needling the guy, but that's what prompts them to spill their true agenda, which is hatred against climate scientists coupled with a strange mix of Malthusian, Luddite, and Cornucopian ideals. They actually lash out angrily when they see somebody trying to scientifically advance the yardstick.

A few days later and in another thread, P.E. wouldn't let go, and went out on a limb:

P.E. | January 4, 2012 at 11:30 am
This scold coming from someone who refuses to download wind data from the Bonneville Power Authority because it’s too much like work.
WebHubTelescope
What a wanker, P.E.
I downloaded data from Ontario and Germany wind farms about a year and a half ago and did a comprehensive analysis of the time series data. I happen to be on vacation over the holiday, so am not going to bend over backwards to work the Bonneville data at your whim and call.
I tell you what; I will analyze the Bonneville data set when I get back, and I will add it to a future revision of my tome.
But then again, maybe not, because trained seal and organ-grinder monkey is not my forte.

As in Judo, we learn to use our opponents force and redirect it back at them. They have no intellectual curiosity, and think that just by asserting something it must be true. Only by applying models to empirical data can we make headway in our understanding of our environment. The last laugh is on P.E. and I have to thank him for providing motivation to do this analysis.

The pushback I see against this model from the skeptical trenches is that it doesn't predict chaotic disturbances. I view chaotic bifurcation as a no-nothing theory. We can't predict anything using it, but know that an external forcing will always overpower the minor chaotic fluctuations that may occur. In other words, chaos is annoying noise in the bigger scheme of things.

The only caveat is whether that noise will also scale with the external forcing -- in other words the fluctuations growing bigger with added heat to the system. This is understood with respect to added heat raising the kinetic energy of the system - (1) increasing the average atmospheric wind speed, and (2) added heat to the ocean allowing for greater occasional releases into the atmosphere. These are likely subtle effects, but they will only go in that direction. Remember, reducing heat leads to lower thermal excitation.

The BesselK wind model has a thermal excitation built-in, as the maximum entropy estimator for average energy is nothing but a thermalized Boltzmann partition function describing the system's entropy.

Entropy measures the thermal excitation of the system and that scales with temperature. Non-linearities and chaos is just a way for the disordered system to explore the ergodic range of states, thus making the statistical mechanical/information theory valid. Look at how a naive pseudo-random number generator is created -- by sticking a handful of non-linear functions together, you will explore the state space (see the topic of iterated maps). Voila, it will approach an ergodic limit. Then add in energy constraints, such as the mean energy of the system and you have the Boltzmann estimator, aka MaxEnt estimation.

See also the topic of K-distribution model which models the disorder, or clutter, in electromagnetic interference. This also uses the BesselK model.

Another interesting observation. In the semiconductor world, the transport regimes analogous to the stochastic and chaotic/deterministic talked about in climate circles are diffusive and ballistic transport. It is very hard to get to ballistic transport because the carriers thermalize so rapidly, which is the essential dissipation route to diffusive transport. The idea of stochastic resonance is all over the place in semiconductor physics and technology, with carrier population inversions, bandgap engineering, etc. but we don't talk about it as stochastic resonance because there is already a well established terminology.

On top of that, there is the concept of dispersive transport, which is diffusive transport with huge amounts of disorder in the mobilities and bandgaps. This is the concept that I have done some recent research on because all the cheap photovoltaic technology is built on amorphous or polycrystalline garbage.

On that level, what I think is happening in the climate system is closer to a mix of dispersion and diffusion than determinism.

For example, relating to wind speed distribution, we need to get a handle on what controls the dispersion. From other posts that I have done on wind speed, remembered that the wind power dissipation is the cube of the wind speed, while the distribution that has a skew (the third statistical moment) equal to the mean is just the BesselK distribution that matches the empirical results. This is to me is somewhat mysterious ... is it just that the mean wind speed matches statistically the dissipation due to friction? Could that be some universal scaling result? In terms of maximum entropy it would be the distribution that has a constraint of a finite mean with the mean equating to cube root of the third moment (i.e Mean = standardDeviation * Skew). This then determines all the fluctuations observed over the larger ensemble.

Thermal Diffusion and the Missing Heat

2012-01-23T06:07:00.000-08:00

I have this documented already (in The Oil Conundrum) but let me put a new spin on it. What I will do is solve the heat equation with initial conditions and boundary conditions for a simple experiment. And then I will add two dimensions of Maximum Entropy priors.

The situation is measuring the temperature of a buried sensor situated at some distance below the surface after an impulse of thermal energy is applied. The physics solution to this problem is the heat kernel function which is the impulse response or Green's function for that variation of the master equation. This is pure diffusion with no convection involved (heat is not sensitive to fields, gravity or electrical, so no convection).

However the diffusion coefficient involved in the solution is not known to any degree of precision. The earthen material that the heat is diffusing through is heterogeneously disordered, and all we can really guess at that it has a mean value for the diffusion coefficient. By inferring through the maximum entropy principle, we can say that the diffusion coefficient has a PDF that is exponentially distributed with a mean value D.

We then work the original heat equation solution with this smeared version of D, and then the kernel simplifies to a exp() solution.
$$ {1\over{2\sqrt{Dt}}}e^{-x/\sqrt{Dt}} $$
But we also don't know the value of x that well and have uncertainty in its value. If we give a Maximum Entropy uncertainty in that value, then the solution simpilfies to
$$ {1\over2}{1\over{x_0+\sqrt{Dt}}} $$
where x0 is a smeared value for x.

This is a valid approximation to the solution of this particular problem and the following Figure 1 is a fit to experimental data. There are two parameters to the model, an asymptotic value that is used to extrapolate a steady state value based on the initial thermal impulse and the smearing value which generates the red line. The slightly noisy blue line is the data, and one can note the good agreement.

Figure 1: Fit of thermal dispersive diffusion model (red) to a heat impulse response (blue).

Notice the long tail on the model fit. The far field response in this case is the probability complement of the near field impulse response. In other words, what diffuses away from the source will show up at the adjacent target. By treating the system as two slabs in this way, we can give it an intuitive feel.

By changing an effective scaled diffusion coefficient from small to large, we can change the tail substantially, see Figure 2. We call it effective because the stochastic smearing on D and Length makes it scale-free and we can longer tell if the mean in D or Length is greater. We could have a huge mean for D and a small mean for Length, or vice versa, but we could not distinguish between the cases, unless we have measurements at more locations.

Figure 2 : Impulse response with increasing diffusion coefficient top to bottom.
The term x represents time, not position .

In practice, we won't have a heat impulse as a stimulus. A much more common situation involves a step input for heat. The unit step response is the integral of the scaled impulse response

The integral shows how the heat sink target transiently draws heat from the source. If the effective diffusion coefficient is very small, an outlet for heat dispersal does not exist and the temperature will continue to rise. If the diffusion coefficient is zero, then the temperature will increase linearly with time, t (again this is without a radiative response to provide an outlet).

Figure 3 : Unit step response of dispersed thermal diffusion. The smaller the effective
thermal diffusion coefficient, the longer the heat can stay near the source.

Eventually the response will attain a square root growth law, indicative of a Fick's law regime of what is often referred to as parabolic growth (somewhat of a misnomer). The larger the diffusion coefficient, the more that the response will diverge from the linear growth. All this means is that the heat is dispersively diffusing to the heat sink.

Application to AGW

This has implications for the "heat in the pipeline" scenario of increasing levels of greenhouse gases and the expected warming of the planet. Since the heat content of the oceans are about 1200 times that of the atmosphere, it is expected that a significant portion of the heat will enter the oceans, where the large volume of water will act as a heat sink. This heat becomes hard to detect because of the ocean's large heat capacity; and it will take time for the climate researchers to integrate the measurements before they can conclusively demonstrate that diffusion path.

In the meantime, the lower atmospheric temperature may not change as much as it could, because the GHG heat gets diverted to the oceans. The heat is therefore "in the pipeline", with the ocean acting as a buffer, capturing the heat that would immediately appear in the atmosphere in the absence of such a large heat sink. The practical evidence for this is a slowing of the atmospheric temperature rise, in accordance with the slower sqrt(t) rise than the linear t. However, this can only go on so long, and when the ocean's heat sink provides a smaller temperature difference than the atmosphere, the excess heat will cause a more immediate temperature rise nearer the source, instead of being spread around.

In terms of AGW, whenever the global temperature measurements start to show divergence from the model, it is likely due to the ocean's heat capacity. Like the atmospheric CO2, the excess heat is not "missing" but merely spread around.

EDIT:
The contents of this post are discussed on The Missing Heat isn't Missing at all.

I mentioned in comments that the analogy is very close to sizing a heat sink for your computer’s CPU. The heat sink works up to a point, then the fan takes over to dissipate that buffered heat via the fins. The problem is that the planet does not have a fan nor fins, but it does have an ocean as a sink. The excess heat then has nowhere left to go. Eventually the heat flow reaches a steady state, and the pipelining or buffering fails to dissipate the excess heat.

What's fittingly apropos is the unification of the two "missing" cases of climate science.

1. The "missing" CO2. Skeptics often complain about the missing CO2 in atmospheric measurements from that anticipated based on fossil fuel emissions. About 40% was missing by most accounts. This lead to confusion between the ideas of residence times versus adjustment times of atmospheric CO2. As it turns out, a simple model of CO2 diffusing to sequestering sites accurately represented the long adjustment times and the diffusion tails account for the missing 40%. I derived this phenomenon using diffusion of trace molecules, while most climate scientists apply a range of time constants that approximate diffusion.

2. The "missing" heat. Concerns also arise about missing heat based on measurements of the average global temperature. When a TCR/ECS* ratio of 0.56 is asserted, 44% of the heat is missing. This leads to confusion about where the heat is in the pipeline. As it turns out, a simple model of thermal energy diffusing to deeper ocean sites may account for the missing 44%. In this post, I derived this using a master heat equation and uncertainty in the parameters. Isaac Held uses a different approach based on time constants.

So that is the basic idea behind modeling the missing quantities of CO2 and of heat -- just apply a mechanism of dispersed diffusion. For CO2, this is the Fokker-Planck equation and for temperature, the heat equation. By applying diffusion principles, the solution arguably comes out much more cleanly and it will lead to better intuition as to the actual physics behind the observed behaviors.

I was alerted to this paper by Hansen et al (1985) which uses a box diffusion model. Hansen’s Figure 2 looks just like my Figure 3 above. This bends over just like Hansen’s does due to the diffusive square root of time dependence. When superimposed, it is not quite as strong a bend as shown in Figure 4 below.

Figure 4: Comparison against Hansen's model of diffusion

This missing heat is now clarified in my mind. In the paper Hansen calls it “unrealized warming”, which is heat entering into the ocean without raising the climate temperature substantially.

EDIT:
The following figure is a guide to the eye which explains the role of the ocean in short- and long-term thermal diffusion, i.e. transient climate response. The data from BEST illustrates the atmospheric-land temperatures, which are part of the fast response to the GHG forcing function. While the GISTEMP temperature data reflects more of the ocean's slow response.

Figure 5: Transient Climate Response explanation

Figure 6: Hansen's original projection of transient climate sensitivity plotted against the GISTEMP data,
which factors in ocean surface temperatures.

*
TCR = Transient Climate Response
ECS = Equilibrium Climate Sensitivity

The belief in Chaos

2012-01-22T09:41:00.000-08:00

The Chief says :

“Nothing just happens randomly in the Earth climate system. Randomness – or stochasticity – is merely a statistical approach to things you haven’t understood yet. ”

One of the unsung achievements in physics, in comparison to the imagination-capturing aspects of relativity and quantum mechanics, is statistical mechanics. This will scale at many levels -- originally intended to bridge the gap between the microscopic theory and macroscopic measurements, such as with the Planck response, scientists have provided statistical explanations to large coarse-grained behaviors as well (wind, ocean wave mechanics, etc). It's not that we don't understand the chaotic underpinnings, more like that we don't always need to, due the near-universal utility of the Boltzmann partition function (see the discussion on the Thermodynamics Climate Etc thread).

Many scientists consider pawning off difficulties to "Chaos" as a common crutch. This is not my original thought, as it is discussed at depth in Science of Chaos or Chaos in Science" by Bricmont. The issue with chaos theories is that they still have to obey some fundamental ideas of energy balance and conservation laws. Since stochastic approaches deal with probabilities, one rarely experiences problems with the fundamental bookkeeping. The basic idea with probability, that it has to integrate to unity probability, making it a slick tool for basic reasoning. That is why I like to use it so much for my own basic understanding of climate science (and all sorts of other things), but unfortunately leads to heated disagreements to the chaos fans and non-linear purists, such as David Young and Chief Hydrologist. They are representative of the opposite side of the debate.

You notice this when Chief states the importance of chaos theory:

"You should try to understand and accept that – along with the reality that my view has considerable support in the scientific literature. You should accept also that I am the future and you are the past.

I think they sould teach the 3 great ideas in 20th centruy physics – relativity, quantum mechanics and chaos theory. They are such fun."

There are only 4 fundamental forces in the universe, gravity, electromagnetism, and the strong and weak nuclear forces. For energy balance of the earth, all that matters is the electromagnetic force, as that is the predominant way that the earth exchanges energy with the rest of the universe.

The 33 degree C warming temperature differential from the earth's gray-body default needs to be completely explained by a photonic mechanism.

The suggestion is that clouds could change the climate. Unfortunately this points it in the incorrect direction of explaining the 33C difference. Water vapor, when not condensed into droplets, acts as a strong GHG and likely does cause a significant fraction of the 33C rise. But when the water vapor starts condensing into droplets and thus forming clouds, the EM radiation begins to partially reflect the incoming radiation, and thus the sun providing even less heat to the earth. So obviously there is a push-pull effect to raising water vapor concentrations in the atmosphere.

Chief is daring us with his statement that "I am the future and you are the past". He evidently thinks that clouds are the feedback that will not be understood unless we drop down to chaos considerations. In other words, that any type of careful statistical considerations of the warming impact of increasing water vapor concentrations with the cooling impact of cloud albedo, will not be explainable unless a full dynamical model is attempted and done correctly.

The divide is between whether one believes as Chief does, that the vague "chaos theory", which is really short-hand for doing a complete dynamical calculation of everything, no exceptions, is the answer. Or is the answer one of energy balance and statistical considerations? I lean toward the latter, along with the great majority of climate scientists, as Andrew Lacis described a while ago here and in his comments. The full dynamics, as Lacis explained is useful for understanding natural variability, and for practical applications such as weather prediction. But it is not the bottom-line, as chaotic natural variability always has to obey the energy balance constraints. And the only practical way to do that is by considering a statistical view.

The bottom-line is that I chuckle at much of the discussion of chaos and non-linearity when it comes to try to understand various natural phenomenon. The classic case is the simplest model of growth described by the logistic differential equation. This is a non-linear equation with a solution described by the so-called logistic function. Huge amounts of work have gone into modeling growth using the logistic equation because of the appearance of an S-shaped curve in some empirical observations. (when it is a logistic difference equation, chaotic solutions result but we will ignore that for this discussion)

Alas, there are trivial ways of deriving the same logistic function without having to assume non-linearity or chaos; instead one only has to assume disorder in the growth parameters and in the growth region. The derivation takes a few lines of math (see the TOC).

Once one considers this picture, the logistic function arguably has a more pragmatic foundation based on stochastics than on non-linear determinism.

That is the essential problem of invoking chaos, in that it precludes (or at least masks) considerations of the much more mundane characteristics of the system. The mundane is that all natural behaviors are smeared out by differences in material properties/characteristics, variation in geometrical considerations, and in thermalization contributing to entropy.

The issue is that obsessives such as the Chief and others think that chaos is the hammer and that they can apply it to every problem that appears to look like a nail.

Certainly, I can easily understand how the disorder in a large system can occasionally trigger tipping points or lead to stochastic resonances, but these are not revealed by analysis of any governing chaotic equations. They simply result from the disorder allowing behaviors to penetrate a wider volume of the state space. When these tickle the right positive feedback modes of the system, then we can observe some of the larger fluctuations. The end result is that the decadal oscillations are of the order of a tenths of degrees in global average temperature.

Of course I am not wedded to this thesis, just that it is a pragmatic result of stochastic and uncertainty considerations that I and a number other people are interested in.

Wave Energy Spectrum

2012-01-17T12:14:00.000-08:00

Ocean waves are just as disordered as the wind. We may not notice this because the scale of waves is usually smaller. In practice, the wind energy distribution relates to an open water wave energy distribution via similar maximum entropy disorder considerations. The following derivation assumes a deep enough water such that troughs do not touch bottom

First, we make a maximum entropy estimation of the energy of a one-dimensional propagating wave driven by a prevailing wind direction. The mean energy of the wave is related to the wave height by the square of the height, H. This makes sense because a taller wave needs a broader base to support that height, leading to a scaled pseudo-triangular shape, as shown in Figure 1 below.

Figure 1: Total energy in a directed wave goes as the square of the height, and the macroscopic fluid
properties suggest that it scales to size. This leads to a dispersive form for the wave size distribution

Since the area of such a scaled triangle goes as H^2, the MaxEnt cumulative probability is:

P(H) = exp(-a*H^2)

where a is related to the mean energy of an ensemble of waves.

We can stop at this point, since this relationship is empirically observed from measurements of ocean wave heights over a sufficient time period. However, we can proceed further and try to derive the dispersion results of wave frequency, which is another very common measure.

So as a next stage, we consider -- based on the energy stored in a specific wave -- the time, t, it will take to drop a height, H, by the Newton's law relation:

t^2 ~ H

and since t goes as 1/f, then we can create a new PDF from the height cumulative as follows:

p(f)*df = dP(H)/dH * dH/df * df

where

H ~ 1/f^2

dH/df ~ -1/f^3

then

p(f) ~ 1/f^5 * exp(-c/f^4)

which is just the Pierson-Moskowitz wave spectra that oceanographers have observed for years (developed first in 1964, variations of this include the Bretschneider and ITTC wave spectra) .

This concise derivation works well despite the correct path of calculating an auto-correlation from p(f) and then deriving a power spectrum from the Fourier Transform of p(f). Yet this convenient shortcut remains useful in understanding the simple physics and probabilities involved.

As we have an interest in using this derived form for an actual potential application, we can seek out public-access stations to obtain and evaluate some real data. The following Figure 2 is data pulled from the first region I accessed -- a pair of measuring stations located off the coast of San Diego. The default data selector picked the first day of this year, 1/1/2012 and the station server provided an averaged wave spectra for the entire day. The red points correspond to best fits from the derived MaxEnt algorithm to the blue data set.

Figure 2: Wave energy spectra from two sites off of the San Diego coastal region.
The Maximum Entropy estimate is in red.

To explore the dataset, here is a link to the interactive page :
http://cdip.ucsd.edu/?nav=historic&sub=data&units=metric&tz=UTC&pub=public&map_stati=1,2,3&stn=167&stream=p1&xyrmo=201201&xitem=product25

Like the wind energy spectrum, the wave spectrum derives simply from maximum entropy conditions.

Refs
Introduction to physical oceanography, RH Stewart

The basic algorithm

2011-12-07T22:01:00.000-08:00

A statement was made by a climate change skeptic:

"The existence of fossil fuels is due to the greater plant productivity in a warm world with lots more CO2."

That is what is so bizarre and surreal about the skeptical AGW discussions. They tend to devolve into either elliptical contradictions or else reveal self-consistent and circular truths.

We have low-levels of CO2 now
Levels of CO2 seem to be rising.
This will lead to a warmer climate due to the GHG effect
Wait a second ! --- heat and extra CO2 may be good after all.
That will promote more plant growth
It also might explain why we have lots of fossil fuels.
Millions of years for the plants to get buried by sedimentation and tectonic activity
That's how much of the CO2 was sequestered, lowering atmospheric concentration
But now we are digging it all up
And burning it, releasing CO2 back into the air
Go to step #2

Now we have to believe in the theory of CO2 assisted climate change. The algorithmic (or is that AlGoreRythmic?) steps have already been laid out. The problem is that the skeptics can't argue their way out of this basic trickbox.

Multiscale Variance Analysis and Ornstein-Uhlenbeck of Temperature Series

2011-11-13T19:41:00.000-08:00

Can we get behavioral information from something as seemingly random and chaotic as the Vostok ice core temperature data? Yes, as it turns out, performing a stochastic analysis of any time series can yield some very fundamental behaviors of the underlying process. The premise is that we can model the historical natural variability that occurs over hundreds of thousands of years by stochastic temperature movements.

I took the approach of modeling a bounded random walk from two different perspectives. In Part 1, I take the classic approach of applying the fundamental random walk hopping formulation.

Part 1 : Ornstein-Uhlenbeck Random Walk

Figure 1 below shows the Vostok data collected at roughly 20 to 600 year intervals and above that a simulated time series of a random walk process that has a similar variance.

Fig. 1: Ornstein-Uhlenbeck Random Walk process (top green) emulates
Vostok temperature variations (below blue)

The basic model I assume is that some type of random walk (red noise) statistics are a result of the earth’s climate inhabiting a very shallow energy well in its quiescent state. See the Figure 2 below. The well is broad on the order of 10 degrees C, and relatively shallow. Within this potential well, the climate hovers in a metastable state whereby small positive feedbacks of CO2, water vapor, or albedo can push it in either direction. The small amounts of perturbation needed to move the global temperature represents this shallowness.

Fig. 2 : Schematic of random walk behavior of a metastable system occupying a shallow quasi-equilibrium

The model of correlated Markovian random walk, whereby the perturbation can go in either a positive or negative direction each year matches this shallow energy well model remarkably well. The random walk can't grow without bound as soft temperature constraints do exist, mainly from the negative feedback of the Stefan-Boltzmann law.

To model these soft constraints, what we do is put a proportional drag on the random walk so the walk produces a noticeable asymmetry on large excursions from the quiescent point, Tq (I don’t call it equilibrium because it is really a quasi-equilibrium). This generates a reversion to the mean on the Markov random walk hopping direction, which is most commonly associated with the general Ornstein-Uhlenbeck process.

The way that random walk works is that if something gets to a certain point proportionately to the square root of time, then you can often interpolate how far the excursion is in a fraction of this time, see Figure 3 below.

$$ \Delta T = \sqrt{Dt} $$


Fig. 3 : Multiscale Variance applied to Monte Carlo simulated random walk. The classic Square Root random walk is shown for reference.

The agreement with the Vostok data is illustrated by Figure 4 below. The multiscale variance essentially takes all possible pairs of points and plots the square root of the average temperature as a function of the time interval, L.

Fig. 4 : Multiscale variance of Vostok data

The Ornstein-Uhlenbeck kernel solution has a simple formula for the variance:

$$Temperature(L) = F(L) = \sqrt{\frac{h}{2 \sigma}(1-e^{-2 \sigma L})}$$

where h is the hop rate and σ is the drag factor. We assume the mean is zero for convenience.

The very simple code for generating a Monte Carlo simulation and the multiscale variance is shown in the awk code snippet below. The BEGIN portion sets up the random walk parameters and then performs a random flip to determine whether the temperature will move high or low. The movements labeled A and C in Figure 2 correspond to random numbers less than 1/2 (positive hops) and the movements labeled B and D correspond to random numbers greater than 1/2 (negative hops). This nearly canonical form captures reversion to the mean well if the mean is set to zero.

The hopping rate, h, corresponds to a diffusion coefficent of 0.0008 degrees^2/year and a drag parameter is set to 0.0012, where the drag causes a reversion to the mean for long times and the resulting hop rate is proportionally reduced according to how large the excursion is away from the mean.

BEGIN { ## Set up the Monte Carlo O-U model
   N = 330000
   hop = 0.02828
   drag = 0.0012
   T[0] = 0

   srand()
   for(i=1; i<N; i++) {
     if(rand()<0.5) {
       T[i] = T[i-1] + hop*(1-T[i-1]*drag)
     } else {
       T[i] = T[i-1] - hop*(1+T[i-1]*drag)
     }
   }
}

{} ## No stream input

END { ## Do the multiscale variance
   L = N/2
   while(L > 1) {
     Sum = 0
     for(i=1; i<=N/2; i++) {
        Val = T[i] - T[i+int(L)]
        Sum += Val * Val
     }
     print int(L) "\t" sqrt(Sum/(N/2))
     L = 0.9*L
   }
}

The END portion calculates the statistics for the time series T(i). The multiscale variance is essentially described by the following correlation measure. All pairs are calculated for a variance and the result is normalized as a square root.

$$ F(L) = [\frac{2}{N}\sum_{j = 1}^{N/2} ( Y_j - Y_{j+L})^2]^{\frac{1}{2}} \quad \forall \quad L \lt N/2 $$

The O-U curve is shown in Figure 4 by the red curve. The long time-scale fluctuations in the variance-saturated region are artifacts based on the relative shortness of the time series. If a longer time series is correlated then these fluctuations will damp out and match the O-U asymptote.

The results from 8 sequential Monte Carlo runs of several hundred thousand years are illustrated in Figure 5 below. Note that on occasion the long term fluctuations coincide with the Vostok data. I am not suggesting that some external forcing factors aren’t triggering some of the movements, but it is interesting how much of the correlation structure shows up at the long time ranges just from random walk fluctuations.

Fig. 5 : Multiscale Variance for Vostok data and
various sampled Ornstein-Uhlenbeck random walk processes

In summary, the Vostok data illustrates the characteristic power-law correlation at shorter times and then it flattens out at longer times starting at about 10,000 years. This behavior results from the climate constraints of the maximum and minimum temperature bounds, set by the negative feedback of S-B at the high end and various factors at the low end. So to model the historical variations over a place like Vostok, we start with a random walk and modify it with an Ornstein-Uhlenbeck drag factor so that as the temperature starts walking closer to the rails, the random walk starts to slow down. By the same token, we give the diffusion coefficient a bit of a positive-feedback push when it starts to go away from the rail.

The question is whether the same random walk extends to the very short time periods ranging from years to a few centuries. Vostok doesn’t give this resolution as is seen from Figure 4 but other data sets do. If we can somehow extend this set of data to a full temporal dynamic range, this would provide an fully stochastic model for local climatic variation.

Figure 6 below shows a Greenland core data set which provides better temporal resolution than the Vostok set.
http://www.ncdc.noaa.gov/paleo/metadata/noaa-icecore-2490.html
ftp://ftp.ncdc.noaa.gov/pub/data/paleo/icecore/greenland/summit/ngrip/ngrip2006d15n-ch4-t.txt

Fig. 6 : Full Greenland time series generated at one year intervals.

Yet, this data is qualified by the caveat that the shorter time intervals are "reconstructed". You can see this by the smoother profile at a shorter time scale shown in Figure 7.

Fig. 7 : At shorter scale note the spline averaging which smooths out the data.

This smoothening reconstruction is a significant hurdle in extracting the actual random walk characteristics at short time scales. The green curve in Figure 8 below illustrates the multiscale variance applied to the set of 50,000 data points. The upper blue curve is a Monte Carlo execution of the excursions assuming O-U diffusion at a resolution of 1 year, and the lower red curve is the multiscale variance applied after processing with a 250 point moving average filter.

Fig. 8 : Multiscale Variance of Ornstein-Uhlenbeck process
does not match Greenland data (blue line)
unless a moving average filter of 250 points is applied (red line).

Higher resolution temperature time series are available from various European cities. The strip chart Figure 9 below shows a time series from Vienna starting from the 18th century.

Fig. 9 : Vienna Temperature time series

The multiscale variance shown in Figure 10 suggests that the data is uncorrelated and consists of white noise.

Fig. 10 : Multiscale Variance on data shows flat variance due to uncorrelated white noise.
Rise at the end due to AGW at current duration.

Sensitivity

This description in no way counters the GHG theory but is used to demonstrate how sensitive the earth’s climate is to small perturbations. The fact that the temperature in places such as Vostok follows this red noise O-U model so well, indicates that it should also be hypersensitive to artificial perturbations such as that caused by adding large amounts of relatively inert CO2 to the atmosphere. The good agreement with a Ornstein-Uhlenbeck process is proof that the climate is sensitive to perturbations, and the CO2 we have introduced is the added forcing function stimulus that the earth is responding to. For the moment it is going in the direction of the stimulus, and the red noise is along for the ride, unfortunately obscuring the observation at the short time scales we are operating under.

Figure 11 below shows a couple of shorter duration simulations using the Vostok parameters. Note that the natural fluctuations can occasionally obscure any warming that may occur. I explored this in the previous post where I calculated cumulative exceedance probabilities here.

Fig. 11 : Simulated time series showing the scale of the red noise

Extrapolating the Vostok data, for an arbitrary 100 year time span, a 0.75 degree or more temperature increase will occur with only a 3% probability. In other words the one we are in right now is a low chance occurrence. It is not implausible, but just that the odds are low.

Alternatively, one could estimate the fraction of the increase due to potential variability with some confidence interval attached to it. This will generate an estimate of uncertainty for the forced fraction versus the natural variability fraction.

The model of Ornstein-Uhlenbeck random walk could include a drift term to represent a forcing function. For a long term data series we keep this low so that it only reaches a maximum value after a specified point. In Figure 12 below, a drift term of 0.0002 degrees/year causes an early variance increase after 5000 years. This may indeed have been a natural forcing function in the historical as it creates an inflection point that matches the data better than the non-forced Ornstein-Uhlenbeck process in Figure 4.

Fig. 12 : Adding a small biased positive forcing function to the hop rate creates a strong early variance
as the temperature shifts to warmer temperatures sooner

Fig. 13 : Another simulation run which matches the early increase in variance seen in the data

The code below illustrates the small bias parameter, force, added to the diffusion term, hop, which adds a small drift term that adds to the variance increase over a long time period.

BEGIN { ## Set up the Monte Carlo O-U model
   N = 330000
   hop = 0.02828
   drag = 0.0012
   force = 0.00015
   T[0] = 0.0

   srand()
   for(i=1; i<N; i++) {
     if(rand()<0.5) {
       T[i] = T[i-1] + (hop+force)*(1-T[i-1]*drag)
     } else {
       T[i] = T[i-1] - (hop-force)*(1+T[i-1]*drag)
     }
     if(i%100 == 0) print T[i]
   }
}

The final caveat is that these historical time series are all localized data proxies or local temperature measurements. They likely give larger excursions than the globally averaged temperature. Alas, that data is hard to generate from paleoclimate data or human records.

Part 2 : Semi-Markov Model

The first part made the assumption that a random walk can generate the long-term temperature time series data. The random-walk generates a power-law of 1/2 in the multiscale variance and the Ornstein-Uhlenbeck master equation provides a reversion-to-the-mean mechanism. One has to look at the model results cross-eyed to see halfway good agreement, so that the validity of this approach for matching the data is ambiguous at best.

Because the temperature excursions don't have orders of magnitude in dynamic range, the power-law may be obscured by other factors. For example, the Vostok data shows large relative amounts of noise over the recent short duration intervals, while the Greenland appears heavily filtered.

That raises the question of whether an alternate stochastic mechanism is operational. If we consider the Greenland data as given, then we need a model that generates smooth and slow natural temperature variations at a short time scale, and not the erratic movements of Brownian motion that the O-U supplies.

One way to potentially proceed is to devise a doubly stochastic model. Along the time dimension, the overall motions are longer term, with the hopping proceeding in one direction until the temperature reverses itself. This is very analogous to a semi-Markov model as described by D.R. Cox and often seen in phenomenon such as the random telegraph signal as described in The Oil ConunDRUM.

A semi-Markov process supplies a very disordered notion of periodicity to the physical behavior. It essentially occupies the space between a more-or-less periodic waveform, with the erratic motion of a pure random walk.

The math behind a semi-Markov process is very straightforward. On the temporal dimension, we apply a probability distribution for a run length, t.
$$p(t) = \alpha e^{-\alpha t}$$
To calculate a variance of the temperature excursion, T, we use the standard definition for the second moment.

$$\sigma_T\,^2 = \int_0^\infty T(t)^2 \cdot p(t) dt$$
Here, p(t) is related to the autocorrelation describing long-range order and a typical PDF is a damped exponential:
$$p(t) = \alpha e^{- \alpha t}$$
The simplification we apply assumes that the short-term temperature excursion is
$$T(t) = \sqrt{D \cdot t}$$
Then to calculate a multiscale variance, we take the square root to convert it to the units of temperature,
$$ var(T) = \sqrt{\sigma_T\,^2}$$
we supply the parameterizations to the integration, and keep a running tally up to the multiscale time of interest, L.
$$var(T(L)) = \sqrt{\int_0^L \alpha e^{-\alpha t} \, Dt \, dt} = \sqrt{D (1-e^{-\alpha L} (\alpha L+1))/\alpha}$$
That is all there is to it. We have a short concise formulation, which has exactly two fitting parameters, one for shifting the time scale
$$ \tau = \frac{1}{\alpha}$$
and the diffusion parameter, D, scaling the temperature excursion.

To test this model, we can do both an analytical fit and Monte Carlo simulation experiments.

We start with the Greenland data as that shows a wide dynamic range in temperature excursions, see Figure 14 below.

Fig. 14 : Greenland multiscale variance with semi-Markov model

The thought was that the Greenland temperature data was reconstructed, with the background white noise removed at the finest resolution. This model reveals the underlying movement of that temperature.
$$D = 0.25 \,\, \mathsf{degrees^2/year}$$
$$ \alpha = 0.005 \,\, \mathsf{year^{-1}}$$
A quick check for asymptotic variance is
$$ var_\infty = \sqrt{D/\alpha} = 7 ^\circ C $$

Fig. 15 : Simulation of Greenland data assuming Semi-Markov with first-order lag bounds.

The simulation in Figure 15 uses the code shown below. Instead of square root excursions, this uses first-order lag movements, which shows generally similar behavior as the analytical expression.

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

main() { // Set up the Monte Carlo model
   #define N 100000
   double char_time = 100.0;
   double hop = 0.02;
   double drag = 0.03;
   int flip = 1;
   double *T = new double[2*N];
   int n = 1;
   int R, L, i;
   double Sum, Val, rnd;

   T[0] = 0.0;
   srand ( time(NULL) );

   while(n < N) {
     flip = -flip;
     rnd = (double)rand()/(double)RAND_MAX;
     R = -(int)(char_time * log(rnd));
     for(i=n; i<=n+R; i++)
        if(flip<0) {
          T[i] = T[i-1] + hop*(1-T[i-1]*drag);
        } else {
          T[i] = T[i-1] - hop*(1+T[i-1]*drag);
        }
     n += R;
   }
   // Do the multiscale variance
   L = N/2;
   while(L > 1) {
     Sum = 0.0;
     for(i=1; i<N/2; i++) {
        Val = T[i] - T[i+L];
        Sum += Val * Val;
     }
     printf("%i\t%g\n", L, sqrt(Sum/(N/2)));
     L = 0.1*L;
   }
}

Greenland shows large temperature volatility when compared against the ice core proxy measurements of the Antarctic. The semi-Markov model when applied to Vostok works well in explaining the data (see Figure 16), with the caveat that a uniform level of white noise is added to the variance profile. The white noise data is uncorrelated so adds the same background level independent of the time scale, i.e. it shows a flat variance. Whether this is due to instrumental uncertainty or natural variation, we can't tell.

Fig. 16 : Semi-Markov model applied to Vostok data

Assuming a white noise level of 0.9 degrees, the parametric data is:
$$D = 0.25 \,\, \mathsf{degrees^2/century}$$
$$ \alpha = 0.018 \,\, \mathsf{century^{-1}}$$
The check for asymptotic variance is
$$ var_\infty = \sqrt{D/\alpha} = 3.7 ^\circ C $$

The remarkable result is that the Vostok temperature shows similar model agreement to Greenland but the semi-Markov changes are on the order of 100's of centuries, instead of 100's of years for Greenland.

Another Antarctic proxy location is EPICA, which generally agrees with Vostok but has a longer temperature record. Figure 17 and Figure 18 also demonstrates the effect of small sample simulations and how that smooths out with over-sampling.

Fig. 17 : EPICA data with model and simulation

Fig. 18 : EPICA data with model and oversampling simulation

Part 3 : Epilog

The semi-Markov model described in Part 2 shows better agreement than the model of random walk with Ornstein-Uhlenbeck process described in Part 1. This better agreement is explained by the likely longer-term order implicit in natural climate change. In other words, the herky-jerky movements of random walk are less well-suited as a process model than longer term stochastic reversing trends appear. The random walk character still exists but likely has a strong bias depending on whether the historical is warming or cooling.

If we can get more of these proxy records for other locations around the world, we should be able to produce better models for natural temperature variability, and thus estimate exceedance probabilities for temperature excursions within well-defined time scales. With that we can estimate whether current temperature trends are out of the ordinary in comparison to historical records.

I found more recent data for Greenland which indicates a likely shift in climate fluctuations starting approximately 10,000 years ago.

Fig. 19 : From http://www.gps.caltech.edu/~jess/ "The figure below shows the record of oxygen isotope variation, a proxy for air temperature, at the Greenland Summit over the past 110,000 years. The last 10,000 years, the Holocene, is marked by relative climatic stability when compared to the preceding glacial period where there are large and very fast transitions between cold and warm times."

The data in Figure 19 above and Figure 20 below should be compared to the Greenland data described earlier, which ranged from 30,000 to 80,000 years ago.

Fig. 20 : Multiscale variance of Greenland data over the last 10,000 years, not converted to temperature.
Note the flat white noise until several thousand years ago, when a noticeable uptick starts.

The recent data definitely shows a more stable warmer regime than the erratic downward fluctuations of the past interglacial years. So the distinction is one between recent times and those of older than approximately 5,000 years ago.

See the post and comments on Climate Etc. regarding climate sensitivity. The following chart (Figure 21) illustrates possible mechanisms for a climate well. Notice the potential for an asymmetry in the potential barriers on the cool versus warm side.

Fig. 21: First Order Climate Well

The following Figure 22 is processed from the 30,000 to 80,000 year range Greenland ice core data. The histogram of cooling and warming rates is tabulated across the extremes (a) on the cool side below average (b) and on the warm side above average (c). Note the fast regime for warming on the warm side.

Fig. 22 : Detection of Skewed Well.
On the high side of paleo temperatures, a fast temperature rate change is observed.
A range of these these rates is equally likely, indicating a fast transition between warm and extreme.

Compare that to Figure 23 over the last 10,000 years (the Holocene era) in Greenland. The split is that rates appear to accelerate as a region transitions from cold to warm and then vice versa. I think we can't make too much out of this effect as the recent Holocene-era data in Greenland is very noisy and often appears to jump between extremes in temperature in consecutive data points. As described earlier this is characteristic of white noise (and perhaps even aliasing) since no long-range order is observed. More recent high resolution research has been published [1] but the raw data is not yet available. The claim is that the climate in Greenland during the Holocene appears highly sensitive to variations in environmental factors, which I find not too surprising.

Fig. 23 : Greenland last 10,000 years.
Note the slope emphasis on the warm side versus the cold side of temperature regime.

[1] Kobashi, T., K. Kawamura, J. P. Severinghaus, J.‐M. Barnola, T. Nakaegawa, B. M. Vinther, S. J. Johnsen, and J. E. Box (2011), High variability of Greenland surface temperature over the past 4000 years estimated from trapped air in an ice core, Geophys. Res. Lett., 38, L21501, doi:10.1029/2011GL049444.

Vostok Ice Cores

2011-10-13T23:42:00.000-07:00

As the partial pressure of CO2 in sea water goes like
$$ c = c_0 * e^{-E/kT}$$
and the climate sensitivity is
$$ T = \alpha * ln(c/c_1) $$
where c is the mean concentration of CO2, then it seems that one could estimate a quasi-equilibrium for the planet’s temperature. Even though they look nasty, these two equations actually solve to a quadratic equation, and one real non-negative value of T will drop out if the coefficients are reasonable.
$$T = \alpha * ln(c/c_1) - \alpha * E / kT $$
For CO2 in fresh water, the activation energy is about 0.23 electron volts. From "Global sea–air CO2 flux basedon climatological surface ocean pCO2, and seasonal biological and temperature effects" by Taro Takahashi,

The pCO2 in surface ocean waters doubles for every 16C temperature increase
(d ln pCO2/ dT=0.0423 C ).

This gives 0.354 eV.

We don’t know what c1 or c0 are and we can use estimates of climate sensitivity for α is (between 1.5/ln(2) and 4.5/ln(2)).

When solving the quadratic the two exponential coefficients can be combined as
$$ ln(w)=ln(c_0/c_1)$$
then the quasi-equilibrium temperature is approximated by this expansion of the quadratic equation.
$$ T = \alpha ln(w) – \frac{E}{k*ln(w)} $$
What the term “w” means is the ratio of CO2 in bulk to that which can effect sensitivity as a GHG.

As a graphical solution of the quadratic consider the following figure. The positive feedback of warming is given by the shallow-sloped violet curve, while the climate sensitivity is given by the strongly exponentially increasing curve. Where the two curves intersect, not enough outgassed CO2 is being produced such that the asymptotically saturated GHG can further act on. The positive feedback essentially has "hit a rail" due to the diminishing return of GHG heat retention.

We can use the Vostok ice core data to map out the rail-to-rail variations. The red curves are rails for the temperature response of CO2 outgassing, given +/- 5% of a nominal coefficient, using the activation energy of 0.354 eV. The green curves are rails for climate sensitivity curves for small variations in α.

This may be an interesting way to look at the problem in the absence of CO2 forcing. The points outside of the slanted parallelogram box are possibly hysteresis terms causes by latencies of in either CO2 sequestering or heat retention. On the upper rail, the concentration drops below the expected value, while as drops to the lower rail, the concentration remains high for awhile.

The cross-correlation of Vostok CO2 with Temperature:

Temperature : ftp://ftp.ncdc.noaa.gov/pub/data/paleo/icecore/antarctica/vostok/deutnat.txt
CO2 core : ftp://ftp.ncdc.noaa.gov/pub/data/paleo/icecore/antarctica/vostok/co2nat.txt

The CO2 data is in approximately 1500 year intervals while the Temperature data is decimated more finely. The ordering of the data is backwards from the current date so the small lead that CO2 shows in the above graph is actually a small lag when the direction of time is considered.

The top chart shows the direction of the CO2:Temperature movements. Lots of noise but a lagged chart will show hints of lissajous figures, which are somewhat noticeable as CCW rotations for a lag. On temperature increase, more of the CO2 is low than high, as you can see it occupying the bottom half of the top curve.

The middle chart shows where both CO2 and T are heading in the same direction. The lower half is more sparsely populated because temperature shoots up more sharply than it cools down.

The bottom chart shows where the CO2 and Temperature are out-of-phase. Again T leads CO2 based on the number you see on the high edge versus the low edge. The lissajous CCW rotations are more obvious as well.

Bottom line is that Temperature will likely lead CO2 because I can’t think of any Paleo events that will spontaneously create 10 to 100 PPM of CO2 quickly, yet Temperature forcings likely occur. Once set in motion, the huge adjustment time of CO2 and the positive feedback outgassing from the oceans will allow it to hit the climate sensitivity rail on the top.

So what is the big deal? We don’t have a historical forcing of CO2 to compare with, yet we have one today that is 100 PPM.

That people is a significant event, and whether it is important or mot we can rely on the models to help.

This is what the changes in temperature look like over different intervals.

The changes follow the MaxEnt estimator of a double sided damped exponential. A 0.2 degree C change per decade(2 degree C per century) is very rare as you can see from the cumulative.

That curve that runs through the cumulative density function (CDF) data is a maximum entropy estimate. The following constraint generated the double-sided exponential or Laplace probability density function (PDF) shown below the cumulative:
$$\int_{I}{|x| p(x)\ dx}=w$$
which when variationally optimized gives
$$p(x)={\beta\over 2}e^{-\beta|x|},\ x\ \in I=(-\infty,\infty)$$
where I fit it to:
$$\beta = 1/0.27$$
which gives a half-width of about +/- 0.27 degrees C.

The Berkeley Earth temperature study shows this kind of dispersion in the spatially separated stations.

Another way to look at the Vostok data is as a random up and down walk of temperature changes. These will occasionally reach high and low excursions corresponding to the interglacial extremes. The following is a Monte Carlo simulation of steps corresponding to 0.0004 deg^2/year.

The trend goes as:
$$ \Delta T \sim \sqrt{Dt}$$
Under maximum entropy this retains its shape:

This can be mapped out with the actual data via a Detrended Fluctuation Analysis.
$$ F(L) = [\frac{1}{L}\sum_{j = 1}^L ( Y_j - aj - b)^2]^{\frac{1}{2}} $$
No trend in this data so the a coefficient was set to 0. This essentially takes all the pairs of points, similar to an autocorrelation function but it shows the Fickian spread in the random walk excursions as opposed to a probability of maintaining the same value.

The intervals are a century apart. Clearly it shows a random walk behavior as the square root fit goes though the data until it hits the long-range correlations.

Sea temperature correlation

2011-10-10T21:05:00.000-07:00

Temperature Induced CO2 Release Adds to the Problem

2011-10-05T20:53:00.000-07:00

As a variable amount of CO2 gets released by decadal global temperature changes, it makes sense that any excess amount would have to follow the same behavior as excess CO2 due to fossil fuel emissions.

From a previous post (Sensitivity of Global Temperature), I was able to detect the differential CO2 sensitivity to global temperature variations. The correlation of temperature anomaly against d[CO2] is very strong with zero lag and a ratio of about 1 PPM change in CO2 per degree temperature change detected per month.

Now, this does not seem like much of a problem, as naively a 1 degree change over a long time span should only add one PPM during the interval. However, two special considerations are involved here. First, the measure being detected is a differential rate of CO2 production and we all know that sustained rates can accumulate into a significant quantities of a substance over time. Secondly, the atmospheric CO2 has a significant adjustment time and the excess isn't immediately reincorporated into sequestering sites. To check this, consider that a slow linear rate of 0.01 degree change per year when accumulated over 100 years will lead to a 50 PPM accumulation, if the excess CO2 is not removed from the system. This is a simple integration where f(T(t)) is the integration function :
$$ [CO2] = f_{co_2}(T(t)) = \int^{100}_0 0.01 t\, dt = \frac{1}{2} 0.01 * 100^2 = 50 $$
The sanity check on this is if you consider that a temperature anomaly of 1 degree change held over 100 years would release 100 PPM into the atmosphere. This is simply a result of Henry's Law applied to the ocean. The ocean has a large heat capacity and so will continue outgassing CO2 at a constant partial-pressure rate as long as the temperature has not reached the new thermal equilibrium. (The CO2 doesn't want to stay in an opened Coke can, and it really doesn't want to stay there when it gets warmed up)

So, if we try the impulse response we derived earlier (Derivation of MaxEnt Diffusion) to this problem, with a characteristic time that matches the IPCC model for Bern CC/TAR, standard:

As another sanity check, the convolution of this with a slow 1 degree change over the course of 100 years will lead to at least a 23 PPM CO2 increase.

Again, this occurs because we are far from any kind of equilibrium, with the ocean releasing the CO2 and the atmosphere retaining what has been released. The slow diffusion into the deep sequestering stores is just too gradual while the biotic carbon cycle is doing just that, cycling the carbon back and forth.

So now we are ready to redo the model of CO2 response to fossil-fuel emissions (Fat-Tail Impulse Response of CO2) with the extra positive feedback term due to temperature changes. This is not too hard as we just need to get temperature data that goes back far enough (the HADCRUT3 series goes back to 1850). So when we do the full combined convolution, we add in the integrated CO2 rate term f(T), which adds in the correction as the earth warms.

$$ [CO2] = FF(t) \otimes R(t) + f_{co_2}(T(t)) \otimes R(t) $$

When we compute the full convolution, the result looks like the following curve (baseline 290 PPM):

The extra CO2 addition is almost 20 PPM just as what we had predicted from the sanity check. The other interesting data feature is that it nearly recreates the cusp around the year 1940. The previous response curve did not pick that up because it is entirely caused by the positive-feedback warming during that time period. The effect is not strong but discernible.

We will continue to watch how this plays out. What is worth looking into is the catastrophic increase of CO2 that will occur as long as the temperature stays elevated and the oceans haven't equilibrated yet.

TOTAL Discovery Data

2011-09-30T23:07:00.000-07:00

From Figure 1 of Laherrere's recent discovery data, he suggests a crude oil URR of 2200 GB.

The oil company Total S.A. also has an accounting of yearly discoveries. I overlaid their data with Laherrere's data in the figure below:

Total must do a bit of a different backdating because their data is consistently above Laherrere's for the majority of the years. As of about 2005, their cumulative is at 2310 GB while Laherrere is at 1930 GB.

This leads to the following asymptotic graph for cumulative oil according to the Dispersive Discovery model. In this case I assigned a URR of 2800 GB to the model, with a model value of 2450 GB as of 2005. In other words, the Total discovery data may hit an asymptote of 2800 GB, which may be a bit generous:

This is really for comparative purposes as I next plotted what Laherrere's discovery data looks like against the same model.

You can see that Laherrere's data likely won't hit that asymptote.

Discovery data for crude oil is hard to come by. Perhaps Laherrere is using 2P probabilities and Total is applying possible reserve growth so that it is >2P? Or perhaps Total is using barrels of oil equivalent (BOE) to inflate the numbers (which Shell oil does)? In the greater scheme of things, does this really matter?

The following chart is the Shock Model applied to the Total discovery data, whereby I tried to follow the historical crude oil (not All Liquids) production rates by varying the extraction rate until the year 2000, then I kept the extraction rate constant at about 3.1% of reserves. This is lower than the currently accepted 4% to 5% extraction rate from reserves.

If Total does use BOE maybe this should actually fit an All Liquids curve, in which case the extraction rates would need to get increased to match the higher levels of All Liquids production.

Bottom line is that the peak plateau might extend for a couple of years and we will have a fatter decline tail if we believe the Total numbers. If it is an All Liquids discovery model, then it is a wash. As if everyone didn't know this by now, peak oil is not about the cumulative, it is about the extractive flow rates, and this is a good example of that.

In general, the URR is incrementally getting pushed up with time. Laherrere had used 2000 GB for a crude oil URR for some time (see this curve from 2005) and now likely because of the deep water oil it is at 2200 GB.

As for going through the trouble of evaluating the Gulf Of Mexico data, that is just noise on the overall curve IMO. It's getting to the point that we have enough historical data that global predictions for crude oil are really starting to stabilize. And the Dispersive Discovery model will anticipate any future discoveries. The longer we wait, the closer all the estimates will start to converge, and any new data will have little effect on the projected asymptotic cumulative.

The following is a curve that takes the Total S.A. discovery data and extrapolates the future discoveries with a dispersive discovery model. The final discovery URR is 2700 billion barrels, which is quite a bit higher than the one Laherrere plots. This is higher because I am making the assumption that Total S.A. is including backdated NGPL and other liquids along with the crude oil. Which means I had to fit against the production data that also used these liquids.

To model the perturbations in production levels, which is necessary to accumulate the reserves properly, I used the Oil Shock Model. In the inset, you can see the changes in extraction rate that occurred over the years. The extraction rate is essentially the same as the Production/Reserve ratio. Notice that the extraction rate was steady until the 1960's at which it ramped up. It started to level off and drop down during the 1970's oil crisis and didn't really start to rise again until the 1990's. I am extrapolating the extraction rate from today to match the peak extraction rate of the 1960's by the year 2050.

This is largely a descriptive working model, which essentially reflects the data that Total S.A. is providing and then reflecting that in terms of what we are seeing in the production numbers. The current plateau could be extended if we try to extract even faster (as in Rockman's words "PO is about oil flow rate") or we can start including other types of fuels to the mix. This latter will happen if the EIA and IEA add biofuels and other sources to the yearly production.

The bottom-line is that it is hard to come up with any scenarios, based on the data that Total and IHS supplies , that can extend this plateau the way that Total suggests it will, peak to 2020. That is what is maddening about this whole business and you wonder why drivel such as what Yergin continues to pump out gets published.

Derivation of MaxEnt Diffusion applied to CO2 adjustment times

2011-09-28T23:07:00.000-07:00

hunter | September 28, 2011 at 11:07 pm |
WHT,
Oil does not give a hoot about your math.
You are an historical illiterate and an imbecile hiding behind your self-declared genius.

The diffusion kernel resulting from solving the Fokker-Planck equation (sans drift term) is:

$$ \large C(t,x|D) = \frac{1}{\sqrt{4 \pi D t}} e^{-x^2/{4 D t}} $$

We place an impulse of concentrate at x=0 and want to watch the evolution of the concentration with time. We have an idea of a mean value for the diffusion coefficient, D, but don't know how much it varies. The remedy for that is to apply a maximum entropy estimate for the variance assuming a mean value D0.

$$ \large p_d(D) = \frac{1}{D_0} e^{-D/D_0} $$

So then we can apply this to the kernel function

$$ \large C(t,x) = \int_0^{\infty} C(t,x|D) p_d(D) dx $$

This actually comes out quite clean (see derivation in The Oil Conundrum book)

$$ \large C(t,x) = \frac{1}{2 \sqrt{D_0 t}} e^{-x/{\sqrt{D_0 t}}} $$

This gives us a result that shows a singularity at t=0 for x=0. In practice, the value of x is not precise, so that we can also place an uncertainty around the value of x.

$$ \large p_x(x) = \frac{1}{x_0} e^{-x/x_0} $$

Once again we can apply this to the concentration, marginalizing x out of the picture:

$$ \large C(t) = \int_0^{\infty} C(t|x) p_x(x) dx $$

This integral is very straightforward

$$ \large C(t) = \frac{1}{2} \frac{1}{x_0 + \sqrt{D_0 t}} $$

which is precisely the form, apart from normalization, used in fitting to the IPCC curves!

See fat-tail-impulse-response-of-co2. This is the impulse response function I used:

$$ \large [CO_2] = f(t) \otimes R(t) = f(t) \otimes \frac{1}{1+0.15\sqrt{t}} $$

So comparing the values, 0.15, is the square root reciprocal of a median diffusion time given by:

$$ \sqrt{t_0} = \frac{x_0}{\sqrt{D_0}} $$

This comes out to 44 years. However, because of the ratio, we lost the ability to figure out either the average diffusion coefficient or the average displacement though. The fat-tail comes about because 44 years is only a median, and a mean diffusion time does not exist.

The chart below is not definitionally correct, as it describes the impulse response as a residence time and not an adjustment time, but you can see how the diffusional model fits on top of the climate science simulations. It falls off rather quickly but then assumes the slow hyperbolic square root growth which matches well over many orders of magnitude. The one key difference is that the diffusional model will asymptotically go to zero, but it will take a long time, as long a geological time as it took for the original fossil fuel to form in all likelihood.

From the textbook "Introduction to Organic Geochemistry", I applied this fit using a value of 0.19 instead of 0.15 for the scaling coefficient. Every chart of the CO2 impulse response I have encountered shows excellent agreement with this very characteristic fat-tailed decline.

We can create a brilliant analogy of CO2 with coins in the transaction system.

Residence time is a coin getting freshly minted and finding out how long it circulates before it gets handed-off as change. When a different coin is exchanged, the person has lost track of its lineage but that does not matter as long as an equivalent coin takes its place. It doesn’t matter if the coin had identifying markers or not, as the excess coins are still circulating.

Adjustment time is a coin cycling through the system long enough that it gets decommissioned or gets lost or destroyed. It then permanently gets removed from the system.

The coin in the financial transaction system is equivalent to the CO2 molecule in the carbon cycle.

Now look at the relative rates of residence time of a coin versus that of the adjustment time. You will find that a typical residence time is measured in weeks, but the adjustment time is years. This is exactly the same rationale of a CO2 residence time in years and a CO2 adjustment time in centuries.

These kinds of studies have been done on circulation of money, look up the Where's George project. I was able to use this data to generate a model for travel patterns of people here. Perhaps not remarkably, this data shows the same fat-tail statistics that the CO2 does. This is all based on standard random-walk and dispersion arguments that come up time and time again.

Fat-Tail Impulse Response of CO2

2011-09-24T00:15:00.000-07:00

This is a combination of the last two posts. The first post called Explaining "The Missing Carbon" reviewed the fat-tail response of atmospheric CO2 adjustment time to a carbon forcing function, and detailed how the Fokker-Planck (FP) diffusion models sequestering of CO2 into deep stores. The first step is to approximate the impulse response, R(t):

$$ \large [CO_2] = f(t) \otimes R(t) = f(t) \otimes \frac{1}{1+0.15\sqrt{t}} $$

Next we convolve this impulse response with the fossil fuel estimates, f(t), from the last 260 years (data from the Carbon Dioxide Information Analysis Center at Oak Ridge National Labs)

Next we convert the carbon emissions in metric tons to a CO2 concentration in ppm, and plot it in comparison to the historically measured CO2 concentrations at Mauna Loa with the NASA GISS global temperature anomaly alongside. We choose a baseline of 290 ppm because that fits better than 280 ppm (and this is in agreement with a previous estimate made).

Now we zero in on the CO2 sensitivity so that we can compare it to the last post "The sensitivity of global temperature to CO2". The BLUE curve below shows the yearly deviations as differential CO2, short-handed as d[CO2], from the convolved impulse response, calculated since 1960. The RED curve below is yearly d[CO2] data from NOAA. The GREEN curve is data directly computed from the fossil fuel emissions, that is no impulse response convolution.

This data removes all the seasonal adjustments and though there is likely a weak correlation with measured d[CO2], most of the CO2 variation is associated with multi-year temperature fluctuations.This link http://www.skepticalscience.com/print.php?r=145 points to work by (Bacastow and Keeling 1981) and Section 7.3.2.4 of the IPCC AR4 Working Group 1 report. The latter has all the details:

[IPCC] — The atmospheric CO2 growth rate exhibits large interannual variations (see Figure 3.3, the TAR and http://lgmacweb.env.uea.ac.uk/lequere/co2/carbon_budget). The variability of fossil fuel emissions and the estimated variability in net ocean uptake are too small to account for this signal, which must be caused by year-to-year fluctuations in land-atmosphere fluxes. Over the past two decades, higher than decadal-mean CO2 growth rates occurred in 1983, 1987, 1994 to 1995, 1997 to 1998 and 2002 to 2003. During such episodes, the net uptake of anthropogenic CO2 (sum of land and ocean sinks) is temporarily weakened. Conversely, small growth rates occurred in 1981, 1992 to 1993 and 1996 to 1997, associated with enhanced uptake.

Those years do indeed match, just odd that no one ever thought to actually plot the numbers and show the cross-correlation. This seems like such obvious scientific book-keeping that I am dumb-founded by the lack of a published plot.

[IPCC] — Since the TAR, many studies have confirmed that the variability of CO2 fluxes is mostly due to land fluxes, and that tropical lands contribute strongly to this signal (Figure 7.9). A predominantly terrestrial origin of the growth rate variability can be inferred from (1) atmospheric inversions assimilating time series of CO2 concentrations from different stations (Bousquet et al., 2000; Rödenbeck et al., 2003b; Baker et al., 2006), (2) consistent relationships between δ13C and CO2 (Rayner et al., 1999), (3) ocean model simulations (e.g., Le Quéré et al., 2003; McKinley et al., 2004a) and (4) terrestrial carbon cycle and coupled model simulations (e.g., C. Jones et al., 2001; McGuire et al., 2001; Peylin et al., 2005; Zeng et al., 2005). Currently, there is no evidence for basin-scale interannual variability of the air-sea CO2 flux exceeding ±0.4 GtC yr–1, but there are large ocean regions, such as the Southern Ocean, where interannual variability has not been well observed.

Take a look at figure 7.9 in particular and one can see how the trends change quite a bit for CO2 measurements over ocean versus land. Curiously Mauna Loa is in the ocean yet it shows the strong correlation of the land stations.

In any case, the long-term atmospheric CO2 concentration is clearly explainable by a fossil-fuel forcing function. The short-term deviations in CO2 are explained by a combination of fluctuating carbon emissions along with a strong dependence on natural short-term temperature fluctuations (shown below).

with a strong cross-correlation between Temperature and the Proportional-Derivative d[CO2] model:

The sensitivity of global temperature to CO2

2011-09-12T23:24:00.000-07:00

Doing exploratory analysis with the CO2 perturbation data at the Wood For Trees data repository and the results are very interesting. The following graph summarizes the results:

Rapid changes in CO2 levels track with global temperatures, with a variance reduction of 30% if d[CO2] derivatives are included in the model. This increase in temperature is not caused by the temperature of the introduced CO2, but is likely due to a reduced capacity for the biosphere to take-up the excess CO2.

This kind of modeling is very easy if you have any experience with engineering controls development. The model is of the type called Proportional-Derivative, and it essentially models a first-order equation
$$\Delta T = k[CO_2] + B \frac{d[CO_2]}{dt}$$

The key initial filter you have to apply is to average the Mauna Loa CO2 data over an entire year. This gets rid of the seasonal changes and the results just pop out.
The numbers I used in the fit are B=1.1 and k=0.0062. The units are months. The B coefficient is large because the CO2 impulse response has got that steep downslope which then tails off:
$$ d[CO_2] \sim \frac{1}{\sqrt{t}}$$

Is this a demonstration of causality, +CO2 => +Temperature?

If another causal chain supports the change in CO2, then likely. We have fairly good records of fossil fuel (FF) emissions over the years. The cross-correlation of the yearly changes, d[CO2] and d[FF], show a zero-lag peak with a significant correlation (below right). The odds of this happening if the two time-series were randomized is about 1 out of 50.

Left chart from Detection of global economic fluctuations in the atmospheric co2 record. This is not as good a cross-correlation as the d[CO2] and dTemperature data -- look at year 1998 in particular, but the zero-lag correlation is clearly visible in the chart..

This is the likely causality chain:
$$d[FF] \longrightarrow d[CO_2] \longrightarrow dTemperature$$

If it was the other way, an increase in temperature would have to lead to both CO2 and carbon emission increases independently. CO2 could happen because of outgassing feedbacks (CO2 in oceans is actually increasing despite the outgassing), but I find it hard to believe that the world economy would increase FF emissions as a result of a warmer climate.

What happens if there is a temperature forcing CO2 ?

With the PD model in place the of d[CO2] against Temperature cross-correlation looks like the following:

The Fourier Transform set looks like the following:

This shows the two curves have the same slope in spectrum and just a scale shift. The upper curve is the ratio between the two curves and is essentially level.

The cross-correlation has zero lag and a strong correlation of 0.9. The model again is
$$ \Delta T = k[CO_2] + B \frac{d[CO_2]}{dt}$$

The first term is a Proportional term and the second is the Derivative term. I chose the coefficients to minimize the variance between the measured Temperature data and the model for [CO2]. In engineering this is a common formulation for a family of feedback control algorithms called PID control (the I stands for integral). The question is what is controlling what.

When I was working with vacuum deposition systems we used PID controllers to control the heat of our furnaces. The difference is that in that situation, the roles are reversed, with the process variable being a temperature reading off a thermocouple and the forcing function is power supplied to a heating coil as a PID combination of T. So it is intuitive for me to immediately think that the [CO2] is the error signal, yet that gives a very strong derivative factor which essentially amplifies the effect. The only way to get a damping factor is by assuming that Temperature is the error signal and then we use a Proportional and an Integral term to model the [CO2] response. Which would then give a similar form and likely an equally good fit.

It is really a question of causality, and the controls community have a couple of terms for this. There is the aspect of Controllability and that of Observability (due to Kalman).

Controllability: In order to be able to do whatever we (Nature) want
with the given dynamic system under control input, the system must be controllable.
Observability: In order to see what is going on inside the system under observation, the system must be observable.

So it gets to the issue of two points of view:
1. The people that think that CO2 is driving the temperature changes have to assume that nature is executing a Proportional/Derivative Controller on observing the [CO2] concentration over time.
2. The people that think that temperature is driving the CO2 changes have to assume that nature is executing a Proportional/Integral Controller on observing the temperature change over time, and the CO2 is simply a side effect.

What people miss is that it can be potentially a combination of the two effects.
Nothing says that we can’t model something more sophisticated like this:
$$ c \Delta T + M \int {\Delta T}dt = k[CO_2] + B \frac{d[CO_2]}{dt}$$

The Laplace transfer function Temperature/CO2 for this is:
$$ \frac {s(k + B s)}{c s + M} $$

Because of the s in the numerator, the derivative is still dominating but the other terms can modulate the effect.

This blog post did the analysis a while ago, the image below is fascinating because the overlay between the dCO2 and Temperature anomaly matches to the point that the noise even looks similar . This doesn't go back to 1960.

If CO2 does follow Temperature, due to the Causius-Clapeyron and the Arrhenius rate law, a positive feedback will occur -- as the released CO2 will provide more GHG which will then potentially increase temperature further. It is a matter of quantifying the effect. It may be subtle or it may be strong.

From the best cross-correlation fit, the perturbation is either around (1) 3.5 ppm change per degree change in a year or (2) 0.3 degree change per ppm change in a year.

(1) makes sense as a Temperature forcing effect as the magnitude doesn’t seem too outrageous and would work as a perturbation playing a minor effect on the 100 ppm change in CO2 that we have observed in the last 100 years.
(2) seems very strong in the other direction as a CO2 forcing effect. You can understand this if we simply made a 100 ppm change in CO2, then we would see a 30 degree change in temperature, which is pretty ridiculous, unless this is a real quick transient effect as the CO2 quickly disperses to generate less of a GHG effect.

Perhaps this explains why the dCO2 versus Temperature data has been largely ignored. Even though the evidence is pretty compelling, it really doesn’t further the argument on either side. On the one side interpretation #1 is pretty small and on the other side interpretation #2 is too large, so #1 may be operational.

One thing I do think this helps with is providing a good proxy for differential temperature measurements. There is a baseline increase of Temperature (and of CO2), and accurate dCO2 measurements can predict at least some of the changes we will see beyond this baseline.

Also, and this is far out, but if #2 is indeed operational, it may give credence to the theory that that we may be seeing the modulation of global temperatures the last 10 years because of a plateauing in oil production. We will no longer see huge excursions in fossil fuel use as it gets too valuable to squander, and so the big transient temperature changes from the baseline no longer occur. That is just a working hypothesis.

I still think that understanding the dCO2 against Temperature will aid in making sense of what is going on. As a piece in the jigsaw puzzle it seems very important although it manifests itself only as a second order effect on the overall trend in temperature. In summary, as a feedback term for Temperature driving CO2 this is pretty small but if we flip it and say it is 3.3 degrees change for every ppm change of CO2 in a month, it looks very significant. I think that order of magnitude effect more than anything else is what is troubling.

One more plot of the alignment. For this one, the periodic portion of the d[CO2] was removed by incorporating a sine wave with an extra harmonic and averaging that with a kernel function for the period. This is the Fourier analysis with t=time starting from the beginning of the year.
$$ 2.78 \cos(2 \pi t - \theta_1) + 0.8 \cos(4 \pi t - \theta_2) $$
$$ \theta_1 = 2 $$
$$ \theta_2 = -0.56 $$
phase shift in radians. The yearly kernel function is calculated from this awk function:

BEGIN {
I=0
}
{
n[I++] = $1
}
END {
Groups = int(I / 12)
## Kernel function
for(i=0; i<=12; i++) {
    x = 0
    for(j=0; j<Groups; j++) {
      x += n[j*12+i]
    }
    G[i] = x/Groups
}
Scale = (G[12]-G[0])
for(i=0; i<=12; i++) {
    Y[i] = (G[i]-G[0]) -i*Scale/12
}

for(j=0; j<Groups; j++) {
    for(i=0; i<12; i++) {
      Diff = n[j*12+i] - Y[i]
      print Diff
    }
}
}

This was then filtered with a 12 month moving average. It looks about the same as the original one from Wood For Trees, with the naive filter applied at the source, and it has the same shape for the cross-correlation. Here it is in any case; I think the fine structure is a bit more apparent(the data near the end points is noisy because I applied the moving average correctly).

How can the derivative of CO2 track the temperature so closely? My working theory assumes that new CO2 is the forcing function. An impulse of CO2 enters the atmosphere and it creates an impulse response function over time. Let's say the impulse response is a damped exponential and the atmospheric temperature responds quickly to this profile.

The CO2 measured at the Mauna Loa station takes some time to disperse over from the original source points. This implies that a smearing function would describe that dispersion, and we can model that as a convolution. The simplest convolution is an exponential with an exponential, as we just need to get the shape right. But what the convolution does is eliminate the strong early impulse response, and thus create a lagged response. As you can see from the Alpha plot below, the way we get the strong impulse back is to take the derivative. What this does is bring the CO2 signal above the sample-and-hold characteristic caused by the fat-tail. The lag disappears and the temperature anomaly now tracks the d[CO2] impulses.

If we believe that CO2 is a forcing function for Temperature, then this behavior must happen as well; the only question is whether the effect is strong enough to be observable.

If you realize that the noisy data below is what we started with, and we had to extract a non-seasonal signal from the green curve, one realizes that detecting that subtle a shift in magnitude is certainly possible.

Explaining the "Missing Carbon"

2011-09-11T08:42:00.000-07:00

Certain laws of physics can't be denied and the model of the carbon cycle is really set in stone. The fundamental law is one of mass balance and every physicist worth his salt is familiar with the master equation, also known as the Fokker-Planck formulation in its continuous form.
What we are really interested in is the detailed mass balance of carbon between the
atmosphere and the earth's surface. The surface can be either land or water, it doesn't matter for argument's sake.

We know that the carbon-cycle between the atmosphere and the biota is relatively fast and the majority of the exchange has a turnover of just a few years. Yet, what we are really interested in is the deep exchange of the carbon with slow-releasing stores. This process is described by diffusion and that is where we can use the Fokker-Planck to represent the flow of CO2.

This part can't be debated because this is the way that the flow of all particles works; they call it the master equation because it invokes the laws of probability and in particular the basic random walk that just about every physical phenomenon displays.

The origin of the master model is best described by considering a flow graph and
drawing edges between compartments of the system. This is often referred to as a
compartment or box model. The flows go both ways and are random, and thus model the random walk between compartments.
The following chart is a Markov model consisting of 50 stages of deeper sequestering with each slab having a constant but small hop rate. The single interface between the atmosphere and the earth has a faster hopping rate corresponding to faster carbon cycling.

That shows how you would solve the system numerically. The basic analytical
solution to the Fokker-Planck assuming a planar source and one-dimensional diffusion is the following:
$$ \frac1{\sqrt{2 \pi t}} exp(-x^2/{2t}) $$

Consider that x=0 near the surface, or at the atmosphere/earth interface.
Because of that, this expression can be approximated by
$$ n(t)=\frac{q}{\sqrt{t}} $$
where n(t) is the concentration evolution over time and q is a scaling factor for
that concentration.
First thing one notices about this expression is that n(t) has a fat tail and
after a rapid initial fall-off only slowly decreases over time. The physical meaning is that, due to diffusion, the concentration randomly walks between the interface and deeper locations in the earth. The square root of time dependence is a classic trait of all random walks and you can't escape seeing this if you have ever watched nature in action. That is just the way particles move around.

For CO2 concentration this in fact describes the evolution of the adjustment
time, and it accurately reflects the infamous IPCC curve for the atmospheric CO2 impulse response. It is called an impulse response because that is the response that one would expect based on an initial impulse of CO2 concentration.

But that is just the first part of the story. As an impulse response, n(t) describes what is called a single point source of initial concentration and its slow evolution. In practice, fossil-fuel emissions generate a continuous stream of CO2 impulses. These have to be incorporated somehow. The way this is done is by the mathematical technique called convolution.

So consider that the incoming stream of new CO2 from fossil fuel emissions is
called F(t). This becomes the forcing function.
Then the system evolution is described by the equation
$$ c(t)=n(t)*F(t) $$
where the operator * is not a multiplication but signifies convolution.

Again, there is no debate over the fundamental correctness of what has been said
so far. This is exactly the way a system will respond.

If we are now to put this into practice and see how well it describes the actual
evolution of CO2, we can understand ever nagging issue that has haunted skeptical observers. It really all becomes very clear.

For the forcing function F(t) we use a growing power law.
$$ F(t) = k t^N $$
where N is the power and k is a scaling constant.

This roughly represents the atmospheric emissions through the industrial era if
we use a power law of N=4. See the following curve:

So all we really want to solve is the convolution of n(t) with F(t). By using
Laplace transforms on the convolution expression, the answer comes out
surprisingly clean and concise. Ignoring the scaling factor :
$$ c(t) \sim t^{N+1/2} $$

With that solved, we can now answer the issue of where the "missing" CO2 went
to. This is an elementary problem of integrating the forcing function, F(t), over time and then comparing the concentration, c(t), to this value. Then this ratio of c(t) to the integral of F(t) is the amount of CO2 that remains in the atmosphere.
Working out the details, this ratio is:
$$ q \sqrt{\frac{\pi}{t}}\frac{(N+1)!}{(N+0.5)!} $$
Plugging in numbers for this expression, q=1, and N=4, then the ratio is about 0.
28 after 200 years of growth. This means that 0.72 of the CO2 is going back into
the deep-stores of the carbon cycle, and 0.28 is remaining in the atmosphere.
If we choose a value of q=2, then 0.56 remains in the atmosphere and 0.44 goes
into the deep store. This ratio is essentially related to the effective
diffusion coefficient of the carbon going into the deep store.

Come up with a good number for the diffusion coefficient and we have an explanation of the evolution of the "missing carbon".

BTW, This is useful for modifying the numbers.
alpha equation

Alexander Harvey comment

R(t,0) = α/Sqrt(t) where α can be derived from the diffusivity and thermal capacity.

The general case R(t,λ) is not easily stated but can be calculated by the deconvolution of R(t,0) with the unit impulse fucntion to which an impulse term 1/λ at t=0 is added and the whole is deconvoluted with the unit impulse function. This is I think the same as can be achieved using the Laplace transforms in the continuous case.

This comment is the closest to my own thinking on the topic. The 1/Sqrt(t) profile is likely the explanation for the "fat-tail" in the residence cum adjustment time that the IPCC has been publishing. Consider this (EVALUATING THE SOCIAL COSTS OF GREENHOUSE GAS)

To represent this carbon cycle, Maier-Reimer and Hasselmann (1987) have suggested that the carbon stock be represented as a series of five boxes, each with a constant, but different atmospheric lifetime. That is, each box is modelled as in equation (7), and atmospheric concentration is a weighted sum of all five boxes. Maier-Reimer and Hasselmann suggested lifetimes of 1.9 years, 17.3 years, 73.6 years, 362.9 years and infinity, to which they attach weights of 0.1, 0.25, 0.32, 0.2 and 0.13, respectively

I punched these numbers in and this is what it looks like in comparison to a disordered diffusion profile:

The boxes they talk about simply describes a model of maximum entropy disorder. If the maximum entropy disorder is governed by force (i.e. drift), the profile will go toward 1/(1+t/a) and if it is governed by diffusion (i.e. random walk) it will go as 1/(1+sqrt(t/a)). Selecting a weighted sum of exponential lifetimes is exactly the same as choosing a maximum entropy distribution of rates, and I believe that whatever their underlying simulation is that it will asymptotically trend toward this result. The lifetime they choose for infinity is needed because otherwise the fat-tail will eventually disappear.

This is the MaxEnt derivation
integral exp(-1/(x*t))*exp(-x)/sqrt(x*t)*dx from x=0 to x=infinity

An example of misapplication of Hubbert Linearization

2011-07-08T05:52:00.000-07:00

This past weekend, one of the principle editors at TOD, Heading Out, wrote a key-post on California oil depletion and an analysis technique known as Hubbert Linearization (HL). I have spent like 5 years on TOD arguing against the scientific credibility and overall applicability of this technique, suggesting really that it is nothing more than a heuristic. In my view of things, a heuristic is not based on science and amounts to nothing more than hand-waving. To understand this long battle that I have slogged through, you really have to realize that most oil depletion analysts are empiricists, which means that they like to look at trends in the data without needing or even caring to understand the fundamentals in a meaningful way. It is kind of like a climate scientist looking at only the trends in temperatures, without doing any extensive modeling to back up what they observe -- I dislike this curve fitting without modeling because it weakens potentially strong arguments one can make against the skeptics and the "Drill,Baby,Drill" crowd. It is really similar to the weather versus climate argument, the knife can cut both ways if you just watch the trends.

So I am taking Heading Out to task for a poor scientific argument and making some real fundamental mistakes while trying to apply the HL technique. See the short TOD post called Tech Talk - California Oil and Hubbert Linearization for some background. The HL portion of the curve is the straight line in the curve below; this line intercepts at a point at which you can estimate the ultimately recoverable cumulative oil. (These were originally plotted by Laherrere).

First note that the "Hubbert Linearization" curves as plotted here is done as Q (cumulative) on the x-axis and P (yearly) on the y-axis. This is not the conventional Hubbert Linearization as applied to a Hubbert Logistic curve. For that curve, the analysts usually plot P/Q (not just P) on the y-axis, which effectively linearizes the equation quoted in the box (i.e. "Stuart's explanation").

(from the TOD post http://www.theoildrum.com/node/8100)
So how are these plots derived?

This is a topic that has been covered almost since the time that The Oil Drum was founded (back in March 2005), since it was in May of that year that Jean (Laherrere) published a paper describing how to “Forecast Production from Discovery.” (His pioneering work, internationally recognized in this field, did not start with this. For example, he co-authored with Colin Campbell the Scientific American article on “The End of Cheap Oil” in 1998.)

Within the pages of TOD Stuart Staniford had briefly explained it while comparing different methods of estimating future production back in September 2005, with a follow-on post looking at specific examples where it might be applied.

The technique derives from a process known as Hubbert Linearization (after King Hubbert, who is largely remembered for predicting the date of peak oil production in the US before it happened). Examining the data from oil production over time, Dr Hubbert postulated that it followed a logistic curve, which as Stuart pointed out, is an accepted model of how exponential growth occurs in a system that is of a certain finite size. It has been used since its original discovery in 1838.

The mathematics of the equation are fairly straightforward, and for consistency I am going to quote Stuart’s explanation:

It's pretty obvious that this change of P vs P/Q makes a huge difference to the interpretation of the plotted curve.

So let me provide the actual mechanics involved.

As a simple premise, consider that one of the most common models of a single field depletion is what one would call a proportional drawdown. This essentially says that the amount you can draw per unit time (the rate) is approximately proportional to the amount left.

Mathematically this relation is expressed as dQ/dt = k(Q0-Q)
which when solved with initial conditions gives the damped exponential:

Q = Q0*(1-exp(-k*t))

where Q0 is the maximum cumulative. This is a basic example of the law of diminishing
returns and occurs in many other natural phenomenon.

Annual production for this curve is simply the slope of the Q curve:

P(t) = dQ/dt

and rewrite it like this:

P = Q0*k*(1-Q/Q0)

Lo and behold, this is clearly a linearized equation as well. Yet, this is not based on a
symmetric Hubbert curve, as the damped exponential is not even close to having the classic peaked shape.

Now, if the parameter k was to increase in direct proportion to Q, we would come up with the Logistic peak, but this would have to be so coincidental and contrived as to be completely absurd. And if this Logistic increase was introduced to these California data sets , then it would make a complete mess of the Laherrere curves as plotted -- they would no longer show linearity with the choice of axes!

Now you can also see how the earlier parts of the plot can be effected. Notice that
the drawdown rate, k, or the ultimate Q0, could potentially change over time. The Kern was one of those cases that benefited from technical advancements over time and ultimately held a large secondary reserve of very heavy oil that was exploited many years after it was originally discovered.

Therefore, everything that is said in this particular TOD post is explained by a much simpler interpretation than the brain-dead Hubbert Logistic. This has nothing to do with HL, and everything to do with first-order drawdown dynamics from a single field and some logical deductions as to what is happening.

All in all, this is another completely flawed analysis published in The Oil Drum. Heading Out included all this useless info on Hubbert Linearization without realizing it had little to do with the problem at hand.

I am now of the opinion that the charter of The Oil Drum doesn't involve wanting to do the math correctly, even if it approaches the trivial. Most of the TOD posters have this long-standing infatuation with heuristics and I have noticed that it always comes back to bite them when they try to "prove" something.

I would advise them to adapt to some current mathematical practices, like say calculus, logic, and probabilities. The old guard of Hubbert and Deffeyes had some insight but clearly didn't have all the answers. Some of the other analysts like Staniford have essentially lead us down blind alleys. Robert Rapier does indeed have it correct in suggesting we drop HL as "The Way" of proving anything. Ultimately this key-post was a case of completely misapplying the technique and then claiming that it serves as an example of where HL actually works.

It does not, however, always apply, and Robert Rapier has explained, in two posts (here, and here) why he has concerns about using the technique. But in terms of giving a ballpark for production (and recognizing that there are always new discoveries and inventions that can be, as they say, game changers) the technique has considerable support. And the consistency with which the Californian fields are following the predictions provide evidence for such an opinion.

Here is the disagreement:
The problem in the TOD piece is that Heading Out thinks he is showing Hubbert Linearization on Laherere's plots whereas he is actually only showing the equivalent of a proportional drawdown from a reservoir.

Here is the important premise:
Consider that someone has the good fortune to earn their income from a trust fund. But the terms of the trust fund are that you can only draw from the trust a percentage per year, k, from what is remaining in the fund. Say that the trust fund starts with a quantity Q0, then the amount you would be expected to draw from it yearly is P, which is the time rate of change of the amount left.

Here is the correct solution:
Mathematically, this gets expressed as

P = dQ/dt = k/100*(Q0-Q)

I don't have to go any further than this because the equation above is all that Heading Out is plotting. The total cumulative amount is the x-axis and the proportional amount is the y-axis. The trust fund analogizes to the oil reservoir and the production draw-down is proportional to the amount of oil remaining. The fit turns out to match proportional draw-down very well. However, this has nothing to do with the Hubbert peak or with Hubbert Linearization of the peak.

Here is the mistake:
Heading Out made the mistake in thinking that he was plotting the completely different

dQ/dt = k*Q*(Q0-Q)

which is the "classic" Hubbert Peak formulation. But that is not what he plotted as the extra multiplicative factor Q is not shown on either of Laherrere's axis. If he plotted it this way, the line through the data would no longer be straight, and Hubbert Linearization would FAIL :(

Here are the implications:
One of the hardest things in technical discussions is to demonstrate how someone did some calculation wrong. You can only try to speculate on why they chose their argument leading to the bad result. And then you run the risk of over-analyzing the situation and confusing someone who can't tell the difference between what is wrong or what is right. School teachers learn how to this endlessly with their students and turn out to be among the most patient people in the world. In the end, I lose patience and just say that the application of Hubbert Linearization that Heading Out performed is WRONG.

What this means is that Heading Out needs to retract that article and start from scratch. You can't put lipstick on a pig, or polish a turd. I wouldn't want to go beyond this and talk at all about how Hubbert Linearization is useful, because this study of California oil production is a case where it is clearly NOT useful. So sorry to say it doesn't work the way that Heading Out is thinking, and what I have written both in this post and in the past substantiates Rapier's views.

For me, these TOD posts are easy to respond to because I only have to refer to what I have written before, where I have tried to tie all the loose ends together:

aka "The Oil Conunundrum"

I stopped posting at my old blog because after I put the book together, I realized that it serves as kind of an archival record of what went into the book. I thought working scientific problems out in the open with some fresh thoughts along with real-time contributions from commenters could lead to a better analysis. We'll see how it turns out. In the meantime, if people have ideas of new places to post on these topics, I would be interested, otherwise I will continue to post here to combat what I think is intellectual laziness amongst the peak oil crowd.