In [1]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
import matplotlib.patches as mpatches
import matplotlib.lines as mlines
from matplotlib import cm
import matplotlib

%matplotlib inline

In [2]:
startingEmissionRate = 10 #e9 #10 GtC/year
startingConcentration = 400 #ppm
PPM_TO_GTC = 2.13 #e9

percentPerYearReduction = 0.015 #1.5%/yr
totalPercentReduction = .99    #until it's been reduced 99% from starting rate

yearsIntoTheFuture = 1000

yearsForAxis = np.linspace(2016,2016+yearsIntoTheFuture-1,yearsIntoTheFuture)

In [3]:
emit = np.zeros(yearsIntoTheFuture)
emit[0] = startingEmissionRate
for i in xrange(len(emit)):
if i>0:
emit[i]=emit[i-1]-(percentPerYearReduction*emit[i-1])

stabilizeEmissionsAt = startingEmissionRate - (totalPercentReduction*startingEmissionRate)
print "stabilize emissions at: ", stabilizeEmissionsAt, ' GtC/year'
indexOfStabilized = len(emit[emit>stabilizeEmissionsAt])
emit[indexOfStabilized+1::] = stabilizeEmissionsAt
#emit
plt.plot(yearsForAxis,emit)
plt.xlabel('Year') #s from Start (Present)')
plt.ylabel('GtC/year emitted')
plt.title('Emissions')

stabilize emissions at:  0.1  GtC/year

Out[3]:
<matplotlib.text.Text at 0x7f6e79c3a810>

## This response function is based on the multi-model mean in Joos et al. (2013). It assumes a 100GtC pulse against a background of 389ppm.¶

### These are the constants needed.¶

In [4]:
# coefficients
a0 = 0.2173 #unitless
a1 = 0.2240
a2 = 0.2824
a3 = 0.2763
# timescales
tau1 = 394.4 #units of years
tau2 = 36.54
tau3 = 4.304

In [5]:
yearsToUse = yearsForAxis #
IRF = np.zeros(len(yearsToUse))
IRFstart = np.zeros(len(yearsToUse))
for year in xrange(len(yearsToUse)):
IRFstart[year] = a0 + (a1*np.exp(-(year+1)/tau1)) #put only the longer timescales to work on the "initial" pulse
IRF[year] = a0 + (a1*np.exp(-(year+1)/tau1)) + (a2*np.exp(-(year+1)/tau2)) + (a3*np.exp(-(year+1)/tau3))
IRFstart = IRFstart/IRFstart[0]
plt.plot(IRF)
plt.plot(IRFstart)
plt.xlabel('Years into Future from Pulse')

Out[5]:
<matplotlib.text.Text at 0x7f6e7c319410>

## Stagger the start of the IRF in respone to new "pulse" each year emissions continue, then sum.¶

In [6]:
yearsIntoTheFuture = len(yearsToUse)
GtCinAtmos = np.zeros((yearsIntoTheFuture,yearsIntoTheFuture))
GtCinAtmos[0,:] = startingConcentration * PPM_TO_GTC #starting from now

for i in xrange(yearsIntoTheFuture):
if i>0:
GtCinAtmos[i,i::] = emit[i]*IRF[i::]
elif i==0:
GtCinAtmos[0,:] = (((startingConcentration-289.0)* PPM_TO_GTC) * IRFstart)

# you have to treat the start as one big pulse above preindustrial, and add preindustrial back onto the end.
GtCinAtmos = GtCinAtmos.sum(0)/PPM_TO_GTC + 289.0 #in ppm
#GtCinAtmos = GtCinAtmos.sum(0)/PPM_TO_GTC

ax = plt.subplot(111)
plt.plot(yearsForAxis,GtCinAtmos)
plt.xlabel('Year')
plt.ylabel('ppm')
plt.title('CO2 Concentration')

ax3 = ax.twinx()
plt.plot(yearsForAxis,emit,'r')
ax3.set_ylabel('GtC/year emitted',color='r')
ax3.tick_params('y', colors='r')


## Where is the maximum?¶

In [7]:
print 'Year ', np.int(yearsForAxis[np.nonzero(GtCinAtmos == GtCinAtmos.max())[0]][0])

