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

(previous, home, next)

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

 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()