Earth Cubed

Distributed Climate Science and Computing

Coriolis Forces in Hopkins and Simmons Vorticity Equation

In the thread vector operations in Hopkins and Simmons, I compute the components of the curl as:

\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]_{\sigma,\mu}+\left[ \nabla \times \vec A \right]_{\sigma,\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]_{\mu,\lambda}+\left[ \nabla \times \vec A \right]_{\mu,\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}=\left[ \nabla \times \vec A \right]_{\lambda,\sigma}+\left[ \nabla \times \vec A \right]_{\lambda,\mu}

In my post Coriolis forces in Hopkins and Simmons I compute the coriolis force as:

\boldsymbol{\Omega \times v} =  \begin{pmatrix} \mu U \\  \pm \sqrt{1-\mu^2} U \\ \mp \sqrt{1-\mu^2} V - \mu  W \end{pmatrix} \  = \begin{pmatrix} A_{\sigma} \\ A_{\mu} \\ A_{\lambda} \end{pmatrix} \

And the partial derivatives are given by:

\frac{\partial A_{\lambda}}{\partial \mu}=\frac{\partial}{\partial \mu}\left(\mp \sqrt{1-\mu^2} V - \mu  W\right)=\pm \frac{2\mu}{\sqrt{1-\mu^2}}V \mp \sqrt{1-\mu^2} \frac{\partial V}{\partial \mu} - W-\frac{\partial W}{\partial \mu}

\frac{\partial A_{\mu}}{\partial \lambda}=\frac{\partial }{\partial \lambda}\left(\pm \sqrt{1-\mu^2} U\right)=\pm \sqrt{1-\mu^2} \frac{\partial U}{\partial \lambda}

\frac{\partial A_{\sigma}}{\partial \lambda}=\frac{\partial }{\partial \lambda}\left( \mu U \right)=\mu \frac{\partial U}{\partial \lambda}

\frac{\partial A_{\lambda}}{\partial \sigma}=\frac{\partial}{\partial \sigma}\left(\mp \sqrt{1-\mu^2} V - \mu  W\right)=\mp \sqrt{1-\mu^2} \frac{\partial V}{\partial \sigma} - \mu  \frac{\partial W}{\partial \sigma}

\frac{\partial A_{\mu}}{\partial \sigma}=\frac{\partial }{\partial \sigma}\left(\pm \sqrt{1-\mu^2} U\right)=\pm \sqrt{1-\mu^2} \frac{\partial U }{\partial \sigma}

\frac{\partial A_{\sigma}}{\partial \mu}=\frac{\partial }{\partial \mu}\left( \mu U \right)=\mu \frac{\partial U}{\partial \mu}

I’ll derive the rest of this later but this doesn’t seem to be the form of the prognostic equation used by Hopkins and Simmons.

Divergence Free Flow

In the post Vector Operations in Hopkins and Simmons I derived the divergence operator as follows:

\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}

If the divergence of the velocity equals zero then:

\nabla \cdot \begin{pmatrix} W \\ V \\ U \end{pmatrix}\ =\frac{1-\mu^2}{r} \frac{\partial V}{\partial \mu} + \frac{1}{r \mu}\frac{\partial U}{\partial \lambda}-\frac{RT}{\sigma mg}\frac{\partial W}{\partial \sigma}=0

Which implies:

\frac{\partial W}{\partial \sigma}=\frac{(1-\mu^2)\sigma mg}{r RT} \frac{\partial V}{\partial \mu} + \frac{\sigma mg}{r \mu \sigma mg}\frac{\partial U}{\partial \lambda}


September 12, 2009 Posted by | GCM (General Circulation Models) | 1 Comment

Lagrangian Mechanics and The Heat Equation

I have been noticing a lot of discussion on Lucia’s blog about an energy Balance Climate Model.

Two Box Model: Algebra for ver. 1 test.
Two Box Model: Now assuming surface temperature are ‘mixed’ values.
Two Box Model: Rough idea how to obtain parameters.
Two Box Models & The 2nd Law of Thermodynamics

