OpenSCM Two Layer Model

OpenSCM two layer model contains implementations of the two layer radiative forcing driven models by Held et al., Geoffroy et al. and Bloch-Johnson et al.

OpenSCM Two Layer Model was designed to provide a clean, modularised, extensible interface for one of the most commonly used simple climate models. It was used in Phase 1 of the Reduced Complexity Model Intercomparison Project as a point of comparison for the other participating models.

The FaIR model implements a mathematically equivalent model (under certain assumptions) but does not provide as clear a conversion between the two-layer model and the two-timescale response as is provided here. We hope that this implementation could interface with other simple climate models like FaIR to allow simpler exploration of the combined behaviour of interacting climate components with minimal coupling headaches.

As implemented here, the “OpenSCM Two Layer Model” interface is target at researchers and as an education tool. Users from other fields are of course encouraged to use it if they find it helpful.

OpenSCM two layer model is free software under a BSD 3-Clause License, see LICENSE.

If you have any issues or would like to discuss a feature request, please raise them in the OpenSCM Two Layer Model issue tracker. If your issue is a feature request or a bug, please use the templates available, otherwise, simply open a normal issue.

Installation

OpenSCM two layer model has only two dependencies:

  • scmdata>=0.9

  • tqdm

OpenSCM two layer model can be installed with pip

pip install openscm-twolayermodel

If you also want to run the example notebooks install additional dependencies using

pip install "openscm-twolayermodel[notebooks]"

Coming soon OpenSCM two layer model can also be installed with conda

conda install -c conda-forge openscm-twolayermodel

We do not ship our tests with the packages. If you wish to run the tests, you must install from source (or download the tests separately and run them on your installation).

Installing from source

To install from source, simply clone the repository and then install it using pip e.g. pip install ".[dev]". Having done this, the tests can be run with pytest tests or using the Makefile (make test will run only the code tests, make checks will run the code tests and all other tests e.g. linting, notebooks, documentation).

For more details, see the development section of the docs.

Usage

Here we provide examples of OpenSCM two layer model’s behaviour and usage. The source code of these usage examples is available in the folder docs/source/usage of the GitHub repository.

Basic demos

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.

[1]:

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.

[2]:

# NBVAL_IGNORE_OUTPUT
openscm_twolayermodel.__version__

[2]:

'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.

[3]:

inspect.getmro(ImpulseResponseModel)

[3]:

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

[4]:

inspect.getmro(TwoLayerModel)

[4]:

(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.

[5]:

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

[5]:

<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

[6]:

# NBVAL_IGNORE_OUTPUT
driver.lineplot()

[6]:

<AxesSubplot:xlabel='time', ylabel='W/m^2'>
_images/usage_getting-started_11_1.png

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

[7]:

# 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()

[7]:
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).

[8]:

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

[8]:

<AxesSubplot:xlabel='time', ylabel='delta_degC'>
_images/usage_getting-started_15_1.png 
[9]:

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

[9]:

<AxesSubplot:xlabel='time', ylabel='W/m^2'>
_images/usage_getting-started_16_1.png

Running scenarios

Here we show how multiple scenarios can be run using the OpenSCM Two Layer Model package.

[1]:

# NBVAL_IGNORE_OUTPUT
import os.path

import numpy as np
import pandas as pd
from openscm_units import unit_registry as ur
import tqdm.autonotebook as tqdman
from scmdata import ScmRun, run_append

from openscm_twolayermodel import TwoLayerModel

import matplotlib.pyplot as plt

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

For this we use RCMIP effective radiative forcing data.

[2]:

DATA_PATH = os.path.join(
    "..",
    "..",
    "..",
    "tests",
    "test-data",
    "rcmip-radiative-forcing-annual-means-v4-0-0.csv",
)
DATA_PATH

[2]:

'../../../tests/test-data/rcmip-radiative-forcing-annual-means-v4-0-0.csv'

[3]:

# NBVAL_IGNORE_OUTPUT
scenarios = ScmRun(DATA_PATH, lowercase_cols=True).filter(
    scenario="historical", keep=False
)
scenarios

[3]:

<scmdata.ScmRun (timeseries: 480, timepoints: 751)>
Time:
        Start: 1750-01-01T00:00:00
        End: 2500-01-01T00:00:00
Meta:
                activity_id mip_era        model region          scenario   unit  \
        0    not_applicable   CMIP5          AIM  World             rcp60  W/m^2
        1    not_applicable   CMIP5          AIM  World             rcp60  W/m^2
        2    not_applicable   CMIP5          AIM  World             rcp60  W/m^2
        3    not_applicable   CMIP5          AIM  World             rcp60  W/m^2
        4    not_applicable   CMIP5          AIM  World             rcp60  W/m^2
        ..              ...     ...          ...    ...               ...    ...
        494  not_applicable   CMIP5  unspecified  World  historical-cmip5  W/m^2
        495  not_applicable   CMIP5  unspecified  World  historical-cmip5  W/m^2
        496  not_applicable   CMIP5  unspecified  World  historical-cmip5  W/m^2
        497  not_applicable   CMIP5  unspecified  World  historical-cmip5  W/m^2
        498  not_applicable   CMIP5  unspecified  World  historical-cmip5  W/m^2

                                                      variable
        0                                    Radiative Forcing
        1                      Radiative Forcing|Anthropogenic
        2             Radiative Forcing|Anthropogenic|Aerosols
        3    Radiative Forcing|Anthropogenic|Aerosols|Aeros...
        4    Radiative Forcing|Anthropogenic|Aerosols|Aeros...
        ..                                                 ...
        494  Radiative Forcing|Anthropogenic|Stratospheric ...
        495  Radiative Forcing|Anthropogenic|Tropospheric O...
        496                          Radiative Forcing|Natural
        497                    Radiative Forcing|Natural|Solar
        498                 Radiative Forcing|Natural|Volcanic

        [480 rows x 7 columns]

We can then run them, for a number of parameter settings, as shown.

[4]:

# NBVAL_IGNORE_OUTPUT
a_values = np.array([0, 0.01]) * ur("W/m^2/delta_degC^2")
a_values

[4]:
Magnitude
[0.0 0.01]
Unitswatt/(delta_degree_Celsius2 meter2)
[5]:

# NBVAL_IGNORE_OUTPUT
runner = TwoLayerModel()
output = []
for a in tqdman.tqdm(a_values, desc="Parameter settings"):
    runner.a = a
    output.append(runner.run_scenarios(scenarios))

output = run_append(output)
output

[5]:

<scmdata.ScmRun (timeseries: 80, timepoints: 751)>
Time:
        Start: 1750-01-01T00:00:00
        End: 2500-01-01T00:00:00
