Predict Your Local Climate with the Dot Product and ODEs
That’s right — with just the dot product, ordinary differential equations, and some physics, you can probably estimate the climate where you are pretty accurately. If you’ve ever been curious about the math behind why it’s warmer near the equator, or why summer days are longer, that’s contained in our dot product today.
We’ll follow a beam of light from the sun across space, and use its journey to set up differential equations to model the temperature of a given location on earth. We’ll then code up a simulation in Python to see it in action. In the end, we’ll take a closer look at the greenhouse effect and how the math we do today is changing to favor a warmer climate.
By the end of this article, hopefully you’ll be able to accurately predict the long-term climate of the area where you live using just a couple differential equations. As an example, here in New England, this is what I got:
Not bad, especially given that the assumptions I made are relatively simple.
Your results may very depending on where you live, but I hope you’ll find a new appreciation for the different factors that determine the diversity of climates on our planet.
You’ll see that while we tend to see climate as a static thing, some of the variables we’re playing with can be manipulated at will. And humanity has started to manipulate them, not just in a Python simulation, but in real life.
Physics background
The sun is very far away from Earth. So far away, in fact, that we can assume the radiation from the sun is constant in both direction and magnitude by the time it gets here. (If you do the trigonometry, the angular variation is about 0.3° at most.)
That greatly simplifies our calculations. So some (roughly) constant amount of energy per unit area hits the Earth every second. We call this value Fₛ: the solar constant.
That’s where most heat on earth comes from. But where does it go from there? Let’s trace the possible paths.
The value of the solar constant turns out to be about 1360 W/m². But not all of this reaches the surface: generally, only about 1000 W/m² does. The rest is absorbed by the atmosphere. If it’s very cloudy, it’s possible that even less makes it to the ground.
Once light reaches the surface, it may get reflected due to the albedo of earth. Albedo represents the fraction of light that is reflected rather than absorbed by the ground. This varies in different places, but a value of around 0.3 is usually reasonable. Generally this value is higher in deserts or icy areas, and lower in forested areas.
So 1360 W/m² comes in from the sun. About 360 W/m² is absorbed by the atmosphere, while the rest goes down to the ground. Of the remaining 1000 W/m², about 30% is reflected back up off the ground, while the other 70% is absorbed by the ground. When heat is reflected back from the ground, it has the potential to be absorbed by the atmosphere.
We will model the ground as an approximate black body: more specifically, we’ll say it absorbs all the incident energy it receives, except for the portion that is reflected due to albedo. This corresponds to the assumption that the ground is perfectly opaque, although it may still reflect some light.
On the other hand, the atmosphere is not modeled as a black body. When it receives energy, only some fraction f is absorbed. The rest simply passes through.
What is this fraction? Given the temperature of the earth, most of its black body radiation is in the infrared spectrum. About 80% of this is absorbed by the atmosphere due to gases in the atmosphere called greenhouse gases. This absorption and the associated heating of earth is what we call the greenhouse effect.
Setting Up ODEs
Now we have almost enough to set up our differential equations. There’s just a couple more things we need. First, a black body (like the earth in our model) radiates energy (per area, per time) proportional to the fourth power of its temperature. The proportionality constant here, σ, is known as the Stefan-Boltzmann constant.
For the atmosphere, we modify this. First, the atmosphere effectively has twice the surface area of the earth, because it can radiate both downwards to earth and upwards to space. Second, because the atmosphere is not a black body, it emits only a fraction f of what a black body at the same temperature would emit.
We’ve done everything in terms of energy flux, but we’re more interested in changes in temperature.
The change in temperature ΔT depends on the net heat energy input Q according to Q = mCΔT, where m is mass and C is a constant called the specific heat of the substance. This represents how difficult it is to change the temperature of a material. Water has a particularly high specific heat.
Mass is equal to density ρ times surface area A times height h. Using this, we can get φ = ρhC dT/dt where φ is the net energy flux coming in. (This comes from φ = (dQ/dt)/A by definition.)
Determining values for ρhC that lead to accurate results can be tricky, and depends on local geography. We’ll get to that at the end in our simulation.
Let’s quickly recap the path of a beam of light from the sun.
- Light comes in from the sun with an intensity of Fₛ. A fraction α is absorbed by the atmosphere, while the other 1000 W/m2 makes it down to the surface. Let F=1000 W/m².
- Given albedo A, an amount FA of the light that reaches the surface is reflected back up, while F(1-A) is absorbed by the earth.
- Of the light that was reflected in the previous step, an amount αFA is absorbed by the atmosphere, while the other (1-α)FA escapes to space.
On top of this, let’s look at the effects of radiation.
- The earth is radiating energy with intensity σT_E⁴ all the time, where T_E is the temperature of earth. An amount fσT_E⁴ is absorbed by the atmosphere, where f is the fraction of infrared radiation that the atmosphere absorbs due to greenhouse gases.
- The atmosphere is radiating energy with an intensity 2fσT_A⁴ all the time, where T_A is the temperature of the atmosphere. The 2 is for two surfaces, one up to space and one down to earth. Half of this (fσT_A⁴) goes down, and we’ll assume that this is all absorbed by the earth.
That was a lot, but we’re finally ready to write our differential equations. Remember that the intensity φ is equal to ρhC dT/dt where ρ is density, h is height, and C is specific heat.
Letting k=1/ρhC, we get dT/dt=kφ. We need two separate constants k_E and k_A for the earth and atmosphere.
From this, the path of light from the sun, and the effects of radiation, we get our final differential equations:
The Dot Product
All this is great, but it has one glaring flaw: there is no such thing as night. Our model has the sun shining happily all day long and all night long, directly overhead, until the earth burns up. If this were the case, temperatures would look something like this:
Luckily, that’s not the case.
Another issue with our model is that we have no impact of latitude or the seasons, arguably the most important factors in regional and temporal climate variation on earth.
Surprisingly, we can model all three of these factors — night, latitude, and seasons — at the same time. We can do this with the humble dot product.
Imagine it’s currently noon along the prime meridian, and we’re at a longitude φ east/west and an angle θ south of the North Pole. (So a latitude of 45° south corresponds to a θ of 135°.) Using spherical coordinates, we can write our location on earth as a vector:
On the other hand, if it’s currently the September or March equinox, we can write the sun’s rays as a vector entirely in the x-direction as shown in this diagram:
If it’s the June solstice, the northern hemisphere tilts 23.5° towards the sun, effectively raising the sun by an angle 23.5° relative to our diagram. We have this vice versa for the December solstice.
The angle by which we “tilt the sun” (from our perspective) is known as the declination. It can be calculated as approximately
where t is the number of days since the December solstice.
Once we have the declination, we can represent the line of the sun’s rays as a vector with components
Now, remember those constants F and Fₛ from before? Those were for a direct plane wave hitting the earth, which we can think of as the equator at noon on the equinox. For all the other cases, we have to multiply F and Fₛ by the dot product r·s.
What if that dot product is negative? This represents a point on the earth that’s opposite the sun, so it’s night there. We just multiply by 0 in that case.
One more thing: the variable φ isn’t actually longitude. It’s actually longitude relative to the place where it’s currently noon, so we’ll want to make φ itself increase linearly by 2π every 24 hours.
And there you have it: night, seasons, and latitude, all represented with a single dot product!
Simulation
The setup is mostly done now: we just need to input the relevant variables into a differential equation solver and get our result.
To recap, these variables are α for initial atmospheric absorption, A for albedo, f for greenhouse gas absorption of infrared radiation, latitude, and what we called “k” separately for both the earth and the atmosphere, depending on density and specific heat.
For albedo, we looked at an albedo map earlier that you can use for your estimate. f is 0.8. α is (1360–1000)/1360 based on how much solar energy is absorbed by the atmosphere.
You can get latitude from a compass app on your phone, or by looking up your coordinates online.
The constant k is the most tricky out of all these. Remember that in our model k=1/ρhC where ρ is density, h is height, and C is specific heat. For the atmosphere, none of these variables are well-defined. But I found ρ=1005 kg/m³ (density of air), h=10km, and C=1.2 J/kg°C (specific heat of air) to work well.
Our model shouldn’t be taken as literally predicting the atmospheric temperature — there’s far too many factors we haven’t considered. The role of the atmosphere is just another layer to more accurately predict ground temperature.
Things get even more difficult when you do this for the earth, and you might have to play around a bit. To start, you can try ρ=1500 kg/m³, h=20 m, and C=1500 J/kg°C. (h=20 m represents the assumption that only the top 20m or so of the ground is really affected by the incoming solar energy or interactions with the atmosphere.)
In general, you can adjust C based on the impact of water where you live. Dryer and more inland areas will have lower values, maybe down to around 1100. More humid and coastal areas will run higher, perhaps up to around 3500 at the very high end. (To any fish reading this: if you actually live underwater, the value is 4184.)
For Southern New England, near the coast, I used ρ=1500 kg/m³, h=20 m, C=1700 J/kg°C, and A=0.25.
Once you have all your constants, you can simply plug them into a solver. Since we’ve already gone over all the concepts, I’ll just give you the full Python code here. (You’ll need numpy, scipy, and matplotlib installed.)
import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt
lat = 42 # latitude in degrees. Use negative for the southern hemisphere
F0 = 1360 # W/m^2
F = 1000 # W/m^2.
A = 0.25 # albedo
f = 0.8 # greenhouse gas factor
alpha = (F0-F)/F0
s = 5.67e-8 # Stefan-Boltzmann constant. W/(m^2 K^4)
k_E = 1/(1500 * 1700 * 20) # kg/m^3, J/(kg K), m
k_A = 1/(1005*1.2 * 10*1000) # J/(kg K), kg/m^3, m
day = 86400 # seconds
year = day*365.2422 # seconds
theta = np.pi/2-(lat*np.pi/180) # radians
def solvr(T, t):
phi = 2*np.pi/day*t # relative longitude in radians
dec = -np.cos(2*np.pi/year*t)*(23.5*np.pi/180) # declination in radians
F_factor = np.sin(theta)*np.cos(phi)*np.cos(dec) + np.sin(dec)*np.cos(theta) # dot product
Ft = np.maximum(0, F*F_factor)
Fa = np.maximum(0, F0*F_factor)
T_E = T[0]
T_A = T[1]
return [k_E*(Ft*(1-A) - s*T_E**4 + f*s*T_A**4), # dT_E/dt
k_A*(Fa*alpha + f*s*T_E**4 + alpha*Ft*A - 2*f*s*T_A**4) #dT_A/dt
]
timesteps = np.arange(0, year*7, day/24) # times to plot in seconds
asol = integrate.odeint(solvr,
[280, 240], # initial T_E, T_A in Kelvin
timesteps, hmax=day/24)
T_E = (asol[:,0]-273.15)*(9/5)+32 # convert from Kelvin to Fahrenheit
T_A = (asol[:,1]-273.15)*(9/5)+32
plt.plot(timesteps/day, T_E)
plt.xlabel("Time (days)")
plt.ylabel("Temperature (°F)")
plt.show()Try playing around with different values for the constants, both where you are and in other places you want to try out. You might find Western Europe doesn’t align well with our model.
While all our math must be done in Kelvin, my code gives a final output in degrees Fahrenheit. If you don’t love freedom units as much as I do, feel free to remove the final *(9/5)+32 to use Celsius.
The Future of Climate
That’s it for the math today. But before I go, I wanted to make a quick note of the variables in this project and why they’re important not only mathematically, but also in physics, society, and even legislation.
By far the most discussed, and most important, is the variable we called f. This represents the effect of greenhouse gases in the atmosphere. f is a very important variable for our survival: try setting f to zero.
But as I’m sure you’ve heard, humans are pumping greenhouse gases in the atmosphere. The effect on f is small in absolute terms, only enough to change temperatures by a few degrees. But as we’re already seeing, this is enough to lead to natural disasters, rising sea levels, wildfires, and more.
As these issues appear in the news and in conversation, I hope this mathematical model gives you a better understanding of what greenhouse gases do.
Even besides the climate change issue, I find the study of climate to be interesting for its own sake. The diversity of climates across the world is incredible, and this model captures just a few of the most important factors. Despite that, I’m amazed by how this is enough to build a broadly accurate model.
I’d like to acknowledge the book Introduction to Atmospheric Chemistry by Daniel J. Jacob for many of the concepts behind this model. If you want a live online demo of this model, with some sacrifice in accuracy, you can check out this link on my website. Thank you for reading!
