From 71e2a4943cfa0a04ad094f70fe656840b20c4e22 Mon Sep 17 00:00:00 2001
From: Fabian Neumann <fabian.neumann@outlook.de>
Date: Tue, 24 Jan 2023 18:44:39 +0100
Subject: [PATCH 01/11] carbon management

---
 Snakefile                                 | 20 +++++++++-
 config.default.yaml                       |  5 ++-
 scripts/build_sequestration_potentials.py | 48 +++++++++++++++++++++++
 scripts/helper.py                         |  2 +-
 scripts/prepare_sector_network.py         | 27 +++++++++----
 5 files changed, 91 insertions(+), 11 deletions(-)
 create mode 100644 scripts/build_sequestration_potentials.py

diff --git a/Snakefile b/Snakefile
index 90fcfb56..cf829322 100644
--- a/Snakefile
+++ b/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="data/complete_map_2020_unit_Mt.geojson",
+            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
diff --git a/config.default.yaml b/config.default.yaml
index ee1c5059..5fa072d4 100644
--- a/config.default.yaml
+++ b/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
diff --git a/scripts/build_sequestration_potentials.py b/scripts/build_sequestration_potentials.py
new file mode 100644
index 00000000..cf8ab6d6
--- /dev/null
+++ b/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)
diff --git a/scripts/helper.py b/scripts/helper.py
index e6ddfd4a..62ae33c0 100644
--- a/scripts/helper.py
+++ b/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))
\ No newline at end of file
diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py
index b6c052be..76b71fc8 100644
--- a/scripts/prepare_sector_network.py
+++ b/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()
@@ -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"]

From 9411cad98b9430fd0ce628915fa1ee401c711a1f Mon Sep 17 00:00:00 2001
From: Fabian Neumann <fabian.neumann@outlook.de>
Date: Tue, 24 Jan 2023 18:47:34 +0100
Subject: [PATCH 02/11] add sequestration potential data link

---
 Snakefile | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Snakefile b/Snakefile
index cf829322..c4680c81 100644
--- a/Snakefile
+++ b/Snakefile
@@ -283,7 +283,7 @@ else:
 if config["sector"].get("sequestration_potential", False):
     rule build_sequestration_potentials:
         input:
-            sequestration_potential="data/complete_map_2020_unit_Mt.geojson",
+            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:

From 697be2d9e3f7cccd9300ebaa089ba77c898d9a46 Mon Sep 17 00:00:00 2001
From: Fabian Neumann <fabian.neumann@outlook.de>
Date: Wed, 25 Jan 2023 08:41:08 +0100
Subject: [PATCH 03/11] fix capacities of biomass transport

The biomass transport links are split into either direction
because they have a marginal cost associated with the
transport. With p_nom_extendable, many decision variables
are created. This commit tries another way where
capacities are non-extendable and set to sufficiently high
values.
---
 scripts/prepare_sector_network.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py
index 76b71fc8..7fcadacb 100644
--- a/scripts/prepare_sector_network.py
+++ b/scripts/prepare_sector_network.py
@@ -1911,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"
         )
 

From 0681cf93ab4641ea081a5b112927e48afc4e2c15 Mon Sep 17 00:00:00 2001
From: Fabian Neumann <fabian.neumann@outlook.de>
Date: Wed, 25 Jan 2023 09:08:51 +0100
Subject: [PATCH 04/11] fix hydrogen network plotting without retrofitting

---
 scripts/plot_network.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/scripts/plot_network.py b/scripts/plot_network.py
index 53c52187..3d0641d5 100644
--- a/scripts/plot_network.py
+++ b/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
 

From a0b18a8fcdb24ef4cae7ded90e140e233cb10794 Mon Sep 17 00:00:00 2001
From: Fabian Neumann <fabian.neumann@outlook.de>
Date: Mon, 6 Feb 2023 09:46:24 +0100
Subject: [PATCH 05/11] manage process emissions in spatial

