Merge pull request #296 from PyPSA/carbon-management

carbon management

.gitignore

@@ -26,6 +26,9 @@ gurobi.log
 /data/Industrial_Database.csv
 /data/retro/tabula-calculator-calcsetbuilding.csv
 /data/nuts*
+data/gas_network/scigrid-gas/
+dask-worker-space/
+publications.jrc.ec.europa.eu/
 
 *.org
 
@@ -48,3 +51,7 @@ doc/_build
 *.xls
 
 *.geojson
+
+*.ipynb
+
+data/costs_*

Snakefile

@@ -288,6 +288,23 @@ else:
     build_biomass_transport_costs_output = {}
 
 
+if config["sector"]["regional_co2_sequestration_potential"]["enable"]:
+    rule build_sequestration_potentials:
+        input:
+            sequestration_potential=HTTP.remote("https://raw.githubusercontent.com/ericzhou571/Co2Storage/main/resources/complete_map_2020_unit_Mt.geojson", keep_local=True),
+            regions_onshore=pypsaeur("resources/regions_onshore_elec_s{simpl}_{clusters}.geojson"),
+            regions_offshore=pypsaeur("resources/regions_offshore_elec_s{simpl}_{clusters}.geojson"),
+        output:
+            sequestration_potential="resources/co2_sequestration_potential_elec_s{simpl}_{clusters}.csv"
+        threads: 1
+        resources: mem_mb=4000
+        benchmark: "benchmarks/build_sequestration_potentials_s{simpl}_{clusters}"
+        script: "scripts/build_sequestration_potentials.py"
+    build_sequestration_potentials_output = rules.build_sequestration_potentials.output
+else:
+    build_sequestration_potentials_output = {}
+
+
 rule build_salt_cavern_potentials:
     input:
         salt_caverns="data/h2_salt_caverns_GWh_per_sqkm.geojson",
@@ -520,7 +537,8 @@ rule prepare_sector_network:
         solar_thermal_rural="resources/solar_thermal_rural_elec_s{simpl}_{clusters}.nc" if config["sector"]["solar_thermal"] else [],
         **build_retro_cost_output,
         **build_biomass_transport_costs_output,
-        **gas_infrastructure
+        **gas_infrastructure,
+        **build_sequestration_potentials_output
     output: RDIR + '/prenetworks/elec_s{simpl}_{clusters}_lv{lv}_{opts}_{sector_opts}_{planning_horizons}.nc'
     threads: 1
     resources: mem_mb=2000

config.default.yaml

@@ -250,10 +250,19 @@ sector:
   coal_cc: false
   dac: true
   co2_vent: false
+  allam_cycle: false
   SMR: true
+  regional_co2_sequestration_potential:
+    enable: false  # enable regionally resolved geological co2 storage potential
+    attribute: 'conservative estimate Mt'
+    include_onshore: false  # include onshore sequestration potentials
+    min_size: 3 # Gt, sites with lower potential will be excluded
+    max_size: 25 # Gt, max sequestration potential for any one site, TODO research suitable value
+    years_of_storage: 25 # years until potential exhausted at optimised annual rate
   co2_sequestration_potential: 200  #MtCO2/a sequestration potential for Europe
   co2_sequestration_cost: 10   #EUR/tCO2 for sequestration of CO2
-  co2_network: false
+  co2_spatial: false
+  co2network: false
   cc_fraction: 0.9  # default fraction of CO2 captured with post-combustion capture
   hydrogen_underground_storage: true
   hydrogen_underground_storage_locations:
@@ -279,7 +288,8 @@ sector:
   gas_network_connectivity_upgrade: 1 # https://networkx.org/documentation/stable/reference/algorithms/generated/networkx.algorithms.connectivity.edge_augmentation.k_edge_augmentation.html#networkx.algorithms.connectivity.edge_augmentation.k_edge_augmentation
   gas_distribution_grid: true
   gas_distribution_grid_cost_factor: 1.0  #multiplies cost in data/costs.csv
-  biomass_transport: false  # biomass transport between nodes
+  biomass_spatial: false  # regionally resolve biomass (e.g. potentials)
+  biomass_transport: false  # allow transport of solid biomass between nodes
   conventional_generation: # generator : carrier
     OCGT: gas
   biomass_to_liquid: false

doc/release_notes.rst

@@ -97,6 +97,30 @@ incorporates retrofitting options to hydrogen.
   as explicit ICE shares for land transport (``land_transport_ice_share``) and
   agriculture machinery (``agriculture_machinery_oil_share``).
 
