Lecture 36: Snowball earth and models of Global Climate (2 days)

(previous, home, next)


References


Model set-up

The dynamic Budyko model for temperature \(T(y,t)\) where \(y\) is the normalized latitude (from 0 at the equator to 1 at the north pole, and \(t\) is time (in years).

\[R \frac{\partial}{\partial t} T{\left (y,t \right )} = Q s(y) \left(1 - \alpha{\left (y,y_{s} \right )} \right) - (A + B T) + C \left( \int_{0}^{1} T{\left (y,t \right )}\, dy - T{\left (y,t \right )} \right)\]

where \(A + B T = 202 + 1.9 T\), \(C = 1.6 B = 3.04\), \(\alpha(y, y _ s) = 0.32 + 0.3 H(y - y_s)\), \(H(x)\) is the Heaviside function, and \(s(y) = 1 - 0.241 (3 y^2 - 1)\).

Model analysis

Steady-states

[Show code]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
#!/usr/bin/env python2
"""
We are working with the dynamic Budyko model and governing law ...

>>> Eq(R*T(y,t).diff(t), Q*s(y)*(1-alpha(y,y_s)) - (A+B*T) + C*(Integral(T(y,t),(y,0,1)) - T(y,t)))


                                                  / 1                     \
      dT                                          | /\                    |
    R -- = Q (1 - a(y, y_s)) s(y) - (A + B T) + C | | T(y, t) dy - T(y, t)|
      dt                                          |\/                     |
                                                  \ 0                     /
"""


from pylab import *

@vectorize
def heaviside(x):
    return (0. if x < 0. else ( 1. if x > 0. else 0.5))

@vectorize
def zap(x):
    # 1 if x >= 0 otherwise NaN
    return ( float('nan') if x < 0. else 1. )


Q_real = 343.
#Q = 330.
A = 202.
B = 1.90
C = 1.6*B
TIce = -10. # average temperature at ice line

def s(y):
    # Latitudinal distribution approximation
    return  1 - 0.482*(3*y**2-1)*0.5

def alpha(y, ys):
    # albedo as a function of latitude and ice-line latitude
    return 0.32 + heaviside(y-ys)*0.3

def mid_alpha(ys):
    # this works only because alpha is a flat Heaviside
    return (alpha(0., ys) + alpha(1., ys))*0.5

def net_alpha(ys):
    # average global albedo, which is Integral(s(y)*alpha(y,ys),(y,0,1))
    return 0.62 - 0.3*ys*(1-0.241*(ys**2-1))

def average_T(ys, Q):
    return (Q*(1 - net_alpha(ys)) - A )/B

def T(y, ys, Q):
    # local steady-state temperature
    return (Q*s(y)*(1 - alpha(y, ys)) - A + C*average_T(ys, Q))/(B + C)


def midpoint_temperature(ys, Q):
    # at the midpoint of the ice line
    return ( Q*s(ys)*(1 - mid_alpha(ys) ) -A + C*average_T(ys, Q))/(B+C)

def critPower(ys):
    return ((B+C)*TIce + A - C*average_T(ys, Q_real))/(s(ys)*(1-mid_alpha(ys)))

def midTemp(ys):
    return average_T(ys, critPower(ys) )

def topTemperature(Q):
    # assume an ice-free world
    ys = 1.
    y = 1.
    return zap(T(y, ys, Q)-TIce)*average_T(ys, Q)

def bottomTemperature(Q):
    # assume an ice-covered world
    ys = 0.
    y = 0.+1e-5
    return zap(TIce-T(y, ys, Q))*average_T(ys, Q)



def main():
    figure(1, figsize=(12, 10))
    y = linspace(0, 1, 1000)
    print T(1e-3, 0., Q_real), T(0.9, 0., Q_real)
    print T(0., 1., Q_real), T(1., 1., Q_real)

    subplot(2, 2, 4)
    for ys in linspace(0, 1, 10):
        plot(y, T(y, ys, Q_real), 'k-')
    #plot(y, -10+0*y, 'r:', y, midpoint_temperature(y, Q_real), 'g--')
    plot(y, -10+0*y, 'r:')
    xlabel('Latitude (normalized)', fontsize=16)
    ylabel('Local temperature for each $y_s$', fontsize=16)
    ylim(-60, 30)

    subplot(2, 2, 3)
    plot(y, s(y), 'b-')
    xlabel('Latitude (normalized)', fontsize=16)
    ylabel('Latitudinal distribution $s(y)$', fontsize=16)

    subplot(2, 2, 1)
    plot(y, alpha(y, 0.5), 'g-')
    xlabel('Latitude (normalized)', fontsize=16)
    ylabel('Albedo if ice line $y_s= 0.5$', fontsize=16)

    subplot(2, 2, 2)
    plot(y, net_alpha(y), 'g-')
    xlabel('Latitude of ice line $y_s$ (normalized)', fontsize=16)
    ylabel(r'Global average albedo $\overline{\alpha}(y_s)$', fontsize=16)

    savefig('snowball.pdf')
    savefig('snowball.png')


    figure(4)
    for ys in linspace(0, 1, 10):
        plot(y, T(y, ys, Q_real), 'k-')
    plot(y, -10+0*y, 'r:', y, midpoint_temperature(y, Q_real), 'g--')
    xlabel('Latitude (normalized)', fontsize=18)
    ylabel('Local temperature for each $y_s$', fontsize=18)
    ylim(-60, 30)

    savefig('icelinesols.pdf')
    savefig('icelinesols.png')



    figure(2)
    ys = linspace(0, 1-1e-5, 100)
    Q = array([ critPower(i) for i in ys ])
    plot(Q, ys, 'k-', [343, 343], [0, 1], 'g:')
    xlim(250, 550)
    xlabel('Solar power (Q watts per meter square)', fontsize=18)
    ylabel('Stationary ice line lattitude ($y_s$)', fontsize=18)
    savefig('power_to_lattitude.png')
    savefig('power_to_lattitude.pdf')


    figure(3)
    mT = array([ midTemp(i) for i in ys])
    plot(Q, mT, 'b-')
    Q = linspace(250, 500, 1000)
    #hT = array([ highTemp(i) for i in Q])
    #lT = array([ lowTemp(i) for i in Q])
    #plot(Q, lT, 'k-', Q, hT, 'k-')
    plot(Q, topTemperature(Q), 'k-', Q, bottomTemperature(Q), 'k-')
    #plot(Q, T(1e-3, 0., Q), 'g-', Q, T(1.-1e-3, 0., Q), 'g-')
    #plot(Q, T(1e-3, 1., Q), 'r-', Q, T(1.-1e-3, 1., Q), 'r-')
    #plot([Q_real]*2, [-80, 100], 'b:')
    xlabel('Solar power (Q watts per meter square)', fontsize=18)
    ylabel('Global average temperature (degrees C)', fontsize=18)
    ylim(-60, 80)
    savefig('power_to_temp.png')
    savefig('power_to_temp.pdf')

    #show()


main()

Stability of steady-states using perturbation analysis

Relaxation oscillators and punctuated dynamics