CO2 simulation

In [8]:
import matplotlib.pyplot as plt
import csv
import numpy as np
from scipy.optimize import curve_fit
import math
%matplotlib inline
plt.rcParams['figure.figsize'] = [16, 8]

Some definitions and constants

In [2]:
##Definitions

def nth_root(num,root):
    answer = num ** (1/root) 
    return answer

def KtoC(In):
    return In-273.15

def SQmeter(Reservoirs):
    a = []
    for i in Reservoirs:
        a.append(i[0]/i[6])
        print(i[5], i[0]/i[6])
    return a

#               km3    ,   GT    ,K ,J/kg*k,J        ,abv , height, area, J/s
Atmosphere = [5.0*10**9,5.0*10**6,260,1.0,1.3*10**21,"At",110000, 540000] # ??km
Surface_upper = [5.0*10**3,2.0*10**4,290,4.0,2.32*10**19,"Su",0.04, 125000] # 4cm
Surface_lower = [1.0*10**5,2.0*10**5,280,4.0,2.24*10**20,"Sl",80, 1250] # 80m
Ocean_upper = [3.8*10**7,3.8*10**7,280,4.0,4.256*10**22,"Ou",100, 380000] # 100m
Ocean_lower = [1.4*10**9,1.4*10**9,275,4.0,1.54*10**24,"Ol",4000, 350000] # 4km

##**~~-ASSUMPTIONS-~~**##
R1 = 6371.0*1000 # radius of earth in meter
H1 = 12.0*1000 # Height of atmosphere to work with
R2 = R1 + H1 # Radius of earth + troposphere
V = ((4*math.pi*R2**3)/3) - ((4*math.pi*R1**3)/3) # m^3 of air
S = 50*2*10**9 # Both ocean and land scavange +/- 50 gigatons of CO2 annualy
PPMinit = 0 # this is set as the first point in the dataset; initial concentration
T = 25 # Averaage 25 degres celsius
Rho = 1.836 # density of CO2 @ 25 C; kg/m3
print(V)
##**~~-ASSUMPTIONS-~~**##


Sfrac = Surface_upper[-1]/(Surface_upper[-1]+Ocean_upper[-1])
Ofrac = 1-Sfrac
Ta = Sfrac*Surface_upper[2]+Ofrac*Ocean_upper[2]


Reservoirs = [Atmosphere,Surface_upper,Surface_lower,Ocean_upper,Ocean_lower]
print(Atmosphere[1]*Atmosphere[2]*Atmosphere[3]*10**12)
print(Surface_upper[1]*Surface_upper[2]*Surface_upper[3]*10**12)
print(Surface_lower[1]*Surface_lower[2]*Surface_lower[3]*10**12)
print(Ocean_upper[1]*Ocean_upper[2]*Ocean_upper[3]*10**12)
print(Ocean_lower[1]*Ocean_lower[2]*Ocean_lower[3]*10**12)
print(1/(Ocean_lower[1]*Ocean_lower[3]*10**12/Ocean_lower[4]))
for i in Reservoirs:
    print(i[5],KtoC(i[2]))

    
areas = SQmeter(Reservoirs)

print(areas[1]+areas[3])
6.132309591141384e+18
1.3e+21
2.32e+19
2.24e+20
4.256e+22
1.54e+24
275.0
At -13.149999999999977
Su 16.850000000000023
Sl 6.850000000000023
Ou 6.850000000000023
Ol 1.8500000000000227
At 45454.545454545456
Su 125000.0
Sl 1250.0
Ou 380000.0
Ol 350000.0
505000.0

Simulation from 1958 to 2014

In [16]:
##**~~-IMPORTANT VARIABLES!-~~**##
Year = [] # 0
Csum = [] # 1
PPM = [] # 8
##**~~-IMPORTANT VARIABLES!-~~**##


FileName = "co2_annmean_mlo.csv"
#File contains:
#   * PPM of CO2 from 1958 to 2014
#   * emission of CO2 in tonnes from 1958 to 2014
#   * more detailed data regarding above posts


