Improve co2 limit

CHANGELOG.md

@@ -1,4 +1,5 @@
 # Changelog
+- Minor improvements to the DE CO2 constraint
 - Bugfix: Enforce stricter H2 derivative import limit to avoid that exports of one type of derivative compensate for imports of another
 - Added an option to source mobility demand from UBA MWMS (Projektionsbericht 2025) for the years 2025-2035
 - Renamed functions and script for exogenous mobility demand

config/config.de.yaml

@@ -4,7 +4,7 @@
 
 # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#run
 run:
-  prefix: 20250901_h2_deriv_fix
+  prefix: 20251013_improve_co2_limit
   name:
   # - ExPol
   - KN2045_Mix

scripts/pypsa-de/additional_functionality.py

@@ -373,6 +373,7 @@ def add_national_co2_budgets(n, snakemake, national_co2_budgets, investment_year
     nyears = nhours / 8760
 
     sectors = determine_emission_sectors(n.config["sector"])
+    energy_totals = pd.read_csv(snakemake.input.energy_totals, index_col=[0, 1])
 
     # convert MtCO2 to tCO2
     co2_totals = 1e6 * pd.read_csv(snakemake.input.co2_totals_name, index_col=0)
@@ -397,8 +398,8 @@ def add_national_co2_budgets(n, snakemake, national_co2_budgets, investment_year
             links = n.links.index[
                 (n.links.index.str[:2] == ct)
                 & (n.links[f"bus{port}"] == "co2 atmosphere")
-                & (
-                    n.links.carrier != "kerosene for aviation"
+                & ~n.links.carrier.str.contains(
+                    "shipping|aviation"
                 )  # first exclude aviation to multiply it with a domestic factor later
             ]
 
@@ -422,27 +423,68 @@ def add_national_co2_budgets(n, snakemake, national_co2_budgets, investment_year
             )
 
         # Aviation demand
-        energy_totals = pd.read_csv(snakemake.input.energy_totals, index_col=[0, 1])
         domestic_aviation = energy_totals.loc[
             (ct, snakemake.params.energy_year), "total domestic aviation"
         ]
         international_aviation = energy_totals.loc[
             (ct, snakemake.params.energy_year), "total international aviation"
         ]
-        domestic_factor = domestic_aviation / (
+        domestic_aviation_factor = domestic_aviation / (
             domestic_aviation + international_aviation
         )
         aviation_links = n.links[
             (n.links.index.str[:2] == ct) & (n.links.carrier == "kerosene for aviation")
         ]
-        lhs.append
-        (
-            n.model["Link-p"].loc[:, aviation_links.index]
-            * aviation_links.efficiency2
-            * n.snapshot_weightings.generators
-        ).sum() * domestic_factor
+        lhs.append(
+            (
+                n.model["Link-p"].loc[:, aviation_links.index]
+                * aviation_links.efficiency2
+                * n.snapshot_weightings.generators
+            ).sum()
+            * domestic_aviation_factor
+        )
+        logger.info(
+            f"Adding domestic aviation emissions for {ct} with a factor of {domestic_aviation_factor}"
+        )
+
+        # Shipping oil
+        domestic_navigation = energy_totals.loc[
+            (ct, snakemake.params.energy_year), "total domestic navigation"
+        ]
+        international_navigation = energy_totals.loc[
+            (ct, snakemake.params.energy_year), "total international navigation"
+        ]
+        domestic_navigation_factor = domestic_navigation / (
+            domestic_navigation + international_navigation
+        )
+        shipping_links = n.links[
+            (n.links.index.str[:2] == ct) & (n.links.carrier == "shipping oil")
+        ]
+        lhs.append(
+            (
+                n.model["Link-p"].loc[:, shipping_links.index]
+                * shipping_links.efficiency2
+                * n.snapshot_weightings.generators
+            ).sum()
+            * domestic_navigation_factor
+        )
+
+        # Shipping methanol
+        shipping_meoh_links = n.links[
+            (n.links.index.str[:2] == ct) & (n.links.carrier == "shipping methanol")
+        ]
+        if not shipping_meoh_links.empty:  # no shipping methanol in 2025
+            lhs.append(
+                (
+                    n.model["Link-p"].loc[:, shipping_meoh_links.index]
+                    * shipping_meoh_links.efficiency2
+                    * n.snapshot_weightings.generators
+                ).sum()
+                * domestic_navigation_factor
+            )
+
         logger.info(
-            f"Adding domestic aviation emissions for {ct} with a factor of {domestic_factor}"
+            f"Adding domestic shipping emissions for {ct} with a factor of {domestic_navigation_factor}"
         )
 
         # Adding Efuel imports and exports to constraint

scripts/pypsa-de/build_scenarios.py

@@ -117,14 +117,16 @@ def get_co2_budget(df, source):
 
     nonco2 = ghg - co2
 
+    # Hard-code values for 2020 and 2025 from UBA reports / projections
+    # Table 1, https://reportnet.europa.eu/public/dataflow/1478, GHG - CO2
+    nonco2[2020] = 733.7 - 647.9
+    nonco2[2025] = 628.8 - 554.6
+
     ## PyPSA disregards nonco2 GHG emissions, but includes bunkers
 
     targets_pypsa = targets_co2 - nonco2
 
     target_fractions_pypsa = targets_pypsa.loc[targets_co2.index] / baseline_pypsa
-    target_fractions_pypsa[2020] = (
-        0.671  # Hard-coded based on REMIND data from ariadne2-internal DB
-    )
 
     return target_fractions_pypsa.round(3)
 

scripts/pypsa-de/export_ariadne_variables.py

@@ -2573,6 +2573,167 @@ def get_final_energy(
 
 
 def get_emissions(n, region, _energy_totals, industry_demand):
+    def get_constraint_emissions(n, ct):
+        lhs = []
+
+        for port in [col[3:] for col in n.links if col.startswith("bus")]:
+            links = n.links.index[
+                (n.links.index.str[:2] == ct)
+                & (n.links[f"bus{port}"] == "co2 atmosphere")
+                & ~n.links.carrier.str.contains(
+                    "shipping|aviation"
+                )  # first exclude aviation to multiply it with a domestic factor later
+            ]
+
+            if port == "0":
+                efficiency = -1.0
+            elif port == "1":
+                efficiency = n.links.loc[links, "efficiency"]
+            else:
+                efficiency = n.links.loc[links, f"efficiency{port}"]
+
+            variables = (
+                n.links_t.p0.loc[:, links]
+                .mul(efficiency)
+                .mul(n.snapshot_weightings.generators, axis=0)
+                .sum()
+            )
+            if not variables.empty:
+                lhs.append(variables)
+
+        # Aviation demand
+        energy_totals = pd.read_csv(snakemake.input.energy_totals, index_col=[0, 1])
+        domestic_aviation = energy_totals.loc[
+            (ct, snakemake.params.energy_totals_year), "total domestic aviation"
+        ]
+        international_aviation = energy_totals.loc[
+            (ct, snakemake.params.energy_totals_year), "total international aviation"
+        ]
+        domestic_aviation_factor = domestic_aviation / (
+            domestic_aviation + international_aviation
+        )
+        aviation_links = n.links[
+            (n.links.index.str[:2] == ct) & (n.links.carrier == "kerosene for aviation")
+        ]
+        lhs.append(
+            (
+                n.links_t.p0.loc[:, aviation_links.index]
+                .mul(aviation_links.efficiency2)
+                .mul(n.snapshot_weightings.generators, axis=0)
+            ).sum()
+            * domestic_aviation_factor
+        )
+
+        # Shipping oil
+        domestic_navigation = energy_totals.loc[
+            (ct, snakemake.params.energy_totals_year), "total domestic navigation"
+        ]
+        international_navigation = energy_totals.loc[
+            (ct, snakemake.params.energy_totals_year), "total international navigation"
+        ]
+        domestic_navigation_factor = domestic_navigation / (
+            domestic_navigation + international_navigation
+        )
+        shipping_links = n.links[
+            (n.links.index.str[:2] == ct) & (n.links.carrier == "shipping oil")
+        ]
+        lhs.append(
+            (
+                n.links_t.p0.loc[:, shipping_links.index].mul(
+                    n.snapshot_weightings.generators, axis=0
+                )
+                * shipping_links.efficiency2
+            ).sum()
+            * domestic_navigation_factor
+        )
+
+        # Shipping methanol
+        shipping_meoh_links = n.links[
+            (n.links.index.str[:2] == ct) & (n.links.carrier == "shipping methanol")
+        ]
+        lhs.append(
+            (
+                n.links_t.p0.loc[:, shipping_meoh_links.index].mul(
+                    n.snapshot_weightings.generators, axis=0
+                )
+                * shipping_meoh_links.efficiency2
+            ).sum()
+            * domestic_navigation_factor
+        )
+
+        # Adding Efuel imports and exports to constraint
+        incoming_oil = n.links.index[n.links.index == f"EU renewable oil -> {ct} oil"]
+        outgoing_oil = n.links.index[n.links.index == f"{ct} renewable oil -> EU oil"]
+
+        lhs.append(
+            (
+                -1
+                * n.links_t.p0.loc[:, incoming_oil].mul(
+                    n.snapshot_weightings.generators, axis=0
+                )
+                * 0.2571
+            ).sum()
+        )
+        lhs.append(
+            (
+                n.links_t.p0.loc[:, outgoing_oil].mul(
+                    n.snapshot_weightings.generators, axis=0
+                )
+                * 0.2571
+            ).sum()
+        )
+
+        incoming_methanol = n.links.index[
+            n.links.index == f"EU methanol -> {ct} methanol"
+        ]
+        outgoing_methanol = n.links.index[
+            n.links.index == f"{ct} methanol -> EU methanol"
+        ]
+
+        lhs.append(
+            (
+                -1
+                * n.links_t.p0.loc[:, incoming_methanol].mul(
+                    n.snapshot_weightings.generators, axis=0
+                )
+                / snakemake.config["sector"]["MWh_MeOH_per_tCO2"]
+            ).sum()
+        )
+
+        lhs.append(
+            (
+                n.links_t.p0.loc[:, outgoing_methanol].mul(
+                    n.snapshot_weightings.generators, axis=0
+                )
+                / snakemake.config["sector"]["MWh_MeOH_per_tCO2"]
+            ).sum()
+        )
+
+        # Methane
+        incoming_CH4 = n.links.index[n.links.index == f"EU renewable gas -> {ct} gas"]
+        outgoing_CH4 = n.links.index[n.links.index == f"{ct} renewable gas -> EU gas"]
+
+        lhs.append(
+            (
+                -1
+                * n.links_t.p0.loc[:, incoming_CH4].mul(
+                    n.snapshot_weightings.generators, axis=0
+                )
+                * 0.198
+            ).sum()
+        )
+
+        lhs.append(
+            (
+                n.links_t.p0.loc[:, outgoing_CH4].mul(
+                    n.snapshot_weightings.generators, axis=0
+                )
+                * 0.198
+            ).sum()
+        )
+
+        return pd.concat(lhs).sum() * t2Mt
+
     energy_totals = _energy_totals.loc[region[0:2]]
 
     industry_DE = industry_demand.filter(
@@ -2587,6 +2748,8 @@ def get_emissions(n, region, _energy_totals, industry_demand):
 
     var = pd.Series()
 
+    var["Emissions|CO2|Model|Constraint"] = get_constraint_emissions(n, region).sum()
+
     co2_emissions = (
         n.statistics.supply(bus_carrier="co2", **kwargs)
         .filter(like=region)