carbon management

Snakefile

@@ -280,6 +280,23 @@ else:
     build_biomass_transport_costs_output = {}
 
 
+if config["sector"].get("sequestration_potential", False):
+    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",
@@ -512,7 +529,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,9 +250,11 @@ sector:
   dac: true
   co2_vent: false
   SMR: true
+  sequestration_potential: true  # geological co2 storage potential
   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:
@@ -275,6 +277,7 @@ 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_spatial: false  # biomass transport between nodes
   biomass_transport: false  # biomass transport between nodes
   conventional_generation: # generator : carrier
     OCGT: gas

scripts/build_sequestration_potentials.py

@@ -0,0 +1,48 @@
+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"
+        )
+
+    # TODO move to config.yaml
+    threshold = 3
+    include_onshore = False
+
+    gdf = gpd.read_file(snakemake.input.sequestration_potential)
+
+    regions = gpd.read_file(snakemake.input.regions_offshore)
+    if include_onshore:
+        onregions = gpd.read_file(snakemake.input.regions_onshore)
+        regions = pd.concat([regions, onregions]).dissolve(by='name').reset_index()
+
+    attr = snakemake.config['sector']["sequestration_potential"]
+    kwargs = dict(attr=attr, threshold=threshold) if isinstance(attr, str) else {}
+
+    s = allocate_sequestration_potential(gdf, regions, **kwargs)
+
+    s = s.where(s>threshold).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/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

@@ -44,7 +44,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"
@@ -61,7 +61,7 @@ def define_spatial(nodes, options):
 
     spatial.co2 = SimpleNamespace()
 
-    if options["co2_network"]:
+    if options.get("co2_spatial", options["co2network"]):
         spatial.co2.nodes = nodes + " co2 stored"
         spatial.co2.locations = nodes
         spatial.co2.vents = nodes + " co2 vent"
@@ -88,8 +88,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)
 
@@ -507,10 +510,18 @@ def add_co2_tracking(n, options):
         unit="t_co2"
     )
 
+    if options["sequestration_potential"]:
+        # TODO make configurable options
+        upper_limit = 25000 # Mt
+        annualiser = 25 # TODO research suitable value
+        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
@@ -1218,7 +1229,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,
@@ -1828,7 +1839,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()
@@ -1900,10 +1911,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"
         )
 
@@ -2052,7 +2063,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
@@ -2775,7 +2786,7 @@ 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)
 
     solver_name = snakemake.config["solving"]["solver"]["name"]