Meta:
            a (watt / delta_degree_Celsius ** 2 / meter ** 2)     activity_id  \
        0                                                0.00  not_applicable
        1                                                0.00  not_applicable
        2                                                0.00  not_applicable
        3                                                0.00  not_applicable
        4                                                0.00  not_applicable
        ..                                                ...             ...
        75                                               0.01  not_applicable
        76                                               0.01  not_applicable
        77                                               0.01  not_applicable
        78                                               0.01  not_applicable
        79                                               0.01  not_applicable

           climate_model  dl (meter)  du (meter)  efficacy (dimensionless)  \
        0      two_layer        1200          50                       1.0
        1      two_layer        1200          50                       1.0
        2      two_layer        1200          50                       1.0
        3      two_layer        1200          50                       1.0
        4      two_layer        1200          50                       1.0
        ..           ...         ...         ...                       ...
        75     two_layer        1200          50                       1.0
        76     two_layer        1200          50                       1.0
        77     two_layer        1200          50                       1.0
        78     two_layer        1200          50                       1.0
        79     two_layer        1200          50                       1.0

            eta (watt / delta_degree_Celsius / meter ** 2)  \
        0                                              0.8
        1                                              0.8
        2                                              0.8
        3                                              0.8
        4                                              0.8
        ..                                             ...
        75                                             0.8
        76                                             0.8
        77                                             0.8
        78                                             0.8
        79                                             0.8

            lambda0 (watt / delta_degree_Celsius / meter ** 2) mip_era          model  \
        0                                            1.246667    CMIP6        AIM/CGE
        1                                            1.246667    CMIP6        AIM/CGE
        2                                            1.246667    CMIP6        AIM/CGE
        3                                            1.246667    CMIP6        AIM/CGE
        4                                            1.246667    CMIP6        AIM/CGE
        ..                                                ...      ...            ...
        75                                           1.246667    CMIP6  REMIND-MAGPIE
        76                                           1.246667    CMIP6  REMIND-MAGPIE
        77                                           1.246667    CMIP6  REMIND-MAGPIE
        78                                           1.246667    CMIP6  REMIND-MAGPIE
        79                                           1.246667    CMIP6  REMIND-MAGPIE

           region  run_idx                   scenario        unit  \
        0   World        0                     ssp370       W/m^2
        1   World        0                     ssp370  delta_degC
        2   World        0                     ssp370  delta_degC
        3   World        0                     ssp370       W/m^2
        4   World        1  ssp370-lowNTCF-aerchemmip       W/m^2
        ..    ...      ...                        ...         ...
        75  World        8                ssp534-over       W/m^2
        76  World        9                     ssp585       W/m^2
        77  World        9                     ssp585  delta_degC
        78  World        9                     ssp585  delta_degC
        79  World        9                     ssp585       W/m^2

                               variable
        0   Effective Radiative Forcing
        1     Surface Temperature|Upper
        2     Surface Temperature|Lower
        3                   Heat Uptake
        4   Effective Radiative Forcing
        ..                          ...
        75                  Heat Uptake
        76  Effective Radiative Forcing
        77    Surface Temperature|Upper
        78    Surface Temperature|Lower
        79                  Heat Uptake

        [80 rows x 15 columns]

[6]:

# NBVAL_IGNORE_OUTPUT
pkwargs = dict(
    hue="scenario", style="a (watt / delta_degree_Celsius ** 2 / meter ** 2)",
)
fig = plt.figure(figsize=(12, 18))

ax = fig.add_subplot(211)
output.filter(variable="Surface Temperature|Upper").lineplot(**pkwargs, ax=ax)

ax = fig.add_subplot(212)
output.filter(variable="Heat Uptake").lineplot(**pkwargs, ax=ax)

[6]:

<AxesSubplot:xlabel='time', ylabel='W/m^2'>
_images/usage_running-scenarios_8_1.png

More detail

Impulse response equivalence

In this notebook we explore the equivalence between the two-layer model and a two-timescale impulse response approach.

Background

Following Geoffroy et al., 2013, Part 2, with notation altered to match our implementation, the two-layer model with efficacy (\(\epsilon\)) and state-dependent climate feedback can be written as

\begin{align} C \frac{dT}{dt} & = F - (\lambda_0 - a T) T - \epsilon \eta (T - T_D) \\ C_D \frac{dT_D}{dt} & = \eta (T - T_D) \end{align}

If the state-dependent feedback factor, \(a\), is non-zero, the two-layer model and impulse response approaches are not equivalent. However, if \(a=0\), they become the same.

Hereafter we assume \(a=0\), however this assumption should not be forgotten. In the case \(a=0\), the two-layer model can be written (adding an \(\epsilon\) for the deep-ocean equation too for simplicity later).

\begin{align} C \frac{dT}{dt} & = F - \lambda_0 T - \epsilon \eta (T - T_D) \\ \epsilon C_D \frac{dT_D}{dt} & = \epsilon \eta (T - T_D) \end{align}

In matrix notation we have

\begin{align} \frac{d\mathbf{X}}{dt} = \mathbf{A} \mathbf{X} + \mathbf{B} \end{align}

where \(\mathbf{X} = \begin{pmatrix}T \\T_D\end{pmatrix}\), \(\mathbf{A} = \begin{bmatrix} - \frac{\lambda_0 + \epsilon \eta}{C} & \frac{\epsilon \eta}{C} \\ \frac{\epsilon \eta}{\epsilon C_d} & -\frac{\epsilon \eta}{\epsilon C_d} \end{bmatrix}\) and \(\mathbf{B} = \begin{pmatrix} \frac{F}{C} \\ 0 \end{pmatrix}\).

As shown in Geoffroy et al., 2013, Part 1, \(\mathbf{A}\) can be diagonalised i.e. written in the form \(\mathbf{A} = \mathbf{\Phi} \mathbf{D} \mathbf{\Phi}^{-1}\), where \(\mathbf{D}\) is a diagonal matrix. Applying the solution given in Geoffroy et al., 2013, Part 1 to our impulse response notation, we have

\begin{align} \mathbf{D} = \begin{bmatrix} - \frac{1}{\tau_1} & 0 \\ 0 & - \frac{1}{\tau_2} \end{bmatrix} \end{align}

and

\begin{align} \mathbf{\Phi} = \begin{bmatrix} 1 & 1 \\ \phi_1 & \phi_2 \end{bmatrix} \end{align}

where

\begin{align} \tau_1 = \frac{C C_D}{2 \lambda_0 \eta} (b - \sqrt{\delta}) \end{align}
\begin{align} \tau_2 = \frac{C C_D}{2 \lambda_0 \eta} (b + \sqrt{\delta}) \end{align}
\begin{align} \phi_1 = \frac{C}{2 \epsilon \eta} (b^* - \sqrt{\delta}) \end{align}
\begin{align} \phi_2 = \frac{C}{2 \epsilon \eta} (b^* + \sqrt{\delta}) \end{align}
\begin{align} b = \frac{\lambda_0 + \epsilon \eta}{C} + \frac{\eta}{C_D} \end{align}
\begin{align} b^* = \frac{\lambda_0 + \epsilon \eta}{C} - \frac{\eta}{C_D} \end{align}
\begin{align} \delta = b^2 - 4 \frac{\lambda_0 \eta}{C C_D} \end{align}

Given this, we can re-write the system as

\begin{align} \frac{d\mathbf{X}}{dt} &= \mathbf{\Phi} \mathbf{D} \mathbf{\Phi}^{-1} \mathbf{X} + \mathbf{B} \\ \mathbf{\Phi}^{-1}\frac{d\mathbf{X}}{dt} &= \mathbf{D} \mathbf{\Phi}^{-1} \mathbf{X} + \mathbf{\Phi}^{-1} \mathbf{B} \\ \frac{d\mathbf{Y}}{dt} &= \mathbf{D} \mathbf{Y} + \mathbf{\Phi}^{-1} \mathbf{B} \\ \end{align}

Defining \(\mathbf{Y} = \begin{pmatrix} T_1 \\ T_2 \end{pmatrix}\), we have

\begin{align} \frac{d}{dt} \begin{pmatrix} T_1 \\ T_2 \end{pmatrix} = \begin{bmatrix} - \frac{1}{\tau_1} & 0 \\ 0 & - \frac{1}{\tau_2} \end{bmatrix} \begin{pmatrix} T_1 \\ T_2 \end{pmatrix} + \frac{1}{\phi_2 - \phi_1}\begin{bmatrix} \phi_2 & -1 \\ -\phi_1 & 1 \end{bmatrix} \begin{pmatrix} \frac{F}{C} \\ 0 \end{pmatrix} \end{align}

or,

\begin{align} \frac{dT_1}{dt} & = \frac{-T_1}{\tau_1} + \frac{\phi_2}{\phi_2 - \phi_1} \frac{F}{C} \\ \frac{dT_2}{dt} & = \frac{-T_2}{\tau_2} - \frac{\phi_1}{\phi_2 - \phi_1} \frac{F}{C} \end{align}

Re-writing, we have,

\begin{align} \frac{dT_1}{dt} & = \frac{1}{\tau_1} \left( \frac{\tau_1 \phi_2}{\phi_2 - \phi_1} \frac{F}{C} - T_1 \right) \\ \frac{dT_2}{dt} & = \frac{1}{\tau_2} \left( \frac{-\tau_2 \phi_1}{\phi_2 - \phi_1} \frac{F}{C} - T_2 \right) \end{align}

We can compare this to the notation of Millar et al., 2017 and see that

