Potential issue with SV (vol inputs)

[ptrscll]

Describe the bug
I was practicing changing parameters in Hector and found that when I decreased the volcanic scaling factor from 1.0 to 0.5, temperatures and CO2 concentrations decreased in the future and across most of the 1900s.

To Reproduce
In the R interface, I tried running the ssp245 scenario and comparing the results using the default parameters with the results after changing only VOLCANIC_SCALE() from 1.0 to 0.5. Graphs below were produced using a date range of 1745-2100 and depicting changes in RF_TOTAL(), GLOBAL_TAS(), VOLCANIC_SO2(), and CONCENTRATIONS_CO2().

Expected behavior
I expected temperatures and CO2 concentrations to remain the same across both scenarios in the future since, to my understanding, Hector doesn't model any volcanic eruptions after the present-day. I also expected temperatures and CO2 concentrations to be higher and radiative forcing during eruptions to be less negative for the volscl=0.5 scenario in the past

Actual behavior
Temperatures and CO2 concentrations behaved mostly as expected before ~1850 but were lower for the volscl=0.5 scenario for most times afterwards. Additionally, radiative forcing during volcanic eruptions appears to have been more negative post-1850 for the volscl=0.5 scenario than volscl=1.0 scenario. Graphs showing these variables and volcanic SO2 are shown below:

Runtime and build environment (please complete the following information):

• OS: Windows

