feat: layered ptes model

config/config.default.yaml

@@ -665,6 +665,9 @@ sector:
       bottom_temperature: 35
       design_top_temperature: 90
       design_bottom_temperature: 35
+      layered:
+        num_layers: 1 
+        heat_transfer_coefficient: 0.01  # H [W/(m²·K)]
     ates:
       enable: false
       suitable_aquifer_types:

data/custom_costs.csv

@@ -10,4 +10,5 @@ all,battery,marginal_cost,0,EUR/MWh,Default value to prevent mathematical degene
 all,battery inverter,marginal_cost,0,EUR/MWh,Default value to prevent mathematical degeneracy in optimisation,
 all,home battery storage,marginal_cost,0,EUR/MWh,Default value to prevent mathematical degeneracy in optimisation,
 all,water tank charger,marginal_cost,0.03,EUR/MWh,Default value to prevent mathematical degeneracy in optimisation,
-all,central water pit charger,marginal_cost,0.025,EUR/MWh,Default value to prevent mathematical degeneracy in optimisation,
+all,central water pit charger,marginal_cost,0.04,EUR/MWh,Default value to prevent mathematical degeneracy in optimisation,
+all,central water pit discharger,marginal_cost,0.05,EUR/MWh,Default value to prevent mathematical degeneracy in optimisation,

rules/build_sector.smk

@@ -529,7 +529,6 @@ def input_heat_source_temperature(
     replace_names: dict[str, str] = {
         "air": "air_total",
         "ground": "soil_total",
-        "ptes": "ptes_top_profiles",
     },
 ) -> dict[str, str]:
     """
@@ -549,7 +548,7 @@ def input_heat_source_temperature(
         for heat sources that require temperature profiles (excludes constant
         temperature sources).
     """