Which started based on a post by Tamino, called “Two Box”

The biggest criticism I’ve scene about these models is that there is very little physics which is used to form these models and consequently they won’t capture much of the nonlinear dynamics.

The basic premises of the model is that the time lags are associated with the heat capacity of the ocean. Consequently, the ocean and some of the atmosphere is partitioned into boxes and the model tries to determine how these boxes are coupled.

The end result though is the time legs end up being associated with eignvectors which, project onto both boxes.  Consequently the modes of the system do not coincide with the boxes.  Moreover, it is desirable to minimize the number of model parameters which need to be fit and therefore principles of physics, should be used to derive a more realistic “Two box model” or should I say more generally “Two mode model” as the modes do not coincide with the boxes.

The most simple equation with regards to the transfer of energy is the heat equation:

1)               \frac{\partial u}{\partial t} -k\left(\frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2}+\frac{\partial^2u}{\partial z^2}\right)=0

The heat equation is essentially a diffusion equation based on Brownian Motion. But we can adjust the constant K to try and account for other heat exchange processes.

In one  dimension the heat equation can be written as:

2)                 \frac{\partial u}{\partial t} -k\frac{\partial^2u}{\partial y^2}=0

In this equation the time derivative only depends on the spatial derivatives.  Keeping a mind on the dynamics, the higher the frequency that the climate forcing is, the more shallow it will penetrate into the ocean because the more ocean you include the greater the thermal inertia (heat capacity) and consequently the greater the damping.  Thus a good fit for the temperature response at a given frequency as a function of depth, may be an exponential function of the depth.

Thus, from this observation and because the mathematics are fairly simple consider a set of basis functions which the special component is exponentially decreasing from the surface.  That is let:

3)                u=\sum a_i exp(-\lambda_i y)

Plugging this into the heat equation one gets:

4)                \frac{du}{dt}=\sum \frac{d ( \ a_i)}{dt} exp(-\lambda_i )=k \lambda_i^2\sum a_i exp(-\lambda_i y)


The constant k(y) represents how quickly heat diffuses in the ocean at a given depth. Equating terms on the heat equation gives:

5)                k\lambda_i^2a_i-\frac{d (a_i)}{dt}=0

Now let there be a function c(y) , and let the quanity being diffused in the heat equation be the temperature. Then the energy constraint is given by:

6)               E=\int_0^{d} c(y)\sum a_i exp(-\lambda_i y)dy

This constraint makes one of the generalized coordinates a_i redundant.

Differentiating (5) with respect to time, gives:

7)             \frac{dE}{dt}=\int_0^{d} \sum c(y)\frac{d(a_i)}{dt} exp(-\lambda_i y)dy

\frac{dE}{dt}   is determined by the forcing. Now some notes, with regards to the coordinate system:

1)I suggest taking the as the redundant basis function, the one that has the smallest value \lambda .
2)Since this model is a linearization (or approximation), I suggest using as the temperature variable at a given depth, the amount by which it exceeds the average temperature at that depth.  That is model the temperature anomaly instead of the actual temperature.
3)Given two then we are looking at the change in the heat transfer from the mean instead of the actual heat transfer.
4)These are only suggestions.

Now the Lagrangian, is essentially a cost function.  In our cost function, we want to minimize the error in the derivatives (the square of equation (5) integrated over y) and as an additional constraint, I want to minimize the change in entropy with respect to y.

When a parcel of water rises because it is hotter then the surrounding and therefore less dense it will expand with constant entropy (adiabatic expansion) if it exchanges no heat with the surroundings. Therefore, the heat induced ocean currents will act in order to try and minimize the entropy change between adjacent regions. (Note the actual effect of this process is to increase entropy).

From thermodynamics:

8 ) dU=T\,dS-p\,dV


9) dS=(dU-p\,dV)/T=\frac{1}{T}(\frac{1}{C(y)}dT-p(y) \frac{dV}{dP}dP)

