Bug with no permafrost run

[kdorheim]

Describe the bug
This was brought to our attention by @cahartin she was having issues with turning the permafrost module off.

To Reproduce

hector_inifile <- file.path(system.file("input", package = "hector"), "hector_ssp245.ini")
name = "No Perm"
hcore <- newcore(hector_inifile, suppresslogging = TRUE, name = name)

# Here I fetch permafrost c and it matches what’s in the ini file. Check.
fetchvars(hcore, NA, PERMAFROST_C())
# scenario year     variable value units
# 1   WPerma   NA permafrost_c   865  Pg C
# matches value in ini file

# Now I set that value to 0.
setvar(hcore, NA, PERMAFROST_C(), 0.0, "Pg C")

# I check again that the value is 0 and it is. Success?
fetchvars(hcore, NA, PERMAFROST_C())
# scenario year     variable value units
# 1  No Perm   NA permafrost_c     0  Pg C

run(hcore)

# I run hector and then pull out the Perma C value and it has carbon in it in 2100! Not what I was expecting
out <- fetchvars(hcore, 1950:2100, PERMAFROST_C())
# scenario year     variable    value units
# 1  No Perm 2100 permafrost_c 654.1661  Pg C

ggplot(data = out) +
    geom_line(aes(year, value)) + 
    labs(title = "Permafrost C Pool", 
         y = "Pg C")


# I also double checked the initial value and it’s back to the ini file value. bummer
fetchvars(hcore, 1745, PERMAFROST_C())
# scenario year     variable value units
# 1  No Perm 1745 permafrost_c   865  Pg C

Expected behavior
Setting the permafrost carbon pool to 0 should mean that the size of the permafrost pool is set at 0 at the beginning on the run and is constant throughout the run.

Actual behavior

Additional context
Hector version 3.2.0 git 1ef0b4df8563a31580ffa94eb11b113e3f317b86

[kdorheim]

Temporary Solution

Manually setting permafrost_c=0 ; permafrost, Pg C in the ini file produces the behavior we expect! See the attached folder for an example.

# LOCAL_DIR should be set to wherever the no_perm_inputs folder is located, note for right now there is only one of the 
# ini files that have set up with the non perm condition 
LOCAL_DIR <- "path/to/no_perm_inputs"  
hector_inifile <- file.path(LOCAL_DIR, "hector_ssp245_noperm.ini")
name = "No Perm"
hcore <- newcore(hector_inifile, suppresslogging = TRUE, name = name)

fetchvars(hcore, NA, vars = PERMAFROST_C())

run(hcore)

out <- fetchvars(hcore, 1900:2100, vars = PERMAFROST_C())

ggplot(data = out) +
    geom_line(aes(year, value)) + 
    labs(title = "Permafrost C Pool", 
         y = "Pg C")

This suggests that there is some problem with how the R interface is setting the initial permafrost value to 0 but not with the actual hector C++ code.

no_perm_inputs.zip

[cahartin]

Thanks @kdorheim!!

[bpbond]

I feel like this should have been caught by an automated test 😞

Thank you @cahartin ! I will fix next week, if that timing is OK for you, or @kdorheim feel free to take a look if your schedule allows.

[bpbond]

After a little investigation, this is a more general problem and has nothing to do with permafrost per se.

> hector_inifile <- file.path(system.file("input", package = "hector"), "hector_ssp245.ini")
> hcore <- newcore(hector_inifile, suppresslogging = TRUE)

> fetchvars(hcore, NA, VEG_C())
  scenario year variable value units
1 Unnamed Hector core   NA    veg_c   562  Pg C

# change the veg C value

> setvar(hcore, NA, VEG_C(), 560, "Pg C") 

> fetchvars(hcore, NA, VEG_C()) 
             scenario year variable value units
1 Unnamed Hector core   NA    veg_c   560  Pg C

# this seems okay, but now try running

> run(hcore)
Auto-resetting core to 1744
Hector core:	Unnamed Hector core
Start date:	1745
End date:	2300
Current date:	2300
Input file:	/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/hector/input/hector_ssp245.ini

> fetchvars(hcore, 1745:1746, VEG_C()) # value is unchanged
  scenario year variable    value units
1 Unnamed Hector core 1745    veg_c 562.0000  Pg C
2 Unnamed Hector core 1746    veg_c 561.9828  Pg C

😕 What's going on? In the C++ code we have a note:

      // For `veg_c`, `detritus_c`, `soil_c`, and `permafrost_c`,
      // if date is not provided, set only the "current" model pool,
      // without touching the time series variable. This is to
      // accommodate the way the INI file is parsed. For
      // interactive use, you will usually want to pass the date
      // -- otherwise, the current value will be overridden by a
      // `reset` (which includes code like `veg_c = veg_c_tv.get(t)`).

[bpbond]

Oh okay, we just need to pass a date?

> hector_inifile <- file.path(system.file("input", package = "hector"), "hector_ssp245.ini")
> hcore <- newcore(hector_inifile, suppresslogging = TRUE)

> setvar(hcore, 1745, VEG_C(), 560, "Pg C")

> run(hcore)
Auto-resetting core to 1744
Hector core:	Unnamed Hector core
Start date:	1745
End date:	2300
Current date:	2300
Input file:	/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/hector/input/hector_ssp245.ini

> fetchvars(hcore, 1745:1746, VEG_C())
             scenario year variable    value units
1 Unnamed Hector core 1745    veg_c 562.0000  Pg C
2 Unnamed Hector core 1746    veg_c 561.9828  Pg C

😞 That didn't work. It reset to 1744 (???) and so picked up the previous value. Resetting to zero doesn't work either, I tried.

What does work is to 'hack' the core by manually changing the reset date.

> hector_inifile <- file.path(system.file("input", package = "hector"), "hector_ssp245.ini")
> hcore <- newcore(hector_inifile, suppresslogging = TRUE)

> setvar(hcore, 1745, VEG_C(), 560, "Pg C")

> hcore$reset_date <- 1745

> run(hcore)
Auto-resetting core to 1745
Hector core:	Unnamed Hector core
Start date:	1745
End date:	2300
Current date:	2300
Input file:	/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/hector/input/hector_ssp245.ini

> fetchvars(hcore, 1745:1746, VEG_C())
             scenario year variable    value units
1 Unnamed Hector core 1745    veg_c 560.0000  Pg C
2 Unnamed Hector core 1746    veg_c 560.0529  Pg C

🎉 and @cahartin this trick works for your permafrost example above

[bpbond]

@kdorheim

This is terrible behavior, right? No one would ever figure it out!

At the very least, I think we need to improve documentation and print a warning if users request to change a carbon pool variable. Thoughts?

[kdorheim]

Yes that is terrible behavior! Why doesn't that happen with the NPP0?

[bpbond]

Because NPP0 isn't a time-varying state variable; it's just a static parameter.

[kdorheim]

hmm, would it make sense for veg_c and others at time 0 to also be considered a static parameter?

[bpbond]

That won't work, I don't think, because the carbon cycle has to be spun up to steady state before we begin history.

[kdorheim]

Could this at all be related to issue #735?

[bpbond]

Not the reset behavior, but the funky indexing of the time series, yeah, related. Thanks for linking.