Couple TES temperatures to district heating supply temperatures

config/config.default.yaml

@@ -471,6 +471,10 @@ sector:
       upper_threshold_ambient_temperature: 10
       rolling_window_ambient_temperature: 72
       relative_annual_temperature_reduction: 0.01
+    ptes:
+      dynamic_capacity: true
+      max_top_temperature: 90 # °C
+      min_bottom_temperature: # °C
     heat_source_cooling: 6 #K
     heat_pump_cop_approximation:
       refrigerant: ammonia

doc/configtables/sector.csv

@@ -17,6 +17,10 @@ district_heating,--,,`prepare_sector_network.py <https://github.com/PyPSA/pypsa-
 -- -- upper_threshold_ambient_temperature,C,float, Assume `min_forward_temperature` if ambient temperature is above this threshold
 -- -- rolling_window_ambient_temperature, h, int, Rolling window size for averaging ambient temperature when approximating supply temperature
 -- -- relative_annual_temperature_reduction,, float, Relative annual reduction of district heating forward and return temperature - defaults to 0.01 (1%)
+-- ptes,,,
+-- -- dynamic_capacity,--,"{true, false}",Add option for dynamic temperature-dependent capacity of pit storage in district heating
+-- -- max_top_temperature,C,float,The maximum top temperature of the pit storage according to DEA technology catalogue (2018)
+-- -- min_bottom_temperature,C,float,The minimum bottom temperature of the pit storage according to DEA technology catalogue (2018)
 -- heat_source_cooling,K,float,Cooling of heat source for heat pumps
 -- heat_pump_cop_approximation,,,
 -- -- refrigerant,--,"{ammonia, isobutane}",Heat pump refrigerant assumed for COP approximation

doc/release_notes.rst

@@ -11,6 +11,10 @@ Release Notes
 Upcoming Release
 ================
 
+* Added option to use dynamic capacity for pit storage using the `e_max_pu` attribute of the store component, which is calculated in the new rule `build_tes_capacity_profiles` and added to the network in `prepare_sector_network`.
+  The dynamic capacity linearly depends on the temperature spread within the storage, which is coupled to the forward and return flow temperatures of the district heating network.
+  The feature can be turned off by setting ``sector: district_heating: ptes: dynamic_capacity: false`` in the configuration file.
+
 * Added rule :mod:`build_co2_sequestration_potentials`, which processes the raw data from `CO2Stop <https://setis.ec.europa.eu/european-co2-storage-
 database_en>`_. Integrated from separate repository (https://github.com/ericzhou571/Co2Storage).
 

rules/build_sector.smk

@@ -380,6 +380,45 @@ rule build_direct_heat_source_utilisation_profiles:
         "../scripts/build_direct_heat_source_utilisation_profiles.py"
 
 
+rule build_tes_capacity_profiles:
+    params:
+        max_top_temperature=config_provider(
+            "sector",
+            "district_heating",
+            "ptes",
+            "max_top_temperature",
+        ),
+        min_bottom_temperature=config_provider(
+            "sector",
+            "district_heating",
+            "ptes",
+            "min_bottom_temperature",
+        ),
+        snapshots=config_provider("snapshots"),
+    input:
+        central_heating_forward_temperature_profiles=resources(
+            "central_heating_forward_temperature_profiles_base_s_{clusters}_{planning_horizons}.nc"
+        ),
+        central_heating_return_temperature_profiles=resources(
+            "central_heating_return_temperature_profiles_base_s_{clusters}_{planning_horizons}.nc"
+        ),
+        regions_onshore=resources("regions_onshore_base_s_{clusters}.geojson"),
+    output:
+        ptes_e_max_pu_profiles=resources(
+            "ptes_e_max_pu_profiles_base_s_{clusters}_{planning_horizons}.nc"
+        ),
+    resources:
+        mem_mb=2000,
+    log:
+        logs("build_tes_capacity_profiles_s_{clusters}_{planning_horizons}.log"),
+    benchmark:
+        benchmarks("build_tes_capacity_profiles/s_{clusters}_{planning_horizons}")
+    conda:
+        "../envs/environment.yaml"
+    script:
+        "../scripts/build_tes_capacity/run.py"
+
+
 def solar_thermal_cutout(wildcards):
     c = config_provider("solar_thermal", "cutout")(wildcards)
     if c == "default":
@@ -1209,6 +1248,15 @@ 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=(
+            resources(
+                "ptes_e_max_pu_profiles_base_s_{clusters}_{planning_horizons}.nc"
+            )
+            if config_provider(
+                "sector", "district_heating", "ptes", "dynamic_capacity"
+            )
+            else []
+        ),
         solar_thermal_total=lambda w: (
             resources("solar_thermal_total_base_s_{clusters}.nc")
             if config_provider("sector", "solar_thermal")(w)

rules/solve_myopic.smk

@@ -74,6 +74,10 @@ rule add_brownfield:
         drop_leap_day=config_provider("enable", "drop_leap_day"),
         carriers=config_provider("electricity", "renewable_carriers"),
         heat_pump_sources=config_provider("sector", "heat_pump_sources"),
+        tes=config_provider("sector", "tes"),
+        dynamic_ptes_capacity=config_provider(
+            "sector", "district_heating", "ptes", "dynamic_capacity"
+        ),
     input:
         unpack(input_profile_tech_brownfield),
         simplify_busmap=resources("busmap_base_s.csv"),

scripts/add_brownfield.py

@@ -292,6 +292,38 @@ def update_heat_pump_efficiency(n: pypsa.Network, n_p: pypsa.Network, year: int)
     )
 
 
+def update_dynamic_ptes_capacity(n, n_p, year):
+    """
+    Updates dynamic pit storage capacity based on district heating temperature changes.
+
+    Parameters
+    ----------
+    n : pypsa.Network
+        Original network.
+    n_p : pypsa.Network
+        Network with updated parameters.
+    year : int
+        Target year for capacity update.
+
+    Returns
+    -------
+    None
+        Updates capacity in-place.
+    """
+    # pit storages in previous iteration
+    dynamic_ptes_idx_previous_iteration = n_p.stores.index[
+        n_p.stores.index.str.contains("water pits")
+    ]
+    # construct names of same-technology dynamic pit storage in the current iteration
+    corresponding_idx_this_iteration = dynamic_ptes_idx_previous_iteration.str[
+        :-4
+    ] + str(year)
+    # update pit storage capacity in previous iteration in-place to capacity in this iteration
+    n_p.stores_t.e_max_pu[dynamic_ptes_idx_previous_iteration] = n.stores_t.e_max_pu[
+        corresponding_idx_this_iteration
+    ].values
+
+
 if __name__ == "__main__":
     if "snakemake" not in globals():
         from _helpers import mock_snakemake
@@ -324,6 +356,9 @@ def update_heat_pump_efficiency(n: pypsa.Network, n_p: pypsa.Network, year: int)
 
     update_heat_pump_efficiency(n, n_p, year)
 
+    if snakemake.params.tes and snakemake.params.dynamic_ptes_capacity:
+        update_dynamic_ptes_capacity(n, n_p, year)
+
     add_brownfield(
         n,
         n_p,

scripts/build_tes_capacity/run.py

@@ -0,0 +1,88 @@
+# SPDX-FileCopyrightText: Contributors to PyPSA-Eur <https://github.com/pypsa/pypsa-eur>
+#
+# SPDX-License-Identifier: MIT
+"""
+Calculate dynamic pit thermal energy storage (PTES) capacity profiles based on
+district heating forward and return flow temperatures. The linear relation between
+temperature difference and capacity is taken from Sorknaes (2018).
+
+The capacity of thermal energy storage systems varies with the temperature difference
+between the forward and return flows in district heating networks assuming a direct
+integration of the storage. This script calculates normalized capacity factors (e_max_pu)
+for PTES systems based on these temperature differences.
+
+Source
+------
+Sorknæs, P. 2018. "Simulation method for a pit seasonal thermal energy storage system with a heat pump in a district heating system", Energy, Volume 152, https://doi.org/10.1016/j.energy.2018.03.152.
+
+Relevant Settings
+-----------------
+
+.. code:: yaml
+    sector:
+        district_heating:
+            dynamic_ptes_capacity: true
+
+Inputs
+------
+- `resources/<run_name>/central_heating_forward_temperature_profiles.nc`: Forward flow temperature profiles
+- `resources/<run_name>/central_heating_return_temperature_profiles.nc`: Return flow temperature profiles
+
+Outputs
+-------
+- `resources/<run_name>/ptes_e_max_pu_profiles.nc`: Normalized PTES capacity profiles
+"""
+
+import logging
+
+import xarray as xr
+from _helpers import set_scenario_config
+from tes_capacity_approximator import TesCapacityApproximator
+
+logger = logging.getLogger(__name__)
+
+if __name__ == "__main__":
+    if "snakemake" not in globals():
+        from _helpers import mock_snakemake
+
+        snakemake = mock_snakemake(
+            "build_tes_parameters",
+            clusters=48,
+            planning_horizons="2050",
+        )
+
+    set_scenario_config(snakemake)
+
+    # Load temperature profiles
+    logger.info("Loading district heating temperature profiles")
+    t_forward = xr.open_dataarray(
+        snakemake.input.central_heating_forward_temperature_profiles
+    )
+    t_return = xr.open_dataarray(
+        snakemake.input.central_heating_return_temperature_profiles
+    )
+
+    # Define operational limits for PTES
+    max_top_temperature = snakemake.params.get("max_top_temperature", 90)  # °C
+    min_bottom_temperature = snakemake.params.get("min_bottom_temperature", 35)  # °C
+
+    logger.info(
+        f"Calculating PTES capacity profiles with max temperature {max_top_temperature}°C"
+    )
+
+    # Create TES capacity approximator
+    tes_capacity_approximator = TesCapacityApproximator(
+        top_temperature=t_forward,
+        bottom_temperature=t_return,
+        max_top_temperature=max_top_temperature,
+        min_bottom_temperature=min_bottom_temperature,
+    )
+
+    # Calculate e_max_pu
+    e_max_pu = tes_capacity_approximator.calculate_e_max_pu()
+
+    # Save output
+    logger.info(
+        f"Saving PTES capacity profiles to {snakemake.output.ptes_e_max_pu_profiles}"
+    )
+    e_max_pu.to_netcdf(snakemake.output.ptes_e_max_pu_profiles)

scripts/build_tes_capacity/tes_capacity_approximator.py

@@ -0,0 +1,104 @@
+# SPDX-FileCopyrightText: Contributors to PyPSA-Eur <https://github.com/pypsa/pypsa-eur>
+#
+# SPDX-License-Identifier: MIT
+
+import xarray as xr
+
+
+class TesCapacityApproximator:
+    """
+    A class to approximate thermal energy storage (TES) capacity based on temperature data.
+
+    This class calculates the normalized temperature difference between the top and bottom
+    layers of a pit thermal energy storage (PTES) system, which determines the available
+    energy storage capacity. The top temperature is clipped to a maximum operational
+    limit of 90°C.
+
+    Attributes
+    ----------
+    top_temperature : xr.DataArray
+        The top temperature data (forward flow temperature).
+    bottom_temperature : xr.DataArray
+        The bottom temperature data (return flow temperature).
+    max_top_temperature : float
+        Maximum operational temperature of top layer in PTES, default 90°C.
+    min_bottom_temperature : float
+        Minimum operational temperature of bottom layer in PTES, default 35°C.
+    """
+
+    def __init__(
+        self,
+        top_temperature: xr.DataArray,
+        bottom_temperature: xr.DataArray,
+        max_top_temperature: float = 90,
+        min_bottom_temperature: float = 35,
+    ):
+        """
+        Initialize TesCapacityApproximator.
+
+        Parameters
+        ----------
+        top_temperature : xr.DataArray
+            The top temperature data (forward flow temperature).
+        bottom_temperature : xr.DataArray
+            The bottom temperature data (return flow temperature).
+        max_top_temperature : float, optional
+            Maximum operational temperature of top layer in PTES, default 90°C.
+        min_bottom_temperature : float, optional
+            Minimum operational temperature of bottom layer in PTES, default 35°C.
+        """
+        self.top_temperature = top_temperature
+        self.bottom_temperature = bottom_temperature
+        self.max_top_temperature = max_top_temperature
+        self.min_bottom_temperature = min_bottom_temperature
+
+    @property
+    def clipped_top_temperature(self) -> xr.DataArray:
+        """
+        Clip top temperature to maximum operational limit.
+
+        Returns
+        -------
+        xr.DataArray
+            Top temperature clipped to maximum operational limit.
+        """
+        return self.top_temperature.where(
+            self.top_temperature <= self.max_top_temperature, self.max_top_temperature
+        )
+
+    @property
+    def e_max_pu(self) -> xr.DataArray:
+        """
+        Calculate the normalized delta T for TES capacity in relation to
+        max and min temperature.
+
+        Returns
+        -------
+        xr.DataArray
+            Normalized delta T values between 0 and 1, representing the
+            available storage capacity as a percentage of maximum capacity.
+        """
+        delta_t = self.clipped_top_temperature - self.bottom_temperature
+        normalized_delta_t = delta_t / (
+            self.max_top_temperature - self.min_bottom_temperature
+        )
+        return normalized_delta_t.clip(min=0)  # Ensure non-negative values
+
+    def calculate_e_max_pu(self, nodes=None) -> xr.DataArray:
+        """
+        Method to compute e_max_pu based on the top and bottom temperatures.
+
+        Parameters
+        ----------
+        nodes : list or pd.Index, optional
+            Subset of nodes to select from the temperature data.
+
+        Returns
+        -------
+        xr.DataArray
+            Computed e_max_pu values, optionally filtered for specific nodes.
+        """
+        result = self.e_max_pu
+        if nodes is not None:
+            result = result.sel(name=nodes)
+        return result

scripts/prepare_sector_network.py

@@ -2760,6 +2760,7 @@ def add_heat(
     cop_profiles_file: str,
     direct_heat_source_utilisation_profile_file: str,
     hourly_heat_demand_total_file: str,
+    ptes_e_max_pu_file: str,
     district_heat_share_file: str,
     solar_thermal_total_file: str,
     retro_cost_file: str,
@@ -3172,18 +3173,30 @@ def add_heat(
                     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
 
+                if options["district_heating"]["ptes"]["dynamic_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)
+                    )
+                else:
+                    e_max_pu = 1
+
                 n.add(
                     "Store",
-                    nodes + f" {heat_system} water pits",
+                    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=1 - np.exp(-1 / 24 / tes_time_constant_days),
                     capital_cost=costs.at["central water pit storage", "capital_cost"],
@@ -6075,6 +6088,7 @@ def add_import_options(
             cop_profiles_file=snakemake.input.cop_profiles,
             direct_heat_source_utilisation_profile_file=snakemake.input.direct_heat_source_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,
             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,