[WIP] Add NVP (Non-Vascular Plant) PFT in CTSM–FATES

[huitang-earth]

Description:

This PR introduces a physically explicit representation of non-vascular plants (NVP; e.g., moss and lichen) in CTSM–FATES and its biophysical effect by assigning them to a dedicated vertical layer in CTSM (ESCOMP/CTSM#3927), and a NVP PFT in FATES. The core architectural decision is to use layer 0 of the soil-snow layer structure to accommodate NVP and its growth (i.e., changes in height) when it is present. This layer integrates consistently into CLM’s radiation (under snow), thermal, and hydrology schemes, while coupling to FATES for canopy radiation (without snow) and photosynthesis (see the following call graph).
Note: The implementation is not fully ready yet and is still open for discussion.

Key implementation details in FATES:

• Introduces NVP allometric subroutine linking cohort-level NVP LAI changes with site(column)-level NVP layer thickness (i.e., height).
• Implements independent NVP radiation, create new absorptance terms (fabd/fabi_nvp_pa) passed to CLM.
• Adds NVP-specific photosynthesis (disabling stomatal conductance and root-dependent processes) driven by NVP conditions (e.g., t_nvp_col, fwet_nvp_col, and absorbed PAR by nvp).

Collaborators:

@rosiealice

Expectation of Answer Changes:

This PR is expected to produce answer changes when NVP is active and NVP layer is present:

• Surface energy balance modified due to new top boundary condition at layer 0 (CLM, FATES)
• Soil temperature profiles altered due to NVP thermal properties (CLM)
• Hydrology affected by additional water storage and evaporation pathway (CLM)
• Carbon fluxes (GPP/NPP) changed due to NVP photosynthesis and modified PAR environment (FATES)

No answer changes are expected when NVP is disabled.

Description of generative AI usage (as necessary)

Claude code has been used for porting changes into the new FATES version and designing the overall structure. The call graph is also generated by Claude code.

Checklist

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
• Describe use of generative AI (if necessary)

Integrator

• FATES PASS/FAIL regression tests were run
• Evaluation of test results for answer changes was performed and results provided
• FATES-CLM6 Code Freeze: satellite phenology regression tests are b4b

Documentation

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

Test Results:

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

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

FATES baseline hash-tag: (TBD)

Test Output: (To be provided)

[samsrabin]

@huitang-earth This is great to see! I'm a software engineer on CTSM and was about to get started in earnest on similar work with @adrifoster, so I might be able to work with you on getting this in. I'll email about scheduling a call. In the meantime, do you have any results you could add here?

[rgknox]

Congratulations on this submission @huitang-earth ! Very exciting stuff! Added myself as a reviewer

[samsrabin]

@huitang-earth I've had a read through this and it looks great. I have a lot of questions and comments, but they can mostly wait until after our call Tuesday. One burning one, though—did you add a moss PFT? I don't see one in the updated parameter file, but I can't figure out how this would work without one. If you have a new parameter file with a moss PFT, could you add it to this PR?

[huitang-earth]

@samsrabin Thanks for your comments. Just a quick answer to your question. I indeed have added moss and lichen as a pft in the parameter file in the old version, with which I have made plenty simulations. But now I am wondering if it is a better idea to have a separate parameter file for NVP? The reason behind this thought is that the traits used for describing NVP are quite different from that for vascular PFT. If we use the same parameter file, there will be a lot of redundant parameters not needed by NVP, while some new parameters (not needed by vascular PFT) will be added in the same file. I am worried that this will make the parameter file quite messy in the future.

I am happy to learn more about your opinions before making such an implementation in the PR.

[samsrabin]

You're right that it will be messy, but I think for now it'd be best to include NVPs in the one existing parameter file. This is analogous to how we have woody and non-woody vascular plants in there already, even though there are parameters used by only one or the other. Well, I guess it's analogous as long as there's ANY overlap—is there?

Either way, I think this is the way to go for now. Our priority should be to get a complete working version so that scientific testing can be carried out with this branch; we can decide later if we want to add the extra infrastructure that would be required for a separate parameter file.

[rgknox]

Speaking specifically to the engineering and coding aspects of this PR, I see two paths that we can follow here. One path could potentially add a lot of technical debt. In the other path we use this set of changes to get even more organized about how to handle different types of plants in different parts of the code. I would also encourage us as we are refining this code going forward, to think about how it will work alongside other axis of differentiation of plants, that we have and maybe don't have, such as woody/non-wood, free-standing/lianas, etc.

[huitang-earth]

You're right that it will be messy, but I think for now it'd be best to include NVPs in the one existing parameter file. This is analogous to how we have woody and non-woody vascular plants in there already, even though there are parameters used by only one or the other. Well, I guess it's analogous as long as there's ANY overlap—is there?

Either way, I think this is the way to go for now. Our priority should be to get a complete working version so that scientific testing can be carried out with this branch; we can decide later if we want to add the extra infrastructure that would be required for a separate parameter file.

I have added NVP (non_vascular_phototroph) into the default parameter json file with some changes in the parameter values when it is relevant. The absent of root for NVP is also added now.

[huitang-earth]

Hi, I have been testing NVP code with our site. Now, we are starting to produce some sensible results. See attached figures. But there are remaining soil energy balance errors (errsoi) which I think is related the surface tiling approach, i.e., the interaction between NVP, soil, snow, and surface water (h2osfc). This is the most tricky part in the work.

I summarized the current implementation below, and hope to discuss with wider community and improve this part so that we could further reduce the errors. @samsrabin @adrifoster @rosiealice @rgknox @wwieder

NVP (moss/lichen) surface restructure: tiled energy split, solar partitioning, and h2osfc

---

1. Layer geometry and notation

The NVP layer occupies a dedicated vertical slot j = 0, between the snow and the soil:

atmosphere
  snow layers      j = snl+1 .. -1     (only when snow present, snl <= -2 with NVP)
  NVP (moss) layer j = 0               (jbot_sno == -1 marks this slot as moss)
  soil layers      j = 1 .. nlevgrnd

Surface fractions. Each column's ground is split into four tiles whose flux-basis weights
sum to 1:

fse              = frac_sno_eff                                   snow-covered fraction (flux partitioning)
frac_nvp_eff     = min(1 - fh2o - fse, max(0, frac_nvp - fse))    exposed moss      
frac_soil        = (1 - fse) - frac_nvp_eff - fh2o              bare soil           
fh2o             = frac_h2osfc         surface-water fraction

fse + frac_nvp_eff + frac_soil + fh2o = 1

Here frac_nvp = col%frac_nvp is the FATES-derived moss ground cover. frac_nvp_eff is what remains exposed after snow and water.

Two exposure bases. Surface water enter the radiation differently from its surface fluxes. fh2o is not considered for solar radiation absorption.

frac_nvp_eff     = min(1 - fh2o - fse, max(0, frac_nvp - fse))
frac_nvp_eff_solar = min(1 - fse, max(0, frac_nvp - fse))     (moss, no-fh2o)
nvp_exp            = frac_nvp_eff       / frac_nvp            (exposed moss / total moss, flux)
nvp_exp_solar      = frac_nvp_eff_solar / frac_nvp            (exposed moss / total moss, solar)

Why two bases? h2osfc affects turbulent/LW fluxes (it is its own tile) but for solar
it is part of the ground albedo, so the solar accounting is on a no-fh2o basis and the
under-water solar is redistributed to the water tile in the deposit step (§3).

---

2. Solar-absorption restructure

2.1 The 3-tile solar ABSORPTION (no-fh2o)

h2osfc is part of the ground albedo (implicitly), so its solar should be inside the soil/moss totals here:

nvp_exp_solar·sabg_nvp  +  (1 − fse)·sabg_soil  +  fse·sabg_snow  =  sabg(p)
Note: sabg_nvp is per-column by default
      sabg_snow is per-area(snow) by default
      sabg_soil is per-area(soil) by default

This partition is adopted in three places consistently:

carve-out (SurfaceRadiationMod):  sabg_soil = [ sabg(p) − fse·sabg_snow − nvp_exp_solar·sabg_nvp ] 
                                              / (1 − fse)   # (guard fse < 1, else 0)
FGR       (SoilFluxesMod):        eflx_soil_grnd += (frac_nvp_eff_solar/frac_nvp)·sabg_nvp = nvp_exp_solar·sabg_nvp
sabg_chk  (SoilTemperatureMod):   sabg_chk       += nvp_exp_solar·sabg_nvp

SABG = FGR = sabg_chk = sabg(p).

2.2 The 4-tile solar DISSIPATION

4-tile solar deposit that places each tile's solar where it can be actually dissipated:

Tile	Solar deposited into the solve	Comment
moss j=0	nvp_exp·sabg_nvp (exposed, inside hs_nvp, with thermostat) for surface; + fse·sabg_lyr_col(c,0) (buried, SNICAR) for layered temperature solver	exposed solar rides the −dhsdT thermostat
soil j=1	(1−fse−fh2o)·sabg_soil (bare + soil-under-exposed-moss) for surface; + fse·sabg_lyr_col(c,1) (under-snow, SNICAR) for layered-temperature solver
h2osfc	fh2o·sabg_soil (soil-under-water) + (nvp_exp_solar−nvp_exp)·sabg_nvp (moss-under-water) for surface only	see §3
snow	fse·sabg_snow (SNICAR, incl. buried-moss sabg_lyr(p,0)) for surface; sabg_lyr(p,-nlevel:-1) for layered-temperature solver

2.3. h2osfc treated together with soil and the NVP layer

h2osfc dissipate soil and nvp absorbed radiation through its own surface budget eflx_gnet_h2osfc (may need further improvement here!):

eflx_gnet_h2osfc (solar) = fh2o·sabg_soil                       (soil-under-water, original term)
                         + (nvp_exp_solar − nvp_exp)·sabg_nvp   (moss-under-water, the new term)

(nvp_exp_solar − nvp_exp) is exactly the moss fraction that sits under water; it is bounded by sabg_nvp/frac_nvp and vanishes when fh2o = 0.

h2osfc is thermally coupled to the soil j=1 (its fn_h2osfc uses t_soisno(c,1)), even though in an NVP column the water physically sits on the moss j=0. This is a pre-existing CLM structure we did not change; it limits how perfectly the under-water energy can balance and is one contributor to the residual errsoi errors (see below).

4. Status and known residuals

Remaining errsoi (mean ~4.3 W/m²): concentrated in partial-snow (snl<0) melt transients, largely pre-existing/structural, not from this rework.