add Biochar as carbon removal technology

config/config.default.yaml

@@ -419,6 +419,13 @@ renewable:
     eia_norm_year: false
     eia_correct_by_capacity: false
     eia_approximate_missing: false
+  # land-eligibility settings for alternative carbon dioxide removal technologies,
+  # feeding determine_availability_matrix.py the same way onwind/solar do above
+  biochar:
+    corine: [12, 13, 14] # CORINE Land Cover codes for non-irrigated arable land (potential biochar feedstock land)
+    natura: false
+    cutout: default
+    excluder_resolution: 250
 
 # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#conventional
 conventional:
@@ -837,6 +844,9 @@ sector:
   methanation: true
   coal_cc: false
   dac: true
+  biochar:
+    enable: false
+    heat_output: false
   co2_vent: false
   heat_vent:
     urban central: true
@@ -1070,6 +1080,11 @@ industry:
   reference_year: 2023
   oil_refining_emissions: 0.013
 
+biochar:
+  application_per_sqkm: 1500 # tonnes of biochar per km² (≈15 t/ha)
+  max_land_usage: 0.2 # maximum fraction of suitable CORINE area used for biochar
+  number_years: 25 # amortisation period (years) for biochar CO2 storage
+
 # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#costs
 costs:
   year: 2050

config/plotting.default.yaml

@@ -738,3 +738,5 @@ plotting:
     import NH3: '#e2ed74'
     import oil: '#93eda2'
     import methanol: '#87d0e6'
+    biochar heat: '#228B22'
+    co2 biochar: '#005000'

rules/build_sector.smk

@@ -920,6 +920,58 @@ rule build_biomass_potentials:
         scripts("build_biomass_potentials.py")
 
 
+rule determine_carbon_dioxide_removal_availability_matrix:
+    input:
+        corine=ancient(rules.retrieve_corine.output["tif_file"]),
+        regions=resources("regions_onshore_base_s_{clusters}.geojson"),
+        cutout=lambda w: input_cutout(
+            w, config_provider("renewable", w.technology, "cutout")(w)
+        ),
+    output:
+        resources(
+            "availability_matrix_carbon_dioxide_removal_{clusters}_{technology}.nc"
+        ),
+    log:
+        logs(
+            "determine_carbon_dioxide_removal_availability_matrix_{clusters}_{technology}.log"
+        ),
+    wildcard_constraints:
+        technology="biochar",
+    threads: config["atlite"].get("nprocesses", 4)
+    resources:
+        mem_mb=config["atlite"].get("nprocesses", 4) * 5000,
+    params:
+        renewable=config_provider("renewable"),
+        plot_availability_matrix=config_provider("atlite", "plot_availability_matrix"),
+    message:
+        "Determining availability matrix for {wildcards.clusters} clusters and {wildcards.technology} carbon dioxide removal technology"
+    script:
+        scripts("determine_availability_matrix.py")
+
+
+rule build_available_land:
+    input:
+        availability_matrix=resources(
+            "availability_matrix_carbon_dioxide_removal_{clusters}_{technology}.nc"
+        ),
+        regions=resources("regions_onshore_base_s_{clusters}.geojson"),
+        cutout=lambda w: input_cutout(
+            w, config_provider("renewable", w.technology, "cutout")(w)
+        ),
+    output:
+        csv_file=resources("{technology}_available_land_s_{clusters}.csv"),
+    log:
+        logs("build_available_land_{clusters}_{technology}.log"),
+    wildcard_constraints:
+        technology="biochar",
+    resources:
+        mem_mb=8000,
+    params:
+        renewable=config_provider("renewable"),
+    script:
+        scripts("build_available_land.py")
+
+
 rule build_biomass_transport_costs:
     input:
         sc1="data/biomass_transport_costs_supplychain1.csv",
@@ -1689,6 +1741,11 @@ rule prepare_sector_network:
             if config_provider("sector", "district_heating", "ates", "enable")(w)
             else []
         ),
