Update default allometry parameters for tree PFTs

[JessicaNeedham]

This PR updates the default global tree PFT allometry parameters for height, crown area, and aboveground biomass.

The current defaults are based on studies with fairly small sample sizes compared to the data now available. It would be nice to have updated versions ahead of no comp calibration. This first pass at updates uses data from BAAD and Tallo databases. The full analysis is here.

Things to note:

• Tallo includes BAAD and has a lot of additional data. But BAAD has a deciduous v evergreen label which we need. I therefore cross reference by species name so we can use the extra Tallo data with the BAAD PFT designations.
• I’m using the Martinez-Cano et al. 2019 height allometry, a power function for crown area, and the Chave et al. 2014 aboveground biomass allometry for all tree PFTs.
• We could use the Chave et al. 2014 height allometry but it depends on the environmental factor E which would need to either be read in for each grid cell or calculated from forcing data.
• The Chave et al. 2014 paper is based on more data than is in BAAD and it therefore makes sense to use parameters from that paper for tropical PFTs AGB allometries.
• For wood density I found the mean value per PFT using @mpaiao's analysis of TRY plus other data sources.
• Crown area fits are surprisingly linear even when data is binned by dbh to address issues of heteroskedasticity.
• Some PFTs contain data from many species with quite different allometries. Ideally we would do some kind of weighting based on real abundances/spatial distributions of species, rather than have the fits be based on biased sampling. For now, I have not tried to correct this. If someone wants to do a more in depth analysis that would be great!
• I’ve set the fates_allom_dbh_maxheight parameter to a large number since the Martinez-Cano function is a Michaelis-Menten and asymptotes and this parameter is therefore redundant (and potentially dangerous!).
• I have not looked at shrubs or grasses. I believe updates for those PFTs are coming from other team members.

A simulation comparing the default allometries with these updates is here. Results are means over 10 years following 100 years of spin up. This is a nocomp simulation run at 4x5 grid. The parameter files are identical apart from allometry parameters, but include Atkin respiration, updates to leaf reflectance parameters, slatop, vcmax and a few other changes compared to the default parameter file. Probably it needs to run for longer than 100 years so I’ll keep updating this notebook.

Collaborators:

Analysis of BAAD data is based on work by @adrifoster. Wood density data from @mpaiao. Conversations with @pollybuotte , @jenniferholm , @ckoven , @rgknox , @glemieux, @mpaiao and @adrifoster.

Expectation of Answer Changes:

This will be answer changing for anyone running with the default parameter file. For runs with custom parameter files it should be bfb. Initial tests show higher biomass of the needleleaf evergreen extratropical tree, and broadleaf evergreen tropical tree, and lower biomass of the broadleaf cold deciduous extra tropical tree. I'll run an ILAMB comparison once it has spun up longer.

Checklist

If this is your first time contributing, please read the CONTRIBUTING document.

All checklist items must be checked to enable merging this pull request:

Contributor

• The in-code documentation has been updated with descriptive comments
• The documentation has been assessed to determine if updates are necessary

Integrator

• FATES PASS/FAIL regression tests were run
• Evaluation of test results for answer changes was performed and results provided

Documentation

• Technical Note update:
• User's Guide update:

This PR does not require a change to the tech note or user's guide.

Test Results:

CTSM (or) E3SM (specify which) test hash-tag:

CTSM (or) E3SM (specify which) baseline hash-tag:

FATES baseline hash-tag:

Test Output:

https://github.com/JessicaNeedham/FATES-plotting-scripts/blob/main/Allometry_test.ipynb

[JessicaNeedham]

FATES nocomp runs with default and updated parameters have run for 250 years now. The ILAMB comparison is here. The updated parameters do better across all the ILAMB biomass datasets. There's a lot of variation within the biomass datasets though.

Some figures of global veg C and LAI by PFT are here. LAI is higher for the needleleaf evergreen extra tropical PFT and the broadleaf cold deciduous extra tropical tree but not that different for any of the others. Veg C is higher for the needleleaf evergreen extratropical PFT and lower for some of the other PFTs.