Year  2140


## And what concentration does it max out at?¶

In [8]:
print GtCinAtmos.max(), ' ppm'

486.516606785  ppm


## The temperature based on the equilibrium climate sensitivity would follow the concentration.¶

In [9]:
def getForcing(time):
return (F_4xCO2/np.log(4)) * np.log(GtCinAtmos[time]/preindustrial)

def getTequilAnom(forcing):
return forcing/feedbackParam

preindustrial = 289.0 #ppm

#model mean:
F_4xCO2 = 6.9 #W/m2
feedbackParam = 1.13

In [10]:
ax = plt.subplot(111)
plt.plot(yearsToUse,getTequilAnom(getForcing(np.arange(yearsIntoTheFuture))),'b')
plt.xlabel('Year')
ax.set_ylabel('Temperature Anom from Preindustrial',color='b')
ax.tick_params('y', colors='b')
plt.title('Equilibrium Temperature based on Concentration')

ax2 = ax.twinx()
plt.plot(yearsToUse,GtCinAtmos,'g')
ax2.set_ylabel('Concentration (ppm)',color='g')
ax2.tick_params('y', colors='g')


## We have to integrate forward from preindustrial.¶

### First get the preindustrial concentrations from RCP site, interpolate to yearly, and stitch together with projection from present forward.¶

In [12]:
preindustrialConc = pd.read_csv('justCO2_1850_2005.csv').transpose()

preindustrialConc.index = pd.to_datetime(preindustrialConc.index)
startYearNow = pd.DataFrame([GtCinAtmos[0]],index=pd.date_range('2016-1-1',periods=1,freq='1AS'))
preindustrialConc = preindustrialConc.append(startYearNow)
preindustrialConc = preindustrialConc.resample('1AS',axis=0).asfreq() #.interpolate('linear')
preindustrialConc = preindustrialConc.interpolate('linear')
preindustrialConc[0].plot()
plt.title('Preindustrial CO2 concentration used in RCPs')
plt.ylabel('ppm')
plt.xlabel('Year')
yearsPreindustrial = np.asarray(preindustrialConc.index.year)
preindustrialConc = np.asarray(preindustrialConc.transpose())[0]

In [13]:
GtCallTime = np.append(preindustrialConc,GtCinAtmos[1::])
yearsAllTime = np.append(yearsPreindustrial,yearsForAxis[1::])
plt.plot(yearsAllTime,GtCallTime)
plt.title('Concentration, historical and future scenario')
plt.ylabel('ppm')
plt.xlabel('Year')

Out[13]:
<matplotlib.text.Text at 0x7f6e795e4ad0>
In [14]:
from scipy.integrate import odeint, ode

#define more model parameters
coupling = 0.74
heatCapacityTop = 7.3
heatCapacityDeep = 91

A = [[-(feedbackParam+coupling)/heatCapacityTop,coupling/heatCapacityTop],
[coupling/heatCapacityDeep,-coupling/heatCapacityDeep]]

def getForcing(time):
return (F_4xCO2/np.log(4)) * np.log(GtCallTime[time]/preindustrial)

def EBM2layer(y,t,args):
b = [getForcing(t)/heatCapacityTop,0]
return np.dot(A,y) + b

y0 = [0,0] #starting from 1850
t = np.arange(len(yearsAllTime))
b = [getForcing(0)/heatCapacityTop,0]
sol = odeint(EBM2layer, y0, t, args=(1,))

print sol.shape
plt.plot(sol[:,0])
plt.plot(sol[:,1])
plt.legend(['"air"','deep ocean'],loc='lower right')

(1166, 2)

/home/disk/grease/naomi/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:12: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future

Out[14]:
<matplotlib.legend.Legend at 0x7f6e7347dc50>
In [15]:
ax = plt.subplot(111)
plt.plot(yearsAllTime,sol[:,0],'b')
plt.xlabel('Year')
ax.set_ylabel('Temperature Anom from Preindustrial',color='b')
ax.tick_params('y', colors='b')
plt.title('Transient Temperature based on Heat Capacities')