\begin{align} d_1 = \tau_1 \\ d_2 = \tau_2 \\ q_1 = \frac{\tau_1 \phi_2}{C(\phi_2 - \phi_1)} \\ q_2 = - \frac{\tau_2 \phi_1}{C(\phi_2 - \phi_1)} \\ \end{align}

Hence we have redemonstrated the equivalence of the two-layer model and a two-timescale impulse response model. Given the parameters of the two-layer model, we can now trivially derive the equivalent parameters of the two-timescale model. Doing the reverse is possible, but requires some more work in order to make a useable route drop out.

The first step is to follow Geoffroy et al., 2013, Part 1, and define two extra constants

\begin{align} a_1 = \frac{\phi_2 \tau_1}{C(\phi_2 - \phi_1)} \lambda_0 \\ a_2 = - \frac{\phi_1 \tau_2}{C(\phi_2 - \phi_1)} \lambda_0 \\ \end{align}

These constants have the useful property that \(a_1 + a_2 = 1\) (proof in Appendix A).

From above, we also see that

\begin{align} a_1 = \lambda_0 q_1 \\ a_2 = \lambda_0 q_2 \\ \end{align}

Hence

\begin{align} a_1 + a_2 = \lambda_0 q_1 + \lambda_0 q_2 = 1 \\ \lambda_0 = \frac{1}{q_1 + q_2} \end{align}

Next we calculate \(C\) via

\begin{align} \frac{q_1}{d_1} + \frac{q_2}{d_2} &= \frac{\phi_2}{C(\phi_2 - \phi_1)} - \frac{\phi_1}{C(\phi_2 - \phi_1)} = \frac{1}{C} \\ C &= \frac{d_1 d_2}{q_1 d_2 + q_2 d_1} \end{align}

We then use further relationships from Table 1 of Geoffroy et al., 2013, Part 1 (proof is left to the reader) to calculate the rest of the constants.

Firstly,

\begin{align} \tau_1 a_1 + \tau_2 a_2 = \frac{C + \epsilon C_D}{\lambda_0} \\ \epsilon C_D = \lambda_0 (\tau_1 a_1 + \tau_2 a_2) - C \end{align}

and then finally,

\begin{align} \tau_1 a_1 + \tau_2 a_2 = \frac{C + \epsilon C_D}{\lambda_0} \\ \epsilon \eta = \frac{\epsilon C_D}{\tau_1 a_2 + \tau_2 a_1} \end{align}

The final thing to notice here is that \(C_D\), \(\epsilon\) and \(\eta\) are not uniquely-defined. This makes sense, as shown by Geoffroy et al., 2013, Part 2, the introduction of the efficacy factor does not alter the behaviour of the system (it is still the same mathematical system) and so it is impossible for simply the two-timescale temperature response to uniquely define all three of these quantities. It can only define the products \(\epsilon C_D\) and \(\epsilon \eta\). Hence when translating from the two-timescale model to the two-layer model with efficacy, an explicit choice for the efficacy must be made. This does not alter the temperature response but it does alter the implied ocean heat uptake of the two-timescale model.

Long story short, when deriving two-layer model parameters from a two-timescale model, one must specify the efficacy.

Given that \(\mathbf{Y} = \mathbf{\Phi}^{-1} \mathbf{X}\) i.e. \(\mathbf{X} = \mathbf{\Phi} \mathbf{Y}\), we can also relate the impulse response boxes to the two layers.

\begin{aligned} \begin{pmatrix} T \\ T_D \end{pmatrix} = \begin{bmatrix} 1 & 1 \\ \phi_1 & \phi_2 \end{bmatrix} \begin{pmatrix} T_1 \\ T_2 \end{pmatrix} \\ \begin{pmatrix} T \\ T_D \end{pmatrix} = \begin{pmatrix} T_1 + T_2 \\ \phi_1 T_1 + \phi_2 T_2 \end{pmatrix} \end{aligned}

Finally, the equivalent of the two-timescale and two-layer models allows us to also calculate the heat uptake of a two-timescale impulse response model. It is given by

\begin{aligned} \text{Heat uptake} &= C \frac{dT}{dt} + C_D \frac{dT_D}{dt} \\ &= F - \lambda_0 T + (1 - \epsilon) \eta (T - T_D) \\ &= F - \lambda_0 (T_1 + T_2) + (1 - \epsilon) \eta ((1 - \phi_1) T_1 + (1 - \phi_2) T_2) \\ &= F - \lambda_0 (T_1 + T_2) - \eta (\epsilon - 1)((1 - \phi_1) T_1 + (1 - \phi_2) T_2) \end{aligned}
Running the code

Here we actually run the two implementations to explore their similarity.

[1]:

import datetime as dt

import numpy as np
import pandas as pd
from openscm_units import unit_registry as ur
from scmdata.run import ScmRun, run_append

from openscm_twolayermodel import ImpulseResponseModel, TwoLayerModel

import matplotlib.pyplot as plt

/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

First we define a scenario to run.

[2]:

time = np.arange(1750, 2501)
forcing = 0.05 * np.sin(time / 15 * 2 * np.pi) + 3.0 * time / time.max()

inp = ScmRun(
    data=forcing,
    index=time,
    columns={
        "scenario": "test_scenario",
        "model": "unspecified",
        "climate_model": "junk input",
        "variable": "Effective Radiative Forcing",
        "unit": "W/m^2",
        "region": "World",
    },
)
inp

[2]:

<scmdata.ScmRun (timeseries: 1, timepoints: 751)>
Time:
        Start: 1750-01-01T00:00:00
        End: 2500-01-01T00:00:00
Meta:
          climate_model        model region       scenario   unit  \
        0    junk input  unspecified  World  test_scenario  W/m^2

                              variable
        0  Effective Radiative Forcing

[3]:

# NBVAL_IGNORE_OUTPUT
inp.lineplot()

[3]:

<AxesSubplot:xlabel='time', ylabel='W/m^2'>
_images/usage_impulse-response-equivalence_6_1.png

Next we run the two-layer model. In order for it to be convertible to a two-timescale model, we must turn state-dependence off (a=0).

[4]:

two_layer_config = {
    "du": 55 * ur("m"),
    "efficacy": 1.2 * ur("dimensionless"),
    #     "efficacy": 1.0 * ur("dimensionless"),
    "a": 0 * ur("W/m^2/delta_degC^2"),
}

[5]:

# NBVAL_IGNORE_OUTPUT
twolayer = TwoLayerModel(**two_layer_config)
res_twolayer = twolayer.run_scenarios(inp)
res_twolayer

[5]:

<scmdata.ScmRun (timeseries: 4, timepoints: 751)>
Time:
        Start: 1750-01-01T00:00:00
        End: 2500-01-01T00:00:00
Meta:
           a (watt / delta_degree_Celsius ** 2 / meter ** 2) climate_model  \
        0                                                0.0     two_layer
        1                                                0.0     two_layer
        2                                                0.0     two_layer
        3                                                0.0     two_layer

           dl (meter)  du (meter)  efficacy (dimensionless)  \
        0        1200          55                       1.2
        1        1200          55                       1.2
        2        1200          55                       1.2
        3        1200          55                       1.2

           eta (watt / delta_degree_Celsius / meter ** 2)  \
        0                                             0.8
        1                                             0.8
        2                                             0.8
        3                                             0.8

           lambda0 (watt / delta_degree_Celsius / meter ** 2)        model region  \
        0                                           1.246667   unspecified  World
        1                                           1.246667   unspecified  World
        2                                           1.246667   unspecified  World
        3                                           1.246667   unspecified  World

           run_idx       scenario        unit                     variable
        0        0  test_scenario       W/m^2  Effective Radiative Forcing
        1        0  test_scenario  delta_degC    Surface Temperature|Upper
        2        0  test_scenario  delta_degC    Surface Temperature|Lower
        3        0  test_scenario       W/m^2                  Heat Uptake

[6]:

# NBVAL_IGNORE_OUTPUT
fig = plt.figure(figsize=(16, 9))

ax = fig.add_subplot(121)
res_twolayer.filter(variable="*Temperature*").lineplot(hue="variable", ax=ax)

