Computational Guided Inquiry for Modeling Earth's Climate (Neshyba, 2023)

Simple Radiative Balance

Absorbed Solar Radiation

The energy responsible for keeping Earth warm comes almost entirely from the Sun. This energy is all "light", but humans can see only part of it; of the rest, some comes in the form of ultraviolet light (which we can detect because it burns our skin) and near-infrared light (which is what you are feeling when you hold your hands out in front of a fire). All three forms of light (uv, visible, and near-IR) are collectively called shortwave radiation.

How intense is shortwave radiation from the sun? To answer that, we need to talk about an instrument for measuring shortwave flux, called a shortwave radiometer. Shortwave radiometers measure the amount of shortwave energy (remember, that's uv, visible, and near-IR) passing through a 1 meter x 1 meter hoop. Real shortwave radiometers aren't actually that big, but they still tell you that key number: the energy intensity of sunlight on a per-meter-squared basis.)

OK, now imagine you were are floating in space, just outside Earth's atmosphere, armed with your shortwave radiometer. It should read a number around $1367 \ W m^{-2}$. Is that a lot? Well, $1367 \ W m^{-2}$ is almost fourteen 100-Watt light bulbs, for every 1x1 square meter! This number is called the solar constant, and symbolized by $S_0$, so we say $S_0=1367 \ W m^{-2}$. $S_0$ varies a little over the course of a year (and of course is smaller for planets farther from the sun, bigger for planets that are closer), but $1367$ is a good average for us.

As it turns out, however, that $1367 \ W m^{-2}$ doesn't actually make it down to every square meter of Earth's surface below you. There are two reasons for this:

1. The sun shines only on one side of the planet at a time, and it usually hits Earth's surface at an oblique angle. Turns out, this reduces the average energy potentially received by earth by a factor of 4, i.e., we're talking $S_0/4=342 \ Wm^{-2}$.
2. Clouds, snowpacks, desert sand, even forests, tend to reflect the shortwave light away, back out to space. The overall fraction reflected is called Earth's albedo, and it's about 30%. The albedo is typically given the symbol $\alpha$, and it's given as a fraction rather than a percent. So the statement "Earth's albedo is 30%" is written mathematically as $\alpha=0.3$.

If we put this all together, we can write down something called the Absorbed Shortwave Radiation, which we'll symbolize here as $ASR$. It can be expressed as

$$ ASR = \dfrac{S_0}{4}(1-\alpha) \ \ \ \ (1) $$

Like $S_0$, the units of $ASR$ are $W m^{-2}$. If you plug in the above numbers for $S_0$ and $\alpha$, you get $ASR=239 \ W m^{-2}$. That's still a lot if you think about it in terms of light bulbs: it's like there are over two 100-W lightbulbs of sunlight being absorbed by Earth's surface, day and night, over every square meter of the planet's surface.

Outgoing Longwave Radiation

Imagine now that you are standing at the earth's surface, armed with a longwave radiometer. That's an instrument that measures the longwave radiant energy passing through a 1 meter x 1 meter hoop. To a good approximation, we'd expect that number to be given by the expression $\sigma T^4$, where $T$ is the temperature of the ground you're standing on. $\sigma$ is known as the Stefan-Boltzmann constant, $5.67 \times 10^{-8} W m^{-2}K^{-4}$. That's just physics, it doesn't depend on climate.

