Merge branch 'master' into scenario-management

.gitignore

@@ -30,6 +30,8 @@ doc/_build
 config/config.yaml
 config/scenarios.yaml
 
+config.yaml
+config/config.yaml
 
 dconf
 /data/links_p_nom.csv

config/config.default.yaml

@@ -317,7 +317,7 @@ pypsa_eur:
 
 # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#energy
 energy:
-  energy_totals_year: 2011
+  energy_totals_year: 2013
   base_emissions_year: 1990
   eurostat_report_year: 2016
   emissions: CO2
@@ -504,6 +504,7 @@ sector:
   SMR_cc: true
   regional_methanol_demand: false
   regional_oil_demand: false
+  regional_coal_demand: false
   regional_co2_sequestration_potential:
     enable: false
     attribute: 'conservative estimate Mt'
@@ -541,8 +542,8 @@ sector:
       efficiency_static: 0.98
       efficiency_per_1000km: 0.977
     H2 pipeline:
-      efficiency_per_1000km: 1 # 0.979
-      compression_per_1000km: 0.019
+      efficiency_per_1000km: 1 # 0.982
+      compression_per_1000km: 0.018
     gas pipeline:
       efficiency_per_1000km: 1 #0.977
       compression_per_1000km: 0.01
@@ -639,6 +640,15 @@ industry:
     2040: 0.12
     2045: 0.16
     2050: 0.20
+  sector_ratios_fraction_future:
+    2020: 0.0
+    2025: 0.1
+    2030: 0.3
+    2035: 0.5
+    2040: 0.7
+    2045: 0.9
+    2050: 1.0
+  basic_chemicals_without_NH3_production_today: 69. #Mt/a, = 86 Mtethylene-equiv - 17 MtNH3
   HVC_production_today: 52.
   MWh_elec_per_tHVC_mechanical_recycling: 0.547
   MWh_elec_per_tHVC_chemical_recycling: 6.9

data/switzerland-new_format-all_years.csv

@@ -0,0 +1,25 @@
+country,item,2010,2011,2012,2013,2014,2015
+CH,total residential,268.2,223.4,243.4,261.3,214.2,229.1
+CH,total residential space,192.2,149.0,168.1,185.5,139.7,154.4
+CH,total residential water,32.2,31.6,31.9,32.2,31.7,31.9
+CH,total residential cooking,9.3,9.3,9.3,9.4,9.5,9.6
+CH,electricity residential,67.9,63.7,65.7,67.6,63.0,64.4
+CH,electricity residential space,15.9,12.8,14.3,15.8,12.3,13.5
+CH,electricity residential water,8.8,8.5,8.5,8.6,8.5,8.6
+CH,electricity residential cooking,4.9,4.9,4.9,4.9,5.0,5.0
+CH,total services,145.9,127.4,136.7,144.0,124.5,132.5
+CH,total services space,80.0,62.2,70.8,77.4,58.3,64.3
+CH,total services water,10.1,10.0,10.1,10.1,10.0,10.0
+CH,total services cooking,2.5,2.4,2.3,2.3,2.4,2.3
+CH,electricity services,60.5,59.2,60.3,61.4,60.3,62.6
+CH,electricity services space,4.0,3.2,3.8,4.2,3.3,3.6
+CH,electricity services water,0.7,0.7,0.7,0.7,0.7,0.7
+CH,electricity services cooking,2.5,2.4,2.3,2.3,2.4,2.3
+CH,total rail,11.5,11.1,11.2,11.4,11.1,11.4
+CH,total road,199.4,200.4,200.4,201.2,202.0,203.1
+CH,electricity road,0.,0.,0.,0.,0.,0.
+CH,electricity rail,11.5,11.1,11.2,11.4,11.1,11.4
+CH,total domestic aviation,3.3,3.2,3.4,3.4,3.5,3.5
+CH,total international aviation,58.0,62.0,63.5,64.2,64.5,66.8
+CH,total domestic navigation,1.6,1.6,1.6,1.6,1.6,1.6
+CH,total international navigation,0.,0.,0.,0.,0.,0.

doc/configtables/industry.csv