-    from scripts.definitions.heat_source import HeatSource
+    from scripts.definitions.heat_source import HeatSource, HeatSourceType
 
     heat_sources = set(
         config_provider("sector", "heat_sources", "urban central")(w)
@@ -562,25 +561,23 @@ def input_heat_source_temperature(
     for heat_source_name in heat_sources:
         heat_source = HeatSource(heat_source_name)
         # Skip heat sources with temperatures defined in config (not from file)
-        if heat_source.temperature_from_config:
+        if (
+            heat_source.temperature_from_config
+            or heat_source.source_type == HeatSourceType.STORAGE
+        ):
             continue
-        if heat_source_name == "ptes":
-            file_names[f"temp_{heat_source_name}"] = resources(
-                f"temp_{replace_names.get(heat_source_name, heat_source_name)}_base_s_{{clusters}}_{{planning_horizons}}.nc"
-            )
         else:
             file_names[f"temp_{heat_source_name}"] = resources(
                 f"temp_{replace_names.get(heat_source_name, heat_source_name)}_base_s_{{clusters}}.nc"
             )
     return file_names
 
 
-def input_ptes_bottom_temperature(w) -> dict[str, str]:
+def input_ptes_operations(w) -> dict[str, str]:
     """
-    Generate conditional input for PTES bottom temperature profiles.
+    Generate conditional input for PTES operations.
 
-    Only includes the input file if PTES is configured as a heat source
-    for urban central heating.
+    Only includes the input file if PTES is enabled.
 
     Parameters
     ----------
@@ -593,11 +590,10 @@ def input_ptes_bottom_temperature(w) -> dict[str, str]:
         Dictionary with "temp_ptes_bottom" key if PTES is a heat source,
         empty dict otherwise.
     """
-    heat_sources = config_provider("sector", "heat_sources", "urban central")(w)
-    if "ptes" in heat_sources:
+    if config_provider("sector", "district_heating", "ptes", "enable")(w):
         return {
-            "temp_ptes_bottom": resources(
-                "temp_ptes_bottom_profiles_base_s_{clusters}_{planning_horizons}.nc"
+            "ptes_operations": resources(
+                "ptes_operations_base_s_{clusters}_{planning_horizons}.nc"
             )
         }
     return {}
@@ -678,9 +674,12 @@ rule build_cop_profiles:
         heat_source_temperatures=config_provider(
             "sector", "district_heating", "heat_source_temperatures"
         ),
+        ptes_enable=config_provider("sector", "district_heating", "ptes", "enable"),
+        ptes_layered=config_provider("sector", "district_heating", "ptes", "layered"),
         snapshots=config_provider("snapshots"),
     input:
         unpack(input_heat_source_temperature),
+        unpack(input_ptes_operations),
         central_heating_forward_temperature_profiles=resources(
             "central_heating_forward_temperature_profiles_base_s_{clusters}_{planning_horizons}.nc"
         ),
@@ -738,6 +737,7 @@ rule build_ptes_operations:
             "ptes",
             "design_bottom_temperature",
         ),
+        layered=config_provider("sector", "district_heating", "ptes", "layered"),
     input:
         central_heating_forward_temperature_profiles=resources(
             "central_heating_forward_temperature_profiles_base_s_{clusters}_{planning_horizons}.nc"
@@ -747,17 +747,8 @@ rule build_ptes_operations:
         ),
         regions_onshore=resources("regions_onshore_base_s_{clusters}.geojson"),
     output:
-        ptes_top_temperature_profiles=resources(
-            "temp_ptes_top_profiles_base_s_{clusters}_{planning_horizons}.nc"
-        ),
-        ptes_bottom_temperature_profiles=resources(
-            "temp_ptes_bottom_profiles_base_s_{clusters}_{planning_horizons}.nc"
-        ),
-        ptes_e_max_pu_profiles=resources(
-            "ptes_e_max_pu_profiles_base_s_{clusters}_{planning_horizons}.nc"
-        ),
-        ptes_boost_per_discharge_profiles=resources(
-            "ptes_boost_per_discharge_profiles_base_s_{clusters}_{planning_horizons}.nc"
+        ptes_operations=resources(
+            "ptes_operations_base_s_{clusters}_{planning_horizons}.nc"
         ),
     resources:
         mem_mb=2000,
@@ -784,7 +775,7 @@ rule build_heat_source_utilisation_profiles:
         ptes_enable=config_provider("sector", "district_heating", "ptes", "enable"),
     input:
         unpack(input_heat_source_temperature),
-        unpack(input_ptes_bottom_temperature),
+        unpack(input_ptes_operations),
         central_heating_forward_temperature_profiles=resources(
             "central_heating_forward_temperature_profiles_base_s_{clusters}_{planning_horizons}.nc"
         ),
@@ -1620,6 +1611,7 @@ rule prepare_sector_network:
     input:
         unpack(input_profile_offwind),
         unpack(input_heat_source_power),
+        unpack(input_ptes_operations),
         **rules.cluster_gas_network.output,
         **rules.build_gas_input_locations.output,
         snapshot_weightings=resources(
@@ -1692,26 +1684,6 @@ rule prepare_sector_network:
         temp_soil_total=resources("temp_soil_total_base_s_{clusters}.nc"),
         temp_air_total=resources("temp_air_total_base_s_{clusters}.nc"),
         cop_profiles=resources("cop_profiles_base_s_{clusters}_{planning_horizons}.nc"),
-        ptes_e_max_pu_profiles=lambda w: (
-            resources(
-                "ptes_e_max_pu_profiles_base_s_{clusters}_{planning_horizons}.nc"
-            )
-            if config_provider("sector", "district_heating", "ptes", "enable")(w)
-            and config_provider(
-                "sector", "district_heating", "ptes", "temperature_dependent_capacity"
-            )(w)
-            else []
-        ),
-        ptes_boost_per_discharge_profiles=lambda w: (
-            resources(
-                "ptes_boost_per_discharge_profiles_base_s_{clusters}_{planning_horizons}.nc"
-            )
-            if config_provider("sector", "district_heating", "ptes", "enable")(w)
-            and config_provider(
-                "sector", "district_heating", "ptes", "discharge_resistive_boosting"
-            )(w)
-            else []
-        ),
         heat_source_direct_utilisation_profiles=lambda w: (
             resources(
                 "heat_source_direct_utilisation_profiles_base_s_{clusters}_{planning_horizons}.nc"

rules/solve_myopic.smk

@@ -116,6 +116,7 @@ rule solve_sector_network_myopic:
         ),
         custom_extra_functionality=input_custom_extra_functionality,
     input:
+        unpack(input_ptes_operations),
         network=resources(
             "networks/base_s_{clusters}_{opts}_{sector_opts}_{planning_horizons}_brownfield.nc"
         ),

rules/solve_overnight.smk

@@ -14,6 +14,7 @@ rule solve_sector_network:
         ),
         custom_extra_functionality=input_custom_extra_functionality,
     input:
+        unpack(input_ptes_operations),
         network=resources(
             "networks/base_s_{clusters}_{opts}_{sector_opts}_{planning_horizons}.nc"
         ),

rules/solve_perfect.smk

@@ -101,6 +101,7 @@ rule solve_sector_network_perfect:
         ),
         custom_extra_functionality=input_custom_extra_functionality,
     input:
+        unpack(input_ptes_operations),
         network=resources(
             "networks/base_s_{clusters}_{opts}_{sector_opts}_brownfield_all_years.nc"
         ),

scripts/build_cop_profiles/run.py

@@ -65,7 +65,7 @@
 from scripts.build_cop_profiles.decentral_heating_cop_approximator import (
     DecentralHeatingCopApproximator,
 )
-from scripts.definitions.heat_source import HeatSource
+from scripts.definitions.heat_source import HeatSource, HeatSourceType
 from scripts.definitions.heat_system_type import HeatSystemType
 
 
@@ -167,6 +167,13 @@ def get_source_temperature(
     heat_source = HeatSource(heat_source_name)
     if heat_source.temperature_from_config:
         return snakemake_params["heat_source_temperatures"][heat_source_name]
+    elif heat_source.source_type == HeatSourceType.STORAGE:
+        # PTES layer temperatures are constants from the ptes_operations dataset
+        if heat_source_name.startswith("ptes layer"):
+            layer_idx = int(heat_source_name.split()[-1])
+            return float(ptes_ds["layer_temperatures"].sel(layer=layer_idx).item())
+        else:
+            return float(ptes_ds.attrs["top_temperature"])
     else:
         return xr.open_dataarray(snakemake_input[f"temp_{heat_source_name}"])
 
@@ -268,6 +275,14 @@ def get_sink_inlet_temperature(
         snakemake.input.central_heating_return_temperature_profiles
     )
 
+    # Load PTES operations dataset if enabled
+    if snakemake.params.ptes_enable:
+        ptes_ds = xr.open_dataset(snakemake.input.ptes_operations)
+        num_ptes_layers = int(ptes_ds.attrs["num_layers"])
+    else:
+        ptes_ds = None
+        num_ptes_layers = 0
+
     cop_all_system_types = []
     for heat_system_type, heat_sources in snakemake.params.heat_sources.items():
         cop_this_system_type = []
@@ -309,6 +324,9 @@ def get_sink_inlet_temperature(
             )
         )
 
+    if ptes_ds is not None:
+        ptes_ds.close()
+
     cop_dataarray = xr.concat(
         cop_all_system_types,
         dim=pd.Index(snakemake.params.heat_sources.keys(), name="heat_system"),

scripts/build_heat_source_utilisation_profiles.py

@@ -61,7 +61,7 @@
 import xarray as xr
 
 from scripts._helpers import configure_logging, set_scenario_config
-from scripts.definitions.heat_source import HeatSource
+from scripts.definitions.heat_source import HeatSource, HeatSourceType
 
 logger = logging.getLogger(__name__)
 
@@ -96,6 +96,13 @@ def get_source_temperature(
     heat_source = HeatSource(heat_source_name)
     if heat_source.temperature_from_config:
         return snakemake_params["heat_source_temperatures"][heat_source_name]
+    elif heat_source.source_type == HeatSourceType.STORAGE:
+        # PTES layer temperatures are constants from the ptes_operations dataset
+        if heat_source_name.startswith("ptes layer"):
+            layer_idx = int(heat_source_name.split()[-1])
+            return float(ptes_ds["layer_temperatures"].sel(layer=layer_idx).item())
+        else:
+            return float(ptes_ds.attrs["top_temperature"])
     else:
         if f"temp_{heat_source_name}" not in snakemake_input.keys():
             raise ValueError(
@@ -182,50 +189,38 @@ def get_preheater_utilisation_profile(
 def get_heat_pump_cooling(
     heat_source_name: str,
     default_heat_source_cooling: float,
-    snakemake_input: dict,
     return_temperature: xr.DataArray = None,
 ) -> float | xr.DataArray:
     """
     Get the additional heat source cooling (temperature drop) through the heat pump for a heat source.
 
     For PTES, this equals the temperature difference between
-    return flow and bottom layers (return_temperature - bottom_temperature), which can
-    be time-varying. For other sources, uses the default constant value.
+    return flow and bottom layer (return_temperature - bottom_temperature).
+    For other sources, uses the default constant value.
 
     Parameters
     ----------
     heat_source_name : str
         Name of the heat source (e.g., 'ptes', 'geothermal', 'air').
     default_heat_source_cooling : float
         Default heat source cooling in Kelvin, from config.
-    snakemake_input : dict
-        Snakemake input files, may contain PTES temperature profiles.
     return_temperature : xr.DataArray, optional
         District heating return temperature profiles in °C. Required for PTES.
 
     Returns
     -------
     float | xr.DataArray
         Heat source cooling in Kelvin. Returns a float for most sources,
-        or a DataArray for PTES when temperatures vary with time.
-
-    Raises
-    ------
-    ValueError
-        If heat source is PTES but bottom temperature profile is not provided.
+        or a DataArray for PTES.
     """
-    if heat_source_name == "ptes":
-        if "temp_ptes_bottom" not in snakemake_input.keys():
-            raise ValueError(
-                "PTES heat source requires bottom temperature profile "
-                "(temp_ptes_bottom) to calculate heat source cooling."
-            )
+    heat_source = HeatSource(heat_source_name)
+    if heat_source.source_type == HeatSourceType.STORAGE:
         if return_temperature is None:
             raise ValueError(
                 "PTES heat source requires return_temperature to calculate heat pump cooling."
             )
-        ptes_bottom_temperature = xr.open_dataarray(snakemake_input["temp_ptes_bottom"])
-        return return_temperature - ptes_bottom_temperature
+        bottom_temperature = float(ptes_ds.attrs["bottom_temperature"])
+        return return_temperature - bottom_temperature
     return default_heat_source_cooling
 
 
@@ -234,7 +229,7 @@ def get_heat_pump_cooling(
         from scripts._helpers import mock_snakemake
 
         snakemake = mock_snakemake(
-            "build_cop_profiles",
+            "build_heat_source_utilisation_profiles",
             clusters=48,
         )
     configure_logging(snakemake)
@@ -243,12 +238,13 @@ def get_heat_pump_cooling(
     heat_sources: list[str] = snakemake.params.heat_sources
     ptes_enable: bool = snakemake.params.ptes_enable
 
-    # Validate PTES configuration
-    if ptes_enable and "ptes" not in heat_sources:
-        raise ValueError(
-            "PTES is enabled (district_heating.ptes.enable=true) but 'ptes' "
-            "is not in heat_sources.urban_central. PTES requires being listed in heat_sources to create the necessary buses and links for heat discharge to the 'urban central heat' bus."
-        )
+    # Load PTES operations dataset if enabled
+    if ptes_enable:
+        ptes_ds = xr.open_dataset(snakemake.input.ptes_operations)
+        num_ptes_layers = int(ptes_ds.attrs["num_layers"])
+    else:
+        ptes_ds = None
+        num_ptes_layers = 0
 
     central_heating_forward_temperature: xr.DataArray = xr.open_dataarray(
         snakemake.input.central_heating_forward_temperature_profiles
@@ -285,11 +281,13 @@ def get_heat_pump_cooling(
                 heat_source_cooling=get_heat_pump_cooling(
                     heat_source_name=heat_source_key,
                     default_heat_source_cooling=snakemake.params.heat_source_cooling,
-                    snakemake_input=snakemake.input,
                     return_temperature=central_heating_return_temperature,
                 ),
             ).assign_coords(heat_source=heat_source_key)
             for heat_source_key in heat_sources
         ],
         dim="heat_source",
     ).to_netcdf(snakemake.output.heat_source_preheater_utilisation_profiles)
+
+    if ptes_ds is not None:
+        ptes_ds.close()