changed implementation to always use 2020 cost

envs/environment.yaml

@@ -20,7 +20,7 @@ dependencies:
 - openpyxl!=3.1.1
 - pycountry
 - seaborn
-- snakemake-minimal>=8.5,<8.6
+- snakemake-minimal>=8.5
 - memory_profiler
 - yaml
 - pytables

rules/build_sector.smk

@@ -248,8 +248,8 @@ rule build_solar_thermal_profiles:
     output:
         solar_thermal=resources("solar_thermal_{scope}_elec_s{simpl}_{clusters}.nc"),
     resources:
-        mem_mb=20000,
+        mem_mb=8000,
-    threads: 16
+    threads: 1
     log:
         logs("build_solar_thermal_profiles_{scope}_s{simpl}_{clusters}.log"),
     benchmark:
@@ -816,6 +816,34 @@ def input_profile_offwind(w):
     }
 
 
+rule build_egs_potentials:
+    params:
+        snapshots=config_provider("snapshots"),
+        sector=config_provider("sector"),
+        costs=config_provider("costs"),
+    input:
+        egs_cost="data/egs_costs.json",
+        regions=resources("regions_onshore_elec_s{simpl}_{clusters}.geojson"),
+        air_temperature=(
+            resources("temp_air_total_elec_s{simpl}_{clusters}.nc")
+            if config_provider("sector", "enhanced_geothermal", "var_cf")
+            else []
+        ),
+    output:
+        egs_potentials=resources("egs_potentials_s{simpl}_{clusters}.csv"),
+        egs_overlap=resources("egs_overlap_s{simpl}_{clusters}.csv"),
+        egs_capacity_factors=resources("egs_capacity_factors_s{simpl}_{clusters}.csv"),
+    threads: 1
+    resources:
+        mem_mb=2000,
+    log:
+        logs("build_egs_potentials_s{simpl}_{clusters}.log"),
+    conda:
+        "../envs/environment.yaml"
+    script:
+        "../scripts/build_egs_potentials.py"
+
+
 rule prepare_sector_network:
     params:
         time_resolution=config_provider("clustering", "temporal", "resolution_sector"),
@@ -931,6 +959,21 @@ rule prepare_sector_network:
             if config_provider("sector", "solar_thermal")(w)
             else []
         ),
+        egs_potentials=lambda w: (
+            resources("egs_potentials_s{simpl}_{clusters}.csv")
+            if config_provider("sector", "enhanced_geothermal", "enable")(w)
+            else []
+        ),
+        egs_overlap=lambda w: (
+            resources("egs_overlap_s{simpl}_{clusters}.csv")
+            if config_provider("sector", "enhanced_geothermal", "enable")(w)
+            else []
+        ),
+        egs_capacity_factors=lambda w: (
+            resources("egs_capacity_factors_s{simpl}_{clusters}.csv")
+            if config_provider("sector", "enhanced_geothermal", "enable")(w)
+            else []
+        ),
     output:
         RESULTS
         + "prenetworks/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}.nc",

scripts/build_egs_potentials.py

@@ -30,6 +30,9 @@ from shapely.geometry import Polygon
 
 
 def prepare_egs_data(egs_file):
+    """
+    Processes the original .json file EGS data to a more human-readable format.
+    """
     with open(egs_file) as f:
         jsondata = json.load(f)
 
@@ -86,6 +89,65 @@ def prepare_egs_data(egs_file):
     return egs_data
 
 
