Skip to content

Fix bugs in the Tree Recruitment Scheme with seedling dynamics #1572

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
r-ward wants to merge 20 commits into
base: main
Choose a base branch
from
+49 −23
Open

Fix bugs in the Tree Recruitment Scheme with seedling dynamics #1572

Show file tree
Hide file tree
Changes from 19 commits
Commits
Show all changes
20 commits
Select commit Hold shift + click to select a range
7d295f2
Fix PAR fraction to use relative area in SeedlingParPatch
Oct 27, 2025
2afcc3f
Skip first-day update of seedling layer smp and moisture deficit days…
Oct 27, 2025
1aa229a
fix typo
Oct 27, 2025
006e8a9
Fix initial value for sdlng_mort_par, change from temp_init_veg to in…
Oct 27, 2025
0e9f7d8
testing update to seedling layer smp
Oct 27, 2025
f2c2233
Fix model day check (must bee > 1.023 to skip 1st model day)
Oct 28, 2025
408c8cd
Add debugging blocks
Oct 29, 2025
e047aab
Fix typo
Oct 29, 2025
8477b83
Fix typo
Oct 29, 2025
5918340
Add print statements for debugging: where and how does seed_germ go n…
Nov 3, 2025
fb4f280
Update seedling layer smp and mdd running means only after first mode…
Nov 13, 2025
735c4cb
Fix negative seed pool issue by capping germination flux. Prevent see…
Nov 16, 2025
1988d06
remove debug messages for seed germ diagnostics
Dec 10, 2025
356b305
Cap seedling h2o mort rate and seed decay rate at 1, add debugging
Dec 10, 2025
c82ba29
Remove debug print statements from TRS seedling layer subroutines
r-ward Apr 27, 2026
ef6ea7e
Make rate capping style consistent in SeedGermination
r-ward Apr 28, 2026
8d7e4ca
Cap TRS seedling-to-sapling recruitment rate at 1 to prevent negative…
r-ward May 22, 2026
c4f5097
Log when TRS seedling mortality/transition rates exceed 1
r-ward May 22, 2026
204dcf4
Fix typos
r-ward May 22, 2026
dc31b22
Remove diagnostic log messages for capped TRS seedling rates
r-ward Aug 5, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
82 changes: 64 additions & 18 deletions biogeochem/EDPhysiologyMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -541,6 +541,7 @@ subroutine PreDisturbanceIntegrateLitter(currentPatch)
integer :: nlevsoil ! number of soil layers
integer :: ilyr ! soil layer loop counter
integer :: dcmpy ! decomposability index
real(r8) :: net_seed_available ! Available seed after inputs and decay

do el = 1, num_elements

Expand All @@ -550,15 +551,20 @@ subroutine PreDisturbanceIntegrateLitter(currentPatch)
! -----------------------------------------------------------------------------------

do pft = 1,numpft
litt%seed(pft) = litt%seed(pft) + &
litt%seed_in_local(pft) + &
litt%seed_in_extern(pft) - &
litt%seed_decay(pft) - &
litt%seed_germ_in(pft)

! Note that the recruitment scheme will use seed_germ
! for its construction costs.
litt%seed_germ(pft) = litt%seed_germ(pft) + &

! Calculate net available seed after inputs and decay
net_seed_available = litt%seed(pft) + litt%seed_in_local(pft) + &
litt%seed_in_extern(pft) - litt%seed_decay(pft)

! Cap germination at available seed to prevent negative seed pool
litt%seed_germ_in(pft) = min(litt%seed_germ_in(pft), net_seed_available)

@rgknox rgknox Jul 28, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@r-ward , I'm wondering if this seed_germ_in could potentially be negative. I don't think it can since seed decay is the only negative term in net_seed_available, and it is by definition a fraction of seed mass. Can you verify this?

