Merge branch 'master' into improve-doc

.github/workflows/ci.yaml

@@ -64,7 +64,7 @@ jobs:
 
       - name: Add solver to environment
         run: |
-          echo -e "  - coincbc\n  - ipopt<3.13.3" >> ../pypsa-eur/envs/environment.yaml
+          echo -e "- coincbc\n- ipopt<3.13.3" >> ../pypsa-eur/envs/environment.yaml
 
       - name: Setup Mambaforge
         uses: conda-incubator/setup-miniconda@v2

.gitignore

@@ -46,3 +46,5 @@ config.yaml
 doc/_build
 
 *.xls
+
+*.geojson

Snakefile

@@ -256,9 +256,9 @@ rule build_biomass_potentials:
         enspreso_biomass=HTTP.remote("https://cidportal.jrc.ec.europa.eu/ftp/jrc-opendata/ENSPRESO/ENSPRESO_BIOMASS.xlsx", keep_local=True),
         nuts2="data/nuts/NUTS_RG_10M_2013_4326_LEVL_2.geojson", # https://gisco-services.ec.europa.eu/distribution/v2/nuts/download/#nuts21
         regions_onshore=pypsaeur("resources/regions_onshore_elec_s{simpl}_{clusters}.geojson"),
-        nuts3_population="../pypsa-eur/data/bundle/nama_10r_3popgdp.tsv.gz",
+        nuts3_population=pypsaeur("data/bundle/nama_10r_3popgdp.tsv.gz"),
-        swiss_cantons="../pypsa-eur/data/bundle/ch_cantons.csv",
+        swiss_cantons=pypsaeur("data/bundle/ch_cantons.csv"),
-        swiss_population="../pypsa-eur/data/bundle/je-e-21.03.02.xls",
+        swiss_population=pypsaeur("data/bundle/je-e-21.03.02.xls"),
         country_shapes=pypsaeur('resources/country_shapes.geojson')
     output:
         biomass_potentials_all='resources/biomass_potentials_all_s{simpl}_{clusters}.csv',
@@ -474,7 +474,7 @@ rule prepare_sector_network:
         co2="data/eea/UNFCCC_v23.csv",
         biomass_potentials='resources/biomass_potentials_s{simpl}_{clusters}.csv',
         heat_profile="data/heat_load_profile_BDEW.csv",
-        costs=CDIR + "costs_{planning_horizons}.csv",
+        costs=CDIR + "costs_{}.csv".format(config['costs']['year']) if config["foresight"] == "overnight" else CDIR + "costs_{planning_horizons}.csv",
         profile_offwind_ac=pypsaeur("resources/profile_offwind-ac.nc"),
         profile_offwind_dc=pypsaeur("resources/profile_offwind-dc.nc"),
         h2_cavern="resources/salt_cavern_potentials_s{simpl}_{clusters}.csv",
@@ -532,6 +532,14 @@ rule copy_config:
     script: "scripts/copy_config.py"
 
 
+rule copy_conda_env:
+    output: SDIR + '/configs/environment.yaml'
+    threads: 1
+    resources: mem_mb=500
+    benchmark: SDIR + "/benchmarks/copy_conda_env"
+    shell: "conda env export -f {output} --no-builds"
+
+
 rule make_summary:
     input:
         overrides="data/override_component_attrs",
@@ -539,7 +547,7 @@ rule make_summary:
             RDIR + "/postnetworks/elec_s{simpl}_{clusters}_lv{lv}_{opts}_{sector_opts}_{planning_horizons}.nc",
             **config['scenario']
         ),
