Climate 3: Insolation and Black-Body Radiation

2-year radiation simulation (first image is corresponding elevation map)

The only way for energy to enter or leave my climate simulation is through radiation. As is discussed more formally below, energy enters via solar radiation, and energy leaves via black-body radiation. I’m starting with this because these two phenomena drive the entire climate system. Without them, winds and weather can’t exist. With only radiation, and without air or water, the simulation will already give recognizable results.

The gif above is the result. The gif shows the ground temperature over 2 simulated years. The timestep was 1 hr and the pictures are 500 hrs apart. Looks like a success to me.

For an in-depth discussion of the results, use the Table of Contents to skip to the Results section. Otherwise, strap in for an in-depth discussion of how this simulation was designed.

Table of Contents

Goals

These are straightforward. Seasons. Day/night. Land vs water temperature differences. Eventually I’d like further albedo differences (e.g. vegetation cover, snow and ice) and effects from clouds and humidity, but that’s a ways down the road.

Ideally, I’d like to generate at least 4 maps, one for each season, which I could rotate between while I’m simulating stuff. I could generate more, say, 12, but the exact number isn’t important. The point is, I want to have a changing seasonal climate which affects the rest of my simulation.

Overall Energy Balance

The conservation of energy gives us the following equation for the energy balance over any system:

Energy In – Energy Out + Energy Created – Energy Destroyed = Accumulation

The system boundary for this climate simulation is the boundary between the [planet + atmosphere] and space. Energy In is solar radiation received by the planet and Energy Out is black-body radiation emitted by the planet.

I’m not going to consider any Energy Created or Energy Destroyed terms. Note that these terms don’t literally mean creating or destroying energy. They represent energy conversions between types of energy the system model tracks and types it doesn’t consider. An example of this is energy generated by radioactive decay in the planet — my simulation tracks heat energy but not radioactive isotopes, so heat generated by radioactive decay would appear as an Energy Created term in the equation. These terms seem unnecessary for the level of simulation I’m working with.

Accumulation is any energy stored in the system. In this simulation, accumulation will consist only of the temperature of the planet surface and atmosphere. When Energy In > Energy Out, the temperature of the planet surface and atmosphere will generally go up. If the simulation ever approaches perfect steady-state, Energy In will equal Energy Out and the temperatures of everything will be constant. So the overall energy balance equation becomes:

Energy In – Energy Out = Accumulation

Or:

Insolation – Black-Body Radiation = Change in Temperature

For convenience, all my equations will be in terms of energy per area, or flux. Since I’m using a tile map, the area is constant and there’s no reason to constantly multiply and divide by the area.

This overall energy balance equation will hold true even after we add winds and other types of convection. All of those other mechanisms will only move energy around within the system (planet + atmosphere), which won’t change the energy balance over the system boundaries.

Now I’ll dive into each term, starting with the most complex: Insolation.

Energy In: Insolation

Note: If you are implementing this, I recommend initially using simplified insolation equations. For reference, on my first iteration of my climate simulator, I used a simple insolation function based on latitude (no seasons, day/night cycle, axis tilt, or spherical mapping). After I the rest of my climate simulator worked, I went back and rewrote the insolation equations to incorporate these phenomena.

The Sun emits energy in the form of electromagnetic radiation. Since the sun is so far from Earth, we can approximate all that electromagnetic radiation as coming in at a single angle. Sunlight hits stuff, and is absorbed, and thus heats stuff up. As shown below, the amount the sun heats the ground is dependent on the angle of incidence.

Parts of the Earth perpendicular to the Sun receive a certain amount of sunlight over a certain area, and parts of the Earth at a steep angle receive the same amount of sunlight spread out over a larger area.

For the parts of Earth perfectly perpendicular to the Sun, the Earth receives (neglecting albedo) a maximum flux $I_{max}$ of 1386 W/m^2. What about non-perpendicular sections of Earth? We need to find the dot product between the surface normal vector of Earth and the negative of the sunlight vector. Here’s a diagram illustrating the planet’s surface normal vector in several places. If we fix Earth in place, we can model the Sun as revolving around the Earth (taking us back to the 1500s).

Hold up. Earth (like most planets) is spherical. How does this work with the cylindrical map I’ve been using?