Cgf = [] # 2; gas fuel
Clf = [] # 3; liquid fuel
Csf = [] # 4; solid fuel
Cc = [] # 5; cement production
Cf = [] # 6; gas flaring
Cpc = [] # 7; per capita
PPMnh = [] # 9; ppm northern hemisphere
PPMsh = [] # 10; ppm southern hemisphere


#### Loading of data-file
with open(FileName, 'r') as csvfile:
    data = csv.reader(csvfile, delimiter=',')
    i = 0
    for row in data:
        if i != 0:
            Year.append(int(row[0]))
            Csum.append(float(row[1])*10**6) # convert from million metric tonnes to tonnes
            Cgf.append(row[2])
            Clf.append(row[3])
            Csf.append(row[4])
            Cc.append(row[5])
            Cf.append(row[6])
            Cpc.append(row[7])
            PPM.append(float(row[8]))
#             PPMnh.append(row[9]) #Data not available in IPCC dataset!
#             PPMsh.append(row[10])
        if i == 1:
            PPMinit = float(row[8])
        i += 1
##Done loading
        
##Test calculation        
m3CO2 = (Csum[-1]*1000)*Rho # Cubic meters of Co2 in 1751 anthropogenic sources
print(m3CO2)
print(Csum[-1]) # kg of CO2 anthropogenic sources
print(S) # scavenging
print((m3CO2/V)*10**6) # ppm of atmosphere


##Setting up plots
fig, ax1 = plt.subplots() # get a plot object
ax1.plot(Year,PPM,'b-',label="Measured PPM") # plot PPM on left axis
ax1.set_xlabel("Year") # set labels
ax1.set_ylabel("ppm CO2")

#CO2 simulator functions
def SimCo2(Csum,V,PPMinit):#Input: An array with emitted CO2 per timestep in Kg; volume of air, initial concentration
    PPMSim = [] # array to store result
    PPMtemp = PPMinit # array to store temporary calculation
    Vupdate = V # initial volume of air
    Err = [] # error based on PPM value from dataset
    Ta = [] # atmosphere temperature
    Tsu = []
    Tsl = []
    Tou = []
    Tol = []
    NR = Reservoirs
    for i,x in enumerate(Csum): # iterate through all values
        M3CO2 = ((x*1000)-S)*Rho # calculate cubic meters of CO2
        PPMtemp += ((M3CO2/Vupdate)*10**6) # calculate how much PPM this is of the total volume
        Vupdate += M3CO2 # update the volume with emitted volume
        try:
            Err.append((PPM[i]-PPMtemp)) # calculate error based on measurements
        except:
            print("err: ",i,PPMtemp)
        PPMSim.append(PPMtemp) # add result to array        
    return PPMSim, Vupdate, Err # return values

### RUN simulation
PPMSim, Vupdate, Err = SimCo2(Csum,V,PPMinit)
##Done

## calculate total air mass growth
print("Vupdate:{} - Vold: {}, %growth: {}".format(Vupdate,V,(Vupdate/V)*100))    

# Plot simulated PPM on the graph, left axis
ax1.plot(Year,PPMSim,'r--',label="Simulated PPM")

## Function to fit on past CO2 emissions. This will be used to generate scenarios
def func(x, a, b, c):
    return a*x**2 + b*x + c


poptemsi, pcov = curve_fit(func, Year[3:], Csum[3:]) # fit on total period
popt, pcov = curve_fit(func, Year, Csum) # fit on first 100 year # changed dataset not so extensive
print(popt)

Emsi = []
Embau = []
Embe = []
Fy = []

## Parabolic function to demonstrate best case scenario
def parab(x,a,b):
    return a*x**2 + b

# manual parabolic fit
poptparab = [-(1/18*10**7),1.2*10**10]


# Get PPM value @1970, when "prediction" is initiated
PPMSimInit = PPM[-(2014-1970)]