For the purpose of radiative balance calculations, however, we need to re-position ourselves back to where we were at the beginning of this story: hovering in space just above Earth's atmosphere. Looking down at the earth's surface, the quantity our longwave radiometer measures now is called the Outgoing Longwave Radiation, or OLR. It will be, it turns out, not quite equal to the $\sigma T^4$ coming off Earth's surface. Why? Well first, something like 85% of the longwave radiation emitted by Earth's surface gets captured by water vapor, $CO_2$ and $CH_4$ in the atmosphere (they're all greenhouse gases). So we'd expect OLR to be reduced to something like 15% of $\sigma T^4$, right? However, it turns out that the atmosphere itself emits longwave radiation. Combining these influences, we can write

$$ OLR = \kappa \sigma T^4 \ \ \ \ (2) $$

where (as it turns out) $\kappa=0.614$.

(Here's a more technical version of what we just said: Atmospheric greenhouse gases absorb 85% of the longwave photons that originate from Earth's surface, but those same gases emit their own longwave photons. Those longwave photons that are emitted by the atmosphere back down toward the earth's surface are responsible for greenhouse warming. Those longwave photons that are emitted by the atmosphere out to space combine with the 15% of Earth's surface photons that manage to get through the atmosphere unscathed, adding up to 64% of the orignal Stefan-Boltzmann ($\sigma T^4$) energy intensity coming off the earth's surface. That leads to the expression for the Outgoing Longwave Radiation in Eq. (2), with $\kappa=0.614$.)

Radiative Balance

All this leads to an algorithm for modeling climate, based on radiative balance. It goes like this: suppose the earth were initially really cold. Then $T$ in Eq. (2) would be small, and that would mean the outgoing longwave radiation, $OLR$, would also be small -- smaller than our ASR of $239 \ W m^{-2}$ coming in. That means the earth's surface would be experiencing a positive flux imbalance -- more energy coming in than going out. What will the earth do then? Well, heat up of course!

Things don't heat up instantly, however. If you're familiar with trying to boil water on a stovetop, you know that the more water in your pot, the longer it takes to heat up. We say that two gallons of water have twice the heat capacity as one gallon. The relevance of heat capacity to our climate model is, even when there's a positive flux imbalance, it would still take time to heat up. We can even write an equation for this:

$$ \Delta T = \dfrac{ASR-OLR}{C_p}\times \Delta t \ \ \ \ (3) $$

where $\Delta T$ is the change in temperature of the earth's surface in the time $\Delta t$, and $C_P$ is our symbol for the heat capacity -- in our model, the heat capacity of a square meter of the average Earth's surface. It turns out that this second number is known: $C_P\approx 2.2\times 10^8 \ J \ m^{-2}K^{-1}$. And notice that it's in the denominator -- which means larger heat capacity means smaller change in temperature. Another term in this equation worth noticing is $ASR-OLR$. $ASR-OLR$ is likely to be positive when the earth is cold (because Eq. (2) says $OLR$ will be small), which means there would be a positive flux difference, which would mean $\Delta T>0$ (i.e. the temperature of the earth will rise).

You can imagine the opposite process when Earth's temperature is initially very high.

We haven't mentioned the factor $\Delta t$ yet. The lower case means we're talking about time, not temperature. And hopefully its appearance in Eq. (3) is fairly intuitive: the longer you heat something (bigger $\Delta t$), the more it will heat up (bigger $\Delta T$).

Equation (3) is the key equation in our model. It allows us to start with a given temperature, $T$, then march through time, at each time step either raising the temperature or lowering it, according to whether $ASR-OLR$ is positive or negative. In "pseudo-code", we could say we have a loop:

while (Delta_T is big):
    calculate OLR from the current temperature using Eq. (2)
    calculate Delta_T from Eq. (3)
    update temperature using T(new) = T(old) + Delta_T

where Delta_T is what we've written here, mathematically, as $\Delta T$.

Why is the condition that we keep going as long as Delta_T is big? Well, there's no point in continuing on if there's no change in temperature, right? In fact, what you'll see is that the temperature calculated in this way evolves toward a steady-state value. If Earth is initially too cold, it gets warmer; if it's too hot, it cools down. That means it's not just any old steady state, it's a stable steady state.

Climate sensitivity to changes in albedo

It turns out that in modeling the climate, a key thing one looks for is what's called the climate sensitivity. For example, we could say

$$ S = \dfrac {\Delta T_{SS}}{\Delta \alpha} \ \ \ \ (4) $$

where $S$ is the sensitivity of Earth's temperature to its albedo ($\alpha$). To find this number from our model, we'd need to run the model once with a "reference" value of the albedo, and then again with a new value of the albedo. We'll call the difference in these two albedos $\Delta \alpha$, and the difference in the resulting steady-state temperatures $\Delta T_{SS}$. From that, we can use Eq. (4) to calculate the climate sensitivity, $S$.

In the code below, you'll do this, and then use the resulting value of $S$ to predict the impact some new albiedo might have on the temperature of the planet -- details to be given below.

Caveats

The model described here ignores a key feature of climate called feedbacks. For example, if lots of snow melts because of higher temperatures (which it will do), Earth's albedo will go down, which will cause Earth's temperature to rise, which will cause more snow to melt: a positive feedback. But that's a challenge for another day.

Learning Goals

1. I'm conversant with fundamental quantities associated with radiative balance, like $\alpha$, $S_o$, $ASR$, $OLR$, $\kappa$, and $\sigma$, and I can supply typical values, with units.
2. I can explain the sense in which the termination of "while" loop in the model is consistent with a stable steady-state temperature.
3. I can explain how to calculate a climate sensitivity (to albedo), and how to use that value as a predictive tool.

import numpy as np
import matplotlib.pyplot as plt

%matplotlib inline

Model Definition

The parameters in this cell don't need to be changed from run to run (unless you want to change the model).

def SRB(temperature_start,alpha):

    # Some constants
    solar_constant=342  # in units W m^-2
    kappa=0.614         # dimensionless
    sigma=5.67e-8       # in units W m^-2 K^-4
    heat_capacity=2.2e8 # in units J m^-2 K^-1
    
    # Time step parameters
    delta_time_years = .1    # Number of years for each time step
    delta_time = delta_time_years*24*60*60*365 # This converts the time step to seconds
    
    # Loop control 
    small_enough = 0.001
    delta_temp = 1000

    # Emergency exit parameter
    maximum_steps = 10000

    # Initializing the arrays that will contain our results
    temperature_list = [temperature_start]
    time_list = [0]
    fluxdiff_list = [0]

    # Looping
    while abs(delta_temp) > small_enough:

        # Let's extract the latest temperature and time
        temperature = temperature_list[-1]
        time = time_list[-1]

        # Outgoing longwave radiation (Wm^{-2})
        OLR=kappa*sigma*temperature**4

        # Absorbed solar radiation (also Wm^{-2})
        ASR=(1-alpha)*solar_constant

        # Flux difference (also Wm^{-2})
        fluxdiff=ASR-OLR

        # What's the temperature adjustment this time period?
        delta_temp = fluxdiff/heat_capacity*delta_time

        # Adjust Earth's temperature for this step
        newtemperature = temperature + delta_temp
        newtime = time + delta_time_years

        # Add to our list of temperatures and times
        temperature_list.append(newtemperature) 
        time_list.append(newtime)
        fluxdiff_list.append(fluxdiff)

        # Emergency exit
        if (len(time_list)>maximum_steps):
            break
            
    return time_list[1:], fluxdiff_list[1:], temperature_list[1:]

Pause for analysis

Spend a little time getting to know this function:

1. What is the physical meaning of the condition that terminates the "while" loop? (Not talking about the emergency exit.)
2. There is one quantity that is calculated in the loop, that doesn't need to be in the loop, since it could easily (and more efficiently) be calculated before the loop begins. Which quantity is that?

Experiment 1

Run the cell below to see how the model responds to an initial temperature of 200 degrees K and an albedo of 0.3.

# Experiment 1

# Run the model
time_list_expt1, fluxdiff_list_expt1, temperature_list_expt1 = SRB(200,0.3)

# Plot the flux differences
plt.figure()
plt.plot(time_list_expt1,fluxdiff_list_expt1,'black')
plt.xlabel('time (years)')
plt.ylabel('Net Flux (W/m^2)')
plt.title('Experiment 1')
plt.grid(True)

# Plot the temperatures
plt.figure()
plt.plot(time_list_expt1,temperature_list_expt1,'black')
plt.xlabel('time (years)')
plt.ylabel('temperature (K)')
plt.title('Experiment 1')
plt.grid(True)

# Report the final (steady-state) temperature
print("The temperature ends up at ", temperature_list_expt1[-1])

The temperature ends up at  287.94734398602714

Pause for analysis

In the cell below, write a few descriptive sentences about this:

1. Is the initial net flux positive or negative?
2. Is the temperature initially rising or falling?
3. What does it mean, physically, that the net flux levels off to a value of zero?

Your turn

Carry out two more numerical experiments:

• In the first cell below, carry out "Experiment 2". This experiment should have a starting temperature of 250 K. In your graph, use a blue line.
• In the following cell, carry out "Experiment 3". This experiment should have a starting temperature of 350 K. In your graph, use a red line.

# Experiment 2

# Run the model with a starting temperature of 250 K
### BEGIN SOLUTION
time_list_expt2, fluxdiff_list_expt2, temperature_list_expt2 = SRB(250,0.3)
### END SOLUTION

# Plot the flux differences
### BEGIN SOLUTION
plt.figure()
plt.plot(time_list_expt2,fluxdiff_list_expt2,'blue')
plt.xlabel('time (years)')
plt.ylabel('Net Flux (W/m^2)')
plt.grid(True)
plt.title('Experiment 2')
### END SOLUTION

# Plot the temperatures
### BEGIN SOLUTION
plt.figure()
plt.plot(time_list_expt2,temperature_list_expt2,'blue')
plt.xlabel('time (years)')
plt.ylabel('temperature (K)')
plt.grid(True)
### END SOLUTION

# Report the final (steady-state) temperature
### BEGIN SOLUTION
print("The temperature ends up at ", temperature_list_expt2[-1])

# Not part of the solution, just some analytical exploration
solar_constant=342  # in units W m^-2
kappa=0.614         # dimensionless
sigma=5.67e-8       # in units W m^-2 K^-4
alpha = 0.3
T_ss_analytical = (solar_constant*(1-alpha)/(kappa*sigma))**(1/4)
print('T_ss_analytical =', T_ss_analytical)
base_term = 1/4*(solar_constant*(1-alpha)/(kappa*sigma))**(1/4-1)
sensitivity_alpha_analytical = base_term*(-solar_constant/(kappa*sigma))
print('sensitivity_alpha_analytical =', sensitivity_alpha_analytical)
sensitivity_solar_constant_analytical = base_term*((1-alpha)/(kappa*sigma))
print('sensitivity_solar_constant_analytical =', sensitivity_solar_constant_analytical)

r_aphelion = 94.5
r_perhelion = 91.4
# The factor of 2 here represents up and down from average distance from the sun
delta_eccentricity = -(1/r_aphelion**2 - 1/r_perhelion**2)*r_perhelion**2/2 # This is a fraction
print('delta_eccentricity =',delta_eccentricity*100)

delta_solar_constant = solar_constant*delta_eccentricity
print('delta_solar_constant =', delta_solar_constant)

delta_T0_Milankovic = sensitivity_solar_constant_analytical*delta_solar_constant
print('delta_T0_Milankovic =',delta_T0_Milankovic)

# 5 here because the range was 10 degrees
delta_T_Pleistocene = 5
gain = delta_T_Pleistocene/delta_T0_Milankovic
print('Pleistocene gain =',gain)
feedback = 1-1/gain
print('Pleistocene feedback parameter =',feedback)

# Sensitivity to CO2 doubling (https://en.wikipedia.org/wiki/Radiative_forcing)
# delta_F_doubling = 3.7/290 # W/m^ for 2x CO2
# delta_F_so_far = delta_F_doubling*(420-290)
# "so far" = 2020
# delta_F_so_far = 5.35 * np.log(423/290)
delta_F_so_far = 2.5
print('delta_F_so_far =', delta_F_so_far, 'W/m^2')
# See https://www.climate.gov/news-features/understanding-climate/climate-change-global-temperature
delta_T_so_far_observed = 1.1
print('delta_T_so_far_observed =', delta_T_so_far_observed)
observed_sensitivity = delta_T_so_far_observed/delta_F_so_far
print('observed_sensitivity =',observed_sensitivity)
eventual_temperature_rise = sensitivity_solar_constant_analytical*delta_F_so_far*gain
print('eventual_temperature_rise =', eventual_temperature_rise)
### END SOLUTION

The temperature ends up at  287.94821469283823
T_ss_analytical = 287.9673030294036
sensitivity_alpha_analytical = -102.84546536764415
sensitivity_solar_constant_analytical = 0.21050241449517806
delta_eccentricity = 3.2266173959295568
delta_solar_constant = 11.035031494079085
delta_T0_Milankovic = 2.3229007735339797
Pleistocene gain = 2.152481094744816
Pleistocene feedback parameter = 0.5354198452932041
delta_F_so_far = 2.5 W/m^2
delta_T_so_far_observed = 1.1
observed_sensitivity = 0.44000000000000006
eventual_temperature_rise = 1.1327561689975196

# Experiment 3

# Run the model with a starting temperature of 350 K
### BEGIN SOLUTION
time_list_expt3, fluxdiff_list_expt3, temperature_list_expt3 = SRB(350,0.3)
### END SOLUTION

# Plot the flux differences
### BEGIN SOLUTION
plt.figure()
plt.plot(time_list_expt3,fluxdiff_list_expt3,'red')
plt.xlabel('time (years)')
plt.ylabel('Net Flux (W/m^2)')
plt.grid(True)
plt.title('Experiment 2')
### END SOLUTION

# Plot the temperatures
### BEGIN SOLUTION
plt.figure()
plt.plot(time_list_expt3,temperature_list_expt3,'red')
plt.xlabel('time (years)')
plt.ylabel('temperature (K)')
plt.grid(True)
### END SOLUTION

# Report the final (steady-state) temperature
### BEGIN SOLUTION
print("The temperature ends up at ", temperature_list_expt3[-1])
### END SOLUTION

The temperature ends up at  287.9865659526379

Your turn (again)

Assuming you've given your results unique names, plot all three temperatures on the same graph, as a function of time.

### BEGIN SOLUTION
# Plot the temperatures
plt.figure()
plt.plot(time_list_expt1,temperature_list_expt1,'black')
plt.plot(time_list_expt2,temperature_list_expt2,'blue')
plt.plot(time_list_expt3,temperature_list_expt3,'red')
plt.xlabel('time (years)')
plt.ylabel('temperature (K)')
plt.grid(True)
plt.xlim([0,18])

### END SOLUTION

(0.0, 18.0)

Your challenge - Climate sensitivity

Here, we'd like to see how the steady-state temperature depends on the albedo. In order to calculate it, we could go about it in the following way. First, specify a new value of the albedo. For example, you could set up Experiment #4, with $\alpha=0.4$. Then, since Experiment #3 used $\alpha=0.3$, we could say

$$ \Delta \alpha = 0.1 $$

When you run Experiment #4, you'll get a different steady-state temperature, of course. If we said $T_4$ is that temperature, and $T_3$ was the steady-state temperature you got in experiment #3, you could find the change in temperature resulting from the change in albedo by calculating

$$ \Delta T_{SS} = T_4-T_3 $$

Of course, the bigger the change in the albedo, the bigger we expect $\Delta T$ to be. What's key, therefore, is the ratio, $\dfrac {\Delta T_{SS}}{\Delta \alpha}$. As described in the introduction, this is called the climate sensitivity to changes in albedo. Doing this in Python, it will be helpful to create a variable holding this value, with a line like the following:

sensitivity = delta_T_SS/delta_alpha

Note that the sign of the sensitivity is important here: if you get a negative number, you're probably on track.

### BEGIN SOLUTION
# Run the models for two different albedo values
time_list_expt3, fluxdiff_list_expt3, temperature_list_expt3 = SRB(200,0.3)
time_list_expt4, fluxdiff_list_expt4, temperature_list_expt4 = SRB(200,0.4)

plt.figure()
plt.plot(time_list_expt3,temperature_list_expt3,'red')
plt.plot(time_list_expt4,temperature_list_expt4,'blue')
plt.xlabel('time (years)')
plt.grid(True)
plt.ylabel('temperature (K)')
plt.legend(['alpha=0.3','alpha=0.4'])

# Calculate differences
delta_T_SS = temperature_list_expt4[-1]-temperature_list_expt3[-1]
delta_alpha = 0.1

# Calculate the climate sensitivity
sensitivity = delta_T_SS/delta_alpha
print("The temperature sensitivity (to changes in albedo) is", sensitivity)
### END SOLUTION

The temperature sensitivity (to changes in albedo) is -108.88877178022824

Pause for analysis

But will Earth's albedo change in a warming climate? That's a big question, which we won't be able to do justice to here. What we can do, is use the foregoing Simple Radiative Balance ("SRB") theory to predict the temperature impact of a hypothetical change. Here's how: if you solve Eq. (4) for $\Delta T_{SS}$, you get

$$ \Delta T_{SS} = S \ \Delta \alpha $$

So what's a reasonable guess for $\Delta \alpha$? A recent study by Goode et al, described in a press release at https://news.agu.org/press-release/Earth-is-dimming-due-to-climate-change, reported that over the past two decades, Earth's average albedo dropped from $0.30$ to $0.295$ -- a change in albedo of $-0.005$. In that case, the corresponding Python code could look like

deltaT_SS = sensitivity*(-0.005)

In the cell below, use SRB theory to predict the change in temperature resulting from this reduction in albedo.

### BEGIN SOLUTION
impact = sensitivity*(-.005)
print(impact)
### END SOLUTION

0.5444438589011412

Pause for analysis

1. What phenomenon are Goode et al actually measuring, that led to the reduction in albedo they report? It will help to take a look at the press release and the paper by Goode et al to answer this. Hint: the phenomenon has been noticed for a long time -- including by Leonardo da Vinci!
2. By comparison with measurements by NASA’s Clouds and the Earth’s Radiant Energy System, Goode et al speculate a reason behind the albedo reduction. What is that reason?

Refresh/save/validate

Almost done! To double-check everything is OK, repeat the "Three steps for refreshing and saving your code," and press the "Validate" button (as usual).

Close/submit/logout

1. Close this notebook using the "File/Close and Halt" dropdown menu
2. Using the Assignments tab, submit this notebook
3. Press the Logout tab of the Home Page