+* Add option to spatially resolve carrier representing stored carbon dioxide
+  (``co2_spatial``). This allows for more detailed modelling of CCUTS, e.g.
+  regarding the capturing of industrial process emissions, usage as feedstock
+  for electrofuels, transport of carbon dioxide, and geological sequestration sites.
+
+* Add option for planning a new carbon dioxide network (``co2network``).
+
+
+* Add option for regionally-resolved geological carbon dioxide sequestration
+  potentials through new rule ``build_sequestration_potentials`` based on
+  `CO2StoP <https://setis.ec.europa.eu/european-co2-storage-database_en>`_. This
+  can be controlled in the section ``regional_co2_sequestration_potential`` of
+  the ``config.yaml``. It includes options to select the level of conservatism,
+  whether onshore potentials should be included, the respective upper and lower
+  limits per region, and an annualisation parameter for the cumulative
+  potential. The defaults are preliminary and will be validated the next
+  release.
+
+* Separate option to regionally resolve biomass (``biomass_spatial``) from
+  option to allow biomass transport (``biomass_transport``).
+
+* Add option to include `Allam cycle gas power plants
+  <https://en.wikipedia.org/wiki/Allam_power_cycle>`_ (``allam_cycle``).
+
 **Bugfixes**
 
 * The CO2 sequestration limit implemented as GlobalConstraint (introduced in the previous version)

scripts/build_sequestration_potentials.py

@@ -0,0 +1,43 @@
+import pandas as pd
+import geopandas as gpd
+
+def area(gdf):
+    """Returns area of GeoDataFrame geometries in square kilometers."""
+    return gdf.to_crs(epsg=3035).area.div(1e6)
+
+
+def allocate_sequestration_potential(gdf, regions, attr='conservative estimate Mt', threshold=3):
+    gdf = gdf.loc[gdf[attr] > threshold, [attr, "geometry"]]
+    gdf["area_sqkm"] = area(gdf)
+    overlay = gpd.overlay(regions, gdf, keep_geom_type=True)
+    overlay["share"] = area(overlay) / overlay["area_sqkm"]
+    adjust_cols = overlay.columns.difference({"name", "area_sqkm", "geometry", "share"})
+    overlay[adjust_cols] = overlay[adjust_cols].multiply(overlay["share"], axis=0)
+    gdf_regions = overlay.groupby("name").sum()
+    gdf_regions.drop(["area_sqkm", "share"], axis=1, inplace=True)
+    return gdf_regions.squeeze()
+
+
+if __name__ == "__main__":
+    if 'snakemake' not in globals():
+        from helper import mock_snakemake
+        snakemake = mock_snakemake(
+            'build_sequestration_potentials',
+            simpl='',
+            clusters="181"
+        )
+
+    cf = snakemake.config["sector"]["regional_co2_sequestration_potential"]
+
+    gdf = gpd.read_file(snakemake.input.sequestration_potential[0])
+
+    regions = gpd.read_file(snakemake.input.regions_offshore)
+    if cf["include_onshore"]:
+        onregions = gpd.read_file(snakemake.input.regions_onshore)
+        regions = pd.concat([regions, onregions]).dissolve(by='name').reset_index()
+
+    s = allocate_sequestration_potential(gdf, regions, attr=cf["attribute"], threshold=cf["min_size"])
+
+    s = s.where(s>cf["min_size"]).dropna()
+
+    s.to_csv(snakemake.output.sequestration_potential)

scripts/helper.py

@@ -138,6 +138,6 @@ def parse(l):
 
 def update_config_with_sector_opts(config, sector_opts):
     for o in sector_opts.split("-"):
-        if o.startswith("CF:"):
+        if o.startswith("CF+"):
             l = o.split("+")[1:]
             update_config(config, parse(l))

scripts/make_summary.py

@@ -273,7 +273,7 @@ def calculate_supply(n, label, supply):
 
             for end in [col[3:] for col in c.df.columns if col[:3] == "bus"]:
 
-                items = c.df.index[c.df["bus" + end].map(bus_map, na_action=None)]
+                items = c.df.index[c.df["bus" + end].map(bus_map).fillna(False)]
 
                 if len(items) == 0:
                     continue
@@ -318,7 +318,7 @@ def calculate_supply_energy(n, label, supply_energy):
 
             for end in [col[3:] for col in c.df.columns if col[:3] == "bus"]:
 
-                items = c.df.index[c.df["bus" + str(end)].map(bus_map, na_action=None)]
+                items = c.df.index[c.df["bus" + str(end)].map(bus_map).fillna(False)]
 
                 if len(items) == 0:
                     continue

scripts/plot_network.py

