Getting Started

This notebook demonstrates the OpenSCM Two Layer Model repository’s basic functionality.

We start with imports, their need will become clearer throughout the notebook.

import inspect

import numpy as np
from openscm_units import unit_registry
from scmdata import ScmRun

import openscm_twolayermodel
from openscm_twolayermodel import ImpulseResponseModel, TwoLayerModel
from openscm_twolayermodel.base import Model

/home/docs/checkouts/readthedocs.org/user_builds/openscm-two-layer-model/envs/stable/lib/python3.7/site-packages/openscm_twolayermodel/base.py:10: TqdmExperimentalWarning: Using `tqdm.autonotebook.tqdm` in notebook mode. Use `tqdm.tqdm` instead to force console mode (e.g. in jupyter console)
  import tqdm.autonotebook as tqdman

As with most Python packages, the version of openscm_twolayermodel being used can always be checked as shown below. This is very helpful for debugging.

# NBVAL_IGNORE_OUTPUT
openscm_twolayermodel.__version__

'0.2.3'

OpenSCM Two Layer Model has two key classes: ImpulseResponseModel and TwoLayerModel. These are implementations of the two major variants of the two-layer model found in the literature. We can see that they both have a common base class using the inspect package.

inspect.getmro(ImpulseResponseModel)

(openscm_twolayermodel.impulse_response_model.ImpulseResponseModel,
 openscm_twolayermodel.base.TwoLayerVariant,
 openscm_twolayermodel.base.Model,
 abc.ABC,
 object)

inspect.getmro(TwoLayerModel)

(openscm_twolayermodel.two_layer_model.TwoLayerModel,
 openscm_twolayermodel.base.TwoLayerVariant,
 openscm_twolayermodel.base.Model,
 abc.ABC,
 object)

These classes can both be used in the same way. We demonstrate the most basic usage here, more comprehensive usage is demonstrated in other notebooks.

The first thing we need is our effective radiative forcing driver. This should be an `ScmRun <https://scmdata.readthedocs.io/en/latest/data.html#the-scmrun-class>`__ instance.

run_length = 200

driver = ScmRun(
    data=np.arange(run_length) * 4 / 70,
    index=1850 + np.arange(run_length),
    columns={
        "unit": "W/m^2",
        "model": "idealised",
        "scenario": "1pctCO2",
        "region": "World",
        "variable": "Effective Radiative Forcing",
    },
)
driver

<scmdata.ScmRun (timeseries: 1, timepoints: 200)>
Time:
        Start: 1850-01-01T00:00:00
        End: 2049-01-01T00:00:00
Meta:
               model region scenario   unit                     variable
        0  idealised  World  1pctCO2  W/m^2  Effective Radiative Forcing

# NBVAL_IGNORE_OUTPUT
driver.lineplot()

<AxesSubplot:xlabel='time', ylabel='W/m^2'>

Then we can initialise instances of our models and run them.

# NBVAL_IGNORE_OUTPUT
two_layer = TwoLayerModel(lambda0=4 / 3 * unit_registry("W/m^2/delta_degC"))
res_two_layer = two_layer.run_scenarios(driver)

impulse_response = ImpulseResponseModel(d1=10 * unit_registry("yr"))
res_impulse_response = impulse_response.run_scenarios(driver)

res = res_two_layer.append(res_impulse_response)
res.head()

																time	1850-01-01	1851-01-01	1852-01-01	1853-01-01	1854-01-01	1855-01-01	1856-01-01	1857-01-01	1858-01-01	1859-01-01	...	2040-01-01	2041-01-01	2042-01-01	2043-01-01	2044-01-01	2045-01-01	2046-01-01	2047-01-01	2048-01-01	2049-01-01
a (watt / delta_degree_Celsius ** 2 / meter ** 2)	climate_model	d1 (a)	d2 (a)	dl (meter)	du (meter)	efficacy (dimensionless)	eta (watt / delta_degree_Celsius / meter ** 2)	lambda0 (watt / delta_degree_Celsius / meter ** 2)	model	q1 (delta_degree_Celsius * meter ** 2 / watt)	q2 (delta_degree_Celsius * meter ** 2 / watt)	region	run_idx	scenario	unit	variable
0.0	two_layer	NaN	NaN	1200.0	50.0	1.0	0.8	1.333333	idealised	NaN	NaN	World	0	1pctCO2	W/m^2	Effective Radiative Forcing	0.0	0.057143	0.114286	0.171429	0.228571	0.285714	0.342857	0.400000	0.457143	0.514286	...	10.857143	10.914286	10.971429	11.028571	11.085714	11.142857	11.200000	11.257143	11.314286	11.371429
															delta_degC	Surface Temperature|Upper	0.0	0.000000	0.008626	0.023100	0.041545	0.062689	0.085676	0.109924	0.135040	0.160761	...	5.710809	5.744627	5.778474	5.812348	5.846251	5.880182	5.914140	5.948127	5.982141	6.016183
																Surface Temperature|Lower	0.0	0.000000	0.000000	0.000043	0.000159	0.000368	0.000681	0.001109	0.001656	0.002328	...	1.937427	1.956414	1.975476	1.994612	2.013823	2.033107	2.052465	2.071897	2.091402	2.110980
															W/m^2	Heat Uptake	0.0	0.000000	0.057143	0.102784	0.140628	0.173178	0.202129	0.228623	0.253435	0.277089	...	3.230641	3.242731	3.254783	3.266797	3.278774	3.290713	3.302615	3.314480	3.326307	3.338098
NaN	two_timescale_impulse_response	10.0	400.0	NaN	NaN	1.0	NaN	NaN	idealised	0.3	0.4	World	0	1pctCO2	W/m^2	Effective Radiative Forcing	0.0	0.057143	0.114286	0.171429	0.228571	0.285714	0.342857	0.400000	0.457143	0.514286	...	10.857143	10.914286	10.971429	11.028571	11.085714	11.142857	11.200000	11.257143	11.314286	11.371429

5 rows × 200 columns

Now we can plot our outputs and compare (of course, we can make these two models the same if we’re clever about how we set the parameters, see the impulse response equivalence notebook).

# NBVAL_IGNORE_OUTPUT
res.filter(variable="Surface Temperature*").lineplot(
    hue="climate_model", style="variable"
)

<AxesSubplot:xlabel='time', ylabel='delta_degC'>

# NBVAL_IGNORE_OUTPUT
res.filter(variable="Heat*").lineplot(hue="climate_model", style="variable")

<AxesSubplot:xlabel='time', ylabel='W/m^2'>