Merge pull request #1175 from PyPSA/update-district-heating-cops

Approximate district heating COPs via Jensen et al. 2018

config/config.default.yaml

@@ -410,6 +410,22 @@ sector:
       2045: 0.8
       2050: 1.0
     district_heating_loss: 0.15
+    forward_temperature: 90 #C
+    return_temperature: 50 #C
+    heat_source_cooling: 6 #K
+    heat_pump_cop_approximation:
+      refrigerant: ammonia
+      heat_exchanger_pinch_point_temperature_difference: 5 #K
+      isentropic_compressor_efficiency: 0.8
+      heat_loss: 0.0
+  heat_pump_sources:
+    urban central:
+    - air
+    urban decentral:
+    - air
+    rural:
+    - air
+    - ground
   cluster_heat_buses: true
   heat_demand_cutout: default
   bev_dsm_restriction_value: 0.75
@@ -492,7 +508,7 @@ sector:
   aviation_demand_factor: 1.
   HVC_demand_factor: 1.
   time_dep_hp_cop: true
-  heat_pump_sink_T: 55.
+  heat_pump_sink_T_individual_heating: 55.
   reduce_space_heat_exogenously: true
   reduce_space_heat_exogenously_factor:
     2020: 0.10  # this results in a space heat demand reduction of 10%

doc/conf.py

@@ -342,5 +342,5 @@ texinfo_documents = [
 
 # Example configuration for intersphinx: refer to the Python standard library.
 intersphinx_mapping = {
-    'https://docs.python.org/': ('https://docs.python.org/3', None),
+    "https://docs.python.org/": ("https://docs.python.org/3", None),
 }

doc/configtables/sector.csv

@@ -9,6 +9,18 @@ district_heating,--,,`prepare_sector_network.py <https://github.com/PyPSA/pypsa-
 -- potential,--,float,maximum fraction of urban demand which can be supplied by district heating
 -- progress,--,Dictionary with planning horizons as keys., Increase of today's district heating demand to potential maximum district heating share. Progress = 0 means today's district heating share. Progress = 1 means maximum fraction of urban demand is supplied by district heating
 -- district_heating_loss,--,float,Share increase in district heat demand in urban central due to heat losses
+-- forward_temperature,°C,float,Forward temperature in district heating
+-- return_temperature,°C,float,Return temperature in district heating. Must be lower than forward temperature
+-- heat_source_cooling,K,float,Cooling of heat source for heat pumps
+-- heat_pump_cop_approximation,,,
+-- -- refrigerant,--,"{ammonia, isobutane}",Heat pump refrigerant assumed for COP approximation
+-- -- 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_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,
+-- -- rural,--,List of heat sources for heat pumps in rural heating,
 cluster_heat_buses,--,"{true, false}",Cluster residential and service heat buses in `prepare_sector_network.py <https://github.com/PyPSA/pypsa-eur-sec/blob/master/scripts/prepare_sector_network.py>`_  to one to save memory.
 ,,,
 bev_dsm_restriction _value,--,float,Adds a lower state of charge (SOC) limit for battery electric vehicles (BEV) to manage its own energy demand (DSM). Located in `build_transport_demand.py <https://github.com/PyPSA/pypsa-eur-sec/blob/master/scripts/build_transport_demand.py>`_. Set to 0 for no restriction on BEV DSM

doc/release_notes.rst

@@ -10,6 +10,8 @@ Release Notes
 Upcoming Release
 ================
 
+* Changed heat pump COP approximation for central heating to be based on `Jensen et al. (2018) <https://backend.orbit.dtu.dk/ws/portalfiles/portal/151965635/MAIN_Final.pdf>`__ and a default forward temperature of 90C. This is more realistic for district heating than the previously used approximation method.
+
 * split solid biomass potentials into solid biomass and municipal solid waste. Add option to use municipal solid waste. This option is only activated in combination with the flag ``waste_to_energy``
 
 * Add option to import solid biomass

rules/build_sector.smk

@@ -217,13 +217,27 @@ rule build_temperature_profiles:
 
 rule build_cop_profiles:
     params:
-        heat_pump_sink_T=config_provider("sector", "heat_pump_sink_T"),
+        heat_pump_sink_T_decentral_heating=config_provider(
+            "sector", "heat_pump_sink_T_individual_heating"
+        ),
+        forward_temperature_central_heating=config_provider(
+            "sector", "district_heating", "forward_temperature"
+        ),
+        return_temperature_central_heating=config_provider(
+            "sector", "district_heating", "return_temperature"
+        ),
+        heat_source_cooling_central_heating=config_provider(
+            "sector", "district_heating", "heat_source_cooling"
+        ),
+        heat_pump_cop_approximation_central_heating=config_provider(
+            "sector", "district_heating", "heat_pump_cop_approximation"
+        ),
+        heat_pump_sources=config_provider("sector", "heat_pump_sources"),
     input:
         temp_soil_total=resources("temp_soil_total_elec_s{simpl}_{clusters}.nc"),
         temp_air_total=resources("temp_air_total_elec_s{simpl}_{clusters}.nc"),
     output:
-        cop_soil_total=resources("cop_soil_total_elec_s{simpl}_{clusters}.nc"),
+        cop_profiles=resources("cop_profiles_elec_s{simpl}_{clusters}.nc"),
-        cop_air_total=resources("cop_air_total_elec_s{simpl}_{clusters}.nc"),
     resources:
         mem_mb=20000,
     log:
@@ -233,7 +247,7 @@ rule build_cop_profiles:
     conda:
         "../envs/environment.yaml"
     script:
-        "../scripts/build_cop_profiles.py"
+        "../scripts/build_cop_profiles/run.py"
 
 
 def solar_thermal_cutout(wildcards):
@@ -941,6 +955,8 @@ rule prepare_sector_network:
         adjustments=config_provider("adjustments", "sector"),
         emissions_scope=config_provider("energy", "emissions"),
         RDIR=RDIR,
+        heat_pump_sources=config_provider("sector", "heat_pump_sources"),
+        heat_systems=config_provider("sector", "heat_systems"),
     input:
         unpack(input_profile_offwind),
         **rules.cluster_gas_network.output,
@@ -1017,8 +1033,7 @@ rule prepare_sector_network:
         ),
         temp_soil_total=resources("temp_soil_total_elec_s{simpl}_{clusters}.nc"),
         temp_air_total=resources("temp_air_total_elec_s{simpl}_{clusters}.nc"),
-        cop_soil_total=resources("cop_soil_total_elec_s{simpl}_{clusters}.nc"),
+        cop_profiles=resources("cop_profiles_elec_s{simpl}_{clusters}.nc"),
-        cop_air_total=resources("cop_air_total_elec_s{simpl}_{clusters}.nc"),
         solar_thermal_total=lambda w: (
             resources("solar_thermal_total_elec_s{simpl}_{clusters}.nc")
             if config_provider("sector", "solar_thermal")(w)

rules/solve_myopic.smk

@@ -9,6 +9,7 @@ rule add_existing_baseyear:
         sector=config_provider("sector"),
         existing_capacities=config_provider("existing_capacities"),
         costs=config_provider("costs"),
+        heat_pump_sources=config_provider("sector", "heat_pump_sources"),
     input:
         network=RESULTS
         + "prenetworks/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}.nc",
@@ -21,8 +22,7 @@ rule add_existing_baseyear:
                 config_provider("scenario", "planning_horizons", 0)(w)
             )
         ),
-        cop_soil_total=resources("cop_soil_total_elec_s{simpl}_{clusters}.nc"),
+        cop_profiles=resources("cop_profiles_elec_s{simpl}_{clusters}.nc"),
-        cop_air_total=resources("cop_air_total_elec_s{simpl}_{clusters}.nc"),
         existing_heating_distribution=resources(
             "existing_heating_distribution_elec_s{simpl}_{clusters}_{planning_horizons}.csv"
         ),
@@ -69,6 +69,7 @@ rule add_brownfield:
         snapshots=config_provider("snapshots"),
         drop_leap_day=config_provider("enable", "drop_leap_day"),
         carriers=config_provider("electricity", "renewable_carriers"),
+        heat_pump_sources=config_provider("sector", "heat_pump_sources"),
     input:
         unpack(input_profile_tech_brownfield),
         simplify_busmap=resources("busmap_elec_s{simpl}.csv"),
@@ -77,8 +78,7 @@ rule add_brownfield:
         + "prenetworks/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}.nc",
         network_p=solved_previous_horizon,  #solved network at previous time step
         costs=resources("costs_{planning_horizons}.csv"),
-        cop_soil_total=resources("cop_soil_total_elec_s{simpl}_{clusters}.nc"),
+        cop_profiles=resources("cop_profiles_elec_s{simpl}_{clusters}.nc"),
-        cop_air_total=resources("cop_air_total_elec_s{simpl}_{clusters}.nc"),
     output:
         RESULTS
         + "prenetworks-brownfield/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}.nc",

rules/solve_perfect.smk

@@ -7,6 +7,7 @@ rule add_existing_baseyear:
         sector=config_provider("sector"),
         existing_capacities=config_provider("existing_capacities"),
         costs=config_provider("costs"),
+        heat_pump_sources=config_provider("sector", "heat_pump_sources"),
     input:
         network=RESULTS
         + "prenetworks/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}.nc",
@@ -19,8 +20,7 @@ rule add_existing_baseyear:
                 config_provider("scenario", "planning_horizons", 0)(w)
             )
         ),
-        cop_soil_total=resources("cop_soil_total_elec_s{simpl}_{clusters}.nc"),
+        cop_profiles=resources("cop_profiles_elec_s{simpl}_{clusters}.nc"),
-        cop_air_total=resources("cop_air_total_elec_s{simpl}_{clusters}.nc"),
         existing_heating_distribution=resources(
             "existing_heating_distribution_elec_s{simpl}_{clusters}_{planning_horizons}.csv"
         ),

