Merge pull request #73 from PyPSA/biomass-transport

Biomass transport

Snakefile

@@ -1,4 +1,7 @@
 
+from snakemake.remote.HTTP import RemoteProvider as HTTPRemoteProvider
+HTTP = HTTPRemoteProvider()
+
 configfile: "config.yaml"
 
 
@@ -180,6 +183,21 @@ rule build_biomass_potentials:
     script: 'scripts/build_biomass_potentials.py'
 
 
+if config["sector"]["biomass_transport"]:
+    rule build_biomass_transport_costs:
+        input:
+            transport_cost_data=HTTP.remote("publications.jrc.ec.europa.eu/repository/bitstream/JRC98626/biomass potentials in europe_web rev.pdf", keep_local=True)
+        output:
+            biomass_transport_costs="resources/biomass_transport_costs.csv",
+        threads: 1
+        resources: mem_mb=1000
+        benchmark: "benchmarks/build_biomass_transport_costs"
+        script: 'scripts/build_biomass_transport_costs.py'
+    build_biomass_transport_costs_output = rules.build_biomass_transport_costs.output
+else:
+    build_biomass_transport_costs_output = {}
+
+
 rule build_ammonia_production:
     input:
         usgs="data/myb1-2017-nitro.xls"
@@ -321,8 +339,8 @@ rule prepare_sector_network:
         energy_totals_name='resources/energy_totals.csv',
         co2_totals_name='resources/co2_totals.csv',
         transport_name='resources/transport_data.csv',
-        traffic_data_KFZ = "data/emobility/KFZ__count",
+        traffic_data_KFZ="data/emobility/KFZ__count",
-        traffic_data_Pkw = "data/emobility/Pkw__count",
+        traffic_data_Pkw="data/emobility/Pkw__count",
         biomass_potentials='resources/biomass_potentials.csv',
         heat_profile="data/heat_load_profile_BDEW.csv",
         costs=CDIR + "costs_{planning_horizons}.csv",
@@ -352,7 +370,8 @@ rule prepare_sector_network:
         solar_thermal_total="resources/solar_thermal_total_elec_s{simpl}_{clusters}.nc",
         solar_thermal_urban="resources/solar_thermal_urban_elec_s{simpl}_{clusters}.nc",
         solar_thermal_rural="resources/solar_thermal_rural_elec_s{simpl}_{clusters}.nc",
-	    **build_retro_cost_output
+        **build_retro_cost_output,
+        **build_biomass_transport_costs_output
     output: RDIR + '/prenetworks/elec_s{simpl}_{clusters}_lv{lv}_{opts}_{sector_opts}_{planning_horizons}.nc'
     threads: 1
     resources: mem_mb=2000

config.default.yaml

@@ -234,6 +234,7 @@ sector:
   electricity_grid_connection: true  # only applies to onshore wind and utility PV
   gas_distribution_grid: true
   gas_distribution_grid_cost_factor: 1.0  #multiplies cost in data/costs.csv
+  biomass_transport: false  # biomass transport between nodes
   conventional_generation: # generator : carrier
     OCGT: gas
 
@@ -275,7 +276,6 @@ industry:
   hotmaps_locate_missing: false
   reference_year: 2015
 
-
 costs:
   lifetime: 25 #default lifetime
   # From a Lion Hirth paper, also reflects average of Noothout et al 2016
@@ -509,5 +509,6 @@ plotting:
     shipping oil: "#6495ED"
     shipping oil emissions: "#6495ED"
     electricity distribution grid: '#333333'
+    solid biomass transport: green
     H2 for industry: "#222222"
     H2 for shipping: "#6495ED"

doc/release_notes.rst

@@ -76,6 +76,12 @@ Future release
 * The share of shipping transformed into hydrogen fuel cell can be now defined for different years in the ``config.yaml`` file. The carbon emission from the remaining share is treated as a negative load on the atmospheric carbon dioxide bus, just like aviation and land transport emissions.
 * The transformation of the Steel and Aluminium production can be now defined for different years in the ``config.yaml`` file.
 * Include the option to alter the maximum energy capacity of a store via the ``carrier+factor`` in the ``{sector_opts}`` wildcard. This can be useful for sensitivity analyses. Example: ``co2 stored+e2`` multiplies the ``e_nom_max`` by factor 2. In this example, ``e_nom_max`` represents the CO2 sequestration potential in Europe.
