Earth Temperature Integrator (Interactive)¶

#Import libraries
%matplotlib widget
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider, Button
import matplotlib.animation as animation
from matplotlib.patches import Rectangle

Integrator¶

Before we discuss the integrator, I’ll go over the adjustments to the equation we make.

Addition of atmosphere¶

To account for the atmosphere, we simply add on another flux term to the right side. This time, it’s possitive as it represents flux being sent back to the Earth from the atmosphere.

We use the term: emissivity * sigma * Temp_Atmosphere^4

For simplicity, we use an estimate for Temp_Atmosphere = 242.4 K.

Removal of emissivity from outgoing flux¶

In this model, we use the emissivity as a property of the atmosphere, so we do not multiply the outgoing flux by it.

Integration over heat content¶

We actually wont integrate directly over the temeprature, but rather over the heat content, descirbed by:

C * dT/dt = dH/dt, where dH/dt is the heat content.

We can then find the temperature at each step by dividing the heat content by C.

Interactivity¶

The code allows the user to input the emissivity, albedo, and solar flux. These initial conditions will change how the planet’s temperature evolves. An animation shows the results and whether the planet will approach a habitable temperature. The integrator starts with a surface temperature T = 0 K.

#Set up integrator

#constants
C = 2.08e8
alpha = float(input("Enter the starting albedo (between 0 and 1.0): ")) #Best = 0.32
epsilon = float(input("Enter the starting emissivity (between 0 and 1.0): ")) #Best = 0.78
sigma = 5.67e-8
S = int(input("Enter the starting solar flux: ")) #Best = 1350

#Initialize
steps = 450
T_now = (np.zeros(steps))
H_now = (np.zeros(steps) +1)
T_next = np.zeros(steps) 
H_next = np.zeros(steps)


dt = 1e6 #seconds

#Differential equation to integrate
#dH = ((1 - alpha)*S)/4 - epsilon*sigma*(T**4)

for i in range(0, steps - 1):
    H_next[i] = H_now[i] + (((1 - alpha)*S)/4 - sigma*(T_now[i]**4) + epsilon*sigma*(242.4**4)) * dt
    
    H_now[i+1] = H_next[i]
    
    T_next[i] = H_now[i] / C                       
    
    T_now[i+1] = T_next[i]
                            

Enter the starting albedo (between 0 and 1.0):  0.9
Enter the starting emissivity (between 0 and 1.0):  0.1
Enter the starting solar flux:  2200

#Set up animation

#initialize arrays
step = []
temp = []

fig = plt.figure()
ax1 = fig.add_subplot(1,1,1)
plt.title('Earth Temperature Evolution')

ax1.set_ylim(0,400)
ax1.set_xlim(0,500)
ax1.set_xlabel('Step')
ax1.set_ylabel('Temperature (K)')
ax1.add_patch((Rectangle((0,280), 500, 20, color = 'green', alpha = 0.2))) #add habitable patch

def animate(i):
    step.append(i)
    temp.append(T_now[i])
    ax1.plot(step, temp, color = 'red')


ani = animation.FuncAnimation(fig, animate, interval = 40)
plt.show()

A Simple Energy Balance Model