@@ -17,6 +17,8 @@ HVC_primary_fraction,--,float,The fraction of high value chemicals (HVC) produce
 HVC_mechanical_recycling _fraction,--,float,The fraction of high value chemicals (HVC) produced using mechanical recycling
 HVC_chemical_recycling _fraction,--,float,The fraction of high value chemicals (HVC) produced using chemical recycling
 ,,,
+sector_ratios_fraction_future,--,Dictionary with planning horizons as keys.,The fraction of total progress in fuel and process switching achieved in the industry sector.
+basic_chemicals_without_NH3_production_today,Mt/a,float,"The amount of basic chemicals produced without ammonia (= 86 Mtethylene-equiv - 17 MtNH3)."
 HVC_production_today,MtHVC/a,float,"The amount of high value chemicals (HVC) produced. This includes ethylene, propylene and BTX. From `DECHEMA (2017) <https://dechema.de/dechema_media/Downloads/Positionspapiere/Technology_study_Low_carbon_energy_and_feedstock_for_the_European_chemical_industry-p-20002750.pdf>`_, Figure 16, page 107"
 Mwh_elec_per_tHVC _mechanical_recycling,MWh/tHVC,float,"The energy amount of electricity needed to produce a ton of high value chemical (HVC) using mechanical recycling. From SI of `Meys et al (2020) <https://doi.org/10.1016/j.resconrec.2020.105010>`_, Table S5, for HDPE, PP, PS, PET. LDPE would be 0.756."
 Mwh_elec_per_tHVC _chemical_recycling,MWh/tHVC,float,"The energy amount of electricity needed to produce a ton of high value chemical (HVC) using chemical recycling. The default value is based on pyrolysis and electric steam cracking. From `Material Economics (2019) <https://materialeconomics.com/latest-updates/industrial-transformation-2050>`_, page 125"

doc/release_notes.rst

@@ -17,6 +17,29 @@ Upcoming Release
   names to run under ``run: name:`` in the configuration file. The latter must
   be a subset of toplevel keys in the scenario file.
 
+* Improved representation of industry transition pathways. A new script was
+  added to interpolate industry sector ratios from today's status quo to future
+  systems (i.e. specific emissions and demands for energy and feedstocks). For
+  each country we gradually switch industry processes from today's specific
+  energy carrier usage per ton material output to the best-in-class energy
+  consumption of tomorrow. This is done on a per-country basis. The ratio of
+  today to tomorrow's energy consumption is set with the ``industry:
+  sector_ratios_fraction_future:`` parameter.
+
+* Fix plotting of retrofitted hydrogen pipelines with pathway optimisation.
+
+* Bugfix: Correct units of subtracted chlorine and methanol demand in
+  :mod:`build_industry_sector_ratios`.
+
+* Include all countries in ammonia production resource. This is so that the full
+  EU28 ammonia demand can be correctly subtracted in the rule
+  :mod:`build_industry_sector_ratios`.
+
+* Regions are assigned to all buses with unique coordinates in the network with
+  a preference given to substations. Previously, only substations had assigned
+  regions, but this could lead to issues when a high spatial resolution was
+  applied.
+
 * The default configuration ``config/config.default.yaml`` is now automatically
   used as a base configuration file and no longer copied to
   ``config/config.yaml`` on first use. The file ``config/config.yaml`` should be
@@ -24,6 +47,10 @@ Upcoming Release
 
 * Merged two OPSD time series data versions into such that the option ``load:
   power_statistics:`` becomes superfluous and was hence removed.
+* Bugfix: The industry coal emissions for industry were not properly tracked.
+
+* Allow industrial coal demand to be regional so its emissions can be included
+  in regional emission limits.
 
 * Add new default to overdimension heating in individual buildings. This allows
   them to cover heat demand peaks e.g. 10% higher than those in the data. The
@@ -120,6 +147,8 @@ Upcoming Release
   workflows with foresight "myopic" and still needs to be added foresight option
   "perfect".
 
+* Switched the energy totals year from 2011 to 2013 to comply with the assumed default weather year.
+
 
 PyPSA-Eur 0.9.0 (5th January 2024)
 ==================================

rules/build_sector.smk