Seasons are created by 1) sphericity and 2) a tilted axis (they have nothing to do with distance to the Sun).

My simulation thus far has used a cylindrical map because it allows me use rectangular coordinates with a hardcoded wraparound boundary. If you imagine replacing Earth in the image above with a tilted cylinder, the angle of incidence will be the same for all latitudes. Thus a cylinder will have no seasons.

So, just for insolation, I will map the cylindrical map onto a sphere.

Spherical Transformation and Calculation of Sunlight Intensity

This transformation is effectively a reversed Mercator map projection, and comes with the same distortions (the polar regions are expanded in the cylindrical version), which I’ll have to account for later.

Remember that in the Mercator projection, the sphere is the true representation, whereas in this project, the cylinder is the true spacial representation.

Note: I use the physics version of spherical coordinates.

The transformation requires developing a function for the surface normal vector in spherical coordinates from the 2D map, which will involve mapping $x \rightarrow \theta$ and $y \rightarrow \phi$ and defining the sunlight vector. Both vectors will be functions of the time of day and time of year (due to the rotation of the planet and revolution around the sun).

The image below shows the important stuff. The axis is tilted by some angle. The planet rotates on its axis. The sunlight vector is horizontal. In reality, the planet revolves around the sun, but for convenience I’ll have the sun revolve around the planet. At any given latitude, there is a surface normal vector which points from the center of the planet to the point on the surface.

Over the course of one day, the planet rotates once on its axis, and the surface normal vector traces a circle at its latitude. Over the course of one year, the sunlight vector travels in a circle around the planet. Remember that the insolation is calculated from the dot product between these two vectors (with a minimum of 0).

Depiction of the surface normal vector and insolation vector

While this transformation is accurate already, calculating the surface normal vectors will be difficult due to the axis tilt. The transformation would be much easier if the planet’s axis aligned with the spherical coordinate vertical axis. I can achieve this by tilting the planet “upright” and tilting the sunlight vector a corresponding amount.

Depiction of insolation vector system, shifted

Now three of the four coordinate transformations are simple linear functions, as you’ll see shortly. I recommend again trying to visualize the movement over the vectors over time. An animation here would be helpful, but what do you think this is? Red Blob Games? Animate in your head.

Here are the coordinate mappings.

Sunlight vector:

$\theta_s : [0, secondsInDay) \rightarrow [0, 2 \pi)$

$\phi_s : [0, daysInYear) \rightarrow [\pi /2 - axisTilt, \pi /2 + axisTilt)$

The sunlight vector component $\phi_s$ is the one that requires a nonlinear function. If you visualize the sun revolving around the planet in the adjusted representation, the sunlight vector will move between $\pi /2 - axisTilt$ and $\pi /2 + axisTilt$ (which are also the tropics) in a sinusoidal path. So the full equation for $\phi_s$ is:

$\phi_s = \pi /2 + axisTilt * sin(2 \pi * dayOfYear / daysInYear)$

Surface normal vector:

$\theta_n : [0, xMax) \rightarrow [0, 2 \pi)$

$\phi_n : [0, yMax) \rightarrow [0, \pi)$

Easy peasy. Now at each simulation tick I can calculate the sunlight vector based on the time. Then at each tile I can calculate the tile surface normal vector and take the dot product. The dot product is best done in Cartesian coordinates. Note that here I am only changing coordinate systems. The map is still in the spherical representation for the surface normal vector.

Transformation of a vector from spherical to Cartesian coordinates:

$<X, Y, Z> = <\rho sin(\phi ) cos(\theta ), \rho sin(\phi ) sin(\theta ), \rho cos(\phi )>$

Note that I’m using unit vectors so $\rho = 1$ . The magnitude will be included as the $I_{max}=1386 \frac{W}{m^2}$ later.

Then take the dot product between the surface normal and the sunlight vectors to find the insolation flux at a given point:

$Insolation = I_{max} * (X_{sunlight} * X_{surface} + Y_{sunlight} * Y_{surface} + Z_{sunlight} * Z_{surface})$

And that’s insolation per area with seasonal and daily variation.

Energy Out: Black-Body Radiation (Stefan-Boltzmann)

Every bit of matter above absolute zero emits energy in the form of electromagnetic radiation (including you, you’re hot stuff).

