Merge pull request #733 from LukasFrankenQ/master

Enhanced Geothermal Systems

config/config.default.yaml

@@ -621,6 +621,13 @@ sector:
       solar: 3
       offwind-ac: 3
       offwind-dc: 3
+  enhanced_geothermal:
+    enable: false
+    flexible: true
+    max_hours: 240
+    max_boost: 0.25
+    var_cf: true
+    sustainability_factor: 0.0025
 
 # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#industry
 industry:
@@ -1184,6 +1191,9 @@ plotting:
     waste: '#e3d37d'
     other: '#000000'
     geothermal: '#ba91b1'
+    geothermal heat: '#ba91b1'
+    geothermal district heat: '#d19D00'
+    geothermal organic rankine cycle: '#ffbf00'
     AC: "#70af1d"
     AC-AC: "#70af1d"
     AC line: "#70af1d"

doc/configtables/sector.csv

@@ -145,3 +145,11 @@ limit_max_growth,,,
 -- -- {carrier},GW,float,The historic maximum growth of a carrier
 -- max_relative_growth,,,
 -- -- {carrier},p.u.,float,The historic maximum relative growth of a carrier
+,,,
+enhanced_geothermal,,,
+-- enable,--,"{true, false}",Add option to include Enhanced Geothermal Systems
+-- flexible,--,"{true, false}",Add option for flexible operation (see Ricks et al. 2024)
+-- max_hours,--,int,The maximum hours the reservoir can be charged under flexible operation
+-- max_boost,--,float,The maximum boost in power output under flexible operation
+-- var_cf,--,"{true, false}",Add option for variable capacity factor (see Ricks et al. 2024)
+-- sustainability_factor,--,float,Share of sourced heat that is replenished by the earth's core (see details in `build_egs_potentials.py <https://github.com/PyPSA/pypsa-eur-sec/blob/master/scripts/build_egs_potentials.py>`_)

doc/release_notes.rst

@@ -7,8 +7,15 @@
 Release Notes
 ##########################################
 
-.. Upcoming Release
+Upcoming Release
-.. ================
+================
+
+* Added Enhanced Geothermal Systems for generation of electricity and district heat.
+  Cost and available capacity assumptions based on `Aghahosseini et al. (2020)
+  <https://www.sciencedirect.com/science/article/pii/S0306261920312551>`__.
+  See configuration ``sector: enhanced_geothermal`` for details; by default switched off.
+
+
 PyPSA-Eur 0.11.0 (25th May 2024)
 =====================================
 
@@ -808,7 +815,7 @@ PyPSA-Eur 0.9.0 (5th January 2024)
 
 * The minimum PyPSA version is now 0.26.1.
 
-* Update to ``tsam>=0.2.3`` for performance improvents in temporal clustering.
+* Update to ``tsam>=0.2.3`` for performance improvements in temporal clustering.
 
 * Pin ``snakemake`` version to below 8.0.0, as the new version is not yet
   supported. The next release will switch to the requirement ``snakemake>=8``.

rules/build_sector.smk

@@ -902,6 +902,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"),
@@ -1022,6 +1050,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

