Add geothermal-sourced central heat pumps

.gitignore

@@ -75,3 +75,5 @@ d1gam3xoknrgr2.cloudfront.net/
 merger-todos.md
 
 *.html
+# private dev folder
+dev/*

config/config.default.yaml

@@ -480,6 +480,18 @@ sector:
       heat_exchanger_pinch_point_temperature_difference: 5 #K
       isentropic_compressor_efficiency: 0.8
       heat_loss: 0.0
+    heat_utilisation_potentials:
+      geothermal:
+        # activate for 85C hydrothermal
+        # key: hydrothermal_85
+        # constant_temperature_celsius: 85
+        key: hydrothermal_65
+        constant_temperature_celsius: 65
+        column_name: Energy_TWh
+        unit: TWh
+        full_load_hours: 4000
+    direct_utilisation_heat_sources:
+    - geothermal
   heat_pump_sources:
     urban central:
     - air
@@ -827,7 +839,7 @@ industry:
 # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#costs
 costs:
   year: 2030
-  version: v0.9.2
+  version: v0.10.0
   social_discountrate: 0.02
   fill_values:
     FOM: 0
@@ -1228,8 +1240,11 @@ plotting:
     services urban decentral air heat pump: '#5af95d'
     services rural air heat pump: '#5af95d'
     urban central air heat pump: '#6cfb6b'
+    urban central geothermal heat pump: '#4f2144'
+    geothermal heat pump: '#4f2144'
+    geothermal heat direct utilisation: '#ba91b1'
     ground heat pump: '#2fb537'
-    residential rural ground heat pump: '#48f74f'
+    residential rural ground heat pump: '#4f2144'
     residential rural air heat pump: '#48f74f'
     services rural ground heat pump: '#5af95d'
     Ambient: '#98eb9d'

doc/configtables/sector.csv

@@ -23,6 +23,15 @@ district_heating,--,,`prepare_sector_network.py <https://github.com/PyPSA/pypsa-
 -- -- heat_exchanger_pinch_point_temperature_difference,K,float,Heat pump pinch point temperature difference in heat exchangers assumed for approximation.
 -- -- isentropic_compressor_efficiency,--,float,Isentropic efficiency of heat pump compressor assumed for approximation. Must be between 0 and 1.
 -- -- heat_loss,--,float,Heat pump heat loss assumed for approximation. Must be between 0 and 1.
+-- heat_utilisation_potentials,--,Dictionary with names of heat sources for which data by Fraunhofer ISI (`Manz et al. 2024 <https://www.sciencedirect.com/science/article/pii/S0960148124001769>) should be used,
+-- -- geothermal,-,Name of the heat source. Must be the same as in ``heat_pump_sources``,
+-- -- -- key,-,string used to complete URL for data download - e.g. `geothermal_65` or `geothermal_85`","i.e file names in `Fordatis <https://fordatis.fraunhofer.de/handle/fordatis/341.3?mode=simple>`,
+-- -- -- constant_temperature_celsius,°C,heat source temperature,
+-- -- -- column_name,-,name of the data column in retrieved GeoDataFrame,
+
+-- -- -- unit,-,unit of heat source potential must be in (K/M/G/T)Wh,
+-- -- -- full_load_hours,h,assumed full-load hours in Manz et al. (used to scale from utilisation to technical potential),
+-- direct_utilisation_heat_sources,--,List of heat sources for direct heat utilisation in district heating. Must be in the keys of `heat_utilisation_potentials` (e.g. ``geothermal``),
 -- heat_pump_sources,--,,
 -- -- urban central,--,List of heat sources for heat pumps in urban central heating,
 -- -- urban decentral,--,List of heat sources for heat pumps in urban decentral heating,

doc/data-retrieval.rst

@@ -9,6 +9,15 @@ Specific retrieval rules
 
 Data in this section is retrieved and extracted in rules specified in ``rules/retrieve.smk``.
 
+
+``data/fraunhofer_heat_source_utilisation_potentials``
+
+- **Source:** Fraunhofer Fordatis
+- **Link:** https://fordatis.fraunhofer.de/handle/fordatis/341.3?mode=simple
+- **License:** `CC BY 4.0 <https://creativecommons.org/licenses/by/4.0/>`__
+- **Description:** Utilisation potentials for different heat sources across Europe, based on Manz et al. 2024<https://doi.org/10.1016/j.renene.2024.120111>.
+
+
 ``data/nuts``
 
 - **Source:** GISCO

doc/release_notes.rst

@@ -11,6 +11,8 @@ Release Notes
 Upcoming Release
 ================
 
+* Feature: Introduce geothermal district heating (direct utilisation and heat pumps). Potentials are based on `Manz et al. 2024: Spatial analysis of renewable and excess heat potentials for climate-neutral district heating in Europe <https://www.sciencedirect.com/science/article/pii/S0960148124001769>`.
+
 * Feature: Allow CHPs to use different fuel sources such as gas, oil, coal, and methanol. Note that the cost assumptions are based on a gas CHP (except for solid biomass-fired CHP).
 
 * Improve `sanitize_carrier`` function by filling in colors of missing carriers with colors mapped after using the function `rename_techs`.

rules/build_sector.smk

@@ -284,6 +284,28 @@ rule build_central_heating_temperature_profiles:
         "../scripts/build_central_heating_temperature_profiles/run.py"
 
 
+rule build_heat_source_potentials:
+    params:
+        heat_utilisation_potentials=config_provider(
+            "sector", "district_heating", "heat_utilisation_potentials"
+        ),
+    input:
+        utilisation_potential="data/heat_source_utilisation_potentials/{heat_source}.gpkg",
+        regions_onshore=resources("regions_onshore_base_s_{clusters}.geojson"),
+    output:
+        resources("heat_source_potential_{heat_source}_base_s_{clusters}.csv"),
+    resources:
+        mem_mb=2000,
+    log:
+        logs("build_heat_source_potentials_{heat_source}_s_{clusters}.log"),
+    benchmark:
+        benchmarks("build_heat_source_potentials/{heat_source}_s_{clusters}")
+    conda:
+        "../envs/environment.yaml"
+    script:
+        "../scripts/build_heat_source_potentials/run.py"
+
+
 rule build_cop_profiles:
     params:
         heat_pump_sink_T_decentral_heating=config_provider(
@@ -296,6 +318,9 @@ rule build_cop_profiles:
             "sector", "district_heating", "heat_pump_cop_approximation"
         ),
         heat_pump_sources=config_provider("sector", "heat_pump_sources"),
+        heat_utilisation_potentials=config_provider(
+            "sector", "district_heating", "heat_utilisation_potentials"
+        ),
         snapshots=config_provider("snapshots"),
     input:
         central_heating_forward_temperature_profiles=resources(
@@ -321,6 +346,39 @@ rule build_cop_profiles:
         "../scripts/build_cop_profiles/run.py"
 
 
+rule build_direct_heat_source_utilisation_profiles:
+    params:
+        direct_utilisation_heat_sources=config_provider(
+            "sector", "district_heating", "direct_utilisation_heat_sources"
+        ),
+        heat_utilisation_potentials=config_provider(
+            "sector", "district_heating", "heat_utilisation_potentials"
+        ),
+        snapshots=config_provider("snapshots"),
+    input:
+        central_heating_forward_temperature_profiles=resources(
+            "central_heating_forward_temperature_profiles_base_s_{clusters}_{planning_horizons}.nc"
+        ),
+    output:
+        direct_heat_source_utilisation_profiles=resources(
+            "direct_heat_source_utilisation_profiles_base_s_{clusters}_{planning_horizons}.nc"
+        ),
+    resources:
+        mem_mb=20000,
+    log:
+        logs(
+            "build_direct_heat_source_utilisation_profiles_s_{clusters}_{planning_horizons}.log"
+        ),
+    benchmark:
+        benchmarks(
+            "build_direct_heat_source_utilisation_profiles/s_{clusters}_{planning_horizons}"
+        )
+    conda:
+        "../envs/environment.yaml"
+    script:
+        "../scripts/build_direct_heat_source_utilisation_profiles.py"
+
+
 def solar_thermal_cutout(wildcards):
     c = config_provider("solar_thermal", "cutout")(wildcards)
     if c == "default":
@@ -1004,6 +1062,20 @@ rule build_egs_potentials:
         "../scripts/build_egs_potentials.py"
 
 
+def input_heat_source_potentials(w):
+
+    return {
+        heat_source_name: resources(
+            "heat_source_potential_" + heat_source_name + "_base_s_{clusters}.csv"
+        )
+        for heat_source_name in config_provider(
+            "sector", "district_heating", "heat_utilisation_potentials"
+        )(w).keys()
+        if heat_source_name
+        in config_provider("sector", "heat_pump_sources", "urban central")(w)
+    }
+
+
 rule prepare_sector_network:
     params:
         time_resolution=config_provider("clustering", "temporal", "resolution_sector"),
@@ -1028,8 +1100,15 @@ rule prepare_sector_network:
         heat_pump_sources=config_provider("sector", "heat_pump_sources"),
         heat_systems=config_provider("sector", "heat_systems"),
         energy_totals_year=config_provider("energy", "energy_totals_year"),
+        heat_utilisation_potentials=config_provider(
+            "sector", "district_heating", "heat_utilisation_potentials"
+        ),
+        direct_utilisation_heat_sources=config_provider(
+            "sector", "district_heating", "direct_utilisation_heat_sources"
+        ),
     input:
         unpack(input_profile_offwind),
+        unpack(input_heat_source_potentials),
         **rules.cluster_gas_network.output,
         **rules.build_gas_input_locations.output,
         snapshot_weightings=resources(
@@ -1119,6 +1198,9 @@ rule prepare_sector_network:
             if config_provider("sector", "enhanced_geothermal", "enable")(w)
             else []
         ),
+        direct_heat_source_utilisation_profiles=resources(
+            "direct_heat_source_utilisation_profiles_base_s_{clusters}_{planning_horizons}.nc"
+        ),
     output:
         RESULTS
         + "prenetworks/base_s_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}.nc",

rules/retrieve.smk

@@ -631,3 +631,21 @@ if config["enable"]["retrieve"] and (
                 "data/osm-raw/{country}/substations_relation.json",
                 country=config_provider("countries"),
             ),
+
+
+if config["enable"]["retrieve"]:
+
+    rule retrieve_heat_source_utilisation_potentials:
+        params:
+            heat_source="{heat_source}",
+            heat_utilisation_potentials=config_provider(
+                "sector", "district_heating", "heat_utilisation_potentials"
+            ),
+        log:
+            "logs/retrieve_heat_source_potentials_{heat_source}.log",
+        resources:
+            mem_mb=500,
+        output:
+            "data/heat_source_utilisation_potentials/{heat_source}.gpkg",
+        script:
+            "../scripts/retrieve_heat_source_utilisation_potentials.py"

scripts/build_cop_profiles/CentralHeatingCopApproximator.py

@@ -144,11 +144,17 @@ def approximate_cop(self) -> Union[xr.DataArray, np.array]:
         """
         Calculate the coefficient of performance (COP) for the system.
 
+        Notes:
+        ------
+        Returns 0 where the source inlet temperature is greater than the sink outlet temperature.
+
         Returns:
         --------
             Union[xr.DataArray, np.array]: The calculated COP values.
         """
-        return (
+        return xr.where(
+            self.t_source_in_kelvin > self.t_sink_out_kelvin,
+            0,
             self.ideal_lorenz_cop
             * (
                 (
@@ -170,7 +176,7 @@ def approximate_cop(self) -> Union[xr.DataArray, np.array]:
             * (1 - self.ratio_evaporation_compression_work)
             + 1
             - self.isentropic_efficiency_compressor_kelvin
-            - self.heat_loss
+            - self.heat_loss,
         )
 
     @property

scripts/build_cop_profiles/run.py

@@ -4,7 +4,7 @@
 # SPDX-License-Identifier: MIT
 """
 Approximate heat pump coefficient-of-performance (COP) profiles for different
-heat sources and systems.
+heat sources and systems. Returns zero where source temperature higher than sink temperature.
 
 For central heating, this is based on Jensen et al. (2018) (c.f. `CentralHeatingCopApproximator <CentralHeatingCopApproximator.py>`_) and for decentral heating, the approximation is based on Staffell et al. (2012) (c.f. `DecentralHeatingCopApproximator <DecentralHeatingCopApproximator.py>`_).
 
@@ -27,11 +27,8 @@
                 urban central:
                 urban decentral:
                 rural:
-    snapshots:
-
 Inputs
 ------
-- `resources/<run_name>/regions_onshore.geojson`: Onshore regions
 - `resources/<run_name>/temp_soil_total`: Ground temperature
 - `resources/<run_name>/temp_air_total`: Air temperature
 
@@ -94,10 +91,6 @@ def get_cop(
         ).approximate_cop()
 
 
-def get_country_from_node_name(node_name: str) -> str:
-    return node_name[:2]
-
-
 if __name__ == "__main__":
     if "snakemake" not in globals():
         from _helpers import mock_snakemake
@@ -109,9 +102,6 @@ def get_country_from_node_name(node_name: str) -> str:
 
     set_scenario_config(snakemake)
 
-    # map forward and return temperatures specified on country-level to onshore regions
-    regions_onshore = gpd.read_file(snakemake.input.regions_onshore)["name"]
-    snapshots = pd.date_range(freq="h", **snakemake.params.snapshots)
     central_heating_forward_temperature: xr.DataArray = xr.open_dataarray(
         snakemake.input.central_heating_forward_temperature_profiles
     )
@@ -123,9 +113,23 @@ def get_country_from_node_name(node_name: str) -> str:
     for heat_system_type, heat_sources in snakemake.params.heat_pump_sources.items():
         cop_this_system_type = []
         for heat_source in heat_sources:
-            source_inlet_temperature_celsius = xr.open_dataarray(
-                snakemake.input[f"temp_{heat_source.replace('ground', 'soil')}_total"]
-            )
+            if heat_source in ["ground", "air"]:
+                source_inlet_temperature_celsius = xr.open_dataarray(
+                    snakemake.input[
+                        f"temp_{heat_source.replace('ground', 'soil')}_total"
+                    ]
+                )
+            elif heat_source in snakemake.params.heat_utilisation_potentials.keys():
+                source_inlet_temperature_celsius = (
+                    snakemake.params.heat_utilisation_potentials[heat_source][
+                        "constant_temperature_celsius"
+                    ]
+                )
+            else:
+                raise ValueError(
+                    f"Unknown heat source {heat_source}. Must be one of [ground, air] or {snakemake.params.heat_sources.keys()}."
+                )
+
             cop_da = get_cop(
                 heat_system_type=heat_system_type,
                 heat_source=heat_source,