@r-ward r-ward Aug 5, 2026

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Correct, seed_germ_in is never negative. Although, I guess it is possible for the flux to go negative if fates_recruit_seed_germination_rate or fates_trs_seedling_a_emerg were erroneously set as negative values. Maybe there should be a check on this in FatesCheckParams (I see for example a check against values of seed_dispersal_fraction < 0). I could add something similar. What do you think?


! Update pools
litt%seed(pft) = net_seed_available - litt%seed_germ_in(pft)

! Note that the recruitment scheme will use seed_germ
! for its construction costs.
litt%seed_germ(pft) = litt%seed_germ(pft) + &
litt%seed_germ_in(pft) - &
litt%seed_germ_decay(pft)

Expand Down Expand Up @@ -2253,6 +2259,7 @@ subroutine SeedDecay( litt , currentPatch )
real(r8) :: seedling_light_mort_rate ! daily seedling mortality rate from light stress
real(r8) :: seedling_h2o_mort_rate ! daily seedling mortality rate from moisture stress
real(r8) :: seedling_mdds ! moisture deficit days accumulated in the seedling layer
real(r8) :: total_seedling_mort_rate ! summed seedling mortality rates

!----------------------------------------------------------------------

Expand Down Expand Up @@ -2317,25 +2324,42 @@ subroutine SeedDecay( litt , currentPatch )

! Get the current seedling moisture deficit days (tracked as a pft-specific exponential
! average)
seedling_mdds = currentPatch%sdlng_mdd(pft)%p%GetMean()
seedling_mdds = currentPatch%sdlng_mdd(pft)%p%GetMean()

! Calculate seedling mortality as a function of moisture deficit days (mdd)
! If the seedling mmd value is below a critical threshold then moisture-based mortality is zero

if (seedling_mdds < EDPftvarcon_inst%seedling_mdd_crit(pft)) then
seedling_h2o_mort_rate = 0.0_r8
else
seedling_h2o_mort_rate = EDPftvarcon_inst%seedling_h2o_mort_a(pft) * seedling_mdds**2 + &
EDPftvarcon_inst%seedling_h2o_mort_b(pft) * seedling_mdds + &
EDPftvarcon_inst%seedling_h2o_mort_c(pft)
! Cap h2o mortality rate at 1, quadratic can produce values >1 for large moisture
! deficit days, if rate exceeds 1 write to fates log

if (seedling_h2o_mort_rate > 1.0_r8) then
write(fates_log(), *) 'TRS seedling_h2o_mort_rate > 1', &

@rgknox rgknox Jul 28, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

logging messages that don't have terminations associated with them are most helpful in my experience while actively debugging. Another alternative is to use FatesWarn() which will count the times this situation triggers and summarize at the end of a run. Another alternative would be create a history variable that tracks the mort rate in this situation.

@r-ward r-ward Aug 5, 2026

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Makes sense, will remove. I'm working on a separate PR to add history variables for some TRS-specific output; I'll add history variables or FatesWarn() messages there.

'pft=', pft, 'uncapped=', seedling_h2o_mort_rate, 'mdds=', seedling_mdds
end if
seedling_h2o_mort_rate = min(1.0_r8, seedling_h2o_mort_rate)
end if ! mdd threshold check

! Step 3. Sum modes of mortality (including background mortality) and send dead seedlings
! to litter
litt%seed_germ_decay(pft) = (litt%seed_germ(pft) * seedling_light_mort_rate) + &
(litt%seed_germ(pft) * seedling_h2o_mort_rate) + &
(litt%seed_germ(pft) * EDPftvarcon_inst%background_seedling_mort(pft) &
* years_per_day)

! Cap total mortality rate at 1 to prevent negative seedling pool
total_seedling_mort_rate = seedling_light_mort_rate + &
seedling_h2o_mort_rate + &
(EDPftvarcon_inst%background_seedling_mort(pft) * years_per_day)
if (total_seedling_mort_rate > 1.0_r8) then
write(fates_log(), *) 'TRS total seedling mortality rate > 1', &

