diff --git a/config/config.default.yaml b/config/config.default.yaml index 1bd6e7d386..0b66f0d7b5 100644 --- a/config/config.default.yaml +++ b/config/config.default.yaml @@ -420,6 +420,10 @@ renewable: eia_correct_by_capacity: false eia_approximate_missing: false +# docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#perennials +perennials: + sequestration_co2: 2 + # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#conventional conventional: unit_commitment: false @@ -837,6 +841,7 @@ sector: methanation: true coal_cc: false dac: true + perennials: false co2_vent: false heat_vent: urban central: true @@ -1325,6 +1330,9 @@ data: instrat_co2_prices: source: primary version: latest + co2_removal_data: + source: primary + version: latest co2stop: source: archive version: latest diff --git a/config/plotting.default.yaml b/config/plotting.default.yaml index 694fc4956b..dcaa6fde47 100644 --- a/config/plotting.default.yaml +++ b/config/plotting.default.yaml @@ -741,3 +741,4 @@ plotting: import NH3: '#e2ed74' import oil: '#93eda2' import methanol: '#87d0e6' + co2 perennials: '#008ffc' diff --git a/config/schema.default.json b/config/schema.default.json index a442c12f17..72ac2f3e37 100644 --- a/config/schema.default.json +++ b/config/schema.default.json @@ -1044,6 +1044,26 @@ } } }, + "co2_removal_data": { + "description": "Configuration for a single data source.", + "properties": { + "source": { + "default": "archive", + "description": "Source of the data. 'archive' retrieves pre-built data, 'primary' retrieves from primary source.", + "enum": [ + "archive", + "primary", + "build" + ], + "type": "string" + }, + "version": { + "default": "latest", + "description": "Version of the data to use. Uses the specific 'version' for the selected 'source' or the dataset tagged 'latest' for this source.", + "type": "string" + } + } + }, "co2stop": { "description": "Configuration for a single data source.", "properties": { @@ -2925,6 +2945,16 @@ } } }, + "PerennialsConfig": { + "description": "Configuration for `perennials` settings.", + "properties": { + "sequestration_co2": { + "default": 2, + "description": "Tonnes of CO2 equivalent sequestered per hectare per year when 1st-generation biofuel cropland is converted to perennial grasses.", + "type": "number" + } + } + }, "PypsaEurConfig": { "description": "Configuration for `pypsa_eur` settings.", "properties": { @@ -4647,6 +4677,11 @@ "description": "Add option for Direct Air Capture (DAC).", "type": "boolean" }, + "perennials": { + "default": false, + "description": "Add option for perennialisation (converting 1st-generation biofuel cropland to perennial grasses) as a carbon dioxide removal (CDR) technology.", + "type": "boolean" + }, "co2_vent": { "default": false, "description": "Add option for vent out CO2 from storages to the atmosphere.", @@ -10091,6 +10126,16 @@ } } }, + "perennials": { + "description": "Configuration for `perennials` settings.", + "properties": { + "sequestration_co2": { + "default": 2, + "description": "Tonnes of CO2 equivalent sequestered per hectare per year when 1st-generation biofuel cropland is converted to perennial grasses.", + "type": "number" + } + } + }, "conventional": { "additionalProperties": true, "description": "Configuration for `conventional` settings.", @@ -11139,6 +11184,11 @@ "description": "Add option for Direct Air Capture (DAC).", "type": "boolean" }, + "perennials": { + "default": false, + "description": "Add option for perennialisation (converting 1st-generation biofuel cropland to perennial grasses) as a carbon dioxide removal (CDR) technology.", + "type": "boolean" + }, "co2_vent": { "default": false, "description": "Add option for vent out CO2 from storages to the atmosphere.", @@ -12901,6 +12951,26 @@ } } }, + "co2_removal_data": { + "description": "Configuration for a single data source.", + "properties": { + "source": { + "default": "archive", + "description": "Source of the data. 'archive' retrieves pre-built data, 'primary' retrieves from primary source.", + "enum": [ + "archive", + "primary", + "build" + ], + "type": "string" + }, + "version": { + "default": "latest", + "description": "Version of the data to use. Uses the specific 'version' for the selected 'source' or the dataset tagged 'latest' for this source.", + "type": "string" + } + } + }, "co2stop": { "description": "Configuration for a single data source.", "properties": { diff --git a/data/versions.csv b/data/versions.csv index 067be135f7..6f7028fc86 100644 --- a/data/versions.csv +++ b/data/versions.csv @@ -11,6 +11,7 @@ bidding_zones_electricitymaps,v1.238.0,primary,latest supported,2026-01-22,,http bidding_zones_electricitymaps,v1.238.0,archive,latest supported,2026-01-22,,https://data.pypsa.org/workflows/eur/bidding_zones_electricitymaps/v1.238.0/world.geojson bidding_zones_entsoepy,v0.6.18,primary,latest supported,2026-01-22,,https://raw.githubusercontent.com/EnergieID/entsoe-py/refs/tags/V0.6.18/entsoe/geo/geojson bidding_zones_entsoepy,v0.6.18,archive,latest supported,2026-01-22,,https://data.pypsa.org/workflows/eur/bidding_zones_entsoepy/v0.6.18 +co2_removal_data,v1.0.1,primary,latest supported,2026-06-22,"Afforestation and perennialisation input data for the CDR sector technologies, archived from https://github.com/BertoGBG/CO2-removal",https://zenodo.org/records/20799337/files/BertoGBG/CO2-removal-v1.0.1.zip co2stop,26-august-2020,primary,latest supported,2025-12-02,,https://setis.ec.europa.eu/document/download/786a884f-0b33-4789-b744-28004b16bd1a_en?filename=co2jrc_openformats.zip co2stop,26-august-2020,archive,latest supported,2026-01-13,,https://data.pypsa.org/workflows/eur/co2stop/26-august-2020/co2jrc_openformats.zip copernicus_land_cover,v3.0.1,primary,latest supported,2025-12-02,"The primary is already from Zenodo, documentation in https://zenodo.org/records/4723921",https://zenodo.org/records/3939050/files/PROBAV_LC100_global_v3.0.1_2019-nrt_Discrete-Classification-map_EPSG-4326.tif diff --git a/doc/release_notes.md b/doc/release_notes.md index 8fcdee22ae..c699ab3e2b 100644 --- a/doc/release_notes.md +++ b/doc/release_notes.md @@ -5,6 +5,7 @@ +* feat: Add perennialisation as a carbon dioxide removal (CDR) technology, converting 1st-generation biofuel cropland to perennial grasses, with node-level potential derived from NUTS2-resolved Eurostat crop yields ([#2143](https://github.com/PyPSA/pypsa-eur/issues/2143)). ## PyPSA-Eur v2026.08.0 (19th August 2026) diff --git a/doc/sector.md b/doc/sector.md index 549a5b98b4..edc48bb26b 100644 --- a/doc/sector.md +++ b/doc/sector.md @@ -35,6 +35,14 @@ Having downloaded the necessary data, ::: build_biomass_potentials +## Rule `build_perennials_yields_eurostat_average` + +::: build_perennials_yields_eurostat_average + +## Rule `build_perennials_yields` + +::: build_perennials_yields + ## Rule `build_egs_potentials` ::: build_egs_potentials diff --git a/rules/build_sector.smk b/rules/build_sector.smk index ce688fafef..6622b6fcb7 100755 --- a/rules/build_sector.smk +++ b/rules/build_sector.smk @@ -953,6 +953,38 @@ rule build_biomass_potentials: scripts("build_biomass_potentials.py") +rule build_perennials_yields_eurostat_average: + input: + nuts2021=rules.retrieve_eu_nuts_2021.output.shapes_level_2, + crops_nuts2=rules.retrieve_co2_removal_data.output.eurostat_crops_nuts2, + crops_nuts0=rules.retrieve_co2_removal_data.output.eurostat_crops_nuts0, + costs=resources(f"costs_{config['costs']['year']}_processed.csv"), + output: + yields_all=resources("perennials_yields_1G_biofuels.csv"), + log: + logs("build_perennials_yields_eurostat_average.log"), + script: + scripts("build_perennials_yields_eurostat_average.py") + + +rule build_perennials_yields: + input: + nuts2=rules.retrieve_eu_nuts_2021.output.shapes_level_2, + country_shapes=resources("country_shapes.geojson"), + perennials_yields_1G_biofuels=resources("perennials_yields_1G_biofuels.csv"), + regions_onshore=resources("regions_onshore_base_s_{clusters}.geojson"), + output: + csv_file=resources("perennials_yields_1G_biofuels_s_{clusters}.csv"), + log: + logs("build_perennials_yields_s_{clusters}.log"), + resources: + mem_mb=8000, + params: + biomass=config_provider("biomass"), + script: + scripts("build_perennials_yields.py") + + rule build_biomass_transport_costs: input: sc1="data/biomass_transport_costs_supplychain1.csv", @@ -1722,6 +1754,11 @@ rule prepare_sector_network: if config_provider("sector", "district_heating", "ates", "enable")(w) else [] ), + perennials_yields_1G_biofuels=lambda w: ( + resources("perennials_yields_1G_biofuels_s_{clusters}.csv") + if config_provider("sector", "perennials")(w) + else [] + ), output: resources( "networks/base_s_{clusters}_{opts}_{sector_opts}_{planning_horizons}.nc" diff --git a/rules/retrieve.smk b/rules/retrieve.smk index 3a96cbc265..968a1d4049 100755 --- a/rules/retrieve.smk +++ b/rules/retrieve.smk @@ -1629,3 +1629,61 @@ if (MOBILITY_PROFILES_DATASET := dataset_version("mobility_profiles"))["source"] run: copy2(input["kfz"], output["kfz"]) copy2(input["pkw"], output["pkw"]) + + +if (CO2_REMOVAL_DATASET := dataset_version("co2_removal_data"))["source"] in [ + "primary", + "archive", +]: + + rule retrieve_co2_removal_data: + input: + zip=storage(CO2_REMOVAL_DATASET["url"]), + output: + afforestation_nuts_biomass_densities=resources( + "afforestation_nuts_biomass_densities.xlsx" + ), + afforestation_nuts2_afforestation_rates=resources( + "afforestation_rates_nuts2_full.csv" + ), + afforestation_nuts2_monthly_weights=resources( + "afforestation_nuts2_monthly_weights.csv" + ), + eurostat_crops_nuts2=resources("eurostat_apro_cpshr_nuts2_raw.csv"), + eurostat_crops_nuts0=resources("eurostat_apro_cpshr_nuts0_raw.csv"), + retries: 2 + message: + "Downloading carbon dioxide removal data (afforestation, perennialisation inputs)" + run: + with ZipFile(input.zip) as z: + # GitHub's release archive nests everything under a single + # top-level "--/" folder whose name + # changes with every release, so resolve it at runtime. + top_dir = z.namelist()[0].split("/")[0] + for src_path, dest in [ + ( + "outputs/afforestation/afforestation_nuts_biomass_densities.xlsx", + output.afforestation_nuts_biomass_densities, + ), + ( + "outputs/afforestation/afforestation_rates_nuts2_full.csv", + output.afforestation_nuts2_afforestation_rates, + ), + ( + "outputs/afforestation/afforestation_nuts2_monthly_weights.csv", + output.afforestation_nuts2_monthly_weights, + ), + ( + "outputs/perennialisation/eurostat_apro_cpshr_nuts2_raw.csv", + output.eurostat_crops_nuts2, + ), + ( + "outputs/perennialisation/eurostat_apro_cpshr_nuts0_raw.csv", + output.eurostat_crops_nuts0, + ), + ]: + with ( + z.open(f"{top_dir}/{src_path}") as src, + open(dest, "wb") as dst, + ): + dst.write(src.read()) diff --git a/scripts/_helpers.py b/scripts/_helpers.py index a606080ad8..cbecd2a153 100644 --- a/scripts/_helpers.py +++ b/scripts/_helpers.py @@ -1093,6 +1093,51 @@ def load_costs(cost_file: str) -> pd.DataFrame: return pd.read_csv(cost_file, index_col=0) +# 1G biofuel crop groups and their target biomass class names. +# In the upstream default config these sit in "not included". +# When sector: perennials=True, resolve_biomass_classes() moves them automatically. +ONE_G_BIOFUEL_CLASSES = { + "Bioethanol barley, wheat, grain maize, oats, other cereals and rye": "biofuels_1G_bioethanol_cereals", + "Sugar from sugar beet": "biofuels_1G_bioethanol_sugar", + "Rape seed": "biofuels_1G_biodiesel", + "Sunflower, soya seed ": "biofuels_1G_biodiesel", +} + + +def resolve_biomass_classes(classes, perennials_enabled): + """ + Auto-reallocate 1G biofuel crop groups from 'not included' into their + biofuels_1G_* target classes when perennials are enabled. + + Cases: + perennials=False → return classes unchanged. + perennials=True, all groups in 'not included' → move them (normal case). + perennials=True, no groups in 'not included' → raise AssertionError. + perennials=True, some groups elsewhere → move available ones, warn about rest. + """ + if not perennials_enabled: + return classes + classes = copy.deepcopy(classes) + not_incl = classes.get("not included", []) + in_ni = [g for g in ONE_G_BIOFUEL_CLASSES if g in not_incl] + elsewhere = [g for g in ONE_G_BIOFUEL_CLASSES if g not in not_incl] + if not in_ni: + raise AssertionError( + "sector: perennials=true but no 1G-biofuel groups are in " + "biomass: classes: 'not included'. Restore upstream defaults " + "so the groups can be reallocated automatically." + ) + if elsewhere: + logger.warning( + "perennials: 1G groups already allocated outside 'not included' " + f"— skipped (running perennials only for available groups): {elsewhere}" + ) + for g in in_ni: + classes["not included"].remove(g) + classes.setdefault(ONE_G_BIOFUEL_CLASSES[g], []).append(g) + return classes + + def _simplify_polys( polys, minarea=100 * 1e6, maxdistance=None, tolerance=None, filterremote=True ): # 100*1e6 = 100 km² if CRS is DISTANCE_CRS diff --git a/scripts/build_biomass_potentials.py b/scripts/build_biomass_potentials.py index 157c06cc4f..aa163338cd 100755 --- a/scripts/build_biomass_potentials.py +++ b/scripts/build_biomass_potentials.py @@ -12,7 +12,11 @@ import numpy as np import pandas as pd -from scripts._helpers import configure_logging, set_scenario_config +from scripts._helpers import ( + configure_logging, + resolve_biomass_classes, + set_scenario_config, +) logger = logging.getLogger(__name__) AVAILABLE_BIOMASS_YEARS = [2010, 2020, 2030, 2040, 2050] @@ -184,17 +188,15 @@ def disaggregate_nuts0(bio): return bio -def build_nuts2_shapes(): +def build_nuts2_shapes(nuts2_fn, country_shapes_fn): """ - load NUTS2 geometries - add RS, AL, BA country shapes (not covered in NUTS 2013) - consistently name ME, MK """ - nuts2 = gpd.GeoDataFrame( - gpd.read_file(snakemake.input.nuts2).set_index("NUTS_ID").geometry - ) + nuts2 = gpd.GeoDataFrame(gpd.read_file(nuts2_fn).set_index("NUTS_ID").geometry) - countries = gpd.read_file(snakemake.input.country_shapes).set_index("name") + countries = gpd.read_file(country_shapes_fn).set_index("name") missing_iso2 = countries.index.intersection(["AL", "RS", "XK", "BA"]) missing = countries.loc[missing_iso2] @@ -382,7 +384,7 @@ def add_unsustainable_potentials(df, input_eurostat): enspreso = disaggregate_nuts0(enspreso) - nuts2 = build_nuts2_shapes() + nuts2 = build_nuts2_shapes(snakemake.input.nuts2, snakemake.input.country_shapes) df_nuts2 = gpd.GeoDataFrame(nuts2.geometry).join(enspreso) @@ -392,7 +394,10 @@ def add_unsustainable_potentials(df, input_eurostat): df.to_csv(snakemake.output.biomass_potentials_all) - grouper = {v: k for k, vv in params["classes"].items() for v in vv} + classes = resolve_biomass_classes( + params["classes"], snakemake.config["sector"].get("perennials", False) + ) + grouper = {v: k for k, vv in classes.items() for v in vv} df = df.T.groupby(grouper).sum().T input_eurostat = snakemake.input.eurostat diff --git a/scripts/build_perennials_yields.py b/scripts/build_perennials_yields.py new file mode 100644 index 0000000000..2188321b22 --- /dev/null +++ b/scripts/build_perennials_yields.py @@ -0,0 +1,242 @@ +# SPDX-FileCopyrightText: Contributors to PyPSA-Eur +# +# SPDX-License-Identifier: MIT +""" +Reproject NUTS2-level 1st-generation (1G) biofuel and perennial crop YIELDS +(from ``build_perennials_yields_eurostat_average.py``) onto the clustered network +regions, for use by ``add_perennials()`` in ``prepare_sector_network.py``. + +Output units, NOT an area or a potential +----------------------------------------- +The output CSV holds per-hectare yield RATES, not land area or a CO2 +potential: MWh/ha/y for the 1G biofuel crop columns (cereals, sugar beet, +rapeseed) and t/ha/y for the ``perennials`` column. The land area available +for conversion to perennial grasses, and the resulting CO2 sequestration +potential, are only derived later - see "Downstream usage" below. + +Method +------ +NUTS2 yields are mapped to clustered regions via an area-weighted average +over the NUTS2/region overlay (NUTS2 geometries intersected with cluster +region geometries, weighted by intersection area) - see +``convert_nuts2_to_regions_yields()``. NUTS2 regions not covered by +Eurostat data (non-EU countries, small islands, city-states) are first +filled from the nearest valid NUTS2 centroid - see +``impute_missing_values()``. Resulting columns are then grouped into +biomass classes (``resolve_biomass_classes()``) to match the class +structure used by ``build_biomass_potentials.py``. + +Downstream usage +----------------- +``add_perennials()`` divides the ENSPRESO 1G biomass potential (MWh/y, an +extensive quantity from ``build_biomass_potentials.py``) by this script's +1G yields (MWh/ha/y) to back out the land area (ha) available for +conversion to perennial grasses, then converts that area into a CO2 +sequestration potential via ``perennials.sequestration_co2`` (tCO2e/ha/y). + +Outputs a single CSV with one column per crop class (cereals, sugar beet, +rapeseed, perennials) indexed by clustered region name. +""" + +import logging + +import geopandas as gpd +import pandas as pd + +from scripts._helpers import ( + configure_logging, + resolve_biomass_classes, + set_scenario_config, +) +from scripts.build_biomass_potentials import build_nuts2_shapes + +logger = logging.getLogger(__name__) + + +def impute_missing_values(df_nuts2, missing_shapes, yield_cols): + """ + Fill missing yield values for NUTS2 regions not covered by Eurostat data + (non-EU countries, small islands, city-states) by copying the values from + the nearest NUTS2 region that does have valid (non-NaN) yields. + + Nearest is determined by centroid distance in an equal-area CRS (EPSG:3035). + + Parameters + ---------- + df_nuts2 : gpd.GeoDataFrame + NUTS2 geometries joined with yield columns; rows for regions without + Eurostat coverage are entirely NaN in ``yield_cols``. + missing_shapes : gpd.GeoDataFrame + Geometries of the NUTS2 (or country-level substitute) regions to + impute, indexed the same way as ``df_nuts2``. + yield_cols : list of str + Columns in ``df_nuts2`` to impute. + + Returns + ------- + gpd.GeoDataFrame + One row per entry in ``missing_shapes``, with ``yield_cols`` filled + from the nearest valid NUTS2 neighbour. + """ + # Keep only rows that have valid yields (drop NaN rows!) + df_valid = df_nuts2.dropna(subset=yield_cols).copy() + + # Project to 3035 + valid_proj = df_valid.to_crs(3035) + missing_proj = missing_shapes.to_crs(3035) + + valid_centroids = valid_proj.centroid + missing_centroids = missing_proj.centroid + + imputed_rows = [] + + for missing_id, c_geom in missing_centroids.items(): + # Distance to only VALID NUTS2 rows + dists = valid_centroids.distance(c_geom) + nearest_nuts2 = dists.idxmin() + + # Copy yields + attrs = df_valid.loc[nearest_nuts2, yield_cols] + + # Create new row + new_row = missing_shapes.copy().loc[[missing_id]] + for col in yield_cols: + new_row[col] = attrs[col] + + imputed_rows.append(new_row) + + return pd.concat(imputed_rows) + + +def convert_nuts2_to_regions_yields(df_nuts2, regions, yield_cols=None): + """ + Convert NUTS2-level yields (intensive, e.g. MWh/ha/y or t/ha/y) to + PyPSA-Eur clustered regions via an area-weighted average over the + NUTS2/region overlay: + + y_n = Σ_i (y_i * A_i∩n) / Σ_i A_i∩n + + Unlike an extensive quantity (a total, e.g. MWh), a yield is a rate and + must be area-weighted-averaged rather than redistributed by area share + - see ``convert_nuts2_to_regions`` in ``build_biomass_potentials.py`` for + the extensive-quantity equivalent. + + Only NUTS2 rows with non-NaN yields are used. Regions with no + overlapping valid NUTS2 get NaN. + + Parameters + ---------- + df_nuts2 : gpd.GeoDataFrame + NUTS2 geometries joined with yield columns. + regions : gpd.GeoDataFrame + PyPSA-Eur clustered onshore regions, with a ``name`` column. + yield_cols : list of str, optional + Columns in ``df_nuts2`` to convert. Defaults to all columns except + ``geometry`` and ``NUTS_ID``. + + Returns + ------- + pd.DataFrame + Area-weighted-average yields indexed by region name, one row per + entry in ``regions``. + """ + + nuts = df_nuts2.copy() + regs = regions.copy() + + # Ensure ID column for NUTS2 + nuts["NUTS_ID"] = nuts.index + + # Identify yield columns if not given + if yield_cols is None: + yield_cols = nuts.columns.difference(["geometry", "NUTS_ID"]) + + # 1) Keep only NUTS2 rows that have at least one non-NaN yield + nuts_valid = nuts.dropna(subset=yield_cols, how="all") + + # 2) Reproject to equal-area CRS for areas + nuts_valid = nuts_valid.to_crs(3035) + regs = regs.to_crs(3035) + + # 3) Overlay: intersection between regions and valid NUTS2 + overlay = gpd.overlay(regs, nuts_valid, keep_geom_type=True) + + # Return empty if nothing overlaps + if overlay.empty: + return pd.DataFrame(index=regs["name"].values, columns=yield_cols, dtype=float) + + # 4) Area of intersections (m² → km²) + overlay["area_intersection"] = overlay.geometry.area + + # Optional: drop absurdly tiny slivers + overlay = overlay[overlay["area_intersection"] > 0] + + # 5) Numerators: y_i * A_i∩n + numerators = overlay[yield_cols].multiply(overlay["area_intersection"], axis=0) + + # 6) Sum per region + numer_by_region = numerators.groupby(overlay["name"]).sum(min_count=1) + denom_by_region = overlay.groupby("name")["area_intersection"].sum() + + # 7) Weighted average + yields_regions = numer_by_region.div(denom_by_region, axis=0) + + # 8) Ensure one row per region (regions with no valid overlaps → NaN) + yields_regions = yields_regions.reindex(regs["name"].values) + + return yields_regions + + +if __name__ == "__main__": + if "snakemake" not in globals(): + from scripts._helpers import mock_snakemake + + snakemake = mock_snakemake( + "build_perennials_yields", + clusters="39", + planning_horizons=2050, + ) + + configure_logging(snakemake) + set_scenario_config(snakemake) + + regions = gpd.read_file(snakemake.input.regions_onshore) + nuts2 = build_nuts2_shapes(snakemake.input.nuts2, snakemake.input.country_shapes) + + yields = pd.read_csv(snakemake.input.perennials_yields_1G_biofuels, index_col=0) + + df_nuts2 = gpd.GeoDataFrame(nuts2.geometry).join(yields) + + # Address the missign countries (df_nuts2 contains all NUTS2 + missing shapes (with NaNs)) + missing_countries = ["AL", "RS", "BA", "XK"] + missing_shapes = df_nuts2.loc[missing_countries, ["geometry"]] + + # Impute yields for missing shapes basest of nearest valid Nuts2 + yield_cols = [ + "Bioethanol barley, wheat, grain maize, oats, other cereals and rye", + "Sugar from sugar beet", + "Rape seed", + "perennials", + ] + imputed_missing = impute_missing_values(df_nuts2, missing_shapes, yield_cols) + df_nuts2 = pd.concat([df_nuts2.drop(index=missing_countries), imputed_missing]) + + # convert nuts2 yields to regions + df = convert_nuts2_to_regions_yields(df_nuts2, regions) + + params = snakemake.params.biomass + + classes = resolve_biomass_classes( + params["classes"], snakemake.config["sector"].get("perennials", False) + ) + grouper = {v: k for k, vv in classes.items() for v in vv} + + # keep the column 'perennials' which is otherwise dropped + for col in df.columns: + if col not in grouper: + grouper[col] = col + + df = df.T.groupby(grouper).sum().T + + df.index.name = "name" + df.to_csv(snakemake.output.csv_file) diff --git a/scripts/build_perennials_yields_eurostat_average.py b/scripts/build_perennials_yields_eurostat_average.py new file mode 100644 index 0000000000..3d80818659 --- /dev/null +++ b/scripts/build_perennials_yields_eurostat_average.py @@ -0,0 +1,399 @@ +# SPDX-FileCopyrightText: Contributors to PyPSA-Eur +# +# SPDX-License-Identifier: MIT +""" +Retrieve crop harvest data from the Eurostat API (dataset ``apro_cpshr``) at +NUTS2 and NUTS0 resolution, and compute production-weighted average yields +(t/ha, then converted to MWh/ha for 1G biofuel crops) per NUTS2 region for +1st-generation biofuel feedstocks and perennial grasses. + +"Eurostat average" in the name refers to this production-weighted averaging +across crop codes and years - not the geometric area-weighted overlay used +downstream in ``build_perennials_yields.py`` to reproject these NUTS2-level +yields onto clustered network regions; do not confuse the two. + +Missing NUTS2 coverage is filled via a three-tier fallback (NUTS2 data -> +NUTS0 country-level data -> spatial neighbor mean or nearest valid NUTS2 for +islands - see ``harmonize_to_nuts2021()``), and results are harmonized to +the NUTS2021 region definitions used by PyPSA-Eur. + +Outputs a single CSV with columns for each crop class (cereals, sugar beet, +rapeseed, perennials) indexed by NUTS2 region. + +Biofuel conversion efficiencies (t_biofuel / t_feedstock) are read from the +``efficiency`` parameter of the ``ethanol from wheat``, ``ethanol from sugar +beet``, and ``biodiesel from rapeseed`` technologies in the technology-data +cost assumptions, sourced from: + + Banja et al. (2013), "Biofuels in the European Union - A general overview", + JRC Technical Report, doi:10.2760/69179, Tables 93, 133, 155, 159. +""" + +import logging +from pathlib import Path + +import geopandas as gpd +import numpy as np +import pandas as pd + +from scripts._helpers import load_costs + +logger = logging.getLogger(__name__) + + +def harmonize_to_nuts2021(df, keep_col, nuts2021_n2): + """ + Df : DataFrame indexed by ['geo', 'TIME_PERIOD', 'mapping'] + contains both NUTS2 and NUTS0 rows + keep_col : column to harmonize (e.g. 'weighted_YL_(t/ha)') + nuts2021_n2 : GeoDataFrame with index=NUTS2_ID and geometry + """ + + # NUTS2 target regions + target = nuts2021_n2.index.sort_values() + + # Split df + df = df.copy() + + # Identify NUTS2 and NUTS0 by index length + df["is_nuts2"] = df.index.str.len() == 4 + + # Split into NUTS2 and NUTS0 + df_nuts2 = df[df["is_nuts2"]] + df_nuts0 = df[~df["is_nuts2"]] + + # Prepare working layer for only NUTS2 + df_work = df_nuts2[[keep_col]].reindex(target) + + # Country code extraction (first 2 chars) + df_work["country"] = df_work.index.str[:2] + + # Fallback 2 — Use NUTS-0 values + df_work = df_work.join( + df_nuts0[[keep_col]].rename(columns={keep_col: "fallback"}), on="country" + ) + df_work[keep_col] = df_work[keep_col].fillna(df_work["fallback"]) + df_work.drop(columns=["fallback"], inplace=True) + + # Join geometry + nuts_proj = nuts2021_n2.to_crs(epsg=3035) + gdf = nuts_proj.join(df_work)[[keep_col, "country", "geometry"]] + + # Fallback 3 — Spatial neighbors mean or nearest region if island + mask_missing = gdf[keep_col].isna() | (gdf[keep_col] <= 0) + if mask_missing.any(): + logger.warning("Spatial fallback required for %d regions", mask_missing.sum()) + + # Pre-calc distance matrix only once + valid = gdf[gdf[keep_col].notna() & (gdf[keep_col] > 0)] + + for idx in gdf[mask_missing].index: + region = gdf.loc[idx, "geometry"] + + # Touching neighbors + neigh_idxs = gdf[gdf.geometry.touches(region)].index.tolist() + neigh_vals = gdf.loc[neigh_idxs, keep_col].dropna() + neigh_vals = neigh_vals[neigh_vals > 0] + + if len(neigh_vals) > 0: + gdf.at[idx, keep_col] = neigh_vals.mean() + else: + # Island fallback: nearest valid NUTS2 region + nearest_idx = valid.distance(region).idxmin() + gdf.at[idx, keep_col] = valid.at[nearest_idx, keep_col] + logger.debug("Spatial fallback for %s: nearest = %s", idx, nearest_idx) + + # Final output tidy + result = gdf[[keep_col]] + result.index.name = "NUTS2" + result.sort_index() + return result + + +def calculate_yields( + filepath_nuts2, filepath_nuts0, crops_sel, crops_mapping, biofuel_yields +): + # filter columns - keep only relevant + df_crops_raw_nuts2 = pd.read_csv(filepath_nuts2) + df_crops_raw_nuts2["TIME_PERIOD"] = df_crops_raw_nuts2["TIME_PERIOD"].astype(int) + + df_crops_raw_nuts0 = pd.read_csv(filepath_nuts0) + df_crops_raw_nuts0["TIME_PERIOD"] = df_crops_raw_nuts0["TIME_PERIOD"].astype(int) + + df_crops_raw = pd.concat( + [df_crops_raw_nuts0, df_crops_raw_nuts2], ignore_index=True + ) + + # drop empty and irrelevant columns + columns_to_drop = [ + "Observation value", + "OBS_FLAG", + "Observation status (Flag) V2 structure", + "CONF_STATUS", + "Confidentiality status (flag)", + "Time", + "STRUCTURE_ID", + "STRUCTURE", + "STRUCTURE_NAME", + "Geopolitical entity (reporting)", + "Time frequency", + ] + + # useful columns: + # Index(['freq', 'crops', 'Crops', 'strucpro', 'Structure of production', 'geo','TIME_PERIOD', 'OBS_VALUE'], + # freq : 'A' meaning annual + # 'crops' : code of the crops e.g. 'G0000', 'G1000', 'G2000', 'G2100', 'G2900' + # 'Crops' : name of the corp e.g. Sugar beet (excluding seed + # 'strucpro' : type of data ['AR', 'MA', 'PR_HU_EU'], where AR is cultivated area, MA is and PR_HU_EU is production at standard EU humidity + # 'OBS_VALUE' : numerical value + + df_crops = df_crops_raw.drop(columns=columns_to_drop, errors="ignore") + df_crops["OBS_VALUE"] = df_crops["OBS_VALUE"].fillna(0) + + # Step 1: Filter to relevant rows for 2023 and strucpro of interest + df_sub = df_crops[ + (df_crops["strucpro"].isin(["AR", "PR_HU_EU"])) + & (df_crops["crops"].isin(crops_sel)) + ][["crops", "geo", "TIME_PERIOD", "strucpro", "OBS_VALUE"]] + + # Pivot so AR and PR_HU_EU are columns for each (crop, geo, year) + df_pivot = ( + df_sub.pivot_table( + index=["crops", "geo", "TIME_PERIOD"], + columns="strucpro", + values="OBS_VALUE", + ) + .dropna(subset=["AR", "PR_HU_EU"]) + .reset_index() + ) + + # Compute yield per year (PR_HU_EU / AR) + df_pivot["YL_(t/ha)"] = np.divide( + df_pivot["PR_HU_EU"], + df_pivot["AR"], + out=np.zeros_like(df_pivot["PR_HU_EU"], dtype=float), + where=df_pivot["AR"] != 0, + ) + + # Average data across years + df_avg_yield = df_pivot.groupby(["crops", "geo"], as_index=False)[ + ["AR", "PR_HU_EU", "YL_(t/ha)"] + ].mean() + + min_year = df_pivot["TIME_PERIOD"].min() + max_year = df_pivot["TIME_PERIOD"].max() + df_avg_yield["TIME_PERIOD"] = f"{min_year}-{max_year}" + + # map crops to categories unsustainable biofuels in + rev_map = { + code: key + for key, val in crops_mapping.items() + for code in (val if isinstance(val, list) else [val]) + } + df_avg_yield["mapping"] = df_avg_yield["crops"].map(rev_map) + + # calculate weighted production per crop within mapping classes + df_avg_yield["PR_share"] = df_avg_yield["PR_HU_EU"] / df_avg_yield.groupby( + ["geo", "TIME_PERIOD", "mapping"] + )["PR_HU_EU"].transform("sum") + df_avg_yield["PR_share"] = df_avg_yield["PR_share"].fillna(0) + + # calculated average weighted yield + df_avg_yield["weighted_YL_(t/ha)"] = ( + df_avg_yield["YL_(t/ha)"] * df_avg_yield["PR_share"] + ) + + # sanity check for very low yields due to small productions + thresholds = { + "MINBIOCRP11": 2.0, # cereals + "MINBIOCRP21": 50.0, # sugar beet + "MINBIORPS1": 1.5, # rapeseed + "PERENNIALS": 5.0, # perennial grasses + } + + # Apply crop-specific minimum threshold + df_avg_yield["weighted_YL_(t/ha)"] = df_avg_yield.apply( + lambda row: ( + row["weighted_YL_(t/ha)"] + if row["weighted_YL_(t/ha)"] >= thresholds.get(row["mapping"], 0) + else 0 + ), + axis=1, + ) + + # weighted yields from current production : applies to unsustainable biofuels + weighted_yields = df_avg_yield.groupby(["geo", "TIME_PERIOD", "mapping"])[ + "weighted_YL_(t/ha)" + ].sum() + + unsustainable_biofuels_yields = pd.DataFrame(weighted_yields) + + # unsustainable biomass yield units from t/ha to MWh/ha + unsustainable_biofuels_yields = unsustainable_biofuels_yields[ + unsustainable_biofuels_yields.index.get_level_values("mapping") != "PERENNIALS" + ] + + unsustainable_biofuels_yields["energy_yields_(MWh/ha)"] = ( + unsustainable_biofuels_yields["weighted_YL_(t/ha)"] + * unsustainable_biofuels_yields.index.get_level_values("mapping").map( + biofuel_yields + ) + ) + + # yields of perennials per hectare in ton/ha + # standard humidity for perennials = 0.65 (tH2O/t_fresh) -> note production is for fresh until 2025 + std_moist_perennials = 0.65 + + perennial_yields = pd.DataFrame(weighted_yields) + perennial_yields = perennial_yields[ + perennial_yields.index.get_level_values("mapping") == "PERENNIALS" + ] * (1 - std_moist_perennials) + + # max yields from current production : applies to perennials for green biorefining + max_yields = df_avg_yield.groupby(["geo", "TIME_PERIOD", "mapping"])[ + "YL_(t/ha)" + ].max() + + perennial_yields_max = pd.DataFrame(max_yields) + perennial_yields_max = perennial_yields_max[ + perennial_yields_max.index.get_level_values("mapping") == "PERENNIALS" + ] * (1 - std_moist_perennials) + + return unsustainable_biofuels_yields, perennial_yields, perennial_yields_max + + +if __name__ == "__main__": + if "snakemake" not in globals(): + from scripts._helpers import mock_snakemake + + snakemake = mock_snakemake("build_perennials_yields_eurostat_average") + + from scripts._helpers import configure_logging, set_scenario_config + + configure_logging(snakemake) + set_scenario_config(snakemake) + + CROPS_CSV_NUTS2 = Path(snakemake.input["crops_nuts2"]) + CROPS_CSV_NUTS0 = Path(snakemake.input["crops_nuts0"]) + NUTS2_2021_GEOJSON = Path(snakemake.input["nuts2021"]) + OUT_CSV_YIELDS_ALL = Path(snakemake.output["yields_all"]) + + CROPS_CSV_NUTS2.parent.mkdir(parents=True, exist_ok=True) + OUT_CSV_YIELDS_ALL.parent.mkdir(parents=True, exist_ok=True) + + # G0000: total green plants (kept – often the only G-code with NUTS2 coverage, e.g. DK) + # G1000: temporary grasses and grazings + # G2000: aggregate legumes (G2100 + G2900) – included for NUTS0 fallback coverage; + # does NOT inflate MAX because it is always ≤ max(G2100, G2900) + # G2100: lucerne/alfalfa; G2900: clover and other leguminous plants + perennial_codes = ["G0000", "G1000", "G2000", "G2100", "G2900"] + + crops_mapping = dict( + MINBIOCRP11=["C0000", "C1000", "C1210", "C1300", "C1310", "C1320"], + MINBIOCRP21="R2000", + MINBIORPS1=["I1110", "I1120", "I1130", "I1110-1130", "I0000"], + PERENNIALS=perennial_codes, + ) + costs = load_costs(snakemake.input.costs) + + # conv = snakemake.params.biofuel_conversion + # LHV values are fixed physical constants, not parameters: JRC Technical Report doi:10.2760/69179 + LHV_fuels = { + "ethanol": 7.447, + "biodiesel": 10.194, + } # MWh/t (26.81 MJ/kg, 36.7 MJ/kg) + + biofuel_yields = { + "MINBIOCRP11": costs.at["ethanol from wheat", "efficiency"] + * LHV_fuels["ethanol"], + "MINBIOCRP21": costs.at["ethanol from sugar beet", "efficiency"] + * LHV_fuels["ethanol"], + "MINBIORPS1": costs.at["biodiesel from rapeseed", "efficiency"] + * LHV_fuels["biodiesel"], + } + + other_crops_codes = [ + item + for v in crops_mapping.values() + for item in (v if isinstance(v, list) else [v]) + ] + crops_sel = perennial_codes + other_crops_codes + + logger.info("Computing crop yields...") + unsustainable_biofuels_yields, perennial_yields, perennial_yields_max = ( + calculate_yields( + filepath_nuts0=CROPS_CSV_NUTS0, + filepath_nuts2=CROPS_CSV_NUTS2, + crops_sel=crops_sel, + crops_mapping=crops_mapping, + biofuel_yields=biofuel_yields, + ) + ) + + yield_MINBIOCRP11 = unsustainable_biofuels_yields[ + unsustainable_biofuels_yields.index.get_level_values("mapping") == "MINBIOCRP11" + ] + yield_MINBIOCRP21 = unsustainable_biofuels_yields[ + unsustainable_biofuels_yields.index.get_level_values("mapping") == "MINBIOCRP21" + ] + yield_MINBIORPS1 = unsustainable_biofuels_yields[ + unsustainable_biofuels_yields.index.get_level_values("mapping") == "MINBIORPS1" + ] + + yield_MINBIOCRP11 = yield_MINBIOCRP11.droplevel(["TIME_PERIOD", "mapping"]) + yield_MINBIOCRP11.index.name = "NUTS2" + yield_MINBIOCRP21 = yield_MINBIOCRP21.droplevel(["TIME_PERIOD", "mapping"]) + yield_MINBIOCRP21.index.name = "NUTS2" + yield_MINBIORPS1 = yield_MINBIORPS1.droplevel(["TIME_PERIOD", "mapping"]) + yield_MINBIORPS1.index.name = "NUTS2" + perennial_yields = perennial_yields.droplevel(["TIME_PERIOD", "mapping"]) + perennial_yields.index.name = "NUTS2" + perennial_yields_max = perennial_yields_max.droplevel(["TIME_PERIOD", "mapping"]) + perennial_yields_max.index.name = "NUTS2" + + logger.info("Harmonizing to NUTS2021 regions...") + nuts2021_n2 = ( + gpd.read_file(NUTS2_2021_GEOJSON) + .loc[:, ["NUTS_ID", "NUTS_NAME", "CNTR_CODE", "geometry"]] + .set_index("NUTS_ID") + ) + + yield_MINBIOCRP11_full = harmonize_to_nuts2021( + yield_MINBIOCRP11, "energy_yields_(MWh/ha)", nuts2021_n2 + ) + yield_MINBIOCRP21_full = harmonize_to_nuts2021( + yield_MINBIOCRP21, "energy_yields_(MWh/ha)", nuts2021_n2 + ) + yield_MINBIORPS1_full = harmonize_to_nuts2021( + yield_MINBIORPS1, "energy_yields_(MWh/ha)", nuts2021_n2 + ) + yields_perennials_max_full = harmonize_to_nuts2021( + perennial_yields_max, "YL_(t/ha)", nuts2021_n2 + ) + yields_perennials_full = harmonize_to_nuts2021( + perennial_yields, "weighted_YL_(t/ha)", nuts2021_n2 + ) + + df_yields_all = pd.concat( + { + "MINBIOCRP11": yield_MINBIOCRP11_full["energy_yields_(MWh/ha)"], + "MINBIOCRP21": yield_MINBIOCRP21_full["energy_yields_(MWh/ha)"], + "MINBIORPS1": yield_MINBIORPS1_full["energy_yields_(MWh/ha)"], + # Max yield across G-codes: models the best-available perennial crop + # choice in each region when substituting 1G biofuel crops. + "PERENNIALS_MAX": yields_perennials_max_full["YL_(t/ha)"], + }, + axis=1, + ) + df_yields_all.columns = [ + "Bioethanol barley, wheat, grain maize, oats, other cereals and rye", + "Sugar from sugar beet", + "Rape seed", + "perennials", + ] + df_yields_all = df_yields_all.sort_index() + + logger.info("Saving output CSV files...") + df_yields_all.to_csv(OUT_CSV_YIELDS_ALL, index=True) + + logger.info("Done.") diff --git a/scripts/lib/validation/config/_schema.py b/scripts/lib/validation/config/_schema.py index ad9e941c00..9f982ba4a7 100644 --- a/scripts/lib/validation/config/_schema.py +++ b/scripts/lib/validation/config/_schema.py @@ -27,6 +27,7 @@ from scripts.lib.validation.config.links import LinksConfig from scripts.lib.validation.config.load import LoadConfig from scripts.lib.validation.config.overpass_api import OverpassApiConfig +from scripts.lib.validation.config.perennials import PerennialsConfig from scripts.lib.validation.config.pypsa_eur import PypsaEurConfig from scripts.lib.validation.config.renewable import RenewableConfig from scripts.lib.validation.config.run import RunConfig @@ -140,6 +141,10 @@ class ConfigSchema(BaseModel): default_factory=RenewableConfig, description="Renewable energy technologies configuration.", ) + perennials: PerennialsConfig = Field( + default_factory=PerennialsConfig, + description="Perennialisation (carbon dioxide removal) configuration.", + ) conventional: ConventionalConfig = Field( default_factory=ConventionalConfig, description="Conventional power plants configuration.", diff --git a/scripts/lib/validation/config/data.py b/scripts/lib/validation/config/data.py index 18477518fd..bdc0dfddf4 100644 --- a/scripts/lib/validation/config/data.py +++ b/scripts/lib/validation/config/data.py @@ -201,6 +201,10 @@ def check_version_files_are_correct_suffix( default_factory=lambda: _DataSourceConfig(source="primary"), description="Instrat CO2 prices data source configuration.", ) + co2_removal_data: _DataSourceConfig = Field( + default_factory=lambda: _DataSourceConfig(source="primary"), + description="Carbon dioxide removal (afforestation, perennialisation) input data source configuration.", + ) co2stop: _DataSourceConfig = Field( default_factory=_DataSourceConfig, description="CO2Stop data source configuration.", diff --git a/scripts/lib/validation/config/perennials.py b/scripts/lib/validation/config/perennials.py new file mode 100644 index 0000000000..a7abc3f631 --- /dev/null +++ b/scripts/lib/validation/config/perennials.py @@ -0,0 +1,22 @@ +# SPDX-FileCopyrightText: Contributors to PyPSA-Eur +# +# SPDX-License-Identifier: MIT + +""" +Perennialisation configuration. + +See docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#perennials +""" + +from pydantic import Field + +from scripts.lib.validation.config._base import ConfigModel + + +class PerennialsConfig(ConfigModel): + """Configuration for `perennials` settings.""" + + sequestration_co2: float = Field( + 2, + description="Tonnes of CO2 equivalent sequestered per hectare per year when 1st-generation biofuel cropland is converted to perennial grasses.", + ) diff --git a/scripts/lib/validation/config/sector.py b/scripts/lib/validation/config/sector.py index fa561f2adf..df7370a8e5 100644 --- a/scripts/lib/validation/config/sector.py +++ b/scripts/lib/validation/config/sector.py @@ -669,6 +669,10 @@ class SectorConfig(BaseModel): False, description="Add option for coal CHPs with carbon capture." ) dac: bool = Field(True, description="Add option for Direct Air Capture (DAC).") + perennials: bool = Field( + False, + description="Add option for perennialisation (converting 1st-generation biofuel cropland to perennial grasses) as a carbon dioxide removal (CDR) technology.", + ) co2_vent: bool = Field( False, description="Add option for vent out CO2 from storages to the atmosphere.", diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index eec3efc6a7..ae1fc07d7b 100755 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -1322,6 +1322,124 @@ def add_dac(n, costs, spatial): ) +def add_perennials(n, costs): + """ + Add perennialisation (CDR via 1st-generation biofuel cropland conversion) + to the network as Bus, Store, and Link components. + + Perennial grasses sequester more soil carbon than the annual 1st-generation + (1G) biofuel crops (cereals, sugar beet, rapeseed) they replace. The land + area available for conversion at each node is backed out from the biomass + potential already allocated to 1G biofuels (``biomass_potentials``, + MWh/y) divided by the 1G crop yield (MWh/ha/y) at that node, giving a + displaced area in ha; multiplying by a fixed CO2 sequestration rate per + hectare (``perennials.sequestration_co2``) gives the store's CO2 potential. + A single "co2 perennials" Link models the harvesting process: CO2 drawn + from the atmosphere (bus0) is converted into biogas (bus3) and stored CO2 + (bus1), with capacity restricted to the April-October harvesting season + via ``p_max_pu``. + + Parameters + ---------- + n : pypsa.Network + The PyPSA network container object + costs : pd.DataFrame + Costs and parameters for different technologies. Must contain a + 'perennials gbr' entry with 'electricity-input', 'biogas-output', + 'capital_cost', 'VOM', and 'lifetime' parameters + + Returns + ------- + None + Modifies the network object in-place by adding the perennials Bus, + Store, and Link + + Notes + ----- + Reads ``snakemake.input.biomass_potentials`` and + ``snakemake.input.perennials_yields_1G_biofuels`` (NUTS2-derived crop + yields aggregated to clustered network regions, see + ``build_perennials_yields_eurostat_average.py`` and + ``build_perennials_yields.py``), and + ``snakemake.config["perennials"]["sequestration_co2"]``. + """ + + logger.info("Adding perennials.") + + # load resources + biomass_potentials = pd.read_csv(snakemake.input.biomass_potentials, index_col=0) + perennials_yields_1G_biofuels = pd.read_csv( + snakemake.input.perennials_yields_1G_biofuels + ).set_index("name") + + # calculate perennials potential based on the conversion on first generation biofuels for equal area + perennials_area = ( + biomass_potentials.filter(regex="biofuels_1G") + / perennials_yields_1G_biofuels.filter(regex="biofuels_1G") + ).sum(axis=1) + # (MWh/y) / (MWh / ha / y) = (ha) returns the area used by sum of the 3 biofuels_1G classes which can be assigned for perennials + perennials_potentials = ( + perennials_area * snakemake.config["perennials"]["sequestration_co2"] + ) # (tCO2seq) = (ha) * (tCO2 seq/ha) + + nodes = pop_layout.index + n.add("Carrier", "co2 perennials") + + n.add( + "Bus", + nodes, + suffix=" co2 perennials", + location=nodes, + carrier="co2 perennials", + unit="t_co2", + ) + + # calculate CO2 sequestration per tDM perennials + perennial_CO2_seq = ( + perennials_yields_1G_biofuels["perennials"] + / snakemake.config["perennials"]["sequestration_co2"] + ) # (tDM/tCO2 seq) + + # calculate biogas production based on harvesting time (in month) + df_harvest = pd.DataFrame(index=n.snapshots, columns=["harvest"]) + df_harvest["harvest"] = df_harvest.index.month.isin([4, 5, 6, 7, 8, 9, 10]).astype( + int + ) + p_max_pu = pd.DataFrame(index=n.snapshots, columns=nodes) + for node in nodes: + p_max_pu[node] = df_harvest["harvest"] + + n.add( + "Link", + nodes, + suffix=" perennials refining", + bus0="co2 atmosphere", + bus1=nodes + " co2 perennials", + bus2=nodes.values, + bus3=spatial.gas.biogas, + efficiency=1, + efficiency2=-costs.at["perennials gbr", "electricity-input"] + * perennial_CO2_seq, + efficiency3=costs.at["perennials gbr", "biogas-output"] * perennial_CO2_seq, + carrier="co2 perennials", + p_nom_extendable=True, + p_max_pu=p_max_pu, + capital_cost=costs.at["perennials gbr", "capital_cost"] * perennial_CO2_seq, + marginal_cost=costs.at["perennials gbr", "VOM"] * perennial_CO2_seq, + lifetime=costs.at["perennials gbr", "lifetime"], + ) + + n.add( + "Store", + nodes, + suffix=" CO2s perennials", + bus=nodes + " co2 perennials", + e_nom=perennials_potentials.values, + carrier="co2 perennials", + e_cyclic=False, + ) + + def add_co2limit(n, options, co2_totals_file, countries, nyears, limit): """ Add a global CO2 emissions constraint to the network. @@ -6634,6 +6752,9 @@ def add_import_options( if options["dac"]: add_dac(n, costs, spatial) + if options.get("perennials"): + add_perennials(n, costs) + if not options["electricity_transmission_grid"]: decentral(n)