Comparisons with plot data from Galbraith et al. 2013 and Piponiot et al. 2022 are here*. The Piponiot et al. data set has AGB, aboveground woody productivity (AWP) and aboveground woody mortality (AWM) by dbh class.. The size distributions don’t look great in the high latitude sites but some tropical sites look good. The Galbraith et al. dataset is total AGB, AWP and carbon residence time across sites. FATES AGB, and AWP are generally too low. Carbon residence time is much too low in FATES and doesn’t vary much across sites (~20-40 years), whereas there is significant variation in the data (~20 - 125 years).

*Three sites (Luquillo, Palamanui and Laupahoehoe) don’t have FATES comparisons because they are islands that don’t appear to have vegetation in 4 x 5 degree runs

[rgknox]

@JessicaNeedham , in the ilamb comparisons, can you confirm that "allom_CA_h_agb" is the set that is proposed in the PR?

[JessicaNeedham]

@JessicaNeedham , in the ilamb comparisons, can you confirm that "allom_CA_h_agb" is the set that is proposed in the PR?

Yep, "allom_CA_h_agb" has all the changes to CA, height and AGB proposed in this PR. The others are showing the default allometry parameters, and the incremental change of just changing CA and changing CA and height. Note though that while these four parameter sets differ only in allometry parameters, they all differ from the default parameter file currently on main in quite a few other parameters that have been updated during calibration - leaf reflectance, leaf longevity, slatop etc.

This ILAMB comparison will be updated as the runs continue - they are currently at about 170 years.
Most recent update is here: https://portal.nersc.gov/cfs/e3sm/jneedham/_build_allom_v2/

[JessicaNeedham]

This PR updates CA by changing the parameter that sets the difference between the leaf biomass and CA exponents - which means LAI will no longer be constant with size. An alternative would be to change the leaf biomass exponent and keep the difference parameter at 0. I'm not sure if we have good empirical data on leaf biomass dbh relationships which could inform this. In the meantime, I'm going to test the alternative and compare them in ILAMB. Results to follow... @rgknox and @glemieux can we hold off merging until I've done this test?

[adamhb]

Hi @JessicaNeedham,
Great work on this! I just wanted to bring up one related point for what its worth. When leaf biomass becomes very large in relation to crown area, LAI/VAI can become so large that it triggers a VAI bin error when the number of required VAI bins becomes larger than the default cap (30). This occurs most readily for tall skinny crowns (e.g. conifers). I'm using the dbh max param to cap leaf biomass in big trees to ensure that the VAI stays below a critical level. I just wanted to bring it up as a way that the dbh max param was being used to avoid other issues, but it doesn't look like you're getting rid of it entirely right?
Adam

[JessicaNeedham]

Thanks @adamhb that's a great point. I wasn't going to get rid of the dbh_maxheight parameter in this PR - just set it very high so it becomes moot. But I was literally just writing up another issue to suggest removing it entirely because it can cause strange behaviour - e.g. a sudden acceleration in radial growth at large sizes.

