Comparative net energy analysis of renewable electricity and carbon capture and storage

Sgouridis, Carbajales-Dale, Csala, Chiesa, Bardi

Use this Jupyter workbook to reproduce all results presented in the paper.

This is document has been created using Jupyter Notebook in the Anaconda distribution and it can be edited and run in active mode by clicking download in top right corner of this page. The code is partitioned into sections, called cells. When you are using this workbook in active mode, double-click on a cell to edit it and then run using Ctrl + Enter. Hitting Shift + Enter runs the code and steps into the next cell, while Alt + Enter runs the code and adds a new, empty cell. If you are running this notebook on a presonal computer, you will need a machine with at least 1GB of memory (2GB recommended) and a processor of 1GHz.

Import dependencies

In [1]:
import pandas as pd, matplotlib.pyplot as plt, numpy as np
from matplotlib.colors import BoundaryNorm
from matplotlib.ticker import MaxNLocator
import matplotlib.patches as patches
import matplotlib.lines as lines
%matplotlib inline
import sys, string
from sympy import *
init_printing() 
print('Running on Python',sys.version)
Running on Python 3.6.6 |Anaconda, Inc.| (default, Jun 28 2018, 11:27:44) [MSC v.1900 64 bit (AMD64)]

Fossil fuels part

Equations

This section reproduces all the equations presented in the paper (supplementary material). This is just for easier visual guidance and for quick checking.

In [2]:
var('EROEI E_out E_in');
eq_S1=Eq(EROEI,E_out/E_in)
eq_S1
Out[2]:
$$EROEI = \frac{E_{out}}{E_{in}}$$
In [3]:
var('E_cap E_f E_OM');
eq_S1i=Eq(EROEI,E_out/(E_cap+E_OM+E_f))
eq_S1i
Out[3]:
$$EROEI = \frac{E_{out}}{E_{OM} + E_{cap} + E_{f}}$$
In [4]:
var('L s_OM');
eq_S2i=Eq(E_in,E_cap*(1+L*s_OM)+E_f)
eq_S2i
Out[4]:
$$E_{in} = E_{cap} \left(L s_{OM} + 1\right) + E_{f}$$
In [5]:
var('P L cf');
eq_S2ii=Eq(E_out,P*cf*L)
eq_S2ii
Out[5]:
$$E_{out} = L P cf$$
In [6]:
var('E_el');
eq_S2iii=Eq(E_out,E_el)
eq_S2iii
Out[6]:
$$E_{out} = E_{el}$$
In [7]:
e=Symbol('eta')
var('E_th');
eq_S2iv=Eq(E_th,P*cf*L/e)
eq_S2iv
Out[7]:
$$E_{th} = \frac{L P cf}{\eta}$$
In [8]:
var('EROEI_th');
eq_S3i=Eq(E_f,E_th/EROEI_th)
eq_S3i
Out[8]:
$$E_{f} = \frac{E_{th}}{EROEI_{th}}$$
In [9]:
#define pretty variable substituion function
def sub(var,eq):
    return [var,solve(eq,var)[0]]
