Adjust potential NPP based on vegetation LUC losses (or gains)

[bpbond]

Per discussion with @dawnlwoodard and @kdorheim yesterday — see #638

Changes

This PR adds a new term to our NPP calculation: vegetation losses (or gains) due to LUC change, with the the aim of avoiding double-counting regrowth fluxes. We save the end-of-spinup vegetation C value; keep a running total of LUC emissions from (or gains to) vegetation; and then adjust NPP:

  npp_luc_adjust = (end_of_spinup_vegc - cum_luc_va) / end_of_spinup_vegc;
  npp = npp * npp_luc_adjust;

Impact (SSP2.45)

> library(hector)
> core <- newcore(system.file("input/hector_ssp245.ini", package = "hector"))
> run(core)
> x <- fetchvars(core, dates = 1750:2300)
> x$which <- "v3_dev"
> [ same with new functionality ]
> y$which <- "This PR"
> z <- rbind(x, y)
> ggplot(z, aes(year, value, color = which)) + facet_wrap(~variable, scales="free") + geom_point() + coord_cartesian(xlim=c(1950, 2030))

[kdorheim]

Cool, I am really interested to hear what @dawnlwoodard says since their work identified this problem.

Something I am slightly confused about though is how the carbon pools are connected with the NPP.

hector/src/simpleNbox-runtime.cpp

Line 353 in d1af1b7

// Update soil, detritus, and atmosphere pools - luc fluxes

Based on my limited understanding of Hector's carbon cycle component it looks like the c pools (veg, detritus, and soil) are updated with luc emissions already right? So does this change now cause a double counting of the luc in these pools? Also does beta not have an effect on NPP anymore?

[kdorheim]

Using the same setup as Ben but focusing on the historical period but including more variables

[bpbond]

Cool—thanks @kdorheim

[bpbond]

it looks like the c pools (veg, detritus, and soil) are updated with luc emissions already right? So does this change now cause a double counting of the luc in these pools? Also does beta not have an effect on NPP anymore?

Old sequence:

1. Calculate beta effect
2. Apply beta effect to NPP
3. Allocate NPP to pools and proceed through carbon cycle

New sequence:

1. (Keep track of the running total of vegetation C lost to LUC)
2. Calculate beta effect
3. Calculate LUC effect as NPP0 * (change of veg_lost_to_luc / preindustrial veg). E.g. if 1% of vegC has been lost to date to LUC, the LUC effect would be 99%
4. Apply beta and LUC effects to NPP: npp = npp0 * beta_effect * luc_effect basically
5. Allocate NPP to pools and proceed through carbon cycle

[dawnlwoodard]

Yay!! Thanks @bpbond! I will test this out tomorrow morning and plan to add the regrowth figures here so we can see how much of an impact this has.

@kdorheim Yes the carbon pools are already updated with LUC emissions, but, unlike respiration, NPP does not itself depend on the size of the current carbon pools. It is purely calculated as a function of initial NPP and a CO2 fertilization multiplier. So this PR actually allows NPP to adjust with LUC emissions. Which, logically, should make some sense. If you cut down half of the forest, the NPP in that forest should not be the exact same as before. What we are trying to do here is something like this:

where if Hector is in equilibrium and then you remove an amount of LUC emissions, Hector's land carbon drops to a new equilibrium. Of course the example I'm showing above doesn't have any nonlinear components, so Hector is not going to look this clean. And to maintain perfect equilibrium you also need to remove those LUC emissions proportional to the current size of each pool, but as Ben pointed out when we discussed this issue, that's not realistically how carbon would move, so we're also not going to do that.

Also, beta should still affect NPP - this is just an additional multiplier.

But at first glance this looks promising - as in, it does what we expect of increasing CO2 and dropping the size of the land carbon pools.

[ssmithClimate]

Note that if you cut down half a forest, NPP doesn't necessarily change. It's just a different vegetation type. Any change would be because the new vegetation type has a lower or higher NPP than the old one. It could go either way. You need to look at the data to see which way that might go.

But perhaps you mean that the NPP of the forest goes down?

[bpbond]

@ssmithClimate thanks, and right, it definitely can get tricky—NPP might not drop. But even if it didn't, the duplication-of-regrowth issue, along with the pave-the-amazon scenario (cf. #531 ), plus the fact that Hector already runs low compared to historical CO2/temp, all argue for addressing this I think.

[dawnlwoodard]

Okay, so running pulse tests with this is interesting.

If we turn off nonlinear land components (by setting Q10 to 1.0 and beta to 0.0) and otherwise run Hector in emissions-driven mode with all emissions zeroed out except a single LUC emissions pulse, we get what I would expect in vegetation and detritus, but soil still has an odd amount of regrowth happening and I'm not sure why:

If we don't turn off the nonlinear pieces, things get extra wonky:

