Canopy aerodynamic scaling fixes

[ckoven]

We've been using the first-PFT values for three parameters related to the aerodynamic calculations: leaf characteristic diameter, displacement height, and canopy roughness length. this does a more sensible weighting.

Description:

per #690, we've been doing a weird thing for some of the scaling-up of turbulence parameters to the patch level. this code does a more proper scaling of the three boundary conditions. for leaf characteristic dimension, we calculate a weighted average, weighted by the total leaf plus stem of all cohorts. for displacement height and roughness length, we weight by the crown area of only cohorts that are in the canopy and multiply by that cohort's height.

fixes #690.

Collaborators:

@rgknox @rosiealice

Expectation of Answer Changes:

because the values of the three parameters are currently all the same, the scaling to of parameters from PFTs to cohorts shouldn't matter. however, because the displacement height and roughness length are being scaled per cohort rather than by the tallest cohort in the patch, this will change answers.

in discussion with @rgknox and @rosiealice we also though that dimensioning the displacement height, and canopy roughness length parameters by PFT doesn't necessarily make sense anyway, and possibly we should make both of these PFT-independent scalars. if so, we may do that as a subsequent change, possibly alongside a more in-depth reworking of how we calculate these quantities.

Checklist:

• My change requires a change to the documentation.
• I have updated the in-code documentation .AND. (the technical note .OR. the wiki) accordingly.
• I have read the CONTRIBUTING document.
• FATES PASS/FAIL regression tests were run
• If answers were expected to change, evaluation was performed and provided

Test Results:

not yet tested formally.

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

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

FATES baseline hash-tag:

Test Output:

[rosiealice]

Why does this part use the stem area for the weighting function? Is that to avoid strange things happening in winter when there are no leaves?

[ckoven]

was following this suggestion here: #690 (comment). but as it should in principle be the LAI rather than the LAI+SAI, and in principle if there is no LAI on the patch then that should be handled by the HLM robustly anyway, perhaps it is better to remove the LAI from the calculation entirely. relatedly: probably we should be doing an ELAI weighting, though that probably involves new science and would need more site-level testing.

[rosiealice]

Did you mean remove LAI, or SAI? I think @rgknox has/had a point in his comment though, that there is still drag if there are no leaves. And rb is used in canopy conductance, i.e. in the efpot line:
https://github.com/ESCOMP/CTSM/blob/6fad071dce82ab4cbf42f6cd47399bd463bd126a/src/biogeophys/CanopyFluxesMod.F90#L1137

So I think we need some sort of value for rb relevant to the stem-only canopy?

[ckoven]

yes, sorry, i meant remove sai. since we're only using the lai (or lai+sai) weighting for dleaf and not the other two terms, it should be safe to remove it, since if there's no leaf there wouldn't be an effect of leaf diameter?

[ckoven]

i see, sorry, you're right. dleaf does go into rb so needs to have a valid value even if lai = 0. so I'll leave as-is then?

[rosiealice]

This was waiting on me to say yes I think, so yes, that's fine...

[rgknox]

@ckoven I notice that we calculate update %c_area here:

https://github.com/NGEET/fates/pull/735/files#diff-315c6d4746f5566a97a47acd7630dc085c4535691dd54fbcbcd2b4f1aa5f8700R1950-R1951

But we use it right before the update here:
https://github.com/NGEET/fates/pull/735/files#diff-315c6d4746f5566a97a47acd7630dc085c4535691dd54fbcbcd2b4f1aa5f8700R1935

[rgknox]

do treelai and treesai need to be re-calculated here?

[ckoven]

Hi All,
A brief update on this. Was curious how the change to using crown-area-weighted-average canopy height rather than maximum-cohort height affected the results, so I did a comparison of four runs: (1) current main branch, (2) this branch, (3) main branch but in FATES-SP mode, and (4) main branch but in "classic" (i.e. big-leaf) SP mode. Plots of the displacement height (variable "DISPLA") are below (results for Z0M are for the most part proportionally identical). For these, I ran for 60 years in 4x5 mode for the first two, and then restarted and ran for five years with the two different codebases, in the full-fates modes, and for the SP modes, I just ran from bare ground for five years.

Summary is that this code does reduce a high bias in displacement height as compared to the SP mode, but a high bias remains as compared to SP mode, especially in the mid- and high-latitude forests. Would be good to benchmark and calibrate against actual data, as in ESCOMP/CTSM#1316, and add capability to use a more sophisticated turbulence scaling model. Both SP modes produce the same results (which is good as they should). But I feel like this PR is (slightly) better than the existing code and should be merged.

1. current main branch:
2. this PR branch:
3. main branch but in FATES-SP mode:
4. main branch but in "classic" SP mode:

[rgknox]

Do want to have a low bias relative to SP? I suppose if we have a lower bias compared to SP, than parameterizations generated through SP are more translatable over to FATES simulations, likewise FATES simulations are more comparable to SP simulations?

[ckoven]

@rgknox its a good question. In principle, the SP roughness height is qualitatively different from the full-FATES roughness height, because the SP is just an input dataset and based on a canopy of uniformly-sized trees, whereas the full-FATES is an emergent outcome. I guess, if we take the input dataset used by SP as being a good metric of the actual canopy height, then we do want them to be close, and we want to adjust the full-FATEs to be more like SP. So I think one interpretation is that we need to change the default height allometry for the non-tropical forest trees to be more representative of those forest types, and which should therefore reduce the bias outside of the tropics. But, yeah, I guess that for the SP-based calibration cascade to work as well as possible, then we do need to match the the aerodynamic properties between the two modes.

[rosiealice]

So the FATES-SP mode is in fact looking the same as the big leaf SP mode? (that's good to know!)
Agree that this does seem to indicate that we need to adjust the allometry of extra-tropical trees.

[glemieux]

Testing is partially complete; due to a known MOSART issue, the single site tests need to be modified, rebuilt and rerun by hand.

All tests returned DIFFs as expected, but two tests failed RUN:
ERP_D_P32x2_Ld3.f19_g17.I2000Clm50FatesCru.cheyenne_intel.clm-FatesColdDef and ERS_Ld30.f45_f45_mg37.I2000Clm50FatesCru.cheyenne_intel.clm-FatesReducedComplexSatellitePhenology

Currently investigating. The test folder is located here on Cheyenne: /glade/u/home/glemieux/scratch/ctsm-tests/tests_canopy_scaling_fix

[glemieux]

All single site tests were manually resubmitted and all returned DIFFs as expected.

I manually resubmitted the ERP_D_P32x2_Ld3 test since the failure mode was on a deallocation of a patch during patch fusion, which struck me as an odd failure mode. For reference, the output just prior to the error line was: 0:MCT::m_Router::initp_: GSMap indices not increasing...Will correct. That said, upon re-submission, this test ran successfully to completion. I'm not sure what this 'instability' suggests.

Reviewing the ERS satellite phenology mode test outputs, I noticed that the run is failing after the restart. We know that the satellite phenology mode has an issue with restarts, so I think this likely just a variation on that failure mode. It should be noted that the SMS version of the satellite phenology mode test runs successfully to completion.