Run the n-layer energy balance model

This notebook shows examples of extending the 3-layer energy balance model to general n.

For the two and three layer cases we’ll take the MLE estimates from Cummins et al. (2020) for HadGEM2-ES, and we’ll use the GISS forcing. Where n > 3 the data is fake.

import matplotlib.pyplot as pl
import numpy as np
import pandas as pd
import pooch

from fair.energy_balance_model import EnergyBalanceModel

df_forcing = pd.read_csv('../tests/test_data/RFMIP_ERF_tier2_GISS-E2-1-G.csv')

ebm3 = EnergyBalanceModel(
    ocean_heat_capacity=[3.62, 9.47, 98.66],
    ocean_heat_transfer=[0.54, 2.39, 0.63],
    deep_ocean_efficacy=1.59,
    gamma_autocorrelation=1.73,
    sigma_xi=0.32,
    sigma_eta=0.43,
    forcing_4co2=6.35,
    stochastic_run=True,
    seed=16
)

ebm3.add_forcing(forcing = df_forcing['GISS-E2-1-G TOT'].values, timestep=1)

ebm3.run()

ebm3.temperature

time = np.arange(1850.5, 2101)

pl.plot(time, ebm3.temperature[:,0], label='surface / top ocean layer')
pl.plot(time, ebm3.temperature[:,1], label='second ocean layer')
pl.plot(time, ebm3.temperature[:,2], label='deep ocean layer')
pl.ylabel('K relative to 1850')
pl.title('SSP2-4.5 temperature change')
pl.legend()

ebm3.emergent_parameters()
ebm3.ecs, ebm3.tcr

ebm2 = EnergyBalanceModel(
    ocean_heat_capacity=[7.73, 89.29],
    ocean_heat_transfer=[0.63, 0.52],
    deep_ocean_efficacy=1.52,
    gamma_autocorrelation=1.58,
    sigma_xi=0.64,
    sigma_eta=0.43,
    stochastic_run=True,
    forcing_4co2=6.86,
    seed=16
)

ebm2.add_forcing(forcing = df_forcing['GISS-E2-1-G TOT'].values, timestep=1)

ebm2.emergent_parameters()
ebm2.ecs, ebm2.tcr

ebm2.run()

pl.plot(time, ebm2.temperature[:,0], label='surface / top ocean layer')
pl.plot(time, ebm2.temperature[:,1], label='deep ocean layer')
pl.ylabel('K relative to 1850')
pl.title('SSP2-4.5 temperature change')
pl.legend()

# this is not based on a tuning to any existing CMIP6 model, but I try to get the TCR close to the
# HadGEM2 2- and 3-layer cases.
ebm4 = EnergyBalanceModel(
    ocean_heat_capacity=[1.3, 9, 20, 80],
    ocean_heat_transfer=[0.54, 3, 3, 0.63],
    deep_ocean_efficacy=1.2,
    gamma_autocorrelation=1.73,
    sigma_xi=0.32,
    sigma_eta=0.43,
    forcing_4co2=6.35,
    stochastic_run=True,
    seed=16
)

ebm4.emergent_parameters()
ebm4.ecs, ebm4.tcr

ebm4.add_forcing(forcing = df_forcing['GISS-E2-1-G TOT'].values, timestep=1)

ebm4.run()

pl.plot(time, ebm4.temperature[:,0], label='surface / top ocean layer')
pl.plot(time, ebm4.temperature[:,1], label='second ocean layer')
pl.plot(time, ebm4.temperature[:,2], label='third ocean layer')
pl.plot(time, ebm4.temperature[:,3], label='deep ocean layer')
pl.ylabel('K relative to 1850')
pl.title('SSP2-4.5 temperature change')
pl.legend()