ax = fig.add_subplot(122)
res_twolayer.filter(variable="Heat Uptake").lineplot(hue="variable", ax=ax)

[6]:

<AxesSubplot:xlabel='time', ylabel='W/m^2'>
_images/usage_impulse-response-equivalence_10_1.png

Next we get the parameters with which we get the equivalent impulse response model.

[7]:

two_timescale_paras = twolayer.get_impulse_response_parameters()
two_timescale_paras

[7]:

{'d1': 103454323.57029569 <Unit('joule / watt')>,
 'd2': 11181891933.114195 <Unit('joule / watt')>,
 'q1': 0.4465999986742509 <Unit('delta_degree_Celsius * meter ** 2 / watt')>,
 'q2': 0.3555390387589074 <Unit('delta_degree_Celsius * meter ** 2 / watt')>,
 'efficacy': 1.2 <Unit('dimensionless')>}

[8]:

# NBVAL_IGNORE_OUTPUT
impulse_response = ImpulseResponseModel(**two_timescale_paras)
res_impulse_response = impulse_response.run_scenarios(inp)
res_impulse_response

[8]:

<scmdata.ScmRun (timeseries: 5, timepoints: 751)>
Time:
        Start: 1750-01-01T00:00:00
        End: 2500-01-01T00:00:00
Meta:
                            climate_model  d1 (joule / watt)  d2 (joule / watt)  \
        0  two_timescale_impulse_response       1.034543e+08       1.118189e+10
        1  two_timescale_impulse_response       1.034543e+08       1.118189e+10
        2  two_timescale_impulse_response       1.034543e+08       1.118189e+10
        3  two_timescale_impulse_response       1.034543e+08       1.118189e+10
        4  two_timescale_impulse_response       1.034543e+08       1.118189e+10

           efficacy (dimensionless)        model  \
        0                       1.2  unspecified
        1                       1.2  unspecified
        2                       1.2  unspecified
        3                       1.2  unspecified
        4                       1.2  unspecified

           q1 (delta_degree_Celsius * meter ** 2 / watt)  \
        0                                         0.4466
        1                                         0.4466
        2                                         0.4466
        3                                         0.4466
        4                                         0.4466

           q2 (delta_degree_Celsius * meter ** 2 / watt) region  run_idx  \
        0                                       0.355539  World        0
        1                                       0.355539  World        0
        2                                       0.355539  World        0
        3                                       0.355539  World        0
        4                                       0.355539  World        0

                scenario        unit                     variable
        0  test_scenario       W/m^2  Effective Radiative Forcing
        1  test_scenario  delta_degC    Surface Temperature|Box 1
        2  test_scenario  delta_degC    Surface Temperature|Box 2
        3  test_scenario  delta_degC          Surface Temperature
        4  test_scenario       W/m^2                  Heat Uptake

[9]:

# NBVAL_IGNORE_OUTPUT
fig = plt.figure(figsize=(16, 9))

ax = fig.add_subplot(121)
res_impulse_response.filter(variable="*Temperature*").lineplot(
    hue="variable", ax=ax
)

ax = fig.add_subplot(122)
res_impulse_response.filter(variable="Heat Uptake").lineplot(
    hue="variable", ax=ax
)

[9]:

<AxesSubplot:xlabel='time', ylabel='W/m^2'>
_images/usage_impulse-response-equivalence_14_1.png

We can compare the two responses as well.

[10]:

# NBVAL_IGNORE_OUTPUT
combined = run_append([res_impulse_response, res_twolayer])
combined

[10]:

<scmdata.ScmRun (timeseries: 9, timepoints: 751)>
Time:
        Start: 1750-01-01T00:00:00
        End: 2500-01-01T00:00:00
Meta:
           a (watt / delta_degree_Celsius ** 2 / meter ** 2)  \
        0                                                NaN
        1                                                NaN
        2                                                NaN
        3                                                NaN
        4                                                NaN
        5                                                0.0
        6                                                0.0
        7                                                0.0
        8                                                0.0

                            climate_model  d1 (joule / watt)  d2 (joule / watt)  \
        0  two_timescale_impulse_response       1.034543e+08       1.118189e+10
        1  two_timescale_impulse_response       1.034543e+08       1.118189e+10
        2  two_timescale_impulse_response       1.034543e+08       1.118189e+10
        3  two_timescale_impulse_response       1.034543e+08       1.118189e+10
        4  two_timescale_impulse_response       1.034543e+08       1.118189e+10
        5                       two_layer                NaN                NaN
        6                       two_layer                NaN                NaN
        7                       two_layer                NaN                NaN
        8                       two_layer                NaN                NaN

           dl (meter)  du (meter)  efficacy (dimensionless)  \
        0         NaN         NaN                       1.2
        1         NaN         NaN                       1.2
        2         NaN         NaN                       1.2
        3         NaN         NaN                       1.2
        4         NaN         NaN                       1.2
        5      1200.0        55.0                       1.2
        6      1200.0        55.0                       1.2
        7      1200.0        55.0                       1.2
        8      1200.0        55.0                       1.2

           eta (watt / delta_degree_Celsius / meter ** 2)  \
        0                                             NaN
        1                                             NaN
        2                                             NaN
        3                                             NaN
        4                                             NaN
        5                                             0.8
        6                                             0.8
        7                                             0.8
        8                                             0.8

           lambda0 (watt / delta_degree_Celsius / meter ** 2)        model  \
        0                                                NaN   unspecified
        1                                                NaN   unspecified
        2                                                NaN   unspecified
        3                                                NaN   unspecified
        4                                                NaN   unspecified
        5                                           1.246667   unspecified
        6                                           1.246667   unspecified
        7                                           1.246667   unspecified
        8                                           1.246667   unspecified

           q1 (delta_degree_Celsius * meter ** 2 / watt)  \
        0                                         0.4466
        1                                         0.4466
        2                                         0.4466
        3                                         0.4466
        4                                         0.4466
        5                                            NaN
        6                                            NaN
        7                                            NaN
        8                                            NaN

           q2 (delta_degree_Celsius * meter ** 2 / watt) region  run_idx  \
        0                                       0.355539  World        0
        1                                       0.355539  World        0
        2                                       0.355539  World        0
        3                                       0.355539  World        0
        4                                       0.355539  World        0
        5                                            NaN  World        0
        6                                            NaN  World        0
        7                                            NaN  World        0
        8                                            NaN  World        0

                scenario        unit                     variable
        0  test_scenario       W/m^2  Effective Radiative Forcing
        1  test_scenario  delta_degC    Surface Temperature|Box 1
        2  test_scenario  delta_degC    Surface Temperature|Box 2
        3  test_scenario  delta_degC          Surface Temperature
        4  test_scenario       W/m^2                  Heat Uptake
        5  test_scenario       W/m^2  Effective Radiative Forcing
        6  test_scenario  delta_degC    Surface Temperature|Upper
        7  test_scenario  delta_degC    Surface Temperature|Lower
        8  test_scenario       W/m^2                  Heat Uptake

To within numerical errors they are equal.

[11]:

# NBVAL_IGNORE_OUTPUT
fig = plt.figure(figsize=(16, 9))

ax = fig.add_subplot(121)
combined.filter(variable="*Temperature*").lineplot(
    hue="variable", style="climate_model", alpha=0.7, linewidth=2, ax=ax
)
ax.legend(loc="upper left")

ax = fig.add_subplot(122)
combined.filter(variable="Heat Uptake").lineplot(
    hue="climate_model", style="climate_model", alpha=0.7, linewidth=2, ax=ax
)

[11]:

<AxesSubplot:xlabel='time', ylabel='W/m^2'>
_images/usage_impulse-response-equivalence_18_1.png
Appendix A

We begin with the definitions of the \(a\) constants,

\begin{align} a_1 = \frac{\phi_2 \tau_1}{C(\phi_2 - \phi_1)} \lambda_0 \\ a_2 = - \frac{\phi_1 \tau_2}{C(\phi_2 - \phi_1)} \lambda_0 \\ \end{align}

We then have