The wavelength of the radiation is correlated with the temperature of the object. This is why room-temperature objects don’t glow, lava glows red, and the hottest stars glow blue.

I don’t care about the wavelength of light emitted by Earth, only the total amount of energy. The amount of energy emitted from a perfect black-body at a temperature T is given by the Stefan-Boltzmann law.

$E = \sigma T^4$

where $\sigma$ is the Stefan-Boltzmann constant with value $5.670367 x 10^{-8} \frac{W}{m^2 K^4}$

I’ll add a unitless term emissivity $\epsilon$ which is a measure of how much energy escapes through the atmosphere. Greenhouse gases warm the planet by reducing the emissivity (no atmosphere would correspond to an emissivity of 1). I doubt I’ll ever use a variable emissivity, but including it here is easy and leaves the possibility open.

$E = \epsilon \sigma T^4$

Since the energy emitted is already in energy per area, and it’s not affected by the map representation, there are no further transformations necessary.

Note: The atmosphere has several additional layers of complexity which are considered by full climate models. The atmosphere is also made of matter warmer than absolute zero, so the air itself is also emitting radiation. A full model integrates the black-body radiation at every altitude (models can use dozens of altitude layers) of the atmosphere and also how much radiation is absorbed or reflected by all above atmosphere layers. It’s convenient to approximate that total escaping energy as black-body radiation from only the ground.

Accumulation: Change in Temperature

The amount of heat energy in an object is directly related to its temperature.

An important coefficient involved in transferring heat is the heat capacity of a substance. The heat capacity is the amount of energy needed to raise the temperature of 1 g of the substance by 1°C.

Substances with high heat capacities can absorb a lot of heat without changing temperature much. In this simulation, we’ll use three heat capacities: dirt, water, and air.

$E = m * c_v * \Delta T$

$E = \rho * V * c_v * \Delta T$

Note: There are two different heat capacities, $c_p$ and $c_v$ . The first is for heating at constant pressure, the second for constant volume. Since solids and liquids are incompressible, for them, $c_p = c_v$ . The difference only matters for gases. I’ll hold the volume constant because my map tiles hav constant volume, so I’ll use $c_v$ and let the pressure vary during the heating of air.

Setting the Temperature

The equations so far are enough to make a complete simulation that approaches steady state with seasons and a day-night cycle. However, there’s an issue: what steady state does the simulation currently approach?

Since I used Earth’s insolation numbers, one might think the simulation should approach an average temperature similar to Earth’s, which is about 10-15 C or 283-288 K. However, the simulation will approach a much lower temperature. Why?

One of the main distortions of the Mercator representation is that the polar regions are much larger on the cylinder (the other has to do with direction).

Okay, so I could bump up the solar intensity $I_{max}$ so that the average temperature is similar to Earth’s. However this would cause the equator to become significantly hotter than Earth’s. I think I would rather have oversized polar regions than unearthly high equatorial temperatures.

So what’s the problem? The problem is that I have to set an initial temperature “guess” for all the layers in my simulation. I want the guess to be as close as possible to the final steady state value in order to minimize the amount of time it takes for my simulation to reach steady state. Additionally, if my “guess” is too far off, it will induce major instabilities in my atmospheric simulation later on.

I would rather flip this around. Instead of guessing what average planet temperature a certain value of $I_{max}$ will give me, it would be convenient to set the average temperature of the planet and then calculate the necessary $I_{max}$ .

Restated: I need to calculate the $I_{max}$ necessary to set the average temperature of the planet to a certain temperature $T_{sp}$ (sp stands for set point).

Let’s see how to calculate this for a spherical Earth-like planet first, then figure out how that transfers to my cylindrical world.

Earth insolation from The Open University

From the above image, you can see that the total amount of insolation received by Earth is equal to $I_{max}$ * (area of profile). The profile is equal to the total area perpendicular to the incoming sunlight (another way to think of this is the area of the shadow cast by the planet).

How could I calculate the area of the shadow cast by my cylindrical world?

I need to integrate my dot product. Currently, it’ll be at least a triple integral (with respect to x, y, and time) with a lot of trig functions. I can simplify it by stopping the rotation and revolution (they affect insolation distribution, not total insolation), which makes t constant. I’ll arbitrarily pick t = 0 because it makes the sunlight vector very convenient.