scripts/add_existing_baseyear.py

@@ -24,6 +24,10 @@ from _helpers import (
 from add_electricity import sanitize_carriers
 from prepare_sector_network import cluster_heat_buses, define_spatial, prepare_costs
 
+from scripts.definitions.heat_sector import HeatSector
+from scripts.definitions.heat_system import HeatSystem
+from scripts.definitions.heat_system_type import HeatSystemType
+
 logger = logging.getLogger(__name__)
 cc = coco.CountryConverter()
 idx = pd.IndexSlice
@@ -416,14 +420,14 @@ def add_power_capacities_installed_before_baseyear(n, grouping_years, costs, bas
 
 
 def add_heating_capacities_installed_before_baseyear(
-    n,
+    n: pypsa.Network,
-    baseyear,
+    baseyear: int,
-    grouping_years,
+    grouping_years: list,
-    ashp_cop,
+    cop: dict,
-    gshp_cop,
+    time_dep_hp_cop: bool,
-    time_dep_hp_cop,
+    costs: pd.DataFrame,
-    costs,
+    default_lifetime: int,
-    default_lifetime,
+    existing_heating: pd.DataFrame,
 ):
     """
     Parameters
@@ -435,141 +439,158 @@ def add_heating_capacities_installed_before_baseyear(
         currently assumed heating capacities split between residential and
         services proportional to heating load in both 50% capacities
         in rural buses 50% in urban buses
+    cop: xr.DataArray
+        DataArray with time-dependent coefficients of performance (COPs) heat pumps. Coordinates are heat sources (see config), heat system types (see :file:`scripts/enums/HeatSystemType.py`), nodes and snapshots.
+    time_dep_hp_cop: bool
+        If True, time-dependent (dynamic) COPs are used for heat pumps
     """
     logger.debug(f"Adding heating capacities installed before {baseyear}")
 
-    existing_heating = pd.read_csv(
+    for heat_system in existing_heating.columns.get_level_values(0).unique():
-        snakemake.input.existing_heating_distribution, header=[0, 1], index_col=0
+        heat_system = HeatSystem(heat_system)
-    )
 
-    for name in existing_heating.columns.get_level_values(0).unique():
+        nodes = pd.Index(
-        name_type = "central" if name == "urban central" else "decentral"
+            n.buses.location[n.buses.index.str.contains(f"{heat_system} heat")]
+        )
 
-        nodes = pd.Index(n.buses.location[n.buses.index.str.contains(f"{name} heat")])
+        if (not heat_system == HeatSystem.URBAN_CENTRAL) and options[
-
+            "electricity_distribution_grid"
-        if (name_type != "central") and options["electricity_distribution_grid"]:
+        ]:
             nodes_elec = nodes + " low voltage"
         else:
             nodes_elec = nodes
 
-        heat_pump_type = "air" if "urban" in name else "ground"
+            too_large_grouping_years = [
-
+                gy for gy in grouping_years if gy >= int(baseyear)
-        # Add heat pumps
-        costs_name = f"decentral {heat_pump_type}-sourced heat pump"
-
-        cop = {"air": ashp_cop, "ground": gshp_cop}
-
-        if time_dep_hp_cop:
-            efficiency = cop[heat_pump_type][nodes]
-        else:
-            efficiency = costs.at[costs_name, "efficiency"]
-
-        too_large_grouping_years = [gy for gy in grouping_years if gy >= int(baseyear)]
-        if too_large_grouping_years:
-            logger.warning(
-                f"Grouping years >= baseyear are ignored. Dropping {too_large_grouping_years}."
-            )
-        valid_grouping_years = pd.Series(
-            [
-                int(grouping_year)
-                for grouping_year in grouping_years
-                if int(grouping_year) + default_lifetime > int(baseyear)
-                and int(grouping_year) < int(baseyear)
             ]
-        )
+            if too_large_grouping_years:
+                logger.warning(
+                    f"Grouping years >= baseyear are ignored. Dropping {too_large_grouping_years}."
+                )
+            valid_grouping_years = pd.Series(
+                [
+                    int(grouping_year)
+                    for grouping_year in grouping_years
+                    if int(grouping_year) + default_lifetime > int(baseyear)
+                    and int(grouping_year) < int(baseyear)
+                ]
+            )
 
-        assert valid_grouping_years.is_monotonic_increasing
+            assert valid_grouping_years.is_monotonic_increasing
 
-        # get number of years of each interval
+            # get number of years of each interval
-        _years = valid_grouping_years.diff()
+            _years = valid_grouping_years.diff()
-        # Fill NA from .diff() with value for the first interval
+            # Fill NA from .diff() with value for the first interval
-        _years[0] = valid_grouping_years[0] - baseyear + default_lifetime
+            _years[0] = valid_grouping_years[0] - baseyear + default_lifetime
-        # Installation is assumed to be linear for the past
+            # Installation is assumed to be linear for the past
-        ratios = _years / _years.sum()
+            ratios = _years / _years.sum()
 
         for ratio, grouping_year in zip(ratios, valid_grouping_years):
+            # Add heat pumps
+            for heat_source in snakemake.params.heat_pump_sources[
+                heat_system.system_type.value
+            ]:
+                costs_name = heat_system.heat_pump_costs_name(heat_source)
 
-            n.madd(
+                efficiency = (
-                "Link",
+                    cop.sel(
-                nodes,
+                        heat_system=heat_system.system_type.value,
-                suffix=f" {name} {heat_pump_type} heat pump-{grouping_year}",
+                        heat_source=heat_source,
-                bus0=nodes_elec,
+                        name=nodes,
-                bus1=nodes + " " + name + " heat",
+                    )
-                carrier=f"{name} {heat_pump_type} heat pump",
+                    .to_pandas()
-                efficiency=efficiency,
+                    .reindex(index=n.snapshots)
-                capital_cost=costs.at[costs_name, "efficiency"]
+                    if time_dep_hp_cop
-                * costs.at[costs_name, "fixed"],
+                    else costs.at[costs_name, "efficiency"]
-                p_nom=existing_heating.loc[nodes, (name, f"{heat_pump_type} heat pump")]
+                )
-                * ratio
+
-                / costs.at[costs_name, "efficiency"],
+                n.madd(
-                build_year=int(grouping_year),
+                    "Link",
-                lifetime=costs.at[costs_name, "lifetime"],
+                    nodes,
-            )
+                    suffix=f" {heat_system} {heat_source} heat pump-{grouping_year}",
+                    bus0=nodes_elec,
+                    bus1=nodes + " " + heat_system.value + " heat",
+                    carrier=f"{heat_system} {heat_source} heat pump",
+                    efficiency=efficiency,
+                    capital_cost=costs.at[costs_name, "efficiency"]
+                    * costs.at[costs_name, "fixed"],
+                    p_nom=existing_heating.loc[
+                        nodes, (heat_system.value, f"{heat_source} heat pump")
+                    ]
+                    * ratio
+                    / costs.at[costs_name, "efficiency"],
+                    build_year=int(grouping_year),
+                    lifetime=costs.at[costs_name, "lifetime"],
+                )
 
             # add resistive heater, gas boilers and oil boilers
             n.madd(
                 "Link",
                 nodes,
-                suffix=f" {name} resistive heater-{grouping_year}",
+                suffix=f" {heat_system} resistive heater-{grouping_year}",
                 bus0=nodes_elec,
-                bus1=nodes + " " + name + " heat",
+                bus1=nodes + " " + heat_system.value + " heat",
-                carrier=name + " resistive heater",
+                carrier=heat_system.value + " resistive heater",
-                efficiency=costs.at[f"{name_type} resistive heater", "efficiency"],
+                efficiency=costs.at[
+                    heat_system.resistive_heater_costs_name, "efficiency"
+                ],
                 capital_cost=(
-                    costs.at[f"{name_type} resistive heater", "efficiency"]
+                    costs.at[heat_system.resistive_heater_costs_name, "efficiency"]
-                    * costs.at[f"{name_type} resistive heater", "fixed"]
+                    * costs.at[heat_system.resistive_heater_costs_name, "fixed"]
                 ),
                 p_nom=(
-                    existing_heating.loc[nodes, (name, "resistive heater")]
+                    existing_heating.loc[nodes, (heat_system.value, "resistive heater")]
                     * ratio
-                    / costs.at[f"{name_type} resistive heater", "efficiency"]
+                    / costs.at[heat_system.resistive_heater_costs_name, "efficiency"]
                 ),
                 build_year=int(grouping_year),
-                lifetime=costs.at[f"{name_type} resistive heater", "lifetime"],
+                lifetime=costs.at[heat_system.resistive_heater_costs_name, "lifetime"],
             )
 
             n.madd(
                 "Link",
                 nodes,
-                suffix=f" {name} gas boiler-{grouping_year}",
+                suffix=f"{heat_system} gas boiler-{grouping_year}",
                 bus0="EU gas" if "EU gas" in spatial.gas.nodes else nodes + " gas",
-                bus1=nodes + " " + name + " heat",
+                bus1=nodes + " " + heat_system.value + " heat",
                 bus2="co2 atmosphere",
-                carrier=name + " gas boiler",
+                carrier=heat_system.value + " gas boiler",
-                efficiency=costs.at[f"{name_type} gas boiler", "efficiency"],
+                efficiency=costs.at[heat_system.gas_boiler_costs_name, "efficiency"],
                 efficiency2=costs.at["gas", "CO2 intensity"],
                 capital_cost=(
-                    costs.at[f"{name_type} gas boiler", "efficiency"]
+                    costs.at[heat_system.gas_boiler_costs_name, "efficiency"]
-                    * costs.at[f"{name_type} gas boiler", "fixed"]
+                    * costs.at[heat_system.gas_boiler_costs_name, "fixed"]
                 ),
                 p_nom=(
-                    existing_heating.loc[nodes, (name, "gas boiler")]
+                    existing_heating.loc[nodes, (heat_system.value, "gas boiler")]
                     * ratio
-                    / costs.at[f"{name_type} gas boiler", "efficiency"]
+                    / costs.at[heat_system.gas_boiler_costs_name, "efficiency"]
                 ),
                 build_year=int(grouping_year),
-                lifetime=costs.at[f"{name_type} gas boiler", "lifetime"],
+                lifetime=costs.at[heat_system.gas_boiler_costs_name, "lifetime"],
             )
 
             n.madd(
                 "Link",
                 nodes,
-                suffix=f" {name} oil boiler-{grouping_year}",
+                suffix=f" {heat_system} oil boiler-{grouping_year}",
                 bus0=spatial.oil.nodes,
-                bus1=nodes + " " + name + " heat",
+                bus1=nodes + " " + heat_system.value + " heat",
                 bus2="co2 atmosphere",
-                carrier=name + " oil boiler",
+                carrier=heat_system.value + " oil boiler",
-                efficiency=costs.at["decentral oil boiler", "efficiency"],
+                efficiency=costs.at[heat_system.oil_boiler_costs_name, "efficiency"],
                 efficiency2=costs.at["oil", "CO2 intensity"],
-                capital_cost=costs.at["decentral oil boiler", "efficiency"]
+                capital_cost=costs.at[heat_system.oil_boiler_costs_name, "efficiency"]
-                * costs.at["decentral oil boiler", "fixed"],
+                * costs.at[heat_system.oil_boiler_costs_name, "fixed"],
                 p_nom=(
-                    existing_heating.loc[nodes, (name, "oil boiler")]
+                    existing_heating.loc[nodes, (heat_system.value, "oil boiler")]
                     * ratio
-                    / costs.at["decentral oil boiler", "efficiency"]
+                    / costs.at[heat_system.oil_boiler_costs_name, "efficiency"]
                 ),
                 build_year=int(grouping_year),
-                lifetime=costs.at[f"{name_type} gas boiler", "lifetime"],
+                lifetime=costs.at[
+                    f"{heat_system.central_or_decentral} gas boiler", "lifetime"
+                ],
             )
 
             # delete links with p_nom=nan corresponding to extra nodes in country
@@ -639,29 +660,22 @@ if __name__ == "__main__":
     )
 
     if options["heating"]:
-        time_dep_hp_cop = options["time_dep_hp_cop"]
+
-        ashp_cop = (
-            xr.open_dataarray(snakemake.input.cop_air_total)
-            .to_pandas()
-            .reindex(index=n.snapshots)
-        )
-        gshp_cop = (
-            xr.open_dataarray(snakemake.input.cop_soil_total)
-            .to_pandas()
-            .reindex(index=n.snapshots)
-        )
-        default_lifetime = snakemake.params.existing_capacities[
-            "default_heating_lifetime"
-        ]
         add_heating_capacities_installed_before_baseyear(
-            n,
+            n=n,
-            baseyear,
+            baseyear=baseyear,
-            grouping_years_heat,
+            grouping_years=grouping_years_heat,
-            ashp_cop,
+            cop=xr.open_dataarray(snakemake.input.cop_profiles),
-            gshp_cop,
+            time_dep_hp_cop=options["time_dep_hp_cop"],
-            time_dep_hp_cop,
+            costs=costs,
-            costs,
+            default_lifetime=snakemake.params.existing_capacities[
-            default_lifetime,
+                "default_heating_lifetime"
+            ],
+            existing_heating=pd.read_csv(
+                snakemake.input.existing_heating_distribution,
+                header=[0, 1],
+                index_col=0,
+            ),
         )
 
     if options.get("cluster_heat_buses", False):

scripts/build_cop_profiles.py

@@ -1,69 +0,0 @@
-# -*- coding: utf-8 -*-
-# SPDX-FileCopyrightText: : 2020-2024 The PyPSA-Eur Authors
-#
-# SPDX-License-Identifier: MIT
-"""
-Build coefficient of performance (COP) time series for air- or ground-sourced
-heat pumps.
-
-The COP is approximated as a quatratic function of the temperature difference between source and
-sink, based on Staffell et al. 2012.
-
-This rule is executed in ``build_sector.smk``.
-
-Relevant Settings
------------------
-
-.. code:: yaml
-    heat_pump_sink_T:
-
-
-Inputs:
--------
-- ``resources/<run_name>/temp_soil_total_elec_s<simpl>_<clusters>.nc``: Soil temperature (total) time series.
-- ``resources/<run_name>/temp_air_total_elec_s<simpl>_<clusters>.nc``: Ambient air temperature (total) time series.
-
-Outputs:
---------
-- ``resources/cop_soil_total_elec_s<simpl>_<clusters>.nc``: COP (ground-sourced) time series (total).
-- ``resources/cop_air_total_elec_s<simpl>_<clusters>.nc``: COP (air-sourced) time series (total).
-
-
-References
-----------
-[1] Staffell et al., Energy & Environmental Science 11 (2012): A review of domestic heat pumps, https://doi.org/10.1039/C2EE22653G.
-"""
-
-import xarray as xr
-from _helpers import set_scenario_config
-
-
-def coefficient_of_performance(delta_T, source="air"):
-    if source == "air":
-        return 6.81 - 0.121 * delta_T + 0.000630 * delta_T**2
-    elif source == "soil":
-        return 8.77 - 0.150 * delta_T + 0.000734 * delta_T**2
-    else:
-        raise NotImplementedError("'source' must be one of  ['air', 'soil']")
-
-
-if __name__ == "__main__":
-    if "snakemake" not in globals():
-        from _helpers import mock_snakemake
-
-        snakemake = mock_snakemake(
-            "build_cop_profiles",
-            simpl="",
-            clusters=48,
-        )
-
-    set_scenario_config(snakemake)
-
-    for source in ["air", "soil"]:
-        source_T = xr.open_dataarray(snakemake.input[f"temp_{source}_total"])
-
-        delta_T = snakemake.params.heat_pump_sink_T - source_T
-
-        cop = coefficient_of_performance(delta_T, source)
-
-        cop.to_netcdf(snakemake.output[f"cop_{source}_total"])

scripts/build_cop_profiles/BaseCopApproximator.py

@@ -0,0 +1,111 @@
+# -*- coding: utf-8 -*-
+# SPDX-FileCopyrightText: : 2020-2024 The PyPSA-Eur Authors
+#
+# SPDX-License-Identifier: MIT
+
+from abc import ABC, abstractmethod
+from typing import Union
+
+import numpy as np
+import xarray as xr
+
+
+class BaseCopApproximator(ABC):
+    """
+    Abstract class for approximating the coefficient of performance (COP) of a
+    heat pump.
+
+    Attributes:
+    ----------
+    forward_temperature_celsius : Union[xr.DataArray, np.array]
+        The forward temperature in Celsius.
+    source_inlet_temperature_celsius : Union[xr.DataArray, np.array]
+        The source inlet temperature in Celsius.
+
+    Methods:
+    -------
+    __init__(self, forward_temperature_celsius, source_inlet_temperature_celsius)
+        Initialize CopApproximator.
+    approximate_cop(self)
+        Approximate heat pump coefficient of performance (COP).
+    celsius_to_kelvin(t_celsius)
+        Convert temperature from Celsius to Kelvin.
+    logarithmic_mean(t_hot, t_cold)
+        Calculate the logarithmic mean temperature difference.
+    """
+
+    def __init__(
+        self,
+        forward_temperature_celsius: Union[xr.DataArray, np.array],
+        source_inlet_temperature_celsius: Union[xr.DataArray, np.array],
+    ):
+        """
+        Initialize CopApproximator.
+
+        Parameters:
+        ----------
+        forward_temperature_celsius : Union[xr.DataArray, np.array]
+            The forward temperature in Celsius.
+        source_inlet_temperature_celsius : Union[xr.DataArray, np.array]
+            The source inlet temperature in Celsius.
+        """
+        pass
+
+    @abstractmethod
+    def approximate_cop(self) -> Union[xr.DataArray, np.array]:
+        """
+        Approximate heat pump coefficient of performance (COP).
+
+        Returns:
+        -------
+        Union[xr.DataArray, np.array]
+            The calculated COP values.
+        """
+        pass
+
+    @staticmethod
+    def celsius_to_kelvin(
+        t_celsius: Union[float, xr.DataArray, np.array]
+    ) -> Union[float, xr.DataArray, np.array]:
+        """
+        Convert temperature from Celsius to Kelvin.
+
+        Parameters:
+        ----------
+        t_celsius : Union[float, xr.DataArray, np.array]
+            Temperature in Celsius.
+
+        Returns:
+        -------
+        Union[float, xr.DataArray, np.array]
+            Temperature in Kelvin.
+        """
+        if (np.asarray(t_celsius) > 200).any():
+            raise ValueError(
+                "t_celsius > 200. Are you sure you are using the right units?"
+            )
+        return t_celsius + 273.15
+
+    @staticmethod
+    def logarithmic_mean(
+        t_hot: Union[float, xr.DataArray, np.ndarray],
+        t_cold: Union[float, xr.DataArray, np.ndarray],
+    ) -> Union[float, xr.DataArray, np.ndarray]:
+        """
+        Calculate the logarithmic mean temperature difference.
+
+        Parameters:
+        ----------
+        t_hot : Union[float, xr.DataArray, np.ndarray]
+            Hot temperature.
+        t_cold : Union[float, xr.DataArray, np.ndarray]
+            Cold temperature.
+
+        Returns:
+        -------
+        Union[float, xr.DataArray, np.ndarray]
+            Logarithmic mean temperature difference.
+        """
+        if (np.asarray(t_hot <= t_cold)).any():
+            raise ValueError("t_hot must be greater than t_cold")
+        return (t_hot - t_cold) / np.log(t_hot / t_cold)

scripts/build_cop_profiles/CentralHeatingCopApproximator.py

@@ -0,0 +1,392 @@
+# -*- coding: utf-8 -*-
+# SPDX-FileCopyrightText: : 2020-2024 The PyPSA-Eur Authors
+#
+# SPDX-License-Identifier: MIT
+
+
+from typing import Union
+
+import numpy as np
+import xarray as xr
+from BaseCopApproximator import BaseCopApproximator
+
+
+class CentralHeatingCopApproximator(BaseCopApproximator):
+    """
+    Approximate the coefficient of performance (COP) for a heat pump in a
+    central heating system (district heating).
+
+    Uses an approximation method proposed by Jensen et al. (2018) and
+    default parameters from Pieper et al. (2020). The method is based on
+    a thermodynamic heat pump model with some hard-to-know parameters
+    being approximated.
+
+    Attributes:
+    ----------
+    forward_temperature_celsius : Union[xr.DataArray, np.array]
+        The forward temperature in Celsius.
+    return_temperature_celsius : Union[xr.DataArray, np.array]
+        The return temperature in Celsius.
+    source_inlet_temperature_celsius : Union[xr.DataArray, np.array]
+        The source inlet temperature in Celsius.
+    source_outlet_temperature_celsius : Union[xr.DataArray, np.array]
+        The source outlet temperature in Celsius.
+    delta_t_pinch_point : float, optional
+        The pinch point temperature difference, by default 5.
+    isentropic_compressor_efficiency : float, optional
+        The isentropic compressor efficiency, by default 0.8.
+    heat_loss : float, optional
+        The heat loss, by default 0.0.
+
+    Methods:
+    -------
+    __init__(
+        forward_temperature_celsius: Union[xr.DataArray, np.array],
+        source_inlet_temperature_celsius: Union[xr.DataArray, np.array],
+        return_temperature_celsius: Union[xr.DataArray, np.array],
+        source_outlet_temperature_celsius: Union[xr.DataArray, np.array],
+        delta_t_pinch_point: float = 5,
+        isentropic_compressor_efficiency: float = 0.8,
+        heat_loss: float = 0.0,
+    ) -> None:
+        Initializes the CentralHeatingCopApproximator object.
+
+    approximate_cop(self) -> Union[xr.DataArray, np.array]:
+        Calculate the coefficient of performance (COP) for the system.
+
+    _approximate_delta_t_refrigerant_source(
+        self, delta_t_source: Union[xr.DataArray, np.array]
+    ) -> Union[xr.DataArray, np.array]:
+        Approximates the temperature difference between the refrigerant and the source.
+
+    _approximate_delta_t_refrigerant_sink(
+        self,
+        refrigerant: str = "ammonia",
+        a: float = {"ammonia": 0.2, "isobutane": -0.0011},
+        b: float = {"ammonia": 0.2, "isobutane": 0.3},
+        c: float = {"ammonia": 0.016, "isobutane": 2.4},
+    ) -> Union[xr.DataArray, np.array]:
+        Approximates the temperature difference between the refrigerant and heat sink.
+
+    _ratio_evaporation_compression_work_approximation(
+        self,
+        refrigerant: str = "ammonia",
+        a: float = {"ammonia": 0.0014, "isobutane": 0.0035},
+    ) -> Union[xr.DataArray, np.array]:
+        Calculate the ratio of evaporation to compression work based on approximation.
+
+    _approximate_delta_t_refrigerant_sink(
+        self,
+        refrigerant: str = "ammonia",
+        a: float = {"ammonia": 0.2, "isobutane": -0.0011},
+        b: float = {"ammonia": 0.2, "isobutane": 0.3},
+        c: float = {"ammonia": 0.016, "isobutane": 2.4},
+    ) -> Union[xr.DataArray, np.array]:
+        Approximates the temperature difference between the refrigerant and heat sink.
+
+    _ratio_evaporation_compression_work_approximation(
+        self,
+        refrigerant: str = "ammonia",
+        a: float = {"ammonia": 0.0014, "isobutane": 0.0035},
+    ) -> Union[xr.DataArray, np.array]:
+        Calculate the ratio of evaporation to compression work based on approximation.
+    """
+
+    def __init__(
+        self,
+        forward_temperature_celsius: Union[xr.DataArray, np.array],
+        source_inlet_temperature_celsius: Union[xr.DataArray, np.array],
+        return_temperature_celsius: Union[xr.DataArray, np.array],
+        source_outlet_temperature_celsius: Union[xr.DataArray, np.array],
+        delta_t_pinch_point: float = 5,
+        isentropic_compressor_efficiency: float = 0.8,
+        heat_loss: float = 0.0,
+    ) -> None:
+        """
+        Initializes the CentralHeatingCopApproximator object.
+
+        Parameters:
+        ----------
+        forward_temperature_celsius : Union[xr.DataArray, np.array]
+            The forward temperature in Celsius.
+        return_temperature_celsius : Union[xr.DataArray, np.array]
+            The return temperature in Celsius.
+        source_inlet_temperature_celsius : Union[xr.DataArray, np.array]
+            The source inlet temperature in Celsius.
+        source_outlet_temperature_celsius : Union[xr.DataArray, np.array]
+            The source outlet temperature in Celsius.
+        delta_t_pinch_point : float, optional
+            The pinch point temperature difference, by default 5.
+        isentropic_compressor_efficiency : float, optional
+            The isentropic compressor efficiency, by default 0.8.
+        heat_loss : float, optional
+            The heat loss, by default 0.0.
+        """
+        self.t_source_in_kelvin = BaseCopApproximator.celsius_to_kelvin(
+            source_inlet_temperature_celsius
+        )
+        self.t_sink_out_kelvin = BaseCopApproximator.celsius_to_kelvin(
+            forward_temperature_celsius
+        )
+
+        self.t_sink_in_kelvin = BaseCopApproximator.celsius_to_kelvin(
+            return_temperature_celsius
+        )
+        self.t_source_out = BaseCopApproximator.celsius_to_kelvin(
+            source_outlet_temperature_celsius
+        )
+
+        self.isentropic_efficiency_compressor_kelvin = isentropic_compressor_efficiency
+        self.heat_loss = heat_loss
+        self.delta_t_pinch = delta_t_pinch_point
+
+    def approximate_cop(self) -> Union[xr.DataArray, np.array]:
+        """
+        Calculate the coefficient of performance (COP) for the system.
+
+        Returns:
+        --------
+            Union[xr.DataArray, np.array]: The calculated COP values.
+        """
+        return (
+            self.ideal_lorenz_cop
+            * (
+                (
+                    1
+                    + (self.delta_t_refrigerant_sink + self.delta_t_pinch)
+                    / self.t_sink_mean_kelvin
+                )
+                / (
+                    1
+                    + (
+                        self.delta_t_refrigerant_sink
+                        + self.delta_t_refrigerant_source
+                        + 2 * self.delta_t_pinch
+                    )
+                    / self.delta_t_lift
+                )
+            )
+            * self.isentropic_efficiency_compressor_kelvin
+            * (1 - self.ratio_evaporation_compression_work)
+            + 1
+            - self.isentropic_efficiency_compressor_kelvin
+            - self.heat_loss
+        )
+
+    @property
+    def t_sink_mean_kelvin(self) -> Union[xr.DataArray, np.array]:
+        """
+        Calculate the logarithmic mean temperature difference between the cold
+        and hot sinks.
+
+        Returns
+        -------
+        Union[xr.DataArray, np.array]
+            The mean temperature difference.
+        """
+        return BaseCopApproximator.logarithmic_mean(
+            t_cold=self.t_sink_in_kelvin, t_hot=self.t_sink_out_kelvin
+        )
+
+    @property
+    def t_source_mean_kelvin(self) -> Union[xr.DataArray, np.array]:
+        """
+        Calculate the logarithmic mean temperature of the heat source.
+
+        Returns
+        -------
+        Union[xr.DataArray, np.array]
+            The mean temperature of the heat source.
+        """
+        return BaseCopApproximator.logarithmic_mean(
+            t_hot=self.t_source_in_kelvin, t_cold=self.t_source_out
+        )
+
+    @property
+    def delta_t_lift(self) -> Union[xr.DataArray, np.array]:
+        """
+        Calculate the temperature lift as the difference between the
+        logarithmic sink and source temperatures.
+
+        Returns
+        -------
+        Union[xr.DataArray, np.array]
+            The temperature difference between the sink and source.
+        """
+        return self.t_sink_mean_kelvin - self.t_source_mean_kelvin
+
+    @property
+    def ideal_lorenz_cop(self) -> Union[xr.DataArray, np.array]:
+        """
+        Ideal Lorenz coefficient of performance (COP).
+
+        The ideal Lorenz COP is calculated as the ratio of the mean sink temperature
+        to the lift temperature difference.
+
+        Returns
+        -------
+        np.array
+            The ideal Lorenz COP.
+        """
+        return self.t_sink_mean_kelvin / self.delta_t_lift
+
+    @property
+    def delta_t_refrigerant_source(self) -> Union[xr.DataArray, np.array]:
+        """
+        Calculate the temperature difference between the refrigerant source
+        inlet and outlet.
+
+        Returns
+        -------
+        Union[xr.DataArray, np.array]
+            The temperature difference between the refrigerant source inlet and outlet.
+        """
+        return self._approximate_delta_t_refrigerant_source(
+            delta_t_source=self.t_source_in_kelvin - self.t_source_out
+        )
+
+    @property
+    def delta_t_refrigerant_sink(self) -> Union[xr.DataArray, np.array]:
+        """
+        Temperature difference between the refrigerant and the sink based on
+        approximation.
+
+        Returns
+        -------
+        Union[xr.DataArray, np.array]
+            The temperature difference between the refrigerant and the sink.
+        """
+        return self._approximate_delta_t_refrigerant_sink()
+
+    @property
+    def ratio_evaporation_compression_work(self) -> Union[xr.DataArray, np.array]:
+        """
+        Calculate the ratio of evaporation to compression work based on
+        approximation.
+
+        Returns
+        -------
+        Union[xr.DataArray, np.array]
+            The calculated ratio of evaporation to compression work.
+        """
+        return self._ratio_evaporation_compression_work_approximation()
+
+    @property
+    def delta_t_sink(self) -> Union[xr.DataArray, np.array]:
+        """
+        Calculate the temperature difference at the sink.
+
+        Returns
+        -------
+        Union[xr.DataArray, np.array]
+            The temperature difference at the sink.
+        """
+        return self.t_sink_out_kelvin - self.t_sink_in_kelvin
+
+    def _approximate_delta_t_refrigerant_source(
+        self, delta_t_source: Union[xr.DataArray, np.array]
+    ) -> Union[xr.DataArray, np.array]:
+        """
+        Approximates the temperature difference between the refrigerant and the
+        source.
+
+        Parameters
+        ----------
+        delta_t_source : Union[xr.DataArray, np.array]
+            The temperature difference for the refrigerant source.
+
+        Returns
+        -------
+        Union[xr.DataArray, np.array]
+            The approximate temperature difference between the refrigerant and heat source.
+        """
+        return delta_t_source / 2
+
+    def _approximate_delta_t_refrigerant_sink(
+        self,
+        refrigerant: str = "ammonia",
+        a: float = {"ammonia": 0.2, "isobutane": -0.0011},
+        b: float = {"ammonia": 0.2, "isobutane": 0.3},
+        c: float = {"ammonia": 0.016, "isobutane": 2.4},
+    ) -> Union[xr.DataArray, np.array]:
+        """
+        Approximates the temperature difference between the refrigerant and
+        heat sink.
+
+        Parameters:
+        ----------
+        refrigerant : str, optional
+            The refrigerant used in the system. Either 'isobutane' or 'ammonia. Default is 'ammonia'.
+        a : float, optional
+            Coefficient for the temperature difference between the sink and source, default is 0.2.
+        b : float, optional
+            Coefficient for the temperature difference at the sink, default is 0.2.
+        c : float, optional
+            Constant term, default is 0.016.
+
+        Returns:
+        -------
+        Union[xr.DataArray, np.array]
+            The approximate temperature difference between the refrigerant and heat sink.
+
+        Notes:
+        ------
+        This function assumes ammonia as the refrigerant.
+
+        The approximate temperature difference at the refrigerant sink is calculated using the following formula:
+        a * (t_sink_out - t_source_out + 2 * delta_t_pinch) + b * delta_t_sink + c
+        """
+        if refrigerant not in a.keys():
+            raise ValueError(
+                f"Invalid refrigerant '{refrigerant}'. Must be one of {a.keys()}"
+            )
+        return (
+            a[refrigerant]
+            * (self.t_sink_out_kelvin - self.t_source_out + 2 * self.delta_t_pinch)
+            + b[refrigerant] * self.delta_t_sink
+            + c[refrigerant]
+        )
+
+    def _ratio_evaporation_compression_work_approximation(
+        self,
+        refrigerant: str = "ammonia",
+        a: float = {"ammonia": 0.0014, "isobutane": 0.0035},
+        b: float = {"ammonia": -0.0015, "isobutane": -0.0033},
+        c: float = {"ammonia": 0.039, "isobutane": 0.053},
+    ) -> Union[xr.DataArray, np.array]:
+        """
+        Calculate the ratio of evaporation to compression work approximation.
+
+        Parameters:
+        ----------
+        refrigerant : str, optional
+            The refrigerant used in the system. Either 'isobutane' or 'ammonia. Default is 'ammonia'.
+        a : float, optional
+            Coefficient 'a' in the approximation equation. Default is 0.0014.
+        b : float, optional
+            Coefficient 'b' in the approximation equation. Default is -0.0015.
+        c : float, optional
+            Coefficient 'c' in the approximation equation. Default is 0.039.
+
+        Returns:
+        -------
+        Union[xr.DataArray, np.array]
+            The approximated ratio of evaporation to compression work.
+
+        Notes:
+        ------
+        This function assumes ammonia as the refrigerant.
+
+        The approximation equation used is:
+        ratio = a * (t_sink_out - t_source_out + 2 * delta_t_pinch) + b * delta_t_sink + c
+        """
+        if refrigerant not in a.keys():
+            raise ValueError(
+                f"Invalid refrigerant '{refrigerant}'. Must be one of {a.keys()}"
+            )
+        return (
+            a[refrigerant]
+            * (self.t_sink_out_kelvin - self.t_source_out + 2 * self.delta_t_pinch)
+            + b[refrigerant] * self.delta_t_sink
+            + c[refrigerant]
+        )

scripts/build_cop_profiles/DecentralHeatingCopApproximator.py

@@ -0,0 +1,110 @@
+# -*- coding: utf-8 -*-
+# SPDX-FileCopyrightText: : 2020-2024 The PyPSA-Eur Authors
+#
+# SPDX-License-Identifier: MIT
+
+
+from typing import Union
+
+import numpy as np
+import xarray as xr
+from BaseCopApproximator import BaseCopApproximator
+
+
+class DecentralHeatingCopApproximator(BaseCopApproximator):
+    """
+    Approximate the coefficient of performance (COP) for a heat pump in a
+    decentral heating system (individual/household heating).
+
+    Uses a quadratic regression on the temperature difference between the source and sink based on empirical data proposed by Staffell et al. 2012.
+
+    Attributes
+    ----------
+    forward_temperature_celsius : Union[xr.DataArray, np.array]
+        The forward temperature in Celsius.
+    source_inlet_temperature_celsius : Union[xr.DataArray, np.array]
+        The source inlet temperature in Celsius.
+    source_type : str
+        The source of the heat pump. Must be either 'air' or 'ground'.
+
+    Methods
+    -------
+    __init__(forward_temperature_celsius, source_inlet_temperature_celsius, source_type)
+        Initialize the DecentralHeatingCopApproximator object.
+    approximate_cop()
+        Compute the COP values using quadratic regression for air-/ground-source heat pumps.
+    _approximate_cop_air_source()
+        Evaluate quadratic regression for an air-sourced heat pump.
+    _approximate_cop_ground_source()
+        Evaluate quadratic regression for a ground-sourced heat pump.
+
+    References
+    ----------
+    [1] Staffell et al., Energy & Environmental Science 11 (2012): A review of domestic heat pumps, https://doi.org/10.1039/C2EE22653G.
+    """
+
+    def __init__(
+        self,
+        forward_temperature_celsius: Union[xr.DataArray, np.array],
+        source_inlet_temperature_celsius: Union[xr.DataArray, np.array],
+        source_type: str,
+    ):
+        """
+        Initialize the DecentralHeatingCopApproximator object.
+
+        Parameters
+        ----------
+        forward_temperature_celsius : Union[xr.DataArray, np.array]
+            The forward temperature in Celsius.
+        source_inlet_temperature_celsius : Union[xr.DataArray, np.array]
+            The source inlet temperature in Celsius.
+        source_type : str
+            The source of the heat pump. Must be either 'air' or 'ground'.
+        """
+
+        self.delta_t = forward_temperature_celsius - source_inlet_temperature_celsius
+        if source_type not in ["air", "ground"]:
+            raise ValueError("'source_type' must be one of ['air', 'ground']")
+        else:
+            self.source_type = source_type
+
+    def approximate_cop(self) -> Union[xr.DataArray, np.array]:
+        """
+        Compute the COP values using quadratic regression for air-/ground-
+        source heat pumps.
+
+        Returns
+        -------
+        Union[xr.DataArray, np.array]
+            The calculated COP values.
+        """
+        if self.source_type == "air":
+            return self._approximate_cop_air_source()
+        elif self.source_type == "ground":
+            return self._approximate_cop_ground_source()
+
+    def _approximate_cop_air_source(self) -> Union[xr.DataArray, np.array]:
+        """
+        Evaluate quadratic regression for an air-sourced heat pump.
+
+        COP = 6.81 - 0.121 * delta_T + 0.000630 * delta_T^2
+
+        Returns
+        -------
+        Union[xr.DataArray, np.array]
+            The calculated COP values.
+        """
+        return 6.81 - 0.121 * self.delta_t + 0.000630 * self.delta_t**2
+
+    def _approximate_cop_ground_source(self) -> Union[xr.DataArray, np.array]:
+        """
+        Evaluate quadratic regression for a ground-sourced heat pump.
+
+        COP = 8.77 - 0.150 * delta_T + 0.000734 * delta_T^2
+
+        Returns
+        -------
+        Union[xr.DataArray, np.array]
+            The calculated COP values.
+        """
+        return 8.77 - 0.150 * self.delta_t + 0.000734 * self.delta_t**2

scripts/build_cop_profiles/run.py

@@ -0,0 +1,95 @@
+# -*- coding: utf-8 -*-
+# SPDX-FileCopyrightText: : 2020-2024 The PyPSA-Eur Authors
+#
+# SPDX-License-Identifier: MIT
+
+import sys
+
+import numpy as np
+import pandas as pd
+import xarray as xr
+from _helpers import set_scenario_config
+from CentralHeatingCopApproximator import CentralHeatingCopApproximator
+from DecentralHeatingCopApproximator import DecentralHeatingCopApproximator
+
+from scripts.definitions.heat_system_type import HeatSystemType
+
+sys.path.append("..")
+
+
+def get_cop(
+    heat_system_type: str,
+    heat_source: str,
+    source_inlet_temperature_celsius: xr.DataArray,
+) -> xr.DataArray:
+    """
+    Calculate the coefficient of performance (COP) for a heating system.
+
+    Parameters
+    ----------
+    heat_system_type : str
+        The type of heating system.
+    heat_source : str
+        The heat source used in the heating system.
+    source_inlet_temperature_celsius : xr.DataArray
+        The inlet temperature of the heat source in Celsius.
+
+    Returns
+    -------
+    xr.DataArray
+        The calculated coefficient of performance (COP) for the heating system.
+    """
+    if HeatSystemType(heat_system_type).is_central:
+        return CentralHeatingCopApproximator(
+            forward_temperature_celsius=snakemake.params.forward_temperature_central_heating,
+            return_temperature_celsius=snakemake.params.return_temperature_central_heating,
+            source_inlet_temperature_celsius=source_inlet_temperature_celsius,
+            source_outlet_temperature_celsius=source_inlet_temperature_celsius
+            - snakemake.params.heat_source_cooling_central_heating,
+        ).approximate_cop()
+
+    else:
+        return DecentralHeatingCopApproximator(
+            forward_temperature_celsius=snakemake.params.heat_pump_sink_T_decentral_heating,
+            source_inlet_temperature_celsius=source_inlet_temperature_celsius,
+            source_type=heat_source,
+        ).approximate_cop()
+
+
+if __name__ == "__main__":
+    if "snakemake" not in globals():
+        from _helpers import mock_snakemake
+
+        snakemake = mock_snakemake(
+            "build_cop_profiles",
+            simpl="",
+            clusters=48,
+        )
+
+    set_scenario_config(snakemake)
+
+    cop_all_system_types = []
+    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"]
+            )
+            cop_da = get_cop(
+                heat_system_type=heat_system_type,
+                heat_source=heat_source,
+                source_inlet_temperature_celsius=source_inlet_temperature_celsius,
+            )
+            cop_this_system_type.append(cop_da)
+        cop_all_system_types.append(
+            xr.concat(
+                cop_this_system_type, dim=pd.Index(heat_sources, name="heat_source")
+            )
+        )
+
+    cop_dataarray = xr.concat(
+        cop_all_system_types,
+        dim=pd.Index(snakemake.params.heat_pump_sources.keys(), name="heat_system"),
+    )
+
+    cop_dataarray.to_netcdf(snakemake.output.cop_profiles)

scripts/definitions/heat_sector.py

@@ -0,0 +1,28 @@
+# -*- coding: utf-8 -*-
+# SPDX-FileCopyrightText: : 2020-2024 The PyPSA-Eur Authors
+#
+# SPDX-License-Identifier: MIT
+
+from enum import Enum
+
+
+class HeatSector(Enum):
+    """
+    Enumeration class representing different heat sectors.
+
+    Attributes:
+        RESIDENTIAL (str): Represents the residential heat sector.
+        SERVICES (str): Represents the services heat sector.
+    """
+
+    RESIDENTIAL = "residential"
+    SERVICES = "services"
+
+    def __str__(self) -> str:
+        """
+        Returns the string representation of the heat sector.
+
+        Returns:
+            str: The string representation of the heat sector.
+        """
+        return self.value

scripts/definitions/heat_system.py

@@ -0,0 +1,267 @@
+# -*- coding: utf-8 -*-
+# SPDX-FileCopyrightText: : 2020-2024 The PyPSA-Eur Authors
+#
+# SPDX-License-Identifier: MIT
+
+from enum import Enum
+
+from scripts.definitions.heat_sector import HeatSector
+from scripts.definitions.heat_system_type import HeatSystemType
+
+
+class HeatSystem(Enum):
+    """
+    Enumeration representing different heat systems.
+
+    Attributes
+    ----------
+    RESIDENTIAL_RURAL : str
+        Heat system for residential areas in rural locations.
+    SERVICES_RURAL : str
+        Heat system for service areas in rural locations.
+    RESIDENTIAL_URBAN_DECENTRAL : str
+        Heat system for residential areas in urban decentralized locations.
+    SERVICES_URBAN_DECENTRAL : str
+        Heat system for service areas in urban decentralized locations.
+    URBAN_CENTRAL : str
+        Heat system for urban central areas.
+
+    Methods
+    -------
+    __str__()
+        Returns the string representation of the heat system.
+    central_or_decentral()
+        Returns whether the heat system is central or decentralized.
+    system_type()
+        Returns the type of the heat system.
+    sector()
+        Returns the sector of the heat system.
+    rural()
+        Returns whether the heat system is for rural areas.
+    urban_decentral()
+        Returns whether the heat system is for urban decentralized areas.
+    urban()
+        Returns whether the heat system is for urban areas.
+    heat_demand_weighting(urban_fraction=None, dist_fraction=None)
+        Calculates the heat demand weighting based on urban fraction and distribution fraction.
+    heat_pump_costs_name(heat_source)
+        Generates the name for the heat pump costs based on the heat source.
+    """
+
+    RESIDENTIAL_RURAL = "residential rural"
+    SERVICES_RURAL = "services rural"
+    RESIDENTIAL_URBAN_DECENTRAL = "residential urban decentral"
+    SERVICES_URBAN_DECENTRAL = "services urban decentral"
+    URBAN_CENTRAL = "urban central"
+
+    def __init__(self, *args):
+        super().__init__(*args)
+
+    def __str__(self) -> str:
+        """
+        Returns the string representation of the heat system.
+
+        Returns
+        -------
+        str
+            The string representation of the heat system.
+        """
+        return self.value
+
+    @property
+    def central_or_decentral(self) -> str:
+        """
+        Returns whether the heat system is central or decentralized.
+
+        Returns
+        -------
+        str
+            "central" if the heat system is central, "decentral" otherwise.
+        """
+        if self == HeatSystem.URBAN_CENTRAL:
+            return "central"
+        else:
+            return "decentral"
+
+    @property
+    def system_type(self) -> HeatSystemType:
+        """
+        Returns the type of the heat system.
+
+        Returns
+        -------
+        str
+            The type of the heat system.
+
+        Raises
+        ------
+        RuntimeError
+            If the heat system is invalid.
+        """
+        if self == HeatSystem.URBAN_CENTRAL:
+            return HeatSystemType.URBAN_CENTRAL
+        elif (
+            self == HeatSystem.RESIDENTIAL_URBAN_DECENTRAL
+            or self == HeatSystem.SERVICES_URBAN_DECENTRAL
+        ):
+            return HeatSystemType.URBAN_DECENTRAL
+        elif self == HeatSystem.RESIDENTIAL_RURAL or self == HeatSystem.SERVICES_RURAL:
+            return HeatSystemType.RURAL
+        else:
+            raise RuntimeError(f"Invalid heat system: {self}")
+
+    @property
+    def sector(self) -> HeatSector:
+        """
+        Returns the sector of the heat system.
+
+        Returns
+        -------
+        HeatSector
+            The sector of the heat system.
+        """
+        if (
+            self == HeatSystem.RESIDENTIAL_RURAL
+            or self == HeatSystem.RESIDENTIAL_URBAN_DECENTRAL
+        ):
+            return HeatSector.RESIDENTIAL
+        elif (
+            self == HeatSystem.SERVICES_RURAL
+            or self == HeatSystem.SERVICES_URBAN_DECENTRAL
+        ):
+            return HeatSector.SERVICES
+        else:
+            "tot"
+
+    @property
+    def is_rural(self) -> bool:
+        """
+        Returns whether the heat system is for rural areas.
+
+        Returns
+        -------
+        bool
+            True if the heat system is for rural areas, False otherwise.
+        """
+        if self == HeatSystem.RESIDENTIAL_RURAL or self == HeatSystem.SERVICES_RURAL:
+            return True
+        else:
+            return False
+
+    @property
+    def is_urban_decentral(self) -> bool:
+        """
+        Returns whether the heat system is for urban decentralized areas.
+
+        Returns
+        -------
+        bool
+            True if the heat system is for urban decentralized areas, False otherwise.
+        """
+        if (
+            self == HeatSystem.RESIDENTIAL_URBAN_DECENTRAL
+            or self == HeatSystem.SERVICES_URBAN_DECENTRAL
+        ):
+            return True
+        else:
+            return False
+
+    @property
+    def is_urban(self) -> bool:
+        """
+        Returns whether the heat system is for urban areas.
+
+        Returns
+        -------
+        bool True if the heat system is for urban areas, False otherwise.
+        """
+        return not self.is_rural
+
+    def heat_demand_weighting(self, urban_fraction=None, dist_fraction=None) -> float:
+        """
+        Calculates the heat demand weighting based on urban fraction and
+        distribution fraction.
+
+        Parameters
+        ----------
+        urban_fraction : float, optional
+            The fraction of urban heat demand.
+        dist_fraction : float, optional
+            The fraction of distributed heat demand.
+
+        Returns
+        -------
+        float
+            The heat demand weighting.
+
+        Raises
+        ------
+        RuntimeError
+            If the heat system is invalid.
+        """
+        if "rural" in self.value:
+            return 1 - urban_fraction
+        elif "urban central" in self.value:
+            return dist_fraction
+        elif "urban decentral" in self.value:
+            return urban_fraction - dist_fraction
+        else:
+            raise RuntimeError(f"Invalid heat system: {self}")
+
+    def heat_pump_costs_name(self, heat_source: str) -> str:
+        """
+        Generates the name for the heat pump costs based on the heat source and
+        system.
+        Used to retrieve data from `technology-data <https://github.com/PyPSA/technology-data>`.
+
+        Parameters
+        ----------
+        heat_source : str
+            The heat source.
+
+        Returns
+        -------
+        str
+            The name for the heat pump costs.
+        """
+        return f"{self.central_or_decentral} {heat_source}-sourced heat pump"
+
+    @property
+    def resistive_heater_costs_name(self) -> str:
+        """
+        Generates the name for the resistive heater costs based on the heat
+        system.
+        Used to retrieve data from `technology-data <https://github.com/PyPSA/technology-data>`.
+
+        Returns
+        -------
+        str
+            The name for the heater costs.
+        """
+        return f"{self.central_or_decentral} resistive heater"
+
+    @property
+    def gas_boiler_costs_name(self) -> str:
+        """
+        Generates the name for the gas boiler costs based on the heat system.
+        Used to retrieve data from `technology-data <https://github.com/PyPSA/technology-data>`.
+
+        Returns
+        -------
+        str
+            The name for the gas boiler costs.
+        """
+        return f"{self.central_or_decentral} gas boiler"
+
+    @property
+    def oil_boiler_costs_name(self) -> str:
+        """
+        Generates the name for the oil boiler costs based on the heat system.
+        Used to retrieve data from `technology-data <https://github.com/PyPSA/technology-data>`.
+
+        Returns
+        -------
+        str
+            The name for the oil boiler costs.
+        """
+        return "decentral oil boiler"

scripts/definitions/heat_system_type.py

@@ -0,0 +1,35 @@
+# -*- coding: utf-8 -*-
+# SPDX-FileCopyrightText: : 2020-2024 The PyPSA-Eur Authors
+#
+# SPDX-License-Identifier: MIT
+
+from enum import Enum
+
+
+class HeatSystemType(Enum):
+    """
+    Enumeration representing different types of heat systems.
+    """
+
+    URBAN_CENTRAL = "urban central"
+    URBAN_DECENTRAL = "urban decentral"
+    RURAL = "rural"
+
+    def __str__(self) -> str:
+        """
+        Returns the string representation of the heat system type.
+
+        Returns:
+            str: The string representation of the heat system type.
+        """
+        return self.value
+
+    @property
+    def is_central(self) -> bool:
+        """
+        Returns whether the heat system type is central.
+
+        Returns:
+            bool: True if the heat system type is central, False otherwise.
+        """
+        return self == HeatSystemType.URBAN_CENTRAL

scripts/prepare_sector_network.py

@@ -37,6 +37,10 @@ from pypsa.geo import haversine_pts
 from pypsa.io import import_components_from_dataframe
 from scipy.stats import beta
 
+from scripts.definitions.heat_sector import HeatSector
+from scripts.definitions.heat_system import HeatSystem
+from scripts.definitions.heat_system_type import HeatSystemType
+
 spatial = SimpleNamespace()
 logger = logging.getLogger(__name__)
 
@@ -1776,7 +1780,7 @@ def build_heat_demand(n):
         .unstack(level=1)
     )
 