\begin{align} a_1 + a_2 &= \frac{\phi_2 \tau_1}{C(\phi_2 - \phi_1)} \lambda_0 - \frac{\phi_1 \tau_2}{C(\phi_2 - \phi_1)} \lambda_0 \\ &= \frac{\lambda_0}{C(\phi_2 - \phi_1)} (\phi_2 \tau_1 - \phi_1 \tau_2) \end{align}

Recalling the definition of the \(\phi\) parameters,

\begin{align} \phi_1 = \frac{C}{2 \epsilon \eta} (b^* - \sqrt{\delta}) \end{align}
\begin{align} \phi_2 = \frac{C}{2 \epsilon \eta} (b^* + \sqrt{\delta}) \end{align}

We have,

\begin{align} \phi_2 - \phi_1 = \frac{C \sqrt{\delta}}{\epsilon \eta} \end{align}

Recalling the definition of the \(\tau\) parameters,

\begin{align} \tau_1 = \frac{C C_D}{2 \lambda_0 \eta} (b - \sqrt{\delta}) \end{align}
\begin{align} \tau_2 = \frac{C C_D}{2 \lambda_0 \eta} (b + \sqrt{\delta}) \end{align}

We have,

\begin{align} \phi_2 \tau_1 - \phi_1 \tau_2 &= \frac{C}{2 \epsilon \eta} (b^* + \sqrt{\delta}) \times \frac{C C_D}{2 \lambda_0 \eta} (b - \sqrt{\delta}) - \frac{C}{2 \epsilon \eta} (b^* - \sqrt{\delta}) \times \frac{C C_D}{2 \lambda_0 \eta} (b + \sqrt{\delta}) \\ &= \frac{C^2 C_d}{4 \epsilon \eta^2 \lambda_0} \left[ (b^* + \sqrt{\delta}) (b - \sqrt{\delta}) - (b^* - \sqrt{\delta}) (b + \sqrt{\delta}) \right] \\ &= \frac{C^2 C_d}{4 \epsilon \eta^2 \lambda_0} \left[ b^*b + b\sqrt{\delta} -b^*\sqrt{\delta} - \delta - bb^* + b\sqrt{\delta} - b^*\sqrt{\delta} + \delta \right] \\ &= \frac{C^2 C_d}{2 \epsilon \eta^2 \lambda_0} \left[ b\sqrt{\delta} -b^*\sqrt{\delta} \right] \\ &= \frac{C^2 C_d \sqrt{\delta}}{2 \epsilon \eta^2 \lambda_0} \left[ b -b^*\right] \end{align}

Recalling the definition of the \(b\) parameters,

\begin{align} b = \frac{\lambda_0 + \epsilon \eta}{C} + \frac{\eta}{C_D} \end{align}
\begin{align} b^* = \frac{\lambda_0 + \epsilon \eta}{C} - \frac{\eta}{C_D} \end{align}

We then have

\begin{align} \phi_2 \tau_1 - \phi_1 \tau_2 &= \frac{C^2 C_d \sqrt{\delta}}{2 \epsilon \eta^2 \lambda_0} \left[ \frac{2\eta}{C_D} \right] \\ &= \frac{C^2 \sqrt{\delta}}{\epsilon \eta \lambda_0} \end{align}

Putting it all back together,

\begin{align} a_1 + a_2 &= \frac{\lambda_0}{C(\phi_2 - \phi_1)} (\phi_2 \tau_1 - \phi_1 \tau_2) \\ &= \frac{\lambda_0}{C} \frac{\epsilon \eta}{C \sqrt{\delta}} \frac{C^2 \sqrt{\delta}}{\epsilon \eta \lambda_0} \\ &= 1 \end{align}

One layer model

Here we show how to run our two-layer model as a single-layer model. There are two different ways to do this, which we present below.

Imports and loading data
[1]:

# NBVAL_IGNORE_OUTPUT
import os.path

import numpy as np
import pandas as pd
from openscm_units import unit_registry as ur
import tqdm.autonotebook as tqdman
from scmdata import ScmRun, run_append

from openscm_twolayermodel import TwoLayerModel

import matplotlib.pyplot as plt

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

For this we use an idealised scenario which is a reasonable representation of the forcing which occurs in response to an abrupt doubling in atmospheric CO\(_2\) concentrations (often referred to as an abrupt-2xCO2 experiment).

[2]:

run_length = 2000

data = np.zeros(run_length)
data[10 : ] = 4.0

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

[2]:

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

[3]:

# NBVAL_IGNORE_OUTPUT
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111)
driver.filter(variable="Effective Radiative Forcing").lineplot()

[3]:

<AxesSubplot:xlabel='time', ylabel='W/m^2'>
_images/usage_one-layer-model_5_1.png
No second layer

The first, and arguably the simplest way, to make a single layer model is to simply remove the connection between the top and second layers. Recalling the equations which define the two-layer model below,

\begin{align} C \frac{dT}{dt} & = F - (\lambda_0 - a T) T - \epsilon \eta (T - T_D) \\ C_D \frac{dT_D}{dt} & = \eta (T - T_D) \end{align}

We see that we can effectively remove the second layer by setting \(\eta = 0\).

[4]:

TwoLayerModel().eta

[4]:
0.8 watt/(delta_degree_Celsius meter2)
[5]:

# NBVAL_IGNORE_OUTPUT
eta_values = np.array([0, 0.8]) * ur("W/m^2/K")
eta_values

[5]:
Magnitude
[0.0 0.8]
Unitswatt/(kelvin meter2)
[6]:

# NBVAL_IGNORE_OUTPUT
du_values = np.array([50, 250, 500]) * ur("m")
du_values

[6]:
Magnitude
[ 50 250 500]
Unitsmeter
[7]:

# NBVAL_IGNORE_OUTPUT
runner = TwoLayerModel()
output = []
equivalent_parameters = []
for eta in tqdman.tqdm(eta_values, desc="eta values", leave=False):
    runner.eta = eta
    for du in tqdman.tqdm(du_values, desc="du values", leave=False):
        runner.du = du
        output.append(runner.run_scenarios(driver))

output = run_append(output)
output.head()

[7]:
time 1850-01-01 00:00:00 1851-01-01 00:00:00 1852-01-01 00:00:00 1853-01-01 00:00:00 1854-01-01 00:00:00 1855-01-01 00:00:00 1856-01-01 00:00:00 1857-01-01 00:00:00 1858-01-01 00:00:00 1859-01-01 00:00:00 ... 3840-01-01 00:00:00 3841-01-01 00:00:00 3842-01-01 00:00:00 3843-01-01 00:00:00 3844-01-01 00:00:00 3845-01-01 00:00:00 3846-01-01 00:00:00 3847-01-01 00:00:00 3848-01-01 00:00:00 3849-01-01 00:00:00
a (watt / delta_degree_Celsius ** 2 / meter ** 2) climate_model dl (meter) du (meter) efficacy (dimensionless) eta (watt / kelvin / meter ** 2) lambda0 (watt / delta_degree_Celsius / meter ** 2) model region run_idx scenario unit variable
0.0 two_layer 1200 50 1.0 0.0 1.246667 idealised World 0 1pctCO2 W/m^2 Effective Radiative Forcing 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 4.000000 4.000000 4.000000 4.000000 4.000000 4.000000 4.000000 4.000000 4.000000 4.000000
delta_degC Surface Temperature|Upper 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 3.208556 3.208556 3.208556 3.208556 3.208556 3.208556 3.208556 3.208556 3.208556 3.208556
Surface Temperature|Lower 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
W/m^2 Heat Uptake 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
250 1.0 0.0 1.246667 idealised World 0 1pctCO2 W/m^2 Effective Radiative Forcing 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 4.000000 4.000000 4.000000 4.000000 4.000000 4.000000 4.000000 4.000000 4.000000 4.000000

5 rows × 2000 columns

As we can see in the plots below, the runs with \(\eta = 0\) only have a single timescale in their response. In contrast, the runs with \(\eta \neq 0\) have two clear, distinct timescales. Notably, because equilibrium warming is independent of ocean heat uptake, the equilibrium warming is the same in all cases.

As expected, we see that the depth of the mixed-layer affects the response time of the mixed-layer (the only response time in the case of \(\eta = 0\)) whilst having a much smaller effect on the response time of the deep ocean.

