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

Euler's Method¶

Overview¶

Here we'll develop some skill at simulating how the atmosphere and land arrive at a steady state over time. The key piece of this is an algorithm called Euler's method.

Fluxes between the atmosphere and the land using 1st-order kinetics¶

We've already talked about expressing the rate (flux) at which carbon is transferred from the atmosphere to land as a process that depends on how much carbon is in the atmosphere. Ignoring for the moment the temperature feedback part of this, we said

$$ F_{atm->land} = k_{al0} + k_{al1} [C_{atm}] \ \ \ (1) $$

where $k_{al0}$ and $k_{al1}$ are constants. Now, say we also have a flux from the land back to the atmosphere,

$$ F_{land->atm} = k_{la} \ \ \ (2) $$

where $k_{la}$ is another constant.

If these fluxes were in balance, namely, $F_{atm->land}=F_{land->atm}$, we would say that the climate is in a state of mass balance. We already know the word for this: homeostasis. And if they are not in balance? The system will respond:

$$ \Delta [C_{atm}] = (F_{land->atm} - F_{atm->land}) \times \Delta t \ \ \ (3) $$

where $\Delta t$ is some period of time we're interested in, and $\Delta [C_{atm}]$ is how much $[C_{atm}]$ changed during that time.

Euler's method¶

Equations 1-3 constitute what's called a set of ordinary differential equations. In general, there are a lot of methods for solving them. Here, we'll learn to use a fairly simple, numerical one called Euler's Method, because it's pretty straightforward to implement in Python. It says, basically, that at each time step, we'll calculate the fluxes based on concentrations in the previous time step. As long as the steps are tiny, it's a pretty accurate method.

For loops and the "+=" construct¶

As you may well have guessed, those tiny steps will be set up using a for loop, which you already know about. The "updating" of $[C_{atm}]$ just referred to, in each iteration of your for loop, will be accomplished by Python's "+=" construct, which you'll catch on to after an example or two.

Learning objectives¶

  1. I can numerically solve 1st order kinetic problems using Euler's method set up in a for loop.
  2. I can use the "+=" construct.
  3. I can interpret the approach to homeostasis in terms of $t_{1/2}$, and predict $t_{1/2}$ from a rate constant $k$ using the "70s rule."
In [5]:
# Bring in resources
import numpy as np # linear algebra
import matplotlib.pyplot as plt
In [6]:
%matplotlib notebook

Implementing Euler's method¶

Now we're going solve for $[C_{atm}]$ as a function of time, according to the theory laid out in the Introduction.

In [7]:
# Constants
k_al0 = 113
k_al1 = 0.0114
k_la = 120

# This is the starting amount of carbon in the atmosphere, in GtC
C_atm = 500

# Here's a time array
nsteps = 500
time = np.linspace(0,250,nsteps)

# Calculating the time step
dt = time[1]-time[0]

# Empty array that will hold GtC carbon in the air over time
C_atm_array = np.empty(0)

# Looping
for i in range(len(time)):
    
    # This calculates the fluxes of this iteration of the loop
    F_al = k_al0 + k_al1*C_atm
    F_la = k_la
    
    # This calculates the change in C_atm
    delta_C_atm = (F_la - F_al)*dt
    
    # Get the new C_atm (and append it to our array)
    C_atm += delta_C_atm
    C_atm_array = np.append(C_atm_array,C_atm)

# Plot the results
plt.figure()
plt.plot(time,C_atm_array)
plt.grid(True)
plt.xlabel('time (years)')
plt.ylabel("[C_atm] (GtC)")
Out[7]:
Text(0, 0.5, '[C_atm] (GtC)')

Pause for analysis¶

You can probably see that this earth has not quite yet achieved homeostasis, but it's getting there! We're going to do a little off-line analysis here, a key part of which is to recognize that the only part of our model that depends on time is the parameter $k_{al1}$, which plays the role of a rate constant. Because of that, the dynamics we're looking at are subject to the "rule of 70s." To explore that, do the following:

  1. Just guestimating, about how many years do you think it takes for $C_{atm}$ to get half-way to homeostasis?
  2. Use the "rule of 70s" to predict a value for $t_{1/2}$. Does your value agree (more or less) with the estimate you just got?

