Earth Cubed

Distributed Climate Science and Computing

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