diff --git a/workflow/config/config.common.yaml b/workflow/config/config.common.yaml
index 8fafacbb8..1d8d0de58 100644
--- a/workflow/config/config.common.yaml
+++ b/workflow/config/config.common.yaml
@@ -1,3 +1,7 @@
+pudl_cache:
+  version: v2025.5.0
+  base_url: https://s3.us-west-2.amazonaws.com/pudl.catalyst.coop
+  directory: repo_data/pudl
 pudl_path: s3://pudl.catalyst.coop/v2025.2.0
 
 foresight: 'perfect'
diff --git a/workflow/repo_data/config/config.common.yaml b/workflow/repo_data/config/config.common.yaml
index a5900e7eb..8577cf638 100644
--- a/workflow/repo_data/config/config.common.yaml
+++ b/workflow/repo_data/config/config.common.yaml
@@ -1,5 +1,10 @@
 pudl_path: s3://pudl.catalyst.coop/v2025.5.0
 
+pudl_cache:
+  version: v2025.5.0
+  base_url: https://s3.us-west-2.amazonaws.com/pudl.catalyst.coop
+  directory: repo_data/pudl
+
 # docs :
 renewable:
   EGS:
diff --git a/workflow/rules/build_electricity.smk b/workflow/rules/build_electricity.smk
index e45204d79..682b1bac6 100644
--- a/workflow/rules/build_electricity.smk
+++ b/workflow/rules/build_electricity.smk
@@ -1,6 +1,107 @@
 ################# ----------- Rules to Build Network ---------- #################
 
 from itertools import chain
+from pathlib import Path
+
+
+PUDL_VERSION = config["pudl_cache"]["version"]
+PUDL_BASE_URL = config["pudl_cache"]["base_url"].rstrip("/")
+PUDL_DIRECTORY = Path(config["pudl_cache"]["directory"]) / PUDL_VERSION
+
+PUDL_POWERPLANT_FILES = [
+    "out_eia__monthly_generators.parquet",
+    "out_eia__yearly_generators.parquet",
+    "core_eia860__scd_generators_energy_storage.parquet",
+    "core_eia860__scd_plants.parquet",
+]
+
+PUDL_COST_FILES = [
+    "core_nrelatb__yearly_projected_financial_cases_by_scenario.parquet",
+    "core_nrelatb__yearly_projected_cost_performance.parquet",
+    "core_eiaaeo__yearly_projected_fuel_cost_in_electric_sector_by_type.parquet",
+]
+
+PUDL_DEMAND_FILES = [
+    "censusdp1tract.state_2010census_dp1.parquet",
+    "core_eiaaeo__yearly_projected_generation_in_electric_sector_by_technology.parquet",
+]
+
+PUDL_FUEL_PRICE_FILES = [
+    "out_eia__monthly_generators.parquet",
+    "out_eia__yearly_generators.parquet",
+    "core_eia860__scd_plants.parquet",
+]
+
+PUDL_FILES = sorted(
+    set(
+        PUDL_POWERPLANT_FILES
+        + PUDL_COST_FILES
+        + PUDL_DEMAND_FILES
+        + PUDL_FUEL_PRICE_FILES
+    )
+)
+
+
+def pudl_paths(filenames):
+    return [str(PUDL_DIRECTORY / filename) for filename in filenames]
+
+
+PUDL_POWERPLANT_PATHS = pudl_paths(PUDL_POWERPLANT_FILES)
+PUDL_COST_PATHS = pudl_paths(PUDL_COST_FILES)
+PUDL_FUEL_PRICE_PATHS = pudl_paths(PUDL_FUEL_PRICE_FILES)
+
+
+def demand_pudl_paths(wildcards):
+    profile = config["electricity"]["demand"]["profile"]
+
+    if profile == "ferc":
+        return pudl_paths(PUDL_DEMAND_FILES)
+
+    if profile == "eia":
+        return pudl_paths(
+            [
+                "core_eiaaeo__yearly_projected_generation_in_electric_sector_by_technology.parquet",
+            ]
+        )
+
+    return []
+
+
+rule retrieve_pudl_file:
+    output:
+        str(PUDL_DIRECTORY / "{pudl_file}.parquet"),
+    params:
+        url=lambda wildcards: (
+            f"{PUDL_BASE_URL}/{PUDL_VERSION}/" f"{wildcards.pudl_file}.parquet"
+        ),
+    wildcard_constraints:
+        pudl_file="|".join(Path(filename).stem for filename in PUDL_FILES),
+    log:
+        "logs/retrieve_pudl/{pudl_file}.log",
+    resources:
+        mem_mb=1000,
+        walltime="00:30:00",
+    shell:
+        r"""
+        mkdir -p "$(dirname {output:q})" "$(dirname {log:q})"
+
+        curl -fL \
+          --retry 5 \
+          --retry-all-errors \
+          --connect-timeout 20 \
+          --max-time 1800 \
+          {params.url:q} \
+          --output {output:q}.part \
+          > {log:q} 2>&1
+
+        test -s {output:q}.part
+
+        python -c \
+          "import duckdb; duckdb.sql(\"SELECT * FROM read_parquet('{output}.part') LIMIT 1\").fetchall()" \
+          >> {log:q} 2>&1
+
+        mv {output:q}.part {output:q}
+        """
 
 
 rule build_shapes:
@@ -101,8 +202,9 @@ rule build_bus_regions:
 rule build_cost_data:
     params:
         aeo=config_provider("costs", "aeo"),
-        pudl_path=config_provider("pudl_path"),
+        pudl_path=str(PUDL_DIRECTORY),
     input:
+        pudl=PUDL_COST_PATHS,
         efs_tech_costs="repo_data/costs/EFS_Technology_Data.xlsx",
         efs_icev_costs="repo_data/costs/efs_icev_costs.csv",
         eia_tech_costs="repo_data/costs/eia_tech_costs.csv",
