Increasing CO2 under 0 NBP constraint

dawnlwoodard commented 1 year ago

Describe the bug If I run Hector in equilibrium (all input emissions are 0) and with an NBP constraint set to 0, there are unexpected increases in atmospheric CO2 (and thus temperature, etc) (figure below).

When running Hector freely in equilibrium (no constraint), we see only small drift and not this larger increase (below).

This increase may seem less significant in this instance, but I suspect this same dynamic is behind much larger CO2 amplification seen when running Hector with a non-zero pulse NBP constraint (e.g. constraint of 10 PgC in 1800, 0 otherwise).

To Reproduce I made the first figure running Hector (latest v3_dev branch as of this morning) with the following R code:

library(hector)
library(dplyr)
library(ggplot2)

scenario_file <- paste0("input/hector_equilibrium_latest.ini")
ini_file <- system.file(scenario_file, package="hector")

year0 <- 1745
core <- hector::newcore(ini_file)
hector::run(core,runtodate = year0)

for (year in 1746:1900){
  setvar(core,year,"NBP_constrain",0.0,"Pg C/yr")
  hector::run(core,runtodate = year)
}

climate_output <- fetchvars(core,dates=1746:1900,vars=c(GLOBAL_TAS(),CONCENTRATIONS_CO2(),FFI_EMISSIONS(), NBP_CONSTRAIN(),NBP(),OCEAN_UPTAKE()))

climate_output %>% tidyr::pivot_wider(id_cols=c("scenario","year","variable","value"),
                                           values_from = value, names_from = variable) %>% mutate(nbp_tot=cumsum(NBP)*-1,ocean_tot=cumsum(ocean_uptake),co2_c=CO2_concentration*2.12) %>%
  tidyr::pivot_longer(cols=c("ffi_emissions","global_tas","co2_c","nbp_tot","ocean_tot","NBP_constrain","CO2_concentration","NBP","ocean_uptake"),names_to="variable",values_to = "value") ->
  climate_out

ggplot(data=climate_out, aes(x=year,y=value)) +
  geom_line() +
  facet_wrap(vars(variable),nrow=3, scales="free_y") +
  theme_classic() ->
  test_figure

Input files

The equilibrium table (below - goes in input/tables): hector_equilibrium_latest.csv

The corresponding ini file (renamed to .txt so Github will upload it, but it is a regular ini file once you change the extension back): hector_equilibrium_latest.txt

bpbond commented 1 year ago

This is interesting, and I'm not sure if we should consider it a bug or not. Note that @dawnlwoodard 's zero-NBP scenario above is weird, because LUC is a part of NBP, so ordering a LUC of 10 Pg C but also 0 NBP is nonsensical. To keep things simpler, let's say that in the year of the pulse, the NBP constraint is -10.

Scenario 1: no constraint versus same-NBP constraint

Consider Hector running unconstrained with no emissions except a LUC pulse in 1800: NPP, RH, and NBP are all flat (because the model is in steady state) until the pulse happens:

This pushes up atmospheric CO2, which raises ocean uptake
NPP drops because of the vegetation loss (#651 ), followed by RH
The model proceeds to a new equilibrium: slightly lower NPP and RH and higher CO2

Reassuringly, if we feed in a NBP constraint identical to the NBP values of the no-constraint run, things look good:

Scenario 2: NBP constraint of zero (except for the pulse)

The 1800 LUC pulse happens.

When under a constraint, NPP and RH fluctuate slightly as the model tries to respond to transient conditions
The green line below is basically Dawn's scenario above: the model gets into a runaway state, driven by an interaction of CO2 fertilization and the constraint-following behavior (see below)
If β is set to zero (blue line), we get better performance

Scenario 3: change constraint behavior

Right now, when the model sees a mismatch between its computed land-atmosphere flux (i.e. NBP) and the ordered NBP constraint, it (i) down/up regulates NPP and RH as needed, and (ii) adjusts the end-of-timestep land and atmosphere C pools correspondingly.

In the LUC pulse scenario, this has the effect of raising atmospheric CO2: the model has computed a CO2-fertilized NPP and thus a positive NBP, but because it's been ordered to hit zero, it lowers NPP, raises RH, and raises atmospheric CO2 (because what was projected to be a land C sink has vanished due to the constraint). The higher CO2 raises NPP in the next time step, etc., leading to the runaway CO2.

Another option would be to take any constraint-driven C change and send it to the deep ocean, as we do when there's a CO2 constraint. This also improves things; NPP and RH are still drifting down, but the changes are small, and the atmosphere is OK.

bpbond commented 1 year ago

Summary

Basically, in this issue we're ordering the model to do several contradictory things simultaneously: constrain to zero NBP, fertilize NPP, and have a LUC pulse.

I see the logic of changing from "if there's a constraint, adjust the atmosphere by any differential" to "if there's a constraint, send any atmosphere differential to the deep ocean". That's consistent with the CO2 constraint approach and will ensure consistency with the solver's end-of-timestep values.

But @dawnlwoodard seems like the simplest thing for you to do is to turn off CO2 fertilization when you're holding an NBP constraint.

nbp_zero.csv

dawnlwoodard commented 1 year ago

Sorry I believe I may have caused some confusion here. If there was a 10 Pg C pulse in my input file that was a mistake on my part. It should be 0 all across the board. I was going back and forth with what values I was experimenting with so I easily could have uploaded the wrong file by mistake. But yeah the figures I showed were from running with no pulse and no emissions. Apologies about that!

dawnlwoodard commented 1 year ago

If I turn off CO2 fertilization & RH and run the same test, I end up with a smaller increase but it's still there:

dawnlwoodard commented 1 year ago

It may not seem significant on this scale but that exact pattern exists (and is amplified substantially, even with CO2 fertilization off) when running with a non-zero NBP constraint (still with emissions left at 0). For example, see the difference between a pulse test of LUC emissions set to 10 Pg C (no constraint on NBP):

and a pulse test of NBP constraint set to 10 Pg C (no LUC emissions) with beta=0 and Q10 = 1:

dawnlwoodard commented 1 year ago

Per our conversation this morning, this may be able to be fixed by shunting the offset carbon to the deep ocean instead of the atmosphere. @bpbond

bpbond commented 1 year ago

@dawnlwoodard Sorry for the slow response. I can confirm that running with a deep ocean shunt fixes the issue:

JGCRI / hector