## Predict from 1970 - untill 2100:
for i in range(2101-1970):
    Embau.append(func(1970+i,*popt)) #Emmissions Business as usual 
    Emsi.append(func(1970+i,*poptemsi)) # Emissions Several improvements
    Embe.append(parab(-120+i,*poptparab)) # Emissions Best case scenario
    Fy.append(1970+i) # Add newly generated years to "Future Year" array
plt.grid()
plt.legend()
18093780000000.0
9855000000.0
100000000000
2.9505653181858142
Vupdate:6.132893543793388e+18 - Vold: 6.132309591141384e+18, %growth: 100.00952255660489
[ 5.96390954e+04 -1.15571477e+08  1.00000000e+00]
/home/jeroen/.local/lib/python3.6/site-packages/scipy/optimize/minpack.py:808: OptimizeWarning: Covariance of the parameters could not be estimated
  category=OptimizeWarning)
Out[16]:
<matplotlib.legend.Legend at 0x7f9965f1c438>

Loading of IPCC Co2 Scenario's (untill 2100)

In [17]:
Fy = []
B1 = []
A1B = []
A1FI = []
with open("CO2 scenarios.csv", 'r') as csvfile:
    data = csv.reader(csvfile, delimiter=',')
    i = 0
    for row in data:
        if i != 0:
            Fy.append(int(row[0]))
            B1.append(float(row[1])*10**8) # convert from million metric tonnes to tonnes
            A1B.append(float(row[2])*10**8)
            A1FI.append(float(row[3])*10**8)
        i += 1

PPMSimInit = PPM[-14]
PPMSimbau, Vupdate, Erra = SimCo2(B1,V,PPMSimInit)
PPMSimsi, Vupdate, Erri = SimCo2(A1B,V,PPMSimInit)
PPMSimbe, Vupdate, Erre = SimCo2(A1FI,V,PPMSimInit)
err:  56 543.1369852355193
err:  57 546.3429285494326
err:  58 549.5192577551521
err:  59 552.6659731327202
err:  60 555.783074965533
err:  61 558.8535908072921
err:  62 561.8775210892652
err:  63 564.8548662392001
err:  64 567.7856266732426
err:  65 570.6698028037201
err:  66 573.5073950355489
err:  67 576.2984037674323
err:  68 579.0428293918606
err:  69 581.7406722942131
err:  70 584.391932853657
err:  71 586.9976740787032
err:  72 589.5578963269397
err:  73 592.0725999500361
err:  74 594.5417852937434
err:  75 596.9654526969967
err:  76 599.3436024928133
err:  77 601.6762350082921
err:  78 603.9633505637169
err:  79 606.2049494734536
err:  80 608.4010320459507
err:  81 610.5543972881696
err:  82 612.6650454778818
err:  83 614.7329768876828
err:  84 616.758191784094
err:  85 618.7406904284614
err:  86 620.6804730769549
err:  87 622.5775399796705
err:  88 624.4318913818284
err:  89 626.2435275225755
err:  90 628.0124486358828
err:  91 629.7444318439718
err:  92 631.4394773398292
err:  93 633.097585312647
err:  94 634.7187559478231
err:  95 636.3029894260628
err:  96 637.8502859242772
err:  97 639.360645615583
err:  98 640.834068668405
err:  99 642.2705552473736
err:  100 643.6701055133253
err:  56 585.3189492412016
err:  57 590.1066490710718
err:  58 594.9041063668756
err:  59 599.7113209881659
err:  60 604.5282927912183
err:  61 609.3354409217245
err:  62 614.1327655183468
err:  63 618.9202667254433
err:  64 623.6979446810946
err:  65 628.465799529077
err:  66 633.2238314068895
err:  67 637.9720404577268
err:  68 642.7104268185063
err:  69 647.4389906318411
err:  70 652.1577320340671
err:  71 656.8642114584509
err:  72 661.5584290751763
err:  73 666.2403850569672
err:  74 670.910079576095
err:  75 675.5675128073716
err:  76 680.2126849191701
err:  77 684.845596082404
err:  78 689.4662464675342
err:  79 694.0746362475622
err:  80 698.6707655890509
err:  81 703.2455054439451
err:  82 707.7988561031651
err:  83 712.3308178622418
err:  84 716.8413910123379
err:  85 721.3305758402482
err:  86 725.7983726373785
err:  87 730.2447816907667
err:  88 734.6698032830826
err:  89 739.0734377016075
err:  90 743.4556852292548
err:  91 747.8125164783268
err:  92 752.14393178458
err:  93 756.4499314788374
err:  94 760.7305158869882
err:  95 764.9856853389667
err:  96 769.2154401597741
err:  97 773.4197806694777
err:  98 777.5987071921907
err:  99 781.7522200470927
err:  100 785.8803195514233
err:  56 616.6768018681867
err:  57 623.8113613148099
err:  58 631.0314554702212
err:  59 638.3370824886811
err:  60 645.7282405025289
err:  61 653.1700159504235
err:  62 660.6624077044806
err:  63 668.2054146291473
err:  64 675.7990355782085
err:  65 683.4432694037675
err:  66 691.1381149472671
err:  67 698.8835710424829
err:  68 706.6796365125315
err:  69 714.5263101788497
err:  70 722.4235908522162
err:  71 730.3655742851453
err:  72 738.3522594170585
err:  73 746.3836451754136
err:  74 754.4597304876768
err:  75 762.5805142693519
err:  76 770.7459954359517
err:  77 778.9561728910279
err:  78 787.2110455351494
err:  79 795.5106122659031
err:  80 803.8548719689159
err:  81 812.2161257512469
err:  82 820.5943731863321
err:  83 828.9896138467483
err:  84 837.4018473012204
err:  85 845.8310731235997
err:  86 854.2772908838858
err:  87 862.7405001512193
err:  88 871.2207004908894
err:  89 879.7178914733117
err:  90 888.2320726650506
err:  91 896.7395448680467
err:  92 905.2403082529559
err:  93 913.7343629932886
err:  94 922.2217092564327
err:  95 930.7023472126303
err:  96 939.1762770319858
err:  97 947.6434988844657
err:  98 956.1040129428907
err:  99 964.5578193739593
err:  100 973.0049183472246