+* Add option to regionally disaggregate biomass potential to individual nodes
+  (currently given per country, then distributed by population density within)
+  and allow the transport of solid biomass.
+  The transport costs are determined based on the `JRC-EU-Times Bioenergy report <http://dx.doi.org/10.2790/01017>`_
+  in the new optional rule ``build_biomass_transport_costs``.
+  Biomass transport can be activated with the setting ``sector: biomass_transport: true``.
 * Compatibility with ``xarray`` version 0.19.
 
 PyPSA-Eur-Sec 0.5.0 (21st May 2021)

doc/spatial_resolution.rst

@@ -44,8 +44,10 @@ Hydrogen network: nodal.
 Methane network: single node for Europe, since future demand is so
 low and no bottlenecks are expected.
 
-Solid biomass:  single node for Europe, until transport costs can be
+Solid biomass: choice between single node for Europe and nodal where biomass
-incorporated.
+potential is regionally disaggregated (currently given per country,
+then distributed by population density within)
+and transport of solid biomass is possible.
 
 CO2:  single node for Europe, but a transport and storage cost is added for
 sequestered CO2. Optionally: nodal, with CO2 transport via pipelines.

scripts/build_biomass_transport_costs.py

@@ -0,0 +1,90 @@
+"""
+Reads biomass transport costs for different countries of the JRC report
+
+    "The JRC-EU-TIMES model.
+    Bioenergy potentials
+    for EU and neighbouring countries."
+    (2015)
+
+converts them from units 'EUR per km/ton' -> 'EUR/ (km MWh)'
+
+assuming as an approximation energy content of wood pellets
+
+@author: bw0928
+"""
+
+import pandas as pd
+import tabula as tbl
+
+ENERGY_CONTENT = 4.8  # unit MWh/t (wood pellets)
+
+def get_countries():
+
+    pandas_options = dict(
+        skiprows=range(6),
+        header=None,
+        index_col=0
+    )
+
+    return tbl.read_pdf(
+        str(snakemake.input.transport_cost_data),
+        pages="145",
+        multiple_tables=False,
+        pandas_options=pandas_options
+    )[0].index
+
+
+def get_cost_per_tkm(page, countries):
+    
+    pandas_options = dict(
+        skiprows=range(6),
+        header=0,
+        sep=' |,',
+        engine='python',
+        index_col=False,
+    )
+
+    sc = tbl.read_pdf(
+        str(snakemake.input.transport_cost_data),
+        pages=page,
+        multiple_tables=False,
+        pandas_options=pandas_options
+    )[0]
+    sc.index = countries
+    sc.columns = sc.columns.str.replace("€", "EUR")
+    
+    return sc
+
+
+def build_biomass_transport_costs():
+
+    countries = get_countries()
+
+    sc1 = get_cost_per_tkm(146, countries)
+    sc2 = get_cost_per_tkm(147, countries)
+
+    # take mean of both supply chains
+    to_concat = [sc1["EUR/km/ton"], sc2["EUR/km/ton"]]
+    transport_costs = pd.concat(to_concat, axis=1).mean(axis=1)
+
+    # convert tonnes to MWh
+    transport_costs /= ENERGY_CONTENT
+    transport_costs.name = "EUR/km/MWh"
+
+    # rename country names
+    to_rename = {
+        "UK": "GB",
+        "XK": "KO",
+        "EL": "GR"
+    }
+    transport_costs.rename(to_rename, inplace=True)
+
+    # add missing Norway with data from Sweden
+    transport_costs["NO"] = transport_costs["SE"]
+
+    transport_costs.to_csv(snakemake.output[0])
+
+
+if __name__ == "__main__":
+
+    build_biomass_transport_costs()

scripts/prepare_sector_network.py

