Plotting workbook

Plotting workbook to accompany the paper Carbon Capture and Storage Energy Systems vs. Dispatchable Renewables for Climate Mitigation: A Bio-physical Comparison by Sgouridis, Dale, Csala, Chiesa & Bardi



In [1]:

    
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.font_manager as font_manager
plt.style.use('seaborn-darkgrid')
%matplotlib inline
prop = font_manager.FontProperties('Segoe UI')
from sympy import *
init_printing()

Eq. 6

Define EROI-CCS function (Eq. 6 in the paper)



In [2]:

    
var('EROI_CCS EROI f_op');



In [3]:

    
eq6=Eq(EROI_CCS,EROI/(1+f_op*EROI))
eq6









    Out[3]:





$$EROI_{CCS} = \frac{EROI}{EROI f_{op} + 1}$$

Eq. 2

Define dispatchable EROI for renewables with storage (Eq. 2 in the paper)



In [4]:

    
phi = Symbol('phi')
eta = Symbol('eta')
var('EROI_disp EROI ESOI phi eta');



In [5]:

    
eq2=Eq(EROI_disp,((1-phi)+(eta*phi))/((1/EROI)+(eta*phi/ESOI)))
eq2









    Out[5]:





$$EROI_{disp} = \frac{\eta \phi - \phi + 1}{\frac{\eta \phi}{ESOI} + \frac{1}{EROI}}$$

Figure 3 & Table 3

Now we generate the contour lines for Figure 3, using the examples from Table 3.

We have two tyes of renewable energy plants, PV and Wind, each with two cases of EROI, taken from the lower and upper boundaries encountered in studies.

We have assumed a storage technology mixes taken from Table 1 of Barnhart et al 2013. Here we have the ESOI values reported for 7 energy storage storage technologies, 5 types of electrochemical storage (Batteries), Compressed air (CAES) and pumped hydro (PHS).

We have considered diffrent energy storage technology mixes for our PV and Wind examples. For PV, we have taken an equal mix of all 5 electrochemical storage technologies, without CAES or PHS (resulting in the pv share vector of [0.2,0.2,0.2,0.2,0.2,0,0]). For Wind, we have considered 50% battery storage (as an equal 10% from each technology), and 25-25% of CAES and PHS, resulting in the wind share vector of [0.1,0.1,0.1,0.1,0.1,0.25,0.25]).



In [6]:

    
etas=[0.90,0.75,0.90,0.75,0.60,0.70,0.85] #from (Barnhart et al 2013)
S=[32,20,5,10,9,797,704]  #from (Barnhart et al 2013)
shr0=[0.2,0.2,0.2,0.2,0.2,0,0] #assumed shares for considered PV exmaples in this paper
shr1=[0.1,0.1,0.1,0.1,0.1,0.25,0.25]  #assumed shares for considered Wind exmaples in this paper
eta_pv=sum([etas[i]*shr0[i] for i in range(len(etas))]) #resulting composite eta for pv
eta_wind=sum([etas[i]*shr1[i] for i in range(len(etas))]) #resulting composite eta for wind
S_pv=sum([S[i]*shr0[i] for i in range(len(etas))])  #resulting composite S for pv
S_wind=sum([S[i]*shr1[i] for i in range(len(etas))])  #resulting composite S for wind
print 'Composite eta for pv',eta_pv
print 'Composite eta for wind',eta_wind
print 'Composite ESOI for pv',S_pv
print 'Composite ESOI for wind',S_wind









    



Composite eta for pv 0.78
Composite eta for wind 0.7775
Composite ESOI for pv 15.2
Composite ESOI for wind 382.85

Construct test cases: Renewables



In [22]:

    
#define pretty output formatter for sympy equation
def format_eq(eq):
    return [solve(eq[i])[0] for i in range(len(eq))]
    #return [float((repr(eq[i])[2:]).split(',')[1].strip()[:-1]) for i in range(len(eq))]



In [23]:

    
eroi_pv=[10,40] #EROEIe for case [F, G] in Table 3 of the paper, PV plant
eroi_disp_pv=format_eq([eq2.subs([[EROI,i],[phi,0.4],[eta,eta_pv],[ESOI,S_pv]]) for i in eroi_pv])
eroi_wind=[10,33] #EROEIe for case [H, I] in Table 3 of the paper, Wind plant
eroi_disp_wind=format_eq([eq2.subs([[EROI,i],[phi,0.4],[eta,eta_wind],[ESOI,S_wind]]) for i in eroi_wind])



In [25]:

    
print 'PV [F,G]'
eroi_disp_pv









    