SessionInfo() output:
R version 4.4.0 (2024-04-24 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 10 x64 (build 19045)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.utf8 LC_CTYPE=English_United States.utf8 LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C LC_TIME=English_United States.utf8

time zone: America/New_York
tzcode source: internal

attached base packages:
[1] stats graphics grDevices utils datasets methods base

other attached packages:
[1] ggplot2_3.5.1 hector_3.2.0

loaded via a namespace (and not attached):
[1] vctrs_0.6.5 cli_3.6.2 rlang_1.1.3 generics_0.1.3 textshaping_0.3.7 glue_1.7.0 labeling_0.4.3
[8] colorspace_2.1-0 ragg_1.3.2 scales_1.3.0 fansi_1.0.6 grid_4.4.0 munsell_0.5.1 tibble_3.2.1
[15] lifecycle_1.0.4 compiler_4.4.0 dplyr_1.1.4 Rcpp_1.0.12 pkgconfig_2.0.3 rstudioapi_0.16.0 farver_2.1.2
[22] systemfonts_1.1.0 R6_2.5.1 tidyselect_1.2.1 utf8_1.2.4 pillar_1.9.0 magrittr_2.0.3 tools_4.4.0
[29] withr_3.0.0 gtable_0.3.5

[bpbond]

This is super clear and useful @ptrscll — thank you! Definitely seems like there's a problem.

[bpbond]

This is happening because all the ini files have this line:

[forcing]
baseyear=1750		; when to start reporting; by definition, all F=0 in this year

When the model hits baseyear, it stores the current forcing values, and then for every subsequent year, subtracts the baseyear values from the calculated ones--i.e., it normalized
radiative forcing relative to the 1750 values.

But...volcanic forcing has a value (0.185746) in that year, and so that's being subtracted from every subsequent value, and so we get this for the current main (with the red dashed line showing where zero should be):

This is an honest-to-goodness, nontrivial bug. Well done @ptrscll !

NB this problem would also apply to RF_misc but, interestingly, it looks like someone thought of that already:

;RF_misc=csv:tables/ssp245_emiss-constraints_rf.csv
RF_misc[1750]=0

[bpbond]

@kdorheim Effect of setting 1750 volcanic forcing value to zero, so that subsequent normalization is correct:

[kdorheim]

Okay yes I think that is totally part of the problem however I think the warming effect seen in RF_vol in the rf_vol_fix seems strange... @ssmithClimate and @claudiatebaldi do you have any insights here?

RF_misc is default set to 0 for the entire run so it is not like that was set up properly either 😬 I guess how do we want to handle this/prevent it from happening again? Add a check when reading in those inputs?

[bpbond]

I think the broader problem is that we have two classes of radiative forcings: human-driven, that we want to scale from 1750 (or whenever), and natural changes (volcanoes, misc) that we don't want to scale--they're absolutes.

[claudiatebaldi]

CMIP6 prescribed a small constant volcanic forcing during piControl. And the same constant is replicated (after a ten year ramp from the end of the historical) for the future scenarios.

[kdorheim]

Okay! So here is some more info on the RF_Vol time series we are using (zenodo archive & GitHub)

Why is volcanic forcing not zero in 1750?
We use a baseline of the last 2500 years of stratospheric aerosol optical depth (SAOD) from Toohey & Sigl (2017), from which volcanic forcing is calculated using a linear relationship. We define zero forcing to be the long-term mean SAOD in this time series. Because of periodic volcanic eruptions that produce spikes of strong negative forcing, the volcanic forcing in quiescent years is a small positive number to make the time average zero. This is done to ensure that long-term equilibrium temperauture anomalies are zero when running with only volcanic forcing. It also ensures a seamless transition to the SSPs, where future volcanic forcing is defined to be zero (i.e. equal to the long term average).

If we used the reference for zero forcing to be either zero SAOD or the background SAOD, we would end up with a spurious long-term cooling trend when this time series is applied to an energy balance model.

[bpbond]

That's extremely helpful--thank you!

So, how to handle the issue of positive forcing in quiescent years? The average of the 1745-2023 input time series we're using is -0.041 W/m2, a small cooling effect overall. So we're not getting the "long-term equilibrium temperauture anomalies are zero when running with only volcanic forcing"!

Options that occur to me:

• Leave it as is
• Accept the above logic, and re-center the time series so that its average is zero (by slightly increasing the quiescent year values, presumably)
• Zero out all positive values; define our time series as being the negative excursions only, and accept that there will be a non-smooth transition to the zero-volcano future

[ssmithClimate]

Shouldn't we should follow the CMIP protocol, using the historical average as given by our time series and scaling factor. That would be applied in the starting year, and an all future years. So that would be something calculated in code.

I believe it's pretty clear we should be doing that for the future - there will be volcanic, forcing in the future, even though we don't know what it is exactly. As for the starting point, I'm a little less clear about that.

Note that there could also be a reference period effect in the output as well, if looking at output that is relative to an average of reference years (1850-1900?) instead of the raw Hector values.

[ssmithClimate]

Seems @kdorheim's post above gives us guidance on the outcome we would want - with only volcanic forcing applied - wouldn't we would want hector to have a long-term zero temperature change (regardless of volcanic scaler chosen by the user) when run well into the future?

[kdorheim]

CMIP6 prescribed a small constant volcanic forcing during piControl. And the same constant is replicated (after a ten year ramp from the end of the historical) for the future scenarios.

Looking at RF_vol plot above it does seems like the its current implementation "main" is consistent with the CMIP6 protocol. Is my interpretation correct @claudiatebaldi?

[claudiatebaldi]

That is what it would appear to me too...

[bpbond]

I'm getting a little confused.

I believe it's pretty clear we should be [applying the historical average] for the future

with only volcanic forcing applied - wouldn't we would want hector to have a long-term zero temperature change (regardless of volcanic scaler chosen by the user) when run well into the future?

Aren't these two things inconsistent?

Looking at RF_vol plot above it does seems like the its current implementation "main" is consistent with the CMIP6 protocol.

Agreed.

[ssmithClimate]

I don't believe these are inconsistent. If the hector temperature "spin up" takes into account (explicitly or otherwise) the average historical volcanic forcing already, then if we did an experiment that kept that average forcing constant over all time - that should result in no temperature signal relative to pre-industrial conditions. I believe that's what we want.

Then the historical eruptions would cause negative temperatures temporarily, but temperatures should ultimately come back to zero over time, particularly after the historical period.