-    sectors = ["residential", "services"]
+    sectors = [sector.value for sector in HeatSector]
     uses = ["water", "space"]
 
     heat_demand = {}
@@ -1804,10 +1808,21 @@ def build_heat_demand(n):
     return heat_demand
 
 
-def add_heat(n, costs):
+def add_heat(n: pypsa.Network, costs: pd.DataFrame, cop: xr.DataArray):
+    """
+    Add heat sector to the network.
+
+    Parameters:
+        n (pypsa.Network): The PyPSA network object.
+        costs (pd.DataFrame): DataFrame containing cost information.
+        cop (xr.DataArray): DataArray containing coefficient of performance (COP) values.
+
+    Returns:
+        None
+    """
     logger.info("Add heat sector")
 
-    sectors = ["residential", "services"]
+    sectors = [sector.value for sector in HeatSector]
 
     heat_demand = build_heat_demand(n)
 
@@ -1826,23 +1841,6 @@ def add_heat(n, costs):
         for sector in sectors:
             heat_demand[sector + " space"] = (1 - dE) * heat_demand[sector + " space"]
 
-    heat_systems = [
-        "residential rural",
-        "services rural",
-        "residential urban decentral",
-        "services urban decentral",
-        "urban central",
-    ]
-
-    cop = {
-        "air": xr.open_dataarray(snakemake.input.cop_air_total)
-        .to_pandas()
-        .reindex(index=n.snapshots),
-        "ground": xr.open_dataarray(snakemake.input.cop_soil_total)
-        .to_pandas()
-        .reindex(index=n.snapshots),
-    }
-
     if options["solar_thermal"]:
         solar_thermal = (
             xr.open_dataarray(snakemake.input.solar_thermal_total)
@@ -1852,31 +1850,34 @@ def add_heat(n, costs):
         # 1e3 converts from W/m^2 to MW/(1000m^2) = kW/m^2
         solar_thermal = options["solar_cf_correction"] * solar_thermal / 1e3
 
-    for name in heat_systems:
+    for (
-        name_type = "central" if name == "urban central" else "decentral"
+        heat_system
+    ) in (
+        HeatSystem
+    ):  # this loops through all heat systems defined in _entities.HeatSystem
 
-        if name == "urban central":
+        if heat_system == HeatSystem.URBAN_CENTRAL:
             nodes = dist_fraction.index[dist_fraction > 0]
         else:
             nodes = pop_layout.index
 
-        n.add("Carrier", name + " heat")
+        n.add("Carrier", f"{heat_system} heat")
 
         n.madd(
             "Bus",
-            nodes + f" {name} heat",
+            nodes + f" {heat_system.value} heat",
             location=nodes,
-            carrier=name + " heat",
+            carrier=f"{heat_system.value} heat",
             unit="MWh_th",
         )
 
-        if name == "urban central" and options.get("central_heat_vent"):
+        if heat_system == HeatSystem.URBAN_CENTRAL and options.get("central_heat_vent"):
             n.madd(
                 "Generator",
-                nodes + f" {name} heat vent",
+                nodes + f" {heat_system} heat vent",
-                bus=nodes + f" {name} heat",
+                bus=nodes + f" {heat_system} heat",
                 location=nodes,
-                carrier=name + " heat vent",
+                carrier=f"{heat_system} heat vent",
                 p_nom_extendable=True,
                 p_max_pu=0,
                 p_min_pu=-1,
@@ -1884,30 +1885,24 @@ def add_heat(n, costs):
             )
 
         ## Add heat load
+        factor = heat_system.heat_demand_weighting(
+            urban_fraction=urban_fraction[nodes], dist_fraction=dist_fraction[nodes]
+        )
+        if not heat_system == HeatSystem.URBAN_CENTRAL:
+            heat_load = (
+                heat_demand[
+                    [
+                        heat_system.sector.value + " water",
+                        heat_system.sector.value + " space",
+                    ]
+                ]
+                .T.groupby(level=1)
+                .sum()
+                .T[nodes]
+                .multiply(factor)
+            )
 
-        for sector in sectors:
+        if heat_system == HeatSystem.URBAN_CENTRAL:
-            # heat demand weighting
-            if "rural" in name:
-                factor = 1 - urban_fraction[nodes]
-            elif "urban central" in name:
-                factor = dist_fraction[nodes]
-            elif "urban decentral" in name:
-                factor = urban_fraction[nodes] - dist_fraction[nodes]
-            else:
-                raise NotImplementedError(
-                    f" {name} not in " f"heat systems: {heat_systems}"
-                )
-
-            if sector in name:
-                heat_load = (
-                    heat_demand[[sector + " water", sector + " space"]]
-                    .T.groupby(level=1)
-                    .sum()
-                    .T[nodes]
-                    .multiply(factor)
-                )
-
-        if name == "urban central":
             heat_load = (
                 heat_demand.T.groupby(level=1)
                 .sum()
@@ -1920,20 +1915,25 @@ def add_heat(n, costs):
         n.madd(
             "Load",
             nodes,
-            suffix=f" {name} heat",
+            suffix=f" {heat_system} heat",
-            bus=nodes + f" {name} heat",
+            bus=nodes + f" {heat_system} heat",
-            carrier=name + " heat",
+            carrier=f"{heat_system} heat",
             p_set=heat_load,
         )
 
         ## Add heat pumps
-
+        for heat_source in snakemake.params.heat_pump_sources[
-        heat_pump_types = ["air"] if "urban" in name else ["ground", "air"]
+            heat_system.system_type.value
-
+        ]:
-        for heat_pump_type in heat_pump_types:
+            costs_name = heat_system.heat_pump_costs_name(heat_source)
-            costs_name = f"{name_type} {heat_pump_type}-sourced heat pump"
             efficiency = (
-                cop[heat_pump_type][nodes]
+                cop.sel(
+                    heat_system=heat_system.system_type.value,
+                    heat_source=heat_source,
+                    name=nodes,
+                )
+                .to_pandas()
+                .reindex(index=n.snapshots)
                 if options["time_dep_hp_cop"]
                 else costs.at[costs_name, "efficiency"]
             )
@@ -1941,10 +1941,10 @@ def add_heat(n, costs):
             n.madd(
                 "Link",
                 nodes,
-                suffix=f" {name} {heat_pump_type} heat pump",
+                suffix=f" {heat_system} {heat_source} heat pump",
                 bus0=nodes,
-                bus1=nodes + f" {name} heat",
+                bus1=nodes + f" {heat_system} heat",
-                carrier=f"{name} {heat_pump_type} heat pump",
+                carrier=f"{heat_system} {heat_source} heat pump",
                 efficiency=efficiency,
                 capital_cost=costs.at[costs_name, "efficiency"]
                 * costs.at[costs_name, "fixed"]
@@ -1954,59 +1954,65 @@ def add_heat(n, costs):
             )
 
         if options["tes"]:
-            n.add("Carrier", name + " water tanks")
+            n.add("Carrier", f"{heat_system} water tanks")
 
             n.madd(
                 "Bus",
-                nodes + f" {name} water tanks",
+                nodes + f" {heat_system} water tanks",
                 location=nodes,
-                carrier=name + " water tanks",
+                carrier=f"{heat_system} water tanks",
                 unit="MWh_th",
             )
 
             n.madd(
                 "Link",
-                nodes + f" {name} water tanks charger",
+                nodes + f" {heat_system} water tanks charger",
-                bus0=nodes + f" {name} heat",
+                bus0=nodes + f" {heat_system} heat",
-                bus1=nodes + f" {name} water tanks",
+                bus1=nodes + f" {heat_system} water tanks",
                 efficiency=costs.at["water tank charger", "efficiency"],
-                carrier=name + " water tanks charger",
+                carrier=f"{heat_system} water tanks charger",
                 p_nom_extendable=True,
             )
 
             n.madd(
                 "Link",
-                nodes + f" {name} water tanks discharger",
+                nodes + f" {heat_system} water tanks discharger",
-                bus0=nodes + f" {name} water tanks",
+                bus0=nodes + f" {heat_system} water tanks",
-                bus1=nodes + f" {name} heat",
+                bus1=nodes + f" {heat_system} heat",
-                carrier=name + " water tanks discharger",
+                carrier=f"{heat_system} water tanks discharger",
                 efficiency=costs.at["water tank discharger", "efficiency"],
                 p_nom_extendable=True,
             )
 
-            tes_time_constant_days = options["tes_tau"][name_type]
+            tes_time_constant_days = options["tes_tau"][
+                heat_system.central_or_decentral
+            ]
 
             n.madd(
                 "Store",
-                nodes + f" {name} water tanks",
+                nodes + f" {heat_system} water tanks",
-                bus=nodes + f" {name} water tanks",
+                bus=nodes + f" {heat_system} water tanks",
                 e_cyclic=True,
                 e_nom_extendable=True,
-                carrier=name + " water tanks",
+                carrier=f"{heat_system} water tanks",
                 standing_loss=1 - np.exp(-1 / 24 / tes_time_constant_days),
-                capital_cost=costs.at[name_type + " water tank storage", "fixed"],
+                capital_cost=costs.at[
-                lifetime=costs.at[name_type + " water tank storage", "lifetime"],
+                    heat_system.central_or_decentral + " water tank storage", "fixed"
+                ],
+                lifetime=costs.at[
+                    heat_system.central_or_decentral + " water tank storage", "lifetime"
+                ],
             )
 
         if options["resistive_heaters"]:
-            key = f"{name_type} resistive heater"
+            key = f"{heat_system.central_or_decentral} resistive heater"
 
             n.madd(
                 "Link",
-                nodes + f" {name} resistive heater",
+                nodes + f" {heat_system} resistive heater",
                 bus0=nodes,
-                bus1=nodes + f" {name} heat",
+                bus1=nodes + f" {heat_system} heat",
-                carrier=name + " resistive heater",
+                carrier=f"{heat_system} resistive heater",
                 efficiency=costs.at[key, "efficiency"],
                 capital_cost=costs.at[key, "efficiency"]
                 * costs.at[key, "fixed"]
@@ -2016,16 +2022,16 @@ def add_heat(n, costs):
             )
 
         if options["boilers"]:
-            key = f"{name_type} gas boiler"
+            key = f"{heat_system.central_or_decentral} gas boiler"
 
             n.madd(
                 "Link",
-                nodes + f" {name} gas boiler",
+                nodes + f" {heat_system} gas boiler",
                 p_nom_extendable=True,
                 bus0=spatial.gas.df.loc[nodes, "nodes"].values,
-                bus1=nodes + f" {name} heat",
+                bus1=nodes + f" {heat_system} heat",
                 bus2="co2 atmosphere",
-                carrier=name + " gas boiler",
+                carrier=f"{heat_system} gas boiler",
                 efficiency=costs.at[key, "efficiency"],
                 efficiency2=costs.at["gas", "CO2 intensity"],
                 capital_cost=costs.at[key, "efficiency"]
@@ -2035,22 +2041,26 @@ def add_heat(n, costs):
             )
 
         if options["solar_thermal"]:
-            n.add("Carrier", name + " solar thermal")
+            n.add("Carrier", f"{heat_system} solar thermal")
 
             n.madd(
                 "Generator",
                 nodes,
-                suffix=f" {name} solar thermal collector",
+                suffix=f" {heat_system} solar thermal collector",
-                bus=nodes + f" {name} heat",
+                bus=nodes + f" {heat_system} heat",
-                carrier=name + " solar thermal",
+                carrier=f"{heat_system} solar thermal",
                 p_nom_extendable=True,
-                capital_cost=costs.at[name_type + " solar thermal", "fixed"]
+                capital_cost=costs.at[
+                    heat_system.central_or_decentral + " solar thermal", "fixed"
+                ]
                 * overdim_factor,
                 p_max_pu=solar_thermal[nodes],
-                lifetime=costs.at[name_type + " solar thermal", "lifetime"],
+                lifetime=costs.at[
+                    heat_system.central_or_decentral + " solar thermal", "lifetime"
+                ],
             )
 
-        if options["chp"] and name == "urban central":
+        if options["chp"] and heat_system == HeatSystem.URBAN_CENTRAL:
             # add gas CHP; biomass CHP is added in biomass section
             n.madd(
                 "Link",
@@ -2107,16 +2117,20 @@ def add_heat(n, costs):
                 lifetime=costs.at["central gas CHP", "lifetime"],
             )
 
-        if options["chp"] and options["micro_chp"] and name != "urban central":
+        if (
+            options["chp"]
+            and options["micro_chp"]
+            and heat_system.value != "urban central"
+        ):
             n.madd(
                 "Link",
-                nodes + f" {name} micro gas CHP",
+                nodes + f" {heat_system} micro gas CHP",
                 p_nom_extendable=True,
                 bus0=spatial.gas.df.loc[nodes, "nodes"].values,
                 bus1=nodes,
-                bus2=nodes + f" {name} heat",
+                bus2=nodes + f" {heat_system} heat",
                 bus3="co2 atmosphere",
-                carrier=name + " micro gas CHP",
+                carrier=heat_system.value + " micro gas CHP",
                 efficiency=costs.at["micro CHP", "efficiency"],
                 efficiency2=costs.at["micro CHP", "efficiency-heat"],
                 efficiency3=costs.at["gas", "CO2 intensity"],
@@ -2152,7 +2166,7 @@ def add_heat(n, costs):
         ) / heat_demand.T.groupby(level=[1]).sum().T
 
         for name in n.loads[
-            n.loads.carrier.isin([x + " heat" for x in heat_systems])
+            n.loads.carrier.isin([x + " heat" for x in HeatSystem])
         ].index:
             node = n.buses.loc[name, "location"]
             ct = pop_layout.loc[node, "ct"]
@@ -3110,27 +3124,23 @@ def add_industry(n, costs):
     if options["oil_boilers"]:
         nodes = pop_layout.index
 
-        for name in [
+        for heat_system in HeatSystem:
-            "residential rural",
+            if not heat_system == HeatSystem.URBAN_CENTRAL:
-            "services rural",
+                n.madd(
-            "residential urban decentral",
+                    "Link",
-            "services urban decentral",
+                    nodes + f" {heat_system} oil boiler",
-        ]:
+                    p_nom_extendable=True,
-            n.madd(
+                    bus0=spatial.oil.nodes,
-                "Link",
+                    bus1=nodes + f" {heat_system} heat",
-                nodes + f" {name} oil boiler",
+                    bus2="co2 atmosphere",
-                p_nom_extendable=True,
+                    carrier=f"{heat_system} oil boiler",
-                bus0=spatial.oil.nodes,
+                    efficiency=costs.at["decentral oil boiler", "efficiency"],
-                bus1=nodes + f" {name} heat",
+                    efficiency2=costs.at["oil", "CO2 intensity"],
-                bus2="co2 atmosphere",
+                    capital_cost=costs.at["decentral oil boiler", "efficiency"]
-                carrier=f"{name} oil boiler",
+                    * costs.at["decentral oil boiler", "fixed"]
-                efficiency=costs.at["decentral oil boiler", "efficiency"],
+                    * options["overdimension_individual_heating"],
-                efficiency2=costs.at["oil", "CO2 intensity"],
+                    lifetime=costs.at["decentral oil boiler", "lifetime"],
-                capital_cost=costs.at["decentral oil boiler", "efficiency"]
+                )
-                * costs.at["decentral oil boiler", "fixed"]
-                * options["overdimension_individual_heating"],
-                lifetime=costs.at["decentral oil boiler", "lifetime"],
-            )
 
     n.madd(
         "Link",
@@ -4186,7 +4196,7 @@ if __name__ == "__main__":
         add_land_transport(n, costs)
 
     if options["heating"]:
-        add_heat(n, costs)
+        add_heat(n=n, costs=costs, cop=xr.open_dataarray(snakemake.input.cop_profiles))
 
     if options["biomass"]:
         add_biomass(n, costs)