---
 scripts/prepare_sector_network.py | 16 +++++++++-------
 1 file changed, 9 insertions(+), 7 deletions(-)

diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py
index 7fcadacb..26d2b2d3 100644
--- a/scripts/prepare_sector_network.py
+++ b/scripts/prepare_sector_network.py
@@ -65,10 +65,12 @@ def define_spatial(nodes, options):
         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)
 
@@ -2408,8 +2410,8 @@ def add_industry(n, costs):
     )
 
     n.add("Bus",
-        "process emissions",
-        location="EU",
+        spatial.co2.process_emissions,
+        location=spatial.co2.locations,
         carrier="process emissions",
         unit="t_co2"
     )
@@ -2417,15 +2419,15 @@ def add_industry(n, costs):
     # 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",
+        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
     )
 
     n.add("Link",
-        "process emissions",
-        bus0="process emissions",
+        spatial.co2.process_emissions,
+        bus0=spatial.co2.process_emissions,
         bus1="co2 atmosphere",
         carrier="process emissions",
         p_nom_extendable=True,
@@ -2436,7 +2438,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",

From ad8ccf3b25d72b3635cc6668b75572319bac539c Mon Sep 17 00:00:00 2001
From: Fabian Neumann <fabian.neumann@outlook.de>
Date: Mon, 6 Feb 2023 09:52:35 +0100
Subject: [PATCH 06/11] fix process emission localisation with co2 (not)
 spatially resolved

---
 scripts/prepare_sector_network.py | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py
index 26d2b2d3..7a7edc25 100644
--- a/scripts/prepare_sector_network.py
+++ b/scripts/prepare_sector_network.py
@@ -2416,13 +2416,19 @@ def add_industry(n, costs):
         unit="t_co2"
     )
 
+    sel = ["process emission", "process emission from feedstock"]
+    if options.get("co2_spatial", 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",
         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",

From fcf527d946506b90db1829799f1e73e1a4748326 Mon Sep 17 00:00:00 2001
From: Fabian <fab.hof@gmx.de>
Date: Tue, 7 Feb 2023 09:50:44 +0100
Subject: [PATCH 07/11] update gitignore

---
 .gitignore | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/.gitignore b/.gitignore
index a7049efe..9371d303 100644
--- a/.gitignore
+++ b/.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
 

From 94fb6c646cdd06fbc3bc5f9276c7968e5bb647bb Mon Sep 17 00:00:00 2001
From: Fabian <fab.hof@gmx.de>
Date: Wed, 8 Feb 2023 22:57:01 +0100
Subject: [PATCH 08/11] prepare sector: add allam cycle build sequestration:
 fix input for new snakemake versions

---
 scripts/build_sequestration_potentials.py |  2 +-
 scripts/make_summary.py                   |  4 +--
 scripts/prepare_sector_network.py         | 34 +++++++++++++++++++----
 3 files changed, 32 insertions(+), 8 deletions(-)

diff --git a/scripts/build_sequestration_potentials.py b/scripts/build_sequestration_potentials.py
index cf8ab6d6..3b06f665 100644
--- a/scripts/build_sequestration_potentials.py
+++ b/scripts/build_sequestration_potentials.py
@@ -31,7 +31,7 @@ def allocate_sequestration_potential(gdf, regions, attr='conservative estimate M
     threshold = 3
     include_onshore = False
 
-    gdf = gpd.read_file(snakemake.input.sequestration_potential)
+    gdf = gpd.read_file(snakemake.input.sequestration_potential[0])
 
     regions = gpd.read_file(snakemake.input.regions_offshore)
     if include_onshore:
diff --git a/scripts/make_summary.py b/scripts/make_summary.py
index f9c74c89..06680cb4 100644
--- a/scripts/make_summary.py
+++ b/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
diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py
index 7a7edc25..4fcdbf84 100644
--- a/scripts/prepare_sector_network.py
+++ b/scripts/prepare_sector_network.py
@@ -61,7 +61,7 @@ def define_spatial(nodes, options):
 
     spatial.co2 = SimpleNamespace()
 
-    if options.get("co2_spatial", options["co2network"]):
+    if options["co2_spatial"]:
         spatial.co2.nodes = nodes + " co2 stored"
         spatial.co2.locations = nodes
         spatial.co2.vents = nodes + " co2 vent"
@@ -433,6 +433,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",
@@ -563,6 +564,28 @@ def add_co2_network(n, costs):
     )
 
 
+def add_allam(n, costs):  
+  
+    logger.info("Adding Allam cycle generators")  
+    
+    nodes = pop_layout.index  
+    
+    n.madd("Link",  
+    nodes + " 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"]
@@ -2409,7 +2432,7 @@ def add_industry(n, costs):
         p_set=industrial_demand.loc[nodes, "electricity"] / 8760
     )
 
