diff --git a/config/config.default.yaml b/config/config.default.yaml
index bbdb363f72..960a6e00e7 100644
--- a/config/config.default.yaml
+++ b/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:
diff --git a/data/custom_costs.csv b/data/custom_costs.csv
index 312cd0b3da..d7b71e3bab 100644
--- a/data/custom_costs.csv
+++ b/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,
diff --git a/rules/build_sector.smk b/rules/build_sector.smk
index 66d591f84a..29143389c2 100755
--- a/rules/build_sector.smk
+++ b/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,12 +561,11 @@ 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"
@@ -575,12 +573,11 @@ def input_heat_source_temperature(
     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"
diff --git a/rules/solve_myopic.smk b/rules/solve_myopic.smk
index d90c332b5d..03f3289eae 100644
--- a/rules/solve_myopic.smk
+++ b/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"
         ),
diff --git a/rules/solve_overnight.smk b/rules/solve_overnight.smk
index 270609ecf2..c0da468064 100644
--- a/rules/solve_overnight.smk
+++ b/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"
         ),
diff --git a/rules/solve_perfect.smk b/rules/solve_perfect.smk
index 59ec549254..4234bbdd0a 100644
--- a/rules/solve_perfect.smk
+++ b/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"
         ),
diff --git a/scripts/build_cop_profiles/run.py b/scripts/build_cop_profiles/run.py
index 6af8b02e8a..c16581138d 100644
--- a/scripts/build_cop_profiles/run.py
+++ b/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"),
diff --git a/scripts/build_heat_source_utilisation_profiles.py b/scripts/build_heat_source_utilisation_profiles.py
index 4d18185fcf..3c92c2833c 100644
--- a/scripts/build_heat_source_utilisation_profiles.py
+++ b/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,15 +189,14 @@ 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
     ----------
@@ -198,8 +204,6 @@ def get_heat_pump_cooling(
         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.
 
@@ -207,25 +211,16 @@ def get_heat_pump_cooling(
     -------
     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,7 +281,6 @@ 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)
@@ -293,3 +288,6 @@ def get_heat_pump_cooling(
         ],
         dim="heat_source",
     ).to_netcdf(snakemake.output.heat_source_preheater_utilisation_profiles)
