Playing with Sunshine

Harys Dalvi

April 2021

Simulation

Average surface temperature: 250 K
Month: March

Help

Hover over a square in the surface temperature or atmospheric temperature layer to see its temperature.

Day/night shows approximately where in the world it is day and where it is night at a moment in the simulation.

Water proportion shows my estimates for the proportion of each square that is water as opposed to land. Due to water's high specific heat, squares with more water tend to resist changes in temperature.

The greenhouse gas level represents the fraction of heat released by the Earth's surface that is later absorbed by the atmosphere. In our present world, this is about 0.77, or 77%. Greenhouse gases tend to make the Earth warmer. [2]

Albedo represents the fraction of the sun's radiation that the Earth reflects rather than absorbing. This is a consequence of the Earth's reflectivity, for example ice that reflects sunlight. In our present world, the albedo is about 0.28, or 28%. Albedo tends to make the Earth colder. [2]

You may notice a gradual oscillation between north and south due to seasons. (Try max speed.) You may also notice an oscilation from east to west due to day and night. (Try speed 1.)

There are many issues associated with cylindrical projections and their inaccurate representation of sizes near the poles. However, I used the equirectangular projection in this simulation because it corresponds best to the spherical/polar coordinate system and helped make the math easy (easier?).

The Science Behind This

Limitations

To make this simulation, I used lots of simplifying assumptions. Here are a few of them to start. The rest will be explained where they are relevant.

The Earth is a perfect sphere of radius 6370 km and mass 5.972 × 10²⁴ kg. The 288 squares of the map total to 0.0005% of this mass, with the rest being contained in the earth's core and mantle or parts of the crust which are assumed unaccessible to the sun and are not part of this simulation. This mass is shared among the squares proportional to the true surface area of each “square”.
The Earth's orbit around the sun is circular with a period of 365.2422 days.
Each "square" on the map (15 degrees of latitude, 15 degrees of longitude) has a constant surface temperature, and a separate constant atmospheric temperature.
The atmosphere is of constant composition, and is 80 km thick.
The initial temperature of the Earth is 275 K on the surface, 250 K in the atmosphere.

There are also some crucial factors that this simulation does not take into account.

Clouds
Ocean currents
Winds

Although it is far from perfect, this simulation is an interesting rough model for the Earth's climate not only globally but also regionally.

Heat Flux

This simulation is primarily based on the Stefan-Boltzmann law and secondarily based on Fourier's law, with the first being much more important as you will see. The Stefan-Boltzmann law states $$\varphi = \sigma T^4$$ Where $\varphi$ is the radiation per unit area emitted by a black body, $T$ is the temperature of that black body, and $\sigma$ is the Stefan-Boltzmann constant, equal to about 5.670374419 × 10^-8 W/m²K⁴. A perfect black body will absorb all radiation that hits it. The Earth is clearly not a perfect black body, as you can tell because it isn't black (meaning it is reflecting rather than absorbing light of some colors). However, for the purposes of this simulation, we can use the Stefan-Boltzmann law for a black body because the infrared wavelengths that make up most of the Earth's radiation emissions are emitted almost like a perfect black body [2].

We then have heat coming in from the sun. The sun heats up the Earth with a heat flux $F_S$ equal to 1370 W/m² However, of this heat, some fraction $a$ is reflected to space rather than absorbed because of the Earth's albedo, leaving $F_S(1-a) $ for Earth. [2].

Now let's consider the greenhouse effect. We will add a layer of atmosphere above the Earth. This atmospheric layer absorbs a fraction $f$ of the heat per area $\sigma T^4$ emitted by the Earth. To make things clear, let's call $T_E$ the temperature of the Earth, and $T_A$ the temperature of the atmosphere. So the heat flux that the atmosphere absorbs from the Earth is $f \sigma T_E^4 $ (with the other $(1-f) \sigma T_E^4$ going to space). This fraction $f$ is based on the amount of greenhouse gases in the atmosphere [2].

Just like Earth, the atmosphere will emit radiation. However, the atmosphere cannot be modeled as a black body. Additionally, the atmosphere emits radiation in two directions: up (to space) and down (to Earth). The heat flux in each direcion is not $\sigma T_A^4$, as you would expect for a black body, but $f \sigma T_A^4$ [2].

To summarize what we have so far, let's follow the path of the heat from the sun.

