## Matlab, Octave and Climate Simulation

Whilst looking about for a simple Mathematica climate simulation, I ran across some .m files dedicated to just that purpose. As I recognized .m as the extension for a Mathematica package, I snapped ’em up and tried installing them. Total, abject failure…

After a bit of research, I figured out that .m was also the extension for Matlab m-files. Since I had an old copy of Matlab from my misspent youth in mechanical engineering, I tried ’em out. Yay! They worked!

On the one hand, I’ve put a lot of time into learning Mathematica and it is a good tool for publishing good-looking mathematical formulae. On the other hand, I’ve found a lot more useful climatic tools for Matlab, like this one and this one here. On the gripping hand, I found a freely available open-source app called Octave that is able to run almost all the Matlab scripts I have found. Please, run out and get Octave if you don’t have Mathematica. If you do have Mathematica, run out and get Octave anyway. I’ve found a few holes in Octave’s coverage of Matlab procedures, mostly in charting tools, but it’s functionality is pretty damn good.

I’ve started making some modifications to temp_gradient.m to increase its flexibility for my purposes. To use my modification, copy the text in temp_gradient2.m, paste it into a new file named temp_gradient2.m. Make sure that the folder where you store the new m-file is either the working folder or on the path for Matlab or Octave.

For Earth, I’m using an obliquity, or axial tilt, of 23.5°, an eccentricity of 0.0172, a latitude of perihelion of 281.37°, an albedo of 0.3, or 30%, and a Solar constant of 1365 W/m2. The highest mean annual temperature calculates to 304.9 K(31.7° C, 89.1° F), and the lowest temperature is 244.7 K(-28.4° C, -19.2° F). The mean temperature is lower than the actual surface temperature of 288 K usually quoted for the Earth, but that makes sense, as the mean is weighting temperatures at the north and south poles exactly the same as it weighs temperatures at the equator even though the equatorial band covers a much, much larger surface area.

The results of applying my modified procedure with Earthlike parameters.

My modified model with obliquity and eccentricity set to zero.

The weaknesses in this model show up pretty quickly if you set the obliquity and eccentricity to zero. Notice that the temperature goes to absolute zero at the poles. This is unrealistic even for a rock ball with no atmosphere.

The same problem exists in the unmodified model. The problem is that there is no spatial or temporal exchange of heat. In the real world, if a very hot place is close to a very cold place heat will flow from the hot end to the cold end. Thus the hot end will be cooler and the cold end will be warmed up. The exchange can be facilitated efficiently by convection in a dense atmosphere, but it also occurs, although much less effectively, due to conduction through the rocky surface of a planet, even without an atmosphere.

The complexity of a procedure required to handle this heat exchange realistically can be pretty doggone extreme. General circulation models can be immensely complex, with realistic modeling of ocean currents air pressure and winds, cloud cover, evaporation and precipitaion. Even with a good fast computer and well-optimized programming, it can take a long time to converge on a realistic global climate pattern without good starting conditions.

In time, I’m hoping to come up with a more middle-of-the-road solution, using an energy balance model with heat exchange. I’ll likely base my model on the following formula for heat exchange.

Haberle et al(1996)

Where ps is the surface pressure of the atmosphere, cp is the specific heat capacity of the atmosphere(at constant pressure), g is the surface gravity of the planet and τ is a transport time scale for the atmosphere. Of these, all except τ can be readily figured out. At the stage of world building where we’re trying to figure out temperature distributions across the planet surface, we should know the surface pressure of the atmosphere and the planetary surface gravity. For cp a value of about 1000 J/kg/K will suffice for dry air(this value was in fact quoted from an engineering site and is the same value as used by Haberle et al for a Venus-like atmosphere, so there’s some wiggle room there). Tla and Tda represent the temperatures of the light-side atmosphere and the dark side atmosphere.

Haberle et al(1996)

The image to the left came from the article from which I sourced the equation. It may clarify how it works in practice. As you can see it is used for a two point energy balance model with the assumption that one point is on the light-side and the other on the dark-side of a face locked planet. In the case of what I am trying to do, a varying heat input S would exist for each latitude band on the planet. Except for the north and south polar bands, there would also be two (ps cp)/(g t) * (Tthis,a – T other,a), one representing exchange with the band to the south and one representing exchange with the band to the north. Because this is used on a much smaller spatial scale(by bands of latitude rather than by hemisphere) the time factor tau would also be smaller. Tau should likely vary with planetary radius all else being equal.

I’m currently trying to adapt the matlab code available here to my purposes. A fuller description of the code and how it works is given here. I still haven’t quite figured out the code itself, but I think instead of figuring out the rate of heat transfer in and out of each band it has an averaging stage where it biases temperatures towards a one-d average determined for the planet as a whole. For most of our purposes, this method should be completely sufficient.

In any case I’ve put a lot of time into this, I hope, useful little digression. I will try to get myself back onto schedule now.

Thank you for your tolerance,

The Astrographer

This entry was posted in Links, Planetary Stuff, Science! and tagged , , , , , , , , , , , , , . Bookmark the permalink.