John Carlos Baez on background
%hide
%latex
In the [Zaliapin-Ghil paper](http://www.nonlin-processes-geophys.net/17/113/2010/npg-17-113-2010.pdf) they take equation
Then in section 3.3.2 they consider the "slow, quasi-adiabatic" dynamics of this equation. That's just a fancy of way of saying they assume
$$ \frac{d T}{d t} = 0$$
and solve for $T$ as a function of $\mu$, holding other things fixed at specific values. Here $\mu Q_0$ is the '"insolation", that is, very roughly, the amount of sunlight reaching the Earth. $Q_0$ is some "standard" amount of insolation, and $\mu$ is some parameter near $1$ which they call the "fractional insolation change". They get this graph:
So, we see that if $\mu$ is low, there is just one solution for $T$, which we could call the "Snowball Earth" solution. As we increase $\mu$ we reach a value there there are three solutions: two stable ones drawn as solid curves (the "Snowball Earth" solution and the "Hot Earth" solution) together with an intermediate unstable one drawn as a dashed curve. As we increase $\mu$ even more there is only one solution: the "Hot Earth" solution.
|
# steady state in ZG10 and their constant values from 3.2
reset()
Q0 = 342.5; c1 = 0.15; c2 = 0.7; m = 0.4
T0= 1.9e-15 ^(-1/6) # this was wrong in ZG10 paper
sigma_sb= 5.6697e-8; Tc = 273.0
#plot ZG10 figure 5 to see if it's correct
var('T')
def zg16_mu(T,kappa=0.05):
alpha = c1 + c2*(1-tanh(kappa*(T-Tc)))/2
g = 1 - m*tanh((T/T0)^6)
return (sigma_sb*T^4 *g)/Q0*(1-alpha)
line([[zg16_mu(t),t] for t in srange(100., 400., 1.)])
|
|
%r
# Nathan Urban - Updated R code
Q0 = 342.5; c1 = 0.15; c2 = 0.7; m = 0.4
sigma = 5.6697e-8; Tc = 273; kappa = 0.05
T0 = (1.9e-15)^(-1/6) # correct
##T0 = (1.9e-15)^(1/6) # ZG10 bug
T = seq(100, 400, by=1)
# fractional insolation change
mu = function(T)
{
alpha = c1 + c2*(1-tanh(kappa*(T-Tc)))/2
lhs = Q0*(1-alpha)
g = 1 - m*tanh((T/T0)^6)
rhs = sigma*T^4 * g
mu = rhs/lhs # mu Q0(1-alpha) = sigma g T^4
}
# EBM bifurcation diagram
plot(mu(T), T, xlim=c(0,3), ylim=c(100,400), type="l", lwd=2, col="blue", xlab="Fractional insolation change (mu)", ylab="Temperature (T)", main="Zaliapin and Ghil (2010), Fig. 5 (corrected)", xaxs="i", yaxs="i")
abline(h=seq(100, 400, by=50), lty="dotted", col="gray")
abline(v=seq(0, 2.5, by=0.5), lty="dotted", col="gray")
# modern Earth state
T = 300
points(mu(T), T, pch=21, col="black", bg="green", cex=2)
Error in png() : X11 is not available Error in png() : X11 is not available Error in png() : X11 is not available Error in png() : X11 is not available Error in png() : X11 is not available Error in png() : X11 is not available Error in png() : X11 is not available Error in png() : X11 is not available |
%r
# fig 4 by Nathan Urban
Q0 = 342.5; c1 = 0.15; c2 = 0.7; m = 0.4
sigma = 5.6697e-8; Tc = 273; kappa = 0.05
T0 = (1.9e-15)^(-1/6) # correct
##T0 = (1.9e-15)^(1/6) # ZG10 bug
T = seq(150, 400, by=1)
# incoming radiation
Ri = function(T, mu)
{
alpha = c1 + c2*(1-tanh(kappa*(T-Tc)))/2
Ri = mu * Q0 * (1-alpha)
}
# outgoing radiation
Ro = function(T)
{
g = 1 - m*tanh((T/T0)^6)
Ro = sigma * g * T^4
}
# plot outgoing radation as a function of temperature
plot(T, Ro(T), xlim=c(150,400), ylim=c(0,630), type="l", lty="dashed", lwd=2, col="darkgreen", xlab="Temperature, T", ylab="Radiation, R", main="Zaliapin and Ghil (2010), Fig. 4", xaxs="i", yaxs="i")
# plot incoming radiation, for different fractional insolations
lines(T, Ri(T, mu=0.5), lwd=2, col="blue")
lines(T, Ri(T, mu=1.0), lwd=2, col="blue")
lines(T, Ri(T, mu=2.0), lwd=2, col="blue")
Error in png() : X11 is not available Error in png() : X11 is not available Error in png() : X11 is not available Error in png() : X11 is not available Error in png() : X11 is not available Error in png() : X11 is not available Error in png() : X11 is not available Error in png() : X11 is not available |
Investigation
%hide
%latex
Next, we might take the same equation, but treat $\mu$ not as a constant but as varying as a function of time. The simplest interesting case would be to let $\mu$ oscillate around an average value $\mu_0$:
$$\mu(t) = \mu_0 + a \cos(\omega t) $$
It'll be interesting if the oscillations are big enough that:
* when $\mu$ is large the only stable solution (from the "adiabatic" analysis in my last comment) is the "Hot Earth" solution, but
* when $\mu$ is small the only stable solution is the "Snowball Earth" solution.
Then if we start with $T(0)$ "hot" and the oscillation frequency $\omega$ is not too fast, the temperature $T(t)$ should stay "Hot" until $\mu$ gets sufficiently small, at which point the Earth should "snap" rather suddenly to the "Snowball" solution... and then when $\mu$ gets big enough it should "snap back" to the "Hot" solution.
The next step would be to add noise, but it would be nice to see this first. In particular, it would be nice if we could agree on some interesting values of all these parameters $c, Q_0, c_1, c_2, \kappa, T_c, \sigma, m, T_0$, and see what this graph looks like for those parameters:
That'll tell us interesting values for the parameters $\mu_0$ and $a$ in this formula:
$$\mu(t) = \mu_0 + a \cos(\omega t) $$
Of course all of this is just an oversimplified model - not much to do with the actual Earth, except in a rough way. Things get a bit more realistic when we reach the [Crucifix-Rougier model].
|
|
|
# Just some draft plots and testing time varying mu. todo:
# Next try to do a bf-plot of this dynamic version
reset()
# Add a function that captures the constants/equations above. returns two functions mu(t) and net radiation
def zg_num(kappa, mu0, tme, ampl, omega):
Q0 = 342.5; c1 = 0.15; c2 = 0.7; m = 0.4; T0= 1.9e-15^(-1/6)
sigma_sb= 5.6697e-8; Tc = 273.0
t,T = var('t,T')
mu = (mu0 + ampl*sin(omega*t))
zgn = mu*(Q0*(1 - (c1 + c2*(1 - tanh(kappa*(T - Tc)))/2.0)) - sigma_sb*T^4 *(1 - m* tanh((T/T0)^6)))
#print "mu(t) = %s"% mu
return mu,zgn
@interact
def sensanalyze(kappa = (0.05,1.0,0.05), mu0 = (0.1,2.5,0.1), tme= (0.1,100.0,.1),omega= (0.0,20.0,.1),ampl = (0.01,5.0,0.05)):
html('<h2 align=center><font color="darkred">Kappa = %g Mu0= %g</font></h2>'%(kappa,mu0))
mu,zg = zg_num(kappa,mu0,tme,ampl,omega)
mup = plot(mu,(t,0,tme),rgbcolor="darkred",ymin=-3,ymax=2,legend_label="Function: $\mu(t)$")
mup.show(axes=True, figsize=[7,2])
# just use line() instead
#bifp = list_plot(zgn(t,T))
mu3d = plot3d(zg,(t,0.,100.),(T,150,370))
mu3d.show(figsize=[7,2])
Click to the left again to hide and once more to show the dynamic interactive window |
# older version by Eric
# updated by staffan as an array computation assuming T and t are numpy arrays
import numpy as np
import scipy
import matplotlib
from pylab import *
sigmab = 5.6697e-8; Q0 = 342.5; kappa= 0.05
c1 = 0.15; c2 = 0.7; Tc = 273.0
m = 0.4; T0 = (1.9e-15)**(-1/6.)
a = .5; omega = .1; c = 1.
mu0 = 2.
def alpha(T):
return c1 + .5*c2*(1 - np.tanh(kappa*(T - Tc)))
def g(T):
return 1 - m*np.tanh((T/T0)**6)
def Ro(T):
return sigmab*g(T)*T**4
def Ri(T):
return mu*Q0*(1 - alpha(T))
def mu_eq(T):
return Ro(T)/Ri(T)
def get_mu(t):
mu = 1.*t;
mu[t < 0] = mu0;
mu[t >= 0] = mu0+a*np.sin(omega*t[t >=0])
return mu
T_init = 300.;
#"mu_init = mu_eq(T_init);
nt = 100.
dt = 2*pi/(nt*omega)
nmin = -nt
nmax = 9*nt
tiv = np.arange(nmin, nmax, dt)
#temporary set mu= 2.
# list_plot(alpha(tiv)) looks good
mu = get_mu(tiv)
list_plot(mu)
Rnet = 1.*tiv
T = tiv.copy()
Rnet = Ri(T)-Ro(T)
list_plot(Rnet)
subplot(311)
plot(tiv, mu,'b-', lw=2)
grid(True)
xlabel('Time')
ylabel('mu')
subplot(312)
plot(tiv, T,'b-', lw=2)
grid(True)
xlabel('Time')
ylabel('Temperature')
subplot(313)
Tp = arange(100.0, 400.0+1.0, 1.0)
plot(Tp, mu_eq(sigma,Q0,kappa,Tp), 'b-', T, mu, 'r-', lw=2)
grid(True)
xlabel('Temperature')
ylabel('mu')
show()
Traceback (click to the left of this block for traceback) ... TypeError: mu_eq() takes exactly 1 argument (4 given) Traceback (most recent call last):
File "", line 1, in <module>
File "/tmp/tmp9D27Fq/___code___.py", line 70, in <module>
plot(Tp, mu_eq(sigma,Q0,kappa,Tp), 'b-', T, mu, 'r-', lw=_sage_const_2 )
TypeError: mu_eq() takes exactly 1 argument (4 given)
|
|
|