-        costs=CDIR + "costs_{}.csv".format(config['scenario']['planning_horizons'][0]),
+        costs=CDIR + "costs_{}.csv".format(config['costs']['year']) if config["foresight"] == "overnight" else CDIR + "costs_{}.csv".format(config['scenario']['planning_horizons'][0]),
         plots=expand(
             RDIR + "/maps/elec_s{simpl}_{clusters}_lv{lv}_{opts}_{sector_opts}-costs-all_{planning_horizons}.pdf",
             **config['scenario']
@@ -589,8 +597,9 @@ if config["foresight"] == "overnight":
         input:
             overrides="data/override_component_attrs",
             network=RDIR + "/prenetworks/elec_s{simpl}_{clusters}_lv{lv}_{opts}_{sector_opts}_{planning_horizons}.nc",
-            costs=CDIR + "costs_{planning_horizons}.csv",
+            costs=CDIR + "costs_{}.csv".format(config['costs']['year']),
-            config=SDIR + '/configs/config.yaml'
+            config=SDIR + '/configs/config.yaml',
+            env=SDIR + '/configs/environment.yaml',
         output: RDIR + "/postnetworks/elec_s{simpl}_{clusters}_lv{lv}_{opts}_{sector_opts}_{planning_horizons}.nc"
         shadow: "shallow"
         log:

config.default.yaml

@@ -33,15 +33,15 @@ scenario:
   # A for agriculture, forestry and fishing
   # solar+c0.5 reduces the capital cost of solar to 50\% of reference value
   # solar+p3 multiplies the available installable potential by factor 3
-  # co2 stored+e2 multiplies the potential of CO2 sequestration by a factor 2
+  # seq400 sets the potential of CO2 sequestration to 400 Mt CO2 per year
   # dist{n} includes distribution grids with investment cost of n times cost in data/costs.csv
   # for myopic/perfect foresight cb states the carbon budget in GtCO2 (cumulative
   # emissions throughout the transition path in the timeframe determined by the
   # planning_horizons), be:beta decay; ex:exponential decay
   # cb40ex0 distributes a carbon budget of 40 GtCO2 following an exponential
   # decay with initial growth rate 0
-  planning_horizons: # investment years for myopic and perfect; or costs year for overnight
+  planning_horizons: # investment years for myopic and perfect; for overnight, year of cost assumptions can be different and is defined under 'costs'
-    - 2030
+    - 2050
   # for example, set to
   # - 2020
   # - 2030
@@ -154,11 +154,11 @@ sector:
     potential: 0.6  # maximum fraction of urban demand which can be supplied by district heating
      # increase of today's district heating demand to potential maximum district heating share
      # progress = 0 means today's district heating share, progress = 1 means maximum fraction of urban demand is supplied by district heating
-    progress: 1
+    progress:
-      # 2020: 0.0
+      2020: 0.0
-      # 2030: 0.3
+      2030: 0.3
-      # 2040: 0.6
+      2040: 0.6
-      # 2050: 1.0
+      2050: 1.0
     district_heating_loss: 0.15
   bev_dsm_restriction_value: 0.75  #Set to 0 for no restriction on BEV DSM
   bev_dsm_restriction_time: 7  #Time at which SOC of BEV has to be dsm_restriction_value
@@ -178,16 +178,16 @@ sector:
   bev_avail_mean: 0.8
   v2g: true #allows feed-in to grid from EV battery
   #what is not EV or FCEV is oil-fuelled ICE
-  land_transport_fuel_cell_share: 0.15 # 1 means all FCEVs
+  land_transport_fuel_cell_share:  # 1 means all FCEVs
-    # 2020: 0
+    2020: 0
-    # 2030: 0.05
+    2030: 0.05
-    # 2040: 0.1
+    2040: 0.1
-    # 2050: 0.15
+    2050: 0.15
-  land_transport_electric_share: 0.85 # 1 means all EVs
+  land_transport_electric_share:  # 1 means all EVs
-    # 2020: 0
+    2020: 0
-    # 2030: 0.25
+    2030: 0.25
-    # 2040: 0.6
+    2040: 0.6
-    # 2050: 0.85
+    2050: 0.85
   transport_fuel_cell_efficiency: 0.5
   transport_internal_combustion_efficiency: 0.3
   agriculture_machinery_electric_share: 0
@@ -195,29 +195,22 @@ sector:
   agriculture_machinery_electric_efficiency: 0.3 # electricity per use
   shipping_average_efficiency: 0.4 #For conversion of fuel oil to propulsion in 2011
   shipping_hydrogen_liquefaction: false # whether to consider liquefaction costs for shipping H2 demands
-  shipping_hydrogen_share: 1 # 1 means all hydrogen FC
+  shipping_hydrogen_share:  0
-    # 2020: 0
-    # 2025: 0
-    # 2030: 0.05
-    # 2035: 0.15
-    # 2040: 0.3
-    # 2045: 0.6
-    # 2050: 1
   time_dep_hp_cop: true #time dependent heat pump coefficient of performance
   heat_pump_sink_T: 55. # Celsius, based on DTU / large area radiators; used in build_cop_profiles.py
    # conservatively high to cover hot water and space heating in poorly-insulated buildings
   reduce_space_heat_exogenously: true  # reduces space heat demand by a given factor (applied before losses in DH)
   # this can represent e.g. building renovation, building demolition, or if
   # the factor is negative: increasing floor area, increased thermal comfort, population growth
-  reduce_space_heat_exogenously_factor: 0.29 # per unit reduction in space heat demand
+  reduce_space_heat_exogenously_factor:  # per unit reduction in space heat demand
   # the default factors are determined by the LTS scenario from http://tool.european-calculator.eu/app/buildings/building-types-area/?levers=1ddd4444421213bdbbbddd44444ffffff11f411111221111211l212221
-    # 2020: 0.10  # this results in a space heat demand reduction of 10%
+    2020: 0.10  # this results in a space heat demand reduction of 10%
-    # 2025: 0.09  # first heat demand increases compared to 2020 because of larger floor area per capita
+    2025: 0.09  # first heat demand increases compared to 2020 because of larger floor area per capita
-    # 2030: 0.09
+    2030: 0.09
-    # 2035: 0.11
+    2035: 0.11
-    # 2040: 0.16
+    2040: 0.16
-    # 2045: 0.21
+    2045: 0.21
-    # 2050: 0.29
+    2050: 0.29
   retrofitting :  # co-optimises building renovation to reduce space heat demand
     retro_endogen: false  # co-optimise space heat savings
     cost_factor: 1.0   # weight costs for building renovation
@@ -252,6 +245,7 @@ sector:
     # - onshore  # more than 50 km from sea
     - nearshore  # within 50 km of sea
     # - offshore
+  ammonia: false # can be false (no NH3 carrier), true (copperplated NH3), "regional" (regionalised NH3 without network)
   use_fischer_tropsch_waste_heat: true
   use_fuel_cell_waste_heat: true
   electricity_distribution_grid: true
@@ -276,36 +270,38 @@ sector:
 
 
 industry:
-  St_primary_fraction: 0.3 # fraction of steel produced via primary route versus secondary route (scrap+EAF); today fraction is 0.6
+  St_primary_fraction:  # fraction of steel produced via primary route versus secondary route (scrap+EAF); today fraction is 0.6
-    # 2020: 0.6
+    2020: 0.6
-    # 2025: 0.55
+    2025: 0.55
-    # 2030: 0.5
+    2030: 0.5
-    # 2035: 0.45
+    2035: 0.45
-    # 2040: 0.4
+    2040: 0.4
-    # 2045: 0.35
+    2045: 0.35
-    # 2050: 0.3
+    2050: 0.3
-  DRI_fraction: 1 # fraction of the primary route converted to DRI + EAF
+  DRI_fraction:  # fraction of the primary route converted to DRI + EAF
-    # 2020: 0
+    2020: 0
-    # 2025: 0
+    2025: 0
-    # 2030: 0.05
+    2030: 0.05
-    # 2035: 0.2
+    2035: 0.2
-    # 2040: 0.4
+    2040: 0.4
-    # 2045: 0.7
+    2045: 0.7
-    # 2050: 1
+    2050: 1
   H2_DRI: 1.7   #H2 consumption in Direct Reduced Iron (DRI),  MWh_H2,LHV/ton_Steel from 51kgH2/tSt in Vogl et al (2018) doi:10.1016/j.jclepro.2018.08.279
   elec_DRI: 0.322   #electricity consumption in Direct Reduced Iron (DRI) shaft, MWh/tSt HYBRIT brochure https://ssabwebsitecdn.azureedge.net/-/media/hybrit/files/hybrit_brochure.pdf
-  Al_primary_fraction: 0.2 # fraction of aluminium produced via the primary route versus scrap; today fraction is 0.4
+  Al_primary_fraction:  # fraction of aluminium produced via the primary route versus scrap; today fraction is 0.4
-    # 2020: 0.4
+    2020: 0.4
-    # 2025: 0.375
+    2025: 0.375
-    # 2030: 0.35
+    2030: 0.35
-    # 2035: 0.325
+    2035: 0.325
-    # 2040: 0.3
+    2040: 0.3
-    # 2045: 0.25
+    2045: 0.25
-    # 2050: 0.2
+    2050: 0.2
+  MWh_NH3_per_tNH3: 5.166 # LHV
   MWh_CH4_per_tNH3_SMR: 10.8 # 2012's demand from https://ec.europa.eu/docsroom/documents/4165/attachments/1/translations/en/renditions/pdf
   MWh_elec_per_tNH3_SMR: 0.7 # same source, assuming 94-6% split methane-elec of total energy demand 11.5 MWh/tNH3
   MWh_H2_per_tNH3_electrolysis: 6.5 # from https://doi.org/10.1016/j.joule.2018.04.017, around 0.197 tH2/tHN3 (>3/17 since some H2 lost and used for energy)
   MWh_elec_per_tNH3_electrolysis: 1.17 # from https://doi.org/10.1016/j.joule.2018.04.017 Table 13 (air separation and HB)
+  MWh_NH3_per_MWh_H2_cracker: 1.46 # https://github.com/euronion/trace/blob/44a5ff8401762edbef80eff9cfe5a47c8d3c8be4/data/efficiencies.csv
   NH3_process_emissions: 24.5 # in MtCO2/a from SMR for H2 production for NH3 from UNFCCC for 2015 for EU28
   petrochemical_process_emissions: 25.5 # in MtCO2/a for petrochemical and other from UNFCCC for 2015 for EU28
   HVC_primary_fraction: 1. # fraction of today's HVC produced via primary route
@@ -327,6 +323,7 @@ industry:
   # Material Economics (2019): https://materialeconomics.com/latest-updates/industrial-transformation-2050
 
 costs:
+  year: 2030
   lifetime: 25 #default lifetime
   # From a Lion Hirth paper, also reflects average of Noothout et al 2016
   discountrate: 0.07
@@ -585,6 +582,12 @@ plotting:
     H2 pipeline retrofitted: '#ba99b5'
     H2 Fuel Cell: '#c251ae'
     H2 Electrolysis: '#ff29d9'
+    # ammonia
+    NH3: '#46caf0'
+    ammonia: '#46caf0'
+    ammonia store: '#00ace0'
+    ammonia cracker: '#87d0e6'
+    Haber-Bosch: '#076987'
     # syngas
     Sabatier: '#9850ad'
     methanation: '#c44ce6'

doc/data.csv

@@ -23,7 +23,7 @@ Floor area missing in hotmaps building stock data,floor_area_missing.csv,unknown
 Comparative level investment,comparative_level_investment.csv,Eurostat,https://ec.europa.eu/eurostat/statistics-explained/index.php?title=Comparative_price_levels_for_investment
 Electricity taxes,electricity_taxes_eu.csv,Eurostat,https://appsso.eurostat.ec.europa.eu/nui/show.do?dataset=nrg_pc_204&lang=en
 Building topologies and corresponding standard values,tabula-calculator-calcsetbuilding.csv,unknown,https://episcope.eu/fileadmin/tabula/public/calc/tabula-calculator.xlsx
-Retrofitting thermal envelope costs for Germany,retro_cost_germany.csv,unkown,https://www.iwu.de/forschung/handlungslogiken/kosten-energierelevanter-bau-und-anlagenteile-bei-modernisierung/
+Retrofitting thermal envelope costs for Germany,retro_cost_germany.csv,unknown,https://www.iwu.de/forschung/handlungslogiken/kosten-energierelevanter-bau-und-anlagenteile-bei-modernisierung/
 District heating most countries,jrc-idees-2015/,CC BY 4.0,https://ec.europa.eu/jrc/en/potencia/jrc-idees,,
-District heating missing countries,district_heat_share.csv,unkown,https://www.euroheat.org/knowledge-hub/country-profiles,,
+District heating missing countries,district_heat_share.csv,unknown,https://www.euroheat.org/knowledge-hub/country-profiles,,
 

doc/installation.rst

@@ -54,7 +54,7 @@ The requirements are the same as `PyPSA-Eur <https://github.com/PyPSA/pypsa-eur>
 xarray version >= 0.15.1, you will need the latest master branch of
 atlite version 0.0.2.
 
-You can create an enviroment using the environment.yaml file in pypsa-eur/envs:
+You can create an environment using the environment.yaml file in pypsa-eur/envs:
 
 .. code:: bash
 

doc/release_notes.rst

@@ -28,7 +28,7 @@ incorporates retrofitting options to hydrogen.
 
 * New rule ``cluster_gas_network`` that clusters the gas transmission network
   data to the model resolution. Cross-regional pipeline capacities are aggregated
-  (while pressure and diameter compability is ignored), intra-regional pipelines
+  (while pressure and diameter compatibility is ignored), intra-regional pipelines
   are dropped. Lengths are recalculated based on the regions' centroids.
 
 * With the option ``sector: gas_network:``, the existing gas network is
@@ -57,6 +57,8 @@ incorporates retrofitting options to hydrogen.
 **New features and functionality**
 
 
+* Add option to aggregate network temporally using representative snapshots or segments (with tsam package)
+
 * Add option for biomass boilers (wood pellets) for decentral heating
 
 * Add option for BioSNG (methane from biomass) with and without CC
@@ -71,10 +73,17 @@ incorporates retrofitting options to hydrogen.
 
 * Add option to sweep the global CO2 sequestration potentials with keyword ``seq200`` in the ``{sector_opts}`` wildcard (for limit of 200 Mt CO2).
 
+* Add option to resolve ammonia as separate energy carrier with Haber-Bosch
+  synthesis, ammonia cracking, storage and industrial demand. The ammonia
+  carrier can be nodally resolved or copperplated across Europe. This feature is
+  controlled by ``sector: ammonia:``.
+
 * Updated `data bundle <https://zenodo.org/record/5824485/files/pypsa-eur-sec-data-bundle.tar.gz>`_ that includes the hydrogan salt cavern storage potentials.
 
 * Updated and extended documentation in <https://pypsa-eur-sec.readthedocs.io/en/latest/>
 
+* Shipping demand now defaults to (synthetic) oil rather than liquefied hydrogen until 2050.
+
 **Bugfixes**
 
 * The CO2 sequestration limit implemented as GlobalConstraint (introduced in the previous version)

doc/spatial_resolution.rst

@@ -6,6 +6,8 @@ Spatial resolution
 
 The default nodal resolution of the model follows the electricity generation and transmission model `PyPSA-Eur <https://github.com/PyPSA/pypsa-eur>`_, which clusters down the electricity transmission substations in each European country based on the k-means algorithm (See `cluster_network <https://pypsa-eur.readthedocs.io/en/latest/simplification/cluster_network.html#rule-cluster-network>`_ for a complete explanation). This gives nodes which correspond to major load and generation centres (typically cities).
 
+The total number of nodes for Europe is set in the ``config.yaml`` file under ``clusters``. The number of nodes can vary between 37, the number of independent countries / synchronous areas, and several hundred. With 200-300 nodes the model needs 100-150 GB RAM to solve with a commercial solver like Gurobi.
+
 Exemplary unsolved network clustered to 512 nodes:
 
 .. image:: ../graphics/elec_s_512.png 

scripts/add_existing_baseyear.py

@@ -165,11 +165,11 @@ def add_power_capacities_installed_before_baseyear(n, grouping_years, costs, bas
     df_agg.loc[biomass_i, 'DateOut'] = df_agg.loc[biomass_i, 'DateOut'].fillna(dateout)
 
 
-    # drop assets which are already phased out / decomissioned
+    # drop assets which are already phased out / decommissioned
     phased_out = df_agg[df_agg["DateOut"]<baseyear].index
     df_agg.drop(phased_out, inplace=True)
 
-    # calculate remaining lifetime before phase-out (+1 because assumming
+    # calculate remaining lifetime before phase-out (+1 because assuming
     # phase out date at the end of the year)
     df_agg["lifetime"] = df_agg.DateOut - df_agg.DateIn + 1
 
@@ -251,7 +251,7 @@ def add_power_capacities_installed_before_baseyear(n, grouping_years, costs, bas
                     # existing capacities are split evenly among regions in every country
                     inv_ind = [i for i in inv_busmap[ind]]
 
-                    # for offshore the spliting only inludes coastal regions
+                    # for offshore the splitting only includes coastal regions
                     inv_ind = [i for i in inv_ind if (i + name_suffix) in n.generators.index]
 
                     p_max_pu = n.generators_t.p_max_pu[[i + name_suffix for i in inv_ind]]

scripts/build_industrial_energy_demand_per_country_today.py

@@ -65,6 +65,8 @@ def industrial_energy_demand_per_country(country):
 
         df = df_dict[sheet][year].groupby(fuels).sum()
 
+        df["ammonia"] = 0.
+
         df['other'] = df['all'] - df.loc[df.index != 'all'].sum()
 
         return df
@@ -89,18 +91,21 @@ def add_ammonia_energy_demand(demand):
     fn = snakemake.input.ammonia_production
     ammonia = pd.read_csv(fn, index_col=0)[str(year)] / 1e3
 
-    def ammonia_by_fuel(x):
+    def get_ammonia_by_fuel(x):
 
         fuels = {'gas': config['MWh_CH4_per_tNH3_SMR'],
                  'electricity': config['MWh_elec_per_tNH3_SMR']}
 
         return pd.Series({k: x*v for k,v in fuels.items()})
 
-    ammonia = ammonia.apply(ammonia_by_fuel).T
+    ammonia_by_fuel = ammonia.apply(get_ammonia_by_fuel).T
+    ammonia_by_fuel = ammonia_by_fuel.unstack().reindex(index=demand.index, fill_value=0.)
+
+    ammonia = pd.DataFrame({"ammonia": ammonia * config['MWh_NH3_per_tNH3']}).T
 
     demand['Ammonia'] = ammonia.unstack().reindex(index=demand.index, fill_value=0.)
 
-    demand['Basic chemicals (without ammonia)'] = demand["Basic chemicals"] - demand["Ammonia"]
+    demand['Basic chemicals (without ammonia)'] = demand["Basic chemicals"] - ammonia_by_fuel
 
     demand['Basic chemicals (without ammonia)'].clip(lower=0, inplace=True)
 

scripts/build_industrial_energy_demand_per_node.py

@@ -9,6 +9,7 @@ if __name__ == '__main__':
             'build_industrial_energy_demand_per_node',
             simpl='',
             clusters=48,
+            planning_horizons=2030,
         )
         
     # import EU ratios df as csv

scripts/build_industry_sector_ratios.py

@@ -60,6 +60,7 @@ index = [
     "hydrogen",
     "heat",
     "naphtha",
+    "ammonia",
     "process emission",
     "process emission from feedstock",
 ]
@@ -432,6 +433,9 @@ def chemicals_industry():
 
     sector = "Ammonia"
     df[sector] = 0.0
+    if snakemake.config["sector"].get("ammonia", False):
+        df.loc["ammonia", sector] = config["MWh_NH3_per_tNH3"]
+    else:
         df.loc["hydrogen", sector] = config["MWh_H2_per_tNH3_electrolysis"]
         df.loc["elec", sector] = config["MWh_elec_per_tNH3_electrolysis"]
 
@@ -614,7 +618,7 @@ def nonmetalic_mineral_products():
     # (c) clinker production (kilns),
     # (d) Grinding, packaging.
     # (b)+(c) represent 94% of fec. So (a) is joined to (b) and (d) is joined to (c).
-    # Temperatures above 1400C are required for procesing limestone and sand into clinker.
+    # Temperatures above 1400C are required for processing limestone and sand into clinker.
     # Everything (except current electricity and heat consumption and existing biomass)
     # is transformed into methane for high T.
 
@@ -1106,7 +1110,7 @@ def non_ferrous_metals():
 
     # Aluminium secondary route
 
-    # All is coverted into secondary route fully electrified.
+    # All is converted into secondary route fully electrified.
 
     sector = "Aluminium - secondary production"
 

scripts/build_retro_cost.py

@@ -33,7 +33,7 @@ The basic equations:
 
         E_space = H_losses - H_gains
 
-    Heat losses constitute from the losses through heat trasmission (H_tr [W/m²K])
+    Heat losses constitute from the losses through heat transmission (H_tr [W/m²K])
     (this includes heat transfer through building elements and thermal bridges)
     and losses by ventilation (H_ve [W/m²K]):
 
@@ -71,7 +71,7 @@ import xarray as xr
 
 # thermal conductivity standard value
 k = 0.035
-# strenght of relative retrofitting depending on the component
+# strength of relative retrofitting depending on the component
 # determined by historical data of insulation thickness for retrofitting
 l_weight = pd.DataFrame({"weight": [1.95, 1.48, 1.]},
                         index=["Roof", "Wall", "Floor"])
@@ -89,8 +89,8 @@ tau_H_0 = 30
 # constant parameter alpha_H_0 [-] according to EN 13790 seasonal method
 alpha_H_0 = 0.8
 
-# paramter for solar heat load during heating season -------------------------
+# parameter for solar heat load during heating season -------------------------
-# tabular standard values table p.8 in documenation
+# tabular standard values table p.8 in documentation
 external_shading = 0.6     # vertical orientation: fraction of window area shaded [-]
 frame_area_fraction = 0.3  # fraction of frame area of window [-]
 non_perpendicular = 0.9    # reduction factor, considering radiation non perpendicular to the glazing[-]
@@ -279,7 +279,7 @@ def prepare_building_stock_data():
 def prepare_building_topology(u_values, same_building_topology=True):
     """
     reads in typical building topologies (e.g. average surface of building elements)
-    and typical losses trough thermal bridging and air ventilation
+    and typical losses through thermal bridging and air ventilation
     """
 
     data_tabula = pd.read_csv(snakemake.input.data_tabula,
@@ -585,7 +585,7 @@ def map_to_lstrength(l_strength, df):
 
 def calculate_heat_losses(u_values, data_tabula, l_strength, temperature_factor):
     """
-    calculates total annual heat losses Q_ht for different insulation thiknesses
+    calculates total annual heat losses Q_ht for different insulation thicknesses
     (l_strength), depening on current insulation state (u_values), standard
     building topologies and air ventilation from TABULA (data_tabula) and
     the accumulated difference between internal and external temperature
@@ -790,7 +790,7 @@ def sample_dE_costs_area(area, area_tot, costs, dE_space, countries,
 
     # drop not considered countries
     cost_dE = cost_dE.reindex(countries,level=0)
-    # get share of residential and sevice floor area
+    # get share of residential and service floor area
     sec_w = area_tot.value / area_tot.value.groupby(level=0).sum()
     # get the total cost-energy-savings weight by sector area
     tot = (cost_dE.mul(sec_w, axis=0).groupby(level="country_code").sum()
@@ -863,7 +863,7 @@ if __name__ == "__main__":
     data_tabula = prepare_building_topology(u_values)
     # costs for retrofitting -------------------------------------------------
     cost_retro, window_assumptions, cost_w, tax_w = prepare_cost_retro(country_iso_dic)
-    # temperature dependend parameters
+    # temperature dependent parameters
     d_heat, temperature_factor = prepare_temperature_data()
 
 

scripts/helper.py

@@ -25,7 +25,7 @@ def override_component_attrs(directory):
 
     Returns
     -------
-    Dictionary of overriden component attributes.
+    Dictionary of overridden component attributes.
     """
 
     attrs = Dict({k : v.copy() for k,v in component_attrs.items()})

scripts/make_summary.py

@@ -273,7 +273,7 @@ def calculate_supply(n, label, supply):
 
             for end in [col[3:] for col in c.df.columns if col[:3] == "bus"]:
 
-                items = c.df.index[c.df["bus" + end].map(bus_map, na_action=False)]
+                items = c.df.index[c.df["bus" + end].map(bus_map, na_action=None)]
 
                 if len(items) == 0:
                     continue
@@ -318,7 +318,7 @@ def calculate_supply_energy(n, label, supply_energy):
 
             for end in [col[3:] for col in c.df.columns if col[:3] == "bus"]:
 
-                items = c.df.index[c.df["bus" + str(end)].map(bus_map, na_action=False)]
+                items = c.df.index[c.df["bus" + str(end)].map(bus_map, na_action=None)]
 
                 if len(items) == 0:
                     continue

scripts/plot_network.py

@@ -23,6 +23,8 @@ def rename_techs_tyndp(tech):
         return "power-to-gas"
     elif tech == "H2":
         return "H2 storage"
+    elif tech in ["NH3", "Haber-Bosch", "ammonia cracker", "ammonia store"]:
+        return "ammonia"
     elif tech in ["OCGT", "CHP", "gas boiler", "H2 Fuel Cell"]:
         return "gas-to-power/heat"
     elif "solar" in tech:

scripts/plot_summary.py

@@ -52,6 +52,7 @@ def rename_techs(label):
         "ror": "hydroelectricity",
         "hydro": "hydroelectricity",
         "PHS": "hydroelectricity",
+        "NH3": "ammonia",
         "co2 Store": "DAC",
         "co2 stored": "CO2 sequestration",
         "AC": "transmission lines",
@@ -107,6 +108,7 @@ preferred_order = pd.Index([
     "natural gas",
     "helmeth",
     "methanation",
+    "ammonia",
     "hydrogen storage",
     "power-to-gas",
     "power-to-liquid",
@@ -255,7 +257,7 @@ def plot_balances():
         df = df / 1e6
 
         #remove trailing link ports
-        df.index = [i[:-1] if ((i != "co2") and (i[-1:] in ["0","1","2","3"])) else i for i in df.index]
+        df.index = [i[:-1] if ((i not in ["co2", "NH3"]) and (i[-1:] in ["0","1","2","3"])) else i for i in df.index]
 
         df = df.groupby(df.index.map(rename_techs)).sum()
 
@@ -399,7 +401,7 @@ def plot_carbon_budget_distribution(input_eurostat):
 
     ax1.plot(emissions, color='black', linewidth=3, label=None)
 
-    #plot commited and uder-discussion targets
+    #plot committed and uder-discussion targets
     #(notice that historical emissions include all countries in the
     # network, but targets refer to EU)
     ax1.plot([2020],[0.8*emissions[1990]],
@@ -425,7 +427,7 @@ def plot_carbon_budget_distribution(input_eurostat):
 
     ax1.plot([2050],[0.125*emissions[1990]],'ro',
                      marker='*', markersize=12, markerfacecolor='black',
-                     markeredgecolor='black', label='EU commited target')
+                     markeredgecolor='black', label='EU committed target')
 
     ax1.legend(fancybox=True, fontsize=18, loc=(0.01,0.01),
                        facecolor='white', frameon=True)

scripts/prepare_sector_network.py

@@ -93,6 +93,19 @@ def define_spatial(nodes, options):
 
     spatial.gas.df = pd.DataFrame(vars(spatial.gas), index=nodes)
 
+    # ammonia
+
+    if options.get('ammonia'):
+        spatial.ammonia = SimpleNamespace()
+        if options.get("ammonia") == "regional":
+            spatial.ammonia.nodes = nodes + " NH3"
+            spatial.ammonia.locations = nodes
+        else:
+            spatial.ammonia.nodes = ["EU NH3"]
+            spatial.ammonia.locations = ["EU"]
+
+        spatial.ammonia.df = pd.DataFrame(vars(spatial.ammonia), index=nodes)
+
     # oil
     spatial.oil = SimpleNamespace()
     spatial.oil.nodes = ["EU oil"]
@@ -664,6 +677,61 @@ def add_generation(n, costs):
         )
 
 
+def add_ammonia(n, costs):
+
+    logger.info("adding ammonia carrier with synthesis, cracking and storage")
+
+    nodes = pop_layout.index
+
+    cf_industry = snakemake.config["industry"]
+
+    n.add("Carrier", "NH3")
+
+    n.madd("Bus",
+        spatial.ammonia.nodes,
+        location=spatial.ammonia.locations,
+        carrier="NH3"
+    )
+
+    n.madd("Link",
+        nodes,
+        suffix=" Haber-Bosch",
+        bus0=nodes,
+        bus1=spatial.ammonia.nodes,
+        bus2=nodes + " H2",
+        p_nom_extendable=True,
+        carrier="Haber-Bosch",
+        efficiency=1 / (cf_industry["MWh_elec_per_tNH3_electrolysis"] / cf_industry["MWh_NH3_per_tNH3"]), # output: MW_NH3 per MW_elec
+        efficiency2=-cf_industry["MWh_H2_per_tNH3_electrolysis"] / cf_industry["MWh_elec_per_tNH3_electrolysis"], # input: MW_H2 per MW_elec
+        capital_cost=costs.at["Haber-Bosch synthesis", "fixed"],
+        lifetime=costs.at["Haber-Bosch synthesis", 'lifetime']
+    )
+
+    n.madd("Link",
+        nodes,
+        suffix=" ammonia cracker",
+        bus0=spatial.ammonia.nodes,
+        bus1=nodes + " H2",
+        p_nom_extendable=True,
+        carrier="ammonia cracker",
+        efficiency=1 / cf_industry["MWh_NH3_per_MWh_H2_cracker"],
+        capital_cost=costs.at["Ammonia cracker", "fixed"] / cf_industry["MWh_NH3_per_MWh_H2_cracker"], # given per MW_H2
+        lifetime=costs.at['Ammonia cracker', 'lifetime']
+    )
+
+    # Ammonia Storage
+    n.madd("Store",
+        spatial.ammonia.nodes,
+        suffix=" ammonia store",
+        bus=spatial.ammonia.nodes,
+        e_nom_extendable=True,
+        e_cyclic=True,
+        carrier="ammonia store",
+        capital_cost=costs.at["NH3 (l) storage tank incl. liquefaction", "fixed"],
+        lifetime=costs.at['NH3 (l) storage tank incl. liquefaction', 'lifetime']
+    )
+
+
 def add_wave(n, wave_cost_factor):
 
     # TODO: handle in Snakefile
@@ -1314,7 +1382,7 @@ def add_land_transport(n, costs):
 
 def build_heat_demand(n):
 
-    # copy forward the daily average heat demand into each hour, so it can be multipled by the intraday profile
+    # copy forward the daily average heat demand into each hour, so it can be multiplied by the intraday profile
     daily_space_heat_demand = xr.open_dataarray(snakemake.input.heat_demand_total).to_pandas().reindex(index=n.snapshots, method="ffill")
 
     intraday_profiles = pd.read_csv(snakemake.input.heat_profile, index_col=0)
@@ -1656,7 +1724,7 @@ def add_heat(n, costs):
 
             # minimum heat demand 'dE' after retrofitting in units of original heat demand (values between 0-1)
             dE = retro_data.loc[(ct, sec), ("dE")]
-            # get addtional energy savings 'dE_diff' between the different retrofitting strengths/generators at one node
+            # get additional energy savings 'dE_diff' between the different retrofitting strengths/generators at one node
             dE_diff = abs(dE.diff()).fillna(1-dE.iloc[0])
             # convert costs Euro/m^2 -> Euro/MWh
             capital_cost =  retro_data.loc[(ct, sec), ("cost")] * floor_area_node / \
@@ -2278,6 +2346,20 @@ def add_industry(n, costs):
         lifetime=costs.at['cement capture', 'lifetime']
     )
 
+    if options.get("ammonia"):
+
+        if options["ammonia"] == 'regional':
+            p_set = industrial_demand.loc[spatial.ammonia.locations, "ammonia"].rename(index=lambda x: x + " NH3") / 8760
+        else:
+            p_set = industrial_demand["ammonia"].sum() / 8760
+
+        n.madd("Load",
+            spatial.ammonia.nodes,
+            bus=spatial.ammonia.nodes,
+            carrier="NH3",
+            p_set=p_set
+        )
+
 
 def add_waste_heat(n):
     # TODO options?
@@ -2417,6 +2499,88 @@ def limit_individual_line_extension(n, maxext):
     hvdc = n.links.index[n.links.carrier == 'DC']
     n.links.loc[hvdc, 'p_nom_max'] = n.links.loc[hvdc, 'p_nom'] + maxext
 
+
+def apply_time_segmentation(n, segments, solver_name="cbc",
+                            overwrite_time_dependent=True):
+    """Aggregating time series to segments with different lengths
+
+    Input:
+        n: pypsa Network
+        segments: (int) number of segments in which the typical period should be
+                  subdivided
+        solver_name: (str) name of solver
+        overwrite_time_dependent: (bool) overwrite time dependent data of pypsa network
+        with typical time series created by tsam
+    """
+    try:
+        import tsam.timeseriesaggregation as tsam
+    except:
+        raise ModuleNotFoundError("Optional dependency 'tsam' not found."
+                                  "Install via 'pip install tsam'")
+
+    # get all time-dependent data
+    columns = pd.MultiIndex.from_tuples([],names=['component', 'key', 'asset'])
+    raw = pd.DataFrame(index=n.snapshots,columns=columns)
+    for c in n.iterate_components():
+        for attr, pnl in c.pnl.items():
+            # exclude e_min_pu which is used for SOC of EVs in the morning
+            if not pnl.empty and attr != 'e_min_pu':
+                df = pnl.copy()
+                df.columns = pd.MultiIndex.from_product([[c.name], [attr], df.columns])
+                raw = pd.concat([raw, df], axis=1)
+
+    # normalise all time-dependent data
+    annual_max = raw.max().replace(0,1)
+    raw = raw.div(annual_max, level=0)
+
+    # get representative segments
+    agg = tsam.TimeSeriesAggregation(raw, hoursPerPeriod=len(raw),
+                                     noTypicalPeriods=1, noSegments=int(segments),
+                                     segmentation=True, solver=solver_name)
+    segmented = agg.createTypicalPeriods()
+
+
+    weightings = segmented.index.get_level_values("Segment Duration")
+    offsets = np.insert(np.cumsum(weightings[:-1]), 0, 0)
+    timesteps = [raw.index[0] + pd.Timedelta(f"{offset}h") for offset in offsets]
+    snapshots =  pd.DatetimeIndex(timesteps)
+    sn_weightings = pd.Series(weightings, index=snapshots, name="weightings", dtype="float64")
+
+    n.set_snapshots(sn_weightings.index)
+    n.snapshot_weightings = n.snapshot_weightings.mul(sn_weightings, axis=0)
+
+    # overwrite time-dependent data with timeseries created by tsam
+    if overwrite_time_dependent:
+        values_t = segmented.mul(annual_max).set_index(snapshots)
+        for component, key in values_t.columns.droplevel(2).unique():
+            n.pnl(component)[key] = values_t[component, key]
+
+    return n
+
+def set_temporal_aggregation(n, opts, solver_name):
+    """Aggregate network temporally."""
+    for o in opts:
+        # temporal averaging
+        m = re.match(r"^\d+h$", o, re.IGNORECASE)
+        if m is not None:
+            n = average_every_nhours(n, m.group(0))
+            break
+        # representative snapshots
+        m = re.match(r"(^\d+)sn$", o, re.IGNORECASE)
+        if m is not None:
+            sn = int(m[1])
+            logger.info(f"use every {sn} snapshot as representative")
+            n.set_snapshots(n.snapshots[::sn])
+            n.snapshot_weightings *= sn
+            break
+        # segments with package tsam
+        m = re.match(r"^(\d+)seg$", o, re.IGNORECASE)
+        if m is not None:
+            segments = int(m[1])
+            logger.info(f"use temporal segmentation with {segments} segments")
+            n = apply_time_segmentation(n, segments, solver_name=solver_name)
+            break
+    return n
 #%%
 if __name__ == "__main__":
     if 'snakemake' not in globals():
@@ -2509,6 +2673,9 @@ if __name__ == "__main__":
     if options['dac']:
         add_dac(n, costs)
 
+    if options['ammonia']:
+        add_ammonia(n, costs)
+
     if "decentral" in opts:
         decentral(n)
 
@@ -2518,11 +2685,8 @@ if __name__ == "__main__":
     if options["co2_network"]:
         add_co2_network(n, costs)
 
-    for o in opts:
+    solver_name = snakemake.config["solving"]["solver"]["name"]
-        m = re.match(r'^\d+h$', o, re.IGNORECASE)
+    n = set_temporal_aggregation(n, opts, solver_name)
-        if m is not None:
-            n = average_every_nhours(n, m.group(0))
-            break
 
     limit_type = "config"
     limit = get(snakemake.config["co2_budget"], investment_year)

scripts/solve_network.py

@@ -227,7 +227,7 @@ def add_co2_sequestration_limit(n, sns):
     limit = n.config["sector"].get("co2_sequestration_potential", 200) * 1e6
     for o in opts:
         if not "seq" in o: continue
-        limit = float(o[o.find("seq")+3:])
+        limit = float(o[o.find("seq")+3:]) * 1e6
         break
 
     name = 'co2_sequestration_limit'

test/config.myopic.yaml

@@ -262,6 +262,9 @@ sector:
   biomass_transport: false  # biomass transport between nodes
   conventional_generation: # generator : carrier
     OCGT: gas
+  biomass_boiler: false
+  biomass_to_liquid: false
+  biosng: false
 
 
 industry:
@@ -316,6 +319,7 @@ industry:
   # Material Economics (2019): https://materialeconomics.com/latest-updates/industrial-transformation-2050
 
 costs:
+  year: 2030
   lifetime: 25 #default lifetime
   # From a Lion Hirth paper, also reflects average of Noothout et al 2016
   discountrate: 0.07

test/config.overnight.yaml

@@ -260,6 +260,9 @@ sector:
   biomass_transport: false  # biomass transport between nodes
   conventional_generation: # generator : carrier
     OCGT: gas
+  biomass_boiler: false
+  biomass_to_liquid: false
+  biosng: false
 
 
 industry:
@@ -314,6 +317,7 @@ industry:
   # Material Economics (2019): https://materialeconomics.com/latest-updates/industrial-transformation-2050
 
 costs:
+  year: 2030
   lifetime: 25 #default lifetime
   # From a Lion Hirth paper, also reflects average of Noothout et al 2016
   discountrate: 0.07