# let's go totally crazy
ebm10 = EnergyBalanceModel(
    ocean_heat_capacity=[0.6, 1.3, 2, 5, 7, 10, 45, 70, 80, 130],
    ocean_heat_transfer=[0.54, 4, 5, 5, 5, 5, 5, 5, 5, 0.63],
    deep_ocean_efficacy=1.2,
    gamma_autocorrelation=1.73,
    sigma_xi=0.32,
    sigma_eta=0.43,
    forcing_4co2=6.35,
    stochastic_run=True,
    seed=16
)

ebm10.emergent_parameters()
ebm10.ecs, ebm10.tcr

ebm10.add_forcing(forcing = df_forcing['GISS-E2-1-G TOT'].values, timestep=1)

ebm10.run()

pl.plot(time, ebm10.temperature[:,0], label='surface / top ocean layer')
pl.plot(time, ebm10.temperature[:,1], label='second ocean layer')
pl.plot(time, ebm10.temperature[:,2], label='third ocean layer')
pl.plot(time, ebm10.temperature[:,3], label='fourth ocean layer')
pl.plot(time, ebm10.temperature[:,4], label='fifth ocean layer')
pl.plot(time, ebm10.temperature[:,5], label='sixth ocean layer')
pl.plot(time, ebm10.temperature[:,6], label='seventh ocean layer')
pl.plot(time, ebm10.temperature[:,7], label='eighth ocean layer')
pl.plot(time, ebm10.temperature[:,8], label='ninth ocean layer')
pl.plot(time, ebm10.temperature[:,9], label='deep ocean layer')
pl.ylabel('K relative to 1850')
pl.title('SSP2-4.5 temperature change')
pl.legend()

pl.plot(time, ebm2.temperature[:,0], label='two layer model')
pl.plot(time, ebm3.temperature[:,0], label='three layer model')
pl.plot(time, ebm4.temperature[:,0], label='four layer model')
pl.plot(time, ebm10.temperature[:,0], label='ten layer model')
pl.ylabel('K relative to 1850')
pl.title('SSP2-4.5 temperature change')
pl.legend()

pl.plot(time, ebm2.toa_imbalance, label='two layer model')
pl.plot(time, ebm3.toa_imbalance, label='three layer model')
pl.plot(time, ebm4.toa_imbalance, label='four layer model')
pl.plot(time, ebm10.toa_imbalance, label='ten layer model')
pl.ylabel('W/m2 relative to 1850')
pl.title('SSP2-4.5 TOA radiation change')
pl.legend()

Repeat everything with stochastic forcing switched off

ebm3 = EnergyBalanceModel(
    ocean_heat_capacity=[3.62, 9.47, 98.66],
    ocean_heat_transfer=[0.54, 2.39, 0.63],
    deep_ocean_efficacy=1.59,
    gamma_autocorrelation=1.73,
    sigma_xi=0.32,
    sigma_eta=0.43,
    forcing_4co2=6.35,
    stochastic_run=False,
    seed=16
)

ebm3.add_forcing(forcing = df_forcing['GISS-E2-1-G TOT'].values, timestep=1)

ebm3.run()

ebm3.temperature

time = np.arange(1850.5, 2101)

pl.plot(time, ebm3.temperature[:,0], label='surface / top ocean layer')
pl.plot(time, ebm3.temperature[:,1], label='second ocean layer')
pl.plot(time, ebm3.temperature[:,2], label='deep ocean layer')
pl.ylabel('K relative to 1850')
pl.title('SSP2-4.5 temperature change')
pl.legend()

ebm3.emergent_parameters()
ebm3.ecs, ebm3.tcr

ebm2 = EnergyBalanceModel(
    ocean_heat_capacity=[7.73, 89.29],
    ocean_heat_transfer=[0.63, 0.52],
    deep_ocean_efficacy=1.52,
    gamma_autocorrelation=1.58,
    sigma_xi=0.64,
    sigma_eta=0.43,
    stochastic_run=False,
    forcing_4co2=6.86,
    seed=16
)

ebm2.add_forcing(forcing = df_forcing['GISS-E2-1-G TOT'].values, timestep=1)

ebm2.emergent_parameters()
ebm2.ecs, ebm2.tcr

ebm2.run()