PV [F,G]






    Out[25]:





$$\left [ 7.56681222707424, \quad 20.0323699421965\right ]$$



In [26]:

    
print 'Wind [H,I]'
eroi_disp_wind









    



Wind [H,I]






    Out[26]:





$$\left [ 9.03659317027671, \quad 29.2781453424333\right ]$$

Construct test cases: Fossils



In [28]:

    
eroi_gas=[9.6, 38.5] #EROEIe for case [A, B] in Table 3 of the paper, gas plant
eroi_ccs_gas=format_eq([eq6.subs([[EROI,i],[f_op,0.147]]) for i in eroi_gas])
eroi_coal=[8.9, 22.6] #EROEIe for case [D, E] in Table 3 of the paper, coal plant
eroi_ccs_coal=format_eq([eq6.subs([[EROI,i],[f_op,0.219]]) for i in eroi_coal])



In [29]:

    
print 'gas [A,B]'
eroi_ccs_gas









    



gas [A,B]






    Out[29]:





$$\left [ 3.9814200398142, \quad 5.78121480591636\right ]$$



In [30]:

    
print 'coal [D,E]'
eroi_ccs_coal









    



coal [D,E]






    Out[30]:





$$\left [ 3.01786985860093, \quad 3.79870239015699\right ]$$

Figure 3

Produce contour lines for multiple $f_{op}$ and $\phi$ levels and place the example plants on the contours. Calculate the $EROI_{CCS}$ and $EROI_{disp}$, respectively, for multiple base $EROI$ values. Likewise we obtain a plot for comparison between dispatchable (i.e. with storage) EROI for renewables and CCS EROI for fossil fuel based power plants.



In [31]:

    
plt.figure(figsize=(5,4))
n=51
erois=range(1,n)

### Renewables contours
colors=['#addd8e',
        '#78c679',
        '#31a354',
        '#006837']
plt.plot(0,0,lw=0,label='${\phi}$ (%)')
#plot EROI_disp contour lines for various EROI values and phi levels
#for eta and S assume a fully battery-based storage (worst case), i.e. the eta_pv and S_pv
pv=1.0 #(pv share in hypothetical renewable technologies)
phis=[0.1,0.2,0.3,0.4] #assumed levels for phi
for r in range(len(phis)):
    phir=phis[r]
    eroi_disps=format_eq([eq2.subs([[EROI,i],[phi,phir],[eta,pv*eta_pv+(1-pv)*eta_wind],[ESOI,pv*S_pv+(1-pv)*S_wind]]) for i in erois])
    plt.plot(erois,eroi_disps,c=colors[r],label=int(phir*100))

plt.xlabel('EROI')
a=plt.gca()
plt.text(-0.07, 0.65, 'EROI$_{disp}$',
        horizontalalignment='right',
        verticalalignment='bottom',rotation=90,#fontproperties=prop,
        size=11,color='#006837',transform=plt.gca().transAxes)
plt.text(-0.07, 0.27, 'EROI$_{CCS}$',
        horizontalalignment='right',
        verticalalignment='bottom',rotation=90,#fontproperties=prop,
        size=11,color='#bd0036',transform=plt.gca().transAxes)

### Fossils contours
colors=['#fed976',
        '#feb24c',
        '#f03b20',
        '#bd0036']
plt.plot(0,0,lw=0,label='$f_{op}$ (%)')
fops=[0.05,0.1,0.15,0.2]  #assumed levels for fop
for r in range(len(fops)):
    fopr=fops[r]
    eroi_ccss=format_eq([eq6.subs([[EROI,i],[f_op,fopr]]) for i in erois])
    plt.plot(erois,eroi_ccss,c=colors[r],label=int(fopr*100))

    
### Table 3 cases
cf="#f03b20"
cr="#31a354"
cg="#aaaaaa"
s=30 #marker size
plt.scatter(eroi_gas,eroi_ccs_gas,s,color=cf,marker='o')
plt.scatter(eroi_coal,eroi_ccs_coal,s,color=cf,marker='^')
plt.scatter(eroi_pv,eroi_disp_pv,s,color=cr,marker='d')
plt.scatter(eroi_wind,eroi_disp_wind,s,color=cr,marker='s')
plt.xlim(0,50)
plt.ylim(0,40)
#[F, G]
plt.text(eroi_pv[1], eroi_disp_pv[1]-1.5, 'PV +\nBattery',
        horizontalalignment='center',
        verticalalignment='top',
        size=9,color=cr)
#[H, I]
plt.text(eroi_wind[1], eroi_disp_wind[1]+1.5, 'Wind + Battery +\n PHS + CAES',
        horizontalalignment='center',
        verticalalignment='bottom',
        size=9,color=cr)