Sunlight vector (Cartesian): $< 1, 0, 0>$

Surface normal vector (Cartesian): $<sin\left( \frac{2\pi}{x_{max}}x \right)cos\left( \frac{\pi}{y_{max}}y \right), sin\left( \frac{2\pi}{x_{max}}x \right)sin\left( \frac{\pi}{y_{max}}y \right), cos\left( \frac{2\pi}{x_{max}}x \right)>$

Taking the dot product:

$InsolationFlux (x,y) = sin\left( \frac{2\pi}{x_{max}}x \right)cos\left( \frac{\pi}{y_{max}}y \right)$

Integrating over the side of the planet exposed to the Sun, the result is:

$totalInsolation = I_{max} * \frac{xMax * yMax}{2\pi ^2}$

Note: This transformation actually can be solved graphically. Thinking back to the last image — what is the area of the shadow cast by my “planet?” — I start with a rectangular map, take half of it (since half of it is day, half night), then bend it into a half circle in both the X and Y directions. The “bending” requires multiplying by 1/pi each time, as it’s mapping a circle’s half-circumference to its diameter.

From the relationship between the totalInsolation and the $I_{max}$ , I can derive the $I_{max}$ from the setPoint temperature from the energy balance equation at steady state (no accumulation so no change in temperature).

$totalInsolation - blackbodyRadiation = temperatureChange = 0$

$totalInsolation = blackbodyRadiation$

$I_{max} \frac{xMax * yMax}{2\pi ^2} = xMax * yMax * \epsilon \sigma T^4$

$I_{max} = 2\epsilon \sigma T^4 \pi^2$

Unfortunately, this formula does not work in practice. This equation assumes the temperature is uniform throughout the map. Because the black-body radiation isn’t linear with respect to temperature (there is a $T^4$ term), a map with a temperature range distributed about the average will have greater radiative emission than a map with uniform temperature. Therefore I end up needing a lot more (about 30% more) insolation than this equation indicates.

To account for the temperature distribution analytically, I’d have to integrate over time, which doesn’t sound like fun. Instead I’ll just use 1.3x the formula as my starting guess and empirically evaluate and update the necessary insolation every tick.

Results

The gif shows the ground temperature over 2 simulated years. The timestep is 1 hr and the images are 500 hrs apart. The temperature range is approximately -60 to 25 C. It won’t be too difficult to shift the temperatures up (this gif was taken before I finished the “Setting the Temperature” section).

When I actually use these climate maps, I’ll linearly interpolate between the grid points to the same tile resolution as the elevation maps (800×400 instead of 50×25). However, for debugging, it’s much more useful for me to view maps without interpolated data.

Overall, it looks like I’m on track to generate seasonal maps and save them.

Some other interesting points:

The land both heats and cools faster than the oceans (as intended). As the “hot band” moves north (approaching northern summer), the northern continents are warmer than the surrounding oceans. As the hot band moves south (approaching northern winter), the northern continents are cooler than the surrounding oceans.
During the first year, the northern hemisphere is warmer than the southern hemisphere. This is because the northern hemisphere’s summer occurs first, so at any point in time except exactly at the beginning of each year, on average, the northern hemisphere has received more insolation than the southern hemisphere. For example, at the end of northern summer 2, the northern hemisphere has experienced 2 summers and 1 winter, whereas the southern hemisphere has experiences 1 summer and 2 winters. This should even out over years of simulation until the difference is negligible.
Day/night is not visible in this gif. The images are 500 hrs apart, so the cycle wouldn’t be visible in this gif, but in any given image I would expect there to be more horizontal variation in temperature. I wonder if the day/night temperature difference is just smaller than the latitude temperature difference or if there is an actual issue with my day/night cycle. I’ll have to investigate.

TODO

When vegetation comes around, I’ll incorporate variable land albedos. I’m also interested in modeling snow/ice cover and cloud cover, but those are stretch goals. I’m pretty happy with the insolation/radiation section as is.

Since terrain generation, I’ve been stuck on the ground. With the overall system energy balance stable and ready to go, it’s time to shatter my earthy shackles and enter the capricious domain of air modeling.