#create multiple substituion method by chaining
def multisub(eq0,msub):
    for i in range(len(msub)//2):
        eq0=eq0.subs([sub(msub[i*2],msub[i*2+1])])
    return eq0
In [10]:
var('EROEI_el');
eq_S3=Eq(EROEI_el,multisub(eq_S1,[E_out,eq_S2ii,E_in,eq_S2i,E_f,eq_S3i,E_th,eq_S2iv]).rhs)
eq_S3
Out[10]:
$$EROEI_{el} = \frac{L P cf}{E_{cap} L s_{OM} + E_{cap} + \frac{L P cf}{EROEI_{th} \eta}}$$
In [11]:
var('CR f_op EROEI_CCS f_cap');
eq_S4=Eq(EROEI_CCS,E_el*(1-f_op(CR))/(E_cap*(1+f_cap)*(1+L*s_OM)+E_f))
eq_S4
Out[11]:
$$EROEI_{CCS} = \frac{E_{el} \left(- \operatorname{f_{op}}{\left (CR \right )} + 1\right)}{E_{cap} \left(f_{cap} + 1\right) \left(L s_{OM} + 1\right) + E_{f}}$$
In [12]:
var('CR f_op EROEI_CCS f_cap');
eq_S5=Eq(EROEI_CCS,(1-f_op(CR))*(E_cap*(1+L*s_OM)+E_f)*EROEI_el/(E_cap*(1+f_cap)*(1+L*s_OM)+E_f))
eq_S5
Out[12]:
$$EROEI_{CCS} = \frac{EROEI_{el} \left(E_{cap} \left(L s_{OM} + 1\right) + E_{f}\right) \left(- \operatorname{f_{op}}{\left (CR \right )} + 1\right)}{E_{cap} \left(f_{cap} + 1\right) \left(L s_{OM} + 1\right) + E_{f}}$$
In [13]:
var('R');
eq_S6i=Eq(R,E_f/(E_cap*(1+L*s_OM)))
eq_S6i
Out[13]:
$$R = \frac{E_{f}}{E_{cap} \left(L s_{OM} + 1\right)}$$
In [14]:
eq_S6=Eq(EROEI_CCS,(1-f_op(CR))*(R+1)*EROEI_el/(R+1+f_cap))
eq_S6
Out[14]:
$$EROEI_{CCS} = \frac{EROEI_{el} \left(R + 1\right) \left(- \operatorname{f_{op}}{\left (CR \right )} + 1\right)}{R + f_{cap} + 1}$$

We have made a major simplifaction step from Eq. S5 to reach Eq. S6. Here we prove by symbolic evaluation of the subtituted expressions that we were correct.

In [15]:
eq_S5.simplify()
Out[15]:
$$EROEI_{CCS} = - \frac{EROEI_{el} \left(E_{cap} \left(L s_{OM} + 1\right) + E_{f}\right) \left(\operatorname{f_{op}}{\left (CR \right )} - 1\right)}{E_{cap} \left(f_{cap} + 1\right) \left(L s_{OM} + 1\right) + E_{f}}$$
In [16]:
eq_S6.subs(R,eq_S6i.rhs).simplify()
Out[16]:
$$EROEI_{CCS} = - \frac{EROEI_{el} \left(E_{cap} \left(L s_{OM} + 1\right) + E_{f}\right) \left(\operatorname{f_{op}}{\left (CR \right )} - 1\right)}{E_{cap} \left(f_{cap} + 1\right) \left(L s_{OM} + 1\right) + E_{f}}$$
In [17]:
#eq_S6 is correct!
eq_S6.subs(R,eq_S6i.rhs).simplify()==eq_S5.simplify()
Out[17]:
True

Simulation

In this part we show you how to reproduce the graphs presented in the paper, which form the basis of our results discussion and conclusions. There is a similar section for the Renewables part. First, we illustrate a representative range for CCS operational parameters for a number of fossil fuel plants. We then take note of the min and max ranges and plot the values alongside those obtained from the Renewables part. The calculations are based on actual power plant data, stored in the Excel file plants.xlsx, found in the project folder and presented in Table S1 in the supplementary material of the paper.

Using the formula in Eq.4 and the plant values from the Table S1 (reproduced from plants.xlsx) calculate the $EROEI_{el}$ for all capacity factors $CF$s between 50 and 90. We then use these values as a proxy to create an a adjustment factor for the $EROEI_{el}$ due to the indirect effect of the change of the capacity factor. The same is true for the operational versus capital energy investment ratio, $R$. We therefore created capacity-factor adjusted values for $R$ as well.

In order to manipulate the Excel file on the fly, we use a python module called xlwings. This will open an Excel application in the background. Be sure to either manually close it after you're done or run the cell just above the Auxiliary plots section.

In [18]:
import xlwings as xw #will open a MS Excel application, let it run in the background
wb = xw.Book('plants.xlsx')

We define a function eroei_r_getter which takes as parameters the capacity factor $CF$ and the $EROEI_{fuel}$ for the two studies fossil fuels, Coal and Natural Gas. Then, this function sets the correct $CF$ and $EROEI_{fuel}$ values in the Excel workbook, which then recalculates the resulting $EROEI_{el}$ and $R$ values. Since we have 3 IGCC Coal, 2 Pulverized Coal and 1 NGCC Gas plant, we average over $EROEI_{el}$ and $R$ values within each group. The resulting pandas dataframe is returned by the function, which now contains 3 columns ['IGCC','PC','NGCC'] and 2 rows ['EROEI_el','R'].

In [19]:
def eroei_r_getter(CF=100,EROEI_fuel_coal=58,EROEI_fuel_gas=87):
    xw.Range('L7').value=CF/100.0 #hard-coded Excel workbook cell ID for capacity factor
    xw.Range('M8').value=EROEI_fuel_coal #hard-coded Excel workbook cell ID for EROEI_fuel for Coal
    xw.Range('O8').value=EROEI_fuel_gas #hard-coded Excel workbook cell ID for EROEI_fuel for Gas
    wb.save() #write new values
    df=pd.read_excel('plants.xlsx',usecols=range(10,24)).loc[[37,38]]
    df=df[df.columns[1::2]]
    df.columns=['IGCC 1','IGCC 3','IGCC 5','PC 9','PC 11','NGCC 13']
    df.index=['EROEI_el','R']
    df['IGCC']=(df['IGCC 1']+df['IGCC 3']+df['IGCC 5'])/3.0
    df['PC']=(df['PC 9']+df['PC 11'])/2.0
    df['NGCC']=df['NGCC 13']
    return df[['IGCC','PC','NGCC']]
In [20]:
eroei_r_getter(95,58,87) #example to get the EROEI_el and R for capacity factor 95% and EROEI_fuels of 58 for Coal and 87 for Gas
Out[20]:
IGCC PC NGCC
EROEI_el 12.0177 14.0583 31.9783
R 1.5745 2.06751 2.91996

Create reasonable operating range for capacity factor $CF$ (%).

In [21]:
CFs=[i*5 for i in range (10,21)]
CFs
Out[21]:
$$\left [ 50, \quad 55, \quad 60, \quad 65, \quad 70, \quad 75, \quad 80, \quad 85, \quad 90, \quad 95, \quad 100\right ]$$

This cell is just an illustration of the pseudo-linear relationship between $CF$ and $EROEI_{el}$ (top row) and $R$ (bottom row), respectively. Note that the relationship is also nonlinear with regards to $EROEI_{fuel}$ (left vs. right column).

In [22]:
EROEI_fuel=[[58,87],[23,20]]
dfe=[]
dfr=[]
for j in range(len(EROEI_fuel)):
    EROEI_fuel_coal=EROEI_fuel[j][0]
    EROEI_fuel_gas=EROEI_fuel[j][1]
    dfs1=[]
    dfs2=[]
    for i in range(len(CFs)):
        df=eroei_r_getter(CFs[i],EROEI_fuel_coal,EROEI_fuel_gas)
        df1=df.loc['EROEI_el']
        df2=df.loc['R']
        df1.name=CFs[i]
        df2.name=CFs[i]
        dfs1.append(df1)
        dfs2.append(df2)
    dfe.append(pd.concat(dfs1,axis=1))
    dfr.append(pd.concat(dfs2,axis=1))

Figure S3. presents two $EROEI_{fuel}$ values for Coal and Gas (columns), respectively, and the dependence of $EROEI_{el}$ (top row) and the Fuel to capital ratio $R$ (bottom row) on the Capacity Factor $CF$.

In [26]:
#Plots init
fig,axes=plt.subplots(2,2,figsize=(8,6),subplot_kw=dict(facecolor='#EEEEEE',axisbelow=True))
for axi in axes:
    for ax in axi:
        ax.grid(color='white', linestyle='solid')                       
plt.subplots_adjust(wspace=0.15)
plt.subplots_adjust(hspace=0.1)

for j in range(len(EROEI_fuel)):
    EROEI_fuel_coal=EROEI_fuel[j][0]
    EROEI_fuel_gas=EROEI_fuel[j][1]
    df=dfe[j]
    ax=axes[0][j]
    df.T.plot(ax=ax,grid=True,legend=False)
    ax.set_title('$EROEI_{fuel}$ Coal$='+str(EROEI_fuel_coal)+'$\n$EROEI_{fuel}$ Gas$='+str(EROEI_fuel_gas)+'$') 
    df=dfr[j]
    ax=axes[1][j]
    df.T.plot(ax=ax,grid=True,legend=False)
    for axi in range(len(axes[j])):
        ax=axes[j][axi]
        ax.text(0.04, 0.95, string.ascii_lowercase[axi*len(axes)+j],
            horizontalalignment='left',
            verticalalignment='top',
            transform=ax.transAxes,size=13,alpha=0.7,
                          bbox=dict(boxstyle="round", fc="w", alpha=0.7))
        

axes[1][1].set_xlabel('Capacity Factor $CF$  (%)',fontsize=13,x=0)
axes[0][0].set_ylabel('$EROEI_{el}$',fontsize=13)
axes[1][0].set_ylabel('Fuel to capital ratio $R$',fontsize=13)
#axes[0][1].set_yticklabels('')   
#axes[1][1].set_yticklabels('')   
axes[0][0].set_xticklabels('')   
axes[0][1].set_xticklabels('')  
axes[0][0].legend(framealpha=0)

#plt.suptitle('Normalized Relationship of $EROEI_{el}$ and $R$ to Capacity Factor $CF$ ',fontsize=13,y=1.01)
plt.savefig('plot/png/figS3.png',dpi=300,bbox_inches = 'tight', pad_inches = 0.1)
plt.savefig('plot/eps/figS3.eps',format='eps',bbox_inches = 'tight', pad_inches = 0.1)
plt.savefig('plot/pdf/figS3.pdf',format='pdf',bbox_inches = 'tight', pad_inches = 0.1)
plt.show()

We need to create a energy penalty conversion ratio for the Capture Ratio $CR$. We use the Energy penalty figure Fig. S1 from the supplementary material. We then use this value to adjust the value of the operational energy penalty $f_{op}$. All this means that a higher Capture Ratio $CR$ lowers $EROEI_{CCS}$ more - and we would like to create a measure about how much. We use the middle solid black line, linear approximation from the plot.

CR energy penalty
60 7.8
90 11.8

We will only test Capture Ratios $CR$ of $60$ and $90$.

In [58]:
CR=[90,60]
CR_ratio=[1, 7.8/11.8]
CR_ratio
Out[58]:
$$\left [ 1, \quad 0.6610169491525423\right ]$$

Create a simple fuel name convereter function.

In [59]:
def fuel_name_converter(fuel):
    if fuel=='NGCC\nGas': return 'NGCC'
    elif fuel=='IGCC\nCoal': return 'IGCC'
    elif fuel=='Pulverized\nCoal': return 'PC'
    else: return None

Define $EROEI_{CCS}$ calculator function eroei_ccs_cycler, including minimums and maximums over $f_{op}$ and $f_{cap}$. This function takes as input:

  • an array names, which holds the fossil fuel type names
  • their $EROEI_{fuel}$s, same length, in the same order
  • $f$ extent envelopes t0s in the format [min $f_{op}$,max $f_{op}$, min $f_{cap}$,max $f_{cap}$], with this array repreated for each name
  • Capacity Factor $CF$ array - of any length, percentage values between $50$ and $100$, increments of $5$.
  • Cpature Ratio $CR$ array - of any length - but only defined for percentage values $60$ and $90$ need to define other levels using Fig. S1 if using differenet values.
  • Ranges of any length for percentage values $f_{op}$ and $f_{cap}$. Values between $0$ and $200$, any increment.

The function will iterate over all inputs and return 4 multidimensional variables. First, uisng the eroei_r_getter, the function will calculate the corresponding $EROEI_{el}$ and $R$ to the $CF$ and $EROEI_{fuel}$ input. Then it will calculate $EROEI_{CCS}$ using Eq. S6, for all variable ranges. The t0x inner loop keeps track of the minimums and maximums. Finally, it returns the calculated 2-dimensional arrays $EROEI_{el}$ and $R$ in eroei_el_cf and R_cf, and multidimensional arrays eroei_ccs, eroei_ccs_min, eroei_ccs_max, holdig the actual values, the minimums and the maximums for $EROEI_{CCS}$.

In [60]:
def eroei_ccs_cycler(names, eroei_fuels, t0s, CF, CR_ratio, fop, fcap, verbose=False):
    eroei_ccs=np.zeros([len(names),len(CF),len(CR_ratio),len(fop),len(fcap)])
    eroei_ccs_min=np.zeros([len(names),len(CF),len(CR_ratio)])+1000
    eroei_ccs_max=np.zeros([len(names),len(CF),len(CR_ratio)])
    eroei_el_cf=np.zeros([len(names),len(CF)])
    R_cf=np.zeros([len(names),len(CF)])
    for i in range(len(names)):
        if verbose: print('Calculating '+names[i].replace('\n',' ')+'...')
        for j in range(len(CF)):
            eroei_r=eroei_r_getter(CF[j],eroei_fuels[i],eroei_fuels[i])[fuel_name_converter(names[i])]
            #convert fuel EROEI to EROEI_el using the CF converstion table
            eroei_el_cf[i][j]=eroei_r.loc['EROEI_el']
            #calculate R using the CF converstion table
            R_cf[i][j]=eroei_r.loc['R']
            for r in range(len(CR_ratio)):
                for k in range(len(fop)):
                    for l in range(len(fcap)):
                        #calculate eoris over maps using eq_S6
                        eroei_ccs[i][j][r][l,k]=(1-fop[k]*CR_ratio[r]/100.0)*\
                            ((R_cf[i][j]+1)/(R_cf[i][j]+1+fcap[l]/100.0))*\
                            eroei_el_cf[i][j]

                        #calculate min-maxes
                        for t0x in range(len(t0s[i])//2):
                            for t0y in range(len(t0s[i])//2,len(t0s[i])):
                                x=t0s[i][t0x]
                                y=t0s[i][t0y]
                                if ((abs(fop[k]-x)<0.5)and(abs(fcap[l]-y)<1.5)):
                                    eroei_ccs_max[i][j][r]=max(eroei_ccs_max[i][j][r],eroei_ccs[i][j][r][l,k])
                                    eroei_ccs_min[i][j][r]=min(eroei_ccs_min[i][j][r],eroei_ccs[i][j][r][l,k])
                if verbose:
                    print('CF='+str(CF[j])+', CR='+str(CR[r])+\
                      ', EROEI_fuel='+str(eroei_fuels[i])+\
                      ', EROEI_el_CF_'+str(CF[j])+'%='+str(np.round(eroei_el_cf[i][j],2))+\
                      ', R_CF_'+str(CF[j])+'%='+str(np.round(R_cf[i][j],2))+\
                      ', EROEI_CCS - min='+str(np.round(eroei_ccs_min[i][j][r],2))+\
                      ', max='+str(np.round(eroei_ccs_max[i][j][r],2)))
    return eroei_el_cf, R_cf, eroei_ccs, eroei_ccs_min, eroei_ccs_max

Having set up the eroei_ccs_cycler function, we define our desired simulation ranges and run it. We also display min and max values for the $EROEI_{CCS}$. These are the bottom left and top right corners of the rectangles in Figure 2.

In [61]:
#Simulation options
fop=np.linspace(5,40,100) #in %
fcap=np.linspace(10,120,100) #in %
CF=[55,85] #in %
title=u'$EROEI_{CCS}$ comparison for Coal and Gas'
#eroei_els=[30,13,13] #deprecated, varies by name
#eroei_els=[32.4,14.3,12.25] #deprecated, 100% CF value from the CF converstion table
eroei_fuels=[87,58,58]
names=['NGCC\nGas','Pulverized\nCoal','IGCC\nCoal'] 
names_=['NGCC Gas','Pulverized Coal','IGCC Coal'] 
t0s=[[12,20,58,110],[17,31,45,95],[15,28,23,40]]  # [fop min, fop max, fcap min fcap max] by each name, in %

eroei_el_cf, R_cf, eroei_ccs, eroei_ccs_min, eroei_ccs_max=eroei_ccs_cycler(names, eroei_fuels, t0s, CF, CR_ratio, fop, fcap, verbose=True)
Calculating NGCC Gas...
CF=55, CR=90, EROEI_fuel=87, EROEI_el_CF_55%=26.97, R_CF_55%=1.69, EROEI_CCS - min=15.23, max=19.67
CF=55, CR=60, EROEI_fuel=87, EROEI_el_CF_55%=26.97, R_CF_55%=1.69, EROEI_CCS - min=16.54, max=20.56
CF=85, CR=90, EROEI_fuel=87, EROEI_el_CF_85%=31.05, R_CF_85%=2.61, EROEI_CCS - min=18.95, max=23.69
CF=85, CR=60, EROEI_fuel=87, EROEI_el_CF_85%=31.05, R_CF_85%=2.61, EROEI_CCS - min=20.57, max=24.76
Calculating Pulverized Coal...
CF=55, CR=90, EROEI_fuel=58, EROEI_el_CF_55%=11.36, R_CF_55%=1.2, EROEI_CCS - min=5.45, max=7.87
CF=55, CR=60, EROEI_fuel=58, EROEI_el_CF_55%=11.36, R_CF_55%=1.2, EROEI_CCS - min=6.29, max=8.41
CF=85, CR=90, EROEI_fuel=58, EROEI_el_CF_85%=13.54, R_CF_85%=1.85, EROEI_CCS - min=6.98, max=9.76
CF=85, CR=60, EROEI_fuel=58, EROEI_el_CF_85%=13.54, R_CF_85%=1.85, EROEI_CCS - min=8.05, max=10.42
Calculating IGCC Coal...
CF=55, CR=90, EROEI_fuel=58, EROEI_el_CF_55%=9.37, R_CF_55%=0.91, EROEI_CCS - min=5.52, max=7.17
CF=55, CR=60, EROEI_fuel=58, EROEI_el_CF_55%=9.37, R_CF_55%=0.91, EROEI_CCS - min=6.26, max=7.58
CF=85, CR=90, EROEI_fuel=58, EROEI_el_CF_85%=11.49, R_CF_85%=1.41, EROEI_CCS - min=7.03, max=8.99
CF=85, CR=60, EROEI_fuel=58, EROEI_el_CF_85%=11.49, R_CF_85%=1.41, EROEI_CCS - min=7.98, max=9.51

Generate Figure 2. This is a 2 by 2 multiplot, with each sub-plot being a filled contour area plot, sometimes also referred to as a colormap or shaded area plot. Here we have plotted the 4-dimensional variable $EROEI_{CCS}$ in the following manner:

  • Columns represent different fossil fuels - with their corresponding $EROEI_{fuel}, EROEI_{el}, R$ values.
  • Rows represent different Capacity Factor $CF$ values.
  • The $x$ and $y$ axes represent ample variation ranges for $f_{op}$ and $f_{cap}$, respectively.
  • We have indicated with a solid rectangle in each sub-plot the reasonable ranges for $f_{op}$ and $f_{cap}$. We have repeated this for an additional Capture Ratio $CR$ value, using a dashed rectangle.
  • The $EROEI_{CCS}$ values are represented by shaded contours. The levels of the contours are indicated in the attached colorbars, per each fossil fuel type (so each column has its own colorbar).
In [62]:
#Plotting options
fig,axes=plt.subplots(2,2,figsize=(10,8))
plt.subplots_adjust(hspace=0.1)
plt.subplots_adjust(wspace=0.25)

eroei_ccs_scale_min2=1.2
eroei_ccs_scale_max2=2.3
ls=['-','--']
ls0='-'
lw0=0.5
t_x_offset=[6.5,12.5,6]
t_y_offset=[47,45,14]
breaker=['\n','\n','\n']
manual_locations=[[[(12,20),(20,20),(30,21),(35,20),(35,60)],\
                  [(7,70),(32,20),(35,60)]],\
                  [[(15,20),(23,19),(30,21),(35,20),(35,60)],\
                  [(7,70),(35,30),(35,70)]]]
levelsi=[[np.arange(11,24,1),np.arange(4,10.5,0.5)],[np.arange(14,27,1),np.arange(6,12.5,0.5)]]
for i in range(len(axes)):
    for j in range(len(axes[i])):
        eroei_el=eroei_el_cf[j][i]
        ax=axes[i][j]
        X, Y = np.meshgrid(fop, fcap)
        z = eroei_ccs[j][i][0]#*100.0/eroei_el #third index is CR, keep at 0, its relative
        levels = levelsi[i][j]
        cmap = plt.get_cmap('coolwarm_r')
        im = ax.contourf(X, Y, z, cmap=cmap, levels=levels,extend='both')
        im2 = ax.contour(X, Y, z, levels=levels,extend='both',colors='k',linewidths=0.4)
        ax.clabel(im2, fmt = '%.i',inline=1, fontsize=11, manual=manual_locations[i][j])
        ax.set_xlim((fop.min(),fop.max()))
        ax.set_ylim((fcap.min(),fcap.max()))
        
        #fig.subplots_adjust(right=0.9)
        cbar_ax = fig.add_axes([0.48+0.43*j, 0.12+0.395*(1-i), 0.015, 0.365])
        fig.colorbar(im, cax=cbar_ax)
        if j==1: cbar_ax.set_ylabel(u'$EROEI_{CCS}$',fontsize=14,labelpad=10)

        for m in range(len(CR_ratio)):
            n=0
            for ti in range(len(t0s)):
                if (((j==1)and(ti==2))or(ti==j)):
                    t0=t0s[ti]
                    t=[t0[0]*CR_ratio[m],t0[1]*CR_ratio[m],t0[2]+n,t0[3]+n]
                    rect = patches.Rectangle((t[0],t[2]),t[1]-t[0],t[3]-t[2],ls=ls[m],
                                             edgecolor='k',facecolor='none')
                    ax.add_patch(rect)
                    '''
                    if m==1:
                        ax.text(t[0]+2,t[3]-14,'CR\n'+str(CR[m]),horizontalalignment='center')
                    else:
                        ax.text(t[0]+t_x_offset[ti],t[3]-t_y_offset[ti],'CR'+breaker[ti]+str(CR[m]),
                                horizontalalignment='center')
                    '''
                    if m==0: 
                        ax.text(t[1]-0.5,t[3]-15,names[ti],horizontalalignment='right',weight='bold')
                    if (ti==1)and(m==1):
                        line = lines.Line2D((t0[0],t[1]),((t0[2]+t0[3])/2,(t[2]+t[3])/2),color='k',
                                            linewidth=lw0,ls=ls0)
                        #ax.add_line(line)  

        '''
        if j==0:
            ax.set_title(names[0].replace('\n',' ')+' - $EROEI_{el}='+\
                     str(round(eroei_el_cf[0][i],1))+\
                         u'$\nCapacity factor $CF='+str(CF[i])+'$%')
        else:
            ax.set_title(names[1].replace('\n',' ')+' - $EROEI_{el}='+\
                     str(round(eroei_el_cf[1][i],1))+'$\n'+\
                         names[2].replace('\n',' ')+' - $EROEI_{el}='+\
                     str(round(eroei_el_cf[2][i],1))+\
                         u'$\nCapacity factor $CF='+str(CF[i])+'$%')
        '''
        ax.text(0.93, 0.95, string.ascii_lowercase[j*len(axes)+i],
            horizontalalignment='left',
            verticalalignment='top',
            transform=ax.transAxes,size=13,alpha=0.7,
                          bbox=dict(boxstyle="round", fc="w", alpha=0.7))
            
axes[1][1].set_xlabel('Operational energy penalty $f_{op}$  (%)',fontsize=13,x=0)
axes[0][0].set_ylabel('Capital energy penalty $f_{cap}$  (%)',fontsize=13,y=0)        
axes[0][1].set_yticklabels('')   
axes[1][1].set_yticklabels('')   
axes[0][0].set_xticklabels('')   
axes[0][1].set_xticklabels('')  
#plt.suptitle(title+'\nframes represent different capture ratios CR (%)',fontsize=13,y=1.04)
plt.savefig('plot/png/fig2.png',dpi=300,bbox_inches = 'tight', pad_inches = 0.1)
plt.savefig('plot/eps/fig2.eps',format='eps',bbox_inches = 'tight', pad_inches = 0.1)
plt.savefig('plot/pdf/fig2.pdf',format='pdf',bbox_inches = 'tight', pad_inches = 0.1)
plt.show()

Dispatchable Renewables Part

Equations

In [63]:
l=Symbol('lambda')
epse=Symbol('epsilon')
var('ESOI D');
eq_S7=Eq(ESOI,l*e*D/epse)
eq_S7
Out[63]:
$$ESOI = \frac{D \eta \lambda}{\epsilon}$$
In [64]:
epsp=Symbol('epsilon_p')
epss=Symbol('epsilon_s')
var('D P C');
eq_S8=Eq(ESOI,C*l*e*D/(P*epsp+C*epss))
eq_S8
Out[64]:
$$ESOI = \frac{C D \eta \lambda}{C \epsilon_{s} + P \epsilon_{p}}$$
In [65]:
f=Symbol('varphi_j')
Ei=Symbol('E_{in_i}')
Ej=Symbol('E_{in_j}')
var('EROEI_disp c j i I J');
eq_S9=Eq(EROEI_disp,(E_out-(1-e)*Sum(f*E_out,(j,1,J))-c*E_out)/(Sum(Ei,(i,1,I))+Sum(Ej,(j,1,J))))
eq_S9
Out[65]:
$$EROEI_{disp} = \frac{- E_{out} c + E_{out} - \left(- \eta + 1\right) \sum_{j=1}^{J} E_{out} \varphi_{j}}{\sum_{i=1}^{I} E_{in_i} + \sum_{j=1}^{J} E_{in_j}}$$
In [66]:
a=Symbol('alpha_i')
ej=Symbol('eta_j')
f=Symbol('varphi_j')
Ei=Symbol('E_{in_i}')
Ej=Symbol('E_{in_j}')
var('EROEI_i u_j ESOI_j');
eq_S10=Eq(EROEI_disp,(E_out*(1-(1-e)*Sum(f,(j,1,J))-c))/(Sum(E_out*a/EROEI_i,(i,1,I))+Sum((E_out*f*ej)/(ESOI_j*u_j),(j,1,J))))
eq_S10
Out[66]:
$$EROEI_{disp} = \frac{E_{out} \left(- c - \left(- \eta + 1\right) \sum_{j=1}^{J} \varphi_{j} + 1\right)}{\sum_{i=1}^{I} \frac{E_{out} \alpha_{i}}{EROEI_{i}} + \sum_{j=1}^{J} \frac{E_{out} \eta_{j} \varphi_{j}}{ESOI_{j} u_{j}}}$$
In [67]:
eq_S11=Eq(EROEI_disp,((1-(1-e)*Sum(f,(j,1,J))-c))/(Sum(a/EROEI_i,(i,1,I))+Sum((f*ej)/(ESOI_j*u_j),(j,1,J))))
eq_S11
Out[67]:
$$EROEI_{disp} = \frac{- c - \left(- \eta + 1\right) \sum_{j=1}^{J} \varphi_{j} + 1}{\sum_{i=1}^{I} \frac{\alpha_{i}}{EROEI_{i}} + \sum_{j=1}^{J} \frac{\eta_{j} \varphi_{j}}{ESOI_{j} u_{j}}}$$

Simulation

Define $EROEI_{disp}$ calcuator function eroei_disp_cycler, including minimums and maximums over Storage fraction $\varphi$. It takes as inputs:

  • an array ESOI, of length how many storage scenarios we are considering, holding $ESOI$ values for each
  • an array phi, of arbitrary length, with values between $0$ and $1$, corresponding to the simulated storage fraction levels
  • an array eta, holding the Roundtrip efficiency for each storage scenario
  • an array c, holding the Curtailment for each storage scenario
  • an array eroei_el, of arbitrary length, with values between $1$ and $100$, corresponding to the $EROEI_{el}$ levels

Return one multidimensional array eroei_disp, corresponding to the Dispatchable Scalable Renewable EROEI $EROEI_{disp}$ for each parameter combination.

In [68]:
def eroei_disp_cycler(ESOI, phi, eta, c, eroei_el):
    #eroei_disp=np.zeros([len(ESOI),len(phi),len(eroei_el)])
    eroei_disp=np.zeros([len(ESOI),len(eroei_el),len(phi)])
    for i in range(len(ESOI)):
        for j in range(len(phi)):
            for k in range(len(eroei_el)):
                #eq_S11
                #Also Eq. 6 from http://pubs.rsc.org/en/content/articlepdf/2013/ee/c3ee41973h
                eroei_disp[i][k,j]=(-c[i]-(+1-eta[i])*phi[j]+1)/\
                                    ((1.0/eroei_el[k])+(eta[i]*phi[j]/ESOI[i]))    
    return eroei_disp

Generate Figure 3. This is a 2 by 2 multiplot, with each sub-plot being a filled contour area plot, sometimes also referred to as a colormap or shaded area plot. Here we have plotted the 4-dimensional variable $EROEI_{disp}$ in the following manner:

  • Sub-plots represent storage scenarios, each with their proper Energy Stored on Energy Invested $ESOI$ values and Roundtrip efficiencies $\eta$. The bottom left case also includes Curtailment $k$, in all other cases $k=0$.
  • The $x$ and $y$ axes represent ample variation ranges for Storage fraction $\varphi$ and $EROEI_{el}$, respectively.
  • The $EROEI_{disp}$ values are represented by shaded contours. The levels of the contours are indicated in the attached colorbar.
In [69]:
segments=800
eroei_el=np.linspace(10,60,segments)
phi=np.linspace(0,35,segments)
ESOI=[11,24,11,249]
ESOI_tech=['Battery storage','Hydrogen storage','Battery + Curtailment','Pumped hydro storage']
eta=np.array([83,30,83,80])
c=[0,0,0.07,0]

eroei_disp=eroei_disp_cycler(ESOI, phi/100.0, eta/100.0, c, eroei_el)

fig,axes=plt.subplots(2,2,figsize=(10,8))
plt.subplots_adjust(wspace=0.1)
plt.subplots_adjust(hspace=0.15)
contour_min=10
contour_max=60
levels = MaxNLocator(nbins=10).tick_values(contour_min, contour_max)
manual_locations=[[(30,10),(7,20),(7,34),(7,50),(3,55)],\
                          [(30,10),(7,20),(7,30),(7,45),(7,55)],\
                          [(30,10),(7,25),(7,42),(7,55)],\
                          [(7,20),(7,30),(7,40),(7,50)]]
for i in range(len(ESOI)):
        ax=axes[i//2][i%2]
        z = eroei_disp[i]
        cmap = plt.get_cmap('coolwarm_r')
        X, Y = np.meshgrid(phi, eroei_el)
        im = ax.contourf(X, Y, z, cmap=cmap, levels=levels,extend='both')
        im2 = ax.contour(X, Y, z, levels=levels,extend='both',colors='k',linewidths=0.4)    
        ax.clabel(im2, fmt = '%i',inline=1, fontsize=11, manual=manual_locations[i])
        if (i==0):
            fig.subplots_adjust(right=0.9)
            cbar_ax = fig.add_axes([0.93, 0.12, 0.015, 0.77])
            fig.colorbar(im, cax=cbar_ax)
            cbar_ax.set_ylabel(u'$EROEI_{disp}$',fontsize=13,labelpad=10)
        ax.set_xlim((phi.min(),phi.max()))
        ax.set_ylim((eroei_el.min(),eroei_el.max()))
        ax.set_title(ESOI_tech[i],fontsize=13)
        ax.text(0.93, 0.95, string.ascii_lowercase[(i%2)*2+i//2],
                horizontalalignment='left',
                verticalalignment='top',
                transform=ax.transAxes,size=13,alpha=0.7,
                              bbox=dict(boxstyle="round", fc="w", alpha=0.7))
            
axes[1][1].set_xlabel('Storage fraction $\\varphi$ (%)',fontsize=13,x=0)
axes[0][0].set_ylabel('$EROEI_{el}$',fontsize=13,y=0)        
axes[0][1].set_yticklabels('')   
axes[1][1].set_yticklabels('')   
axes[0][0].set_xticklabels('')   
axes[0][1].set_xticklabels('')  
#plt.suptitle(r'$EROEI_{disp}$ for Scalable Dispatchable Renewables - PV and Wind',fontsize=13,y=0.98)
plt.savefig('plot/png/fig3.png',dpi=150,bbox_inches = 'tight', pad_inches = 0.1)
plt.savefig('plot/eps/fig3.eps',format='eps',bbox_inches = 'tight', pad_inches = 0.1)
plt.savefig('plot/pdf/fig3.pdf',format='pdf',bbox_inches = 'tight', pad_inches = 0.1)
plt.show()

Comparison

Sensitivity analysis

In this section we conduct a comparative sensitivity analysis over a range of parameters presented in the Fossils and Renewables sections, in order to be able to represent $EROEI_{CCS}$ and $EROEI_{disp}$ in one common axis, for a range of $EROEI_{el}$ values. In essence this means replicating the rectangles of Figure 2. subplots and taking their minimum and maximum ranges, over many starting $EROEI_{fuel}$ - and thus varying $EROEI_{el}$ and $R$ - values, to result in trapezoids. This is then compared with the 4 energy storage cases presented in Figure 3., with separate $EROEI_{el}$ range assumptions for Wind and PV.

First, we define a convex hull function to get outermost points of trapezoids in the sensitivity analysis area.

In [39]:
from scipy.spatial import ConvexHull
In [40]:
def convexhull(p):
    p = np.array(p)
    hull = ConvexHull(p)
    return p[hull.vertices,:]

Define wieldy color functions.

In [41]:
def cbrew(i):
    colors=['#984ea3',
            '#e41a1c',
            '#377eb8',
            '#018571',
            '#e66101']
    return colors[i]
def cbrew2(i):
    colors=['#984ea3',
            '#e41a1c',
            '#377eb8',
            '#4daf4a',
            '#ff7f00']
    return colors[i]

This is main sensitivity analysis calculation. We calculate envelopes over a range of parameters for Fossils and Renewables.

Fossils parameters:

  • $EROEI_{fuel}$ varies over $20-87$ for Gas and $23-58$ for Coal, using $800$ steps in between. This is stored in the variable eroei_el_ns with the number of steps indicated in segments.
  • The other parameters have been already defined in the generation of Figure 2. We re-utilize their variation ranges here.
In [42]:
#Fossil part
segments=800
eroei_el_ns=[np.linspace(20,87,segments),np.linspace(23,58,segments),np.linspace(23,58,segments)]   
eroei_el_cf_N, eroei_ccs_N, eroei_ccs_min_N, eroei_ccs_max_N=[],[],[],[]
for i in range(len(names)):
    print('calculating '+names[i].replace('\n',' ')+' area...')
    eroei_el_cf_n, eroei_ccs_n, eroei_ccs_min_n, eroei_ccs_max_n=[],[],[],[]
    for eroi in range(len(eroei_el_ns[i])):
        eroei_el_cf, R_cf, eroei_ccs, eroei_ccs_min, eroei_ccs_max=\
                eroei_ccs_cycler([names[i]], [eroei_el_ns[i][eroi]], [t0s[i]], CF, CR_ratio, fop, fcap)
        eroei_el_cf_n.append(eroei_el_cf)
        eroei_ccs_n.append(eroei_ccs)
        eroei_ccs_min_n.append(eroei_ccs_min)
        eroei_ccs_max_n.append(eroei_ccs_max)
    eroei_el_cf_N.append(eroei_el_cf_n)
    eroei_ccs_N.append(eroei_ccs_n)
    eroei_ccs_min_N.append(eroei_ccs_min_n)
    eroei_ccs_max_N.append(eroei_ccs_max_n)
calculating NGCC Gas area...
calculating Pulverized Coal area...
calculating IGCC Coal area...

OPTIONAL. Save generated data for easier loading later. No need to run these two cells if you have run the previous cell. GitHub does not hold large files, so you would have to generate your own data by running the sensitivity cell as least once.

In [ ]:
#optional SAVE
np.save('eroei_el_cf_N',eroei_el_cf_N)
np.save('eroei_ccs_N',eroei_ccs_N)
np.save('eroei_ccs_min_N',eroei_ccs_min_N)
np.save('eroei_ccs_max_N',eroei_ccs_max_N)
In [68]:
#optional LOAD
eroei_el_cf_N=np.load('eroei_el_cf_N.npz')
eroei_ccs_N=np.load('eroei_ccs_N.npz')
eroei_ccs_min_N=np.load('eroei_ccs_min_N.npz')
eroei_ccs_max_N=np.load('eroei_ccs_max_N.npz')

Renewables parameters:

  • $EROEI_{el}$ varies over $9-30$ for PV and $20-40$ for Wind, using segments number of steps (set to $800$ in the Fossils cell. This is stored in eroei_el_ren.
  • Higher resolution areas for PV between $9-17$ for moderate insolation and $12-23$ for high insolation systems. We also indicate a future expected range of $20-40$ by the year 2022.
  • Both current and future feasible wind $EROEI_{el}$ range is defined as $20-40$.
  • We also indicate a likely future $EROEI_{disp}$ example for the European electricity grid, using a mix of PV and Wind technologies, with shares $0.33$ and $0.67$, respectively, corresponding $EROEI_{el}$s of $18$ and $35$, and a resulting composite $EROEI_{disp}$ of $21.9$.
  • We re-utilize the energy storage cases and their respective parameters $ESOI, \eta, k, \varphi$ from Figure 3.
In [43]:
#RE part
eroei_el_ren=[np.linspace(9,30,segments),np.linspace(20,40,segments)]
eroei_el_pv=[[9,17],[12,23],[20,40]] #current medium and high PV EROEI_el ranges, future 2022 range
eroei_el_wind=[20,40] #current medium and high PV EROEI_el ranges, future 2022 range
eroei_el_eu=[[18,35],[0.33,0.67],21.9] #Table S3, current EU solar PV and wind EROEI_el, shares and resulting EROEI_disp
eroei_ren_names=['PV','Wind']
eroei_disp_n=[]
for i in range(len(eroei_ren_names)):
    print('calculating '+eroei_ren_names[i]+' area...')
    eroei_disp=eroei_disp_cycler(ESOI, phi/100.0, eta/100.0, c, eroei_el_ren[i])
    eroei_disp_n.append(eroei_disp)  
calculating PV area...
calculating Wind area...
In [142]:
#optional SAVE
np.save('eroei_disp_n',eroei_disp_n)
In [ ]:
#optional LOAD
eroei_disp_n=np.load('eroei_disp_n.npz')

We convert the sensitivity analysis mutlidimensional arrays into points on the canvas.

In [44]:
#Calculate points for convex hull
all_points={}

#Renewables
for i in range(len(eroei_disp_n)):
    print('calculating '+eroei_ren_names[i]+' points...')
    if eroei_ren_names[i] not in all_points:all_points[eroei_ren_names[i]]=[[] for ii in range(len(eroei_disp_n[i]))]
    for j in range(len(eroei_disp_n[i])):
        for eroi in range(len(eroei_disp_n[i][j])):
            x0=eroei_el_ren[i][eroi]
            y0min=min(eroei_disp_n[i][j][eroi].flatten())
            y0max=max(eroei_disp_n[i][j][eroi].flatten())
            all_points[eroei_ren_names[i]][j].append([x0,y0min])
            all_points[eroei_ren_names[i]][j].append([x0,y0max])

#Fossil part
for i in range(len(names)):
    print('calculating '+names[i].replace('\n',' ')+' points...')
    if names[i] not in all_points:all_points[names[i]]=[]
    for eroi in range(len(eroei_el_ns[i])):
        for j in range(len(CF)):
            x0=eroei_el_cf_N[i][eroi][0][j]
            y0min=min(eroei_ccs_min_N[i][eroi][0][j])
            y0max=max(eroei_ccs_max_N[i][eroi][0][j])
            all_points[names[i]].append([x0,y0min])
            all_points[names[i]].append([x0,y0max])
calculating PV points...
calculating Wind points...
calculating NGCC Gas points...
calculating Pulverized Coal points...
calculating IGCC Coal points...

Intersections

Define a simple function for finding nearest values in an array, in order to be able to identify intersections between Fossil $EROEI_{CCS}$ and Renewable $EROEI_{disp}$ cases.

In [45]:
def find_nearest(array, value):
    array = np.asarray(array)
    idx = (np.abs(array - value)).argmin()
    return idx

We define the $EROEI_{fuel}$ values, for each Fossil fuel source to calculate the intersections for. We reasonably use a Capacity Factor of $CF=85$% and a Capture Ratio of $CR=90$%. For the EROEI we select the highest value and the mid-point of the defined sensitvity range, for each fuel.

Once this is done, we cycle through all Fossil-Renewable combinations and take note of the matching $(EROEI_{el},\varphi)$ combinations, for all cases when $EROEI_{CCS}=EROEI_{disp}$.

In [46]:
eroei_ids=[-1,segments//2] #Fossil EROEI range IDs to calculate intersections for. Default: LAST and MID
eroei_ida=['high','mid'] #corresponding scenario names
CF_id=1 #85%
CR_id=0 #90%      
all_intersections=[]
for eroei_idi in range(len(eroei_ids)):
    eroei_id=eroei_ids[eroei_idi]
    case_id=eroei_ida[eroei_idi]
    for name_id in range(len(names)):
        name=names[name_id].replace('\n',' ')
        fs_x=eroei_el_cf_N[name_id][eroei_id][0][CF_id]
        fs_y_max=eroei_ccs_max_N[name_id][eroei_id][0][CF_id][CR_id]
        fs_y_min=eroei_ccs_min_N[name_id][eroei_id][0][CF_id][CR_id]
        fs_y=np.round((fs_y_max+fs_y_min)/2.0,2)
        for i in range(len(eroei_disp_n))[::-1]:
            for j in range(len(eroei_disp_n[i])):
                minimums=all_points[eroei_ren_names[i]][j][0::2] #bottom line on RE area plots
                maximums=all_points[eroei_ren_names[i]][j][1::2] #top line on RE area plots    
                intersections=[]
                for eroi in range(len(eroei_el_ren[i])):
                    if minimums[eroi][1]<fs_y<maximums[eroi][1]:
                        intersections.append(eroi)
                if len(intersections)>0:
                    imin=min(intersections)
                    rmin=np.round(eroei_el_ren[i][imin],1)
                    emin=eroei_disp_n[i][j][imin]
                    phimin=np.round(phi[find_nearest(emin,fs_y)],1)
                    imax=max(intersections)
                    rmax=np.round(eroei_el_ren[i][imax],1)
                    emax=eroei_disp_n[i][j][imax]
                    phimax=np.round(phi[find_nearest(emax,fs_y)],1)
                else: 
                    imax=np.nan
                    emax=np.nan
                    rmax=np.nan
                    phimax=np.nan
                    
                if phimax>34:phimax=35.0
                if phimin<0.31:phimin=0.0
                
                all_intersections.append([case_id,name,eroei_ren_names[i],ESOI_tech[j],
                    np.round(eroei_el_ns[name_id][eroei_id],1),np.round(fs_x,1),np.round(fs_y,1),
                    rmin,phimin,rmax,phimax])

all_intersections=pd.DataFrame(all_intersections)
all_intersections.columns=['EROEI case','Fossil','Renewable','Storage',
                           'Fossil $EROEI_{fuel}$','Fossil $EROEI_{el}$','Fossil $EROEI_{CCS}$',
                           'RE $EROEI_{el}$ min',r'$\varphi$ min (%)',
                           'RE $EROEI_{el}$ max',r'$\varphi$ max (%)']
all_intersections['Case ID']=[i+1 for i in range(len(all_intersections))]
all_intersections2=all_intersections.dropna()
all_intersections2['Case ID']=[i+1 for i in range(len(all_intersections2))]
C:\Users\csala\AppData\Local\Continuum\anaconda2\envs\python3\lib\site-packages\ipykernel_launcher.py:52: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy

We present the obtained intersections in Table S4, sorting by EROEI case, Fossil fuel, Renewable type and Energy storage scenario.

In [47]:
all_intersections2.set_index(['EROEI case','Fossil','Renewable','Storage','Case ID'])
Out[47]:
Fossil $EROEI_{fuel}$ Fossil $EROEI_{el}$ Fossil $EROEI_{CCS}$ RE $EROEI_{el}$ min $\varphi$ min (%) RE $EROEI_{el}$ max $\varphi$ max (%)
EROEI case Fossil Renewable Storage Case ID
high NGCC Gas Wind Battery storage 1 87.0 31.0 21.3 21.3 0.0 40.0 26.2
Hydrogen storage 2 87.0 31.0 21.3 21.3 0.0 32.2 35.0
Battery + Curtailment 3 87.0 31.0 21.3 22.9 0.0 40.0 22.3
Pumped hydro storage 4 87.0 31.0 21.3 21.3 0.0 23.5 35.0
PV Battery storage 5 87.0 31.0 21.3 21.3 0.0 30.0 16.3
Hydrogen storage 6 87.0 31.0 21.3 21.3 0.0 30.0 29.9
Battery + Curtailment 7 87.0 31.0 21.3 22.9 0.0 30.0 12.4
Pumped hydro storage 8 87.0 31.0 21.3 21.3 0.0 23.5 35.0
Pulverized Coal PV Battery storage 9 58.0 13.5 8.4 9.0 8.7 11.6 35.0
Hydrogen storage 10 58.0 13.5 8.4 9.0 8.7 11.6 35.0
Battery + Curtailment 11 58.0 13.5 8.4 9.0 0.4 12.9 35.0
Pumped hydro storage 12 58.0 13.5 8.4 9.0 30.8 9.1 35.0
IGCC Coal PV Battery storage 13 58.0 11.5 8.0 9.0 14.2 11.0 35.0
Hydrogen storage 14 58.0 11.5 8.0 9.0 13.8 11.1 35.0
Battery + Curtailment 15 58.0 11.5 8.0 9.0 5.2 12.2 35.0
mid NGCC Gas Wind Battery storage 16 53.5 21.4 15.6 20.0 16.5 29.4 35.0
Hydrogen storage 17 53.5 21.4 15.6 20.0 24.8 22.6 35.0
Battery + Curtailment 18 53.5 21.4 15.6 20.0 11.3 33.8 35.0
PV Battery storage 19 53.5 21.4 15.6 15.6 0.0 29.4 35.0
Hydrogen storage 20 53.5 21.4 15.6 15.6 0.0 22.6 35.0
Battery + Curtailment 21 53.5 21.4 15.6 16.8 0.0 30.0 30.6
Pumped hydro storage 22 53.5 21.4 15.6 15.6 0.0 17.0 35.0
Pulverized Coal PV Hydrogen storage 23 40.5 10.6 6.8 9.0 31.0 9.4 35.0
Battery + Curtailment 24 40.5 10.6 6.8 9.0 25.4 9.8 35.0
IGCC Coal PV Hydrogen storage 25 40.5 9.2 6.6 9.0 34.8 9.0 35.0
Battery + Curtailment 26 40.5 9.2 6.6 9.0 30.4 9.4 35.0

Area plots

In this section we represent the sensitivity analysis results from both the Fossils and the Renewables in one graph.

For aesthetics, we define a custom shaded area plot, which fades towards the edges.

In [48]:
def alpharater(g,n,t):
    if g<t*n:
        return (g+1)*1.0/(t*n)
    elif g+1>(1-t)*n:
        return ((n-g))*1.0/(t*n)
    else:
        return 1
In [49]:
#Renewables plotter
def plot_renewables(axes,segment_wise=False,more_labels=False,edgeless=False,blow_up=False):
    for i in range(len(eroei_disp_n))[::-1]:
        print('plotting '+eroei_ren_names[i]+' area...')
        for j in range(len(eroei_disp_n[i])):
            points=all_points[eroei_ren_names[i]][j]
            ax=axes[j//2][j%2]
            colr=cbrew(4-i)
            if segment_wise:
                nfaces=400
                if not blow_up:
                    tfaces=0.2
                    basealpha=0.15 #0.2
                    lws=0.95
                else:
                    tfaces=0.2
                    basealpha=0.3 #0.4
                    lws=1.7
                for g,p in enumerate(np.split(np.array(points),nfaces)):
                    psegment=convexhull(p)
                    poly = plt.Polygon(psegment, edgecolor=colr,
                                       alpha=basealpha*alpharater(g,nfaces,tfaces),fill=False,lw=lws)
                    ax.add_patch(poly)
                    if not edgeless:
                        ax.plot(psegment.T[0][:2],psegment.T[1][:2],c=colr,alpha=0.9,lw=0.5)
                        ax.plot(psegment.T[0][2:],psegment.T[1][2:],c=colr,alpha=0.9,lw=0.5)
                        if g==0:
                            ax.plot(psegment.T[0][0::2],psegment.T[1][0::2],c=colr,alpha=0.9,lw=0.5)
                        if g==len(np.split(np.array(points),nfaces))-1:
                            ax.plot(psegment.T[0][1::2],psegment.T[1][1::2],c=colr,alpha=0.9,lw=0.5)
            else:
                poly = plt.Polygon(convexhull(points), edgecolor="none",facecolor=colr,alpha=0.4)
                poly2 = plt.Polygon(convexhull(points), edgecolor=colr,alpha=0.8,fill=False,lw=1)
                ax.add_patch(poly)
                ax.add_patch(poly2)
            legendpoly = plt.Polygon([[0,0],[0,-1],[-1,-1]], edgecolor=colr,facecolor=colr,alpha=legendalpha,lw=1,label=eroei_ren_names[i])
            ax.add_patch(legendpoly)
            if (i==0):
                ax.set_title(ESOI_tech[j],fontsize=13)
                ax.text(0.04, 0.12, string.ascii_lowercase[(j%2)*2+j//2],
                        horizontalalignment='left',
                        verticalalignment='top',
                        transform=ax.transAxes,size=13,alpha=0.7,
                        bbox=dict(boxstyle="round", fc="w", alpha=0.7))
                ax.plot([0,100],[0,100],c='grey',lw=1,ls='--')
                if more_labels:
                    if (j==2): 
                        ax.text(0.22, 0.27, 'no storage',
                                horizontalalignment='center',
                                verticalalignment='center',
                                transform=ax.transAxes,size=11,rotation=40, color='grey')
                        ax.arrow(35, 35, 0, -1.5, head_width=0.5, head_length=0.5, fc='k', ec='k',ls='-',lw=0.5)
                        ax.arrow(35, 32, 0, -15, head_width=0.5, head_length=0.5, fc='k', ec='k',ls='-',lw=0.5)
                        ax.text(34.5, 34, 'curtailment',
                                horizontalalignment='right',
                                verticalalignment='center', rotation=0, color='k')
                        ax.text(34.5, 20, 'more\nstorage',
                                horizontalalignment='right',
                                verticalalignment='center', rotation=0, color='k')
                    if (j==3):
                        ax.arrow(31.5, 18.7, 0, 6, head_width=0.5, head_length=0.5, fc='k', ec='k',ls='-',lw=0.5)
                        ax.arrow(31.5, 18.7, -4, -4, head_width=0.5, head_length=0.5, fc='k', ec='k',ls='-',lw=0.5)
                        ax.text(32.5, 21.7, 'lower\ncapture\nratio',
                                horizontalalignment='left',
                                verticalalignment='center', rotation=0, color='k')
                        ax.text(30, 12.2, 'lower\ncapacity\nfactor',
                                horizontalalignment='left',
                                verticalalignment='center', rotation=0, color='k')
                    if (j==1):
                        arrowoffset=[1.5,0.75,1.5]
                        arrowtextoffset=[1.5,4.5,1]
                        arrowtext=['PV moderate       ','   PV high','PV in 2022']
                        arrowtextoffsetx=[1,1,1]
                        arrowtextoffsety=[35.75,32.5,29.25]
                        #From main text + Breyer estimate
                        for a in range(len(eroei_el_pv)):
                            x1,y1=eroei_el_pv[a][0], arrowoffset[a]
                            w1=eroei_el_pv[a][1]-eroei_el_pv[a][0]
                            ax.arrow(x1,y1,w1,0,length_includes_head=True,
                                     head_width=0.5, head_length=0.5, fc=colr, ec=colr,ls='-',lw=0.5)
                            ax.arrow(x1+w1,y1,-w1,0,length_includes_head=True,
                                     head_width=0.5, head_length=0.5, fc=colr, ec=colr,ls='-',lw=0.5)
                            ax.text(arrowtextoffset[a]+x1+w1/2, 3, arrowtext[a],
                                horizontalalignment='center', 
                                verticalalignment='center', rotation=0, color=colr)
                            ax.text(arrowtextoffsetx[a], arrowtextoffsety[a], arrowtext[a].strip()+\
                                ': $'+str(eroei_el_pv[a][0])+'-'+str(eroei_el_pv[a][1])+'$',
                                horizontalalignment='left', 
                                verticalalignment='top', rotation=0, color=colr)

                        x2=eroei_el_eu[0][0]*eroei_el_eu[1][0]+eroei_el_eu[0][1]*eroei_el_eu[1][1]
                        y2=eroei_el_eu[2]
                        ax.scatter([x2],[y2],color='k',zorder=100)
                        #From Table S2 in supplementary
                        ax.text(x2*1.04, y2*1.03, 'European\nsystem\nexample',
                                horizontalalignment='left',
                                verticalalignment='center', rotation=0, color='k')
                        
            #From Davidsson, S., Höök, M. & Wall, G. A review of life cycle assessments on wind energy systems
            elif (j==1):
                if more_labels:
                    xk=3
                    x1,y1=eroei_el_wind[0], 6
                    w1=eroei_el_wind[1]-eroei_el_wind[0]
                    ax.arrow(x1+w1-xk,y1,xk,0,length_includes_head=True,
                                 head_width=0.5, head_length=0.5, fc=colr, ec=colr,ls='--',lw=0.5)
                    ax.arrow(x1+w1-xk,y1,-w1+xk,0,length_includes_head=True,
                                 head_width=0.5, head_length=0.5, fc=colr, ec=colr,ls='-',lw=0.5)
                    ax.text(x1+w1/2, y1+1.5, 'Wind',
                            horizontalalignment='center', 
                            verticalalignment='center', rotation=0, color=colr)
                    ax.text(1, 39, 
                        'Wind: $EROEI_{el}='+str(eroei_el_wind[0])+'-'+str(eroei_el_wind[1])+'$+',
                        horizontalalignment='left', 
                        verticalalignment='top', rotation=0, color=colr)
In [50]:
#Fossil plotter
def plot_fossils(axes,segment_wise=False,more_labels=False,edgeless=False,blow_up=False):
    fossil_minmax=[[],[]]
    for i in range(len(names)):
        points=all_points[names[i]]
        print('plotting '+names[i].replace('\n',' ')+' area...')
        for axi in range(len(axes)):
            for axj in range(len(axes[axi])):
                ax=axes[axi][axj]
                colr=cbrew(i)
                if segment_wise:
                    nfaces=100
                    if blow_up:
                        tfaces=0.2
                        basealpha=0.15 #0.2
                        lws=1.05
                    else:
                        tfaces=0.24
                        basealpha=0.11 #0.15
                        lws=1
                    for g,p in enumerate(np.split(np.array(points),nfaces)):
                        psegment=convexhull(p)
                        poly = plt.Polygon(psegment, edgecolor=colr, facecolor=colr,
                                           alpha=basealpha*alpharater(g,nfaces,tfaces),lw=lws)
                        ax.add_patch(poly)
                        if not edgeless: ax.plot([psegment.T[0][-1],psegment.T[0][-2]],[psegment.T[1][-1],psegment.T[1][-2]],
                                    c=colr,alpha=0.9,lw=0.5)
                        if g==len(np.split(np.array(points),nfaces))-1:
                            if not edgeless: ax.plot([psegment.T[0][0],psegment.T[0][-1]],[psegment.T[1][0],psegment.T[1][-1]],
                                    c=colr,alpha=0.9,lw=0.5)
                            if axi==0 and axj==0:
                                fossil_minmax[1].append(np.round(psegment.T[0][0],1))
                    psegment=convexhull(points)
                    if axi==0 and axj==0:
                        fossil_minmax[0].append(np.round(points[0][0],1))
                    for g in range(len(psegment.T[0])-2):
                        if not edgeless: ax.plot([psegment.T[0][g],psegment.T[0][g+1]],[psegment.T[1][g],psegment.T[1][g+1]],
                                    c=colr,alpha=0.9,lw=0.5)
                else:
                    poly = plt.Polygon(convexhull(points), edgecolor="k",facecolor=colr,alpha=0.4)
                    poly2 = plt.Polygon(convexhull(points), edgecolor=colr,alpha=0.8,fill=False,lw=1)
                    ax.add_patch(poly)
                    ax.add_patch(poly2)
                legendpoly = plt.Polygon([[0,0],[0,-1],[-1,-1]], edgecolor=colr,facecolor=colr,alpha=legendalpha,lw=1,label=names[i].replace('\n',' '))
                ax.add_patch(legendpoly)
        
        if more_labels:
            ax=axes[1][1]
            arrowoffset=[0.75,2.25,1.5]
            arrowtextoffsetx=[1,1,1]
            arrowtextoffsety=[39,32.5,26]
            arrowtextoffsetx2=[23,15,14.5]
            arrowtextoffsety2=[2.25,5,2.25]
            sep=['\t  ','\t ','\t    ']
            x1,y1=fossil_minmax[0][i], arrowoffset[i]
            w1=fossil_minmax[1][i]-fossil_minmax[0][i]
            ax.arrow(x1,y1,w1,0,length_includes_head=True,
                         head_width=0.5, head_length=0.5, fc=colr, ec=colr,ls='-',lw=0.5)
            ax.arrow(x1+w1,y1,-w1,0,length_includes_head=True,
                         head_width=0.5, head_length=0.5, fc=colr, ec=colr,ls='-',lw=0.5)
            ax.text(arrowtextoffsetx[i], arrowtextoffsety[i], 
                    '$EROEI_{th}='+str(int(min(eroei_el_ns[i])))+'-'+str(int(max(eroei_el_ns[i])))+'$'+\
                    '\n$EROEI_{el}='+str(fossil_minmax[0][i])+'-'+str(fossil_minmax[1][i])+'$',
                    horizontalalignment='left', 
                    verticalalignment='top', rotation=0, color=colr)
            ax.text(arrowtextoffsetx2[i], arrowtextoffsety2[i], names[i].replace('\n',''),
                    horizontalalignment='left', verticalalignment='center', rotation=0, color=colr)
In [51]:
#Intersections plotter
def plot_intersections(axes,blow_up=False):
    print('plotting Intersections...')
    for i in all_intersections2.index:
        j=ESOI_tech.index(all_intersections2.loc[i]['Storage'])
        z=names.index(all_intersections2.loc[i]['Fossil'].replace(' ','\n'))
        w=eroei_ren_names.index(all_intersections2.loc[i]['Renewable'])
        ax=axes[j//2][j%2]
        x=all_intersections2.loc[i]['Fossil $EROEI_{el}$']
        y=all_intersections2.loc[i]['Fossil $EROEI_{CCS}$']
        if (i-1)%4==0:
            for axi in axes:
                for axj in axi:
                    axj.scatter([x],[y],c='k',s=20,zorder=100)
                    axj.scatter([x],[y],c=cbrew(z),s=12,zorder=100)

        x1=float(all_intersections2.loc[i]['RE $EROEI_{el}$ min'])
        x2=float(all_intersections2.loc[i]['RE $EROEI_{el}$ max'])
        colr=cbrew(4-w)
        if (not blow_up and x1>13) or (blow_up and x2<maxlabel):
            ax.plot([x1,x2],[y,y],c='k',lw=2,alpha=0.3)
            ax.plot([x1,x2],[y,y],c=colr,lw=1,alpha=0.5)
            if blow_up: z=0.1
            else: z=1
            ax.text(x1+(x2-x1)/2, y+z+z*(w), all_intersections2.loc[i]['Case ID'],
                        horizontalalignment='center',
                        verticalalignment='bottom', fontsize=8,
                        bbox=dict(boxstyle="round", fc=colr, alpha=0.7))

We plot the Renewables and Fossils on a common axis, with $EROEI_{el}$ on the $x$ and $EROEI_{CCS}$ or $EROEI_{disp}$ on the $y$ axis. We split the figure into 4 sub-plots, split by enegeyr storage cases, exactly as in Figure 3. Each sub-plot will have 5 shaded trapezoid areas, for Wind, PV, NGCC Gas, Pulverized Coal and IGCC Coal, representing the extent ranges of the sensitivity analysis.

We produce a set of these area plots, the first of which is the explicit representation of the Intersection cases from Table S4.

In [52]:
#Plots init
fig,axes=plt.subplots(2,2,figsize=(9,8),subplot_kw=dict(facecolor='#EEEEEE',axisbelow=True))
for axi in axes:
    for ax in axi:
        ax.grid(color='white', linestyle='solid')                       
plt.subplots_adjust(wspace=0.1)
plt.subplots_adjust(hspace=0.15)

legendalpha=0.7

#Renewables part
plot_renewables(axes,segment_wise=True)
                                  
#Fossils part
plot_fossils(axes,segment_wise=True)

#Intersections
plot_intersections(axes)

#Format labels
axes[0][0].set_ylabel('$EROEI_{CCS}$ | $EROEI_{disp}$',fontsize=14,y=0)
axes[1][1].set_xlabel('$EROEI_{el}$',fontsize=14,x=0)        
maxlabel=max(np.array(eroei_el_ren).flatten())
labelrange=[i if i%2==0 else '' for i in range(int(maxlabel)+1) ]
for axi in axes:
    for ax in axi:
        ax.set_xlim(0,maxlabel)
        ax.set_ylim(0,maxlabel)
        ax.set_xticks(range(int(maxlabel)+1))
        ax.set_yticks(range(int(maxlabel)+1))
        ax.set_yticks(range(int(maxlabel)+1))
        ax.set_xticklabels(labelrange,fontsize=8)   
        ax.set_yticklabels(labelrange,fontsize=8) 

axes[0][0].legend(loc=(0.03,0.52),framealpha=0,fontsize=11)  
        
#plt.suptitle('Comparison of adjusted EROEIs for \nfossil fuels with carbon-capture and storage $EROEI_{CCS}$ and \ndispatchable renewables with energy storage $EROEI_{disp}$',y=1.03,fontsize=13)
plt.savefig('plot/extra/figS4.png',dpi=150,bbox_inches = 'tight', pad_inches = 0.1)
plt.show()
plotting Wind area...
plotting PV area...
plotting NGCC Gas area...
plotting Pulverized Coal area...
plotting IGCC Coal area...
plotting Intersections...

Next, we produce a blow-up of the same plot of the $EROEI_{el}$ range of $3-14$ for better visuals.

In [53]:
#Plots init
fig,axes=plt.subplots(2,2,figsize=(9,8),subplot_kw=dict(facecolor='#EEEEEE',axisbelow=True))
for axi in axes:
    for ax in axi:
        ax.grid(color='white', linestyle='solid')                       
plt.subplots_adjust(wspace=0.1)
plt.subplots_adjust(hspace=0.15)

legendalpha=0.7

#Renewables part
plot_renewables(axes,segment_wise=True,blow_up=True)

#Fossil part
plot_fossils(axes,segment_wise=True,blow_up=True)

minlabel=3
maxlabel=14
            
#Intersections
plot_intersections(axes,blow_up=True)

#Format labels
axes[0][0].set_ylabel('$EROEI_{CCS}$ | $EROEI_{disp}$',fontsize=14,y=0)
axes[1][1].set_xlabel('$EROEI_{el}$',fontsize=14,x=0)        
labelrange=[i for i in range(minlabel,int(maxlabel)+1) ]
for axi in axes:
    for ax in axi:
        ax.set_xlim(minlabel,maxlabel)
        ax.set_ylim(minlabel,maxlabel)
        ax.set_xticks(range(minlabel,int(maxlabel)+1))
        ax.set_yticks(range(minlabel,int(maxlabel)+1))
        ax.set_yticks(range(minlabel,int(maxlabel)+1))
        ax.set_xticklabels(labelrange,fontsize=8)   
        ax.set_yticklabels(labelrange,fontsize=8) 

axes[0][0].legend(loc=(0.03,0.52),framealpha=0,fontsize=11)  
        
#plt.suptitle('Comparison of adjusted EROEIs for \nfossil fuels with carbon-capture and storage $EROEI_{CCS}$ and \ndispatchable renewables with energy storage $EROEI_{disp}$',y=1.03,fontsize=13)
plt.savefig('plot/extra/figS5.png',dpi=150,bbox_inches = 'tight', pad_inches = 0.1)
plt.show()
plotting Wind area...
plotting PV area...
plotting NGCC Gas area...
plotting Pulverized Coal area...
plotting IGCC Coal area...
plotting Intersections...

And finally, we produce an annotated plots, with the $EROEI_{fuel}$s and $EROEI_{el}$ indicated as on-figure texts, as well as teh European system example presented in Table S3. This becomes Figure 4. in the paper.

In [54]:
#Plots init
fig,axes=plt.subplots(2,2,figsize=(9,8),subplot_kw=dict(facecolor='#FFFFFF',axisbelow=True))
for axi in axes:
    for ax in axi:
        ax.grid(color='white', linestyle='solid')                       
plt.subplots_adjust(wspace=0.1)
plt.subplots_adjust(hspace=0.15)
legendalpha=0.7

#Renewables part
plot_renewables(axes,segment_wise=True,more_labels=True)

#Fossil part
plot_fossils(axes,segment_wise=True,more_labels=True)

#Intersections
#plot_intersections(axes)

#Format labels
for axi in axes:
    for ax in axi:
        ax.set_xlim(0,max(np.array(eroei_el_ren).flatten()))
        ax.set_ylim(0,max(np.array(eroei_el_ren).flatten()))
axes[0][0].set_ylabel('$EROEI_{CCS}$ | $EROEI_{disp}$',fontsize=14,y=0)
axes[1][1].set_xlabel('$EROEI_{el}$',fontsize=14,x=-0.05)        
axes[0][1].set_yticklabels('')   
axes[1][1].set_yticklabels('')   
axes[0][0].set_xticklabels('')   
axes[0][1].set_xticklabels('')  
axes[0][0].legend(loc=(0.03,0.52),framealpha=0,fontsize=11)       
        
#plt.suptitle('Comparison of adjusted EROEIs for \nfossil fuels with carbon-capture and storage $EROEI_{CCS}$ and \ndispatchable renewables with energy storage $EROEI_{disp}$',y=1.03,fontsize=13)
plt.savefig('plot/png/fig4.png',dpi=150,bbox_inches = 'tight', pad_inches = 0.1)
plt.savefig('plot/eps/fig4.eps',format='eps',bbox_inches = 'tight', pad_inches = 0.1)
plt.savefig('plot/pdf/fig4.pdf',format='pdf',bbox_inches = 'tight', pad_inches = 0.1)
plt.show()
plotting Wind area...
plotting PV area...
plotting NGCC Gas area...
plotting Pulverized Coal area...
plotting IGCC Coal area...

An alternative iteration of Figure 4. with the area edges removed.

In [55]:
#Plots init
fig,axes=plt.subplots(2,2,figsize=(9,8),subplot_kw=dict(facecolor='#FFFFFF',axisbelow=True))
for axi in axes:
    for ax in axi:
        ax.grid(color='white', linestyle='solid')                       
plt.subplots_adjust(wspace=0.1)
plt.subplots_adjust(hspace=0.15)
legendalpha=0.7

#Renewables part
plot_renewables(axes,segment_wise=True,more_labels=True,edgeless=True)

#Fossil part
plot_fossils(axes,segment_wise=True,more_labels=True,edgeless=True)

#Intersections
#plot_intersections(axes)

#Format labels
for axi in axes:
    for ax in axi:
        ax.set_xlim(0,max(np.array(eroei_el_ren).flatten()))
        ax.set_ylim(0,max(np.array(eroei_el_ren).flatten()))
axes[0][0].set_ylabel('$EROEI_{CCS}$ | $EROEI_{disp}$',fontsize=14,y=0)
axes[1][1].set_xlabel('$EROEI_{el}$',fontsize=14,x=-0.05)        
axes[0][1].set_yticklabels('')   
axes[1][1].set_yticklabels('')   
axes[0][0].set_xticklabels('')   
axes[0][1].set_xticklabels('')  
axes[0][0].legend(loc=(0.03,0.52),framealpha=0,fontsize=11)       
        
#plt.suptitle('Comparison of adjusted EROEIs for \nfossil fuels with carbon-capture and storage $EROEI_{CCS}$ and \ndispatchable renewables with energy storage $EROEI_{disp}$',y=1.03,fontsize=13)
plt.savefig('plot/extra/fig4_edgeless.png',dpi=150,bbox_inches = 'tight', pad_inches = 0.1)
plt.show()
plotting Wind area...
plotting PV area...
plotting NGCC Gas area...
plotting Pulverized Coal area...
plotting IGCC Coal area...