@@ -344,7 +446,6 @@ def demand_raw_data(wildcards):
     elif profile == "ferc":
         return [
             DATA + "pudl/out_ferc714__hourly_estimated_state_demand.parquet",
-            DATA + "pudl/censusdp1tract.sqlite",
         ]
     elif profile == "eulp":
         return [
@@ -422,8 +523,9 @@ rule build_electrical_demand:
         planning_horizons=config["scenario"]["planning_horizons"],
         renewable_weather_years=config["renewable_weather_years"],
         snapshots=config["snapshots"],
-        pudl_path=config_provider("pudl_path"),
+        pudl_path=str(PUDL_DIRECTORY),
     input:
+        pudl=demand_pudl_paths,
         network=RESOURCES + "{interconnect}/elec_base_network.nc",
         demand_files=demand_raw_data,
         demand_scaling_file=demand_scaling_data,
@@ -457,10 +559,12 @@ rule build_service_demand:
         dissagregate_files=demand_dissagregate_data,
         demand_scaling_file=demand_scaling_data,
     output:
-        electricity=RESOURCES + "{interconnect}/demand/{end_use}_electricity.pkl",
-        space_heat=RESOURCES + "{interconnect}/demand/{end_use}_space-heating.pkl",
-        water_heat=RESOURCES + "{interconnect}/demand/{end_use}_water-heating.pkl",
-        cool=RESOURCES + "{interconnect}/demand/{end_use}_cooling.pkl",
+        electricity=RESOURCES + "{interconnect}/demand/sector/{end_use}_electricity.pkl",
+        space_heat=RESOURCES
+        + "{interconnect}/demand/sector/{end_use}_space-heating.pkl",
+        water_heat=RESOURCES
+        + "{interconnect}/demand/sector/{end_use}_water-heating.pkl",
+        cool=RESOURCES + "{interconnect}/demand/sector/{end_use}_cooling.pkl",
     log:
         LOGS + "{interconnect}/demand/{end_use}_build_demand.log",
     benchmark:
@@ -480,15 +584,14 @@ rule build_industry_demand:
         profile_year=pd.to_datetime(config["snapshots"]["start"]).year,
         eia_api=config_provider("api", "eia"),
         snapshots=config_provider("snapshots"),
-        pudl_path=config_provider("pudl_path"),
     input:
         network=RESOURCES + "{interconnect}/elec_base_network.nc",
         demand_files=demand_raw_data,
         dissagregate_files=demand_dissagregate_data,
         demand_scaling_file=demand_scaling_data,
     output:
-        electricity=RESOURCES + "{interconnect}/demand/{end_use}_electricity.pkl",
-        heat=RESOURCES + "{interconnect}/demand/{end_use}_heating.pkl",
+        electricity=RESOURCES + "{interconnect}/demand/sector/{end_use}_electricity.pkl",
+        heat=RESOURCES + "{interconnect}/demand/sector/{end_use}_heating.pkl",
     log:
         LOGS + "{interconnect}/demand/{end_use}_build_demand.log",
     benchmark:
@@ -515,10 +618,10 @@ rule build_transport_road_demand:
         dissagregate_files=demand_dissagregate_data,
         demand_scaling_file=demand_scaling_data,
     output:
-        light_duty=RESOURCES + "{interconnect}/demand/{end_use}_light-duty.pkl",
-        med_duty=RESOURCES + "{interconnect}/demand/{end_use}_med-duty.pkl",
-        heavy_duty=RESOURCES + "{interconnect}/demand/{end_use}_heavy-duty.pkl",
-        bus=RESOURCES + "{interconnect}/demand/{end_use}_bus.pkl",
+        light_duty=RESOURCES + "{interconnect}/demand/sector/{end_use}_light-duty.pkl",
+        med_duty=RESOURCES + "{interconnect}/demand/sector/{end_use}_med-duty.pkl",
+        heavy_duty=RESOURCES + "{interconnect}/demand/sector/{end_use}_heavy-duty.pkl",
+        bus=RESOURCES + "{interconnect}/demand/sector/{end_use}_bus.pkl",
     log:
         LOGS + "{interconnect}/demand/{end_use}_build_demand.log",
     benchmark:
@@ -546,7 +649,7 @@ rule build_transport_other_demand:
         demand_files=demand_raw_data,
         dissagregate_files=demand_dissagregate_data,
     output:
-        RESOURCES + "{interconnect}/demand/{end_use}_{vehicle}.pkl",
+        RESOURCES + "{interconnect}/demand/sector/{end_use}_{vehicle}.pkl",
     log:
         LOGS + "{interconnect}/demand/{end_use}_{vehicle}_build_demand.log",
     benchmark:
@@ -558,53 +661,80 @@ rule build_transport_other_demand:
         "../scripts/build_demand.py"
 
 
-def demand_to_add(wildcards):
+def sector_demand_files(wildcards):
+    """
+    Return compact sector-demand inputs for the clustered network.
+
+    Parameters
+    ----------
+    wildcards : snakemake.io.Wildcards
+        Wildcards supplied by Snakemake. Returned paths retain the existing
+        interconnect placeholder for workflow expansion.
+
+    Returns
+    -------
+    list or itertools.chain
+        An empty list for electricity-only studies. For sector studies,
+        an iterator over residential, commercial, industrial, road-transport,
+        and other transport demand files.
+
+    Notes
+    -----
+    These files are consumed by add_extra_components after spatial
+    clustering. They are not inputs to the initial add_demand rule.
+
+    The existing service-sector configuration determines whether heating
+    demand is split into space and water heating. No new configuration
+    options are introduced.
+    """
+    if config["scenario"]["sector"] in ("E", ""):
+        return []
 
-    if config["scenario"]["sector"] == "E":
-        return RESOURCES + "{interconnect}/demand/power_electricity.csv"
+    services = ["residential", "commercial"]
+    if config["sector"]["service_sector"]["split_space_water_heating"]:
+        fuels = ["electricity", "cooling", "space-heating", "water-heating"]
     else:
-        # service demand
-        services = ["residential", "commercial"]
-        if config["sector"]["service_sector"]["split_space_water_heating"]:
-            fuels = ["electricity", "cooling", "space-heating", "water-heating"]
-        else:
-            fuels = ["electricity", "cooling", "heating"]
-        service_demands = [
-            RESOURCES + "{interconnect}/demand/" + service + "_" + fuel + ".pkl"
-            for service in services
-            for fuel in fuels
-        ]
-        # industrial demand
-        fuels = ["electricity", "heating"]
-        industrial_demands = [
-            RESOURCES + "{interconnect}/demand/industry_" + fuel + ".pkl"
-            for fuel in fuels
-        ]
-        # road transport demands
-        vehicles = ["light-duty", "med-duty", "heavy-duty", "bus"]
-        road_demand = [
-            RESOURCES + "{interconnect}/demand/transport_" + vehicle + ".pkl"
-            for vehicle in vehicles
-        ]
-
-        # other transport demands
-        vehicles = ["boat-shipping", "rail-shipping", "rail-passenger", "air"]
-        non_road_demand = [
-            RESOURCES + "{interconnect}/demand/transport_" + vehicle + ".pkl"
-            for vehicle in vehicles
-        ]
-
-        return chain(service_demands, industrial_demands, road_demand, non_road_demand)
+        fuels = ["electricity", "cooling", "heating"]
+
+    service_demands = [
+        RESOURCES + "{interconnect}/demand/sector/" + service + "_" + fuel + ".pkl"
+        for service in services
+        for fuel in fuels
+    ]
+
+    fuels = ["electricity", "heating"]
+    industrial_demands = [
+        RESOURCES + "{interconnect}/demand/sector/industry_" + fuel + ".pkl"
+        for fuel in fuels
+    ]
+
+    vehicles = ["light-duty", "med-duty", "heavy-duty", "bus"]
+    road_demand = [
+        RESOURCES + "{interconnect}/demand/sector/transport_" + vehicle + ".pkl"
+        for vehicle in vehicles
+    ]
+
+    vehicles = ["boat-shipping", "rail-shipping", "rail-passenger", "air"]
+    non_road_demand = [
+        RESOURCES + "{interconnect}/demand/sector/transport_" + vehicle + ".pkl"
+        for vehicle in vehicles
+    ]
+
+    return chain(
+        service_demands,
+        industrial_demands,
+        road_demand,
+        non_road_demand,
+    )
 
 
 rule add_demand:
     params:
-        sectors=config["scenario"]["sector"],
         planning_horizons=config_provider("scenario", "planning_horizons"),
         snapshots=config_provider("snapshots"),
     input:
         network=RESOURCES + "{interconnect}/elec_base_network.nc",
-        demand=demand_to_add,
+        demand=RESOURCES + "{interconnect}/demand/power_electricity.csv",
     output:
         network=RESOURCES + "{interconnect}/elec_base_network_dem.nc",
     log:
@@ -629,8 +759,9 @@ rule build_fuel_prices:
     params:
         snapshots=config["snapshots"],
         api_eia=config["api"]["eia"],
-        pudl_path=config_provider("pudl_path"),
+        pudl_path=str(PUDL_DIRECTORY),
     input:
+        pudl=PUDL_FUEL_PRICE_PATHS,
         gas_balancing_area=ba_gas_dynamic_fuel_price_files,
     output:
         state_ng_fuel_prices=RESOURCES + "{interconnect}/state_ng_power_prices.csv",
@@ -665,8 +796,10 @@ def dynamic_fuel_price_files(wildcards):
 
 rule build_powerplants:
     params:
-        pudl_path=config_provider("pudl_path"),
+        pudl_path=str(PUDL_DIRECTORY),
+        renewable_weather_year=config_provider("renewable_weather_years"),
     input:
+        pudl=PUDL_POWERPLANT_PATHS,
         wecc_ads="repo_data/WECC_ADS_public",
         eia_ads_generator_mapping="repo_data/WECC_ADS_public/eia_ads_generator_mapping_updated.csv",
         fuel_costs="repo_data/plants/fuelCost22.csv",
@@ -790,6 +923,7 @@ rule simplify_network:
         + "{interconnect}/Geospatial/regions_offshore.geojson",
     output:
         network=RESOURCES + "{interconnect}/elec_s{simpl}.nc",
+        busmap=RESOURCES + "{interconnect}/busmap_s{simpl}.csv",
         regions_onshore=RESOURCES
         + "{interconnect}/Geospatial/regions_onshore_s{simpl}.geojson",
         regions_offshore=RESOURCES
@@ -881,6 +1015,17 @@ rule add_extra_components:
             if hour.isdigit()
         },
         network=RESOURCES + "{interconnect}/elec_s{simpl}_c{clusters}.nc",