from which we can get:

10) \frac{dS}{dy}=\frac{1}{T}(\frac{1}{C(y)}\frac{dT}{dy}-p(y) \frac{dV}{dP}\frac{dP}{dy})

Substituting into (10) the expression for temperature (equation (3)) equation (10) gives:

11) \frac{dS}{dy}=-\frac{1}{T}(\frac{1}{C(y)}\sum a_i\frac{-1}{\lambda_i} exp(\lambda_i y)+p(y) \frac{dV}{dP}\frac{dP}{dy})

The entropy cost function is:

12) \int_o^d (\frac{dS}{dy})^2dy=\int_0^d\frac{1}{T^2}(\frac{1}{C(y)}\sum a_i\frac{1}{\lambda_i} exp(\lambda_i y)+p(y) \frac{dV}{dP}\frac{dP}{dy})^2dy

The Lagrangian is given by:

13) \mathcal L(a_1,...a_n, ...;\dot a_1,..,\dot a_n, ...;t)\,.= \alpha \int_o^d (\frac{dS}{dy})^2dy+\sum(k\lambda_i^2a_i-\frac{d (a_i)}{dt})^2

Where \alpha is a stiffness parameter for the entropy.

Equation 13) is subject to the constraints for some a_j  (These are the energy constraints mentioned above.

14) \frac{d(a_j)}{dt}=\frac{\frac{dE}{dt}-\int_0^{d} \sum_{i \neq j} c(y)\frac{d(a_i)}{dt} exp(-\lambda_i y)dy}{\int_0^{d}c(y)exp(-\lambda_j y)dy}

15) a_j=\frac{E-\int_0^{d} c(y)\sum_{i \neq j} a_i exp(-\lambda_i y)dy}{\int_0^{d} c(y) exp(-\lambda_i y)dy}

Now the Euler Lagrange Equation is used to obtain the dynamics (differential equations) from the Lagrangian

16) \left(\frac{\mathcal{L}(a_1,\ldots,a_n,\ldots,\dot a_1,\ldots,\dot a_n,\ldots;t)}{\partial a_i}\right) - \left(\frac{d}{dy} \frac{\mathcal{L}(a_1,\ldots,a_n, \ldots,;\dot a_1,\ldots,\dot a_n, \ldots;t)}{\partial \dot a_i}\right) = 0 \quad \textrm{for } i = 1, \dots, n, i \neq j

This will give a second order differential equation and it might be necessary to use some linear algebra to rearrange the equation. You can convert this into a first order differential equation by using the Hamiltonian form.  If only three basis  functions are used then one gets “two mode model”. There are no restrictions on the number of bais functions used.

September 1, 2009 Posted by | GCM (General Circulation Models) | 7 Comments

Coriolis Forces

A derivation  for coriolis forces can be found be found on the wikipedia page for fictitious forces.

In general for an accelerating reference frame in rectangular coordinates the factious forces are given by:

\mathbf{F}_{\mbox{fictitious}} =-m\ \mathbf{a}_{AB} -2m\ \sum_{j=1}^3 v_j \ \frac{d \mathbf{u}_j}{dt} - m\ \sum_{j=1}^3 x_j \ \frac{d^2 \mathbf{u}_j}{dt^2}\ .


he first term is the Coriolis force, the second term is the centrifugal force, and the third term is the Euler force. When the rate of rotation doesn’t change, as is typically the case for a planet, the Euler force is zero.

Looking specifically at the Coriolis force:

- 2 m \boldsymbol \omega  \times \mathbf{v}

which gives in (east-west,north-sourth, height) coordinates:

\boldsymbol{ \Omega} = \omega \begin{pmatrix} 0 \\ \cos \varphi \\ \sin \varphi \end{pmatrix}\ , \boldsymbol{ v} = \begin{pmatrix} v_e \\ v_n \\ v_u \end{pmatrix}\ ,