+
+    if ptes_ds is not None:
+        ptes_ds.close()
diff --git a/scripts/build_ptes_operations/ptes_approximator.py b/scripts/build_ptes_operations/ptes_approximator.py
new file mode 100644
index 0000000000..05a7a43e47
--- /dev/null
+++ b/scripts/build_ptes_operations/ptes_approximator.py
@@ -0,0 +1,217 @@
+# SPDX-FileCopyrightText: Contributors to PyPSA-Eur <https://github.com/pypsa/pypsa-eur>
+#
+# SPDX-License-Identifier: MIT
+
+import numpy as np
+import xarray as xr
+
+
+class PtesApproximator:
+    def __init__(
+        self,
+        forward_temperature: xr.DataArray,
+        return_temperature: xr.DataArray,
+        top_temperature: float,
+        bottom_temperature: float,
+        num_layers: int,
+        design_top_temperature: float,
+        design_bottom_temperature: float,
+        design_standing_losses: float,
+        interlayer_heat_transfer_coefficient: float,
+    ):
+        self.forward_temperature = forward_temperature
+        self.return_temperature = return_temperature
+        self.top_temperature = top_temperature
+        self.bottom_temperature = bottom_temperature
+        self.num_layers = num_layers
+        self.design_top_temperature = design_top_temperature
+        self.design_bottom_temperature = design_bottom_temperature
+        self.design_standing_losses = design_standing_losses
+        self.interlayer_heat_transfer_coefficient = interlayer_heat_transfer_coefficient
+
+        self._time_dim = [d for d in forward_temperature.dims if d != "name"][0]
+
+    @property
+    def layer_temperatures(self) -> np.ndarray:
+        """Calculate the fixed temperature of each layer [°C], hottest first."""
+        return np.linspace(
+            self.top_temperature, self.bottom_temperature, self.num_layers + 1
+        )[:-1]
+
+    @property
+    def volume_weights(self) -> np.ndarray:
+        """
+        Volume weight W_l per layer.
+
+        W_l = (T_top - T_bottom) / (T_l - T_bottom).
+        A unit of energy in a colder layer occupies more volume.
+        W_1 = 1 for the hottest layer; W_l ≥ 1 for colder layers.
+        """
+        weights = (self.design_top_temperature - self.design_bottom_temperature) / (
+            self.layer_temperatures - self.bottom_temperature
+        )
+        return weights
+
+    @property
+    def charging_availability(self) -> xr.DataArray:
+        """
+        Binary charging availability Φ⁺_{l,t}.
+
+        For each timestep, only the layer whose temperature is closest to the forward temperature can receive charge. Returns a DataArray with dimensions (snapshot, name, layer).
+        """
+        # |T_l - T_fwd_t| for each layer
+        # Broadcast: T_fwd is (snapshot, name), T_layers is (layer,)
+        layer_dim = xr.DataArray(
+            self.layer_temperatures,
+            dims=["layer"],
+            coords={"layer": np.arange(len(self.layer_temperatures))},
+        )
+        abs_diff = np.abs(
+            layer_dim - self.forward_temperature
+        )  # (snapshot, name, layer)
+
+        # Find the layer index with minimum distance at each (snapshot, name)
+        closest_layer = abs_diff.argmin(dim="layer")  # (snapshot, name)
+
+        # Build binary availability: 1 where layer == closest_layer
+        layer_indices = xr.DataArray(
+            np.arange(len(self.layer_temperatures)),
+            dims=["layer"],
+            coords={"layer": np.arange(len(self.layer_temperatures))},
+        )
+
+        binary_availability = (layer_indices == closest_layer).astype(float)
+        binary_availability = binary_availability.transpose(
+            self._time_dim, "name", "layer"
+        )
+
+        return binary_availability
+
+    @property
+    def interlayer_heat_transfer_coefficients(self) -> np.ndarray:
+        """."""
+        return np.full(self.num_layers - 1, self.interlayer_heat_transfer_coefficient)
+
+    def _interlayer_heat_transfer_coefficients(
+        self,
+        conductivity=0.6,
+        density=1000,
+        heat_capacity=4184,
+        storage_height=15,
+        thermocline=1,
+        seconds_per_hour=3600,
+    ) -> np.ndarray:
+        """
+        Interlayer heat transfer coefficient relative to layer state-of-charge. Requires "nominal" layer height, which is assumed to be storage_height distributed equally across layers. Returns array of length num_layers - 1, where each entry corresponds to the coefficient between layer l and l+1.
+
+        Args:
+        ----
+        conductivity: Thermal conductivity of the storage medium [W/m/K]
+        density: Density of the storage medium [kg/m3]
+        heat_capacity: Specific heat capacity of the storage medium [J/kg/K]
+        storage_height: Total height of the storage medium [m]
+        seconds_per_hour: Number of seconds in an hour (for unit conversion)
+        """
+        layer_height = storage_height / (
+            self.num_layers + 1
+        )  # assuming equal layer heights
+
+        return (
+            conductivity
+            / (density * heat_capacity * layer_height)
+            * seconds_per_hour
+            * (self.layer_temperatures[:-1] - self.layer_temperatures[1:])
+            / (self.layer_temperatures[:-1] - self.bottom_temperature)
+            / thermocline
+        )
+
+    @property
+    def boost_ratios(self) -> xr.DataArray:
+        """
+        Resistive boost ratio α^RH_{l,t} per layer per timestep.
+
+        α^RH_{l,t} = max(0, (T_fwd_t - T_l) / (T_l - T_bottom))
+
+        Returns DataArray with dimensions (snapshot, name, layer).
+        """
+        layer_dim = xr.DataArray(
+            self.layer_temperatures,
+            dims=["layer"],
+            coords={"layer": np.arange(len(self.layer_temperatures))},
+        )
+
+        numerator = self.forward_temperature - layer_dim  # (snapshot, name, layer)
+        denominator = layer_dim - self.bottom_temperature  # (layer,)
+
+        alpha = (numerator / denominator).clip(min=0)
+        return alpha
+
+    @property
+    def e_max_pu(self) -> float:
+        """
+        Normalized storage capacity as fraction of design capacity.
+
+        e_max_pu = (T_top - T_bottom) / (T_design_top - T_design_bottom),
+        clipped to non-negative.
+        """
+        delta_t = self.top_temperature - self.bottom_temperature
+        max_delta_t = self.design_top_temperature - self.design_bottom_temperature
+        return max(0, delta_t / max_delta_t)
+
+    @property
+    def standing_losses(self) -> float:
+        """Tbd."""
+        return np.full(self.num_layers, self.design_standing_losses)
+
+    def to_dataset(self) -> xr.Dataset:
+        """
+        Export all pre-computed parameters as a single xr.Dataset.
+
+        Variables:
+        - layer_temperatures: (layer,)
+        - volume_weights: (layer,)
+        - charging_availability: (snapshot, name, layer)
+        - interlayer_transfer_coefficients: (layer_pair,)
+        - boost_ratios: (snapshot, name, layer)
+        - standing_losses: (layer,)
+        """
+        layer_coord = np.arange(self.num_layers)
+        pair_coord = np.arange(self.num_layers - 1)
+
+        ds = xr.Dataset(
+            {
+                "layer_temperatures": xr.DataArray(
+                    self.layer_temperatures,
+                    dims=["layer"],
+                    coords={"layer": layer_coord},
+                ),
+                "volume_weights": xr.DataArray(
+                    self.volume_weights,
+                    dims=["layer"],
+                    coords={"layer": layer_coord},
+                ),
+                "charging_availability": self.charging_availability,
+                "interlayer_heat_transfer_coefficients": xr.DataArray(
+                    self.interlayer_heat_transfer_coefficients,
+                    dims=["layer_pair"],
+                    coords={"layer_pair": pair_coord},
+                ),
+                "boost_ratios": self.boost_ratios,
+                "standing_losses": xr.DataArray(
+                    self.standing_losses,
+                    dims=["layer"],
+                    coords={"layer": layer_coord},
+                ),
+                "e_max_pu": self.e_max_pu,
+            },
+            attrs={
+                "num_layers": self.num_layers,
+                "top_temperature": self.top_temperature,
+                "bottom_temperature": self.bottom_temperature,
+                "storage_height": 15,  # m, assumed for interlayer coefficient calculation
+                "design_top_temperature": self.design_top_temperature,
+                "design_bottom_temperature": self.design_bottom_temperature,
+                "design_standing_losses": self.design_standing_losses,
+            },
+        )
+        return ds
diff --git a/scripts/build_ptes_operations/ptes_temperature_approximator.py b/scripts/build_ptes_operations/ptes_temperature_approximator.py
deleted file mode 100644
index 61140ef31d..0000000000
--- a/scripts/build_ptes_operations/ptes_temperature_approximator.py
+++ /dev/null
@@ -1,315 +0,0 @@
-# SPDX-FileCopyrightText: Contributors to PyPSA-Eur <https://github.com/pypsa/pypsa-eur>
-#
-# SPDX-License-Identifier: MIT
-
-import logging
-
-import xarray as xr
-
-logger = logging.getLogger(__name__)
-
-
-class PtesTemperatureApproximator:
-    """
-    A unified class to handle pit thermal energy storage (PTES) temperature-related calculations.
-
-    It calculates top temperature profiles, determines when charge or discharge boosting is needed,
-    and approximates storage capacity based on temperature differences.
-
-    Attributes
-    ----------
-    forward_temperature : xr.DataArray
-        The forward temperature profile from the district heating network.
-    return_temperature : xr.DataArray
-        The return temperature profile from the district heating network.
-    top_temperature : float | str
-        Operational temperature specification for top layer in PTES.
-    bottom_temperature : float | str
-        Operational temperature specification for bottom layer in PTES.
-    charge_boosting_required : bool
-        Whether charge boosting is required/allowed.
-    discharge_boosting_required : bool
-        Whether discharge boosting is required/allowed.
-    temperature_dependent_capacity : bool
-        Whether storage capacity varies with temperature. If False, assumes constant capacity.
-    design_top_temperature : float
-        Maximum design temperature for the top layer of PTES, used for capacity normalization.
-    design_bottom_temperature : float
-        Minimum design temperature for the bottom layer of PTES, used for capacity normalization.
-    """
-
-    def __init__(
-        self,
-        forward_temperature: xr.DataArray,
-        return_temperature: xr.DataArray,
-        top_temperature: float | str,
-        bottom_temperature: float | str,
-        charge_boosting_required: bool,
-        discharge_boosting_required: bool,
-        temperature_dependent_capacity: bool,
-        design_top_temperature: float,
-        design_bottom_temperature: float,
-    ):
-        """
-        Initialize PtesTemperatureApproximator.
-
-        Parameters
-        ----------
-        forward_temperature : xr.DataArray
-            The forward temperature profile from the district heating network.
-        return_temperature : xr.DataArray
-            The return temperature profile from the district heating network.
-        top_temperature : float | str
-            Operational temperature of top layer in PTES. Either a float value or 'forward' for dynamic profiles.
-        bottom_temperature : float | str
-            Operational temperature of bottom layer in PTES. Either a float value or 'return' for dynamic profiles.
-        charge_boosting_required : bool
-            Whether charge boosting is required/allowed.
-        discharge_boosting_required : bool
-            Whether discharge boosting is required/allowed.
-        temperature_dependent_capacity : bool
-            Whether storage capacity varies with temperature. If False, assumes constant capacity.
-        design_top_temperature : float
-            Maximum design temperature for the top layer of PTES, used for capacity normalization
-            and clipping dynamic top temperature profiles.
-        design_bottom_temperature : float
-            Minimum design temperature for the bottom layer of PTES, used for capacity normalization.
-        """
-        self.forward_temperature = forward_temperature
-        self.return_temperature = return_temperature
-        self.top_temperature = top_temperature
-        self.bottom_temperature = bottom_temperature
-        self.charge_boosting_required = charge_boosting_required
-        self.discharge_boosting_required = discharge_boosting_required
-        self.temperature_dependent_capacity = temperature_dependent_capacity
-        self.design_top_temperature = design_top_temperature
-        self.design_bottom_temperature = design_bottom_temperature
-
-        if self.charge_boosting_required:
-            raise NotImplementedError(
-                "Charge boosting for PTES is currently not supported but might be retintroduced in the future."
-            )
-
-    @property
-    def top_temperature_profile(self) -> xr.DataArray:
-        """
-        PTES top layer temperature profile.
-
-        Returns either the forward temperature (if top_temperature == 'forward'),
-        clipped to the design_top_temperature, or a constant temperature profile
-        (if top_temperature is a numeric value).
-
-        Returns
-        -------
-        xr.DataArray
-            The resulting top temperature profile for PTES.
-        """
-        if self.top_temperature == "forward":
-            logger.info(
-                f"PTES top temperature profile: Using dynamic forward temperature from district heating network, clipped by design top temperature to {self.design_top_temperature}°C "
-                f"Forward temperature range: {float(self.forward_temperature.min().values):.1f}°C to {float(self.forward_temperature.max().values):.1f}°C)"
-            )
-            return self.forward_temperature.where(
-                self.forward_temperature <= self.design_top_temperature,
-                self.design_top_temperature,
-            )
-        elif isinstance(self.top_temperature, (int, float)):
-            logger.info(
-                f"PTES top temperature profile: Using constant temperature of {self.top_temperature}°C "
-                f"for all {self.forward_temperature.size} snapshots and nodes"
-            )
-            return xr.full_like(self.forward_temperature, self.top_temperature)
-        else:
-            raise ValueError(
-                f"Invalid top_temperature: {self.top_temperature}. "
-                "Must be 'forward' or a numeric value."
-            )
-
-    @property
-    def bottom_temperature_profile(self) -> xr.DataArray:
-        """
-        PTES bottom layer temperature profile.
-
-        Returns either the return temperature (if bottom_temperature == 'return')
-        or a constant temperature profile (if bottom_temperature is a numeric value).
-
-        Returns
-        -------
-        xr.DataArray
-            The resulting bottom temperature profile for PTES.
-        """
-        if self.bottom_temperature == "return":
-            logger.info(
-                f"PTES bottom temperature profile: Using dynamic return temperature from district heating network "
-                f"(shape: {self.return_temperature.shape}, range: {float(self.return_temperature.min().values):.1f}°C to {float(self.return_temperature.max().values):.1f}°C)"
-            )
-            return self.return_temperature
-        elif isinstance(self.bottom_temperature, (int, float)):
-            logger.info(
-                f"PTES bottom temperature profile: Using constant temperature of {self.bottom_temperature}°C "
-                f"for all {self.return_temperature.size} snapshots and nodes"
-            )
-            return xr.full_like(self.return_temperature, self.bottom_temperature)
-        else:
-            raise ValueError(
-                f"Invalid bottom_temperature: {self.bottom_temperature}. "
-                "Must be 'return' or a numeric value."
-            )
-
-    @property
-    def e_max_pu(self) -> xr.DataArray:
-        """
-        Calculate e_max_pu for PTES as design_temperature_delta / actual_temperature_delta.
-
-        Returns
-        -------
-        xr.DataArray
-            Normalized delta T values between 0 and 1, representing the
-            available storage capacity as a fraction of maximum design capacity.
-            If temperature_dependent_capacity is False, returns constant capacity of 1.0.
-        """
-        if self.temperature_dependent_capacity:
-            delta_t = self.top_temperature_profile - self.bottom_temperature_profile
-            # Get max possible delta_t for normalization
-            max_delta_t = self.design_top_temperature - self.design_bottom_temperature
-            normalized_delta_t = delta_t / max_delta_t
-            result = normalized_delta_t.clip(min=0)  # Ensure non-negative values
-            logger.info(
-                f"PTES capacity (e_max_pu): Calculating temperature-dependent capacity. "
-                f"Normalization: max_delta_t={max_delta_t:.2f}K (max_top={self.design_top_temperature:.2f}°C, min_bottom={self.design_bottom_temperature:.2f}°C). "
-                f"Resulting capacity range: {float(result.min().values):.3f} to {float(result.max().values):.3f} p.u."
-            )
-            return result
-        else:
-            logger.info(
-                f"PTES capacity (e_max_pu): Using constant capacity of 1.0 p.u. (temperature-independent) "
-                f"for all {self.forward_temperature.size} snapshots and nodes"
-            )
-            return xr.ones_like(self.forward_temperature)
-
-    @property
-    def boost_per_discharge(self) -> xr.DataArray:
-        """
-        Calculate the additional lift required between the store's
-        current top temperature and the forward temperature with the lift
-        already achieved inside the store.
-
-        Notes
-        -----
-        The total thermal output required to reach the forward temperature is:
-
-            Q_total = Q_discharge + Q_boost
-
-        The total heat transfer is partitioned into:
-
-            Q_discharge = Ṽ·ρ·cₚ·(T_top − T_bottom)
-            Q_boost  = Ṽ·ρ·cₚ·(T_forward − T_top)
-
-        Solving this to constant Ṽ gives α as the ratio of required boost to available store energy:
-
-            α = Q_boost / Q_discharge
-              = (T_forward − T_top) / (T_top − T_bottom)
-
-        This expression quantifies the share of PTES output that is covered
-        by stored energy relative to the additional heating needed to meet
-        the desired forward temperature.
-
-        Returns
-        -------
-        xr.DataArray
-            The resulting fraction of PTES charge that must be further heated.
-        """
-        if self.discharge_boosting_required:
-            result = (
-                (self.forward_temperature - self.top_temperature_profile)
-                / (self.top_temperature_profile - self.bottom_temperature_profile)
-            ).where(self.forward_temperature > self.top_temperature_profile, 0)
-
-            # Count how many snapshots require boosting
-            boosting_needed = (
-                (self.forward_temperature > self.top_temperature_profile).sum().values
-            )
-            total_snapshots = self.forward_temperature.size
-
-            logger.info(
-                f"Discharge boosting (boost_per_discharge): Enabled. "
-                f"Boosting required for {int(boosting_needed)}/{total_snapshots} snapshot-node combinations "
-                f"(ratio range: {float(result.min().values):.3f} to {float(result.max().values):.3f})"
-            )
-            return result
-        else:
-            logger.info(
-                f"Discharge boosting (boost_per_discharge): Not required. "
-                f"Returning boost_per_discharge=0 for all {self.forward_temperature.size} snapshots and nodes"
-            )
-            return xr.zeros_like(self.forward_temperature)
-
-    @property
-    def boost_per_charge(self) -> xr.DataArray:
-        """
-        Calculate the additional boost energy ratio required during charging.
-
-        .. note::
-            Charge boosting is currently not implemented and will raise
-            NotImplementedError if charge_boosting_required is True.
-
-        This calculates how much additional energy is needed to raise the PTES
-        top layer from the forward temperature to the design top temperature,
-        relative to the energy delivered by charging.
-
-        Notes
-        -----
-        To fill the storage from the return temperature all the way up to its
-        design top temperature, the total thermal energy required is split into:
-
-            Q_charge = Ṽ·ρ·cₚ·(T_forward − T_bottom)
-            Q_boost  = Ṽ·ρ·cₚ·(T_top − T_forward)
-
-        - Q_charge is the energy delivered by charging to the forward setpoint.
-        - Q_boost is the extra boost energy needed to reach the design top temperature.
-
-        Defining α as the ratio of boost energy to charge energy:
-
-            α = Q_boost / Q_charge
-              = (T_top − T_forward) / (T_forward − T_return)
-
-        Wherever the forward temperature meets or exceeds the design top temperature,
-        α is set to zero since no further boost is needed.
-
-        Returns
-        -------
-        xr.DataArray
-            The ratio of additional boost energy needed to the energy delivered by charging.
-        """
-        if self.charge_boosting_required:
-            # Get the max top temperature value
-            max_top = (
-                self.top_temperature
-                if isinstance(self.top_temperature, (int, float))
-                else self.top_temperature_profile
-            )
-            result = (
-                (
-                    (max_top - self.forward_temperature)
-                    / (self.forward_temperature - self.return_temperature)
-                )
-                .where(self.forward_temperature < max_top, 0)
-                .clip(max=1)
-            )
-
-            # Count how many snapshots require boosting
-            boosting_needed = (self.forward_temperature < max_top).sum().values
-            total_snapshots = self.forward_temperature.size
-
-            logger.info(
-                f"Charge boosting (boost_per_charge): Enabled. "
-                f"Boosting required for {int(boosting_needed)}/{total_snapshots} snapshot-node combinations "
-                f"(ratio range: {float(result.min().values):.3f} to {float(result.max().values):.3f})"
-            )
-            return result
-        else:
-            logger.info(
-                f"Charge boosting (boost_per_charge): Not required. "
-                f"Returning boost_per_charge=0 for all {self.forward_temperature.size} snapshots and nodes"
-            )
-            return xr.zeros_like(self.forward_temperature)
diff --git a/scripts/build_ptes_operations/run.py b/scripts/build_ptes_operations/run.py
index 347c5e5efe..d16ddb3565 100644
--- a/scripts/build_ptes_operations/run.py
+++ b/scripts/build_ptes_operations/run.py
@@ -4,24 +4,12 @@
 """
 Build PTES (Pit Thermal Energy Storage) operational profiles.
 
-This script calculates temperature and capacity profiles for pit thermal energy
-storage systems integrated with district heating networks. It determines:
+This script pre-computes temperature-dependent parameters for pit thermal
+energy storage systems integrated with district heating networks, using the
+unified ``PtesApproximator`` class.
 
-1. **Top temperature profiles**: The operational top layer temperature, either
-   following the district heating forward temperature (clipped to design limits)
-   or a constant value.
-
-2. **Capacity profiles (e_max_pu)**: Normalized storage capacity based on the
-   temperature difference between top and bottom layers, relative to the design
-   temperature difference. This captures how storage capacity varies when
-   operating temperatures deviate from design conditions.
-
-3. **Discharge boosting profiles**: When the PTES top temperature is below the
-   required forward temperature, additional heating (boosting) is needed during
-   discharge. This profile quantifies the ratio of boost energy to stored energy.
-
-The outputs are used by ``prepare_sector_network.py`` to configure PTES storage
-components, charger/discharger links, and optional resistive boosting infrastructure.
+Outputs are consumed by ``prepare_sector_network.py`` and by the COP / heat
+source utilisation profile builders.
 
 Relevant Settings
 -----------------
@@ -31,31 +19,23 @@
         district_heating:
             ptes:
                 enable: true
-                temperature_dependent_capacity: false # if true, e_max_pu varies with temperature difference (static but scaled if top/bottom are constant)
-                charge_boosting_required: false # currently not supported
-                discharge_resistive_boosting: false # if true, adds resistive boosting link for discharge boosting and disables heat pump boosting
-                top_temperature: 90  # or "forward" for dynamic
-                bottom_temperature: 35  # or "return" for dynamic
-                design_top_temperature: 90 # used to compute design temperature difference for e_max_pu if temperature_dependent_capacity is true
-                design_bottom_temperature: 35 # used to compute design temperature difference for e_max_pu if temperature_dependent_capacity is true
+                top_temperature: 90
+                bottom_temperature: 35
+                design_top_temperature: 90
+                design_bottom_temperature: 35
+                discharge_resistive_boosting: false
+                layered:
+                    num_layers: 1
 
 Inputs
 ------
 - ``resources/<run_name>/central_heating_forward_temperature_profiles_base_s_{clusters}_{planning_horizons}.nc``
-    Forward temperature profiles for district heating networks (°C).
 - ``resources/<run_name>/central_heating_return_temperature_profiles_base_s_{clusters}_{planning_horizons}.nc``
-    Return temperature profiles for district heating networks (°C).
 
 Outputs
 -------
-- ``resources/<run_name>/temp_ptes_top_profiles_base_s_{clusters}_{planning_horizons}.nc``
-    PTES top layer temperature profile (°C), clipped to design_top_temperature.
-- ``resources/<run_name>/ptes_e_max_pu_profiles_base_s_{clusters}_{planning_horizons}.nc``
-    Normalized PTES capacity profiles (0-1 p.u.).
-- ``resources/<run_name>/ptes_boost_per_discharge_profiles_base_s_{clusters}_{planning_horizons}.nc``
-    Discharge boosting ratio profiles. Values represent the ratio of boost energy
-    to discharge energy needed when T_forward > T_top. Only used when resistive
-    boosting is enabled.
+- ``resources/<run_name>/ptes_layered_params_base_s_{clusters}_{planning_horizons}.nc``
+    Pre-computed parameter dataset for the PTES model.
 
 References
 ----------
@@ -69,9 +49,7 @@
 import xarray as xr
 
 from scripts._helpers import set_scenario_config
-from scripts.build_ptes_operations.ptes_temperature_approximator import (
-    PtesTemperatureApproximator,
-)
+from scripts.build_ptes_operations.ptes_approximator import PtesApproximator
 
 logger = logging.getLogger(__name__)
 
@@ -87,43 +65,32 @@
 
     set_scenario_config(snakemake)
 
-    logger.info(
-        "Loading district heating temperature profiles and calculating PTES operational profiles"
-    )
+    logger.info("Building PTES operational profiles")
     logger.info(f"PTES configuration: {snakemake.params}")
 
-    # Discharge boosting profiles are calculated when resistive boosting is enabled.
-    # For heat pump-based boosting, add "ptes" to central heating heat sources instead and disable resistive boosting.
-    discharge_boosting_required: bool = snakemake.params.discharge_resistive_boosting
-
-    ptes_temperature_approximator = PtesTemperatureApproximator(
-        forward_temperature=xr.open_dataarray(
-            snakemake.input.central_heating_forward_temperature_profiles
-        ),
-        return_temperature=xr.open_dataarray(
-            snakemake.input.central_heating_return_temperature_profiles
-        ),
-        top_temperature=snakemake.params.top_temperature,
-        bottom_temperature=snakemake.params.bottom_temperature,
-        charge_boosting_required=snakemake.params.charge_boosting_required,
-        discharge_boosting_required=discharge_boosting_required,
-        temperature_dependent_capacity=snakemake.params.temperature_dependent_capacity,
-        design_bottom_temperature=snakemake.params.design_bottom_temperature,
-        design_top_temperature=snakemake.params.design_top_temperature,
+    forward_temperature = xr.open_dataarray(
+        snakemake.input.central_heating_forward_temperature_profiles
     )
-
-    ptes_temperature_approximator.top_temperature_profile.to_netcdf(
-        snakemake.output.ptes_top_temperature_profiles
+    return_temperature = xr.open_dataarray(
+        snakemake.input.central_heating_return_temperature_profiles
     )
 
-    ptes_temperature_approximator.bottom_temperature_profile.to_netcdf(
-        snakemake.output.ptes_bottom_temperature_profiles
-    )
+    # Resolve top/bottom temperature to constant floats
+    top_temp = float(snakemake.params.top_temperature)
+    bottom_temp = float(snakemake.params.bottom_temperature)
 
-    ptes_temperature_approximator.e_max_pu.to_netcdf(
-        snakemake.output.ptes_e_max_pu_profiles
-    )
+    layered_config = snakemake.params.layered
 
-    ptes_temperature_approximator.boost_per_discharge.to_netcdf(
-        snakemake.output.ptes_boost_per_discharge_profiles
+    approximator = PtesApproximator(
+        forward_temperature=forward_temperature,
+        return_temperature=return_temperature,
+        top_temperature=top_temp,
+        bottom_temperature=bottom_temp,
+        num_layers=layered_config["num_layers"],
+        design_top_temperature=snakemake.params.design_top_temperature,
+        design_bottom_temperature=snakemake.params.design_bottom_temperature,
+        design_standing_losses=0.0,  # placeholder; actual value set from costs data
     )
+
+    # Write single dataset with all pre-computed PTES parameters
+    approximator.to_dataset().to_netcdf(snakemake.output.ptes_operations)
diff --git a/scripts/definitions/heat_source.py b/scripts/definitions/heat_source.py
index 6aba369787..a9ff972cf6 100644
--- a/scripts/definitions/heat_source.py
+++ b/scripts/definitions/heat_source.py
@@ -78,6 +78,11 @@ class HeatSource(Enum):
     AIR = "air"
     GROUND = "ground"
     PTES = "ptes"
+    PTES_LAYER_0 = "ptes layer 0"
+    PTES_LAYER_1 = "ptes layer 1"
+    PTES_LAYER_2 = "ptes layer 2"
+    PTES_LAYER_3 = "ptes layer 3"
+    PTES_LAYER_4 = "ptes layer 4"
     # PTX excess heat sources
     ELECTROLYSIS_WASTE = "electrolysis_waste"
     FISCHER_TROPSCH_WASTE = "fischer_tropsch_waste"
@@ -111,7 +116,14 @@ def source_type(self) -> HeatSourceType:
             return HeatSourceType.INEXHAUSTIBLE
         elif self in [HeatSource.GEOTHERMAL, HeatSource.RIVER_WATER]:
             return HeatSourceType.SUPPLY_LIMITED
-        elif self == HeatSource.PTES:
+        elif self in [
+            HeatSource.PTES,
+            HeatSource.PTES_LAYER_0,
+            HeatSource.PTES_LAYER_1,
+            HeatSource.PTES_LAYER_2,
+            HeatSource.PTES_LAYER_3,
+            HeatSource.PTES_LAYER_4,
+        ]:
             return HeatSourceType.STORAGE
         else:
             return HeatSourceType.PROCESS_WASTE
@@ -289,7 +301,10 @@ def requires_heat_pump(self, ptes_discharge_resistive_boosting: bool) -> bool:
         bool
             False for PTES with resistive boosting, True otherwise.
         """