Like, what is happening with that sharp spike in NPP/detritus?? And the vegetation has a weird regrowth pattern.

[kdorheim]

Hmmm @dawnlwoodard agreed that the soil C pool dynamics look funky @bpbond any ideas as to where the growth in soil C is coming from?

[dawnlwoodard]

As I'm looking into this further, another piece I don't understand is why the pools each have 2 large drops (in 1800 and 1801) rather than just one in 1800 when the pulse happens.

[bpbond]

@dawnlwoodard

soil still has an odd amount of regrowth happening

Probably because we're not taking out the LUC pulse proportionately? I.e., we're still following Hector's default (more realistic, I said to you last week) fractions.

another piece I don't understand is why the pools each have 2 large drops (in 1800 and 1801)

I suspect that this is a time-smudge issue: half is being taken out in 1800-1800.5, and half in 1800.5-1801, and because "1800.0" is January 1, 1800, it gets charged to two separate years. But yeah that's not what people will expect...

Seems like I should make the first change, at least, so you have a clean basis for testing.

[bpbond]

@dawnlwoodard Temporarily shifted to proportional LUC losses/gains in 3b40525 ; haven't tested except yes it runs, but should do what you want.

[dawnlwoodard]

Thanks, @bpbond! So running with the changes and with nonlinear pieces turned off we would expect to produce a clean drop to a new equilibrium with no regrowth. But instead it now has consistent, small regrowth in every land C pool. It's a small amount, so not super concerned, but I am curious what's driving it. With the beta effect off, there's no actual dependence on the atmospheric carbon concentration, so it's not a response to that pool changing.

With nonlinear back on (below) the jerkiness is fixed but the regrowth is still amplified somewhat (which we would expect). The question is just whether this is too much regrowth or not. Which, as @bpbond and I discussed when we first talked about this, is not necessarily a clear line to draw. My thought is that this is maybe small enough to not worry about if we can fix the regrowth that's happening in the purely linear one ^. Because this amplification doesn't seem like much if that weren't already happening. Mind you, none of what I've said here is necessarily a reason to not merge this in, as long as we can at least explain why that bit of regrowth is happening still in the purely linear version above. Like, less double counting is still an improvement. But imho if we're going to make a change this significant (in terms of impact on model results, necessitating recalibration, etc), then it does seem better to get it all the way how we want it.

[dawnlwoodard]

My other thought on this is that before you changed the luc fractions, @bpbond, the soil vs vegetation were not actually behaving in a more realistic way, unless I'm not thinking about this right. Based on what you and I discussed, what we talked about as being more realistic is that the soil should lose relatively less carbon due to LUC emissions, compared to vegetation and detritus. So what you'd expect is the soil carbon pool to drop further over time as it comes back into equilibrium, right? Whereas instead we saw the biggest drop in soil, such that it was regrowing, and vegetation and detritus were the ones dropping further to come into equilibrium. Am I missing something? Is that actually the behavior we'd expect?

[bpbond]

Fixed every test issue except for the old-new test; that'll be the last change, when everything else is resolved.

Thanks @dawnlwoodard — interesting to see the graphs above! On the one hand, like you I am curious as to why, with concentration-driven mode and no beta and Q10=1, why there's adjustment over time in the carbon pools after your LUC pulse. I would guess that it takes time for the NPP-LUC effect to work its way through the system in some way we're not expecting.

On the other, as you note, it's a small adjustment: <0.1%. Tough for me to prioritize this.

So what you'd expect is the soil carbon pool to drop further over time as it comes back into equilibrium, right? Whereas instead we saw the biggest drop in soil, such that it was regrowing, and vegetation and detritus were the ones dropping further to come into equilibrium. Am I missing something? Is that actually the behavior we'd expect?

Hmm, let's touch base about this in person? I'm finding it hard to walk this through in my head.

[dawnlwoodard]

Sounds good, @bpbond. Let's discuss Wednesday.

Re: "On the other, as you note, it's a small adjustment: <0.1%. Tough for me to prioritize this." - yes, the remaining effect here is small and alone I agree would not be a priority. But the intended default way we'd like to run Hector is of course with feedbacks on, which regrows a larger fraction of the carbon lost, and likely also with LUC fractions fixed, which regrows even more. I'm not sure we can just dismiss that amount of regrowth as insignificant. And it's unclear to me whether the fixes are unrelated or not. But, as I said before, less regrowth is still better. So we could certainly justify letting this be a good start and leaving it at that.

[dawnlwoodard]

@bpbond reminder to circle back to this. Per our discussion, we believe we are pulling too much carbon from all three land C pools relative to new NPP. Possibly an issue with calculations that happen outside of the solver versus fractions that get calculated in the solver.

[bpbond]