\boldsymbol{ a}_C =-2\boldsymbol{\Omega \times v}= 2\,\omega\, \begin{pmatrix} v_n \sin \varphi-v_u \cos \varphi \\ -v_e \sin \varphi \\ v_e \cos\varphi\end{pmatrix}\ .

Where \phi is the latitudinal coordinate (equator=zero latitude).

In general the cross product for a coordinate with orthonormal direction vectors is given by:

\boldsymbol{\Omega \times v} = \begin{vmatrix} \boldsymbol{i}&\boldsymbol{j}&\boldsymbol{k} \\ \Omega_x & \Omega_y & \Omega_z \\ v_x & v_y & v_z \end{vmatrix}\ = \begin{pmatrix} \Omega_y v_z - \Omega_z v_y \\ \Omega_z v_x - \Omega_x v_z \\ \Omega_x v_y - \Omega_y v_x \end{pmatrix}\ ,

since the basis direction vectors are orthogonal in hopkins and simmons coordinates write write:

\boldsymbol{\Omega \times v} = \begin{vmatrix} \boldsymbol{e_{\sigma}}&\boldsymbol{e_{\mu}}&\boldsymbol{e_{\lambda}} \\ \Omega_{\sigma} & \Omega_{\mu} & \Omega_{\lambda} \\ W & V & U \end{vmatrix}\ = \begin{pmatrix} \Omega_{\mu} U - \Omega_{\lambda} V \\ \Omega_{\lambda} W - \Omega_{\sigma} U \\ \Omega_{\sigma} V - \Omega_{\mu} W \end{pmatrix}\ ,

(note with regards to weather the system is right handed we can choose the direction of the logitudanal cordinate \lambda to make it right handed.)

Just to recall from the post (Hoskins and Simmons (1974) Coordinate System):

\mu = sin( \theta ) where theta is the latitude.
\sigma = pressure/P_* Where P_*  is the surface pressure and \sigma is the vertical coordinate.
\lambda is the longitude.


U is the longitudinal component of the velocity
V is the latitudinal component of the velocity
W is the vertical component of the veolicty (not used in Hopkins and Simmons 1974)

Now the angular velocity of the earth in Hopkins and Simmons is given by:

\boldsymbol{ \Omega} = \omega \begin{pmatrix} \cos \varphi \\ \sin \varphi \\ 0 \end{pmatrix}\ = \omega \begin{pmatrix} \mp \sqrt{1-\mu^2} \\ \mu \\ 0 \end{pmatrix}\, \boldsymbol{ v} = \begin{pmatrix} W \\ V \\ U \end{pmatrix}\ ,

Where the sign of \pm is positive for the northern hemisphere and negative for the southern hemisphere.


\boldsymbol{\Omega \times v} =  \begin{pmatrix} \mu U \\  \pm \sqrt{1-\mu^2} U \\ \mp \sqrt{1-\mu^2} V - \mu  W \end{pmatrix} \

Some comments:

The result obtained is essentialy the same result that one would get if, they took the(east west,north south, altitude) coordinate system and replaced \phi with asin(\mu) .

The only differences are the order and sign of the components.  These are the only differences because both coordinate contain the same unit vectors.  In my example of a Hopkins and Simon’s like coordinate system I used a different order for the components then was used in my example for the (east-west, north south altitude) coordinate system. This will effect the sign in the cross product.

I wrote the z component of the angular velocity as \mp \sqrt{1-\mu^2} to emphasize that the positive direction for the z component in Simpons coordinate system  is downward. However, the actual angular rotation of the earth in simons coordinate system still have a postive \mu component depending on the which direction is defined as positive for the longitudinal  coordinate.

The order which we specify the coordinates determines the right handedness of the coordinate system.  Therefore, righthandedness is not inherently a geometric property because it depends on the order of the coordinates. For instance, in standard Cartesian coordinates e_z \times e_y = -e_x

In our case the first coordinat,e \sigma , was specified in the downward direction, our second coordinate,  \mu , points south, now using the right hand  rule means that gives the positive direction for the third coordinate \lambda in the east direction.