[8]:

# NBVAL_IGNORE_OUTPUT
scenario_to_plot = "1pctCO2"
xlim = [1850, 3500]
pkwargs = dict(
    hue="du (meter)",
    style="eta (watt / kelvin / meter ** 2)",
    time_axis="year"
)
fig = plt.figure(figsize=(9, 9))

ax = fig.add_subplot(211)
output.filter(scenario=scenario_to_plot, variable="Surface Temperature|Upper").lineplot(**pkwargs, ax=ax)
ax.set_title("Surface Temperature|Upper")

ax = fig.add_subplot(212, sharex=ax)
output.filter(scenario=scenario_to_plot, variable="Heat Uptake").lineplot(**pkwargs, ax=ax)
ax.set_title("Heat Uptake")
ax.set_xlim(xlim)

plt.tight_layout()
_images/usage_one-layer-model_12_0.png
Infinite reservoir second layer

If we make the deep ocean component of the two-layer model infinitely deep, then we also have a single layer model. The concept is described by Equation 4 of Geoffroy et al. 2013, Part 1.

\begin{align} C \frac{dT}{dt} & = F - (\lambda_0 - a T) T - \epsilon \eta (T - T_D) \\ C_D \frac{dT_D}{dt} & = \eta (T - T_D) \end{align}

In short, if \(C_D \rightarrow \infty\), then \(T_D = 0\) and the equation governing the mixed layer response becomes

\begin{align} C \frac{dT}{dt} & = F - (\lambda_0 - a T) T - \epsilon \eta T \end{align}

In effect, we alter the climate feedback factor from \(\lambda_0 - a T\) to \(\lambda_0 - a T + \epsilon \eta\) we increase the climate feedback factor and hence lower the equilibrium climate sensitivity.

[9]:

# NBVAL_IGNORE_OUTPUT
dl_values = np.array([10 ** 3, 10 ** 4, 10 ** 5, 10 ** 15]) * ur("m")
dl_values

[9]:
Magnitude
[            1000            10000           100000 1000000000000000]
Unitsmeter
[10]:

# NBVAL_IGNORE_OUTPUT
runner = TwoLayerModel()
output = []
equivalent_parameters = []
for dl in tqdman.tqdm(dl_values, desc="dl values", leave=False):
    runner.dl = dl
    output.append(runner.run_scenarios(driver))
    equivalent_parameters.append(({"two-layer deep ocean depth": runner.dl}, runner.get_impulse_response_parameters()))

output = run_append(output)
output.head()

[10]:
time 1850-01-01 00:00:00 1851-01-01 00:00:00 1852-01-01 00:00:00 1853-01-01 00:00:00 1854-01-01 00:00:00 1855-01-01 00:00:00 1856-01-01 00:00:00 1857-01-01 00:00:00 1858-01-01 00:00:00 1859-01-01 00:00:00 ... 3840-01-01 00:00:00 3841-01-01 00:00:00 3842-01-01 00:00:00 3843-01-01 00:00:00 3844-01-01 00:00:00 3845-01-01 00:00:00 3846-01-01 00:00:00 3847-01-01 00:00:00 3848-01-01 00:00:00 3849-01-01 00:00:00
a (watt / delta_degree_Celsius ** 2 / meter ** 2) climate_model dl (meter) du (meter) efficacy (dimensionless) eta (watt / delta_degree_Celsius / meter ** 2) lambda0 (watt / delta_degree_Celsius / meter ** 2) model region run_idx scenario unit variable
0.0 two_layer 1000 50 1.0 0.8 1.246667 idealised World 0 1pctCO2 W/m^2 Effective Radiative Forcing 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 4.000000 4.000000 4.000000 4.000000 4.000000 4.000000 4.000000 4.000000 4.000000 4.000000
delta_degC Surface Temperature|Upper 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 3.207635 3.207638 3.207642 3.207645 3.207648 3.207652 3.207655 3.207658 3.207661 3.207665
Surface Temperature|Lower 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 3.206227 3.206236 3.206244 3.206252 3.206261 3.206269 3.206278 3.206286 3.206294 3.206302
W/m^2 Heat Uptake 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.001153 0.001149 0.001144 0.001140 0.001136 0.001132 0.001128 0.001124 0.001120 0.001115
10000 50 1.0 0.8 1.246667 idealised World 0 1pctCO2 W/m^2 Effective Radiative Forcing 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 4.000000 4.000000 4.000000 4.000000 4.000000 4.000000 4.000000 4.000000 4.000000 4.000000

5 rows × 2000 columns

As we can see below, as the deep ocean becomes deeper and deeper, its equivalent timescale increases. This demonstrates that the deep ocean is becoming increasingly close to being an infinite reservoir.

[11]:

for v in equivalent_parameters:
    v[1]["d1"] = v[1]["d1"].to("yr")
    v[1]["d2"] = v[1]["d2"].to("yr")

equivalent_parameters

[11]:

[({'two-layer deep ocean depth': 1000 <Unit('meter')>},
  {'d1': 3.211845269334279 <Unit('a')>,
   'd2': 273.9854219906419 <Unit('a')>,
   'q1': 0.4810875417166762 <Unit('delta_degree_Celsius * meter ** 2 / watt')>,
   'q2': 0.32105149571648217 <Unit('delta_degree_Celsius * meter ** 2 / watt')>,
   'efficacy': 1.0 <Unit('dimensionless')>}),
 ({'two-layer deep ocean depth': 10000 <Unit('meter')>},
  {'d1': 3.2342013046041056 <Unit('a')>,
   'd2': 2720.915300585768 <Unit('a')>,
   'q1': 0.4878523568580023 <Unit('delta_degree_Celsius * meter ** 2 / watt')>,
   'q2': 0.3142866805751838 <Unit('delta_degree_Celsius * meter ** 2 / watt')>,
   'efficacy': 1.0 <Unit('dimensionless')>}),
 ({'two-layer deep ocean depth': 100000 <Unit('meter')>},
  {'d1': 3.2364276918618766 <Unit('a')>,
   'd2': 27190.435420502465 <Unit('a')>,
   'q1': 0.4885246922441889 <Unit('delta_degree_Celsius * meter ** 2 / watt')>,
   'q2': 0.31361434518885845 <Unit('delta_degree_Celsius * meter ** 2 / watt')>,
   'efficacy': 1.0 <Unit('dimensionless')>}),
 ({'two-layer deep ocean depth': 1000000000000000 <Unit('meter')>},
  {'d1': 3.238959349323281 <Unit('a')>,
   'd2': 271883581625601.66 <Unit('a')>,
   'q1': 0.4889441930744013 <Unit('delta_degree_Celsius * meter ** 2 / watt')>,
   'q2': 0.31555972744518723 <Unit('delta_degree_Celsius * meter ** 2 / watt')>,
   'efficacy': 1.0 <Unit('dimensionless')>})]

As shown in the plots below, as the deep ocean becomes bigger, it can uptake more heat and hence mixed-layer warming is reduced. However, whilst the deep ocean is finite, the mixed-layer warming does eventually reach the same equilibrium (independent of deep ocean depth), it just takes longer to do so. Once the deep ocean becomes infinite, as discussed above, we effectively have a single-layer model with an increased climate feedback factor (lower equilibrium climate sensitivity) which can uptake heat forever.

[12]:

# NBVAL_IGNORE_OUTPUT
scenario_to_plot = "1pctCO2"
xlim = [1850, 3500]
pkwargs = dict(
    hue="dl (meter)",
    style="variable",
    time_axis="year"
)
fig = plt.figure(figsize=(9, 9))

ax = fig.add_subplot(211)
output.filter(scenario=scenario_to_plot, variable="Surface Temperature|Upper").lineplot(**pkwargs, ax=ax)

ax = fig.add_subplot(212, sharex=ax)
output.filter(scenario=scenario_to_plot, variable="Heat Uptake").lineplot(**pkwargs, ax=ax)
ax.set_xlim(xlim)

[12]:

(1850.0, 3500.0)
_images/usage_one-layer-model_19_1.png

Development