pl.plot(time, ebm2.temperature[:,0], label='surface / top ocean layer')
pl.plot(time, ebm2.temperature[:,1], label='deep ocean layer')
pl.ylabel('K relative to 1850')
pl.title('SSP2-4.5 temperature change')
pl.legend()

# this is not based on a tuning to any existing CMIP6 model, but I try to get the TCR close to the
# HadGEM2 2- and 3-layer cases.
ebm4 = EnergyBalanceModel(
    ocean_heat_capacity=[1.3, 9, 20, 80],
    ocean_heat_transfer=[0.54, 3, 3, 0.63],
    deep_ocean_efficacy=1.2,
    gamma_autocorrelation=1.73,
    sigma_xi=0.32,
    sigma_eta=0.43,
    forcing_4co2=6.35,
    stochastic_run=False,
    seed=16
)

ebm4.emergent_parameters()
ebm4.ecs, ebm4.tcr

ebm4.add_forcing(forcing = df_forcing['GISS-E2-1-G TOT'].values, timestep=1)

ebm4.run()

pl.plot(time, ebm4.temperature[:,0], label='surface / top ocean layer')
pl.plot(time, ebm4.temperature[:,1], label='second ocean layer')
pl.plot(time, ebm4.temperature[:,2], label='third ocean layer')
pl.plot(time, ebm4.temperature[:,3], label='deep ocean layer')
pl.ylabel('K relative to 1850')
pl.title('SSP2-4.5 temperature change')
pl.legend()

# let's go totally crazy
ebm10 = EnergyBalanceModel(
    ocean_heat_capacity=[0.6, 1.3, 2, 5, 7, 10, 45, 70, 80, 130],
    ocean_heat_transfer=[0.54, 4, 5, 5, 5, 5, 5, 5, 5, 0.63],
    deep_ocean_efficacy=1.2,
    gamma_autocorrelation=1.73,
    sigma_xi=0.32,
    sigma_eta=0.43,
    forcing_4co2=6.35,
    stochastic_run=False,
    seed=16
)

ebm10.emergent_parameters()
ebm10.ecs, ebm10.tcr

ebm10.add_forcing(forcing = df_forcing['GISS-E2-1-G TOT'].values, timestep=1)

ebm10.run()

pl.plot(time, ebm10.temperature[:,0], label='surface / top ocean layer')
pl.plot(time, ebm10.temperature[:,1], label='second ocean layer')
pl.plot(time, ebm10.temperature[:,2], label='third ocean layer')
pl.plot(time, ebm10.temperature[:,3], label='fourth ocean layer')
pl.plot(time, ebm10.temperature[:,4], label='fifth ocean layer')
pl.plot(time, ebm10.temperature[:,5], label='sixth ocean layer')
pl.plot(time, ebm10.temperature[:,6], label='seventh ocean layer')
pl.plot(time, ebm10.temperature[:,7], label='eighth ocean layer')
pl.plot(time, ebm10.temperature[:,8], label='ninth ocean layer')
pl.plot(time, ebm10.temperature[:,9], label='deep ocean layer')
pl.ylabel('K relative to 1850')
pl.title('SSP2-4.5 temperature change')
pl.legend()

pl.plot(time, ebm2.temperature[:,0], label='two layer model')
pl.plot(time, ebm3.temperature[:,0], label='three layer model')
pl.plot(time, ebm4.temperature[:,0], label='four layer model')
pl.plot(time, ebm10.temperature[:,0], label='ten layer model')
pl.ylabel('K relative to 1850')
pl.title('SSP2-4.5 temperature change')
pl.legend()

pl.plot(time, ebm2.toa_imbalance, label='two layer model')
pl.plot(time, ebm3.toa_imbalance, label='three layer model')
pl.plot(time, ebm4.toa_imbalance, label='four layer model')
pl.plot(time, ebm10.toa_imbalance, label='ten layer model')
pl.ylabel('W/m2 relative to 1850')
pl.title('SSP2-4.5 TOA radiation change')
pl.legend()