It is for these reasons that differences can arrise, and therefore it is very important when doing cross products to clearly express the postive direction of the coordinate unit vectors and the order of the coordinates.

August 29, 2009 Posted by | GCM (General Circulation Models) | 2 Comments

Hoskins and Simmons (1974) Coordinate System

sHoskins and Simmons (1974) Coordinate System

Hoskin and Simon used the following coordinate system:

\mu = sin( \theta ) where theta is the latitude.

\sigma = pressure/P_* Where P_*  is the surface pressure and \sigma is the verticle coordinate.

\lambda is the longitude.

This coordinate system is also used in the following model which was recomended to me as a simple model:

As the complexity of general circulation models has been and still is growing considerably, it is not surprising that, for both education and research, models simpler than those comprehensive GCMs at the cutting edge of the development, are becoming more and more attractive. These medium complexity models do not simply enhance the climate model hierarchy. They support understanding atmospheric or climate phenomena by simplifying the system gradually to reveal
the key mechanisms. They also provide an ideal tool kit for students to be educated and to teach themselves, gaining practice in model building or modeling. Our aim is to provide such a model of intermediate complexity for the university environment: the PlanetSimulator. It
can be used for training the next GCM developers, to support scientists to understand climate processes, and to do fundamental research.

In order to use the coordinates given by Hoskins and Simmons in the Navier Stokes equations it is necissary to find the operations for Div Grad and Curl in this coordinate system. You can get a free PDF here:

Which shows how to derive these operations for any coordinate system. These derivations express the apove operations in terms of scaling quantities. The scaling quantitys are the magnitude of the derivative of the change in postion, with respect to a change in one of the coordinates in the new coordinate system.

In our case for the coordinate \mu the scalling quanity is:

h_{\mu}=\left| \frac{ \partial \vec r}{\partial \mu} \right| = \left| \frac{ \partial (r  \theta)}{\partial \mu} \right|

because when you change latitude the distance moved is proportional to the radius multiplied by the change in latitudinal (measured in radians).

Now expressing the latitude in terms of \mu
h_{\mu}=\left| \frac{ \partial (r \ sin^{-1}( \mu ))}{\partial \mu} \right| = \frac{r}{\sqrt{1-\mu^2}}

Now for the \lambda which is the longitude coordinate, h_{\lambda} depends on how you define zero latitude. The distance changed for a given change in longitude is the product of the change in logitude (measured in radians) multiplied by the distance to the earths axis of rotation. If zero lattitude is the center of the earth then this distance is r cos(\theta) if zero latitude is the north pole (typical in spherical coordinates) then this distance is r sin(\theta) .

We will see that taking zero latitude as the north pole will be simpler.

h_{\lambda}=\left| \frac{ \partial \vec r}{\partial \lambda} \right| = \left| \frac{ \partial (r  \lambda sin(\theta))}{\partial \lambda} \right|=r \ sin(\theta)

But \theta=sin^{-1}(\mu) Therefore:

h_{\lambda}=r \mu

The vertical coordinate \sigma is the ratio of the hydrostatic pressure to the hydrostatic pressure at sea level.

In Hoskins & Simmons the hydrostatic pressure used in the coordinate ignores the lapse rate. In some models variations on this are used. For instance in a model used by Judith Curry this coordinate is somewhat adjusted based on the surface topology of the earth. (see my post: hybird Coordinate System)

hybird coordinate system

In Curries model this is done to try and keep wind flow on the horizontal plane. Anyway, the equation for hydrostatic equilibrium is:

\frac{dp}{dh} = - \frac {mpg}{RT}


\frac{\partial \sigma}{\partial h}=\frac{\partial \sigma}{\partial P}\frac{\partial P}{\partial h}=\frac{1}{P_*}\frac{\partial P}{h}=-\frac{1}{P_*} \frac {mpg}{RT} =-\sigma \frac{mg}{RT}

Hmmmmmmm: From what I’ve done the only way I can see to make the constants \frac{mg}{RT} disappear is to set T=\frac{mg}{R}