#[A, B]
plt.text(eroi_gas[1], eroi_ccs_gas[1]-2.2, 'Gas + CCS',
        horizontalalignment='center',
        verticalalignment='top',
        size=9,color=cf)
#[D, E]
plt.text(eroi_coal[1], eroi_ccs_coal[1]-1, 'Coal + CCS',
        horizontalalignment='center',
        verticalalignment='top',
        size=9,color=cf)

#Annotate cases
ap=0.7
plt.text(eroi_coal[0]+2, eroi_ccs_coal[0]-0.5, 'D',
        horizontalalignment='center',
        verticalalignment='top',
        size=10,color=cf,
        bbox=dict(facecolor='white', edgecolor=cf,  alpha=ap, boxstyle='round'))
plt.text(eroi_gas[0]+2, eroi_ccs_gas[0]+2, 'B',
        horizontalalignment='center',
        verticalalignment='top',
        size=10,color=cf,
        bbox=dict(facecolor='white', edgecolor=cf,  alpha=ap, boxstyle='round'))
plt.text(eroi_coal[1], eroi_ccs_coal[1]+3.3, 'E',
        horizontalalignment='center',
        verticalalignment='top',
        size=10,color=cf,
        bbox=dict(facecolor='white', edgecolor=cf,  alpha=ap, boxstyle='round'))
plt.text(eroi_gas[1], eroi_ccs_gas[1]+3.3, 'A',
        horizontalalignment='center',
        verticalalignment='top',
        size=10,color=cf,
        bbox=dict(facecolor='white', edgecolor=cf, alpha=ap, boxstyle='round'))

plt.text(eroi_pv[1], eroi_disp_pv[1]+3.3, 'G',
        horizontalalignment='center',
        verticalalignment='top',
        size=10,color=cr,
        bbox=dict(facecolor='white', edgecolor=cr, alpha=ap, boxstyle='round'))
plt.text(eroi_wind[0]+2.6, eroi_disp_wind[0], 'F',
        horizontalalignment='center',
        verticalalignment='top',
        size=10,color=cr,
        bbox=dict(facecolor='white', edgecolor=cr,  alpha=ap, boxstyle='round'))
plt.text(eroi_pv[0], eroi_disp_pv[0]+5.3, 'H',
        horizontalalignment='center',
        verticalalignment='top',
        size=10,color=cr,
        bbox=dict(facecolor='white', edgecolor=cr,  alpha=ap, boxstyle='round'))
plt.text(eroi_wind[1], eroi_disp_wind[1]-2, 'I',
        horizontalalignment='center',
        verticalalignment='top',
        size=10,color=cr,
        bbox=dict(facecolor='white', edgecolor=cr,  alpha=ap, boxstyle='round'))

plt.legend(framealpha=0,loc=2,fontsize=9)
plt.savefig('fig3.png',bbox_inches = 'tight', facecolor='w', pad_inches = 0.1, dpi=200)
plt.show()

Figure 4

For this figure we are considering various levels of powerplant efficiency, and its effect on the EROI , under various capture ratio assumptions. This is described in Eq. 1 and Eq. 4 in the paper.

Define parameters

Electricity penalty (b) curve from Sanpasternich (2009) Figure 8 [kWh/kg CO2]



In [33]:

    
#CCS_op #operating CCS penalty
#CR #capture ratio
var('CCS_op CR');
CCS_op_eq=Eq(CCS_op,25.957**6 - 85.031*CR**5 + \
114.5*CR**4 - 80.385*CR**3 + \
31.47*CR**2 - 6.7725*CR  + 1.1137)
CCS_op_eq









    Out[33]:





$$CCS_{op} = - 85.031 CR^{5} + 114.5 CR^{4} - 80.385 CR^{3} + 31.47 CR^{2} - 6.7725 CR + 305863028.421403$$

Operational cost of coal fuel cycle over plant lifetime output [kWh]



In [34]:

    
#PP_op_L #operational cost of coal fuel cycle over plant lifetime
#PP_op #operational cost of coal fuel cycle
#PP_CF #powerplant capacity factor
#PP_L #powerplant lifetime
#eta #powerplant efficiency
var('PP_op_L PP_op PP_CF eta PP_L');
PP_op_L_eq=Eq(PP_op_L,8760*PP_L*PP_CF*PP_op/eta)
PP_op_L_eq









    Out[34]:





$$PP_{op L} = \frac{8760 PP_{CF}}{\eta} PP_{L} PP_{op}$$