+        sector_demand=lambda w: list(sector_demand_files(w)),
+        busmap_s=(
+            RESOURCES + "{interconnect}/busmap_s{simpl}.csv"
+            if config["scenario"]["sector"] not in ("E", "")
+            else []
+        ),
+        busmap_c=(
+            RESOURCES + "{interconnect}/busmap_s{simpl}_{clusters}.csv"
+            if config["scenario"]["sector"] not in ("E", "")
+            else []
+        ),
         tech_costs=lambda wildcards: expand(
             RESOURCES + "costs/costs_{year}.csv",
             year=config["scenario"]["planning_horizons"],
diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk
index 3ddae2eb6..07295d0eb 100644
--- a/workflow/rules/common.smk
+++ b/workflow/rules/common.smk
@@ -7,8 +7,11 @@ from functools import partial, lru_cache
 
 import os, sys, glob
 
-path = workflow.source_path("../scripts/_helpers.py")
-sys.path.insert(0, os.path.dirname(path))
+# Cache local modules before importing them from Snakemake's source cache.
+for source in ("../scripts/_helpers.py", "../scripts/constants.py"):
+    source_dir = os.path.dirname(workflow.source_path(source))
+    if source_dir not in sys.path:
+        sys.path.insert(0, source_dir)
 
 from _helpers import validate_checksum, update_config_from_wildcards
 from constants import HOURS_PER_YEAR
diff --git a/workflow/scripts/_helpers.py b/workflow/scripts/_helpers.py
index 3f33c662d..d5aae598b 100644
--- a/workflow/scripts/_helpers.py
+++ b/workflow/scripts/_helpers.py
@@ -17,6 +17,8 @@
 
 REGION_COLS = ["geometry", "name", "x", "y", "country"]
 
+logger = logging.getLogger(__name__)
+
 
 def configure_logging(snakemake, skip_handlers=False):
     """
@@ -866,3 +868,272 @@ def get_multiindex_snapshots(
             get_snapshots(sns_config).map(lambda x: x.replace(year=year)),
         )
     return pd.MultiIndex.from_arrays([sns.year, sns])
+
+
+def prepare_sector_demand(profiles, allocation, growth, unit_conversion=1.0):
+    """
+    Package zonal sector demand and its original spatial allocation factors.
+
+    Parameters
+    ----------
+    profiles : pd.DataFrame
+        Zonal demand profiles indexed by the timestamps of all required
+        investment periods. Columns identify source zones. Demand growth
+        and unit conversion must not have been applied yet.
+    allocation : pd.DataFrame
+        Allocation table indexed by unique original network bus identifiers.
+        Required columns are ``zone``, identifying the source demand zone,
+        and ``laf``, containing the original nodal load allocation factor.
+    growth : pd.Series
+        Multiplicative demand growth factor for each timestamp. Its index
+        must match ``profiles.index`` exactly, including ordering.
+    unit_conversion : float, default 1.0
+        Conversion applied after nodal allocation and demand growth.
+        Transport demand may require a conversion; other demand categories
+        normally use 1.0.
+
+    Returns
+    -------
+    dict
+        Compact demand representation with ``profiles``, ``allocation``,
+        ``growth``, and ``unit_conversion`` entries. Allocation rows contain
+        only buses passing the existing nodal cutoff and having a matching
+        source-zone profile.
+
+    Raises
+    ------
+    ValueError
+        If profiles are empty, profile and growth indices differ, bus or
+        profile-zone identifiers are duplicated, or numerical inputs contain
+        non-finite values.
+
+    Notes
+    -----
+    Allocation factors below 1e-6 are discarded. Buses whose source zones
+    have no demand profile are skipped with a warning, as can occur near
+    interconnection boundaries. Retained factors are not renormalized.
+
+    No hourly nodal demand matrix is constructed here. Allocation, growth,
+    conversion, and nodal rounding are evaluated when demand is aggregated
+    onto the clustered network.
+    """
+    if profiles.empty or not profiles.index.equals(growth.index):
+        raise ValueError(
+            "Sector profiles and growth factors must have matching, nonempty snapshots",
+        )
+    if not allocation.index.is_unique or not profiles.columns.is_unique:
+        raise ValueError("Sector demand buses and profile zones must be unique")
+    if not np.isfinite(profiles.to_numpy()).all() or not np.isfinite(growth.to_numpy()).all():
+        raise ValueError("Sector demand profiles and growth factors must be finite")
+    if not np.isfinite(allocation.laf.to_numpy()).all() or not np.isfinite(unit_conversion):
+        raise ValueError(
+            "Sector demand allocation factors and unit conversion must be finite",
+        )
+
+    # Apply the original nodal cutoff before any spatial aggregation.
+    allocation = allocation.loc[
+        allocation.laf >= 0.000001,
+        ["zone", "laf"],
+    ].copy()
+
+    missing = allocation.loc[
+        ~allocation.zone.isin(profiles.columns),
+        "zone",
+    ].unique()
+    if len(missing):
+        logger.warning("No demand found for %s", missing.tolist())
+
+    allocation = allocation.loc[allocation.zone.isin(profiles.columns)]
+
+    return {
+        "profiles": profiles,
+        "allocation": allocation,
+        "growth": growth,
+        "unit_conversion": unit_conversion,
+    }
+
+
+def read_busmap(path):
+    """
+    Read a bus mapping without changing textual bus identifiers.
+
+    Parameters
+    ----------
+    path : str or path-like
+        CSV file with exactly two columns: original source bus identifiers
+        followed by destination bus identifiers. Empty destinations indicate
+        buses explicitly removed from the network.
+
+    Returns
+    -------
+    pd.Series
+        Mapping indexed by source bus identifiers. Identifiers are preserved
+        as strings, including leading zeros. Empty destinations become NaN.
+
+    Raises
+    ------
+    ValueError
+        If the file does not contain exactly two columns or source bus
+        identifiers are duplicated.
+
+    Notes
+    -----
+    Both columns are read as ordinary string columns before setting the
+    index. This prevents CSV index inference from converting identifiers
+    such as ``"001"`` into integers.
+
+    Default NA parsing is disabled so identifiers such as ``"NA"`` and
+    ``"nan"`` remain literal strings. Only empty destination fields represent
+    removed buses.
+    """
+    frame = pd.read_csv(path, dtype=str, keep_default_na=False)
+
+    if frame.shape[1] != 2 or not frame.iloc[:, 0].is_unique:
+        raise ValueError(
+            f"Expected source and destination columns with unique source buses in {path}",
+        )
+
+    frame = frame.set_index(frame.columns[0])
+    return frame.iloc[:, 0].replace("", np.nan)
+
+
+def compose_busmaps(first, second):
+    """
+    Compose successive spatial mappings while retaining explicit removals.
+
+    Parameters
+    ----------
+    first : pd.Series
+        Mapping from original buses to intermediate buses, indexed by unique
+        original bus identifiers.
+    second : pd.Series
+        Mapping from intermediate buses to final buses, indexed by unique
+        intermediate bus identifiers.
+
+    Returns
+    -------
+    pd.Series
+        Mapping from original buses to final buses, retaining the index of
+        ``first``. Null destinations propagate through the composition.
+
+    Raises
+    ------
+    ValueError
+        If either mapping has duplicate source identifiers or a non-null
+        destination in ``first`` is absent from the index of ``second``.
+
+    Notes
+    -----
+    A null destination represents an intentional removal. A missing
+    intermediate mapping entry is treated as an error rather than silently
+    interpreted as a removed bus.
+    """
+    if not first.index.is_unique or not second.index.is_unique:
+        raise ValueError("Bus maps must have unique source buses")
+
+    missing = pd.Index(first.dropna().unique()).difference(second.index)
+    if len(missing):
+        raise ValueError(
+            f"Bus map is missing intermediate buses: {missing.tolist()[:10]}",
+        )
+
+    return first.map(second)
+
+
+def aggregate_sector_demand(demand, busmap, block_size=128):
+    """
+    Construct sector demand directly on the surviving clustered buses.
+
+    Parameters
+    ----------
+    demand : dict
+        Compact demand representation produced by ``prepare_sector_demand``.
+        Profiles contain investment-period timestamps, while allocation rows
+        refer to original network buses.
+    busmap : pd.Series
+        Mapping from original buses to final clustered buses. Every retained
+        allocation bus must appear in the mapping index. Null destinations
+        identify buses explicitly removed during network processing.
+    block_size : int, default 128
+        Maximum number of original buses reconstructed simultaneously within
+        each source zone. Must be positive.
+
+    Returns
+    -------
+    pd.DataFrame
+        Demand indexed by ``demand["profiles"].index``, with sorted columns
+        identifying surviving destination buses receiving demand. Values
+        retain the units of the demand category after conversion.
+
+    Raises
+    ------
+    ValueError
+        If block_size is below one or an allocation bus is missing from the
+        bus-map index.
+
+    Notes
+    -----
+    Original nodal values are calculated in this order:
+
+    ``zonal_profile * allocation_factor * growth * unit_conversion``
+
+    Values are rounded to four decimal places at the original bus level
+    before summation onto clustered buses. Allocation factors are not
+    renormalized or aggregated before rounding.
+
+    Demand assigned to explicit null destinations is discarded with a
+    warning. Missing source mapping entries instead raise an error.
+
+    Only the final clustered matrix and one bounded nodal block are
+    materialized. Floating-point summation order can differ from the original
+    network aggregation.
+    """
+    if block_size < 1:
+        raise ValueError("Demand block size must be positive")
+
+    allocation = demand["allocation"]
+
+    missing = allocation.index.difference(busmap.index)
+    if len(missing):
+        raise ValueError(
+            f"Bus map is missing demand buses: {missing.tolist()[:10]}",
+        )
+
+    destinations = busmap.reindex(allocation.index)
+    removed = destinations.isna()
+
+    if removed.any():
+        logger.warning(
+            "Dropping sector demand on %s removed buses",
+            removed.sum(),
+        )
+
+    allocation = allocation.loc[~removed]
+    destinations = destinations.loc[~removed]
+
+    result = pd.DataFrame(
+        0.0,
+        index=demand["profiles"].index,
+        columns=pd.Index(sorted(destinations.unique())),
+    )
+    growth = demand["growth"].to_numpy()[:, None]
+
+    for zone, group in allocation.groupby("zone", sort=False):
+        profile = demand["profiles"][zone].to_numpy(dtype=float)[:, None]
+
+        for start in range(0, len(group), block_size):
+            buses = group.iloc[start : start + block_size]
+
+            # Preserve nodal multiplication and rounding before aggregation.
+            values = profile * buses.laf.to_numpy()[None, :]
+            values *= growth
+            values *= demand["unit_conversion"]
+
+            block = pd.DataFrame(
+                np.round(values, 4).T,
+                index=buses.index,
+            )
+            grouped = block.groupby(destinations.loc[buses.index]).sum().T
+            result.loc[:, grouped.columns] += grouped.to_numpy()
+
+    return result
diff --git a/workflow/scripts/add_demand.py b/workflow/scripts/add_demand.py
index a7c5a4afc..34bb41702 100644
--- a/workflow/scripts/add_demand.py
+++ b/workflow/scripts/add_demand.py
@@ -1,8 +1,15 @@
 """
-Adds demand to the network.
+Attach electricity demand and provide clustered sector-demand attachment.
 
-Depending on study, the load will all be aggregated to a single load
-type, or distributed to different sectors and end use fuels.
+The add_demand workflow rule always attaches the configured electricity
+demand to the original network. Network simplification and clustering
+therefore operate on electricity loads for both electricity-only and
+sector studies.
+
+The attach_sector_demand function is called by add_extra_components after
+spatial clustering. It replaces the initial electricity loads with the
+sector-specific demand profiles, before subsequent demand-dependent
+network preparation.
 """
 
 import logging
@@ -10,7 +17,14 @@
 
 import pandas as pd
 import pypsa
-from _helpers import configure_logging, get_multiindex_snapshots, mock_snakemake
+from _helpers import (
+    aggregate_sector_demand,
+    compose_busmaps,
+    configure_logging,
+    get_multiindex_snapshots,
+    mock_snakemake,
+    read_busmap,
+)
 from constants_sector import (
     TRANSPORT_FUELS,
     SecCarriers,
@@ -22,14 +36,45 @@
 
 def attach_demand(n: pypsa.Network, df: pd.DataFrame, carrier: str, suffix: str):
     """
-    Add demand to network from specified configuration setting.
-
-    Returns network with demand added.
+    Add bus-indexed demand profiles as network Load components.
+
+    Parameters
+    ----------
+    n : pypsa.Network
+        Network receiving the loads. Its snapshots must already be configured.
+    df : pd.DataFrame
+        Demand with timestamps as rows and existing bus identifiers as columns.
+        Rows must correspond positionally to the network snapshots.
+    carrier : str
+        Carrier assigned to all loads created by this call.
+    suffix : str
+        Suffix appended to each bus identifier to construct the Load name.
+        Bus assignments use the unmodified column identifiers.
+
+    Returns
+    -------
+    None
+        Loads are added to the supplied network in place.
+
+    Raises
+    ------
+    AssertionError
+        If the number of demand rows differs from the network snapshot count.
+
+    Notes
+    -----
+    The input DataFrame index is converted to timestamps and then replaced
+    by the network snapshot index in place. This supports investment-period
+    MultiIndices while preserving the existing positional attachment.
+
+    This helper validates row count, not timestamp equality. Callers requiring
+    timestamp validation must perform it before attachment.
     """
     df.index = pd.to_datetime(df.index)
     assert len(df.index) == len(
         n.snapshots,
     ), "Demand time series length does not match network snapshots"
+
     df.index = n.snapshots
     n.madd(
         "Load",
@@ -41,62 +86,134 @@ def attach_demand(n: pypsa.Network, df: pd.DataFrame, carrier: str, suffix: str)
     )
 
 
-if __name__ == "__main__":
-    if "snakemake" not in globals():
-        snakemake = mock_snakemake("add_demand", interconnect="western")
-    configure_logging(snakemake)
-
-    demand_files = snakemake.input.demand
-    n = pypsa.Network(snakemake.input.network)
-
-    sectors = snakemake.params.sectors
-
-    # add snapshots
-    sns_config = snakemake.params.snapshots
-    planning_horizons = snakemake.params.planning_horizons
-
-    n.snapshots = get_multiindex_snapshots(sns_config, planning_horizons)
-    n.set_investment_periods(periods=planning_horizons)
-
+def attach_sector_demand(n, demand_files, busmap_s, busmap_c):
+    """
+    Replace initial electricity loads with clustered sector demand.
+
+    Parameters
+    ----------
+    n : pypsa.Network
+        Spatially clustered network before temporal aggregation, trimming,
+        or demand-response preparation. Existing loads must have carrier
+        ``AC``, or the network must contain no loads.
+    demand_files : str or list[str]
+        Compact sector-demand pickle files produced by build_demand.
+        Filenames follow ``{sector}_{end-use}.pkl``, for example
+        ``residential_electricity.pkl`` or ``transport_light-duty.pkl``.
+    busmap_s : str or path-like
+        CSV mapping original buses to simplified buses, including transformer
+        removal, substation aggregation, and optional simplification clustering.
+    busmap_c : str or path-like
+        CSV mapping simplified buses to final clustered buses. Null
+        destinations represent buses explicitly removed during processing.
+
+    Returns
+    -------
+    None
+        Existing electricity Load components are removed and sector-specific
+        loads are added to the supplied network in place.
+
+    Raises
+    ------
+    ValueError
+        If no demand files are supplied, existing loads have unexpected
+        carriers, mappings are invalid or incomplete, final destinations are
+        absent from the network, filenames are malformed, or timestamps differ.
+    KeyError
+        If a sector or end-use label is unknown, or a compact demand file
+        lacks a required field.
+
+    Notes
+    -----
+    The initial electricity demand is replaced, not added to the electrical
+    components of sector demand. This avoids double counting.
+
+    The AC Carrier itself is retained because network buses and other
+    electrical components continue to use it.
+
+    Original bus allocation factors, demand growth, unit conversion, and
+    nodal rounding are applied before aggregation onto final buses. Load
+    names use ``"{bus} {carrier}"``, matching the spatial aggregation naming
+    convention expected by downstream sector builders.
+
+    This function does not export the network or alter its snapshots.
+    Changes are not transactional: an error while processing a later file
+    can occur after electricity loads have been removed.
+    """
     if isinstance(demand_files, str):
         demand_files = [demand_files]
 
-    if sectors == "E" or sectors == "":  # electricity only
-        assert len(demand_files) == 1
-
-        suffix = ""
-        carrier = "AC"
+    if not demand_files:
+        raise ValueError("No sector demand files supplied")
 
-        df = pd.read_csv(demand_files[0], index_col=0)
-        attach_demand(n, df, carrier, suffix)
-        logger.info("Electricity demand added to network")
+    if not n.loads.carrier.eq("AC").all():
+        raise ValueError(
+            "Expected an electricity-only network before attaching sector demand; "
+            "rebuild add_demand, simplify_network and cluster_network",
+        )
 
-    else:  # sector files
-        for demand_file in demand_files:
-            parsed_name = Path(demand_file).name.split("_")
-            parsed_name[-1] = parsed_name[-1].split(".pkl")[0]
-
-            if len(parsed_name) == 2:
-                sector = parsed_name[0].upper()
-                end_use = parsed_name[1].upper().replace("-", "_")
+    busmap = compose_busmaps(
+        read_busmap(busmap_s),
+        read_busmap(busmap_c),
+    )
 
-                sec_name = SecNames[sector].value
-                if sector.lower() == "transport":  # hack for now to get names to work
-                    sec_car = TRANSPORT_FUELS[end_use.lower()]
-                else:
-                    sec_car = SecCarriers[end_use].value
+    missing = pd.Index(busmap.dropna().unique()).difference(n.buses.index)
+    if len(missing):
+        raise ValueError(
+            f"Demand maps to absent clustered buses: {missing.tolist()[:10]}",
+        )
+
+    # Replace the initial electricity demand to avoid double counting.
+    n.mremove("Load", n.loads.index.copy())
+
+    for demand_file in demand_files:
+        sector, end_use = Path(demand_file).stem.split("_")
+        sec_name = SecNames[sector.upper()].value
+        end_use = end_use.upper().replace("-", "_")
+
+        if sector == "transport":
+            sec_car = TRANSPORT_FUELS[end_use.lower()]
+        else:
+            sec_car = SecCarriers[end_use].value
+
+        carrier = f"{sec_name}-{sec_car}"
+        df = aggregate_sector_demand(
+            pd.read_pickle(demand_file),
+            busmap,
+        )
+
+        if not pd.DatetimeIndex(df.index).equals(
+            n.snapshots.get_level_values("timestep"),
+        ):
+            raise ValueError(
+                f"Demand timestamps do not match network snapshots: {demand_file}",
+            )
+
+        attach_demand(n, df, carrier, suffix=f" {carrier}")
+        logger.info(
+            "%s %s demand added to network",
+            sector.upper(),
+            end_use,
+        )
 
-                carrier = f"{sec_name}-{sec_car}"
 
-                log_statement = f"{sector} {end_use} demand added to network"
+if __name__ == "__main__":
+    if "snakemake" not in globals():
+        snakemake = mock_snakemake("add_demand", interconnect="western")
 
-            else:
-                raise NotImplementedError
+    configure_logging(snakemake)
 
-            suffix = f"-{carrier}"
+    n = pypsa.Network(snakemake.input.network)
+    n.snapshots = get_multiindex_snapshots(
+        snakemake.params.snapshots,
+        snakemake.params.planning_horizons,
+    )
+    n.set_investment_periods(
+        periods=snakemake.params.planning_horizons,
+    )
 
-            df = pd.read_pickle(demand_file)
-            attach_demand(n, df, carrier, suffix)
-            logger.info(log_statement)
+    df = pd.read_csv(snakemake.input.demand, index_col=0)
+    attach_demand(n, df, carrier="AC", suffix="")
 
+    logger.info("Electricity demand added to network")
     n.export_to_netcdf(snakemake.output.network)
diff --git a/workflow/scripts/add_extra_components.py b/workflow/scripts/add_extra_components.py
index 29f6b3ac1..43a609f41 100644
--- a/workflow/scripts/add_extra_components.py
+++ b/workflow/scripts/add_extra_components.py
@@ -7,6 +7,7 @@
 import pandas as pd
 import pypsa
 from _helpers import calculate_annuity, configure_logging, load_costs
+from add_demand import attach_sector_demand
 from add_electricity import add_missing_carriers
 from constants import HOURS_PER_YEAR
 from eia import FuelCosts
@@ -1545,6 +1546,19 @@ def add_dac(n: pypsa.Network, config: dict, sector: bool):
     configure_logging(snakemake)
 
     n = pypsa.Network(snakemake.input.network)
+
+    # Introduce sector demand before demand-dependent network preparation.
+    if snakemake.config["scenario"]["sector"] not in ("E", ""):
+        attach_sector_demand(
+            n,
+            snakemake.input.sector_demand,
+            snakemake.input.busmap_s,
+            snakemake.input.busmap_c,
+        )
+
+        # Register sector carriers before assigning their metadata.
+        add_missing_carriers(n, n.loads.carrier)
+
     elec_config = snakemake.config["electricity"]
 
     costs_dict = {
diff --git a/workflow/scripts/build_cost_data.py b/workflow/scripts/build_cost_data.py
index 43a865df9..14a6a8ef2 100644
--- a/workflow/scripts/build_cost_data.py
+++ b/workflow/scripts/build_cost_data.py
@@ -63,17 +63,8 @@
 ]  # https://github.com/NREL/ReEDS-2.0/blob/e65ed5ed4ffff973071839481309f77d12d802cd/inputs/plant_characteristics/maxage.csv#L4
 
 
-def create_duckdb_instance():
-    """Set up DuckDB to read parquet files directly."""
-    duckdb.connect(database=":memory:", read_only=False)
-    # Install httpfs extension to access remote files if needed
-    duckdb.query("INSTALL httpfs;")
-
-
 def load_pudl_atb_data(parquet_path: str):
     """Loads ATB data directly from parquet files."""
-    create_duckdb_instance()
-
     query = f"""
     WITH finance_cte AS (
         SELECT
diff --git a/workflow/scripts/build_cutout.py b/workflow/scripts/build_cutout.py
index 37e3a62f4..085068897 100644
--- a/workflow/scripts/build_cutout.py
+++ b/workflow/scripts/build_cutout.py
@@ -139,4 +139,8 @@
     logging.info(f"Preparing cutout with parameters {cutout_params}.")
     features = cutout_params.pop("features", None)
     cutout = atlite.Cutout(snakemake.output[0], **cutout_params)
-    cutout.prepare(features=features)
+    cutout.prepare(
+        features=features,
+        monthly_requests=True,
+        concurrent_requests=False,
+    )
diff --git a/workflow/scripts/build_demand.py b/workflow/scripts/build_demand.py
index 035c8766b..9b5e2673d 100644
--- a/workflow/scripts/build_demand.py
+++ b/workflow/scripts/build_demand.py
@@ -17,7 +17,7 @@
 import pandas as pd
 import pypsa
 import xarray as xr
-from _helpers import configure_logging, get_multiindex_snapshots
+from _helpers import configure_logging, get_multiindex_snapshots, prepare_sector_demand
 from constants_sector import FIPS_2_STATE, NAICS, VMT_UNIT_CONVERSION, DemandFuels
 from eia import EnergyDemand, TransportationDemand
 
@@ -63,7 +63,7 @@ def _read(self) -> pd.DataFrame:
 
     def _write(self, demand: pd.DataFrame, zone: str, **kwargs) -> pd.DataFrame:
         """Delegate writing to the strategy."""
-        return self._write_strategy.dissagregate_demand(demand, zone, **kwargs)
+        return self._write_strategy.prepare_profiles(demand, zone, **kwargs)
 
     def prepare_demand(self, **kwargs) -> pd.DataFrame:
         """Read in and dissagregate demand."""
@@ -255,9 +255,6 @@ def _read_data(self) -> pd.DataFrame:
 
     def _read_census_data(self) -> pd.DataFrame:
         """Reads in census data for population weighting using parquet."""
-        duckdb.connect(database=":memory:", read_only=False)
-        duckdb.query("INSTALL httpfs;")
-
         parquet_path = snakemake.params.pudl_path
 
         sql = f"""
@@ -699,7 +696,7 @@ class ReadCliu(ReadStrategy):
     - MECS is used to scale county level data
 
     Sources:
-    - https://data.nrel.gov/submissions/97
+    - https://data.nlr.gov/submissions/97
     - https://www.eia.gov/consumption/manufacturing/data/2014/#r3 (Table 3.2)
     - https://loadshape.epri.com/enduse
     - https://github.com/NREL/Industry-Energy-Tool/
@@ -1590,7 +1587,7 @@ def _get_load_allocation_factor(
         """
         pass
 
-    def dissagregate_demand(
+    def prepare_profiles(
         self,
         df: pd.DataFrame,
         zone: str,
@@ -1600,53 +1597,77 @@ def dissagregate_demand(
         sns: pd.DatetimeIndex | None = None,
     ) -> pd.DataFrame:
         """
-        Public load dissagregation method.
-
-        df: pd.DataFrame
-            Demand dataframe
-        zone: str
-            Zones of demand ('ba', 'state', 'reeds')
-        sector: Optional[str | List[str]] = None,
-            Sectors to group
-        subsector: Optional[str | List[str]] = None,
-            Subsectors to group
-        fuel: Optional[str | List[str]] = None,
-            End use fules to group
-        sns: Optional[pd.DatetimeIndex] = None
-            Filter data over this period. If not provided, use network snapshots
-
-        Data is returned in the format of:
+        Prepare filtered zonal profiles and the original bus allocation table.
+
+        Parameters
+        ----------
+        df : pd.DataFrame
+            Source demand with source-zone columns and a MultiIndex containing
+            snapshot, sector, subsector, and fuel levels.
+        zone : {"ba", "state", "reeds"}
+            Spatial classification associating original network buses with
+            source zones. The selected allocation strategy must support it.
+        sector : str, list[str], or None, default None
+            Sector labels to retain. If None, no sector filter is applied.
+        subsector : str, list[str], or None, default None
+            Subsector labels to retain. If None, no subsector filter is applied.
+        fuel : str, list[str], or None, default None
+            End-use labels to retain. If None, no fuel filter is applied.
+        sns : pd.DatetimeIndex or None, default None
+            Source timestamps to retain. If None, no explicit timestamp filter
+            is applied.
 
-        |                     | BusName_1 | BusName_2 | ... | BusName_n |
-        |---------------------|-----------|-----------|-----|-----------|
-        | 2019-01-01 00:00:00 |    ###    |    ###    |     |    ###    |
-        | 2019-01-01 01:00:00 |    ###    |    ###    |     |    ###    |
-        | 2019-01-01 02:00:00 |    ###    |    ###    |     |    ###    |
-        | ...                 |           |           |     |    ###    |
-        | 2019-12-31 23:00:00 |    ###    |    ###    |     |    ###    |
+        Returns
+        -------
+        pd.DataFrame
+            Selected demand summed by timestamp, with source-zone columns:
+
+            |                     | Zone 1 | Zone 2 | ... | Zone n |
+            |---------------------|--------|--------|-----|--------|
+            | 2019-01-01 00:00:00 |  ###   |  ###   |     |  ###   |
+            | 2019-01-01 01:00:00 |  ###   |  ###   |     |  ###   |
+            | ...                 |        |        |     |        |
+
+        Raises
+        ------
+        AssertionError
+            If the zone classification is invalid or the source demand fails
+            the existing input-structure checks.
+
+        Notes
+        -----
+        The method sets ``self.allocation`` to a DataFrame indexed by original
+        buses, with ``zone`` and ``laf`` columns. Allocation factors are
+        calculated using the existing population or industrial strategy.
+
+        If filtering produces no demand, the existing zero-demand replacement
+        based on network snapshots is used.
+
+        Returned values have not been multiplied by nodal allocation factors.
+        Investment-year projection, unit conversion, and final rounding are
+        handled by subsequent processing steps.
         """
-        # 'state' is states based on power regions
-        # 'full_state' is actual geographic boundaries
         assert zone in ("ba", "state", "reeds")
         self._check_datastructure(df)
 
-        # get zone area demand for specific sector and fuel
+        # Select and aggregate source demand without expanding it to buses.
         demand = self._filter_demand(df, sector, subsector, fuel, sns)
         demand = self._group_demand(demand)
         if demand.empty:
             demand = self._make_empty_demand(columns=df.columns)
 
-        # assign buses to dissagregation zone
+        # Preserve the allocation calculated on the original network.
         dissagregation_zones = self._get_load_dissagregation_zones(zone)
-
-        # get implementation specific dissgregation factors
-        laf = self._get_load_allocation_factor(df=dissagregation_zones, zone=zone)
-
-        # disaggregate load to buses
+        laf = self._get_load_allocation_factor(
+            df=dissagregation_zones,
+            zone=zone,
+        )
         zone_data = dissagregation_zones.to_frame(name="zone").join(
             laf.to_frame(name="laf"),
         )
-        return self._disaggregate_demand_to_buses(demand, zone_data)
+
+        self.allocation = zone_data
+        return demand
 
     def _get_load_dissagregation_zones(self, zone: str) -> pd.Series:
         """Map each bus to the load dissagregation zone (states, ba, ...)."""
@@ -1814,7 +1835,7 @@ class WriteIndustrial(WriteStrategy):
     """
     Based on county level energy use from 2014.
 
-    https://data.nrel.gov/submissions/97
+    https://data.nlr.gov/submissions/97
     """
 
     def __init__(self, n: pypsa.Network, filepath: str) -> None:
@@ -1978,37 +1999,134 @@ def __init__(
         else:
             self.scaler = None
 
-    def format_demand(self, df: pd.DataFrame, sector: str, **kwargs) -> pd.DataFrame:
-        """Public method to format demand ready to be ingested into the model."""
+    def format_profiles(
+        self,
+        df: pd.DataFrame,
+        sector: str,
+        **kwargs,
+    ) -> tuple[pd.DataFrame, pd.Series]:
+        """
+        Align demand profiles to investment years and calculate growth factors.
+
+        Parameters
+        ----------
+        df : pd.DataFrame
+            Demand indexed by timestamps from the available source years.
+            Columns may identify source zones or original network buses.
+        sector : str
+            Category passed to the configured growth scaler. Transport demand
+            uses the corresponding vehicle or transport category.
+        **kwargs
+            Additional keyword arguments retained for call compatibility.
+            They are not used by this implementation.
+
+        Returns
+        -------
+        profiles : pd.DataFrame
+            Profiles reindexed to the configured investment years, with
+            unchanged source values and unchanged columns.
+        growth : pd.Series
+            Multiplicative growth factors indexed exactly like ``profiles``.
+            Profiles already supplied for an investment year receive 1.0.
+
+        Raises
+        ------
+        AssertionError
+            If the existing scaling check requires a DemandScaler but none
+            is available, no source year precedes or matches the first
+            investment year, or the assembled profile length differs from
+            the snapshot count.
+
+        Notes
+        -----
+        Each missing investment year uses the latest available source year
+        that does not exceed it. For example, a 2020 profile can provide a
+        missing 2030 profile, while an available 2040 profile can provide a
+        missing 2045 profile.
+
+        Growth factors are returned separately so sector demand can retain
+        the multiplication order of nodal allocation followed by demand
+        growth. This method performs neither unit conversion nor rounding.
+
+        If source years already match the investment periods, the supplied
+        profiles are returned directly with unit growth factors.
+        """
         if self.need_scaling(df):
             assert isinstance(self.scaler, DemandScaler)
 
         demand_periods = df.index.year.unique().to_list()
         if demand_periods == self.investment_periods:
             logger.info("No demand formatting required")
-            return df
+            return df, pd.Series(1.0, index=df.index, name="growth")
 
-        # need a starting reference year for scaling
         assert min(demand_periods) <= min(self.investment_periods)
-        demand_per_period = []
+
+        profiles_per_period = []
+        growth_per_period = []
+
         for investment_year in self.investment_periods:
-            formatted_demand = df[df.index.year == investment_year]
-            if not formatted_demand.empty:
-                demand_per_period.append(formatted_demand)
-            else:
-                nearest_year = max([x for x in demand_periods if x <= investment_year])
-                formatted_demand = df[df.index.year == nearest_year]
-                demand_per_period.append(
-                    self.scaler.scale(
-                        formatted_demand,
-                        nearest_year,
-                        investment_year,
-                        sector,
-                    ),
+            profile = df[df.index.year == investment_year]
+            factor = 1.0
+
+            if profile.empty:
+                nearest_year = max(x for x in demand_periods if x <= investment_year)
+                profile = self.scaler.reindex(
+                    df[df.index.year == nearest_year],
+                    investment_year,
                 )
-        demand = pd.concat(demand_per_period)
-        assert len(demand) == len(self.sns)
-        return demand
+                factor = self.scaler.get_growth(
+                    nearest_year,
+                    investment_year,
+                    sector,
+                )
+
+            profiles_per_period.append(profile)
+            growth_per_period.append(
+                pd.Series(factor, index=profile.index, name="growth"),
+            )
+
+        profiles = pd.concat(profiles_per_period)
+        growth = pd.concat(growth_per_period)
+
+        assert len(profiles) == len(self.sns)
+        return profiles, growth
+
+    def format_demand(
+        self,
+        df: pd.DataFrame,
+        sector: str,
+        **kwargs,
+    ) -> pd.DataFrame:
+        """
+        Return demand projected onto the configured investment periods.
+
+        Parameters
+        ----------
+        df : pd.DataFrame
+            Demand indexed by source timestamps. Columns identify the spatial
+            units whose demand is being projected.
+        sector : str
+            Demand category used by the configured growth scaler.
+        **kwargs
+            Additional arguments forwarded to ``format_profiles``.
+
+        Returns
+        -------
+        pd.DataFrame
+            Investment-year demand after applying the corresponding growth
+            factors. Columns are unchanged and rows use a DatetimeIndex.
+
+        Notes
+        -----
+        Electricity demand calls this method after nodal disaggregation,
+        preserving the existing allocation-before-growth order.
+
+        Sector demand instead consumes the separate outputs of
+        ``format_profiles`` so nodal allocation and growth can be applied
+        together during aggregation onto the clustered network.
+        """
+        profiles, growth = self.format_profiles(df, sector, **kwargs)
+        return profiles.mul(growth, axis=0)
 
     def need_scaling(self, df: pd.DataFrame) -> bool:
         """Checks if any demand needs to be scaled."""
@@ -2022,9 +2140,7 @@ def assign_scaler(self):  # type DemandScaler
             assert self.api, "Must provide eia api key"
             return AeoEnergyScaler(self.api)
         elif self.scaling_method == "aeo_electricity":
-            assert self.filepath.startswith(
-                "s3://pudl.catalyst.coop/",
-            ), "Must provide pudl S3 URL (s3://pudl.catalyst.coop/...)"
+            assert self.filepath, "Must provide a local PUDL directory"
             return AeoElectricityScaler(self.filepath)
         elif self.scaling_method == "efs":
             assert self.filepath.endswith(".csv"), "Must provide EFS.csv data"
@@ -2126,9 +2242,6 @@ def get_projections(self) -> pd.DataFrame:
         | 2049 |  ###  |  ###  |
         | 2050 |  ###  |  ###  |
         """
-        duckdb.connect(database=":memory:", read_only=False)
-        duckdb.query("INSTALL httpfs;")
-
         query = f"""
         SELECT
             projection_year,
@@ -2610,8 +2723,14 @@ def _get_sector_fuels(end_use: str, vehicle: str | None = None) -> list[str]:
     # extract demand based on strategies
     # this is raw demand, not scaled or garunteed to align to network snapshots
 
-    if end_use == "power":  # only one demand for electricity only studies
-        demand = demand_converter.prepare_demand(sns=sns)  # pd.DataFrame
+    if end_use == "power":
+        demand = demand_converter.prepare_demand(sns=sns)
+
+        # Electricity demand is attached before network simplification.
+        demand = writer._disaggregate_demand_to_buses(
+            demand,
+            writer.allocation,
+        )
         demands = {"electricity": demand}
     else:
         fuels = _get_sector_fuels(end_use, vehicle)
@@ -2636,13 +2755,30 @@ def _get_sector_fuels(end_use: str, vehicle: str | None = None) -> list[str]:
     )
 
     formatted_demand = {}
-    if end_use == "transport":  # transport fuel is actually vehicle types
-        for fuel, demand in demands.items():
-            vmt_conversion = VMT_UNIT_CONVERSION.get(fuel, 1)  # for buses
-            formatted_demand[fuel] = demand_formatter.format_demand(demand, fuel).mul(vmt_conversion)
-    else:
-        for fuel, demand in demands.items():
-            formatted_demand[fuel] = demand_formatter.format_demand(demand, end_use)
+
+    for fuel, demand in demands.items():
+        if end_use == "power":
+            formatted_demand[fuel] = demand_formatter.format_demand(
+                demand,
+                end_use,
+            )
+            continue
+
+        # Transport growth projections are selected by transport category.
+        scaling_sector = fuel if end_use == "transport" else end_use
+
+        profiles, growth = demand_formatter.format_profiles(
+            demand,
+            scaling_sector,
+        )
+        unit_conversion = VMT_UNIT_CONVERSION.get(fuel, 1) if end_use == "transport" else 1
+
+        formatted_demand[fuel] = prepare_sector_demand(
+            profiles,
+            writer.allocation,
+            growth,
+            unit_conversion,
+        )
 
     # electricity sector study
     if end_use == "power":
@@ -2659,6 +2795,4 @@ def _get_sector_fuels(end_use: str, vehicle: str | None = None) -> list[str]:
                     out_f = snakemake.output[0]
                 else:
                     raise KeyError(e)
-            formatted_demand[fuel].round(4).to_pickle(
-                out_f,
-            )
+            pd.to_pickle(formatted_demand[fuel], out_f)
diff --git a/workflow/scripts/build_natural_gas.py b/workflow/scripts/build_natural_gas.py
index 10cde0b04..f2669c1e4 100644
--- a/workflow/scripts/build_natural_gas.py
+++ b/workflow/scripts/build_natural_gas.py
@@ -853,7 +853,7 @@ def _expand_costs(self, n: pypsa.Network, costs: pd.DataFrame) -> pd.DataFrame:
         expanded_costs = []
         for invesetment_period in n.investment_periods:
             # reindex to match any tsa
-            cost = costs.copy()
+            cost = costs.loc[costs.index.year == self.year].copy()
             cost.index = cost.index.map(lambda x: x.replace(year=invesetment_period))
             cost = cost.reindex(n.snapshots.get_level_values(1), method="nearest")
             # set investment periods
diff --git a/workflow/scripts/build_powerplants.py b/workflow/scripts/build_powerplants.py
index 29c509554..aa64c46f7 100644
--- a/workflow/scripts/build_powerplants.py
+++ b/workflow/scripts/build_powerplants.py
@@ -12,11 +12,6 @@
 logger = logging.getLogger(__name__)
 
 
-def initialize_duckdb():
-    duckdb.connect(database=":memory:", read_only=False)
-    duckdb.query("INSTALL httpfs;")
-
-
 def load_eia_operable_data(parquet_path: str):
     """Queries the parquet files directly for operable plant data."""
     return duckdb.query(
@@ -694,9 +689,9 @@ def apply_cems_heat_rates(plants, crosswalk_fn, cems_fn):
     start_date = f"{data_year}-01-01"
     end_date = f"{data_year + 1}-01-01"
 
-    initialize_duckdb()
-    eia_data_operable = load_eia_operable_data(snakemake.params.pudl_path)
-    heat_rates = load_heat_rates_data(snakemake.params.pudl_path, start_date, end_date)
+    pudl_path = snakemake.params.pudl_path
+    eia_data_operable = load_eia_operable_data(pudl_path)
+    heat_rates = load_heat_rates_data(pudl_path, start_date, end_date)
 
     eia_data_operable = merge_fc_hr_data(
         eia_data_operable,
diff --git a/workflow/scripts/build_sector_costs.py b/workflow/scripts/build_sector_costs.py
index 083931ca2..264f141cc 100644
--- a/workflow/scripts/build_sector_costs.py
+++ b/workflow/scripts/build_sector_costs.py
@@ -2,7 +2,7 @@
 Builds costs data from EFS studies.
 
 https://www.nrel.gov/docs/fy18osti/70485.pdf
-https://data.nrel.gov/submissions/78
+https://data.nlr.gov/submissions/78
 """
 
 from abc import ABC, abstractmethod
@@ -22,7 +22,7 @@
 
 # fixed maintenance costs in %/year
 # From https://www.nrel.gov/docs/fy18osti/71500.pdf
-# Data in Sheet "13" at https://data.nrel.gov/submissions/93
+# Data in Sheet "13" at https://data.nlr.gov/submissions/93
 # manually claculated as fixed/capex
 ICE_MAINTENANCE_COSTS = {  # units of % / year
     "light_duty_cars": 11.34,
@@ -43,7 +43,7 @@
 # maintenance costs in $/mile. This gave very high maintenance costs
 # when using exogenous lifetime miles assumption.
 # From https://www.nrel.gov/docs/fy18osti/71500.pdf
-# Data in Sheet "A2" at https://data.nrel.gov/submissions/93
+# Data in Sheet "A2" at https://data.nlr.gov/submissions/93
 """
 ICE_MAINTENANCE_COSTS = {  # units of $ / miles
     "light_duty_cars": 0.033,
@@ -96,7 +96,7 @@ class EfsTechnologyData:
 
     args:
         file_path: str
-            Path to this file https://data.nrel.gov/submissions/78 which is the
+            Path to this file https://data.nlr.gov/submissions/78 which is the
             data from this report https://www.nrel.gov/docs/fy18osti/70485.pdf
 
     Buildings sector will return both commercial and residential data
@@ -254,7 +254,7 @@ def get_capex(self):  # noqa: D102
         df = df[
             (df.Sector == "Transportation") & (df["EFS Case"] == self.efs_case) & (df.Metric == "Capital Cost")
         ].copy()
-        source = "NREL EFS at https://data.nrel.gov/submissions/78"
+        source = "NREL EFS at https://data.nlr.gov/submissions/78"
         description = "Supplemented with NREL ATB"
         df["Metric"] = "investment"
         df = self._format_data_structure(df, source=source, description=description)
@@ -278,7 +278,7 @@ def get_efficiency(self):
         df = df[
             (df.Sector == "Transportation") & (df["EFS Case"] == self.efs_case) & (df.Metric == "Main Efficiency")
         ].copy()
-        source = "NREL EFS at https://data.nrel.gov/submissions/78"
+        source = "NREL EFS at https://data.nlr.gov/submissions/78"
         df["Metric"] = "efficiency"
         df = self._format_data_structure(df, source=source, description="")
         df = self._correct_efficiency_units(df)
@@ -329,7 +329,7 @@ class EfsIceTransportationData:
     Only contains ICE vehicles, as this data is manually scraped.
 
     - See Table 5 in https://www.nrel.gov/docs/fy18osti/70485.pdf for the sources
-    - See this file for the raw data https://data.nrel.gov/submissions/93
+    - See this file for the raw data https://data.nlr.gov/submissions/93
 
     Args:
         file_path: str
@@ -511,7 +511,7 @@ def get_capex(self):  # noqa: D102
             & (df.Metric == "Installed Cost")
             & (df.Units.str.startswith("2016$/kBtu/hr"))
         ].copy()
-        source = "NREL EFS at https://data.nrel.gov/submissions/78"
+        source = "NREL EFS at https://data.nlr.gov/submissions/78"
         description = ""
         df["Metric"] = "investment"
         df = self._format_data_structure(df, source=source, description=description)
@@ -529,7 +529,7 @@ def get_lifetime(self):  # noqa: D102
     def get_efficiency(self):  # noqa: D102
         df = self.data.copy()
         df = df[(df.Sector == "Buildings") & (df["EFS Case"] == self.efs_case) & (df.Metric == "Efficiency")].copy()
-        source = "NREL EFS at https://data.nrel.gov/submissions/78"
+        source = "NREL EFS at https://data.nlr.gov/submissions/78"
         df["Metric"] = "efficiency"
         df["Units"] = "per unit"
         df = self._format_data_structure(df, source=source, description="")
diff --git a/workflow/scripts/build_stock_data.py b/workflow/scripts/build_stock_data.py
index c97bb6f07..5a3b966f5 100644
--- a/workflow/scripts/build_stock_data.py
+++ b/workflow/scripts/build_stock_data.py
@@ -832,7 +832,7 @@ def add_brownfield_lpg(
         # existing stock efficiencies taken from 2016 EFS Technology data
         # This is consistent with where future efficiencies are taken from
         # Historical uses lowest EIA 2017 case where available
-        # https://data.nrel.gov/submissions/93
+        # https://data.nlr.gov/submissions/93
         # https://www.nrel.gov/docs/fy18osti/70485.pdf
 
         ## Efficiencies of existing stock are quite sensitive! ##
diff --git a/workflow/scripts/cluster_network.py b/workflow/scripts/cluster_network.py
index 0d29bd714..4ca419cfa 100644
--- a/workflow/scripts/cluster_network.py
+++ b/workflow/scripts/cluster_network.py
@@ -985,6 +985,8 @@ def calibrate_tamu_transmission_capacity(
 
     n = pypsa.Network(snakemake.input.network)
 
+    input_buses = n.buses.index.copy()
+
     n.set_investment_periods(
         periods=snakemake.params.planning_horizons,
     )
@@ -1265,11 +1267,9 @@ def calibrate_tamu_transmission_capacity(
 
     clustering.network.export_to_netcdf(snakemake.output.network)
 
-    for attr in (
-        "busmap",
-        "linemap",
-    ):  # also available: linemap_positive, linemap_negative
-        getattr(clustering, attr).to_csv(snakemake.output[attr])
+    # Retain explicit null destinations for buses removed during processing.
+    clustering.busmap.reindex(input_buses).to_csv(snakemake.output.busmap)
+    clustering.linemap.to_csv(snakemake.output.linemap)
 
     cluster_regions((clustering.busmap,), snakemake.input, snakemake.output)
     n.consistency_check()
diff --git a/workflow/scripts/eia.py b/workflow/scripts/eia.py
index 5d27c9145..3a264e77b 100644
--- a/workflow/scripts/eia.py
+++ b/workflow/scripts/eia.py
@@ -1607,7 +1607,7 @@ def extract_state(description: str) -> str:
     @staticmethod
     def map_state_names(state: str) -> str:
         """Maps state name to code."""
-        return "U.S." if state == "U.S. Total" else STATE_CODES[state]
+        return "U.S." if state in {"U.S.", "U.S. Total"} else STATE_CODES[state]
 
 
 class _GasProduction(DataExtractor):
diff --git a/workflow/scripts/opts/sector.py b/workflow/scripts/opts/sector.py
index 186ff19a7..1a1f05e9f 100644
--- a/workflow/scripts/opts/sector.py
+++ b/workflow/scripts/opts/sector.py
@@ -420,7 +420,8 @@ def add_export_limits(n, data, constraint, multiplier=None):
     # add domestic limits
 
     trade = Trade("gas", False, "exports", year, api).get_data()
-    trade = _format_data(trade, " trade")
+    trade = trade.loc[trade.index == year]
+    trade = _format_data(trade, " trade").groupby(level=0).sum()
 
     add_import_limits(n, trade, "min", import_min)
     add_export_limits(n, trade, "min", export_min)
@@ -433,7 +434,8 @@ def add_export_limits(n, data, constraint, multiplier=None):
     # add international limits
 
     trade = Trade("gas", True, "exports", year, api).get_data()
-    trade = _format_data(trade, " trade")
+    trade = trade.loc[trade.index == year]
+    trade = _format_data(trade, " trade").groupby(level=0).sum()
 
     add_import_limits(n, trade, "min", import_min)
     add_export_limits(n, trade, "min", export_min)
@@ -595,7 +597,7 @@ def add_ev_generation_constraint(n, config, snakemake):
 
     Default limits taken from:
     - (Fig ES2) https://www.nrel.gov/docs/fy18osti/71500.pdf
-    - (Sheet 6.3 - high case) https://data.nrel.gov/submissions/90
+    - (Sheet 6.3 - high case) https://data.nlr.gov/submissions/90
     """
     mode_mapper = {
         "light_duty": "lgt",
diff --git a/workflow/scripts/plot_network_maps.py b/workflow/scripts/plot_network_maps.py
index 3e6c5c9a5..55b54bab9 100644
--- a/workflow/scripts/plot_network_maps.py
+++ b/workflow/scripts/plot_network_maps.py
@@ -466,9 +466,22 @@ def plot_capacity_map_by_horizon(
         bus_scale = get_bus_scale(interconnect) if interconnect else 1
         line_scale = get_line_scale(interconnect) if interconnect else 1
 
-        if bus_values.empty and kind == "new":
+        has_new_capacity = (
+            bus_values.fillna(0).gt(0).any() or line_values.fillna(0).gt(0).any() or link_values.fillna(0).gt(0).any()
+        )
+
+        if kind == "new" and not has_new_capacity:
+            logger.info("No new capacity found; creating an empty new-capacity map.")
             fig, ax = plt.subplots(figsize=(10, 10))
-            ax.text(0.5, 0.5, "No new capacity built", ha="center", va="center", fontsize=14)
+            ax.text(
+                0.5,
+                0.5,
+                "No new capacity built",
+                ha="center",
+                va="center",
+                fontsize=14,
+                transform=ax.transAxes,
+            )
             ax.set_title(title, fontsize=TITLE_SIZE, pad=20)
             ax.axis("off")
         else:
@@ -551,6 +564,25 @@ def plot_capacity_map_by_horizon(
             remove_sector_buses(bus_values).groupby(["bus", "carrier"]).sum() if not bus_values.empty else bus_values
         )
 
+        has_new_capacity = (
+            bus_values.fillna(0).gt(0).any() or line_values.fillna(0).gt(0).any() or link_values.fillna(0).gt(0).any()
+        )
+
+        if kind == "new" and not has_new_capacity:
+            logger.info("No new capacity found for horizon %s.", horizon)
+            ax.text(
+                0.5,
+                0.5,
+                "No new capacity built",
+                ha="center",
+                va="center",
+                fontsize=14,
+                transform=ax.transAxes,
+            )
+            ax.set_title(f"{horizon}", fontsize=TITLE_SIZE)
+            ax.axis("off")
+            continue
+
         artifacts = _plot_capacity_on_ax(
             n,
             bus_values,
diff --git a/workflow/scripts/simplify_network.py b/workflow/scripts/simplify_network.py
index 44543a0eb..6b74acda7 100644
--- a/workflow/scripts/simplify_network.py
+++ b/workflow/scripts/simplify_network.py
@@ -297,6 +297,9 @@ def assign_line_lengths(n, line_length_factor):
         )
         n = clustering.network
 
+        # Extend the original bus mapping through simplification clustering.
+        busmaps = busmaps.map(clustering.busmap)
+
         cluster_regions((clustering.busmap,), snakemake.input, snakemake.output)
     else:
         for which in ("regions_onshore", "regions_offshore"):  # pass through regions
@@ -305,4 +308,6 @@ def assign_line_lengths(n, line_length_factor):
 
     update_p_nom_max(n)
 
+    # Preserve original bus assignments for sector-demand aggregation.
+    busmaps.rename("bus").to_csv(snakemake.output.busmap)
     n.export_to_netcdf(snakemake.output[0])
diff --git a/workflow/scripts/zenodo_downloader.py b/workflow/scripts/zenodo_downloader.py
index 4d7df7d9f..1a73d4ad4 100644
--- a/workflow/scripts/zenodo_downloader.py
+++ b/workflow/scripts/zenodo_downloader.py
@@ -1,6 +1,8 @@
 """Download scenarios from Zenodo."""
 
+import shutil
 from pathlib import Path
+from zipfile import BadZipFile, ZipFile
 
 import requests
 
@@ -42,7 +44,10 @@ def get_record_metadata(self, record_id):
         url = f"https://zenodo.org/api/records/{record_id}"
 
         try:
-            response = requests.get(url)
+            response = requests.get(
+                url,
+                timeout=(30, 120),
+            )
             response.raise_for_status()
             metadata = response.json()
             self._metadata_cache[record_id] = metadata
@@ -68,7 +73,10 @@ def download_scenario_file(self, scenario_final, scenario, filename, force_redow
         force_redownload : bool, optional
             If True, re-download the file even if it exists locally. Default is False.
         """
-        (self.download_dir / "zenodo" / scenario).mkdir(exist_ok=True)
+        (self.download_dir / "zenodo" / scenario).mkdir(
+            parents=True,
+            exist_ok=True,
+        )
         local_filepath = f"{self.download_dir}/zenodo/{scenario}/{filename}"
 
         # Check if file already exists locally and skip Zenodo
@@ -82,12 +90,9 @@ def download_scenario_file(self, scenario_final, scenario, filename, force_redow
             record_id = self.scenario_records.get(scenario_final)
 
             if not record_id:
-                print(f"No record ID found for scenario: {scenario_final}")
-                print("Available scenarios with record IDs:")
-                for scenario, rec_id in self.scenario_records.items():
-                    if rec_id is not None:
-                        print(f"  - {scenario} (ID: {rec_id})")
-                return None
+                raise ValueError(
+                    f"No Zenodo record ID configured for scenario: {scenario_final}",
+                )
 
             return self._download_file(record_id, filename, Path(local_filepath), force_redownload)
 
@@ -144,71 +149,126 @@ def download_by_record_id(self, record_id, filename, force_redownload=False):
         # Only proceed with download if needed
         return self._download_file(record_id, filename, Path(local_filepath), force_redownload)
 
-    def _download_file(self, record_id, filename, local_filepath, force_redownload=False):
-        """
-        Internal method to download a file from Zenodo.
-
-        This is only called after confirming the file doesn't exist locally.
-        """
-        # Ensure directory exists
+    def _download_file(
+        self,
+        record_id,
+        filename,
+        local_filepath,
+        force_redownload=False,
+    ):
+        """Download a file directly or extract it from a Zenodo ZIP archive."""
         local_filepath.parent.mkdir(parents=True, exist_ok=True)
 
-        # Get record metadata
         metadata = self.get_record_metadata(record_id)
         if not metadata:
-            return None
-
-        # Find the specific file
-        target_file = None
-        for file_info in metadata.get("files", []):
-            if file_info["key"] == filename:
-                target_file = file_info
-                break
-
-        if not target_file:
-            print(f"File '{filename}' not found in record {record_id}")
-            print("Available files:")
-            for file_info in metadata.get("files", []):
-                print(f"  - {file_info['key']}")
-            return None
-
-        # Download the file
-        download_url = target_file["links"]["self"]
-        file_size_mb = target_file["size"] / (1024 * 1024)
+            raise RuntimeError(
+                f"Could not retrieve metadata for Zenodo record {record_id}",
+            )
 
-        print(f"Downloading {filename} from record {record_id}...")
-        print(f"Size: {file_size_mb:.1f} MB")
-        print(f"Saving to: {local_filepath}")
+        files = metadata.get("files", [])
 
-        try:
-            response = requests.get(download_url, stream=True)
-            response.raise_for_status()
+        # Some Zenodo records expose the requested NetCDF file directly.
+        target_file = next(
+            (file_info for file_info in files if file_info["key"] == filename),
+            None,
+        )
 
-            total_size = int(response.headers.get("content-length", 0))
-            downloaded_size = 0
+        if target_file is not None:
+            self._download_url(
+                target_file["links"]["self"],
+                local_filepath,
+                target_file["size"],
+            )
+            return str(local_filepath)
 
-            with open(local_filepath, "wb") as f:
-                for chunk in response.iter_content(chunk_size=8192):
-                    if chunk:
-                        f.write(chunk)
-                        downloaded_size += len(chunk)
+        # Future renewable datasets are published as ZIP archives.
+        zip_file = next(
+            (file_info for file_info in files if file_info["key"].lower().endswith(".zip")),
+            None,
+        )
+
+        if zip_file is None:
+            available_files = ", ".join(file_info["key"] for file_info in files)
+            raise FileNotFoundError(
+                f"Neither '{filename}' nor a ZIP archive was found in "
+                f"Zenodo record {record_id}. Available files: {available_files}",
+            )
 
-                        # Show progress for large files
-                        if total_size > 10 * 1024 * 1024:  # Show progress for files > 10MB
-                            progress = (downloaded_size / total_size) * 100
-                            print(f"\rProgress: {progress:.1f}%", end="", flush=True)
+        archive_path = local_filepath.parent / zip_file["key"]
 
-            if total_size > 10 * 1024 * 1024:
-                print()  # New line after progress
+        if force_redownload or not archive_path.exists():
+            self._download_url(
+                zip_file["links"]["self"],
+                archive_path,
+                zip_file["size"],
+            )
 
-            print(f"Successfully downloaded {filename}")
-            return str(local_filepath)
+        try:
+            with ZipFile(archive_path) as archive:
+                matches = [member for member in archive.namelist() if Path(member).name == filename]
+
+                if len(matches) != 1:
+                    raise FileNotFoundError(
+                        f"Expected exactly one '{filename}' in '{archive_path.name}', found {len(matches)}",
+                    )
+
+                with (
+                    archive.open(matches[0]) as source,
+                    local_filepath.open("wb") as destination,
+                ):
+                    shutil.copyfileobj(
+                        source,
+                        destination,
+                        length=1024 * 1024,
+                    )
+
+        except BadZipFile:
+            local_filepath.unlink(missing_ok=True)
+            archive_path.unlink(missing_ok=True)
+            raise
+
+        except OSError:
+            local_filepath.unlink(missing_ok=True)
+            raise
+
+        return str(local_filepath)
+
+    def _download_url(self, download_url, destination, expected_size):
+        """Download a URL to a local path using a temporary partial file."""
+        partial_path = destination.with_suffix(
+            destination.suffix + ".part",
+        )
+
+        print(
+            f"Downloading {destination.name} ({expected_size / 1024**3:.1f} GiB)...",
+        )
 
-        except requests.exceptions.RequestException as e:
-            print(f"Download failed: {e}")
-            if Path(local_filepath).exists():
-                Path(local_filepath).unlink()  # Remove partial file
-            return None
+        try:
+            with requests.get(
+                download_url,
+                stream=True,
+                timeout=(30, 300),
+            ) as response:
+                response.raise_for_status()
+
+                with partial_path.open("wb") as output:
+                    for chunk in response.iter_content(
+                        chunk_size=1024 * 1024,
+                    ):
+                        if chunk:
+                            output.write(chunk)
+
+            actual_size = partial_path.stat().st_size
+            if expected_size and actual_size != expected_size:
+                raise OSError(
+                    f"Incomplete download for {destination.name}: expected {expected_size} bytes, got {actual_size}",
+                )
+
+            partial_path.replace(destination)
+
+        except Exception:
+            partial_path.unlink(missing_ok=True)
+            raise
 
     def list_available_files(self, scenario_name):
         """List all available files in a scenario dataset."""
@@ -252,16 +312,32 @@ def get_available_scenarios(self):
         return available
 
 
-def download_scenario_file(scenario_final, scenario, filename, download_dir="./data/zenodo"):
-    """Quick function to download a single file from a scenario."""
+def download_scenario_file(
+    scenario_final,
+    scenario,
+    filename,
+    download_dir="./data",
+):
+    """Download a single file from a configured scenario."""
     downloader = ZenodoScenarioDownloader(download_dir)
-    return downloader.download_scenario_file(scenario_final, scenario, filename)
-
-
-def download_by_record_id(record_id, filename, download_dir="./data/zenodo"):
-    """Quick function to download a file directly by record ID."""
+    return downloader.download_scenario_file(
+        scenario_final,
+        scenario,
+        filename,
+    )
+
+
+def download_by_record_id(
+    record_id,
+    filename,
+    download_dir="./data",
+):
+    """Download a file directly using its Zenodo record ID."""
     downloader = ZenodoScenarioDownloader(download_dir)
-    return downloader.download_by_record_id(record_id, filename)
+    return downloader.download_by_record_id(
+        record_id,
+        filename,
+    )
 
 
 def list_available_scenarios():