@@ -31,12 +31,31 @@ def define_spatial(nodes):
     ----------
     nodes : list-like
     """
-    
+
     global spatial
     global options
-    
+
     spatial.nodes = nodes
     
+    # biomass
+
+    spatial.biomass = SimpleNamespace()
+
+    if options["biomass_transport"]:
+        spatial.biomass.nodes = nodes + " solid biomass"
+        spatial.biomass.locations = nodes
+        spatial.biomass.industry = nodes + " solid biomass for industry"
+        spatial.biomass.industry_cc = nodes + " solid biomass for industry CC"
+    else:
+        spatial.biomass.nodes = ["EU solid biomass"]
+        spatial.biomass.locations = ["EU"]
+        spatial.biomass.industry = ["solid biomass for industry"]
+        spatial.biomass.industry_cc = ["solid biomass for industry CC"]
+
+    spatial.biomass.df = pd.DataFrame(vars(spatial.biomass), index=nodes)
+    
+    # co2
+
     spatial.co2 = SimpleNamespace()
     
     if options["co2_network"]:
@@ -53,7 +72,7 @@ def define_spatial(nodes):
 
 def emission_sectors_from_opts(opts):
 
-    sectors = ["electricity"]
+    sectors = ["electricity"]   
     if "T" in opts:
         sectors += [
             "rail non-elec",
@@ -206,6 +225,53 @@ def add_lifetime_wind_solar(n, costs):
         n.generators.loc[gen_i, "lifetime"] = costs.at[carrier, 'lifetime']
 
 
+def create_network_topology(n, prefix, connector=" -> ", bidirectional=True):
+    """
+    Create a network topology like the power transmission network.
+    
+    Parameters
+    ----------
+    n : pypsa.Network
+    prefix : str
+    connector : str
+    bidirectional : bool, default True
+        True: one link for each connection
+        False: one link for each connection and direction (back and forth)
+        
+    Returns
+    -------
+    pd.DataFrame with columns bus0, bus1 and length
+    """
+
+    ln_attrs = ["bus0", "bus1", "length"]
+    lk_attrs = ["bus0", "bus1", "length", "underwater_fraction"]
+
+    candidates = pd.concat([
+        n.lines[ln_attrs],
+        n.links.loc[n.links.carrier == "DC", lk_attrs]
+    ]).fillna(0)
+
+    positive_order = candidates.bus0 < candidates.bus1
+    candidates_p = candidates[positive_order]
+    swap_buses = {"bus0": "bus1", "bus1": "bus0"}
+    candidates_n = candidates[~positive_order].rename(columns=swap_buses)
+    candidates = pd.concat([candidates_p, candidates_n])
+    
+    def make_index(c):
+        return prefix + c.bus0 + connector + c.bus1
+
+    topo = candidates.groupby(["bus0", "bus1"], as_index=False).mean()
+    topo.index = topo.apply(make_index, axis=1)
+    
+    if not bidirectional:
+        topo_reverse = topo.copy()
+        topo_reverse.rename(columns=swap_buses, inplace=True)
+        topo_reverse.index = topo_reverse.apply(make_index, axis=1)
+        topo = topo.append(topo_reverse)
+    
+    return topo
+
+
 # TODO merge issue with PyPSA-Eur
 def update_wind_solar_costs(n, costs):
     """
@@ -1623,8 +1689,16 @@ def add_biomass(n, costs):
 
     biomass_potentials = pd.read_csv(snakemake.input.biomass_potentials, index_col=0)
 
-    n.add("Carrier", "biogas")
+    if options["biomass_transport"]:
+        # potential per node distributed within country by population
+        biomass_potentials_spatial = (biomass_potentials.loc[pop_layout.ct]
+                            .set_index(pop_layout.index)
+                            .mul(pop_layout.fraction, axis="index")
+                            .rename(index=lambda x: x + " solid biomass"))
+    else:
+        biomass_potentials_spatial = biomass_potentials.sum()
 
+    n.add("Carrier", "biogas")
     n.add("Carrier", "solid biomass")
 
     n.add("Bus",
@@ -1633,9 +1707,9 @@ def add_biomass(n, costs):
         carrier="biogas"
     )
 
-    n.add("Bus",
+    n.madd("Bus",
-        "EU solid biomass",
+        spatial.biomass.nodes,
-        location="EU",
+        location=spatial.biomass.locations,
         carrier="solid biomass"
     )
 
@@ -1648,13 +1722,13 @@ def add_biomass(n, costs):
         e_initial=biomass_potentials.loc[countries, "biogas"].sum()
     )
 
