Hello, If anyone calculated global mean surface temperture anomaly caused by volcanic ERF, give me a hand please.

We want to calculate surface temperature anomaly caused by volcanic aerosol forcing. We have the montly ERF anomaly (i.e., relative to novolcano simulation). How should we calculate the surface temperature anomaly caused by the ERF anomaly but relative to a reference state?

Below is my code as my unserstanding. I don't know if it is right. Also, this temperature anomaly is from a reference state, I want to get a temperature anomaly only from ERF anomaly (i.e., volcanic minus novolcanic simulation). What should I do?

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

def yyyymm(ys,ye,dtype=int): year = np.arange(ys, ye+1, 1, dtype = dtype) month = np.arange(1,13,1) ym = [] for iy in year: for im in month: ym.append(iy*100 + im) return ym

def cal_temp_ano(scenarios, configs, time, df_volcanic): f = FAIR() f.define_time(time[0], time[1], time[2])

f.timebounds[:5], f.timebounds[-5:]
f.timepoints[:5], f.timepoints[-5:]

f.define_configs(configs)
f.define_scenarios(scenarios)

# 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)
species, properties = read_properties()
f.define_species(species, properties)
f.allocate()
f.emissions
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')
f.fill_from_rcmip()

volcanic_forcing = np.zeros(f.timebounds.shape)
volcanic_forcing[:] = df_volcanic[time[0]:time[1]].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)
df = pd.read_csv("./tests/test_data/4xCO2_cummins_ebm3.csv")
models = df['model'].unique()

seed = 30

# remember we rescale sigma
if (np.array(np.array(configs).shape) > 1):
    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] / np.sqrt(f.timestep),
             config=config)
        fill(f.climate_configs['sigma_xi'], df.loc[condition, 'sigma_xi'].values[0] / np.sqrt(f.timestep),
             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)
else:
    config = configs
    condition = (df['model'] == 'CESM2-WACCM-FV2') & (df['run'] == 'r1i1p1f1')
    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] / np.sqrt(f.timestep),
         config=config)
    fill(f.climate_configs['sigma_xi'], df.loc[condition, 'sigma_xi'].values[0] / np.sqrt(f.timestep),
         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)

f.run()

tempano = f.temperature.loc[dict(scenario='ssp585', layer=0)]
forcano = f.forcing.loc[dict(scenario='ssp585', specie='Volcanic')]
return tempano, forcano,

============== user defined ================

ys = 1971 ye = 1993

ym = np.array(yyyymm(ys, ye, dtype=int))

df_volcanic = pd.read_csv('./tests/test_data/volcanic_ERF-PINATUBO.csv', index_col='year') temp_volc, forc_volc = cal_temp_ano(['ssp585'], ['CESM2-WACCM-FV2'] ,[ys, ye+0.917, 1/12], df_volcanic=df_volcanic)

OMS-NetZero / FAIR

how to calculate GMST anomaly from volcanic ERF #128

============== user defined ================