BEGIN SOLUTION¶

The time constant is about 61 years, which is agreement with what the graph says.

END SOLUTION¶

Your turn¶

Now you're going do almost the exact same thing, but with a few wrinkles. Here's what to do:

  1. Extend the time range to go out to 10 times $t_{1/2}$ years;

  2. Inside your loop, append the concentration of carbon in the atmosphere to your array as before, but also append the fluxes to their own arrays; and

  3. Once the loop is over, make three plots:

    -In one window, plot the concentration as a function of time, as before;
    -In a second window, plot the two flux arrays together. Annotate these using the label/legend method; and
    -In a third window, plot the net flux to the atmosphere, defined as $F_{net}=F_{la}-F_{al}$.

In [8]:
# We'll use the same starting amount of carbon in the atmosphere, in GtC
C_atm = 500

# Initializing three arrays (for atmospheric carbon, flux in, and flux out)
C_atm_array = np.empty(0)
F_la_array = np.empty(0)
F_al_array = np.empty(0)

# Now, the loop ...
#    Start with creating the time array using np.linspace, this time going up to 3 x t_1/2
#    You'll need to calculate the time step of your new array by subtracting the 2nd time from the 1st
#    Loop over each element in the time array, calculate the fluxes and new concentrations, and append to arrays

### BEGIN SOLUTION
thalf = 70/k_al1/100
print('Time constant = ', thalf)
tstop = 10*thalf
time = np.linspace(0,tstop,nsteps)
dt = time[1]-time[0]
for i in range(len(time)):
    
    # This calculates the fluxes of this iteration of the loop
    F_al = k_al0 + k_al1*C_atm
    F_la = k_la
    
    # This calculates the change in C_atm
    delta_C_atm = (F_la - F_al)*dt
    
    # Get the new C_atm (and append it to our array)
    C_atm += delta_C_atm
    C_atm_array = np.append(C_atm_array,C_atm)
    
    # Same with the new fluxes
    F_la_array = np.append(F_la_array,F_la)
    F_al_array = np.append(F_al_array,F_al)
### END SOLUTION

# Now plot all the concentration results over time
plt.figure()
### BEGIN SOLUTION
plt.plot(time,C_atm_array)
plt.grid(True)
plt.xlabel('time (years)')
plt.ylabel("[C_atm] (GtC)")
### END SOLUTION

# The two flux results (in and out) on the same graph, labeling using the "label/legend"
plt.figure()
### BEGIN SOLUTION
plt.plot(time,F_la_array,label='F_la')
plt.plot(time,F_al_array,label='F_al')
plt.grid(True)
plt.xlabel('time (years)')
plt.ylabel("Fluxes (GtC/year)")
plt.legend()
### END SOLUTION

# The net flux to the atmosphere
plt.figure()
### BEGIN SOLUTION
plt.plot(time,F_la_array-F_al_array,label='Net flux to atmosphere')
plt.grid(True)
plt.xlabel('time (years)')
plt.ylabel("Fluxes (GtC/year)")
plt.legend()
### END SOLUTION

Time constant =  61.40350877192982
Out[8]:
<matplotlib.legend.Legend at 0x7efffa270970>

Pause for analysis¶

Looking at the beginning of your flux graphs, does the flux going into the atmosphere ($F_{la}$) exceed the flux from the land to the atmosphere ($F_{al}$), or the other way round? Hopefully, your answer is consistent with the fact that the concentration of carbon in the atmosphere is rising as it approaches homeostasis.

BEGIN SOLUTION¶

END SOLUTION¶

Refresh/save/validate¶

Double-check everything is OK, and press the "Validate" button (as usual).

Close/submit/logout¶

Close, submit, and log out.