@dawnlwoodard OK, I have fixed and merged ( #644 ) ANOTHER bug that your work has uncovered! 👏 The LUC pulse test now matches your spreadsheet I think — see graph in #643 .

@kdorheim This PR now has some significant changes that, combined, definitely affect model behavior:

• Interpolation bugfix
• LUC loss fractions changed to proportional
• NPP adjustment for vegetation LUC losses

At your convenience, could you take a hard look at the performance of this PR, e.g. by running it through your v3 assessment RMarkdown you showed me yesterday?

[bpbond]

Two other notes before I go offline for a while.

1. The old-new test is currently disabled, just so we can quickly and easily see if everything else passes.
2. What about the LUC fractions?

LUC fractions

To date Hector has let users set how LUC emissions (uptake) pull from (flow into) veg, detritus, and soil:

f_lucv=0.1			; Fraction of land use change flux from vegetation
f_lucd=0.01			; Fraction of land use change flux from detritus (balance from soil)

But I'm persuaded by Dawn's argument here: this approach is both ad hoc and dependent on what kind of LUC occurs, and it's simpler and more defensible to just pull proportionately from the three pools. That what was implemented in 3b40525 on a test basis, and I'm proposing we make permanent. Two questions then:

• Does this seem reasonable?
• If we do this, do we rip out the f_lucv and f_lucd settings entirely? Or let users set them if they want to override the default behavior?

[kdorheim]

@bpbond will do taking a look at it now!

[dawnlwoodard]

@bpbond so exciting!! So glad you figured out what was behind the weird behavior I was seeing and that these tests can help with future diagnosis as well.

Re: LUC fractions, I think it could be useful to have users be able to override them especially if they're biome-specific. Then someone could break out a particular biome that should have specific LUC-related behavior if they want. But if it's going to be a pain to keep them around, may not be worth it.

[kdorheim]

Okay I ran it through the manuscript three calibration/figures we were looking at yesterday.

Look at how great the Hector v3 with the new luc update fits [co2]!

The comparison with the other observations GISS temp and GCP are not as good but I am not that concerned.

[dawnlwoodard]

@kdorheim glad CO2 is fitting better! Also, on these figures, is old luc just Hector v3 without Ben's latest fixes? I guess I'm particularly wondering why old and new luc have such different land sink behavior.

[kdorheim]

@dawnlwoodard so there are two differences in the Hector output, the hand full of Ben's latest commits and also the parameters used. The difference between the land sink is a combination of these two factors. In the plot below I ran the same luc branch but with different parameter combinations and we see the large change in land sink behavior.

old params

     diff    q10_rh      beta 
1.7263138 0.9330690 0.2890119 

new params

     diff    q10_rh      beta 
1.1759621 2.1558709 0.5529147 

[bpbond]

I took an exploratory look at keeping f_lucv and f_lucd around — see the luc_option branch.

Making these fractions biome-specific, and allowing for two different behaviors (default proportional, but user could override on a biome-specific basis), that has to be implemented in two different places (calc_derivs() and stashCValues()), would introduce a fair amount of complexity. And imho that's not something we need more of at the moment, at least not without a strong reason!

So, my vote would be to remove them entirely, which is what #645 does. We can discuss next week, but if this makes sense, @kdorheim can merge #645 into this PR and then re-evaluate changed from v3_dev please.

🙏

[kdorheim]

@bpbond in your figure generated above comparing Hector's outputs could you please include NBP, NPP, and Rh?

[bpbond]

@kdorheim

[kdorheim]

Thanks Ben I will take a look at the PR in a bit. The terrestrial carbon cycle changes seem large do these changes seem reasonable/along the lines of what you were expecting?

[bpbond]

I don't know. Those NBP lines are pretty distressingly different.

[kdorheim]

Yea... it seems like these changes have made Hector more sensitive to parameter values... What does the output looks like when
aero = 1
volscl = 1
ECS = 3
NPP0 = 56.2
diff = 1.18
q10 = 2.16
beta = 0.55

[bpbond]

[kdorheim]

hmmm interesting the developments seemed to have really increased Hector's parameter sensitivity

[kdorheim]

@bpbond this is only a partial review but may be relevant to the interpolation devlopment

[kdorheim]

Why start at 1745? Is this a new requirement based on your recent developments?

[bpbond]

The model starts at 1745 🤷 otherwise no reason

[kdorheim]

hahah that is a pretty good reason, I will update hectordata to make sure the input files start at year 1745

[kdorheim]

@bpbond did you intend to make these changes in the NAMESPACE file? I get the impression that it is a consequence of changing the formatting in the src/rcpp_constants.cpp. Anyhow the overall effect removes all of these helper functions from the R code base which is not something I think we want to do. Unless do you have a reason why we would want to drop these functions then we should probably talk about it.

[bpbond]

Uh no whoops

[bpbond]

Thank you

[bpbond]

Superseded by #651