refining existing Haber-Bosch plants + Documentation

Author: toniseibold

scripts/add_existing_baseyear.py

@@ -717,12 +717,25 @@ def add_heating_capacities_installed_before_baseyear(
             )
 
 
-def prepare_plant_data():
+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(snakemake.input.regions_onshore).set_index("name")
+    regions = gpd.read_file(regions_fn).set_index("name")
 
-    isi_data = pd.read_excel(snakemake.input.isi_database, sheet_name="Database", index_col=1)
+    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")
@@ -744,19 +757,27 @@ def prepare_plant_data():
     return plant_data, regions
 
 
-def add_existing_ammonia_plants(n):
-
+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()
+    plant_data, regions = prepare_plant_data(
+        snakemake.input.regions_onshore,
+        snakemake.input.isi_database,
+    )
 
-    ammonia_plants = plant_data[plant_data.Product=="Ammonia"]
+    fh_ammonia = plant_data[plant_data.Product=="Ammonia"]
 
-    ammonia_plants = ammonia_plants.groupby(['bus', 'Country', 'grouping_year', 'Product'], as_index=False)['Production in tons (calibrated)'].sum()
+    fh_ammonia = fh_ammonia.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
+    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)
@@ -765,22 +786,45 @@ def add_existing_ammonia_plants(n):
     gdf = gpd.sjoin(gdf, regions, how="inner", predicate="within")
 
     gdf.rename(columns={"name": "bus"}, inplace=True)
-    gdf["country"] = gdf.bus.str[:2]
+    gdf["Country"] = gdf.bus.str[:2]    
     # filter for countries that are missing
-    gdf[~gdf.country.isin(ammonia_plants.Country.unique())]
-
+    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=ammonia_plants.bus,
+        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=[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,
+        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 / costs.at["Haber-Bosch", "electricity-input"],
-        efficiency2=-costs.at["Haber-Bosch", "hydrogen-input"]
-        / costs.at["Haber-Bosch", "electricity-input"],
+        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"]
@@ -873,7 +917,7 @@ def add_existing_ammonia_plants(n):
         cluster_heat_buses(n)
 
     # add existing industry plants
-    if snakemake.sector.ammonia:
+    if snakemake.params.sector["ammonia"]:
         add_existing_ammonia_plants(n)
 
     n.meta = dict(snakemake.config, **dict(wildcards=dict(snakemake.wildcards)))