A heat flux $F_S$ comes in from the sun. Note that unlike the other heat flux values, to get the actual heat from the sun, we multiply not by the surface area of the Earth but by the area perpendicular to the sun's rays.
Of this heat flux $F_S$, a fraction $a$ is reflected back to space due to the Earth's albedo. $F_S (1-a)$ remains, and this goes to heat up the Earth.
The Earth emits a heat flux $\sigma T_E^4$. Of this, a fraction $f$ is absorbed by the atmosphere, while the other $(1-f) \sigma T_E^4$ is lost to space.
The atmosphere radiates a heat flux $f \sigma T_A^4$ down to Earth, and an equal heat flux up into space.

Note that the ultimate source of all heat in this cycle is the sun, from step 1. Heat can escape to space in steps 2, 3, and 4. Looking at steps 3 and 4, we see some amount of heat is "bouncing" between the Earth and the atmosphere rather than being radiated into space. This is the greenhouse effect.

Now we can see the general direction of the flow of heat, but not the magnitude. In addition, we still haven't considered the different squares on the map. Since heat flux is heat per unit area, in general we can find the rate of heat flow by multiplying the heat flux by the area. The question is, which area? This is where some tricky math comes in.

Areas

First we consider the heat flux $F_S$ coming from the sun's rays. The rate of heat flow will depend on this heat flux as well as the area perpendicular to the sun's rays. For the entire Earth, this perpendicular cross-sectional area is simply $\pi R^2$, where $R$ is the radius of Earth [2]. However, this simulation is not looking at the whole Earth as one, but is splitting it into 288 squares. Therefore, we have to calculate the perpendicular area of each square. The perpendicular area for each square is

$$R^2 \cdot \frac{2 \sin (\delta) (\lambda_1 - \lambda_2) (\cos^2 \theta_2 - \cos^2 \theta_1) + \cos(\delta) (\sin \lambda_1 - \sin \lambda_2) (\sin (2 \theta_2) - 2 \theta_2 - \sin(2 \theta_1) + 2 \theta_1)}{4}$$

where $\lambda_2$ and $\lambda_1$ are the east and west longitudes of the square relative not to the prime meridian, but to the longitude where it is solar noon; $\theta_2$ and $\theta_1$ are the north and south latitudes of the square relative to the north pole (so $\theta = \pi - \phi$ in radians); $\delta$ is the declination; and $R$ is the radius of the Earth. The declination is a value that varies sinusoidally between 23.5° and -23.5° (from the axial tilt of the Earth) to represent either the northern or southern hemisphere tilting towards the sun.

Explanation of this result

We can solve this by using a surface integral. We can model the sun's rays as a vector field

$$\mathbf{F_S} = \langle 0, \cos \delta, \sin \delta \rangle$$

and we can model the Earth as a surface parametrized by

$$\mathbf{r} = \langle R \sin \varphi \sin \lambda, R \sin \varphi \cos \lambda, R \cos \varphi \rangle.$$

Now to find the area $A_{\perp}$ of a given square perpendicular to the sun, we take the section of the sphere with $\theta$ from $\theta_1$ to $\theta_2$, and similarly for $\lambda$. If we call this section a surface $S$, we can take the surface integral

$$A_{\perp} = \iint \mathbf{F_S} \cdot d \mathbf{S} = \int_{\lambda_1}^{\lambda_2} \int_{\theta_1}^{\theta_2} \mathbf{F_S} \cdot \frac{\mathbf{r}}{|\mathbf{r}|} \cdot R^2 \sin \theta \ d \theta \, d \lambda.$$

Unfortunately, this integral turns out to be extremely tedious when you plug in values. I won't go through it here, but if you like torturing yourself with this stuff, study multivariable calculus.

If this value is negative, it means that this square is not in the sunlight, and so will receive no heat from the sun. Notice that at higher latitudes, this area becomes smaller because $\theta \approx 0$, $\sin (2 \theta) \approx 2 \theta$ near 90° N, and 90° S is similar with a phase shift. This means that areas farther from the equator will be colder.

Now let's consider the heat flux emitted by the Earth and the heat flux emitted by the atmosphere. For these, the area is equal to the surface area of the part of the spherical Earth on this section of the map. This area is $$\frac{1}{12} \pi R^2 (\sin \phi_2 - \sin \phi_1)$$

Explanation of this result