@rgknox rgknox Jul 28, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same comment as above

'pft=', pft, 'uncapped=', total_seedling_mort_rate, &
'seedling_light_mort_rate=', seedling_light_mort_rate, &
'seedling_h2o_mort_rate=', seedling_h2o_mort_rate
end if
litt%seed_germ_decay(pft) = litt%seed_germ(pft) * &
min(1.0_r8, total_seedling_mort_rate)

else

litt%seed_germ_decay(pft) = litt%seed_germ(pft) * &
Expand Down Expand Up @@ -2440,6 +2464,15 @@ subroutine SeedGermination( litt, cold_stat, drought_stat, currentPatch )
if ( seedling_layer_smp .GE. EDPftvarcon_inst%seedling_psi_emerg(pft) ) then
seedling_emerg_rate = photoblastic_germ_modifier * EDPftvarcon_inst%a_emerg(pft) * &
wetness_index**EDPftvarcon_inst%b_emerg(pft)
! Cap emergence rate at 1, rate can exceed 1 for some parameter combinations,
! leading to negative seed bank, write to log if this is the case
if (seedling_emerg_rate > 1.0_r8) then
write(fates_log(),*) 'TRS seedling_emerg_rate > 1', &
'pft=', pft, 'uncapped=', seedling_emerg_rate, &
'photoblastic_germ_modifier=', photoblastic_germ_modifier, &
'wetness_index=', wetness_index
end if
seedling_emerg_rate = min(1.0_r8, seedling_emerg_rate )
else

seedling_emerg_rate = 0.0_r8
Expand Down Expand Up @@ -2522,6 +2555,7 @@ subroutine recruitment(currentSite, currentPatch, bc_in)
real(r8) :: stem_drop_fraction !
real(r8) :: fnrt_drop_fraction !
real(r8) :: sdlng2sap_par ! running mean of PAR at the seedling layer [MJ/m2/day]
real(r8) :: sdlng2sap_rate ! seedling-to-sapling transition rate factor (unitless)
real(r8) :: seedling_layer_smp ! soil matric potential at seedling rooting depth [mm H2O suction]
integer, parameter :: recruitstatus = 1 ! whether the newly created cohorts are recruited or initialized
integer :: ilayer_seedling_root ! the soil layer at seedling rooting depth
Expand Down Expand Up @@ -2667,11 +2701,22 @@ subroutine recruitment(currentSite, currentPatch, bc_in)
sdlng2sap_par = currentPatch%sdlng2sap_par%GetMean()* &
sec_per_day*megajoules_per_joule

sdlng2sap_rate = EDPftvarcon_inst%seedling_light_rec_a(ft)* &
sdlng2sap_par**EDPftvarcon_inst%seedling_light_rec_b(ft)

! If the seedling to sapling transition rate exceeds 1,
! cap at 1 (prevent seed_germ from going negative) and write to fates log

if (sdlng2sap_rate > 1.0_r8) then
write(fates_log(), *) 'TRS seedling-to-sapling transition rate > 1', &
'pft=', ft, 'uncapped=', sdlng2sap_rate, &
'sdlng2sap_par=', sdlng2sap_par
end if

mass_avail = currentPatch%area* &
currentPatch%litter(el)%seed_germ(ft)* &
EDPftvarcon_inst%seedling_light_rec_a(ft)* &
sdlng2sap_par**EDPftvarcon_inst%seedling_light_rec_b(ft)

min(1.0_r8, sdlng2sap_rate)