+        biochar_potentials=lambda w: (
+            resources("biochar_available_land_s_{clusters}.csv")
+            if config_provider("sector", "biochar", "enable")(w)
+            else []
+        ),
     output:
         resources(
             "networks/base_s_{clusters}_{opts}_{sector_opts}_{planning_horizons}.nc"

scripts/build_available_land.py

@@ -0,0 +1,71 @@
+# SPDX-FileCopyrightText: Contributors to PyPSA-Eur <https://github.com/pypsa/pypsa-eur>
+#
+# SPDX-License-Identifier: MIT
+"""
+Aggregate a per-bus land-eligibility availability matrix (as produced by
+``determine_availability_matrix.py`` / ``determine_rock_weathering_availability_matrix.py``)
+into eligible land area per network node, in km².
+
+Shared by biochar, afforestation and rock weathering: the per-bus fractional
+eligibility (0-1 per cutout grid cell) is multiplied by the cutout grid cell
+area and summed per bus, recovering the same eligible-area-in-km² number the
+previous per-technology scripts (``build_corine_potentials.py``,
+``build_EW_potentials.py``) each computed independently via
+``atlite.gis.shape_availability``.
+
+Output: CSV with columns ``area [sqkm]`` (total node area) and
+``potential [sqkm]`` (eligible area after exclusions), indexed by node name.
+The downstream consumers (``add_biochar``, ``build_afforestation_potentials``,
+``add_rock_weathering``) apply their own technology-specific rate (t/km²,
+growth rate, etc.) to the ``potential [sqkm]`` column.
+"""
+
+import logging
+
+import geopandas as gpd
+import pandas as pd
+import xarray as xr
+
+from scripts._helpers import configure_logging, load_cutout, set_scenario_config
+
+logger = logging.getLogger(__name__)
+
+
+if __name__ == "__main__":
+    if "snakemake" not in globals():
+        from scripts._helpers import mock_snakemake
+
+        snakemake = mock_snakemake(
+            "build_available_land", clusters="39", technology="biochar"
+        )
+    configure_logging(snakemake)
+    set_scenario_config(snakemake)
+
+    technology = snakemake.wildcards.technology
+
+    availability = xr.open_dataarray(snakemake.input.availability_matrix)
+    cutout = load_cutout(snakemake.input.cutout)
+    regions = (
+        gpd.read_file(snakemake.input.regions).set_index("name").rename_axis("bus")
+    )
+
+    cell_area = cutout.grid.to_crs(3035).area / 1e6
+    cell_area = xr.DataArray(
+        cell_area.values.reshape(cutout.shape),
+        [cutout.coords["y"], cutout.coords["x"]],
+    )
+
+    potential = (availability * cell_area).sum(dim=["x", "y"]).to_pandas()
+    area = regions.geometry.to_crs(3035).area / 1e6
+
+    df = pd.DataFrame({"area [sqkm]": area, "potential [sqkm]": potential})
+    df.index.name = "node"
+
+    logger.info(
+        "Total %s potential: %.0f sqkm (of %.0f sqkm total node area)",
+        technology,
+        df["potential [sqkm]"].sum(),
+        df["area [sqkm]"].sum(),
+    )
+
+    df.to_csv(snakemake.output.csv_file)

scripts/prepare_sector_network.py

@@ -1247,6 +1247,150 @@ def add_dac(n, costs):
     )
 
 