Here is what Hoskins and Simmons have to say about temperature scaling and such.

all non-dimensionalized using as length scale the radius of the planet, a, ; as time scale the reciprocal of its angular velocity \Omega ; temperature scale a^2\Omega^2/R (R as the gas constant) and pressure scale P_o=1 Bar.

In T_o was taken to be 250K. Some numbers of rellevence:

Gas constant of air: 287.058 J kg−1 K−1

and for the molar mass of air:

Standard temperature =273.16k

Source for the molar mass of air:

Standard pressure=1 Atm=1013.250*10^2 Pa.
The molar mass of air at standard temperature and pressure is 28.964*10^-3 kg/mole. Molar

Standard Gravity: 9.80665 m/s2

August 27, 2009 Posted by | GCM (General Circulation Models), Math | 3 Comments

A Comparison of GCMs and Some History (Two Links)

I found two links which seem to be a good place to start with understanding the differences between the current GCM, and the history of GCM.

The GCM-Reality Intercomparison
Project for SPARC (GRIPS):
Scientific Issues and Initial Results
S. Pawson,a,b,c K. Kodera,d K. Hamilton,e T. G. Shepherd,f S. R. Beagley,g B. A. Boville,h
J. D. Farrara,i T. D. A. Fairlie,j A. Kitoh,d W. A. Lahoz,k U. Langematz,c E. Manzini,l
D. H. Rind,m A. A. Scaife,n K. Shibata,e P. Simon,o R. Swinbank,n L. Takacs,p
R. J. Wilson,e J. A. Al-Saadi,q M. Amodei,o M. Chiba,r L. Coy,p J. de Grandpré,g
R. S. Eckman,q M. Fiorino,s,t W. L. Grose,q H. Koide,d J. N. Koshyk,f D. Li,k,n
J. Lerner,m J. D. Mahlman,e N. A. McFarlane,u C. R. Mechoso,i A. Molod,p
A. O’Neill,k R. B. Pierce,q W. J. Randel,h R. B. Rood,b and F. Wuh
To investigate the effects of the middle atmosphere on climate, the World Climate Research Programme is supporting the project “Stratospheric Processes and their Role in Climate” (SPARC). A  entral theme of SPARC, to examine model simulations of the coupled troposphere–middle atmosphere system, is being  erformed through the initiative called GRIPS (GCM-Reality Intercomparison Project for SPARC). In this paper, an overview of the objectives of GRIPS is given. Initial activities include an assessment of the performance of middle atmosphere climate models, and preliminary results from this evaluation are presented here. It is shown that although all 13 models evaluated represent most major features of the mean atmospheric state, there are deficiencies in the magnitude and location of the features, which cannot easily be traced to the formulation (resolution or the parameterizations included) of the models. Most models show a cold bias in all locations, apart from the tropical tropopause region where they can be either too warm or too cold. The strengths and locations of the major jets are often misrepresented in the models. Looking at three-dimensional fields reveals, for some models, more severe deficiencies in the magnitude and positioning of the dominant structures (such as the Aleutian high in the stratosphere), although undersampling might explain some of these differences from observations. All the models have shortcomings in their simulations of the present-day climate, which might limit the accuracy of predictions of the climate response to ozone change and other anomalous forcing.

1965-75: Spread of AGCMs

By 1965, three groups in the United States had established ongoing efforts in general circulation modeling:

  • Geophysical Fluid Dynamics Laboratory
  • UCLA Dept. of Meteorology
  • National Center for Atmospheric Research

At this point, AGCMs and modeling techniques began to spread by a variety of means. Commonly, new modeling groups began with some version of another group’s model. Some new groups were started by post-docs or graduate students from one of the three original AGCM groups. Others built new models from scratch. An AGCM family tree offers a visual map of these relationships, with links to some of the models and modeling groups.

I plan to have more to say on this later.

August 14, 2009 Posted by | GCM (General Circulation Models) | 1 Comment