@@ -310,7 +310,7 @@ def plot_h2_map(network, regions):
 
     else:
 
-        h2_total = h2_new
+        h2_total = h2_new.p_nom_opt
 
     link_widths_total = h2_total / linewidth_factor
 

scripts/prepare_sector_network.py

@@ -48,7 +48,7 @@ def define_spatial(nodes, options):
 
     spatial.biomass = SimpleNamespace()
 
-    if options["biomass_transport"]:
+    if options.get("biomass_spatial", options["biomass_transport"]):
         spatial.biomass.nodes = nodes + " solid biomass"
         spatial.biomass.locations = nodes
         spatial.biomass.industry = nodes + " solid biomass for industry"
@@ -65,14 +65,16 @@ def define_spatial(nodes, options):
 
     spatial.co2 = SimpleNamespace()
 
-    if options["co2_network"]:
+    if options["co2_spatial"]:
         spatial.co2.nodes = nodes + " co2 stored"
         spatial.co2.locations = nodes
         spatial.co2.vents = nodes + " co2 vent"
+        spatial.co2.process_emissions = nodes + " process emissions"
     else:
         spatial.co2.nodes = ["co2 stored"]
         spatial.co2.locations = ["EU"]
         spatial.co2.vents = ["co2 vent"]
+        spatial.co2.process_emissions = ["process emissions"]
 
     spatial.co2.df = pd.DataFrame(vars(spatial.co2), index=nodes)
 
@@ -92,8 +94,11 @@ def define_spatial(nodes, options):
         spatial.gas.locations = ["EU"]
         spatial.gas.biogas = ["EU biogas"]
         spatial.gas.industry = ["gas for industry"]
-        spatial.gas.industry_cc = ["gas for industry CC"]
         spatial.gas.biogas_to_gas = ["EU biogas to gas"]
+        if options.get("co2_spatial", options["co2network"]):
+            spatial.gas.industry_cc = nodes + " gas for industry CC"
+        else:
+            spatial.gas.industry_cc = ["gas for industry CC"]
 
     spatial.gas.df = pd.DataFrame(vars(spatial.gas), index=nodes)
 
@@ -432,6 +437,7 @@ def add_carrier_buses(n, carrier, nodes=None):
         e_nom_extendable=True,
         e_cyclic=True,
         carrier=carrier,
+        capital_cost=0.2 * costs.at[carrier, "discount rate"] # preliminary value to avoid zeros
     )
 
     n.madd("Generator",
@@ -511,10 +517,17 @@ def add_co2_tracking(n, options):
         unit="t_co2"
     )
 
+    if options["regional_co2_sequestration_potential"]["enable"]:
+        upper_limit = options["regional_co2_sequestration_potential"]["max_size"] * 1e3 # Mt
+        annualiser = options["regional_co2_sequestration_potential"]["years_of_storage"]
+        e_nom_max = pd.read_csv(snakemake.input.sequestration_potential, index_col=0).squeeze()
+        e_nom_max = e_nom_max.reindex(spatial.co2.locations).fillna(0.).clip(upper=upper_limit).mul(1e6) / annualiser # t
+        e_nom_max = e_nom_max.rename(index=lambda x: x + " co2 stored")
+
     n.madd("Store",
         spatial.co2.nodes,
         e_nom_extendable=True,
-        e_nom_max=np.inf,
+        e_nom_max=e_nom_max,
         capital_cost=options['co2_sequestration_cost'],
         carrier="co2 stored",
         bus=spatial.co2.nodes
@@ -554,6 +567,29 @@ def add_co2_network(n, costs):
     )
 
 
+def add_allam(n, costs):  
+
+    logger.info("Adding Allam cycle gas power plants.")  
+
+    nodes = pop_layout.index  
+
+    n.madd("Link",  
+        nodes,
+        suffix=" allam",  
+        bus0=spatial.gas.df.loc[nodes, "nodes"].values,  
+        bus1=nodes,  
+        bus2=spatial.co2.df.loc[nodes, "nodes"].values,  
+        carrier="allam",  
+        p_nom_extendable=True,  
+        # TODO: add costs to technology-data
+        capital_cost=0.6*1.5e6*0.1, # efficiency * EUR/MW * annuity  
+        marginal_cost=2,  
+        efficiency=0.6,  
+        efficiency2=costs.at['gas', 'CO2 intensity'],  
+        lifetime=30.,
+    )
+
+
 def add_dac(n, costs):
 
     heat_carriers = ["urban central heat", "services urban decentral heat"]