+def add_biochar(n, costs):
+    """
+    Add biochar (CDR via biomass pyrolysis and stable-carbon soil storage)
+    to the network as Bus, Store, and Link components, all sharing a single
+    "co2 biochar" carrier; optionally exports pyrolysis waste heat to urban
+    central heating.
+
+    A single pyrolysis Link per node converts solid biomass (bus2) and
+    electricity (bus3) into stored biochar-carbon (bus1, drawn from the
+    atmosphere via bus0) and, if enabled, district heat (bus4, with an
+    overflow Store for heat that cannot be absorbed by the local heat
+    network). The store's capacity is dimensioned from the eligible CORINE
+    land area (see ``determine_availability_matrix.py`` and
+    ``build_available_land.py``) multiplied by a per-km2 application rate
+    and the stable-carbon fraction (``co2_per_tonne``), amortised over an
+    assumed storage horizon.
+
+    Parameters
+    ----------
+    n : pypsa.Network
+        The PyPSA network container object
+    costs : pd.DataFrame
+        Costs and parameters for different technologies. Must contain a
+        'biochar pyrolysis' entry with 'biomass-input', 'yield-biochar',
+        'electricity-input', 'heat-output', 'capital_cost', and 'VOM'
+        parameters
+
+    Returns
+    -------
+    None
+        Modifies the network object in-place by adding the biochar Bus,
+        Store, and Link (and, if enabled, heat-output Bus/Link/Store)
+
+    Notes
+    -----
+    Reads ``snakemake.input.biochar_potentials`` (eligible area in km2 per
+    node) and ``snakemake.config["biochar"]["application_per_sqkm"]`` /
+    ``["max_land_usage"]`` / ``["number_years"]``. Heat output is gated by
+    ``sector.heating`` and ``sector.biochar.heat_output``.
+    """
+    logger.info("Adding biochar.")
+
+    biochar_potentials = pd.read_csv(snakemake.input.biochar_potentials).set_index(
+        "node"
+    )
+
+    n.add("Carrier", "co2 biochar")
+
+    n.add(
+        "Bus",
+        spatial.nodes + " co2 biochar",
+        location=spatial.nodes,
+        carrier="co2 biochar",
+        unit="t_co2",
+    )
+
+    co2_per_tonne = (
+        1
+        / costs.at["biochar pyrolysis", "biomass-input"]
+        * 1
+        / costs.at["biochar pyrolysis", "yield-biochar"]
+    )  # tCO2 / t_biochar
+
+    n.add(
+        "Store",
+        spatial.nodes + " co2 biochar",
+        bus=spatial.nodes + " co2 biochar",
+        carrier="co2 biochar",
+        e_nom_extendable=True,
+        e_nom_max=(
+            biochar_potentials["potential [sqkm]"].values
+            * co2_per_tonne
+            * snakemake.config["biochar"]["application_per_sqkm"]
+            * snakemake.config["biochar"]["max_land_usage"]
+            / snakemake.config["biochar"]["number_years"]
+        ),
+    )
+
+    if len(spatial.biomass.nodes) == 1:
+        biomass_buses = spatial.biomass.nodes
+    else:
+        biomass_buses = spatial.nodes + " solid biomass"
+
+    # heat output into urban central heat system
+    if (
+        snakemake.config["sector"]["heating"]
+        and snakemake.config["sector"]["biochar"]["heat_output"]
+    ):
+        logger.info("Adding biochar heat output to urban central heat system.")
+        biochar_heat_buses = []
+        for node in spatial.nodes:
+            if (node + " urban central heat") in n.buses.index:
+                biochar_heat_bus = node + " biochar heat"
+                biochar_heat_buses.append(biochar_heat_bus)
+                n.add("Bus", biochar_heat_bus, carrier="biochar heat")
+                n.add(
+                    "Link",
+                    biochar_heat_bus,
+                    bus0=biochar_heat_bus,
+                    bus1=node + " urban central heat",
+                    p_nom_extendable=True,
+                    carrier="biochar heat",
+                )
+                biochar_heat_bus_waste = node + " biochar heat waste"
+                n.add("Bus", biochar_heat_bus_waste, carrier="biochar heat")
+                n.add(
+                    "Store",
+                    biochar_heat_bus_waste,
+                    bus=biochar_heat_bus_waste,
+                    e_nom_extendable=True,
+                    carrier="biochar heat",
+                )
+                n.add(
+                    "Link",
+                    biochar_heat_bus_waste,
+                    bus0=biochar_heat_bus,
+                    bus1=biochar_heat_bus_waste,
+                    p_nom_extendable=True,
+                    carrier="biochar heat",
+                )
+            else:
+                biochar_heat_buses.append(None)
+    else:
+        biochar_heat_buses = [None]
+
+    n.add(
+        "Link",
+        spatial.nodes + " biochar",
+        bus0="co2 atmosphere",
+        bus1=spatial.nodes + " co2 biochar",
+        bus2=biomass_buses,
+        bus3=spatial.nodes,
+        bus4=biochar_heat_buses,
+        carrier="co2 biochar",
+        capital_cost=costs.at["biochar pyrolysis", "capital_cost"],
+        marginal_cost=costs.at["biochar pyrolysis", "VOM"],
+        efficiency=1.0,
+        efficiency2=-costs.at["biochar pyrolysis", "biomass-input"],
+        efficiency3=-costs.at["biochar pyrolysis", "electricity-input"],
+        efficiency4=costs.at["biochar pyrolysis", "heat-output"],
+        p_nom_extendable=True,
+    )
+
+
 def add_co2limit(n, options, co2_totals_file, countries, nyears, limit):
     """
     Add a global CO2 emissions constraint to the network.
@@ -6518,6 +6662,9 @@ def add_import_options(
     if options["dac"]:
         add_dac(n, costs)
 
+    if options.get("biochar", {}).get("enable"):
+        add_biochar(n, costs)
+
     if not options["electricity_transmission_grid"]:
         decentral(n)