Carbon intensity of electricity [kg CO2/kWh]



In [35]:

    
#C_CO2 #carbon dioxide content of coal [kg/MJ]
#Elec_CO2 #carbon intensity of electricity [kg CO2/kWh]
var('Elec_CO2 C_CO2');
Elec_CO2_eq=Eq(Elec_CO2,C_CO2*3.6/eta)
Elec_CO2_eq









    Out[35]:





$$Elec_{CO2} = \frac{3.6 C_{CO2}}{\eta}$$

Energy penalty of CCS [dmnl]



In [36]:

    
#b #energy penalty of CCS
var('b');
b_eq=Eq(b,CCS_op*Elec_CO2*CR)
b_eq









    Out[36]:





$$b = CCS_{op} CR Elec_{CO2}$$

Emissions



In [37]:

    
#E #emissions
var('E');
E_eq=Eq(E,CCS_op*(1-CR))
E_eq









    Out[37]:





$$E = CCS_{op} \left(- CR + 1\right)$$

Energy cost of constructing CCS plant



In [40]:

    
#CCS_cons_energy #energy cost of constructing CCS plant
#f_cap #energy cost share of constructing CCS plant
#PP_cons_energy #energy cost of power plant construction [MJ/MW], does not include energy embodied in materials
var('CCS_cons_energy PP_cons_energy f_cap');
CCS_cons_energy_eq=Eq(CCS_cons_energy,f_cap*PP_cons_energy)
CCS_cons_energy_eq









    Out[40]:





$$CCS_{cons energy} = PP_{cons energy} f_{cap}$$

EROIs

Energy cost of constructing CCS plant



In [45]:

    
#EROI_1 #EROI where CCS construction energy is a ratio of power plant construction energy [MJ_e/MW/(MJ_p/MW) = dmnl]
#CCS_L #CCS plant lifetime
var('EROI_1 CCS_L');
EROI_1_eq=Eq(EROI_1,8760*PP_L*PP_CF*(1-b)*3.6/\
                    (PP_cons_energy + PP_op_L + CCS_cons_energy*PP_L/CCS_L))
EROI_1_eq









    Out[45]:





$$EROI_{1} = \frac{31536.0 PP_{CF} PP_{L} \left(- b + 1\right)}{PP_{cons energy} + PP_{op L} + \frac{CCS_{cons energy}}{CCS_{L}} PP_{L}}$$



In [42]:

    
#EROI_1_adj #makes adjustment for electrical output [MJ_p/MJ_p = dmnl]
var('EROI_1_adj');
EROI_1_adj_eq=Eq(EROI_1_adj,EROI_1/0.3)
EROI_1_adj_eq









    Out[42]:





$$EROI_{1 adj} = 3.33333333333333 EROI_{1}$$

Energy cost of operating CCS plant



In [46]:

    
#EROI_2 #where CCS construction energy is a ratio of CCS operation energy [MJ_e/MJ_p = dmnl]
var('EROI_2');
EROI_2_eq=Eq(EROI_2,8760*PP_L*PP_CF*(1-b)/\
                    (PP_cons_energy + PP_op_L + f_op*CCS_op + CCS_cons_energy*PP_L/CCS_L))
EROI_2_eq









    Out[46]:





$$EROI_{2} = \frac{8760 PP_{CF} PP_{L} \left(- b + 1\right)}{CCS_{op} f_{op} + PP_{cons energy} + PP_{op L} + \frac{CCS_{cons energy}}{CCS_{L}} PP_{L}}$$



In [47]:

    
#EROI_2_adj #makes adjustment for electrical output [MJ_p/MJ_p = dmnl]
var('EROI_2_adj');
EROI_2_adj_eq=Eq(EROI_2_adj,EROI_1/0.3)
EROI_2_adj_eq









    Out[47]:





$$EROI_{2 adj} = 3.33333333333333 EROI_{1}$$

Translation to quantities as per defined in the paper

Eq. 5



In [48]:

    
#R #E_op/E_cap
var('R');
eq5=Eq(EROI_CCS,(R+1)*EROI/(R+1+f_cap+f_op*(R+1)*EROI))
eq5









    Out[48]:





$$EROI_{CCS} = \frac{EROI \left(R + 1\right)}{EROI f_{op} \left(R + 1\right) + R + f_{cap} + 1}$$



In [51]:

    
var('E_op E_cap');
R_eq=Eq(R,E_op/E_cap)
R_eq









    Out[51]:





$$R = \frac{E_{op}}{E_{cap}}$$

Eq. 4