We can imagine the square on the map as a solid of rotation. Think of a circle with the radius of the Earth and equation $x^2 + y^2 = R^2$. From this, we get $y = \sqrt{R^2-x^2}$. We can take the section of this line between latitudes $\phi_1$ and $\phi_2$ and then revolve it 360° around the $x$-axis; this will be the total surface area of all squares between latitudes $\phi_1$ and $\phi_2$. Since there are 24 of these squares, we must divide by 24 to get the area of just one square. The formula for surface area of a solid revolved around the $x$-axis is [6] $$2 \pi \int y \sqrt{1 + \bigg( \frac{dy}{dx} \bigg)^2} \, dx$$ Since we have $y = \sqrt{R^2-x^2}$, we find $$\frac{dy}{dx} = \frac{-x}{\sqrt{R^2-x^2}}$$ Plugging in for $y$ and $\frac{dy}{dx}$, and dividing by 24 to find the area of just one square, we get

$$A = \frac{\pi}{12} \int_{R \sin \phi_1}^{R \sin \phi_2} \sqrt{R^2-x^2} \sqrt{1 + \frac{x^2}{R^2-x^2}} \, dx$$

Interestingly, this becomes the integral of a constant and does not depend on $x$. $$A = \frac{\pi}{12} \int_{R \sin \phi_1}^{R \sin \phi_2} \sqrt{R^2-x^2+x^2} \, dx$$ $$A = \frac{\pi}{12} R^2 (\sin \phi_2 - \sin \phi_1)$$

Now that we have all the areas, we can combine these with the heat flux values to get the rate of heat flow into and out of the Earth and the atmosphere for each map square. This is the core of the simulation.

Fourier's Law

Fourier's law describes heat conductivity, the motion of heat. The equation states $$\varphi = -k \nabla T$$ where $\varphi$ is the heat flux, $k$ is the thermal conductivity of the material, and $\nabla T$ is the gradient of temperature, a vector measure of the direction in which temperature is increasing and by how much. Since $\varphi = \frac{dQ/dt}{A}$ by definition, where $\frac{dQ}{dt}$ is the rate of heat flow and $A$ is the area, we can write $$\frac{dQ}{dt} = -k A | \nabla T |.$$ However, when implementing this equation, I found this effect to be negligible compared to the effects discussed earlier. This means that there is very little heat transfer between squares in the simulation, and the squares are mostly independent. In the real world, I would imagine that there is much more heat transfer than predicted by this model; if not because of Fourier's law, it would exist as a significant factor in climate due to winds and ocean currents.

Heat and Temperature

A general formula relating heat and temperature is $$Q = mC \Delta T$$ where $Q$ is the amount of heat (positive for added or negative for removed), $m$ is the mass of a body, $C$ is the specific heat of the body, and $\Delta T$ is the temperature change experienced by the body. This simulation uses $\Delta T = \frac{Q}{mC}$ to relate our previous work with heat to the temperatures shown in the map above.

Constants

For specific heat, I used (in J/kg K) 1005 for air, 800 for soil, and 4182 for water [4]. For thermal conductivity, I used (in W/m K) 0.020 for air, 0.5 for soil, and 0.606 for water [5]. However, many map squares are a mix of land and water. For these tiles, I estimated the proportion of the tile that was water, and then used a weighted average using this proportion between water and soil for both specific heat and thermal conductivity of the surface. For the atmosphere, I used the values for air.

When finding the thermal conductivity between the surface and atmosphere, I used the value for air, because most of this distance is air. For thermal conductivity between two squares on the ground, I took the average of the thermal conductivities of both squares.

For $dt$ in expressions with $\frac{dQ}{dt}$ (the time increment of this simulation), I used 25 hours. Why 25 and not 24 (one day)? If I had used 24 hours, the same part of the Earth would be in the sun all the time. 25 hours means in addition to a day passing, the sun moves by one hour to light the world evenly.

For the intensity of the sun's radiation $F_S$ I used 1370 W/m² [2].

It can be shown with calculus that if we let the density of the Earth be proportional to $A-Br$, then using the assumption that the relevant crust contains 0.0005% of the Earth's mass, it is around 2 km thick.

References

Area With Polar Coordinates from Paul's Online Math Notes (Lamar University)

Introduction to Atmospheric Chemistry by Daniel J. Jacob (Princeton University Press) (Chapter 7: "The Greenhouse Effect")

Integral Calculator by David Scherfgen

Specific Heat of some common Subtances from Engineering ToolBox

Thermal Conductivity of some selected Materials and Gases from Engineering ToolBox

Surface Area from Paul's Online Math Notes (Lamar University)

Resources

R. Stöckli, E. Vermote, N. Saleous, R. Simmon and D. Herring (2005). The Blue Marble Next Generation - A true color earth dataset including seasonal dynamics from MODIS. Published by the NASA Earth Observatory. Corresponding author: rstockli@climate.gsfc.nasa.gov