-    n.add("Bus",
+    n.madd("Bus",
         spatial.co2.process_emissions,
         location=spatial.co2.locations,
         carrier="process emissions",
@@ -2417,21 +2440,21 @@ def add_industry(n, costs):
     )
 
     sel = ["process emission", "process emission from feedstock"]
-    if options.get("co2_spatial", options["co2network"]):
+    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",
+    n.madd("Load",
         spatial.co2.process_emissions,
         bus=spatial.co2.process_emissions,
         carrier="process emissions",
         p_set=p_set,
     )
 
-    n.add("Link",
+    n.madd("Link",
         spatial.co2.process_emissions,
         bus0=spatial.co2.process_emissions,
         bus1="co2 atmosphere",
@@ -2796,6 +2819,7 @@ def set_temporal_aggregation(n, opts, solver_name):
 
     if options["co2network"]:
         add_co2_network(n, costs)
+        add_allam(n, costs)
 
     solver_name = snakemake.config["solving"]["solver"]["name"]
     n = set_temporal_aggregation(n, opts, solver_name)

From ccde2e7f931ee1a6a74a10ccb0b32f80039a90b9 Mon Sep 17 00:00:00 2001
From: Fabian Neumann <fabian.neumann@outlook.de>
Date: Thu, 16 Feb 2023 17:21:58 +0100
Subject: [PATCH 09/11] make sequestration options configurable

---
 .gitignore                                |  4 +++
 Snakefile                                 |  2 +-
 config.default.yaml                       | 13 ++++++--
 scripts/build_sequestration_potentials.py | 13 +++-----
 scripts/prepare_sector_network.py         | 38 ++++++++++++-----------
 5 files changed, 39 insertions(+), 31 deletions(-)

diff --git a/.gitignore b/.gitignore
index 9371d303..e3f625c3 100644
--- a/.gitignore
+++ b/.gitignore
@@ -51,3 +51,7 @@ doc/_build
 *.xls
 
 *.geojson
+
+*.ipynb
+
+data/costs_*
\ No newline at end of file
diff --git a/Snakefile b/Snakefile
index ef1ed840..b14adacf 100644
--- a/Snakefile
+++ b/Snakefile
@@ -288,7 +288,7 @@ else:
     build_biomass_transport_costs_output = {}
 
 
-if config["sector"].get("sequestration_potential", False):
+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),
diff --git a/config.default.yaml b/config.default.yaml
index d09ff950..c1172bd1 100644
--- a/config.default.yaml
+++ b/config.default.yaml
@@ -250,8 +250,15 @@ sector:
   coal_cc: false
   dac: true
   co2_vent: false
+  allam_cycle: false
   SMR: true
-  sequestration_potential: true  # geological co2 storage potential
+  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_spatial: false
@@ -281,8 +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_spatial: false  # biomass transport between nodes
-  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
diff --git a/scripts/build_sequestration_potentials.py b/scripts/build_sequestration_potentials.py
index 3b06f665..4983640b 100644
--- a/scripts/build_sequestration_potentials.py
+++ b/scripts/build_sequestration_potentials.py
@@ -27,22 +27,17 @@ def allocate_sequestration_potential(gdf, regions, attr='conservative estimate M
             clusters="181"
         )
 
-    # TODO move to config.yaml
-    threshold = 3
-    include_onshore = False
+    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 include_onshore:
+    if cf["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, attr=cf["attribute"], threshold=cf["min_size"])
 
-    s = allocate_sequestration_potential(gdf, regions, **kwargs)
-
-    s = s.where(s>threshold).dropna()
+    s = s.where(s>cf["min_size"]).dropna()
 
     s.to_csv(snakemake.output.sequestration_potential)
diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py
index 5348a79e..a47c5710 100644
--- a/scripts/prepare_sector_network.py
+++ b/scripts/prepare_sector_network.py
@@ -517,10 +517,9 @@ 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
+    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")
@@ -570,24 +569,25 @@ def add_co2_network(n, costs):
 
 def add_allam(n, costs):  
   
-    logger.info("Adding Allam cycle generators")  
+    logger.info("Adding Allam cycle gas power plants.")  
     
     nodes = pop_layout.index  
     
     n.madd("Link",  
-    nodes + " 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.,
-    )  
+        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):
@@ -2923,6 +2923,8 @@ def set_temporal_aggregation(n, opts, solver_name):
 
     if options["co2network"]:
         add_co2_network(n, costs)
+
+    if options["allam_cycle"]:
         add_allam(n, costs)
 
     solver_name = snakemake.config["solving"]["solver"]["name"]

From b2ea0576cdbe7b034b20ce6b819a72bfde12a1e4 Mon Sep 17 00:00:00 2001
From: Fabian Neumann <fabian.neumann@outlook.de>
Date: Thu, 16 Feb 2023 17:32:36 +0100
Subject: [PATCH 10/11] add release note

---
 doc/release_notes.rst | 24 ++++++++++++++++++++++++
 1 file changed, 24 insertions(+)

diff --git a/doc/release_notes.rst b/doc/release_notes.rst
index 210de55a..4700f27b 100644
--- a/doc/release_notes.rst
+++ b/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)

From 126a5be729516d4c3e61dbf48d147346e7c91969 Mon Sep 17 00:00:00 2001
From: Fabian Neumann <fabian.neumann@outlook.de>
Date: Thu, 16 Feb 2023 17:34:13 +0100
Subject: [PATCH 11/11] add new options to test config files

---
 test/config.myopic.yaml    | 2 ++
 test/config.overnight.yaml | 2 ++
 2 files changed, 4 insertions(+)

diff --git a/test/config.myopic.yaml b/test/config.myopic.yaml
index 3719921e..255d6734 100644
--- a/test/config.myopic.yaml
+++ b/test/config.myopic.yaml
@@ -235,6 +235,8 @@ sector:
   dac: true
   co2_vent: true
   SMR: true
+  regional_co2_sequestration_potential:
+    enable: false
   co2_sequestration_potential: 200  #MtCO2/a sequestration potential for Europe
   co2_sequestration_cost: 10   #EUR/tCO2 for sequestration of CO2
   co2_network: false
diff --git a/test/config.overnight.yaml b/test/config.overnight.yaml
index 6663d393..a68ad540 100644
--- a/test/config.overnight.yaml
+++ b/test/config.overnight.yaml
@@ -233,6 +233,8 @@ sector:
   dac: true
   co2_vent: true
   SMR: true
+  regional_co2_sequestration_potential:
+    enable: false
   co2_sequestration_potential: 200  #MtCO2/a sequestration potential for Europe
   co2_sequestration_cost: 10   #EUR/tCO2 for sequestration of CO2
   co2_network: false