Earth Cubed

Distributed Climate Science and Computing

Laplace Transform Via Limits

Quite a while back I had a brief interest in wheather the Laplace transform made numeric sense. I was interested in this for its own sake:

http://www.physicsforums.com/showthread.php?t=334191

and because if we could numericaly understand the Laplace transform it might tell us more about a dynamic system then a fourier transform.

https://earthcubed.wordpress.com/2009/08/30/using-the-fft-to-calculate-the-laplace-transform/

Memory of my failures in obtaining the desired results led me to believe that the Laplace transform was not Riemann integratable. However when chalanged on this in another forum:

http://forums.philosophyforums.com/threads/notitle-54250-12.html

I rethought a suggested solution and decided that I was perhaps partly incorrect. We can derive the integral from sums & limits bit it converges in a very strange way. The sum converges exactly where the geometric series diverges. All that is left prior to cancellation is a ratio of infinitesimally small quantities.

The summation suggested by andrewk for the laplace transform was:

(1) – F(s)= \int_0^{oo} f(t) dt = lim_{m->oo} lim_{n->oo} \sum_{k=1}^{n} f(km/n) \frac{m}{n}

It looks simple enough that I may have or should have considered it but it is not entirely obvious for how the double sum should work. Moreover, wikipedia discusses some problems with trying to use Riemann integration with improper integrals.

http://en.wikipedia.org/wiki/Riemann_integral#Generalizations

Needless to I am able to use a Riemann like result and get the right answer. Whether the approach is mathematically sound is an entirely different question. In the above equation if we are trying to integrate an exponentail function of the form

(2) – f(x)=exp(- \alpha t)

(3) – L( f(x)) = lim_{m->oo} lim_{n->oo} \sum_{k=1}^{n} e^{-skm/n} exp(- \alpha km/n) (m/n) 

combining terms

(4) – L[f(x)] =lim_{m->oo}lim_{n->oo}\sum_{k-1}^n exp(-(s+alpha)km/n) (m/n) 

Expressing as a geometric series:

(5) – L[f(x)]=lim_{m->oo}lim_{n->oo}\sum_{k-1}^n \left( e^{-(s+alpha)m/n} \right)^k (m/n) 

Using the geometric series:

(6) – L[f(x)]=lim_{m->oo}lim_{n->oo} \frac{1- \left( e^{-(s+alpha)m/n \right)^n}}{1-e^{-(s+alpha)m/n}} \frac{m}{n} 

Now here is where the Voodoo comes in we need to let m/n approach zero and m approach infinity. Rather them attempting to do this formally lets just cut to the chase.

when m/n approaches zero we can we can approximate:$latex e^{-(s+alpha)m/n} $ as 1 - (s+alpha) (m/n) 

and as m approach infinity \left( e^{-(s+alpha)m/n \right}^n   approaches zero.

Making these substitution into equation (6) we get:

(7) – L[f(x)]=lim_{m->oo}lim_{n->oo} \frac{1}{1-(1+(s+alpha)m/n)} \frac{m}{n}

which further simplifies to

(8) – L{f(x)}=lim_{m->oo}lim_{n->oo} \frac{1}{(s+alpha)m/n)} x m/n

cancelling m/n from the numerator and denominator (

(9) – L{f(x)}=lim_{m->oo}lim_{n->oo} \frac{1}{(s+alpha))}

and now the limits are irrelivant:

(10) – L{f(x)}=lim_{m->oo}lim_{n->oo} \frac{1}{(s+alpha))}

July 14, 2012 Posted by | Uncategorized | Leave a comment

log(CO2) and Scary Graphs

After reading Aurther’s blog post “New Congressional Budget Office Report on Climate Change” I got curious as to what the temperature response to CO2 would look like if it is actually logarithmic.

I descovered (what I should have realized from the start) that the curve is concave up but will converge to a linear curve for large T

It is commonly believed that the response of the earth to greenhouse gasses is logarithmic. I heard people suggest on other forums that this view was obtained empirically though climate models and perhaps more specifically radiative transfer models.

