Add industry sites

config/config.default.yaml

@@ -435,6 +435,7 @@ solar_thermal:
 existing_capacities:
   grouping_years_power: [1920, 1950, 1955, 1960, 1965, 1970, 1975, 1980, 1985, 1990, 1995, 2000, 2005, 2010, 2015, 2020, 2025]
   grouping_years_heat: [1980, 1985, 1990, 1995, 2000, 2005, 2010, 2015, 2019] # heat grouping years >= baseyear will be ignored
+  grouping_years_industry: [1995, 2000, 2005, 2010, 2015, 2020, 2025]
   threshold_capacity: 10
   default_heating_lifetime: 20
   conventional_carriers:
@@ -770,7 +771,7 @@ sector:
     var_cf: true
     sustainability_factor: 0.0025
   solid_biomass_import:
-    enable: false
+    enable: true
     price: 54 #EUR/MWh
     max_amount: 1390 # TWh
     upstream_emissions_factor: .1 #share of solid biomass CO2 emissions at full combustion
@@ -937,7 +938,7 @@ clustering:
       ramp_limit_down: max
   temporal:
     resolution_elec: false
-    resolution_sector: false
+    resolution_sector: 8760H
 
 # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#adjustments
 adjustments:

rules/solve_myopic.smk

@@ -11,6 +11,8 @@ rule add_existing_baseyear:
         costs=config_provider("costs"),
         heat_pump_sources=config_provider("sector", "heat_pump_sources"),
         energy_totals_year=config_provider("energy", "energy_totals_year"),
+        countries=config_provider("countries"),
+        MWh_NH3_per_tNH3=config_provider("industry", "MWh_NH3_per_tNH3"),
     input:
         network=resources(
             "networks/base_s_{clusters}_{opts}_{sector_opts}_{planning_horizons}.nc"
@@ -29,6 +31,9 @@ rule add_existing_baseyear:
             "existing_heating_distribution_base_s_{clusters}_{planning_horizons}.csv"
         ),
         heating_efficiencies=resources("heating_efficiencies.csv"),