@@ -279,7 +279,7 @@ rule build_energy_totals:
     input:
         nuts3_shapes=resources("nuts3_shapes.geojson"),
         co2="data/bundle-sector/eea/UNFCCC_v23.csv",
-        swiss="data/bundle-sector/switzerland-sfoe/switzerland-new_format.csv",
+        swiss="data/switzerland-new_format-all_years.csv",
         idees="data/bundle-sector/jrc-idees-2015",
         district_heat_share="data/district_heat_share.csv",
         eurostat=input_eurostat,
@@ -406,8 +406,6 @@ rule build_salt_cavern_potentials:
 
 
 rule build_ammonia_production:
-    params:
-        countries=config_provider("countries"),
     input:
         usgs="data/bundle-sector/myb1-2017-nitro.xls",
     output:
@@ -447,6 +445,31 @@ rule build_industry_sector_ratios:
         "../scripts/build_industry_sector_ratios.py"
 
 
+rule build_industry_sector_ratios_intermediate:
+    params:
+        industry=config["industry"],
+    input:
+        industry_sector_ratios=RESOURCES + "industry_sector_ratios.csv",
+        industrial_energy_demand_per_country_today=RESOURCES
+        + "industrial_energy_demand_per_country_today.csv",
+        industrial_production_per_country=RESOURCES
+        + "industrial_production_per_country.csv",
+    output:
+        industry_sector_ratios=RESOURCES
+        + "industry_sector_ratios_{planning_horizons}.csv",
+    threads: 1
+    resources:
+        mem_mb=1000,
+    log:
+        LOGS + "build_industry_sector_ratios_{planning_horizons}.log",
+    benchmark:
+        BENCHMARKS + "build_industry_sector_ratios_{planning_horizons}"
+    conda:
+        "../envs/environment.yaml"
+    script:
+        "../scripts/build_industry_sector_ratios_intermediate.py"
+
+
 rule build_industrial_production_per_country:
     params:
         industry=config_provider("industry"),
@@ -560,7 +583,7 @@ rule build_industrial_production_per_node:
 
 rule build_industrial_energy_demand_per_node:
     input:
-        industry_sector_ratios=resources("industry_sector_ratios.csv"),
+        industry_sector_ratios=resources("industry_sector_ratios_{planning_horizons}.csv"),
         industrial_production_per_node=resources(
             "industrial_production_elec_s{simpl}_{clusters}_{planning_horizons}.csv"
         ),
@@ -596,7 +619,6 @@ rule build_industrial_energy_demand_per_country_today:
         industry=config_provider("industry"),
     input:
         jrc="data/bundle-sector/jrc-idees-2015",
-        ammonia_production=resources("ammonia_production.csv"),
         industrial_production_per_country=resources(
             "industrial_production_per_country.csv"
         ),

scripts/base_network.py

@@ -559,6 +559,7 @@ def _set_countries_and_substations(n, config, country_shapes, offshore_shapes):
     for b, df in product(("bus0", "bus1"), (n.lines, n.links)):
         has_connections_b |= ~df.groupby(b).under_construction.min()
 
+    buses["onshore_bus"] = onshore_b
     buses["substation_lv"] = (
         lv_b & onshore_b & (~buses["under_construction"]) & has_connections_b
     )

scripts/build_ammonia_production.py

@@ -28,13 +28,14 @@ if __name__ == "__main__":
         header=0,
         index_col=0,
         skipfooter=19,
+        na_values=["--"],
     )
 
     ammonia.index = cc.convert(ammonia.index, to="iso2")
 
     years = [str(i) for i in range(2013, 2018)]
-    countries = ammonia.index.intersection(snakemake.params.countries)
-    ammonia = ammonia.loc[countries, years].astype(float)
+
+    ammonia = ammonia[years]
 
     # convert from ktonN to ktonNH3
     ammonia *= 17 / 14

scripts/build_bus_regions.py

@@ -136,7 +136,13 @@ if __name__ == "__main__":
         c_b = n.buses.country == country
 
         onshore_shape = country_shapes[country]