For analytic justification one would derive an expression for how the spectral width of an absorption peak grows with CO2 concentration (as I’ve done here). However, I am not convinced that this is sufficient justification as there will always be some radiative transfer at a given wave length even if the majority of it is absorbed over a very short distance. This is because the temperature gradient produced by the lapse rate produces an outward radiative flux that would exceed radiative feedback. I discuss this in more detail in my post:
Tropospheric Feedback

The logarithmic response is important because it is a type of saturation, in other words the more CO2 that is added to the atmosphere the less effective the next unit of CO2 will be in contributing to the warming. What I learned from reading Arther’s blog is that the current CO2 levels have not yet overwhelmed the natural levels of CO2.  This can be seen in the following graph:

More specificity from about 1000-1800 the CO2 concentration in the atmosphere stayed around 280 ppm. The following graph is more useful for measuring the current growth in CO2 concentration:

This graph is surprisingly very linear. If the growth in CO2 is truely exponential then it should be possible to estimate in from the slope on this graph which is given as 1.4203 PPM per year. For an exponential function:

y=y_o+Aexp(\lambda (t-t_o))
The derivative is:
y'=A\lambda exp( \lambda (t-t_o))
And the second derivative is:
y''=A\lambda^2 exp( \lambda (t-t_o))
The second derivative was taken because two equations are needed to find both A \lambda can be found.

The site also where I obtained the above figures gives a quadratic fit which can be used to estimate the first and second derivatives:

44690.5-46.1486x+0.0119942x^2

Therefore at year 2007 the first derivative is given by:

-46.1486+2*0.0119942*2007=1.9961

and the second derivative is 0.0119942

Giving:

A exp( \lambda (t-t_o)) \lambda =1.9961
A exp( \lambda (t-t_o)) \lambda^2 =0.0119942

Dividing the second equation by the first:

\lambda=0.0119942/1.9961=0.006008817

From this the doubling time can be obtained as follows:

ln(2)/\lambda=ln(2)/0.0060=115.5245 years.

y_o is taken to be the base level of CO2. That is:

y_0=280

A third equation can now be found as:

384-280=Aexp(0.006(1997-t_0))

giving:

A=\frac{384-280}{exp(0.006(1997-t_0))}

If t_o is taken to be 1800 this gives:

A=31.8931

which suggests the CO2 growth rate has decreased over the last 200 years.

The CO2 is estimated to follow this function:

CO_2(t)=280+31.8931 exp(0.006 (t-1800))

The question now is how does this growth rate in CO2 effect the temperature. There are several estimates for the sensitivity of the climate to changes in CO2. Lucia’s one box model “lumpy” suggest a sensitivity of:

1.7 degrees Celsius per CO2 doubling.
The IPCC estimates the lower bound for sesitivity to be:
1.5 degrees Celsius per CO2 doubling (see CO2 Climate Sensitivity)
Isaac M. Held suggests a climate sensitivity of about 2.8C/CO2 doubling. See:

http://www.gfdl.gov/isaac-held-homepage

Selected recent papers on climate sensitivity:

Here is what wikipedia has to say:

In Intergovernmental Panel on Climate Change (IPCC) reports, equilibrium climate sensitivity refers to the equilibrium change in global mean near-surface air temperature that would result from a sustained doubling of the atmospheric (equivalent) CO2 concentration. This value is estimated, by the IPCC Fourth Assessment Report (AR4) as likely to be in the range 2 to 4.5°C with a best estimate of about 3°C, and is very unlikely to be less than 1.5°C. Values substantially higher than 4.5°C cannot be excluded, but agreement of models with observations is not as good for those values. This is a slight change from the IPCC Third Assessment Report (TAR), which said it was “likely to be in the range of 1.5 to 4.5°C” [1]. AR3 defined climate sensitivity alternatively in systematic units, equilibrium climate sensitivity refers to the equilibrium change in surface air temperature following a unit change in radiative forcing and is expressed in units of °C/(W/m2) or equivalently K/(W/m2). In practice, the evaluation of the equilibrium climate sensitivity from models requires very long simulations with coupled global climate models, or it may be deduced from observations. Therefore the 2007 AR4 renamed the alternative climate sensitivity to climate sensitivity parameter adding a new definition of effective climate sensitivity which is “a measure of the strengths of the climate feedbacks at a particular time and may vary with forcing history and climate state”.