If you’re interested in contributing to OpenSCM Two Layer Model, we’d love to have you on board! This section of the docs details how to get setup to contribute and how best to communicate.

Contributing

All contributions are welcome, some possible suggestions include:

  • tutorials (or support questions which, once solved, result in a new tutorial :D)

  • blog posts

  • improving the documentation

  • bug reports

  • feature requests

  • pull requests

Please report issues or discuss feature requests in the OpenSCM Two Layer Model issue tracker. If your issue is a feature request or a bug, please use the templates available, otherwise, simply open a normal issue.

As a contributor, please follow a couple of conventions:

  • Create issues in the OpenSCM Two Layer Model issue tracker for changes and enhancements, this ensures that everyone in the community has a chance to comment

  • Be welcoming to newcomers and encourage diverse new contributors from all backgrounds: see the Python Community Code of Conduct

  • Only push to your own branches, this allows people to force push to their own branches as they need without fear or causing others headaches

  • Start all pull requests as draft pull requests and only mark them as ready for review once they’ve been rebased onto master, this makes it much simpler for reviewers

  • Try and make lots of small pull requests, this makes it easier for reviewers and faster for everyone as review time grows exponentially with the number of lines in a pull request

Getting setup

To get setup as a developer, we recommend the following steps (if any of these tools are unfamiliar, please see the resources we recommend in Development tools):

  1. Install conda and make

  2. Run make virtual-environment, if that fails you can try doing it manually

    1. Change your current directory to OpenSCM Two Layer Model’s root directory (i.e. the one which contains README.rst), cd openscm-twolayermodel

    2. Create a virtual environment to use with OpenSCM Two Layer Model python3 -m venv venv

    3. Activate your virtual environment source ./venv/bin/activate

    4. Upgrade pip pip install --upgrade pip

    5. Install the development dependencies (very important, make sure your virtual environment is active before doing this) pip install -e .[dev]

  3. Make sure the tests pass by running make checks, if that fails the commands can be read out of the Makefile

Getting help

Whilst developing, unexpected things can go wrong (that’s why it’s called ‘developing’, if we knew what we were doing, it would already be ‘developed’). Normally, the fastest way to solve an issue is to contact us via the issue tracker. The other option is to debug yourself. For this purpose, we provide a list of the tools we use during our development as starting points for your search to find what has gone wrong.

Development tools

This list of development tools is what we rely on to develop OpenSCM Two Layer Model reliably and reproducibly. It gives you a few starting points in case things do go inexplicably wrong and you want to work out why. We include links with each of these tools to starting points that we think are useful, in case you want to learn more.

Other tools

We also use some other tools which aren’t necessarily the most familiar. Here we provide a list of these along with useful resources.

  • Regular expressions

    • we use regex101.com to help us write and check our regular expressions, make sure the language is set to Python to make your life easy!

Formatting

To help us focus on what the code does, not how it looks, we use a couple of automatic formatting tools. These automatically format the code for us and tell use where the errors are. To use them, after setting yourself up (see Getting setup), simply run make format. Note that make format can only be run if you have committed all your work i.e. your working directory is ‘clean’. This restriction is made to ensure that you don’t format code without being able to undo it, just in case something goes wrong.

Buiding the docs

After setting yourself up (see Getting setup), building the docs is as simple as running make docs (note, run make -B docs to force the docs to rebuild and ignore make when it says ‘… index.html is up to date’). This will build the docs for you. You can preview them by opening docs/build/html/index.html in a browser.

For documentation we use Sphinx. To get ourselves started with Sphinx, we started with this example then used Sphinx’s getting started guide.

Gotchas

To get Sphinx to generate pdfs (rarely worth the hassle), you require Latexmk. On a Mac this can be installed with sudo tlmgr install latexmk. You will most likely also need to install some other packages (if you don’t have the full distribution). You can check which package contains any missing files with tlmgr search --global --file [filename]. You can then install the packages with sudo tlmgr install [package].

Docstring style

For our docstrings we use numpy style docstrings. For more information on these, here is the full guide and the quick reference we also use.

Releasing

First step

  1. Test installation with dependencies make test-install

  2. Update CHANGELOG.rst

    • add a header for the new version between master and the latest bullet point

    • this should leave the section underneath the master header empty

  3. git add .

  4. git commit -m "Prepare for release of vX.Y.Z"

  5. git tag vX.Y.Z

  6. Test version updated as intended with make test-install

Push to repository

To do the release, push the tags and the repository state.

  1. git push

  2. git push --tags

Assuming all the checks pass, this automatically triggers a release on PyPI via the .github/workflows/ci-cd-workflow.yml action.

Why is there a Makefile in a pure Python repository?

Whilst it may not be standard practice, a Makefile is a simple way to automate general setup (environment setup in particular). Hence we have one here which basically acts as a notes file for how to do all those little jobs which we often forget e.g. setting up environments, running tests (and making sure we’re in the right environment), building docs, setting up auxillary bits and pieces.

Base API

Module containing the base for model implementations

class openscm_twolayermodel.base.Model

Bases: abc.ABC

Base class for model implementations

reset()

Reset everything so that a new run can be performed.

Called as late as possible before run().

run()

Run the model.

abstract set_drivers(*args, **kwargs)

Set the model’s drivers

step()

Do a single time step.

class openscm_twolayermodel.base.TwoLayerVariant

Bases: openscm_twolayermodel.base.Model

Base for variations of implementations of the two-layer model

property delta_t

pint.Quantity Time step for forward-differencing approximation

property erf

pint.Quantity Effective radiative forcing

reset()

Reset everything so that a new run can be performed.

Called as late as possible before run().

run()

Run the model.

run_scenarios(scenarios, driver_var='Effective Radiative Forcing', progress=True)

Run scenarios.

The model timestep is automatically adjusted based on the timestep used in scenarios. The timestep used in scenarios must be constant because this implementation has a constant timestep. Pull requests to upgrade the implementation to support variable timesteps are welcome https://github.com/openscm/openscm-twolayermodel/pulls.

Parameters

  • scenarios (ScmDataFrame or ScmRun or pyam.IamDataFrame or pd.DataFrame or np.ndarray or str) – Scenarios to run. The input will be converted to an ScmRun before the run takes place.

  • driver_var (str) – The variable in scenarios to use as the driver of the model

  • progress (bool) – Whether to display a progress bar

Returns

Results of the run (including drivers)

Return type

ScmRun

Raises

ValueError – No data is available for driver_var in the "World" region in scenarios.

set_drivers(erf)

Set drivers for a model run

Parameters

erf (pint.Quantity) – Effective radiative forcing (W/m^2) to use to drive the model

Raises

AssertionError – erf is not one-dimensional

step()

Do a single time step.

Impulse Response Model API

Module containing the impulse response model

The 2-timescale impulse response model is mathematically equivalent to the two-layer model without state dependence.

class openscm_twolayermodel.impulse_response_model.ImpulseResponseModel(q1=<Quantity(0.3, 'delta_degree_Celsius * meter ** 2 / watt')>, q2=<Quantity(0.4, 'delta_degree_Celsius * meter ** 2 / watt')>, d1=<Quantity(9.0, 'a')>, d2=<Quantity(400.0, 'a')>, efficacy=<Quantity(1.0, 'dimensionless')>, delta_t=<Quantity(0.0833333333, 'a')>)

Bases: openscm_twolayermodel.base.TwoLayerVariant

TODO: top line and paper references

This implementation uses a forward-differencing approach. This means that temperature and ocean heat uptake values are start of timestep values. For example, temperature[i] is only affected by drivers from the i-1 timestep. In practice, this means that the first temperature and ocean heat uptake values will always be zero and the last value in the input drivers has no effect on model output.

property d1

pint.Quantity Response timescale of first box

property d2

pint.Quantity Response timescale of second box

property delta_t

pint.Quantity Time step for forward-differencing approximation

property efficacy

pint.Quantity Efficacy factor

property erf

pint.Quantity Effective radiative forcing

get_two_layer_parameters()

Get equivalent two-layer model parameters

For details on how the equivalence is calculated, please see the notebook impulse-response-equivalence.ipynb in the OpenSCM Two Layer model repository.

