implementing existing ammonia plants

Author: toniseibold

config/config.default.yaml

@@ -472,6 +472,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:
@@ -812,7 +813,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
@@ -968,7 +969,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

@@ -12,6 +12,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_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

@@ -10,6 +10,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"
@@ -27,6 +29,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

@@ -13,6 +13,7 @@
 import country_converter as coco
 import numpy as np
 import pandas as pd
+import geopandas as gpd
 import powerplantmatching as pm
 import pypsa
 import xarray as xr
@@ -716,6 +717,80 @@ def add_heating_capacities_installed_before_baseyear(
             )
 
 
+def prepare_plant_data():
+
+    # add existing industry
+    regions = gpd.read_file(snakemake.input.regions_onshore).set_index("name")
+    
+    isi_data = pd.read_excel(snakemake.input.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):
+
+    logger.info("Adding existing ammonia plants.")
+
+    plant_data, regions = prepare_plant_data()
+
+    ammonia_plants = plant_data[plant_data.Product=="Ammonia"]
+    
+    ammonia_plants = ammonia_plants.groupby(['bus', 'Country', 'grouping_year', 'Product'], as_index=False)['Production in tons (calibrated)'].sum()
+    
+    ammonia_plants.index = ammonia_plants['bus'] + " Haber-Bosch-" + ammonia_plants['grouping_year'].astype(str)
+
+    # add dataset
+    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.country.isin(ammonia_plants.Country.unique())]
+
+    n.add(
+        "Link",
+        ammonia_plants.index,
+        bus0=ammonia_plants.bus,
+        bus1=[bus + " NH3" for bus in ammonia_plants.bus] if snakemake.params.sector["ammonia"] else "EU NH3",
+        bus2=[bus + " gas" for bus in ammonia_plants.bus] if snakemake.params.sector["gas_network"] else "EU gas",
+        p_nom=ammonia_plants["Production in tons (calibrated)"].mul(snakemake.params.MWh_NH3_per_tNH3).div(costs.at["Haber-Bosch", "electricity-input"]).div(8760).values,
+        p_nom_extendable=False,
+        carrier="Haber-Bosch",
+        efficiency=1 / costs.at["Haber-Bosch", "electricity-input"],
+        efficiency2=-costs.at["Haber-Bosch", "hydrogen-input"]
+        / costs.at["Haber-Bosch", "electricity-input"],
+        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
@@ -797,6 +872,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.sector.ammonia:
+        add_existing_ammonia_plants(n)
+
     n.meta = dict(snakemake.config, **dict(wildcards=dict(snakemake.wildcards)))
 
     sanitize_custom_columns(n)