In [52]:

    
var('E_CCS E_out E_red');
eq4=Eq(EROI_CCS,E_out/((1+f_cap)*E_cap+E_op+f_op*E_out))
eq4









    Out[52]:





$$EROI_{CCS} = \frac{E_{out}}{E_{cap} \left(f_{cap} + 1\right) + E_{op} + E_{out} f_{op}}$$



In [53]:

    
f_op_eq=Eq(f_op,E_red/E_out)
f_op_eq









    Out[53]:





$$f_{op} = \frac{E_{red}}{E_{out}}$$



In [103]:

    
f_cap_eq=Eq(f_cap,E_CCS/E_cap)
f_cap_eq









    Out[103]:





$$f_{cap} = \frac{E_{CCS}}{E_{cap}}$$

Eq. 3



In [92]:

    
eq3=Eq(EROI_CCS,E_out/(E_cap+E_op+E_red+E_CCS))
eq3









    Out[92]:





$$EROI_{CCS} = \frac{E_{out}}{E_{CCS} + E_{cap} + E_{op} + E_{red}}$$

Eq. 1



In [93]:

    
var('E_in');
eq1=Eq(EROI,E_out/E_in)
eq1









    Out[93]:





$$EROI = \frac{E_{out}}{E_{in}}$$



In [94]:

    
E_in_eq=Eq(E_in,E_cap+E_op)
E_in_eq









    Out[94]:





$$E_{in} = E_{cap} + E_{op}$$



In [95]:

    
solve(E_in_eq,E_in)[0]









    Out[95]:





$$E_{cap} + E_{op}$$



In [96]:

    
eq1.subs([[E_in,solve(E_in_eq,E_in)[0]]])









    Out[96]:





$$EROI = \frac{E_{out}}{E_{cap} + E_{op}}$$

Verification

Bactkrack to Eq. 3 from Eq. 5



In [97]:

    
eq5









    Out[97]:





$$EROI_{CCS} = \frac{EROI \left(R + 1\right)}{EROI f_{op} \left(R + 1\right) + R + f_{cap} + 1}$$

Get back Eq. 4 from Eq. 5



In [98]:

    
eq5.subs([[R,solve(R_eq,R)[0]],\
          [EROI,solve(eq1.subs([[E_in,solve(E_in_eq,E_in)[0]]]),EROI)[0]]]).simplify()









    Out[98]:





$$EROI_{CCS} = \frac{E_{out}}{E_{cap} \left(f_{cap} + 1\right) + E_{op} + E_{out} f_{op}}$$



In [99]:

    
eq4









    Out[99]:





$$EROI_{CCS} = \frac{E_{out}}{E_{cap} \left(f_{cap} + 1\right) + E_{op} + E_{out} f_{op}}$$

Get back Eq. 3 from Eq. 4



In [104]:

    
eq4.subs([[f_op,solve(f_op_eq,f_op)[0]],\
         [f_cap,solve(f_cap_eq,f_cap)[0]]]).simplify()









    Out[104]:





$$EROI_{CCS} = \frac{E_{out}}{E_{CCS} + E_{cap} + E_{op} + E_{red}}$$



In [105]:

    
eq3









    Out[105]:





$$EROI_{CCS} = \frac{E_{out}}{E_{CCS} + E_{cap} + E_{op} + E_{red}}$$

Translate

Translate example notation into paper equation notations



In [108]:

    
eq3









    Out[108]:





$$EROI_{CCS} = \frac{E_{out}}{E_{CCS} + E_{cap} + E_{op} + E_{red}}$$



In [109]:

    
E_out_eq=Eq(E_out,8760*PP_L*PP_CF*(1-b))
E_out_eq









    Out[109]:





$$E_{out} = 8760 PP_{CF} PP_{L} \left(- b + 1\right)$$



In [115]:

    
E_out_eq.subs([[b,solve(b_eq,b)[0]]])\
        .subs([[CCS_op,solve(CCS_op_eq,CCS_op)[0]]])\
        .subs([[Elec_CO2,solve(Elec_CO2_eq,Elec_CO2)[0]]])









    Out[115]:





$$E_{out} = 8760 PP_{CF} PP_{L} \left(- \frac{3.6 CR}{\eta} C_{CO2} \left(- 85.031 CR^{5} + 114.5 CR^{4} - 80.385 CR^{3} + 31.47 CR^{2} - 6.7725 CR + 305863028.421403\right) + 1\right)$$

... ?

Construct from bottom up



In [125]:

    
eq3









    Out[125]:





$$EROI_{CCS} = \frac{E_{out}}{E_{CCS} + E_{cap} + E_{op} + E_{red}}$$