-    n.add("Store",
+    n.madd("Store",
-        "EU solid biomass",
+        spatial.biomass.nodes,
-        bus="EU solid biomass",
+        bus=spatial.biomass.nodes,
         carrier="solid biomass",
-        e_nom=biomass_potentials.loc[countries, "solid biomass"].sum(),
+        e_nom=biomass_potentials_spatial["solid biomass"],
         marginal_cost=costs.at['solid biomass', 'fuel'],
-        e_initial=biomass_potentials.loc[countries, "solid biomass"].sum()
+        e_initial=biomass_potentials_spatial["solid biomass"]
     )
 
     n.add("Link",
@@ -1669,6 +1743,32 @@ def add_biomass(n, costs):
         p_nom_extendable=True
     )
 
+    if options["biomass_transport"]:
+
+        transport_costs = pd.read_csv(
+            snakemake.input.biomass_transport_costs,
+            index_col=0,
+            squeeze=True
+        )
+        
+        # add biomass transport
+        biomass_transport = create_network_topology(n, "biomass transport ", bidirectional=False)
+
+        # costs
+        bus0_costs = biomass_transport.bus0.apply(lambda x: transport_costs[x[:2]])
+        bus1_costs = biomass_transport.bus1.apply(lambda x: transport_costs[x[:2]])
+        biomass_transport["costs"] = pd.concat([bus0_costs, bus1_costs], axis=1).mean(axis=1)
+
+        n.madd("Link",
+            biomass_transport.index,
+            bus0=biomass_transport.bus0 + " solid biomass",
+            bus1=biomass_transport.bus1 + " solid biomass",
+            p_nom_extendable=True,
+            length=biomass_transport.length.values,
+            marginal_cost=biomass_transport.costs * biomass_transport.length.values,
+            capital_cost=1,
+            carrier="solid biomass transport"
+        )
 
     #AC buses with district heating
     urban_central = n.buses.index[n.buses.carrier == "urban central heat"]
@@ -1679,7 +1779,7 @@ def add_biomass(n, costs):
 
         n.madd("Link",
             urban_central + " urban central solid biomass CHP",
-            bus0="EU solid biomass",
+            bus0=spatial.biomass.df.loc[urban_central, "nodes"].values,
             bus1=urban_central,
             bus2=urban_central + " urban central heat",
             carrier="urban central solid biomass CHP",
@@ -1693,7 +1793,7 @@ def add_biomass(n, costs):
 
         n.madd("Link",
             urban_central + " urban central solid biomass CHP CC",
-            bus0="EU solid biomass",
+            bus0=spatial.biomass.df.loc[urban_central, "nodes"].values,
             bus1=urban_central,
             bus2=urban_central + " urban central heat",
             bus3="co2 atmosphere",
@@ -1719,35 +1819,37 @@ def add_industry(n, costs):
     # 1e6 to convert TWh to MWh
     industrial_demand = pd.read_csv(snakemake.input.industrial_demand, index_col=0) * 1e6
 
-    solid_biomass_by_country = industrial_demand["solid biomass"].groupby(pop_layout.ct).sum()
+    n.madd("Bus",
-
+        spatial.biomass.industry,
-    n.add("Bus",
+        location=spatial.biomass.locations,
-        "solid biomass for industry",
-        location="EU",
         carrier="solid biomass for industry"
     )
 
-    n.add("Load",
+    if options["biomass_transport"]:
-        "solid biomass for industry",
+        p_set = industrial_demand.loc[spatial.biomass.locations, "solid biomass"].rename(index=lambda x: x + " solid biomass for industry") / 8760
-        bus="solid biomass for industry",
+    else:
+        p_set = industrial_demand["solid biomass"].sum() / 8760
+
+    n.madd("Load",
+        spatial.biomass.industry,
+        bus=spatial.biomass.industry,
         carrier="solid biomass for industry",
-        p_set=solid_biomass_by_country.sum() / 8760
+        p_set=p_set
     )
 
-    n.add("Link",
+    n.madd("Link",
-        "solid biomass for industry",
+        spatial.biomass.industry,
-        bus0="EU solid biomass",
+        bus0=spatial.biomass.nodes,
-        bus1="solid biomass for industry",
+        bus1=spatial.biomass.industry,
         carrier="solid biomass for industry",
         p_nom_extendable=True,
         efficiency=1.
     )
 
     n.madd("Link",
-        spatial.co2.locations,
+        spatial.biomass.industry_cc,
-        suffix=" solid biomass for industry CC",
+        bus0=spatial.biomass.nodes,
-        bus0="EU solid biomass",
+        bus1=spatial.biomass.industry,
-        bus1="solid biomass for industry",
         bus2="co2 atmosphere",
         bus3=spatial.co2.nodes,
         carrier="solid biomass for industry CC",
@@ -2156,6 +2258,8 @@ if __name__ == "__main__":
         if o[:4] == "dist":
             options['electricity_distribution_grid'] = True
             options['electricity_distribution_grid_cost_factor'] = float(o[4:].replace("p", ".").replace("m", "-"))
+        if o == "biomasstransport":
+            options["biomass_transport"] = True
 
     nodal_energy_totals, heat_demand, ashp_cop, gshp_cop, solar_thermal, transport, avail_profile, dsm_profile, nodal_transport_data = prepare_data(n)