The max vai bins being hard coded at 30 seems like a problem (also see issue #1081). I've been manually changing nlevleaf to larger values to get around those errors. Are you using the default vai parameters? The need for thinner and more numerous leaf layers should go away with the new two stream radiation scheme. But if you're getting those errors even with default vai parameters that suggests the current limit of 30 is not high enough to accommodate reasonable PFT allometries.

[adamhb]

Hi @JessicaNeedham ,
Good timing! Yes, I'm using default vai parameters, so I agree that the limit of 30 could be raised, but I'm not sure how this would affect computation time so I stayed away from raising this for now. Another option would be to have an asymptotic leaf biomass equation similar to the Martinez-Cano form you're implementing for height (then you could get rid of dbh max height). As you allude to above, the leaf biomass data is far more sparse than the height data, but I suppose fitting the same data to a different functional form would be fine right? I bet leaf biomass data is particularly limited for large trees where we would want to verify this asymptotic relationship, but it would essentially be the same assumption as having a cap on leaf biomass like what I'm doing now. Does that make sense?
Adam

[JessicaNeedham]

Hi @adamhb . The current dbh to leaf biomass function is a power law, so an asymptotic function would be a completely different relationship. It would be good to look at some data or literature on this - wouldn't any increase in size require an increase in leaves in order to meet increased maintenance costs?

I'm also a bit confused about how the dbh max height parameter caps leaf biomass. That parameter stops height growth which has the effect of increasing dbh growth - so given the power relationship between dbh and leaf biomass wouldn't the dbh max height parameter actually lead to more leaf biomass?

If the vai bins are being maxed out because leaf biomass is too high I think it makes sense to increase the number of leaf layer bins if you think the leaf biomass and crown area allometries are reasonable. Or if leaf biomass actually is too high, then you could reduce the exponent in the dbh to leaf biomass relationship (but keep the crown area allometry exponent the same).

[adamhb]

Hi @JessicaNeedham,
When allom_lmode = 3 leaf biomass is capped based on fates_allom_dbh_maxheight in this section of code. But I agree that when running without the cap then an uncapped power function would be very different from anything asymptotic. The question of whether leaves should scale forever with dbh is interesting. I think that if maintenance costs keeping increasing then your intuition that leaves should too makes sense. I wonder though if the live tissues (and therefore maintenance) can stop growing, while the tree continues to put down heartwood and therefore keeps increasing in size? Will see if I can find anything on this in the literature.
Adam

[JoshuaRady]

@JessicaNeedham,
Thank you for putting together this detailed PR. I think this looks pretty good. I'm glad you kept the old values as a comment in the parameter file. Since these changes are going to impact the majority of users I would like to see this well advertised as a release or a least a well documented tag (if we do that).

I have some questions related to the fates_allom_dbh_maxheight discussion above but I see issue #1102 and will add my comments there.

[rgknox]

Great idea @JoshuaRady , I'll send out an email to the community when we integrate

[rgknox]

This PR is pretty clean at this point. There is a clear trail of the changes that were made, and the changes have been substantiated. This PR has had plenty of visibility as well.

[JessicaNeedham]

I tried adding in leaf biomass updates but the BAAD data is very limited for some of the PFTs, especially at large sizes. I was getting very high exponents, and when I tried to do a dbh initialisation with 50 cm trees the vai bins maxed out, even with nlevleaf set to 50. This seems to be the same problem @adamhb was encountering. For now, I'm going to say this PR is good to go, with just CA, height, and agb updates. Future changes could address the leaf biomass allometries.

Updated ILAMB results are here.

[adamhb]

@JessicaNeedham. It makes sense that if your leaf biomass exponents were large that is would max out the vai bins, especially for narrow crowned trees. I agree that the leaf biomass question can be tackled in a separate PR. Nice work!

[JessicaNeedham]

I think I’m done updating these parameters now. There are still issues with leaf biomass parameters that should be addressed in a future PR. The problem is that there is very little data on leaf biomass, but by updating crown area but not leaf biomass we end up with very strange individual tree LAIs. Height, crown area and AGB parameters are now all based on data, but leaf biomass parameters I hand tuned to give what looked like reasonable individual tree LAIs but we should do a more formal calibration at some point.

Estimating parameters from data is here.

A comparison of default v updated allometries is here.

Results from a FATES run (250 years) comparing them is here.

An ILAMB comparison is here but note that runs were only 100 years in when I ran this. Perlmutter is down now but i can update the comparison with the spun up versions.

[glemieux]

Regression testing on derecho against fates-sci.1.72.3_api.34.0.0-ctsm5.1.dev174 is complete. All expected tests passed b4b, with the exception of the two stream tests. While these tests have a known issue with not passing exact restart, this is the first time I've seen them not pass baseline comparison. I'm rerunning just those tests now to see if I can replicate the issue.

I also ran the full fates regression testing suite with the parameter file updates in this PR to confirm that all expected tests would run successfully. No unexpected errors were seen. These defaults will be synchronized with an HLM-side API update at a later date via the PR corresponding to #1128.

[glemieux]

I figured out what was happening here. The two stream tests auto build the default parameter file on the fly so while the other test mods are using the current default nc file listed in the HLM, this test mod actually reflects the actual state of the updated parameter. Note that the seed dispersal tests do this as well, but those are breaking due to a different known issue (ESCOMP/CTSM#2423).

This is ready to integrate.