James Lovelock
Part of the charm in solving a differential equation is in the feeling that we are getting something for nothing. So little information appears to go into the solution that there is a sense of surprise over the extensive results that are derived.
George Robert Stibitz
As James Lovelock points out when it comes to climate change timescale is everything. I am involved in a study of the Pleistocene epoch, a period from about 2.6 million years ago to about 12 thousand years ago. This was a period of cycles of advancing and retreating glaciers and ice sheets, especially in the Northern Hemisphere. Obviously, there was extensive change in the global climate during this period. However, the geological record shows that the climate changes occurred over periods of thousands to tens of thousands of years rather than hundreds of years.
I'm not a climate scientist. I'm more of a data science hacker who is now involved in looking at paleoclimate data.. In order to get some background, I have been reading reading a few papers about climate models and the proxies climate scientists use to estimate quantities like $C{O_2}$ levels.
I recently read a paper entitled Modeling the Dynamics of Glacial Cycles by Hans Engler, Hans G. Kaper, Tasso J. Kaper, and Theodore Vo and a companion paper Dynamical systems analysis of the Maasch–Saltzman model for glacial cycles by the same authors. In the papers they present a simple climate model based on a model originally presented by Maasch and Saltzman [J. Geophys. Res., 95, D2 (1990), pp. 1955-1963]. Since, I sometimes feel that I don't understand a subject unless I can write a program that implements it, I decided to try and reproduce one of the results from the paper.
From Modeling the Dynamics of Glacial Cycles paper:
To reconstruct the Pleistocene climate, geoscientists rely on geological proxies, particularly a dimensionless quantity denoted by $\delta {{\rm O}^{18}}$, which is measured in parts per mille. This quantity measures the deviation of the ratio $\frac{{\delta {{\rm O}^{18}}}}{{\delta {{\rm O}^{16}}}}$ of the stable oxygen isotopes $\delta {{\rm O}^{18}}$ and $\delta {{\rm O}^{16}}$ in a given marine sediment sample from the same ratio in a universally accepted standard sample. The relative amount of the isotope $\delta {{\rm O}^{18}}$ in ocean water is known to be higher at tropical latitudes than near the poles, since water with the heavier oxygen isotope is slightly less likely to evaporate and more likely to precipitate rst. Similarly, water with the lighter isotope $\delta {{\rm O}^{16}}$ is more likely to be found in ice sheets and in rain water at high latitudes, since it is favored in atmospheric transport across latitudes. The global distribution of $\delta {{\rm O}^{18}}$ in ocean water therefore varies in a known way between glacial and interglacial periods. A record of these variations is preserved in the calcium carbonate shells of foraminifera, a class of common single cell marine organisms. These fossil records may be sampled from deep sea sediment cores, and their age and $\delta {{\rm O}^{18}}$ may be determined. Details are described in lisiecki2005pliocene.
I won't go into the all of the gory details of the model. The paper is pretty clear about the details of its development. The Maasch and Saltzman model is based on five state variables:
- The total ice mass.
- The concentration of $\delta {{\rm O}^{18}}$ in the atmosphere.
- The volume of the North Atlantic Deep Water (NADW).
- The global mean sea surface temperature.
- The mean volume of summer sea ice
The volume of the NADW is a measure of the strength of the global oceanic circulation. Notice that there are no external inputs to this model. An important climate input is the level of solar activity. The earth's orbit has variations from it's basic elliptical path. These variations affect the amount of solar radiation reaching the earth's atmosphere and thus climate. The Maasch and Saltzman model is a conceptual model based on physical arguments. It describes internal climate dynamics, but is far from the whole story. We can think of it as the starting point for a more complex model.
The Maasch and Saltzman Model
The basic model consists of three coupled ordinary differential equations (ODEs).
${\frac{{dx}}{{dt}} = - (x + y)\\}$
${\frac{{dy}}{{dt}} = (r - {z^2})y - pz + s{z^2}\\}$
${\frac{{dz}}{{dt}} = - q(x + z)\\}$
x, y, z are dimensionless quantities that represent the deviations from the global average of the amount of ice, the atmospheric $C{O_2}$ level, and the volume of the NADW (the oceanic $C{O_2}$ pump). Time is measured in 10K years. p, q, r, and s are positive dimensionless parameters. See the paper for the details of their definition.
My interest is in reproducing figure 3.1 from the paper. To do this, I wrote a simple Python program to solve the dynamic system and plot the results.
"""
climate3.py - A simple climate model based on
https://open.bu.edu/bitstream/handle/2144/28962/1705.07387v1.pdf?sequence=1&isAllowed=y
author: Bill Thompson
license: GPL
copyright: May 28, 2018
"""
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
"""
The Maasch and Saltzman model
This function is called by the ODE solver.
arguments:
t - time
u - a 3 element vector representing x,y,z from the paper
args - a tuple of the p, q, r, s arguments to the model
returns:
f - a list of the current values of u
note:
The form of this function is dictated by the ODE solver, in this case
scipy.integrate.solve_ivp.
"""
def maasch_saltzman(t, u, args):
p, q, r, s = args
f = [-(u[0] + u[1]),
(r - u[2]**2) * u[1] - p * u[2] + s * u[2]**2,
-q * (u[0] + u[2])]
return f
"""
plot_it - use matplotlib to plot the results of the ODE solver
arguments:
sol - the solver object returned by solve_ivp
"""
def plot_it(sol):
fig, ax = plt.subplots()
ax.plot(sol.t, sol.y[0], label="Ice", color="black")
ax.plot(sol.t, sol.y[1], label="CO2", color="red")
ax.plot(sol.t, sol.y[2], label="NADW", color="blue")
ax.legend(loc="upper left")
plt.xlabel("Time [10K years]")
plt.ylabel("x,y,z")
plt.show()
def main():
# parameter values from figure 3.1
p = 1.0
q = 1.2
r = 0.8
s = 0.8
# initial conditions - I took these from the figure 3.1 so may not exactly
# match those used in the paper
u0 = [-0.4, 0.1, 0.9]
abserr = 1.0e-8
relerr = 1.0e-6
stoptime = 40.0
# ODE solver
# use lambda to define a warpper because solve_ivp doesn't pass
# arguments directly
sol = solve_ivp(fun = lambda t, y: maasch_saltzman(t, y, (p, q, r, s)),
t_span = (0, stoptime),
y0 = u0,
rtol=relerr,
atol=abserr)
plot_it(sol)
if __name__ == "__main__":
main()
Here are the results.
This plot matched figure 3.1 from the paper. It shows cycles in the ice mass, level of $C{O_2}$, and volume of NADW. The variables are correlated. As the concentration of atmospheric $C{O_2}$ increases, the climate gets warmer and the total ice mass decreases; as the volume of NADW increases, the strength of the North Atlantic circulation increases, more atmospheric $C{O_2}$ is absorbed by the ocean and, consequently, the atmospheric $C{O_2}$ concentration decreases.
The Same Thing in Julia
I have been experimenting with the Julia language. Here is the same set of equations solved using Julia's DifferentialEquations module.
#=
climate1.jl - A simple climate model based on
https://open.bu.edu/bitstream/handle/2144/28962/1705.07387v1.pdf?sequence=1&isAllowed=y
author: Bill Thompson
license: GPL
copyright: May 28, 2018
=#
using DifferentialEquations
using PyPlot
#=
The Maasch and Saltzman model
This function is called by the ODE solver.
arguments:
du - the derivatives of x,y,,z
u - a 3 element vector representing x,y,z from the paper
p - a tuple of the p, q, r, s arguments to the model
t - time
returns:
du - a list of the current values of the derivative of u
note:
The form of this function is dictated by the ODE solver, in this case
ODEProblem.
=#
function maasch_saltzman(du, u, p, t)
p1, q, r, s = p
du[1] = -(u[1] + u[2])
du[2] = (r - u[3]^2) * u[2] - p1 * u[3] + s * u[3]^2
du[3] = -q * (u[1] + u[3])
end
#=
plot_it - use matplotlib to plot the results of the ODE solver
arguments:
sol - the solver object returned by solve
=#
function plot_it(sol)
u1 = [sol.u[i][1] for i in 1:length(sol.u)]
u2 = [sol.u[i][2] for i in 1:length(sol.u)]
u3 = [sol.u[i][3] for i in 1:length(sol.u)]
fig, ax = subplots()
ax[:plot](sol.t, u1, label="Ice", color="black")
ax[:plot](sol.t, u2, label="CO2", color="red")
ax[:plot](sol.t, u3, label="NADW", color="blue")
ax[:legend](loc="lower left")
xlabel("Time [10K years]")
ylabel("x,y,z")
show()
end
function main()
# parameter values from figure 3.1
p = 1.0
q = 1.2
r = 0.8
s = 0.8
# initial conditions - I took these from the figure 3.1 so may not exactly
# match those used in the paper
u0 = [-0.4, 0.1, 0.9]
tspan = (0.0, 50.0) # time span
prob = ODEProblem(maasch_saltzman, u0, tspan, (p, q, r, s))
sol = solve(prob, AutoTsit5(Rosenbrock23()), reltol=1e-8, abstol=1e-8)
plot_it(sol)
end
main()
This code produces the same plot as the Python version.
I still think Julia has great potential as a scientific and statistical programming tool. However, I have been running into problems with it. I use Windows Subsystem for Linux (WSL) and have had trouble getting plotting to work. With WSL, I cannot get PyPlot to work at all despite trying some of the fixes I found online. I found colors don't work properly with GR or plotly. Some of these problems may come from WSL which can be a bit cranky.
Using Linux, PyPlot works properly, but GR and plotly don't display when used from a script rather than the Julia REPL despite using display(). I'm sure I am doing something wrong in each of these cases, but it's hard to spend too much time tracking down implementation problems when you can get results quickly in Python.
Module load time in Julia is a bigger problem. Runtime for an app like the one above is dominated by the amount of time it takes to load DifferentialEquations and PyPlot. This has been a problem since the early versions of Julia.. While Julia is fast once it gets loaded, small apps can't can't compete with Python or R in terms of speed if load times dominate the runtime.
The full Python and Julia code is available for download here.
No comments:
Post a Comment