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

Now what does this concentration/forcing trajectory mean for temperature (with no coupling between Carbon and climate) using Geoffrey et al. (2013) two-layer EBM to account for the deep ocean?

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

So, for ~65 more years it would rise perceptibly if we say that less than 0.5˚C/century is imperceptible.

Importantly, temperature does not come down after that, and continues to creep up for quite a while.

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