+        regions_onshore=resources("regions_onshore_base_s_{clusters}.geojson"),
+        ammonia="data/ammonia_plants.csv",
+        isi_database="data/1-s2.0-S0196890424010586-mmc2.xlsx",
     output:
         resources(
             "networks/base_s_{clusters}_{opts}_{sector_opts}_{planning_horizons}_brownfield.nc"

rules/solve_perfect.smk

@@ -9,6 +9,8 @@ rule add_existing_baseyear:
         costs=config_provider("costs"),
         heat_pump_sources=config_provider("sector", "heat_pump_sources"),
         energy_totals_year=config_provider("energy", "energy_totals_year"),
+        countries=config_provider("countries"),
+        MWh_NH3_per_tNH3=config_provider("industry", "MWh_NH3_per_tNH3"),
     input:
         network=resources(
             "networks/base_s_{clusters}_{opts}_{sector_opts}_{planning_horizons}.nc"
@@ -28,6 +30,9 @@ rule add_existing_baseyear:
         ),
         existing_heating="data/existing_infrastructure/existing_heating_raw.csv",
         heating_efficiencies=resources("heating_efficiencies.csv"),
+        regions_onshore=resources("regions_onshore_base_s_{clusters}.geojson"),
+        ammonia="data/ammonia_plants.csv",
+        isi_database="data/1-s2.0-S0196890424010586-mmc2.xlsx",
     output:
         resources(
             "networks/base_s_{clusters}_{opts}_{sector_opts}_{planning_horizons}_brownfield.nc"

scripts/add_existing_baseyear.py

@@ -11,6 +11,7 @@
 from types import SimpleNamespace
 
 import country_converter as coco
+import geopandas as gpd
 import numpy as np
 import pandas as pd
 import powerplantmatching as pm
@@ -706,6 +707,161 @@ def add_heating_capacities_installed_before_baseyear(
             )
 
 
+def prepare_plant_data(
+    regions_fn: str,
+    isi_database: str,
+) -> tuple[pd.DataFrame, gpd.GeoDataFrame]:
+    """
+    Reads in the Fraunhofer ISI database with high resolution plant data and maps them to the bus regions.
+    Returns the database as df as well as the regions as gdf.
+
+    Parameters
+    ----------
+    regions_fn : str
+        path to the onshore regions file
+    isi_database: str
+        path to the fraunhofer isi database
+    """
+    # add existing industry
+    regions = gpd.read_file(regions_fn).set_index("name")
+
+    isi_data = pd.read_excel(isi_database, sheet_name="Database", index_col=1)
+    # assign bus region to each plant
+    geometry = gpd.points_from_xy(isi_data["Longitude"], isi_data["Latitude"])
+    plant_data = gpd.GeoDataFrame(isi_data, geometry=geometry, crs="EPSG:4326")
+    plant_data = gpd.sjoin(plant_data, regions, how="inner", predicate="within")
+    plant_data.rename(columns={"name": "bus"}, inplace=True)
+    # filter for countries in model scope
+    plant_data = plant_data[plant_data.Country.isin(snakemake.params.countries)]
+    # replace UK with GB in Country column
+    plant_data["Country"] = plant_data["Country"].replace("UK", "GB")
+    # assign industry grouping year
+    grouping_years = snakemake.params.existing_capacities["grouping_years_industry"]
+    plant_data.loc[:, "Year of last modernisation"] = plant_data[
+        "Year of last modernisation"
+    ].replace("x", np.nan)
+    plant_data["grouping_year"] = 0
+    valid_mask = plant_data["Year of last modernisation"].notna()
+    valid_years = plant_data.loc[valid_mask, "Year of last modernisation"]
+    indices = np.searchsorted(grouping_years, valid_years, side="right")
+    plant_data.loc[valid_years.index, "grouping_year"] = np.array(grouping_years)[
+        indices
+    ]
+
+    return plant_data, regions
+
+
+def add_existing_ammonia_plants(
+    n: pypsa.Network,
+) -> None:
+    """
+    Adds existing Haber-Bosch plants.
+    The plants are running on natural gas only since the retrofitting to hydrogen would be associated with costs. Plants are not expected to run at a minimal part load to avoid forcing the use of natural gas in planning horizons with climate targets.
+    Exhaust heat is not integrated since assuming that heat is integrated to make the current process more efficient.
+    """
+    logger.info("Adding existing ammonia plants.")
+
+    plant_data, regions = prepare_plant_data(
+        snakemake.input.regions_onshore,
+        snakemake.input.isi_database,
+    )
+
+    fh_ammonia = plant_data[plant_data.Product == "Ammonia"]
+
+    fh_ammonia = fh_ammonia.groupby(
+        ["bus", "Country", "grouping_year", "Product"], as_index=False
+    )["Production in tons (calibrated)"].sum()
+
+    fh_ammonia.index = (
+        fh_ammonia["bus"]
+        + " Haber-Bosch-SMR-"
+        + fh_ammonia["grouping_year"].astype(str)
+    )
+    # add dataset for Non EU27 countries
+    df = pd.read_csv(snakemake.input.ammonia, index_col=0)
+
+    geometry = gpd.points_from_xy(df.Longitude, df.Latitude)
+    gdf = gpd.GeoDataFrame(df, geometry=geometry, crs="EPSG:4326")
+
+    gdf = gpd.sjoin(gdf, regions, how="inner", predicate="within")
+
+    gdf.rename(columns={"name": "bus"}, inplace=True)
+    gdf["Country"] = gdf.bus.str[:2]
+    # filter for countries that are missing
+    gdf = gdf[
+        (~gdf.Country.isin(fh_ammonia.Country.unique()))
+        & (gdf.Country.isin(snakemake.params.countries))
+    ]
+    # following approach from build_industrial_distribution_key.py
+    for country in gdf.Country:
+        facilities = gdf.query("Country == @country")
+        production = facilities["Ammonia [kt/a]"]
+        # assume 50% of the minimum production for missing values
+        production = production.fillna(0.5 * facilities["Ammonia [kt/a]"].min())
+
+    # missing data
+    gdf.drop(gdf[gdf["Ammonia [kt/a]"].isna()].index, inplace=True)
+
+    # get average plant age:
+    avg_age = plant_data[plant_data.Product == "Ammonia"][
+        "Year of last modernisation"
+    ].mean()
+    gdf["grouping_year"] = min(
+        y
+        for y in snakemake.params.existing_capacities["grouping_years_industry"]
+        if y > avg_age
+    )
+    # match database
+    gdf.index = (
+        gdf["bus"] + " Haber-Bosch-SMR-" + gdf["grouping_year"].values.astype(str)
+    )
+    gdf.rename(
+        columns={"Ammonia [kt/a]": "Production in tons (calibrated)"}, inplace=True
+    )
+    gdf["Production in tons (calibrated)"] *= 1e3
+
+    ammonia_plants = pd.concat(
+        [
+            fh_ammonia,
+            gdf[["bus", "Country", "grouping_year", "Production in tons (calibrated)"]],
+        ]
+    )
+
+    # https://dechema.de/dechema_media/Downloads/Positionspapiere/Technology_study_Low_carbon_energy_and_feedstock_for_the_European_chemical_industry.pdf
+    # page 56: 1.83 t_CO2/t_NH3
+    ch4_per_nh3 = (
+        1.83 / costs.at["gas", "CO2 intensity"] / snakemake.params["MWh_NH3_per_tNH3"]
+    )
+    n.add(
+        "Link",
+        ammonia_plants.index,
+        bus0=[bus + " gas" for bus in ammonia_plants.bus]
+        if snakemake.params.sector["gas_network"]
+        else "EU gas",
+        bus1=[bus + " NH3" for bus in ammonia_plants.bus]
+        if snakemake.params.sector["ammonia"]
+        else "EU NH3",
+        bus2=ammonia_plants.bus,
+        bus3="co2 atmosphere",
+        p_nom=ammonia_plants["Production in tons (calibrated)"]
+        .mul(snakemake.params.MWh_NH3_per_tNH3)
+        .div(ch4_per_nh3)
+        .div(8760)
+        .values,
+        p_nom_extendable=False,
+        carrier="Haber-Bosch",
+        efficiency=1 / ch4_per_nh3,
+        efficiency1=-costs.at["Haber-Bosch", "electricity-input"] / ch4_per_nh3,
+        efficiency2=costs.at["gas", "CO2 intensity"],
+        capital_cost=costs.at["Haber-Bosch", "capital_cost"]
+        / costs.at["Haber-Bosch", "electricity-input"],
+        marginal_cost=costs.at["Haber-Bosch", "VOM"]
+        / costs.at["Haber-Bosch", "electricity-input"],
+        build_year=ammonia_plants["grouping_year"],
+        lifetime=costs.at["Haber-Bosch", "lifetime"],
+    )
+
+
 if __name__ == "__main__":
     if "snakemake" not in globals():
         from scripts._helpers import mock_snakemake
@@ -789,6 +945,10 @@ def add_heating_capacities_installed_before_baseyear(
     if options.get("cluster_heat_buses", False):
         cluster_heat_buses(n)
 
+    # add existing industry plants
+    if snakemake.params.sector["ammonia"]:
+        add_existing_ammonia_plants(n)
+
     n.meta = dict(snakemake.config, **dict(wildcards=dict(snakemake.wildcards)))
 
     sanitize_custom_columns(n)