SSP emissions runs using CMIP6 calibrations
This notebook gives an example of running the SSP scenarios in FaIR using pre-calculated climate response calibrations from 66 CMIP6 models for a total of 8 × 66 = 528 scenarios in parallel. We will run in stochastic mode to demonstrate internal variability.
This demonstrates the speed in which FaIR can run scenarios and configurations in parallel. This might be a good template notebook to use for SSP emissions-driven runs.
Refer to the basic_emissions_run_example
for more detail on setting
up a run step by step; this notebook will be a little lighter on the
detail, but hopefully still enough to be useful.
0. import required modules
import numpy as np
import matplotlib.pyplot as pl
import pandas as pd
from fair import FAIR
from fair.io import read_properties
from fair.interface import fill, initialise
from fair.earth_params import seconds_per_year
1. initialise FaIR
We want to enable the methane lifetime routine that is a function of SLCFs and reactive gases, i.e. using the coefficients and feedbacks from Thornhill et al. 2021 and Skeie et al. 2020. We set this option in the initialiser this time.
f = FAIR(ch4_method='thornhill2021')
2. Define time horizon
create world running from 1750 to 2100, at 1-year intervals
f.define_time(1750, 2100, 1)
3. Define scenarios
We’ll take the SSPs; important that the names are consistent with those in the RCMIP database
# Define SSP scenarios
scenarios = ['ssp119', 'ssp126', 'ssp245', 'ssp370', 'ssp434', 'ssp460', 'ssp534-over', 'ssp585']
f.define_scenarios(scenarios)
4. Define configs
Our list of configs are going to be each CMIP6 climate model’s 4xCO2 response, which has been pre-calculated in the calibration notebooks.
We could also modify the response for different aerosol, ozone, methane lifetime tunings etc., but not every model has this data available.
df = pd.read_csv("../tests/test_data/4xCO2_cummins_ebm3.csv")
models = df['model'].unique()
configs = []
for imodel, model in enumerate(models):
for run in df.loc[df['model']==model, 'run']:
configs.append(f"{model}_{run}")
f.define_configs(configs)
5. Define species and properties
FaIR contains a few helper functions that populate the model with
sensible defaults. One is the read_properties
function that obtains
default species (the kitchen sink) and their properties for an
emissions-driven run
species, properties = read_properties()
#species = list(properties.keys())
species[:5]
properties['CO2 FFI']
f.define_species(species, properties)
~6. Modify run options~
Not required, because we set our run options in step 1.
7. Create input and output data
f.allocate()
8. Fill in the data
8a. get default species configs
Again we read in a default list of species configs that will apply to each config. If you want to change specific configs then you can still use this function to set defaults and tweak what you need. We will do this with the methane lifetime, which has a different value calibrated for the Thornhill 2021 lifetime option.
I’m also going to subtract the RCMIP 1750 emissions from CH4 and N2O. This is not in the default configs.
f.fill_species_configs()
fill(f.species_configs['unperturbed_lifetime'], 10.8537568, specie='CH4')
fill(f.species_configs['baseline_emissions'], 19.01978312, specie='CH4')
fill(f.species_configs['baseline_emissions'], 0.08602230754, specie='N2O')
8b. fill emissions
grab emissions (+solar and volcanic forcing) from RCMIP datasets using
the fill_from_rcmip
helper function. This function automatically
selects the emissions, concentration or forcing you want depending on
the properties
for each of the SSP scenarios defined.
I’m then going to make one change: replace the volcanic dataset with the AR6 volcanic dataset, as I want to compare the impact of monthly volcanic forcing in the monthly comparison.
We also need to initialise the first timestep of the run in terms of its per-species forcing, temperature, cumulative and airborne emissions. We set these all to zero. The concentration in the first timestep will be set to the baseline concentration, which are the IPCC AR6 1750 values.
df_volcanic = pd.read_csv('../tests/test_data/volcanic_ERF_monthly_175001-201912.csv', index_col='year')
df_volcanic[1750:].head()
f.fill_from_rcmip()
# overwrite volcanic
volcanic_forcing = np.zeros(351)
volcanic_forcing[:271] = df_volcanic[1749:].groupby(np.ceil(df_volcanic[1749:].index) // 1).mean().squeeze().values
fill(f.forcing, volcanic_forcing[:, None, None], specie="Volcanic") # sometimes need to expand the array
initialise(f.concentration, f.species_configs['baseline_concentration'])
initialise(f.forcing, 0)
initialise(f.temperature, 0)
initialise(f.cumulative_emissions, 0)
initialise(f.airborne_emissions, 0)
8c. fill climate configs
Take pre-calculated values from the Cummins et al. three layer model. We will use a reproducible random seed to define the stochastic behaviour.
df = pd.read_csv("../tests/test_data/4xCO2_cummins_ebm3.csv")
models = df['model'].unique()
seed = 1355763
for config in configs:
model, run = config.split('_')
condition = (df['model']==model) & (df['run']==run)
fill(f.climate_configs['ocean_heat_capacity'], df.loc[condition, 'C1':'C3'].values.squeeze(), config=config)
fill(f.climate_configs['ocean_heat_transfer'], df.loc[condition, 'kappa1':'kappa3'].values.squeeze(), config=config)
fill(f.climate_configs['deep_ocean_efficacy'], df.loc[condition, 'epsilon'].values[0], config=config)
fill(f.climate_configs['gamma_autocorrelation'], df.loc[condition, 'gamma'].values[0], config=config)
fill(f.climate_configs['sigma_eta'], df.loc[condition, 'sigma_eta'].values[0], config=config)
fill(f.climate_configs['sigma_xi'], df.loc[condition, 'sigma_xi'].values[0], config=config)
fill(f.climate_configs['stochastic_run'], True, config=config)
fill(f.climate_configs['use_seed'], True, config=config)
fill(f.climate_configs['seed'], seed, config=config)
seed = seed + 399
9. Run FaIR
look at it go.
You can turn off the progress bar with progress=False
.
f.run()
10. Make some nice plots
Presently this is accessed using the xarray
notation; perhaps we can
write a nice filter function like I did with fill
and
initialise
.
The output attributes of FAIR
of interest are - temperature (layer=0
is surface) - emissions (an output for GHGs driven with concentration) -
concentration (as above, vice versa) - forcing: the per-species
effective radiative forcing - forcing_sum: the total forcing -
airborne_emissions: total emissions of a GHG remaining in the atmosphere
- airborne_fraction: the fraction of GHG emissions remaining in the
atmosphere - alpha_lifetime: the scaling factor to unperturbed lifetime.
Mutiply the two values to get the atmospheric lifetime of a greenhouse
gas (see methane example below) - cumulative_emissions -
ocean_heat_content_change - toa_imbalance - stochastic_forcing: if
stochastic variability is activated, the non-deterministic part of the
forcing
pl.plot(f.timebounds, f.temperature.loc[dict(scenario='ssp119', layer=0)], label=f.configs);
pl.title('ssp119: temperature')
pl.xlabel('year')
pl.ylabel('Temperature anomaly (K)')
pl.plot(f.timebounds, f.species_configs['unperturbed_lifetime'].loc[dict(specie='CH4', gasbox=0)].data * f.alpha_lifetime.loc[dict(scenario='ssp119', specie='CH4')], label=f.configs);
pl.title('ssp119: methane lifetime')
pl.xlabel('year')
pl.ylabel('methane lifetime (yr)')
pl.plot(f.timebounds, f.forcing.loc[dict(scenario='ssp119', specie='CO2')], label=f.configs);
pl.title('ssp119: CO2 forcing')
pl.xlabel('year')
pl.ylabel('W/m2')
pl.plot(f.timebounds, f.forcing.loc[dict(scenario='ssp119', specie='CH4')], label=f.configs);
pl.title('ssp119: methane forcing')
pl.xlabel('year')
pl.ylabel('W/m2')
pl.plot(f.timebounds, f.concentration.loc[dict(scenario='ssp119', specie='CH4')], label=f.configs);
pl.title('ssp119: methane concentration')
pl.xlabel('year')
pl.ylabel('Temperature anomaly (K)')
pl.plot(f.timebounds, f.concentration.loc[dict(scenario='ssp119', specie='Equivalent effective stratospheric chlorine')], label=f.configs);
pl.title('ssp119: EESC')
pl.xlabel('year')
pl.ylabel('ppt')
pl.plot(f.timebounds, f.forcing.loc[dict(scenario='ssp119', specie='N2O')], label=f.configs);
pl.title('ssp119: N2O concentration')
pl.xlabel('year')
pl.ylabel('ppb')
pl.plot(f.timebounds, f.forcing.loc[dict(scenario='ssp119', specie='N2O')], label=f.configs);
pl.title('ssp119: N2O forcing')
pl.xlabel('year')
pl.ylabel('W/m2')
pl.plot(f.timebounds, f.concentration.loc[dict(scenario='ssp119', specie='CH3Cl')], label=f.configs);
pl.title('ssp119: Halon-1211 concentration')
pl.xlabel('year')
pl.ylabel('ppt')
pl.plot(f.timebounds, f.forcing.loc[dict(scenario='ssp119', specie='Aerosol-radiation interactions')], label=f.configs);
pl.title('ssp119: ERFari')
pl.xlabel('year')
pl.ylabel('W/m2')
pl.plot(f.timebounds, f.forcing.loc[dict(scenario='ssp119', specie='Aerosol-cloud interactions')], label=f.configs);
pl.title('ssp119: ERFaci')
pl.xlabel('year')
pl.ylabel('W/m2')
pl.plot(f.timebounds, f.forcing.loc[dict(scenario='ssp119', specie='Ozone')], label=f.configs);
pl.title('ssp119: Ozone forcing')
pl.xlabel('year')
pl.ylabel('W/m2')
pl.plot(f.timebounds, f.forcing.loc[dict(scenario='ssp119', specie='Contrails')], label=f.configs);
pl.title('ssp119: Contrails')
pl.xlabel('year')
pl.ylabel('W/m2')
pl.plot(f.timebounds, f.forcing.loc[dict(scenario='ssp119', specie='Light absorbing particles on snow and ice')], label=f.configs);
pl.title('ssp119: LAPSI')
pl.xlabel('year')
pl.ylabel('W/m2')
pl.plot(f.timebounds, f.forcing.loc[dict(scenario='ssp119', specie='Land use')], label=f.configs);
pl.title('ssp119: land use forcing')
pl.xlabel('year')
pl.ylabel('W/m2')
pl.plot(f.timebounds, f.forcing.loc[dict(scenario='ssp119', specie='Solar')], label=f.configs);
pl.title('ssp119: solar forcing')
pl.xlabel('year')
pl.ylabel('W/m2')
pl.plot(f.timebounds, f.forcing.loc[dict(scenario='ssp119', specie='Volcanic')], label=f.configs);
pl.title('ssp119: volcanic forcing')
pl.xlabel('year')
pl.ylabel('W/m2')
pl.plot(f.timebounds, f.forcing.loc[dict(scenario='ssp119', specie='Stratospheric water vapour')], label=f.configs);
pl.title('ssp119: Stratospheric water vapour forcing')
pl.xlabel('year')
pl.ylabel('W/m2')
pl.plot(f.timebounds, f.temperature.loc[dict(scenario='ssp126', layer=0)], label=f.configs);
pl.title('ssp126: temperature')
pl.xlabel('year')
pl.ylabel('Temperature anomaly (K)')
#pl.legend()
pl.plot(f.timebounds, f.temperature.loc[dict(scenario='ssp245', layer=0)], label=f.configs);
pl.title('ssp245: temperature')
pl.xlabel('year')
pl.ylabel('Temperature anomaly (K)')
#pl.legend()
pl.plot(f.timebounds, f.temperature.loc[dict(scenario='ssp370', layer=0)], label=f.configs);
pl.title('ssp370: temperature')
pl.xlabel('year')
pl.ylabel('Temperature anomaly (K)')
#pl.legend()
pl.plot(f.timebounds, f.temperature.loc[dict(scenario='ssp585', layer=0)], label=f.configs);
pl.title('ssp585: temperature')
pl.xlabel('year')
pl.ylabel('Temperature anomaly (K)')
#pl.legend()
pl.plot(f.timebounds, f.forcing.loc[dict(scenario='ssp126', specie='Ozone')], label=f.configs);
pl.title('ssp126: Ozone forcing')
pl.xlabel('year')
pl.ylabel('W/m2')
pl.plot(f.timebounds, f.forcing.loc[dict(scenario='ssp245', specie='Ozone')], label=f.configs);
pl.title('ssp245: Ozone forcing')
pl.xlabel('year')
pl.ylabel('W/m2')
pl.plot(f.timebounds, f.forcing.loc[dict(scenario='ssp370', specie='Ozone')], label=f.configs);
pl.title('ssp370: Ozone forcing')
pl.xlabel('year')
pl.ylabel('W/m2')
pl.plot(f.timebounds, f.forcing.loc[dict(scenario='ssp585', specie='Ozone')], label=f.configs);
pl.title('ssp585: Ozone forcing')
pl.xlabel('year')
pl.ylabel('W/m2')
pl.plot(f.timebounds, f.airborne_emissions.loc[dict(scenario='ssp585', specie='CO2')], label=f.configs);
pl.title('ssp585: Airborne emissions of CO2')
pl.xlabel('year')
pl.ylabel('GtCO2')
pl.plot(f.timebounds, f.airborne_fraction.loc[dict(scenario='ssp585', specie='CO2')], label=f.configs);
pl.title('ssp585: Airborne fraction of CO2')
pl.xlabel('year')
pl.ylabel('[1]')
pl.plot(f.timebounds, f.cumulative_emissions.loc[dict(scenario='ssp585', specie='CO2')], label=f.configs);
pl.title('ssp585: Cumulative emissions of CO2')
pl.xlabel('year')
pl.ylabel('GtCO2')
pl.plot(f.timebounds, f.ocean_heat_content_change.loc[dict(scenario='ssp585')], label=f.configs);
pl.title('ssp585: Ocean heat content change')
pl.xlabel('year')
pl.ylabel('J')
pl.plot(f.timebounds, f.toa_imbalance.loc[dict(scenario='ssp585')], label=f.configs);
pl.title('ssp585: Top of atmosphere energy imbalance')
pl.xlabel('year')
pl.ylabel('W m$^{-2}$')
pl.plot(f.timebounds, f.stochastic_forcing.loc[dict(scenario='ssp585')], label=f.configs);
pl.title('ssp585: Total forcing')
pl.xlabel('year')
pl.ylabel('W m$^{-2}$')
pl.plot(f.timebounds, f.forcing_sum.loc[dict(scenario='ssp585')], label=f.configs);
pl.title('ssp585: Deterministic forcing')
pl.xlabel('year')
pl.ylabel('W m$^{-2}$')
pl.plot(f.timebounds, f.stochastic_forcing.loc[dict(scenario='ssp585')]-f.forcing_sum.loc[dict(scenario='ssp585')], label=f.configs);
pl.title('ssp585: Stochastic forcing component')
pl.xlabel('year')
pl.ylabel('W m$^{-2}$')