+def prepare_capex(prepared_data):
+    """
+    The source paper provides only data for year and regions where LCOE <
+    100Euro/MWh. However, this implementations starts with the costs for 2020
+    for all regions and then adjusts the costs according to the user's chosen
+    setting in the config file.
+
+    As such, for regions where cost data is available only from, say,
+    2035, we need to reverse-engineer the costs for 2020. This is done
+    in the following (unfortunately verbose) function.
+    """
+
+    default_year = 2020
+
+    # obtains all available CAPEX data
+    capex_df = pd.DataFrame(columns=prepared_data.keys())
+
+    for year in capex_df.columns:
+
+        year_data = prepared_data[year].groupby("geometry").mean().reset_index()
+
+        for g in year_data.geometry:
+
+            if not g in year_data.geometry.tolist():
+                # weird but apparently necessary
+                continue
+
+            capex_df.loc[g, year] = year_data.loc[
+                year_data.geometry == g, "CAPEX"
+            ].values[0]
+
+    capex_df = capex_df.loc[:, default_year:]
+
+    # fill up missing values assuming cost reduction factors similar to existing values
+    for sooner, later in zip(capex_df.columns[::-1][1:], capex_df.columns[::-1]):
+
+        missing_mask = capex_df[sooner].isna()
+        cr_factor = (
+            capex_df.loc[~missing_mask, later] / capex_df.loc[~missing_mask, sooner]
+        )
+
+        capex_df.loc[missing_mask, sooner] = (
+            capex_df.loc[missing_mask, later] / cr_factor.mean()
+        )
+
+    # harmonice capacity and CAPEX
+    p_nom_max = prepared_data[2050].groupby("geometry")["PowerSust"].mean()
+    p_nom_max = p_nom_max.loc[p_nom_max > 0]
+
+    capex_df = capex_df.loc[p_nom_max.index]
+
+    data = (
+        pd.concat((capex_df[default_year], p_nom_max), axis=1)
+        .reset_index()
+        .rename(columns={2020: "CAPEX"})
+    )
+    return gpd.GeoDataFrame(data, geometry=data.geometry)
+
+
 def get_capacity_factors(network_regions_file, air_temperatures_file):
     """
     Performance of EGS is higher for lower temperatures, due to more efficient
@@ -146,28 +208,9 @@ if __name__ == "__main__":
     egs_config = snakemake.params["sector"]["enhanced_geothermal"]
     costs_config = snakemake.params["costs"]
 
-    sustainability_factor = egs_config["sustainability_factor"]
-    # the share of heat that is replenished from the earth's core.
-    # we are not constraining ourselves to the sustainable share, but
-    # inversely apply it to our underlying data, which refers to the
-    # sustainable heat. Source: Relative magnitude of sustainable heat vs
-    # nonsustainable heat in the paper "From hot rock to useful energy..."
-
     egs_data = prepare_egs_data(snakemake.input.egs_cost)
+    egs_data = prepare_capex(egs_data)
 
-    if egs_config["optimism"]:
-        egs_data = egs_data[(year := costs_config["year"])]
-        logger.info(
-            f"EGS optimism! Building EGS potentials with costs estimated for {year}."
-        )
-
-    else:
-        egs_data = egs_data[(default_year := 2020)]
-        logger.info(
-            f"No EGS optimism! Building EGS potentials with {default_year} costs."
-        )
-
-    egs_data = egs_data.loc[egs_data["PowerSust"] > 0].reset_index(drop=True)
     egs_regions = egs_data.geometry
 
     network_regions = (
@@ -188,7 +231,12 @@ if __name__ == "__main__":
 
     overlap_matrix.to_csv(snakemake.output["egs_overlap"])
 
-    # consider not only replenished heat
+    # the share of heat that is replenished from the earth's core.
+    # we are not constraining ourselves to the sustainable share, but
+    # inversely apply it to our underlying data, which refers to the
+    # sustainable heat. Source: Relative magnitude of sustainable heat vs
+    # nonsustainable heat in the paper "From hot rock to useful energy..."
+    sustainability_factor = egs_config["sustainability_factor"]
     egs_data["p_nom_max"] = egs_data["PowerSust"] / sustainability_factor
 
     egs_data[["p_nom_max", "CAPEX"]].to_csv(snakemake.output["egs_potentials"])

scripts/prepare_sector_network.py

@@ -190,7 +190,7 @@ def define_spatial(nodes, options):
 
     # deep geothermal
     spatial.geothermal_heat = SimpleNamespace()
-    spatial.geothermal_heat.nodes = ["EU deep geothermal"]
+    spatial.geothermal_heat.nodes = ["EU enhanced geothermal systems"]
     spatial.geothermal_heat.locations = ["EU"]
 
     return spatial
@@ -3599,15 +3599,11 @@ def add_enhanced_geothermal(
 
     # under egs optimism, the expected cost reductions also cover costs for ORC
     # hence, the ORC costs are no longer taken from technology-data
-    orc_capex = (
+    orc_capex = costs.at["organic rankine cycle", "investment"]
-        costs.at["organic rankine cycle", "investment"]
-        if not snakemake.params.sector["enhanced_geothermal_optimism"]
-        else 0.0
-    )
 
     # cost for ORC is subtracted, as it is already included in the geothermal cost.
     # The orc cost are attributed to a separate link representing the ORC.
-    # also capital_cost conversion Eurofalse/kW -> Euro/MW
+    # also capital_cost conversion Euro/kW -> Euro/MW
 
     egs_potentials["capital_cost"] = (
         (egs_annuity + FOM / (1.0 + FOM))
@@ -3619,8 +3615,8 @@ def add_enhanced_geothermal(
         egs_potentials["capital_cost"] > 0
     ).all(), "Error in EGS cost, negative values found."
 
-    plant_annuity = calculate_annuity(costs.at["organic rankine cycle", "lifetime"], dr)
+    orc_annuity = calculate_annuity(costs.at["organic rankine cycle", "lifetime"], dr)
-    plant_capital_cost = (plant_annuity + FOM / (1 + FOM)) * orc_capex * Nyears
+    orc_capital_cost = (orc_annuity + FOM / (1 + FOM)) * orc_capex * Nyears
 
     efficiency_orc = costs.at["organic rankine cycle", "efficiency"]
     efficiency_dh = costs.at["geothermal", "efficiency residential heat"]
@@ -3647,19 +3643,22 @@ def add_enhanced_geothermal(
         p_nom_extendable=True,
     )
 
-    if snakemake.params.sector["enhanced_geothermal_var_cf"]:
+    if egs_config["var_cf"]:
         efficiency = pd.read_csv(
             snakemake.input.egs_capacity_factors, parse_dates=True, index_col=0
         )
         logger.info("Adding Enhanced Geothermal with time-varying capacity factors.")
     else:
-        efficiency = pd.Series(1, overlap.index)
+        efficiency = 1.0
 
     # if urban central heat exists, adds geothermal as CHP
     as_chp = "urban central heat" in n.loads.carrier.unique()
 
     if as_chp:
-        logger.info("Adding Enhanced Geothermal as Combined Heat and Power.")
+        logger.info("Adding EGS as Combined Heat and Power.")
+
+    else:
+        logger.info("Adding EGS for Electricity Only.")
 
     for bus, bus_overlap in overlap.iterrows():
         if not bus_overlap.sum():
@@ -3674,34 +3673,36 @@ def add_enhanced_geothermal(
         bus_egs["p_nom_max"] = bus_egs["p_nom_max"].multiply(bus_overlap)
         bus_egs = bus_egs.loc[bus_egs.p_nom_max > 0.0]
 
-        if egs_config["performant"]:
+        appendix = " " + pd.Index(np.arange(len(bus_egs)).astype(str))
-            bus_egs = bus_egs.sort_values(by="capital_cost").iloc[:1]
-            appendix = pd.Index([""])
-        else:
-            appendix = " " + pd.Index(np.arange(len(bus_egs)).astype(str))
 
         # add surface bus
-        n.add(
+        n.madd(
             "Bus",
-            f"{bus} geothermal heat surface",
+            pd.Index([f"{bus} geothermal heat surface"]),
+            location=bus,
+            unit="MWh_th",
+            carrier="geothermal heat",
         )
 
         bus_egs.index = np.arange(len(bus_egs)).astype(str)
         well_name = f"{bus} enhanced geothermal" + appendix
 
-        bus_eta = pd.concat(
+        if egs_config["var_cf"]:
-            (efficiency[bus].rename(idx) for idx in well_name),
+            bus_eta = pd.concat(
-            axis=1,
+                (efficiency[bus].rename(idx) for idx in well_name),
-        )
+                axis=1,
+            )
+        else:
+            bus_eta = efficiency
 
         p_nom_max = bus_egs["p_nom_max"]
         capital_cost = bus_egs["capital_cost"]
         bus1 = pd.Series(f"{bus} geothermal heat surface", well_name)
 
+        # adding geothermal wells as multiple generators to represent supply curve
         n.madd(
             "Link",
             well_name,
-            location=bus,
             bus0=spatial.geothermal_heat.nodes,
             bus1=bus1,
             carrier="geothermal heat",
@@ -3711,6 +3712,7 @@ def add_enhanced_geothermal(
             efficiency=bus_eta,
         )
 
+        # adding Organic Rankine Cycle as a single link
         n.add(
             "Link",
             bus + " geothermal organic rankine cycle",
@@ -3718,7 +3720,7 @@ def add_enhanced_geothermal(
             bus1=bus,
             p_nom_extendable=True,
             carrier="geothermal organic rankine cycle",
-            capital_cost=plant_capital_cost * efficiency_orc,
+            capital_cost=orc_capital_cost * efficiency_orc,
             efficiency=efficiency_orc,
         )
 
@@ -3729,7 +3731,7 @@ def add_enhanced_geothermal(
                 bus0=f"{bus} geothermal heat surface",
                 bus1=bus + " urban central heat",
                 carrier="geothermal district heat",
-                capital_cost=plant_capital_cost
+                capital_cost=orc_capital_cost
                 * efficiency_orc
                 * costs.at["geothermal", "district heating cost"],
                 efficiency=efficiency_dh,
@@ -3745,8 +3747,8 @@ def add_enhanced_geothermal(
             # Hence, it is counter-intuitive to install it at the surface bus,
             # this is however the more lean and computationally efficient solution.
 
-            max_hours = egs_config["reservoir_max_hours"]
+            max_hours = egs_config["max_hours"]
-            boost = egs_config["reservoir_max_boost"]
+            boost = egs_config["max_boost"]
 
             n.add(
                 "StorageUnit",
@@ -3756,6 +3758,7 @@ def add_enhanced_geothermal(
                 p_nom_extendable=True,
                 p_min_pu=-boost,
                 max_hours=max_hours,
+                cyclic_state_of_charge=True,
             )
 
 
@@ -3884,6 +3887,12 @@ if __name__ == "__main__":
     if options["electricity_distribution_grid"]:
         insert_electricity_distribution_grid(n, costs)
 
+    if options["enhanced_geothermal"].get("enable", False):
+        logger.info("Adding Enhanced Geothermal Systems (EGS).")
+        add_enhanced_geothermal(
+            n, snakemake.input["egs_potentials"], snakemake.input["egs_overlap"], costs
+        )
+
     maybe_adjust_costs_and_potentials(n, snakemake.params["adjustments"])
 
     if options["gas_distribution_grid"]:
@@ -3911,12 +3920,6 @@ if __name__ == "__main__":
     if options.get("cluster_heat_buses", False) and not first_year_myopic:
         cluster_heat_buses(n)
 
-    if options["enhanced_geothermal"].get("enable", False):
-        logger.info("Adding Enhanced Geothermal Potential.")
-        add_enhanced_geothermal(
-            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)

scripts/solve_network.py

@@ -801,7 +801,7 @@ def add_geothermal_chp_constraint(n):
     elec_index = n.links.loc[
         n.links.carrier == "geothermal organic rankine cycle"
     ].index
-    heat_index = n.links.loc[n.links.carrier == "geothermal heat district heat"].index
+    heat_index = n.links.loc[n.links.carrier == "geothermal district heat"].index
 
     p_nom_lhs = (
         n.model["Link-p_nom"].loc[heat_index] - n.model["Link-p_nom"].loc[elec_index]
@@ -824,7 +824,7 @@ def add_flexible_egs_constraint(n):
 
     n.model.add_constraints(
         p_nom_lhs <= p_nom_rhs,
-        name="Upper bounds the charging capacity of the storage unit",
+        name="Upper bounds the charging capacity of the geothermal reservoir according to the well capacity",
     )
 
 
@@ -954,28 +954,6 @@ def solve_network(n, config, solving, **kwargs):
         n.model.print_infeasibilities()
         raise RuntimeError("Solving status 'infeasible'")
 
-    # check if enhanced_geothermal_performant might have changed model results
-    if (
-        snakemake.config["sector"]["enhanced_geothermal"]
-        and snakemake.config["sector"]["enhanced_geothermal_performant"]
-    ):
-        mask = (mask := n.links.carrier == "geothermal heat") & (
-            n.links.loc[mask, "p_nom_max"] > 0.0
-        )
-
-        saturated = n.links.loc[mask, "p_nom_max"] == n.links.loc[mask, "p_nom_opt"]
-
-        if saturated.any():
-            logger.warning(
-                (
-                    "Potential for enhanced geothermal heat is saturated at bus(es):\n"
-                    f"{', '.join(n.links.loc[saturated.loc[saturated.astype(bool)].index, 'location'].tolist())}.\n"
-                    "Consider setting config['sector']['enhanced_geothermal_performant'] to False."
-                )
-            )
-
-    return n
-
 
 if __name__ == "__main__":
     if "snakemake" not in globals():