@@ -0,0 +1,249 @@
+# -*- coding: utf-8 -*-
+# SPDX-FileCopyrightText: : 2023 @LukasFranken, The PyPSA-Eur Authors
+#
+# SPDX-License-Identifier: MIT
+"""
+This rule extracts potential and cost for electricity generation through
+enhanced geothermal systems.
+
+For this, we use data from "From hot rock to useful energy..." by Aghahosseini, Breyer (2020)
+'https://doi.org/10.1016/j.apenergy.2020.115769'
+Note that we input data used here is not the same as in the paper, but was passed on by the authors.
+
+The data provides a lon-lat gridded map of Europe (1° x 1°), with each grid cell assigned
+a heat potential (in GWh) and a cost (in EUR/MW).
+
+This scripts overlays that map with the network's regions, and builds a csv with CAPEX, OPEX and p_nom_max
+"""
+
+import logging
+
+logger = logging.getLogger(__name__)
+
+import json
+
+import geopandas as gpd
+import numpy as np
+import pandas as pd
+import xarray as xr
+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)
+
+    def point_to_square(p, lon_extent=1.0, lat_extent=1.0):
+        try:
+            x, y = p.coords.xy[0][0], p.coords.xy[1][0]
+        except IndexError:
+            return p
+
+        return Polygon(
+            [
+                [x - lon_extent / 2, y - lat_extent / 2],
+                [x - lon_extent / 2, y + lat_extent / 2],
+                [x + lon_extent / 2, y + lat_extent / 2],
+                [x + lon_extent / 2, y - lat_extent / 2],
+            ]
+        )
+
+    years = [2015, 2020, 2025, 2030, 2035, 2040, 2045, 2050]
+    lcoes = ["LCOE50", "LCOE100", "LCOE150"]
+
+    egs_data = dict()
+
+    for year in years:
+        df = pd.DataFrame(columns=["Lon", "Lat", "CAPEX", "HeatSust", "PowerSust"])
+
+        for lcoe in lcoes:
+            for country_data in jsondata[lcoe]:
+                try:
+                    country_df = pd.DataFrame(
+                        columns=df.columns,
+                        index=range(len(country_data[0][years.index(year)]["Lon"])),
+                    )
+                except TypeError:
+                    country_df = pd.DataFrame(columns=df.columns, index=range(0))
+
+                for col in df.columns:
+                    country_df[col] = country_data[0][years.index(year)][col]
+
+                if country_df.dropna().empty:
+                    continue
+                elif df.empty:
+                    df = country_df.dropna()
+                else:
+                    df = pd.concat((df, country_df.dropna()), ignore_index=True)
+
+        gdf = gpd.GeoDataFrame(
+            df.drop(columns=["Lon", "Lat"]), geometry=gpd.points_from_xy(df.Lon, df.Lat)
+        ).reset_index(drop=True)
+
+        gdf["geometry"] = gdf.geometry.apply(lambda geom: point_to_square(geom))
+        egs_data[year] = gdf
+
+    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
+    air cooling Data from Ricks et al.: The Role of Flexible Geothermal Power
+    in Decarbonized Elec Systems.
+    """
+
+    # these values are taken from the paper's
+    # Supplementary Figure 20 from https://zenodo.org/records/7093330
+    # and relate deviations of the ambient temperature from the year-average
+    # ambient temperature to EGS capacity factors.
+    delta_T = [-15, -10, -5, 0, 5, 10, 15, 20]
+    cf = [1.17, 1.13, 1.07, 1, 0.925, 0.84, 0.75, 0.65]
+
+    x = np.linspace(-15, 20, 200)
+    y = np.interp(x, delta_T, cf)
+
+    upper_x = np.linspace(20, 25, 50)
+    m_upper = (y[-1] - y[-2]) / (x[-1] - x[-2])
+    upper_y = upper_x * m_upper - x[-1] * m_upper + y[-1]
+
+    lower_x = np.linspace(-20, -15, 50)
+    m_lower = (y[1] - y[0]) / (x[1] - x[0])
+    lower_y = lower_x * m_lower - x[0] * m_lower + y[0]
+
+    x = np.hstack((lower_x, x, upper_x))
+    y = np.hstack((lower_y, y, upper_y))
+
+    network_regions = gpd.read_file(network_regions_file).set_crs(epsg=4326)
+    index = network_regions["name"]
+
+    air_temp = xr.open_dataset(air_temperatures_file)
+
+    snapshots = pd.date_range(freq="h", **snakemake.params.snapshots)
+    capacity_factors = pd.DataFrame(index=snapshots)
+
+    # bespoke computation of capacity factors for each bus.
+    # Considering the respective temperatures, we compute
+    # the deviation from the average temperature and relate it
+    # to capacity factors based on the data from above.
+    for bus in index:
+        temp = air_temp.sel(name=bus).to_dataframe()["temperature"]
+        capacity_factors[bus] = np.interp((temp - temp.mean()).values, x, y)
+
+    return capacity_factors
+
+
+if __name__ == "__main__":
+    if "snakemake" not in globals():
+        from _helpers import mock_snakemake
+
+        snakemake = mock_snakemake(
+            "build_egs_potentials",
+            simpl="",
+            clusters=37,
+        )
+
+    egs_config = snakemake.params["sector"]["enhanced_geothermal"]
+    costs_config = snakemake.params["costs"]
+
+    egs_data = prepare_egs_data(snakemake.input.egs_cost)
+    egs_data = prepare_capex(egs_data)
+
+    egs_regions = egs_data.geometry
+
+    network_regions = (
+        gpd.read_file(snakemake.input.regions)
+        .set_index("name", drop=True)
+        .set_crs(epsg=4326)
+    )
+
+    overlap_matrix = pd.DataFrame(
+        index=network_regions.index,
+        columns=egs_data.index,
+    )
+
+    for name, polygon in network_regions.geometry.items():
+        overlap_matrix.loc[name] = (
+            egs_regions.intersection(polygon).area
+        ) / egs_regions.area
+
+    overlap_matrix.to_csv(snakemake.output["egs_overlap"])
+
+    # 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"])
+
+    capacity_factors = get_capacity_factors(
+        snakemake.input["regions"],
+        snakemake.input["air_temperature"],
+    )
+
+    capacity_factors.to_csv(snakemake.output["egs_capacity_factors"])

scripts/prepare_sector_network.py

@@ -196,6 +196,11 @@ def define_spatial(nodes, options):
     spatial.lignite.nodes = ["EU lignite"]
     spatial.lignite.locations = ["EU"]
 
+    # deep geothermal
+    spatial.geothermal_heat = SimpleNamespace()
+    spatial.geothermal_heat.nodes = ["EU enhanced geothermal systems"]
+    spatial.geothermal_heat.locations = ["EU"]
+
     return spatial
 
 
@@ -976,7 +981,7 @@ def insert_electricity_distribution_grid(n, costs):
         .get("efficiency_static")
     ):
         logger.info(
-            f"Deducting distribution losses from electricity demand: {100*(1-efficiency)}%"
+            f"Deducting distribution losses from electricity demand: {np.around(100*(1-efficiency), decimals=2)}%"
         )
         n.loads_t.p_set.loc[:, n.loads.carrier == "electricity"] *= efficiency
 
@@ -3726,6 +3731,210 @@ def lossy_bidirectional_links(n, carrier, efficiencies={}):
         )
 
 
+def add_enhanced_geothermal(n, egs_potentials, egs_overlap, costs):
+    """
+    Adds EGS potential to model.
+
+    Built in scripts/build_egs_potentials.py
+    """
+
+    if len(spatial.geothermal_heat.nodes) > 1:
+        logger.warning(
+            "'add_enhanced_geothermal' not implemented for multiple geothermal nodes."
+        )
+    logger.info(
+        "[EGS] implemented with 2020 CAPEX from Aghahosseini et al 2021: 'From hot rock to...'."
+    )
+    logger.info(
+        "[EGS] Recommended usage scales CAPEX to future cost expectations using config 'adjustments'."
+    )
+    logger.info("[EGS] During this the relevant carriers are:")
+    logger.info("[EGS] drilling part -> 'geothermal heat'")
+    logger.info(
+        "[EGS] electricity generation part -> 'geothermal organic rankine cycle'"
+    )
+    logger.info("[EGS] district heat distribution part -> 'geothermal district heat'")
+
+    egs_config = snakemake.params["sector"]["enhanced_geothermal"]
+    costs_config = snakemake.config["costs"]
+
+    # matrix defining the overlap between gridded geothermal potential estimation, and bus regions
+    overlap = pd.read_csv(egs_overlap, index_col=0)
+    overlap.columns = overlap.columns.astype(int)
+    egs_potentials = pd.read_csv(egs_potentials, index_col=0)
+
+    Nyears = n.snapshot_weightings.generators.sum() / 8760
+    dr = costs_config["fill_values"]["discount rate"]
+    lt = costs.at["geothermal", "lifetime"]
+    FOM = costs.at["geothermal", "FOM"]
+
+    egs_annuity = calculate_annuity(lt, dr)
+
+    # 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 = costs.at["organic rankine cycle", "investment"]
+
+    # 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 Euro/kW -> Euro/MW
+
+    egs_potentials["capital_cost"] = (
+        (egs_annuity + FOM / (1.0 + FOM))
+        * (egs_potentials["CAPEX"] * 1e3 - orc_capex)
+        * Nyears
+    )
+
+    assert (
+        egs_potentials["capital_cost"] > 0
+    ).all(), "Error in EGS cost, negative values found."
+
+    orc_annuity = calculate_annuity(costs.at["organic rankine cycle", "lifetime"], dr)
+    orc_capital_cost = (orc_annuity + FOM / (1 + FOM)) * orc_capex * Nyears
+
+    efficiency_orc = costs.at["organic rankine cycle", "efficiency"]
+    efficiency_dh = costs.at["geothermal", "district heat-input"]
+
+    # p_nom_max conversion GW -> MW
+    egs_potentials["p_nom_max"] = egs_potentials["p_nom_max"] * 1000.0
+
+    # not using add_carrier_buses, as we are not interested in a Store
+    n.add("Carrier", "geothermal heat")
+
+    n.madd(
+        "Bus",
+        spatial.geothermal_heat.nodes,
+        carrier="geothermal heat",
+        unit="MWh_th",
+    )
+
+    n.madd(
+        "Generator",
+        spatial.geothermal_heat.nodes,
+        bus=spatial.geothermal_heat.nodes,
+        carrier="geothermal heat",
+        p_nom_extendable=True,
+    )
+
+    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 = 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 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():
+            continue
+
+        overlap = bus_overlap.loc[bus_overlap > 0.0]
+        bus_egs = egs_potentials.loc[overlap.index]
+
+        if not len(bus_egs):
+            continue
+
+        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]
+
+        appendix = " " + pd.Index(np.arange(len(bus_egs)).astype(str))
+
+        # add surface bus
+        n.madd(
+            "Bus",
+            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
+
+        if egs_config["var_cf"]:
+            bus_eta = pd.concat(
+                (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,
+            bus0=spatial.geothermal_heat.nodes,
+            bus1=bus1,
+            carrier="geothermal heat",
+            p_nom_extendable=True,
+            p_nom_max=p_nom_max.set_axis(well_name) / efficiency_orc,
+            capital_cost=capital_cost.set_axis(well_name) * efficiency_orc,
+            efficiency=bus_eta,
+        )
+
+        # adding Organic Rankine Cycle as a single link
+        n.add(
+            "Link",
+            bus + " geothermal organic rankine cycle",
+            bus0=f"{bus} geothermal heat surface",
+            bus1=bus,
+            p_nom_extendable=True,
+            carrier="geothermal organic rankine cycle",
+            capital_cost=orc_capital_cost * efficiency_orc,
+            efficiency=efficiency_orc,
+        )
+
+        if as_chp and bus + " urban central heat" in n.buses.index:
+            n.add(
+                "Link",
+                bus + " geothermal heat district heat",
+                bus0=f"{bus} geothermal heat surface",
+                bus1=bus + " urban central heat",
+                carrier="geothermal district heat",
+                capital_cost=orc_capital_cost
+                * efficiency_orc
+                * costs.at["geothermal", "district heat surcharge"]
+                / 100.0,
+                efficiency=efficiency_dh,
+                p_nom_extendable=True,
+            )
+        elif as_chp and not bus + " urban central heat" in n.buses.index:
+            n.links.at[bus + " geothermal organic rankine cycle", "efficiency"] = (
+                efficiency_orc
+            )
+
+        if egs_config["flexible"]:
+            # this StorageUnit represents flexible operation using the geothermal reservoir.
+            # 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["max_hours"]
+            boost = egs_config["max_boost"]
+
+            n.add(
+                "StorageUnit",
+                bus + " geothermal reservoir",
+                bus=f"{bus} geothermal heat surface",
+                carrier="geothermal heat",
+                p_nom_extendable=True,
+                p_min_pu=-boost,
+                max_hours=max_hours,
+                cyclic_state_of_charge=True,
+            )
+
+
 # %%
 if __name__ == "__main__":
     if "snakemake" not in globals():
@@ -3857,6 +4066,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"]:

scripts/solve_network.py

@@ -948,6 +948,25 @@ def add_pipe_retrofit_constraint(n):
     n.model.add_constraints(lhs == rhs, name="Link-pipe_retrofit")
 
 
+def add_flexible_egs_constraint(n):
+    """
+    Upper bounds the charging capacity of the geothermal reservoir according to
+    the well capacity.
+    """
+    well_index = n.links.loc[n.links.carrier == "geothermal heat"].index
+    storage_index = n.storage_units.loc[
+        n.storage_units.carrier == "geothermal heat"
+    ].index
+
+    p_nom_rhs = n.model["Link-p_nom"].loc[well_index]
+    p_nom_lhs = n.model["StorageUnit-p_nom"].loc[storage_index]
+
+    n.model.add_constraints(
+        p_nom_lhs <= p_nom_rhs,
+        name="upper_bound_charging_capacity_of_geothermal_reservoir",
+    )
+
+
 def add_co2_atmosphere_constraint(n, snapshots):
     glcs = n.global_constraints[n.global_constraints.type == "co2_atmosphere"]
 
@@ -1013,6 +1032,9 @@ def extra_functionality(n, snapshots):
     else:
         add_co2_atmosphere_constraint(n, snapshots)
 
+    if config["sector"]["enhanced_geothermal"]["enable"]:
+        add_flexible_egs_constraint(n)
+
     if snakemake.params.custom_extra_functionality:
         source_path = snakemake.params.custom_extra_functionality
         assert os.path.exists(source_path), f"{source_path} does not exist"