-        onshore_locs = n.buses.loc[c_b & n.buses.substation_lv, ["x", "y"]]
+        onshore_locs = (
+            n.buses.loc[c_b & n.buses.onshore_bus]
+            .sort_values(
+                by="substation_lv", ascending=False
+            )  # preference for substations
+            .drop_duplicates(subset=["x", "y"], keep="first")[["x", "y"]]
+        )
         onshore_regions.append(
             gpd.GeoDataFrame(
                 {

scripts/build_industrial_energy_demand_per_country_today.py

@@ -74,7 +74,7 @@ def industrial_energy_demand_per_country(country, year, jrc_dir):
     def get_subsector_data(sheet):
         df = df_dict[sheet][year].groupby(fuels).sum()
 
-        df["ammonia"] = 0.0
+        df["hydrogen"] = 0.0
 
         df["other"] = df["all"] - df.loc[df.index != "all"].sum()
 
@@ -95,51 +95,50 @@ def industrial_energy_demand_per_country(country, year, jrc_dir):
     return df
 
 
-def add_ammonia_energy_demand(demand):
-    # MtNH3/a
-    fn = snakemake.input.ammonia_production
-    ammonia = pd.read_csv(fn, index_col=0)[str(year)] / 1e3
+def separate_basic_chemicals(demand, production):
 
-    def get_ammonia_by_fuel(x):
-        fuels = {
-            "gas": params["MWh_CH4_per_tNH3_SMR"],
-            "electricity": params["MWh_elec_per_tNH3_SMR"],
+    ammonia = pd.DataFrame(
+        {
+            "hydrogen": production["Ammonia"] * params["MWh_H2_per_tNH3_electrolysis"],
+            "electricity": production["Ammonia"]
+            * params["MWh_elec_per_tNH3_electrolysis"],
         }
-
-        return pd.Series({k: x * v for k, v in fuels.items()})
-
-    ammonia_by_fuel = ammonia.apply(get_ammonia_by_fuel).T
-    ammonia_by_fuel = ammonia_by_fuel.unstack().reindex(
-        index=demand.index, fill_value=0.0
-    )
-
-    ammonia = pd.DataFrame({"ammonia": ammonia * params["MWh_NH3_per_tNH3"]}).T
+    ).T
+    chlorine = pd.DataFrame(
+        {
+            "hydrogen": production["Chlorine"] * params["MWh_H2_per_tCl"],
+            "electricity": production["Chlorine"] * params["MWh_elec_per_tCl"],
+        }
+    ).T
+    methanol = pd.DataFrame(
+        {
+            "gas": production["Methanol"] * params["MWh_CH4_per_tMeOH"],
+            "electricity": production["Methanol"] * params["MWh_elec_per_tMeOH"],
+        }
+    ).T
 
     demand["Ammonia"] = ammonia.unstack().reindex(index=demand.index, fill_value=0.0)
+    demand["Chlorine"] = chlorine.unstack().reindex(index=demand.index, fill_value=0.0)
+    demand["Methanol"] = methanol.unstack().reindex(index=demand.index, fill_value=0.0)
 
-    demand["Basic chemicals (without ammonia)"] = (
-        demand["Basic chemicals"] - ammonia_by_fuel
+    demand["HVC"] = (
+        demand["Basic chemicals"]
+        - demand["Ammonia"]
+        - demand["Methanol"]
+        - demand["Chlorine"]
     )
 
-    demand["Basic chemicals (without ammonia)"].clip(lower=0, inplace=True)
-
     demand.drop(columns="Basic chemicals", inplace=True)
 
+    demand["HVC"].clip(lower=0, inplace=True)
+
     return demand
 
 
-def add_non_eu28_industrial_energy_demand(countries, demand):
+def add_non_eu28_industrial_energy_demand(countries, demand, production):
     non_eu28 = countries.difference(eu28)
     if non_eu28.empty:
         return demand
-    # output in MtMaterial/a
-    fn = snakemake.input.industrial_production_per_country
-    production = pd.read_csv(fn, index_col=0) / 1e3
-
-    # recombine HVC, Chlorine and Methanol to Basic chemicals (without ammonia)
-    chemicals = ["HVC", "Chlorine", "Methanol"]
-    production["Basic chemicals (without ammonia)"] = production[chemicals].sum(axis=1)
-    production.drop(columns=chemicals, inplace=True)
 
     eu28_production = production.loc[countries.intersection(eu28)].sum()
     eu28_energy = demand.groupby(level=1).sum()
@@ -184,9 +183,15 @@ if __name__ == "__main__":
 
     demand = industrial_energy_demand(countries.intersection(eu28), year)
 
-    demand = add_ammonia_energy_demand(demand)
+    # output in MtMaterial/a
+    production = (
+        pd.read_csv(snakemake.input.industrial_production_per_country, index_col=0)
+        / 1e3
+    )
 
-    demand = add_non_eu28_industrial_energy_demand(countries, demand)
+    demand = separate_basic_chemicals(demand, production)
+
+    demand = add_non_eu28_industrial_energy_demand(countries, demand, production)
 
     # for format compatibility
     demand = demand.stack(dropna=False).unstack(level=[0, 2])

scripts/build_industrial_energy_demand_per_node.py

@@ -21,23 +21,31 @@ if __name__ == "__main__":
         )
     set_scenario_config(snakemake)
 
-    # import EU ratios df as csv
+    # import ratios
     fn = snakemake.input.industry_sector_ratios
-    industry_sector_ratios = pd.read_csv(fn, index_col=0)
+    sector_ratios = pd.read_csv(fn, header=[0, 1], index_col=0)
 
-    # material demand per node and industry (kton/a)
+    # material demand per node and industry (Mton/a)
     fn = snakemake.input.industrial_production_per_node
-    nodal_production = pd.read_csv(fn, index_col=0)
+    nodal_production = pd.read_csv(fn, index_col=0) / 1e3
 
     # energy demand today to get current electricity
     fn = snakemake.input.industrial_energy_demand_per_node_today
     nodal_today = pd.read_csv(fn, index_col=0)
 
-    # final energy consumption per node and industry (TWh/a)
-    nodal_df = nodal_production.dot(industry_sector_ratios.T)
+    nodal_sector_ratios = pd.concat(
+        {node: sector_ratios[node[:2]] for node in nodal_production.index}, axis=1
+    )
 
-    # convert GWh to TWh and ktCO2 to MtCO2
-    nodal_df *= 0.001
+    nodal_production_stacked = nodal_production.stack()
+    nodal_production_stacked.index.names = [None, None]
+
+    # final energy consumption per node and industry (TWh/a)
+    nodal_df = (
+        (nodal_sector_ratios.multiply(nodal_production_stacked))
+        .T.groupby(level=0)
+        .sum()
+    )
 
     rename_sectors = {
         "elec": "electricity",

scripts/build_industrial_production_per_country.py

@@ -261,7 +261,11 @@ def separate_basic_chemicals(demand, year):
     demand["Basic chemicals"].clip(lower=0.0, inplace=True)
 
     # assume HVC, methanol, chlorine production proportional to non-ammonia basic chemicals
-    distribution_key = demand["Basic chemicals"] / demand["Basic chemicals"].sum()
+    distribution_key = (
+        demand["Basic chemicals"]
+        / params["basic_chemicals_without_NH3_production_today"]
+        / 1e3
+    )
     demand["HVC"] = params["HVC_production_today"] * 1e3 * distribution_key
     demand["Chlorine"] = params["chlorine_production_today"] * 1e3 * distribution_key
     demand["Methanol"] = params["methanol_production_today"] * 1e3 * distribution_key

scripts/build_industry_sector_ratios.py

@@ -408,15 +408,15 @@ def chemicals_industry():
     df.loc["methane", sector] -= ammonia_total * params["MWh_CH4_per_tNH3_SMR"]
     df.loc["elec", sector] -= ammonia_total * params["MWh_elec_per_tNH3_SMR"]
 
-    # subtract chlorine demand
+    # subtract chlorine demand (in MtCl/a)
     chlorine_total = params["chlorine_production_today"]
-    df.loc["hydrogen", sector] -= chlorine_total * params["MWh_H2_per_tCl"]
-    df.loc["elec", sector] -= chlorine_total * params["MWh_elec_per_tCl"]
+    df.loc["hydrogen", sector] -= chlorine_total * params["MWh_H2_per_tCl"] * 1e3
+    df.loc["elec", sector] -= chlorine_total * params["MWh_elec_per_tCl"] * 1e3
 
-    # subtract methanol demand
+    # subtract methanol demand (in MtMeOH/a)
     methanol_total = params["methanol_production_today"]
-    df.loc["methane", sector] -= methanol_total * params["MWh_CH4_per_tMeOH"]
-    df.loc["elec", sector] -= methanol_total * params["MWh_elec_per_tMeOH"]
+    df.loc["methane", sector] -= methanol_total * params["MWh_CH4_per_tMeOH"] * 1e3
+    df.loc["elec", sector] -= methanol_total * params["MWh_elec_per_tMeOH"] * 1e3
 
     # MWh/t material
     df.loc[sources, sector] = df.loc[sources, sector] / s_out

scripts/build_industry_sector_ratios_intermediate.py

@@ -0,0 +1,81 @@
+# -*- coding: utf-8 -*-
+# SPDX-FileCopyrightText: : 2020-2024 The PyPSA-Eur Authors
+#
+# SPDX-License-Identifier: MIT
+"""
+Build specific energy consumption by carrier and industries and by country,
+that interpolates between the current average energy consumption (from
+2015-2020) and the ideal future best-in-class consumption.
+"""
+
+import pandas as pd
+from prepare_sector_network import get
+
+
+def build_industry_sector_ratios_intermediate():
+
+    # in TWh/a
+    demand = pd.read_csv(
+        snakemake.input.industrial_energy_demand_per_country_today,
+        header=[0, 1],
+        index_col=0,
+    )
+
+    # in Mt/a
+    production = (
+        pd.read_csv(snakemake.input.industrial_production_per_country, index_col=0)
+        / 1e3
+    ).stack()
+    production.index.names = [None, None]
+
+    # in MWh/t
+    future_sector_ratios = pd.read_csv(
+        snakemake.input.industry_sector_ratios, index_col=0
+    )
+
+    today_sector_ratios = demand.div(production, axis=1)
+
+    today_sector_ratios.dropna(how="all", axis=1, inplace=True)
+
+    rename = {
+        "waste": "biomass",
+        "electricity": "elec",
+        "solid": "coke",
+        "gas": "methane",
+        "other": "biomass",
+        "liquid": "naphtha",
+    }
+    today_sector_ratios = today_sector_ratios.rename(rename).groupby(level=0).sum()
+
+    fraction_future = get(params["sector_ratios_fraction_future"], year)
+
+    intermediate_sector_ratios = {}
+    for ct, group in today_sector_ratios.T.groupby(level=0):
+        today_sector_ratios_ct = (
+            group.droplevel(0)
+            .T.reindex_like(future_sector_ratios)
+            .fillna(future_sector_ratios)
+        )
+        intermediate_sector_ratios[ct] = (
+            today_sector_ratios_ct * (1 - fraction_future)
+            + future_sector_ratios * fraction_future
+        )
+    intermediate_sector_ratios = pd.concat(intermediate_sector_ratios, axis=1)
+
+    intermediate_sector_ratios.to_csv(snakemake.output.industry_sector_ratios)
+
+
+if __name__ == "__main__":
+    if "snakemake" not in globals():
+        from _helpers import mock_snakemake
+
+        snakemake = mock_snakemake(
+            "build_industry_sector_ratios_intermediate",
+            planning_horizons="2030",
+        )
+
+    year = int(snakemake.wildcards.planning_horizons[-4:])
+
+    params = snakemake.params.industry
+
+    build_industry_sector_ratios_intermediate()

scripts/plot_hydrogen_network.py

@@ -24,6 +24,7 @@ def group_pipes(df, drop_direction=False):
     """
     Group pipes which connect same buses and return overall capacity.
     """
+    df = df.copy()
     if drop_direction:
         positive_order = df.bus0 < df.bus1
         df_p = df[positive_order]
@@ -32,12 +33,13 @@ def group_pipes(df, drop_direction=False):
         df = pd.concat([df_p, df_n])
 
     # there are pipes for each investment period rename to AC buses name for plotting
+    df["index_orig"] = df.index
     df.index = df.apply(
         lambda x: f"H2 pipeline {x.bus0.replace(' H2', '')} -> {x.bus1.replace(' H2', '')}",
         axis=1,
     )
     return df.groupby(level=0).agg(
-        {"p_nom_opt": "sum", "bus0": "first", "bus1": "first"}
+        {"p_nom_opt": "sum", "bus0": "first", "bus1": "first", "index_orig": "first"}
     )
 
 
@@ -95,17 +97,18 @@ def plot_h2_map(n, regions):
             )
 
     if not h2_retro.empty:
-        positive_order = h2_retro.bus0 < h2_retro.bus1
-        h2_retro_p = h2_retro[positive_order]
-        swap_buses = {"bus0": "bus1", "bus1": "bus0"}
-        h2_retro_n = h2_retro[~positive_order].rename(columns=swap_buses)
-        h2_retro = pd.concat([h2_retro_p, h2_retro_n])
+        if snakemake.params.foresight != "myopic":
+            positive_order = h2_retro.bus0 < h2_retro.bus1
+            h2_retro_p = h2_retro[positive_order]
+            swap_buses = {"bus0": "bus1", "bus1": "bus0"}
+            h2_retro_n = h2_retro[~positive_order].rename(columns=swap_buses)
+            h2_retro = pd.concat([h2_retro_p, h2_retro_n])
 
-        h2_retro["index_orig"] = h2_retro.index
-        h2_retro.index = h2_retro.apply(
-            lambda x: f"H2 pipeline {x.bus0.replace(' H2', '')} -> {x.bus1.replace(' H2', '')}",
-            axis=1,
-        )
+            h2_retro["index_orig"] = h2_retro.index
+            h2_retro.index = h2_retro.apply(
+                lambda x: f"H2 pipeline {x.bus0.replace(' H2', '')} -> {x.bus1.replace(' H2', '')}",
+                axis=1,
+            )
 
         retro_w_new_i = h2_retro.index.intersection(h2_new.index)
         h2_retro_w_new = h2_retro.loc[retro_w_new_i]

scripts/prepare_sector_network.py

@@ -176,6 +176,13 @@ def define_spatial(nodes, options):
     spatial.coal.nodes = ["EU coal"]
     spatial.coal.locations = ["EU"]
 
+    if options["regional_coal_demand"]:
+        spatial.coal.demand_locations = nodes
+        spatial.coal.industry = nodes + " coal for industry"
+    else:
+        spatial.coal.demand_locations = ["EU"]
+        spatial.coal.industry = ["EU coal for industry"]
+
     # lignite
     spatial.lignite = SimpleNamespace()
     spatial.lignite.nodes = ["EU lignite"]
@@ -3052,19 +3059,40 @@ def add_industry(n, costs):
 
         mwh_coal_per_mwh_coke = 1.366  # from eurostat energy balance
         p_set = (
-            industrial_demand["coal"].sum()
-            + mwh_coal_per_mwh_coke * industrial_demand["coke"].sum()
+            industrial_demand["coal"]
+            + mwh_coal_per_mwh_coke * industrial_demand["coke"]
         ) / nhours
 
+        if not options["regional_coal_demand"]:
+            p_set = p_set.sum()
+
+        n.madd(
+            "Bus",
+            spatial.coal.industry,
+            location=spatial.coal.demand_locations,
+            carrier="coal for industry",
+            unit="MWh_LHV",
+        )
+
         n.madd(
             "Load",
-            spatial.coal.nodes,
-            suffix=" for industry",
-            bus=spatial.coal.nodes,
+            spatial.coal.industry,
+            bus=spatial.coal.industry,
             carrier="coal for industry",
             p_set=p_set,
         )
 
+        n.madd(
+            "Link",
+            spatial.coal.industry,
+            bus0=spatial.coal.nodes,
+            bus1=spatial.coal.industry,
+            bus2="co2 atmosphere",
+            carrier="coal for industry",
+            p_nom_extendable=True,
+            efficiency2=costs.at["coal", "CO2 intensity"],
+        )
+
 
 def add_waste_heat(n):
     # TODO options?

scripts/simplify_network.py

@@ -612,6 +612,7 @@ if __name__ == "__main__":
         "symbol",
         "tags",
         "under_construction",
+        "onshore_bus",
         "substation_lv",
         "substation_off",
         "geometry",