! If soil moisture is below pft-specific seedling moisture stress threshold the
! recruitment does not occur.
ilayer_seedling_root = minloc(abs(bc_in%z_sisl(:) - &
Expand Down Expand Up @@ -2780,6 +2825,7 @@ subroutine recruitment(currentSite, currentPatch, bc_in)
currentPatch%litter(el)%seed_germ(ft) = &
currentPatch%litter(el)%seed_germ(ft) - cohort_n / currentPatch%area * &
(m_struct + m_leaf + m_fnrt + m_sapw + m_store + m_repro)

end if

end do
Expand Down
2 changes: 1 addition & 1 deletion biogeochem/FatesPatchMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -666,7 +666,7 @@ subroutine InitRunningMeans(this, current_tod, regeneration_model, numpft)
call this%seedling_layer_par24%InitRMean(fixed_24hr, &
init_value=init_seedling_par, init_offset=real(current_tod, r8))
call this%sdlng_mort_par%InitRMean(ema_sdlng_mort_par, &
init_value=temp_init_veg)
init_value=init_seedling_par)
call this%sdlng2sap_par%InitRMean(ema_sdlng2sap_par, &
init_value=init_seedling_par)

Expand Down
12 changes: 8 additions & 4 deletions main/FatesInterfaceMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2335,8 +2335,12 @@ subroutine UpdateFatesRMeansTStep(sites,bc_in, bc_out)
endif

! Update the seedling layer smp and mdd running means
call cpatch%sdlng_emerg_smp(pft)%p%UpdateRMean(new_seedling_layer_smp)
call cpatch%sdlng_mdd(pft)%p%UpdateRMean(new_seedling_mdd)
! Only update after first model day to avoid adding
! initialization spike in smp to running means
if (hlm_model_day > 2.0_r8) then
call cpatch%sdlng_mdd(pft)%p%UpdateRMean(new_seedling_mdd)
call cpatch%sdlng_emerg_smp(pft)%p%UpdateRMean(new_seedling_layer_smp)

@rgknox rgknox Jul 28, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We want to make sure that hlm_model_day is continuous across restarts. If for some reason someone is running two day restarts, and hlm_model_day is resetting each time, this running mean will never increment with meaningful data.

While I don't think SMP spikes are unexpected at the beginning of a simulation, and your move to ignore them makes sense. But I also wanted to check in regarding frozen soils. I think there was some discussion of it here:

#1500

SHould we apply some filters to prevent calculations when soils are frozen?

@r-ward r-ward Aug 5, 2026

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, I see. Re: frozen soils - I think you're right that we need some conditionality, but I don't think we should stop calculations when soils are frozen. We could do something analogous to the phenology code (https://github.com/r-ward/fates/blob/3b58472a513e39fcb8d0dafc7bc0ae31b91f6e95/biogeochem/EDPhysiologyMod.F90#L1222)) and use a lower bound value similar to smp_lwr_bound for cases where check_layer_water finds no liquid water. Frozen conditions are physiologically dry, so I think that makes sense (although really we should implement behavior dependent on moisture and chilling degree days for temperate/boreal forest uses of TRS).

In the tropics, though, using an smp_lwr_bound during the SMP initialization spike suppresses germination and inflates seedling mortality for weeks/months, so this fix alone doesn't work.

Can we differentiate real frozen soils from initialization/erroneous freezing by adding a check on soil temperature? If soil temperature is warm and there is no liquid water, skip the calculations. If soil temperature is freezing, and there is no liquid water, use an smp_lwr_bound. @lmkueppers - what do you think?

endif

enddo !end pft loop

Expand Down Expand Up @@ -2410,8 +2414,8 @@ subroutine SeedlingParPatch(cpatch, &

! Start with the assumption that there is a single canopy layer
seedling_par_high = atm_par_dir+atm_par_dif
par_high_frac = 1._r8-cpatch%total_canopy_area
par_low_frac = cpatch%total_canopy_area
par_high_frac = 1._r8 - (cpatch%total_canopy_area / cpatch%area)
par_low_frac = cpatch%total_canopy_area / cpatch%area

! Work up through the canopy layers from the bottom layer
do cl = cpatch%NCL_p,max(1,cpatch%NCL_p-1),-1
Expand Down