pre-testrun version added to prepare_sector_network

rules/build_sector.smk

@@ -783,6 +783,8 @@ rule prepare_sector_network:
         cop_air_total=RESOURCES + "cop_air_total_elec_s{simpl}_{clusters}.nc",
         cop_air_rural=RESOURCES + "cop_air_rural_elec_s{simpl}_{clusters}.nc",
         cop_air_urban=RESOURCES + "cop_air_urban_elec_s{simpl}_{clusters}.nc",
+        egs_potentials=RESOURCES + "egs_potentials_s{simpl}_{clusters}.csv",
+        egs_overlap=RESOURCES + "egs_overlap_s{simpl}_{clusters}.csv",
         solar_thermal_total=RESOURCES
         + "solar_thermal_total_elec_s{simpl}_{clusters}.nc"
         if config["sector"]["solar_thermal"]

scripts/build_egs_potentials.py

@@ -23,7 +23,7 @@ import json
 import pandas as pd
 import geopandas as gpd
 
-from shapely.geometry import Polygon, Point
+from shapely.geometry import Polygon
 
 
 def prepare_egs_data(egs_file):
@@ -108,13 +108,23 @@ if __name__ == "__main__":
     egs_data.index = egs_data.geometry.astype(str)
     egs_shapes = egs_data.geometry
     
-    network_shapes = gpd.read_file(snakemake.input.shapes).set_index("name", drop=True).set_crs(epsg=4326)
-    egs_shapes = egs_data.geometry
+    network_shapes = (
+        gpd.read_file(snakemake.input.shapes)
+        .set_index("name", drop=True)
+        .set_crs(epsg=4326)
+    )
 
-    overlap_matrix = pd.DataFrame(index=network_shapes.index, columns=egs_shapes.index)
+    overlap_matrix = (
+        pd.DataFrame(
+            index=network_shapes.index,
+            columns=(egs_shapes := egs_data.geometry).astype(str).values)
+    )
 
     for name, polygon in network_shapes.geometry.items():
-        overlap_matrix.loc[name] = egs_shapes.intersection(polygon).area / egs_shapes.area
+        overlap_matrix.loc[name] = (
+            egs_shapes
+            .intersection(polygon).area
+        ) / egs_shapes.area
 
     overlap_matrix.to_csv(snakemake.output["egs_overlap"])
 

scripts/prepare_sector_network.py

@@ -3280,6 +3280,121 @@ def set_temporal_aggregation(n, opts, solver_name):
     return n
 
 
+def add_egs_potential(
+    n,
+    egs_potentials,
+    egs_overlap,
+    costs,
+):
+    """
+    Adds EGS potential to model.
+
+    Built in scripts/build_egs_potentials.py
+    """
+
+    n.add(
+        "Carrier",
+        "geothermal heat",
+        nice_name="Geothermal Heat",
+        color=snakemake.config["plotting"]["tech_colors"]["geothermal"],
+        co2_emissions=costs.at["geothermal", "CO2 intensity"],
+    )
+
+    config = snakemake.config
+
+    overlap = pd.read_csv(egs_overlap, index_col=0)
+    indicator_matrix = (
+        pd.DataFrame(
+            np.ceil(overlap.values),
+            index=overlap.index,
+            columns=overlap.columns
+        )
+    )
+
+    egs_potentials = pd.read_csv(egs_potentials, index_col=0)
+    
+    Nyears = n.snapshot_weightings.generators.sum() / 8760
+    dr = config["costs"]["fill_values"]["discount rate"]
+    lt = costs.at["geothermal", "lifetime"]
+
+    egs_annuity = calculate_annuity(lt, dr)
+
+    egs_potentials["capital_cost"] = (
+        (egs_annuity + 0.02 / 1.02) * egs_potentials["CAPEX"] * Nyears 
+    )
+
+    efficiency = costs.at["geothermal", "efficiency electricity"]
+
+    # p_nom conversion GW -> MW
+    # capital_cost conversion Euro/kW -> Euro/MW
+    egs_potentials["p_nom_max"] = egs_potentials["p_nom_max"] * 1000.0
+    egs_potentials["capital_cost"] = egs_potentials["capital_cost"] * 1000.0
+
+    n.add(
+        "Bus",
+        "EU geothermal heat",
+        carrier="geothermal heat",
+        unit="MWh_th",
+        location="EU",
+    )
+
+    n.add(
+        "Generator",
+        "EU geothermal heat",
+        bus="EU geothermal heat",
+        carrier="geothermal heat",
+        p_nom_max=np.inf,
+        p_nom_extendable=True,
+    )
+
+    for bus, overlaps in overlap.iterrows():
+        if not overlaps.sum():
+            continue
+
+        overlaps.index = np.arange(len(overlaps))
+    
+        overlaps = overlaps.loc[overlaps > 0.]
+        indicators = indicator_matrix.loc[bus]
+        indicators.index = np.arange(len(indicators))
+
+        bus_egs = egs_potentials.loc[overlaps.loc[overlaps > 0.].index]
+
+        if not len(bus_egs):
+            continue
+        
+        bus_egs["p_nom_max"] = bus_egs["p_nom_max"].multiply(overlaps)
+        bus_egs = bus_egs.loc[bus_egs.p_nom_max > 0.]
+
+        if config["sector"]["enhanced_geothermal_best_only"]:
+            bus_egs = bus_egs.sort_values(by="capital_cost").iloc[:1]
+            appendix = ""
+        else:
+            appendix = " " + bus_egs.index
+
+        bus_egs.index = np.arange(len(bus_egs)).astype(str)
+
+        p_nom_max = bus_egs["p_nom_max"]
+        p_nom_max.index = f"{bus} enhanced geothermal" + appendix
+
+        capital_cost = bus_egs["capital_cost"]
+        capital_cost.index = f"{bus} enhanced geothermal" + appendix
+
+        n.madd(
+            "Link",
+            f"{bus} enhanced geothermal" + appendix,
+            bus0="EU geothermal heat",
+            bus1=bus,
+            bus2="co2 atmosphere",
+            carrier="geothermal heat",
+            p_nom_extendable=True,
+            p_nom_max=p_nom_max / efficiency,
+            capital_cost=capital_cost * efficiency,
+            efficiency=efficiency,
+            efficiency2=costs.at["geothermal", "CO2 intensity"] * efficiency,
+        )
+
+
+
 if __name__ == "__main__":
     if "snakemake" not in globals():
         from _helpers import mock_snakemake
@@ -3453,6 +3568,17 @@ if __name__ == "__main__":
     if options.get("cluster_heat_buses", False) and not first_year_myopic:
         cluster_heat_buses(n)
 
+    if options.get("enhanced_geothermal"):
+        logger.info("Adding Enhanced Geothermal Potential.")
+
+        add_egs_potential(
+            n,
+            snakemake.input["egs_potentials"],
+            snakemake.input["egs_overlap"],
+            costs
+        )
+
+
     n.meta = dict(snakemake.config, **dict(wildcards=dict(snakemake.wildcards)))
 
     sanitize_carriers(n, snakemake.config)