@@ -1222,7 +1258,7 @@ def add_storage_and_grids(n, costs):
             bus0=spatial.coal.nodes,
             bus1=spatial.nodes,
             bus2="co2 atmosphere",
-            bus3="co2 stored",
+            bus3=spatial.co2.nodes,
             marginal_cost=costs.at['coal', 'efficiency'] * costs.at['coal', 'VOM'], #NB: VOM is per MWel
             capital_cost=costs.at['coal', 'efficiency'] * costs.at['coal', 'fixed'] + costs.at['biomass CHP capture', 'fixed'] * costs.at['coal', 'CO2 intensity'], #NB: fixed cost is per MWel
             p_nom_extendable=True,
@@ -1832,7 +1868,7 @@ def add_biomass(n, costs):
     else:
         biogas_potentials_spatial = biomass_potentials["biogas"].sum()
 
-    if options["biomass_transport"]:
+    if options.get("biomass_spatial", options["biomass_transport"]):
         solid_biomass_potentials_spatial = biomass_potentials["solid biomass"].rename(index=lambda x: x + " solid biomass")
     else:
         solid_biomass_potentials_spatial = biomass_potentials["solid biomass"].sum()
@@ -1904,10 +1940,10 @@ def add_biomass(n, costs):
             biomass_transport.index,
             bus0=biomass_transport.bus0 + " solid biomass",
             bus1=biomass_transport.bus1 + " solid biomass",
-            p_nom_extendable=True,
+            p_nom_extendable=False,
+            p_nom=5e4,
             length=biomass_transport.length.values,
             marginal_cost=biomass_transport.costs * biomass_transport.length.values,
-            capital_cost=1,
             carrier="solid biomass transport"
         )
 
@@ -2056,7 +2092,7 @@ def add_industry(n, costs):
         unit="MWh_LHV"
     )
 
-    if options["biomass_transport"]:
+    if options.get("biomass_spatial", options["biomass_transport"]):
         p_set = industrial_demand.loc[spatial.biomass.locations, "solid biomass"].rename(index=lambda x: x + " solid biomass for industry") / 8760
     else:
         p_set = industrial_demand["solid biomass"].sum() / 8760
@@ -2402,25 +2438,31 @@ def add_industry(n, costs):
         p_set=industrial_demand.loc[nodes, "electricity"] / 8760
     )
 
-    n.add("Bus",
-        "process emissions",
-        location="EU",
+    n.madd("Bus",
+        spatial.co2.process_emissions,
+        location=spatial.co2.locations,
         carrier="process emissions",
         unit="t_co2"
     )
 
+    sel = ["process emission", "process emission from feedstock"]
+    if options["co2_spatial"] or options["co2network"]:
+        p_set = -industrial_demand.loc[nodes, sel].sum(axis=1).rename(index=lambda x: x + " process emissions") / 8760
+    else:
+        p_set = -industrial_demand.loc[nodes, sel].sum(axis=1).sum() / 8760
+
     # this should be process emissions fossil+feedstock
     # then need load on atmosphere for feedstock emissions that are currently going to atmosphere via Link Fischer-Tropsch demand
-    n.add("Load",
-        "process emissions",
-        bus="process emissions",
+    n.madd("Load",
+        spatial.co2.process_emissions,
+        bus=spatial.co2.process_emissions,
         carrier="process emissions",
-        p_set=-industrial_demand.loc[nodes,["process emission", "process emission from feedstock"]].sum(axis=1).sum() / 8760
+        p_set=p_set,
     )
 
-    n.add("Link",
-        "process emissions",
-        bus0="process emissions",
+    n.madd("Link",
+        spatial.co2.process_emissions,
+        bus0=spatial.co2.process_emissions,
         bus1="co2 atmosphere",
         carrier="process emissions",
         p_nom_extendable=True,
@@ -2431,7 +2473,7 @@ def add_industry(n, costs):
     n.madd("Link",
         spatial.co2.locations,
         suffix=" process emissions CC",
-        bus0="process emissions",
+        bus0=spatial.co2.process_emissions,
         bus1="co2 atmosphere",
         bus2=spatial.co2.nodes,
         carrier="process emissions CC",
@@ -2879,9 +2921,12 @@ def set_temporal_aggregation(n, opts, solver_name):
     if "noH2network" in opts:
         remove_h2_network(n)
 
-    if options["co2_network"]:
+    if options["co2network"]:
         add_co2_network(n, costs)
 
+    if options["allam_cycle"]:
+        add_allam(n, costs)
+
     solver_name = snakemake.config["solving"]["solver"]["name"]
     n = set_temporal_aggregation(n, opts, solver_name)