-        if self == HeatSource.PTES and ptes_discharge_resistive_boosting:
+        if (
+            self.source_type == HeatSourceType.STORAGE
+            and ptes_discharge_resistive_boosting
+        ):
             logging.info(
                 "PTES configured with resistive boosting during discharge; "
                 "heat pump not built for PTES."
diff --git a/scripts/definitions/heat_system.py b/scripts/definitions/heat_system.py
index b7ad87994b..2f22d2aa13 100644
--- a/scripts/definitions/heat_system.py
+++ b/scripts/definitions/heat_system.py
@@ -231,15 +231,18 @@ def heat_pump_costs_name(self, heat_source: HeatSource | str) -> str:
         # Check if this is an excess-heat-sourced heat pump
         if heat_source in [
             HeatSource.PTES,
+            HeatSource.PTES_LAYER_0,
+            HeatSource.PTES_LAYER_1,
+            HeatSource.PTES_LAYER_2,
             HeatSource.GEOTHERMAL,
             HeatSource.SEA_WATER,
             HeatSource.RIVER_WATER,
-            HeatSource.ELECTROLYSIS_waste,
-            HeatSource.FISCHER_TROPSCH_waste,
-            HeatSource.SABATIER_waste,
-            HeatSource.HABER_BOSCH_waste,
-            HeatSource.METHANOLISATION_waste,
-            HeatSource.FUEL_CELL_waste,
+            HeatSource.ELECTROLYSIS_WASTE,
+            HeatSource.FISCHER_TROPSCH_WASTE,
+            HeatSource.SABATIER_WASTE,
+            HeatSource.HABER_BOSCH_WASTE,
+            HeatSource.METHANOLISATION_WASTE,
+            HeatSource.FUEL_CELL_WASTE,
         ]:
             return f"{self.central_or_decentral} excess-heat-sourced heat pump"
         else:
diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py
index 1302f10fb3..ed80d9cd63 100755
--- a/scripts/prepare_sector_network.py
+++ b/scripts/prepare_sector_network.py
@@ -2755,8 +2755,7 @@ def add_heat(
     heat_source_direct_utilisation_profile_file: str,
     heat_source_preheater_utilisation_profile_file: str,
     hourly_heat_demand_total_file: str,
-    ptes_e_max_pu_file: str,
-    ptes_boost_per_discharge_profile_file: str,
+    ptes_operations_file: str,
     ates_e_nom_max: str,
     ates_capex_as_fraction_of_geothermal_heat_source: float,
     ates_recovery_factor: float,
@@ -3098,84 +3097,138 @@ def add_heat(
         ):
             n.add("Carrier", f"{heat_system} water pits")
 
-            n.add(
-                "Bus",
-                nodes + f" {heat_system} water pits",
-                location=nodes,
-                carrier=f"{heat_system} water pits",
-                unit="MWh_th",
-            )
+            ptes_ds = xr.open_dataset(ptes_operations_file)
+            num_layers = int(ptes_ds.attrs["num_layers"])
 
             energy_to_power_ratio_water_pit = costs.at[
                 "central water pit storage", "energy to power ratio"
             ]
 
-            n.add(
-                "Link",
-                nodes,
-                suffix=f" {heat_system} water pits charger",
-                bus0=nodes + f" {heat_system} heat",
-                bus1=nodes + f" {heat_system} water pits",
-                efficiency=costs.at[
-                    "central water pit charger",
-                    "efficiency",
-                ],
-                carrier=f"{heat_system} water pits charger",
-                p_nom_extendable=True,
-                lifetime=costs.at["central water pit storage", "lifetime"],
-                marginal_cost=costs.at["central water pit charger", "marginal_cost"],
-            )
+            e_max_pu = float(ptes_ds["e_max_pu"])
 
-            n.add(
-                "Bus",
-                HeatSource.PTES.resource_bus(nodes, heat_system),
-                location=nodes,
-                carrier=f"{heat_system} ptes heat",
-                unit="MWh_th",
-            )
+            for layer in range(num_layers):
+                layer_suffix = f" layer {layer}" if num_layers > 1 else ""
+                n.add(
+                    "Bus",
+                    nodes + f" {heat_system} water pits{layer_suffix}",
+                    location=nodes,
+                    carrier=f"{heat_system} water pits",
+                    unit="MWh_th",
+                )
 
-            n.add(
-                "Link",
-                nodes,
-                suffix=f" {heat_system} water pits discharger",
-                bus0=nodes + f" {heat_system} water pits",
-                bus1=HeatSource.PTES.resource_bus(nodes, heat_system),
-                carrier=f"{heat_system} water pits discharger",
-                efficiency=costs.at[
-                    "central water pit discharger",
-                    "efficiency",
-                ],
-                p_nom_extendable=True,
-                lifetime=costs.at["central water pit storage", "lifetime"],
-            )
-            n.links.loc[
-                nodes + f" {heat_system} water pits charger",
-                "energy to power ratio",
-            ] = energy_to_power_ratio_water_pit
+                n.add(
+                    "Link",
+                    nodes,
+                    suffix=f" {heat_system} water pits charger{layer_suffix}",
+                    bus0=nodes + f" {heat_system} heat",
+                    bus1=nodes + f" {heat_system} water pits{layer_suffix}",
+                    efficiency=costs.at[
+                        "central water pit charger",
+                        "efficiency",
+                    ],
+                    carrier=f"{heat_system} water pits charger",
+                    p_nom_extendable=True,
+                    lifetime=costs.at["central water pit storage", "lifetime"],
+                    marginal_cost=costs.at[
+                        "central water pit charger", "marginal_cost"
+                    ],
+                    p_max_pu=ptes_ds["charging_availability"]
+                    .sel(layer=layer)
+                    .to_pandas(),
+                )
 
-            if options["district_heating"]["ptes"]["temperature_dependent_capacity"]:
-                # Load pre-calculated e_max_pu profiles
-                e_max_pu_data = xr.open_dataarray(ptes_e_max_pu_file)
-                e_max_pu = (
-                    e_max_pu_data.sel(name=nodes).to_pandas().reindex(index=n.snapshots)
+                heat_source = HeatSource(f"ptes{layer_suffix}")
+                n.add(
+                    "Bus",
+                    heat_source.resource_bus(nodes, heat_system),
+                    location=nodes,
+                    carrier=f"{heat_system} ptes heat",
+                    unit="MWh_th",
                 )
-            else:
-                e_max_pu = 1
 
-            n.add(
-                "Store",
-                nodes,
-                suffix=f" {heat_system} water pits",
-                bus=nodes + f" {heat_system} water pits",
-                e_cyclic=True,
-                e_nom_extendable=True,
-                e_max_pu=e_max_pu,
-                carrier=f"{heat_system} water pits",
-                standing_loss=costs.at["central water pit storage", "standing_losses"]
-                / 100,  # convert %/hour into unit/hour
-                capital_cost=costs.at["central water pit storage", "capital_cost"],
-                lifetime=costs.at["central water pit storage", "lifetime"],
-            )
+                n.add(
+                    "Link",
+                    nodes,
+                    suffix=f" {heat_system} water pits discharger{layer_suffix}",
+                    bus0=nodes + f" {heat_system} water pits{layer_suffix}",
+                    bus1=heat_source.resource_bus(nodes, heat_system),
+                    carrier=f"{heat_system} water pits discharger",
+                    efficiency=costs.at[
+                        "central water pit discharger",
+                        "efficiency",
+                    ],
+                    marginal_cost=costs.at[
+                        "central water pit discharger", "marginal_cost"
+                    ],
+                    p_nom_extendable=True,
+                    lifetime=costs.at["central water pit storage", "lifetime"],
+                )
+                n.links.loc[
+                    nodes + f" {heat_system} water pits charger{layer_suffix}",
+                    "energy to power ratio",
+                ] = energy_to_power_ratio_water_pit
+                n.add(
+                    "Store",
+                    nodes,
+                    suffix=f" {heat_system} water pits{layer_suffix}",
+                    bus=nodes + f" {heat_system} water pits{layer_suffix}",
+                    # e_cyclic=True,
+                    e_initial=0,
+                    e_nom_extendable=True,
+                    e_max_pu=e_max_pu,
+                    carrier=f"{heat_system} water pits",
+                    standing_loss=costs.at[
+                        "central water pit storage", "standing_losses"
+                    ]
+                    / 100,
+                    capital_cost=costs.at["central water pit storage", "capital_cost"]
+                    if num_layers == 1
+                    else 0,
+                    lifetime=costs.at["central water pit storage", "lifetime"],
+                )
+
+            # Inter-layer links: bidirectional flow between adjacent layers
+            for pair_idx in range(num_layers - 1):
+                upper_layer = f"layer {pair_idx}"
+                lower_layer = f"layer {pair_idx + 1}"
+
+                n.add(
+                    "Link",
+                    nodes,
+                    suffix=f" {heat_system} water pits inter {upper_layer}-{lower_layer}",
+                    bus0=nodes + f" {heat_system} water pits {upper_layer}",
+                    bus1=nodes + f" {heat_system} water pits {lower_layer}",
+                    carrier=f"{heat_system} water pits interlayer",
+                    p_nom_extendable=True,
+                    p_min_pu=-1,
+                )
+
+            if num_layers > 1:
+                # Main bus for the aggregate Store
+                n.add(
+                    "Bus",
+                    nodes + f" {heat_system} water pits",
+                    location=nodes,
+                    carrier=f"{heat_system} water pits",
+                    unit="MWh_th",
+                )
+
+                n.add(
+                    "Store",
+                    nodes,
+                    suffix=f" {heat_system} water pits",
+                    bus=nodes + f" {heat_system} water pits",
+                    e_cyclic=True,
+                    e_nom_extendable=True,
+                    e_max_pu=e_max_pu,
+                    carrier=f"{heat_system} water pits",
+                    standing_loss=costs.at[
+                        "central water pit storage", "standing_losses"
+                    ]
+                    / 100,
+                    capital_cost=costs.at["central water pit storage", "capital_cost"],
+                    lifetime=costs.at["central water pit storage", "lifetime"],
+                )
 
         if enable_ates and heat_system == HeatSystem.URBAN_CENTRAL:
             n.add("Carrier", f"{heat_system} aquifer thermal energy storage")
@@ -3238,18 +3291,6 @@ def add_heat(
 
             costs_name_heat_pump = heat_system.heat_pump_costs_name(heat_source)
 
-            cop_heat_pump = (
-                cop.sel(
-                    heat_system=heat_system.system_type.value,
-                    heat_source=heat_source.value,
-                    name=nodes,
-                )
-                .transpose("time", "name")
-                .to_pandas()
-                if options["time_dep_hp_cop"]
-                else costs.loc[[costs_name_heat_pump], ["efficiency"]]
-            )
-
             heat_carrier = heat_source.heat_carrier(heat_system)
 
             if heat_source.requires_bus:
@@ -3376,6 +3417,38 @@ def add_heat(
                     "ptes"
                 ]["discharge_resistive_boosting"]
             ):
+                cop_heat_pump = (
+                    cop.sel(
+                        heat_system=heat_system.system_type.value,
+                        heat_source=heat_source.value,
+                        name=nodes,
+                    )
+                    .transpose("time", "name")
+                    .to_pandas()
+                    if options["time_dep_hp_cop"]
+                    else costs.loc[[costs_name_heat_pump], ["efficiency"]]
+                )
+
+                p_min_pu = (
+                    0
+                    if heat_source == HeatSource.PTES
+                    and params.sector["district_heating"]["ptes"]["layered"][
+                        "num_layers"
+                    ]
+                    > 1
+                    else -(cop_heat_pump > 0).astype(float)
+                )
+                capital_cost = (
+                    0
+                    if heat_source
+                    in [
+                        HeatSource.PTES_LAYER_0,
+                        HeatSource.PTES_LAYER_1,
+                        HeatSource.PTES_LAYER_2,
+                    ]
+                    else costs.at[costs_name_heat_pump, "capital_cost"] * overdim_factor
+                )
+
                 n.add(
                     "Link",
                     nodes,
@@ -3384,18 +3457,16 @@ def add_heat(
                     bus1=nodes,
                     bus2=heat_source.get_heat_pump_input_bus(nodes, heat_system),
                     carrier=f"{heat_system} {heat_source} heat pump",
-                    efficiency=1 / cop_heat_pump.clip(lower=0.001).squeeze(),
+                    efficiency=1 / cop_heat_pump.clip(lower=0.001),
                     efficiency2=heat_source.get_heat_pump_efficiency2(cop_heat_pump),
-                    capital_cost=costs.at[costs_name_heat_pump, "capital_cost"]
-                    * overdim_factor,
-                    p_min_pu=-(cop_heat_pump > 0).squeeze().astype(float),
+                    capital_cost=capital_cost,
+                    p_min_pu=p_min_pu,
                     p_max_pu=0,
                     p_nom_extendable=True,
                     lifetime=costs.at[costs_name_heat_pump, "lifetime"],
                 )
 
         if options["resistive_heaters"]:
-            ptes_heat_source = HeatSource.PTES
             key = f"{heat_system.central_or_decentral} resistive heater"
 
             if (
@@ -3404,16 +3475,7 @@ def add_heat(
                 and params.sector["district_heating"]["ptes"][
                     "discharge_resistive_boosting"
                 ]
-                and ptes_heat_source.value
-                in params.heat_sources[heat_system.system_type.value]
             ):
-                ptes_boost_per_discharge_profiles = (
-                    xr.open_dataarray(ptes_boost_per_discharge_profile_file)
-                    .sel(name=nodes)
-                    .to_pandas()
-                    .reindex(index=n.snapshots)
-                )
-
                 n.add(
                     "Bus",
                     nodes,
@@ -3422,24 +3484,34 @@ def add_heat(
                     carrier=f"{heat_system} resistive heat",
                 )
 
-                n.add(
-                    "Link",
-                    nodes,
-                    suffix=f" {heat_system} water pits resistive booster",
-                    bus0=nodes + f" {heat_system} heat",
-                    bus1=nodes + f" {heat_system} resistive heat",
-                    bus2=ptes_heat_source.preheater_input_bus(nodes, heat_system),
-                    # eff = 1 - eff2 (energy conservation)
-                    efficiency=ptes_boost_per_discharge_profiles
-                    / (ptes_boost_per_discharge_profiles + 1),
-                    # Use 1 unit of medium-temperature heat to produce (ptes_boost_per_discharge_profiles + 1) units of district heating
-                    # (similar to HP balance: p_el x COP = p_source + p_el )
-                    efficiency2=1 / (ptes_boost_per_discharge_profiles + 1),
-                    p_nom_extendable=True,
-                    p_min_pu=-(ptes_boost_per_discharge_profiles > 0).astype(float),
-                    carrier=f"{heat_system} water pits resistive booster",
-                    p_max_pu=0,
-                )
+                for layer in range(num_layers):
+                    layer_suffix = f" layer {layer}" if num_layers > 1 else ""
+                    ptes_heat_source = HeatSource(f"ptes{layer_suffix}")
+                    layer_boost_per_discharge_profile = (
+                        ptes_ds["boost_ratios"]
+                        .sel(name=nodes, layer=layer)
+                        .to_pandas()
+                        .reindex(index=n.snapshots)
+                    )
+
+                    n.add(
+                        "Link",
+                        nodes,
+                        suffix=f" {heat_system} water pits resistive booster{layer_suffix}",
+                        bus0=nodes + f" {heat_system} heat",
+                        bus1=nodes + f" {heat_system} resistive heat",
+                        bus2=ptes_heat_source.preheater_input_bus(nodes, heat_system),
+                        # eff = 1 - eff2 (energy conservation)
+                        efficiency=layer_boost_per_discharge_profile
+                        / (layer_boost_per_discharge_profile + 1),
+                        # Use 1 unit of medium-temperature heat to produce (ptes_boost_per_discharge_profiles + 1) units of district heating
+                        # (similar to HP balance: p_el x COP = p_source + p_el )
+                        efficiency2=1 / (layer_boost_per_discharge_profile + 1),
+                        p_nom_extendable=True,
+                        p_min_pu=-(layer_boost_per_discharge_profile > 0).astype(float),
+                        carrier=f"{heat_system} water pits resistive booster",
+                        p_max_pu=0,
+                    )
 
                 n.add(
                     "Link",
@@ -6391,7 +6463,7 @@ def add_import_options(
             heat_source_direct_utilisation_profile_file=snakemake.input.heat_source_direct_utilisation_profiles,
             heat_source_preheater_utilisation_profile_file=snakemake.input.heat_source_preheater_utilisation_profiles,
             hourly_heat_demand_total_file=snakemake.input.hourly_heat_demand_total,
-            ptes_e_max_pu_file=snakemake.input.ptes_e_max_pu_profiles,
+            ptes_operations_file=getattr(snakemake.input, "ptes_operations", ""),
             ates_e_nom_max=snakemake.input.ates_potentials,
             ates_capex_as_fraction_of_geothermal_heat_source=snakemake.params.sector[
                 "district_heating"
@@ -6403,7 +6475,6 @@ def add_import_options(
                 "recovery_factor"
             ],
             enable_ates=snakemake.params.sector["district_heating"]["ates"]["enable"],
-            ptes_boost_per_discharge_profile_file=snakemake.input.ptes_boost_per_discharge_profiles,
             district_heat_share_file=snakemake.input.district_heat_share,
             solar_thermal_total_file=snakemake.input.solar_thermal_total,
             retro_cost_file=snakemake.input.retro_cost,
diff --git a/scripts/solve_network.py b/scripts/solve_network.py
index 9bfc69f543..e8bbfcb533 100644
--- a/scripts/solve_network.py
+++ b/scripts/solve_network.py
@@ -864,10 +864,12 @@ def add_TES_energy_to_power_ratio_constraints(n: pypsa.Network) -> None:
     """
     indices_charger_p_nom_extendable = n.links.index[
         n.links.index.str.contains("water tanks charger|water pits charger")
+        & ~n.links.index.str.contains("layer")
         & n.links.p_nom_extendable
     ]
     indices_stores_e_nom_extendable = n.stores.index[
         n.stores.index.str.contains("water tanks|water pits")
+        & ~n.stores.index.str.contains("layer")
         & n.stores.e_nom_extendable
     ]
 
@@ -930,12 +932,14 @@ def add_TES_charger_ratio_constraints(n: pypsa.Network) -> None:
         n.links.index.str.contains(
             "water tanks charger|water pits charger|aquifer thermal energy storage charger"
         )
+        & ~n.links.index.str.contains("layer")
         & n.links.p_nom_extendable
     ]
     indices_discharger_p_nom_extendable = n.links.index[
         n.links.index.str.contains(
             "water tanks discharger|water pits discharger|aquifer thermal energy storage discharger"
         )
+        & ~n.links.index.str.contains("layer")
         & n.links.p_nom_extendable
     ]
 
@@ -970,6 +974,225 @@ def add_TES_charger_ratio_constraints(n: pypsa.Network) -> None:
     n.model.add_constraints(lhs == 0, name="TES_charger_ratio")
 
 
+def _has_layered_ptes(n: pypsa.Network) -> bool:
+    """Return True if the network contains layered PTES stores."""
+    return n.stores.index.str.contains("water pits layer").any()
+
+
+def add_layered_ptes_volume_capacity_constraint(
+    n: pypsa.Network, ptes_ds: xr.Dataset
+) -> None:
+    """
+    Enforce volume-weighted layer SOC ≤ aggregate e_nom.
+
+    For each node with layered PTES:
+        Σ_l  W_l · e_{l,t}  ≤  ē   ∀ t
+
+    where W_l is the volume weight of layer l and ē is the aggregate
+    store's optimal capacity (e_nom variable).
+    """
+    layer_stores = n.stores.index[
+        n.stores.index.str.contains("water pits layer") & n.stores.e_nom_extendable
+    ]
+    agg_stores = n.stores.index[
+        n.stores.index.str.contains("water pits")
+        & ~n.stores.index.str.contains("layer")
+        & n.stores.e_nom_extendable
+    ]
+
+    if layer_stores.empty or agg_stores.empty:
+        logger.warning(
+            "No valid layered or aggregate PTES stores found. "
+            "Not enforcing volume capacity constraints."
+        )
+        return
+
+    constraints = []
+    for agg in agg_stores:
+        # Find corresponding layer stores (same prefix before "water pits")
+        prefix = agg.split("water pits")[0] + "water pits layer"
+        layers = layer_stores[layer_stores.str.startswith(prefix)]
+        if layers.empty:
+            continue
+
+        weighted_sum = None
+        for layer_name in layers:
+            layer_idx = int(layer_name.split("layer")[-1].strip())
+            w = float(ptes_ds["volume_weights"].sel(layer=layer_idx).item())
+            layer_e = n.model["Store-e"].loc[:, layer_name]
+            term = w * layer_e
+            weighted_sum = term if weighted_sum is None else weighted_sum + term
+
+        agg_e_nom = n.model["Store-e_nom"].loc[agg]
+        # Σ W_l · e_{l,t} - ē ≤ 0
+        constraints.append(weighted_sum - agg_e_nom)
+
+    if not constraints:
+        return
+
+    merged = linopy.expressions.merge(
+        constraints, dim="Store-ext" if PYPSA_V1 else "name"
+    )
+    n.model.add_constraints(merged <= 0, name="layered_ptes_volume_capacity")
+
+
+def add_layered_ptes_heat_pump_capacity_constraint(
+    n: pypsa.Network,
+) -> None:
+    """
+    Tbd
+    """
+    layer_heat_pumps = n.links.index[
+        n.links.index.str.contains("ptes layer") & n.links.p_nom_extendable
+    ]
+    agg_heat_pumps = n.links.index[
+        n.links.index.str.contains("ptes")
+        & ~n.links.index.str.contains("layer")
+        & n.links.p_nom_extendable
+    ]
+
+    if layer_heat_pumps.empty or agg_heat_pumps.empty:
+        logger.warning(
+            "No valid layered or aggregate PTES heat pumps found. "
+            "Not enforcing volume capacity constraints."
+        )
+        return
+
+    constraints = []
+    for agg in agg_heat_pumps:
+        # Find corresponding layer stores (same prefix before "ptes")
+        prefix = agg.split("ptes")[0] + "ptes layer"
+        layers = layer_heat_pumps[layer_heat_pumps.str.startswith(prefix)]
+        if layers.empty:
+            continue
+
+        layer_p_sum = None
+        for layer_name in layers:
+            layer_p = n.model["Link-p"].loc[:, layer_name]
+            layer_p_sum = layer_p if layer_p_sum is None else layer_p_sum + layer_p
+
+        agg_p_nom = n.model["Link-p_nom"].loc[agg]
+        # Σ W_l · e_{l,t} - ē ≤ 0
+        constraints.append(layer_p_sum - agg_p_nom)
+
+    if not constraints:
+        return
+
+    merged = linopy.expressions.merge(
+        constraints, dim="Link-ext" if PYPSA_V1 else "name"
+    )
+    n.model.add_constraints(merged <= 0, name="layered_ptes_heat_pump_capacity")
+
+
+def add_layered_ptes_interlayer_flow_constraint(
+    n: pypsa.Network, ptes_ds: xr.Dataset
+) -> None:
+    """
+    Enforce interlayer heat diffusion proportional to upper-layer SOC.
+
+    For each interlayer link:
+        p_{l→l+1, t}  =  -Ψ_{l,l+1} · e_{l, t}   ∀ t
+
+    where Ψ is the interlayer transfer coefficient and e_{l,t} is the
+    state-of-energy of the upper layer store.
+    """
+    interlayer_links = n.links.index[
+        n.links.carrier.str.contains("water pits interlayer")
+    ]
+
+    if interlayer_links.empty:
+        return
+
+    constraints = []
+    for link_name in interlayer_links:
+        # Extract pair index from link name (e.g., "... inter layer 0-layer 1" → 0)
+        pair_str = link_name.split("inter ")[-1]
+        pair_idx = int(pair_str.split("layer ")[1].split("-")[0])
+        heat_transfer_coef = float(
+            ptes_ds["interlayer_heat_transfer_coefficients"]
+            .sel(layer_pair=pair_idx)
+            .item()
+        )
+
+        upper_bus = n.links.at[link_name, "bus0"]
+        # Find the store on the upper bus
+        upper_store = n.stores.index[n.stores.bus == upper_bus]
+        if upper_store.empty:
+            continue
+        upper_store = upper_store[0]
+
+        link_p = n.model["Link-p"].loc[:, link_name]
+        store_e = n.model["Store-e"].loc[:, upper_store]
+        # p = Ψ · e
+        constraints.append(link_p - heat_transfer_coef * store_e)
+
+    if not constraints:
+        return
+
+    merged = linopy.expressions.merge(
+        constraints, dim="Link, Store" if PYPSA_V1 else "name"
+    )
+    n.model.add_constraints(merged == 0, name="layered_ptes_interlayer_flow")
+
+
+def add_layered_ptes_aggregate_throughput_constraint(
+    n: pypsa.Network, ptes_ds: xr.Dataset
+) -> None:
+    """
+    Limit total volume-weighted charging + discharging power by aggregate capacity.
+
+    For each node with layered PTES:
+        R · Σ_l W_l · (p_ch_{l,t} + p_dis_{l,t})  ≤  ē_nom   ∀ t
+
+    where R is the energy-to-power ratio, W_l the volume weight, p_ch/p_dis
+    the layer charger/discharger operational power, and ē_nom the aggregate
+    store capacity.
+    """
+    layer_chargers = n.links.index[
+        n.links.index.str.contains("water pits charger")
+        & n.links.index.str.contains("layer")
+    ]
+    agg_stores = n.stores.index[
+        n.stores.index.str.contains("water pits")
+        & ~n.stores.index.str.contains("layer")
+        & n.stores.e_nom_extendable
+    ]
+
+    if layer_chargers.empty or agg_stores.empty:
+        return
+
+    constraints = []
+    for agg in agg_stores:
+        prefix = agg.split("water pits")[0] + "water pits"
+
+        chargers = layer_chargers[layer_chargers.str.startswith(prefix + " charger")]
+        dischargers = pd.Index(
+            [c.replace(" charger ", " discharger ") for c in chargers]
+        ).intersection(n.links.index)
+        if chargers.empty or dischargers.empty:
+            continue
+
+        R = n.links.at[chargers[0], "energy to power ratio"]
+
+        weighted_sum = None
+        for ch, dis in zip(chargers, dischargers):
+            layer_idx = int(ch.rsplit("layer ", 1)[-1])
+            w = float(ptes_ds["volume_weights"].sel(layer=layer_idx).item())
+            term = w * (n.model["Link-p"].loc[:, ch] + n.model["Link-p"].loc[:, dis])
+            weighted_sum = term if weighted_sum is None else weighted_sum + term
+
+        agg_e_nom = n.model["Store-e_nom"].loc[agg]
+        constraints.append(R * weighted_sum - agg_e_nom)
+
+    if not constraints:
+        return
+
+    merged = linopy.expressions.merge(
+        constraints, dim="Store-ext" if PYPSA_V1 else "name"
+    )
+    n.model.add_constraints(merged <= 0, name="layered_ptes_aggregate_throughput")
+
+
 def add_battery_constraints(n):
     """
     Add constraint ensuring that charger = discharger, i.e.
@@ -1216,6 +1439,13 @@ def extra_functionality(
             add_TES_energy_to_power_ratio_constraints(n)
             add_TES_charger_ratio_constraints(n)
 
+    if _has_layered_ptes(n):
+        ptes_ds = xr.open_dataset(snakemake.input.ptes_operations)
+        add_layered_ptes_volume_capacity_constraint(n, ptes_ds)
+        add_layered_ptes_interlayer_flow_constraint(n, ptes_ds)
+        # add_layered_ptes_aggregate_throughput_constraint(n, ptes_ds)
+        ptes_ds.close()
+
     add_battery_constraints(n)
     add_lossy_bidirectional_link_constraints(n)
     add_pipe_retrofit_constraint(n)