Plotting of past simulation + measurement and predicted scenario CO2 concentrations

In [15]:
plt.plot(Year,PPM,'b-',label="Measured")
plt.plot(Year,PPMSim,'r--',label="Calculated based on emissions")
plt.plot(Fy,PPMSimbau,"g-.",label="B1 Scenario")
plt.plot(Fy,PPMSimsi,"y-.",label="A1B Scenario")
plt.plot(Fy,PPMSimbe,"m-.",label="A1FI Scenario")
plt.fill_between(Fy,PPMSimbau,PPMSimsi, facecolor='green',alpha = 0.5)
plt.fill_between(Fy,PPMSimsi,PPMSimbe, facecolor='red',alpha = 0.5)
plt.gca().set_xlabel("Year") # set labels
plt.gca().set_ylabel("PPM CO2")
plt.grid()
#plt.yscale('log')
plt.legend()
plt.show()

Emissions in kg CO2 for the 3 scenario's

In [13]:
plt.plot(Year,Csum,label="Measured CO2")
plt.plot(Fy,B1,"g-.",label="B1 Scenario")
plt.plot(Fy,A1B,"y-.",label="A1B Scenario")
plt.plot(Fy,A1FI,"m-.",label="A1FI Scenario")
plt.fill_between(Fy,B1,A1B, facecolor='green',alpha = 0.5)
plt.fill_between(Fy,A1B,A1FI, facecolor='red',alpha = 0.5)
plt.gca().set_xlabel("Year") # set labels
plt.gca().set_ylabel("kg CO2")
plt.grid()
plt.legend()
plt.show()

Prediction error over period: 1958 - 2014 per year

In [20]:
plt.plot(Year,Err)
Out[20]:
[<matplotlib.lines.Line2D at 0x7f9963705400>]