The logarithmic law of CO2 forcing is given as:

\Delta T=k*log_2(CO_2(t_2)/CO_2(t_1))=\frac{k}{log_n(2)}log_n((CO_2(t_2)/CO_2(t_1))

where k is the CO2 sensitivity for doubling CO2

I plotted this function for several values of the doubling sensitivity k

The labels on the right hand side of the plot are the climate sensitivities for each curve. This is actually a considerably smaller response then one would expect given the doubling time is around 100 years. However, while this is a nearly sufficient time for the exponential part of the curve to double,  the CO2 only increases by a factor of 1.2 since at 1997 the CO2 concentration is 384 ppm and in 2100 the CO2 concentration is projected by this fit to be  473ppm. This reduces the expected response by a factor of:

log(473/384)/log(2)=0.3

Notice that if the sensitivity ranges given by the IPCC are used then with this fit to CO2 emission growth is around 0.5-1.2 degrees which is hardly the doomsday scenario shown in the following graph which was posted on Aurther’s blog.

As a final note, the MATLAB code I used to produce the above graph is:

clear all
CO2= @(t)280+31.8931*exp(0.006*(t-1800))
CO2_0=CO2(1997)
t=linspace(1997,2100)
K=[1.5 4 7 11 15 21 25 29 33]
CO2s=CO2(t);
DT=@(k)(k/log(2))*log(CO2s/CO2_0)
for i=1:length(K)
  DTs=DT(K(i))
  plot(t,DTs)
  AXIS([1997 2100 0 10])
  gtext(num2str(K(i)))
  hold on;
end
xlabel('Year')
ylabel('Temperature Change in Degrees Celcius')
hold off,

September 19, 2009 Posted by | Uncategorized | 1 Comment

Vector Operations in Hoskins and Simmons Coordinates

In my post Hoskins and Simmons (1974) Coordinate System, I derived the following scaling quantities which will be used to derive the vector operations of Grad, Div and Curl in Hoskins Coordinate system.

h_{\mu}= \left| \frac{r}{\sqrt{1-\mu^2}} \right|
h_{\lambda}=r \mu
h_{\sigma} =-\sigma \frac{mg}{RT}

The coordinates in Hoskins coordinate system are dimensionless . (see nondimentionalization of Navier Stokes).

The gradient is defined as (see lectures on coordinate transforms):

\nabla \phi = \frac{\vec e_{\mu}}{h_{\mu}} \frac{\partial \phi}{\partial \mu} + \frac{\vec e_{\lambda}}{h_{\lambda}}\frac{\partial \phi}{\partial \lambda}+\frac{\vec e_{\sigma}}{h_{\sigma}}\frac{\partial \phi}{\partial \sigma}=\left(\vec e_{\mu}\right)\frac{1-\mu^2}{r} \frac{\partial \phi}{\partial \mu} + \left(\vec e_{\lambda}\right)\frac{1}{r \mu}\frac{\partial \phi}{\partial \lambda}-\left(\vec e_{\sigma}\right)\frac{RT}{\sigma mg}\frac{\partial \phi}{\partial \sigma}

The divergence is defined as:

\nabla \cdot = \frac{1}{h_{\mu}} \frac{\partial }{\partial \mu} + \frac{1}{h_{\lambda}}\frac{\partial}{\partial \lambda}+\frac{1}{h_{\sigma}}\frac{\partial}{\partial \sigma}=\frac{1-\mu^2}{r} \frac{\partial}{\partial \mu} + \frac{1}{r \mu}\frac{\partial}{\partial \lambda}-\frac{RT}{\sigma mg}\frac{\partial}{\partial \sigma}

The curl is defined by:

\nabla \times \vec A = \frac{1}{h_{\sigma}h_{\mu}h_{\lambda}} \left| \begin{array}{ccc} \vec e_{\sigma} & \vec e_{\mu} & \vec e_{\lambda} \\ \frac{\partial}{\partial \sigma} & \frac{\partial}{\partial \mu} & \frac{\partial}{\partial \lambda} \\ A_{\sigma}h_{\sigma} & A_{\mu}h_{\mu} &  A_{\lambda}h_{\lambda} \end{array} \right|

(note, the direction of the longitudinal coordinate \lambda is defined to obey the right hand rule)

This gives for the components

\left[ \nabla \times \vec A \right]_{\sigma}=\frac{RT}{r \mu \sigma mg}\left(\sigma \frac{mg}{RT}\frac{\partial A_{\lambda}}{\partial \mu}+ r \mu\frac{\partial A_{\mu}}{\partial \lambda} \right)

\left[ \nabla \times \vec A \right]_{\mu}= -\frac{\sqrt{1-\mu^2}RT}{r\sigma mg}\left( \frac{r}{\sqrt{1-\mu^2}}\frac{\partial A_{\sigma}}{\partial \lambda}+\sigma \frac{mg}{RT} \frac{\partial A_{\lambda}}{\partial \sigma} \right)

\left[ \nabla \times \vec A \right]_{\lambda}= \frac{\sqrt{1-\mu^2}RT}{r^2 \mu }\left( r \mu\frac{\partial A_{\mu}}{\partial \sigma}-\frac{r}{\sqrt{1-\mu^2}} \frac{\partial A_{\sigma}}{\partial \mu} \right)

Which Simplifies to:

\left[ \nabla \times \vec A \right]_{\sigma}= \frac{1}{r\mu}\frac{\partial A_{\lambda}}{\partial \mu}+ \frac{RT}{\sigma m g}\frac{\partial A_{\mu}}{\partial \lambda}

\left[ \nabla \times \vec A \right]_{\mu}=  -\frac{RT}{\sigma m g}\frac{\partial A_{\sigma}}{\partial \lambda}- \frac{\sqrt{1-\mu^2}}{r} \frac{\partial A_{\lambda}}{\partial \sigma}

\left[ \nabla \times \vec A \right]_{\lambda}=  \frac{\sqrt{1-\mu^2}RT}{r}\frac{\partial A_{\mu}}{\partial \sigma}-\frac{RT}{r\mu} \frac{\partial A_{\sigma}}{\partial \mu}

August 29, 2009 Posted by | Uncategorized | 1 Comment

The Coordinate System

Given that the surface of the eath is curved one might think that it would be difficult to construct a coordinate system that behaves like a rectangular coordinate system on a small scale but curves around the earth on a large scale.

To approach this problem we ask. How is a sphere like a cube? To a topologist they are nearly identical as it doesn’t take much effort to deform one into the other.

Imagine we had a cubic balloon. We blow it up a little and it is a cube. We draw a separate grid on each face of the cube shaped balloon. Then we blow it up more until it becomes a sphere. If our grid was fine enough then even with the balloon inflated to a spherical shape it still might look like a grid of squares is covering the balloon.

This should give confidence that we can divide the surface of a sphere up into sections that look like squares even though they cover a round object. So to do this we start thinking about what a square balloon might look like if it was inflated into a sphere.  The center of the top face will be mapped to the north pole. The center of the bottom face will be mapped to the south pole, the four faces on the side of the cube will divide the earth into four more sections. And the center point of each of these sections will be mapped to somewhere on the equator.

There are 360 degrees of longitude in the earth. Thus each section covers 90 degrees of longitude. If  we take the prime meridian as one edge then each face adjacent to the prime merridian will be centered at +/-45 degrees longitude.  Additionally the center points of the remain faces will be at +/-135 degrees of longitude.

These results are summarized in the following table. The remaining columns will be discussed in another post:

gridPoints

id

parent

level

left

right

top

bottom

up

down

lat

long

time

P_Theory

1

0

1

5

4

6

3

0

0

90

0

12

0.5

2

0

1

5

4

3

6

0

0

-90

0

12

0.5

3

0

1

5

4

1

2

0

0

0

45

12

0.5

4

0

1

3

6

1

2

0

0

0

135

12

0.5

5

0

1

6

3

1

2

0

0

0

-45

12

0.5

6

0

1

4

5

2

1

0

0

0

-135

12

0.5

March 1, 2009 Posted by | Computer Science and Modeling, Uncategorized | Leave a comment