ax2 = ax.twinx()
plt.plot(yearsAllTime,GtCallTime,'g')
ax2.set_ylabel('Concentration (ppm)',color='g')
ax2.tick_params('y', colors='g')


## temperature maximum is:¶

In [16]:
print "%.2f" %np.max(sol[:,0]), '˚C above preindustrial'
print 'maximum in year ', np.int(yearsAllTime[np.nonzero(sol[:,0] == sol[:,0].max())[0]][0])

1.97 ˚C above preindustrial
maximum in year  2321


## To be quantitative about it, we could say that when the rate becomes imperceptible, that's when it has stabilized.¶

In [20]:
threshold = 0.5 #degrees/century
threshold = threshold/100.0 #degrees/year

plt.figure(figsize=[15,5])
ax = plt.subplot(111)
justAir = sol[:,0]
deriv = [justAir[i+1]-justAir[i] for i in xrange(len(justAir)-1)]
deriv = np.asarray(deriv)
plt.plot(yearsAllTime[1::],deriv)
ax.set_ylabel('First Derivative of Temperature Anom',color='b')
ax.tick_params('y', colors='b')

print "Until the year ", int(yearsAllTime[np.nonzero(deriv > threshold)][-1])

ax2 = ax.twinx()
plt.plot(yearsAllTime,sol[:,0],'g')
ax2.set_ylabel('Temperature Anom from Preindustrial',color='g')
ax2.tick_params('y', colors='g')

Until the year  2081


## Finally, also include low and high climate sensitivity estimates.¶

### Using parameters that fit the two individual models with highest and lowest climate sensitivity. This is only a lower bound of what the range could be.¶

In [21]:
#multi-model mean, as above:
F_4xCO2 = 6.9 #W/m2
feedbackParam = 1.13
coupling = 0.74
heatCapacityTop = 7.3
heatCapacityDeep = 91

#low: GISS-E2R
feedbackParam = 1.7
F_4xCO2 = 7.3 #W/m2

A = [[-(feedbackParam+coupling)/heatCapacityTop,coupling/heatCapacityTop],
[coupling/heatCapacityDeep,-coupling/heatCapacityDeep]]

y0 = [0,0] #starting from 1850
t = np.arange(len(yearsAllTime))
low = odeint(EBM2layer, y0, t, args=(1,))

#high: CSIRO-Mk3.6.0
feedbackParam = 0.61
F_4xCO2 = 5.1 #W/m2

A = [[-(feedbackParam+coupling)/heatCapacityTop,coupling/heatCapacityTop],
[coupling/heatCapacityDeep,-coupling/heatCapacityDeep]]

y0 = [0,0] #starting from 1850
t = np.arange(len(yearsAllTime))
high = odeint(EBM2layer, y0, t, args=(1,))

plt.figure(figsize=[15,5])
ax = plt.subplot(111)
plt.plot(yearsAllTime,sol[:,0],'b')
plt.xlabel('Year')
ax.set_ylabel('Temperature Anom from Preindustrial') #,color='b')
#ax.tick_params('y', colors='b')
plt.title('Transient Temperature with High and Low Climate Sensitivity')
plt.plot(yearsAllTime,low[:,0],'b--')
plt.plot(yearsAllTime,high[:,0],'b--')

/home/disk/grease/naomi/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:12: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future

Out[21]:
[<matplotlib.lines.Line2D at 0x7f6e72e5b6d0>]
In [22]:
def getFarhrenheit(c):
return c*(9.0/5.0) #+ 32.0

print 'low climate sensitivity: ', "%.2f" %low[:,0].max(), '˚C or ', "%.2f" %getFarhrenheit(low[-1,0]), '˚F'
print 'mean climate sensitivity: ', "%.2f" %sol[:,0].max(), '˚C or ', "%.2f" %getFarhrenheit(sol[-1,0]), '˚F'
print 'high climate sensitivity: ', "%.2f" %high[:,0].max(), '˚C or ', "%.2f" %getFarhrenheit(high[-1,0]), '˚F'

low climate sensitivity:  1.44 ˚C or  2.20 ˚F
mean climate sensitivity:  1.97 ˚C or  3.15 ˚F
high climate sensitivity:  2.52 ˚C or  4.32 ˚F