Returns

dict of str – Input arguments to initialise an openscm_twolayermodel.TwoLayerModel with the same temperature response as self

Return type

pint.Quantity

property q1

pint.Quantity Sensitivity of first box response to radiative forcing

property q2

pint.Quantity Sensitivity of second box response to radiative forcing

reset()

Reset everything so that a new run can be performed.

Called as late as possible before run().

run()

Run the model.

run_scenarios(scenarios, driver_var='Effective Radiative Forcing', progress=True)

Run scenarios.

The model timestep is automatically adjusted based on the timestep used in scenarios. The timestep used in scenarios must be constant because this implementation has a constant timestep. Pull requests to upgrade the implementation to support variable timesteps are welcome https://github.com/openscm/openscm-twolayermodel/pulls.

Parameters

  • scenarios (ScmDataFrame or ScmRun or pyam.IamDataFrame or pd.DataFrame or np.ndarray or str) – Scenarios to run. The input will be converted to an ScmRun before the run takes place.

  • driver_var (str) – The variable in scenarios to use as the driver of the model

  • progress (bool) – Whether to display a progress bar

Returns

Results of the run (including drivers)

Return type

ScmRun

Raises

ValueError – No data is available for driver_var in the "World" region in scenarios.

set_drivers(erf)

Set drivers for a model run

Parameters

erf (pint.Quantity) – Effective radiative forcing (W/m^2) to use to drive the model

Raises

AssertionError – erf is not one-dimensional

step()

Do a single time step.

Two Layer Model API

Module containing the two-layer model

class openscm_twolayermodel.two_layer_model.TwoLayerModel(du=<Quantity(50, 'meter')>, dl=<Quantity(1200, 'meter')>, lambda0=<Quantity(1.24666667, 'watt / delta_degree_Celsius / meter ** 2')>, a=<Quantity(0.0, 'watt / delta_degree_Celsius ** 2 / meter ** 2')>, efficacy=<Quantity(1.0, 'dimensionless')>, eta=<Quantity(0.8, 'watt / delta_degree_Celsius / meter ** 2')>, delta_t=<Quantity(31557600.0, 'second')>)

Bases: openscm_twolayermodel.base.TwoLayerVariant

TODO: top line and paper references

This implementation uses a forward-differencing approach. This means that temperature and ocean heat uptake values are start of timestep values. For example, temperature[i] is only affected by drivers from the i-1 timestep. In practice, this means that the first temperature and ocean heat uptake values will always be zero and the last value in the input drivers has no effect on model output.

property a

pint.Quantity Dependence of climate feedback factor on temperature

property delta_t

pint.Quantity Time step for forward-differencing approximation

property dl

pint.Quantity Depth of lower layer

property du

pint.Quantity Depth of upper layer

property efficacy

pint.Quantity Efficacy factor

property erf

pint.Quantity Effective radiative forcing

property eta

pint.Quantity Heat transport efficiency

get_impulse_response_parameters()

Get equivalent two-timescale impulse response model parameters

For details on how the equivalence is calculated, please see the notebook impulse-response-equivalence.ipynb in the OpenSCM Two Layer model repository.

Returns

dict of str – Input arguments to initialise an openscm_twolayermodel.ImpulseResponseModel with the same temperature response as self

Return type

pint.Quantity

Raises

ValueError – self.a is non-zero, the two-timescale model does not support state-dependence.

property heat_capacity_lower

pint.Quantity Heat capacity of lower layer

property heat_capacity_upper

pint.Quantity Heat capacity of upper layer

property lambda0

pint.Quantity Initial climate feedback factor

reset()

Reset everything so that a new run can be performed.

Called as late as possible before run().

run()

Run the model.

run_scenarios(scenarios, driver_var='Effective Radiative Forcing', progress=True)

Run scenarios.

The model timestep is automatically adjusted based on the timestep used in scenarios. The timestep used in scenarios must be constant because this implementation has a constant timestep. Pull requests to upgrade the implementation to support variable timesteps are welcome https://github.com/openscm/openscm-twolayermodel/pulls.

Parameters

  • scenarios (ScmDataFrame or ScmRun or pyam.IamDataFrame or pd.DataFrame or np.ndarray or str) – Scenarios to run. The input will be converted to an ScmRun before the run takes place.

  • driver_var (str) – The variable in scenarios to use as the driver of the model

  • progress (bool) – Whether to display a progress bar

Returns

Results of the run (including drivers)

Return type

ScmRun

Raises

ValueError – No data is available for driver_var in the "World" region in scenarios.

set_drivers(erf)

Set drivers for a model run

Parameters

erf (pint.Quantity) – Effective radiative forcing (W/m^2) to use to drive the model

Raises

AssertionError – erf is not one-dimensional

step()

Do a single time step.

Constants API

Physical constants used in calculations

openscm_twolayermodel.constants.DENSITY_WATER = <Quantity(1000.0, 'kilogram / meter ** 3')>

density of water

Type

pint.Quantity

openscm_twolayermodel.constants.HEAT_CAPACITY_WATER = <Quantity(4181.0, 'joule / delta_degree_Celsius / kilogram')>

heat capacity of water

Type

pint.Quantity

Errors API

Exceptions raised within openscm_twolayermodel

exception openscm_twolayermodel.errors.ModelStateError

Bases: ValueError

Exception raised if a model’s state is incompatible with the action

args
with_traceback()

Exception.with_traceback(tb) – set self.__traceback__ to tb and return self.

exception openscm_twolayermodel.errors.UnitError

Bases: ValueError

Exception raised if something has the wrong units

args
with_traceback()

Exception.with_traceback(tb) – set self.__traceback__ to tb and return self.

Utils API

Utility functions

openscm_twolayermodel.utils.convert_lambda_to_ecs(lambda_val, f2x=<Quantity(3.74, 'watt / meter ** 2')>)

Convert a lambda value to equilibrium climate sensitivity (ECS)

Parameters

  • lambda_val (pint.Quantity) – Value of lambda to convert to ECS

  • f2x (pint.Quantity) – Value of the forcing due to a doubling of atmospheric CO2 to assume during the conversion

Returns

ECS value

Return type

pint.Quantity

Raises

TypeError – lambda_val or f2x is not a pint.Quantity.

Changelog

All notable changes to this project will be documented in this file.

The format is based on Keep a Changelog, and this project adheres to Semantic Versioning.

The changes listed in this file are categorised as follows:

  • Added: new features

  • Changed: changes in existing functionality

  • Deprecated: soon-to-be removed features

  • Removed: now removed features

  • Fixed: any bug fixes

  • Security: in case of vulnerabilities.

v0.2.3 - 2021-04-27

Fixed

v0.2.2 - 2021-04-27

Added

  • (#33) Information in README and testing for conda install

Changed

  • (#32) Include LICENSE, README.rst and CHANGELOG in package

  • (#30) Require scmdata>=0.9

  • (#27) Fixed the discussion (in the relevant notebook) of how a one-layer model can be made from the two-layer implementation here

Fixed

  • (#30) Incorrect call to scmdata.ScmRun() in tests

v0.2.1 - 2020-12-23

Added

Changed

  • (#26) Updated to new scmdata version (and hence new openscm-units API)

  • (#25) JOSS paper following JOSS review 1

  • (#23) Moved notebooks into full documentation following JOSS review (closes #17)

  • (#21) Quoted pip install instructions to ensure cross-shell compatibility following JOSS review (closes #16)

  • (#20) Option to remove tqdm progress bar by passing progress=False

v0.2.0 - 2020-10-09

Added

  • (#7) JOSS paper draft

Changed

  • (#7) Require scmdata>=0.7

v0.1.2 - 2020-31-07

Changed

  • (#12) Upgrade to scmdata>=0.6.2 so that package can be installed

v0.1.1 - 2020-06-29

Added

  • (#8) Add notebook showing how to run a single-layer model

Changed

  • (#11) Re-wrote the getting started notebook

  • (#10) Re-wrote CHANGELOG

  • (#9) Update to scmdata 0.5.Y

v0.1.0 - 2020-05-15

Added

  • (#3) Add first implementation of the models

  • (#1) Setup repository