Seed pool can become negative when using TRS with seedling dynamics, eventually crashing fates

### Describe the issue


When running FATES with the Tree Recruitment Scheme (TRS) and seedling dynamics, the seed pool can become negative, eventually crashing the model. I think this is due to two problems in the calculation of `seed_germ_in` (the mass germinating out of the seed pool). TRS then uses `seed_germ` for its construction costs. 

**Background:** The seed pool is updated in the PreDisturbanceIntegrateLitter Mod: 
https://github.com/NGEET/fates/blob/75a465c0b650fad9e40059204903f6017f22a1ee/biogeochem/EDPhysiologyMod.F90#L544-L552

When TRS with seedling dynamics is enabled, `seed_decay` includes non-reproductive biomass from in-fluxing litter and decay from the existing seed bank.  

And `seed_germ_in` = an environmentally sensitive fraction of the seed pool: `seed` * `seedling_emerg_rate` (Step 4 below):

https://github.com/NGEET/fates/blob/75a465c0b650fad9e40059204903f6017f22a1ee/biogeochem/EDPhysiologyMod.F90#L2402-L2445

**There are two problems here:** 
**Problem 1.** Order of operations: when `seedling_emerg_rate` is high, the updated seed pool can become negative if (seed pool * seed_emerg_rate)  > (seed pool - seed_decay)   
- Potential fix: calculate in two steps to make sure that the material in seed_germ_in does not exceed the available seed material after decay e.g.
```
! Calculate 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
   litt%seed_germ_in(pft) = min(litt%seed_germ_in(pft), net_seed_available))
   
   ! Update pools
   litt%seed(pft) = net_seed_available - litt%seed_germ_in(pft)
   litt%seed_germ(pft) = litt%seed_germ(pft) + litt%seed_germ_in(pft) - litt%seed_germ_decay(pft)
```
   
**Problem 2.** The `seedling_emerg_rate` is currently not capped, and can be > 1,  which contributes to problem 1 above.   
- Potential fix: cap seedling emergence at 1?

As an example, I'm pasting log output from a day 353 of single-point run at BCI with 2-pfts.
PFT1 = light demanding (LD) PFT 2 = shade tolerant (ST), I used TRS parameter values from the tech documentation when LD/ST specific values are given (Table 1.7 , https://fates-users-guide.readthedocs.io/projects/tech-doc/en/latest/fates_tech_note.html). 

The eq for the seedling_emerg_rate is : 
```
 seedling_emerg_rate = photoblastic_germ_modifier * EDPftvarcon_inst%a_emerg(pft) * &
        wetness_index**EDPftvarcon_inst%b_emerg(pft)
```
In this case the parameter of interest is b_emerg (a_emerg is 0.0003 for both PFTs, from log output: photoblastic_germ_modifier = 0.93753000329937553, wetness index = 390.85645048765576).      

When b_emerg is 1.2 (ST pft) everything works fine: seedling emergence rate = 0.36268297569017249     
When b_emerg is 1.6 (LD pft) , seedling emergence rate = 3.9476085760981108     

The model crashes at day 646 with `Patch area not closing` (error message pasted below).

I'm surprised the parameters tested for the LD PFT are causing >1 emergence rates, and I’m not sure if the formulation itself needs tweaking, or if capping seedling_emergence at 1 is sufficient.

Advice would be appreciated! 

### Relevant log output

```shell
===SEED DECAY DIADNOSTICS====
 Model Day =    353.02083333333331     
 PFT:           1
 Seedling par (raw):   83.540806467338086     
 Seedling par (converted):   230.97362172089632     
 Seedling mdds:   0.0000000000000000     
 Seedling mdd crit:   1400000.0000000000     
 Seedling light mort rate:   1.7994349006037014E-101
 Seedling h2o mort rate:   0.0000000000000000     
 Seed_germ pool (before decay):   2.6708814696409877E-006
 Seed_germ_decay (calculated):   1.2805596087319802E-009
 Background_seedling_mort:   4.7945205479452054E-004
 Total mort rate:   4.7945205479452054E-004
 Does seed decay exceeds pool?: F
 ===SEED DECAY DIADNOSTICS====
 Model Day =    353.02083333333331     
 PFT:           2
 Seedling par (raw):   83.540806467338086     
 Seedling par (converted):   230.97362172089632     
 Seedling mdds:   0.0000000000000000     
 Seedling mdd crit:   1400000.0000000000     
 Seedling light mort rate:   1.7994349006037014E-101
 Seedling h2o mort rate:   0.0000000000000000     
 Seed_germ pool (before decay):   2.5402458337975548E-006
 Seed_germ_decay (calculated):   1.0439366440263924E-009
 Background_seedling_mort:   4.1095890410958901E-004
 Total mort rate:   4.1095890410958901E-004
 Does seed decay exceeds pool?: F
 === SeedGermination DIAGNOSTICS ===
 Model Day:   353.02083333333331     
 PFT:           1
 seed pool:  -1.3374934093787974E-006
 seedling_layer_par:   9.8450410540560043     
 photoblastic_germ_modifier:  0.93753000329937553     
 seedling_layer_smp:  -261.06980247919256     
 wetness_index:   390.85645048765576     
 seedling_emerg_rate:   3.9476085760981108     
 seed_germ_in (calculated):  -5.2799004533384421E-006
 === SeedGermination DIAGNOSTICS ===
 Model Day:   353.02083333333331     
 PFT:           2
 seed pool:   3.5904396208826199E-007
 seedling_layer_par:   9.8450410540560043     
 photoblastic_germ_modifier:  0.93753000329937553     
 seedling_layer_smp:  -261.06980247919256     
 wetness_index:   390.85645048765576     
 seedling_emerg_rate:  0.36268297569017249     
 seed_germ_in (calculated):   1.3021913257376033E-007
 === PreDisturbanceIntegrateLitter DIAGNOSTICS ===
 Model Day:   353.02083333333331     
 PFT:           1
 seed_germ (before update):   2.6708814696409877E-006
 seed_germ_in:  -5.2799004533384421E-006
 seed_germ_decay:   1.2805596087319802E-009
 seed_germ (after update):  -2.6102995433061863E-006
 Net change:  -5.2811810129471739E-006
 Decay exceeds pool?: T
 === PreDisturbanceIntegrateLitter DIAGNOSTICS ===
 Model Day:   353.02083333333331     
 PFT:           2
 seed_germ (before update):   2.5402458337975548E-006
 seed_germ_in:   1.3021913257376033E-007
 seed_germ_decay:   1.0439366440263924E-009
 seed_germ (after update):   2.6694210297272885E-006
 Net change:   1.2917519592973395E-007
 Decay exceeds pool?: F




 PATCH AREA CHECK NOT CLOSING
 patch area:   7.3448618706999813E-002
 fraction that is imperfect (unclosed):   0.0000000000000000     
 lat:   9.1530000000000005     
 lon:   280.15390000000002     
 spread:  0.84999999999999987     
 -----------------------------------------
 layer:            1  area:    7.3448618706999980E-002
 bias [m2] (layer-patch):    1.6653345369377348E-016
 -----------
  co area:   9.0581374081781301E-010
  co dbh:   0.76451299099786252     
  co pft:            2
  co n:    4.9653109197993802E-009
 -----------
  co area:   8.9155179987784708E-010
  co dbh:   0.76407217793207960     
  co pft:            2
  co n:    4.8913116373993979E-009
 -----------
  co area:   7.3448616909634445E-002
  co dbh:   0.76407217793207960     
  co pft:            1
  co n:   0.40296040531824168     
 -----------------------------------------
 layer:            2  area:    7.3448618706999536E-002
 bias [m2] (layer-patch):   -2.7755575615628914E-016
 -----------
  co area:   9.0581374081780639E-010
  co dbh:   0.76451299099786252     
  co pft:            2
  co n:    4.9653109197993438E-009
 -----------
  co area:   8.9155179987785338E-010
  co dbh:   0.76407217793207960     
  co pft:            2
  co n:    4.8913116373994327E-009
 -----------
  co area:   7.3448616909634001E-002
  co dbh:   0.76407217793207960     
  co pft:            1
  co n:   0.40296040531823929     
 -----------------------------------------
 layer:            3  area:    7.3448618706999536E-002
 bias [m2] (layer-patch):   -2.7755575615628914E-016
 -----------
  co area:   9.0581374081780660E-010
  co dbh:   0.76451299099786252     
  co pft:            2
  co n:    4.9653109197993447E-009
 -----------
  co area:   8.9155179987785039E-010
  co dbh:   0.76407217793207960     
  co pft:            2
  co n:    4.8913116373994161E-009
 -----------
  co area:   7.3448616909634001E-002
  co dbh:   0.76407217793207960     
  co pft:            1
  co n:   0.40296040531823923     
 -----------------------------------------
 layer:            4  area:    7.3448618707000438E-002
 bias [m2] (layer-patch):    6.2450045135165055E-016
 -----------
  co area:   9.0581374081781301E-010
  co dbh:   0.76451299099786252     
  co pft:            2
  co n:    4.9653109197993802E-009
 -----------
  co area:   8.9155179987785359E-010
  co dbh:   0.76407217793207960     
  co pft:            2
  co n:    4.8913116373994335E-009
 -----------
  co area:   7.3448616909634903E-002
  co dbh:   0.76407217793207960     
  co pft:            1
  co n:   0.40296040531824417     
 -----------------------------------------
 layer:            5  area:    7.3448618706999536E-002
 bias [m2] (layer-patch):   -2.7755575615628914E-016
 -----------
  co area:   9.0581374081780329E-010
  co dbh:   0.76451299099786252     
  co pft:            2
  co n:    4.9653109197993265E-009
 -----------
  co area:   8.9155179987784356E-010
  co dbh:   0.76407217793207960     
  co pft:            2
  co n:    4.8913116373993789E-009
 -----------
  co area:   7.3448616909634001E-002
  co dbh:   0.76407217793207960     
  co pft:            1
  co n:   0.40296040531823923     
 -----------------------------------------
 layer:            6  area:    7.3448618706999536E-002
 bias [m2] (layer-patch):   -2.7755575615628914E-016
 -----------
  co area:   9.0581374081780639E-010
  co dbh:   0.76451299099786252     
  co pft:            2
  co n:    4.9653109197993438E-009
 -----------
  co area:   8.9155179987785028E-010
  co dbh:   0.76407217793207960     
  co pft:            2
  co n:    4.8913116373994153E-009
 -----------
  co area:   7.3448616909634001E-002
  co dbh:   0.76407217793207960     
  co pft:            1
  co n:   0.40296040531823923     
 -----------------------------------------
 layer:            7  area:    7.3448618706999980E-002
 bias [m2] (layer-patch):    1.6653345369377348E-016
 -----------
  co area:   9.0581374081780970E-010
  co dbh:   0.76451299099786252     
  co pft:            2
  co n:    4.9653109197993620E-009
 -----------
  co area:   8.9155179987784708E-010
  co dbh:   0.76407217793207960     
  co pft:            2
  co n:    4.8913116373993979E-009
 -----------
  co area:   7.3448616909634445E-002
  co dbh:   0.76407217793207960     
  co pft:            1
  co n:   0.40296040531824168     
 -----------------------------------------
 layer:            8  area:    7.3448618706999980E-002
 bias [m2] (layer-patch):    1.6653345369377348E-016
 -----------
  co area:   9.0581374081780970E-010
  co dbh:   0.76451299099786252     
  co pft:            2
  co n:    4.9653109197993620E-009
 -----------
  co area:   8.9155179987784697E-010
  co dbh:   0.76407217793207960     
  co pft:            2
  co n:    4.8913116373993971E-009
 -----------
  co area:   7.3448616909634445E-002
  co dbh:   0.76407217793207960     
  co pft:            1
  co n:   0.40296040531824173     
 -----------------------------------------
 layer:            9  area:                        NaN
 bias [m2] (layer-patch):                        NaN
 -----------
  co area:   9.0581374081781322E-010
  co dbh:   0.76451299099786252     
  co pft:            2
  co n:    4.9653109197993810E-009
 -----------
  co area:   8.9155179987785369E-010
  co dbh:   0.76407217793207960     
  co pft:            2
  co n:    4.8913116373994343E-009
 -----------
  co area:   7.3448616909634903E-002
  co dbh:   0.76407217793207960     
  co pft:            1
  co n:   0.40296040531824417     
 -----------------------------------------
 layer:           10  area:    7.2269511633305839E-317
 bias [m2] (layer-patch):   -7.3448618706999813E-002
 -----------
  co area:   9.0581374081780639E-010
  co dbh:   0.76451299099786252     
  co pft:            2
  co n:    4.9653109197993438E-009
 -----------
  co area:   8.9155179987785039E-010
  co dbh:   0.76407217793207960     
  co pft:            2
  co n:    4.8913116373994161E-009
 -----------
  co area:   7.3448616909634223E-002
  co dbh:   0.76407217793207960     
  co pft:            1
  co n:   0.40296040531824051     
 -----------------------------------------
 layer:           11  area:    7.3448598384861494E-002
 bias [m2] (layer-patch):   -2.0322138291617442E-008
 -----------
  co area:   9.0581374081780639E-010
  co dbh:   0.76451299099786252     
  co pft:            2
  co n:    4.9653109197993438E-009
 -----------
  co area:   8.9155179987785039E-010
  co dbh:   0.76407217793207960     
  co pft:            2
  co n:    4.8913116373994161E-009
 -----------
  co area:   7.3448616909634223E-002
  co dbh:   0.76407217793207960     
  co pft:            1
  co n:   0.40296040531824051     
 -----------------------------------------
 layer:           12  area:    1.8257451319375437     
 bias [m2] (layer-patch):    1.7522965132305439     
 -----------
  co area:   2.2516216871790394E-008
  co dbh:   0.76451299099786252     
  co pft:            2
  co n:    1.2342495202726042E-007
 -----------
  co area:   2.2161701433631037E-008
  co dbh:   0.76407217793207960     
  co pft:            2
  co n:    1.2158551880187157E-007
 -----------
  co area:   1.8257450872596255     
  co dbh:   0.76407217793207960     
  co pft:            1
  co n:    10.016566837127510     
 ENDRUN: 
 ERROR in /glade/u/home/rward/CTSM/src/fates/biogeochem/EDCanopyStructureMod.F90 at line 327
```

### FATES tag

sci.1.87.5_api.41.0.0-16-g59183409

### Host land model tag

ctsm5.3.077

### Machine

derecho

### Other supported machine name

_No response_

### Additional context

_No response_

	! Update the bank of viable seeds
	! -----------------------------------------------------------------------------------

	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)

	! If TRS seedling dynamics is switched on we calculate seedling emergence (i.e. germination)
	! as a pft-specific function of understory light and soil moisture.
	else if ( hlm_regeneration_model == TRS_regeneration .and. &
	prt_params%allom_dbh_maxheight(pft) > min_max_dbh_for_trees ) then

	! Step 1. Calculate how germination rate is modified by understory light
	! This applies to photoblastic germinators (e.g. many tropical pioneers)

	! Calculate mean PAR at the seedling layer (MJ m-2 day-1) over the prior 24 hours
	seedling_layer_par = currentPatch%seedling_layer_par24%GetMean() * sec_per_day * megajoules_per_joule

	! Calculate the photoblastic germination rate modifier (Eq. 3 Hanbury-Brown et al., 2022)
	photoblastic_germ_modifier = seedling_layer_par / &
	(seedling_layer_par + EDPftvarcon_inst%par_crit_germ(pft))

	! Step 2. Calculate how germination rate is modified by soil moisture in the rooting zone of
	! the seedlings. This is a pft-specific running mean based on pft-specific seedling rooting
	! depth.

	! Get running mean of soil matric potential (mm of H2O suction) at the seedling rooting depth
	! This running mean based on pft-specific seedling rooting depth.
	seedling_layer_smp = currentPatch%sdlng_emerg_smp(pft)%p%GetMean()

	! Calculate a soil wetness index (1 / -soil matric pontential (MPa) ) used by the TRS
	! to calculate seedling mortality from moisture stress.
	wetness_index = 1.0_r8 / (seedling_layer_smp * (-1.0_r8) * mpa_per_mm_suction)

	! Step 3. Calculate the seedling emergence rate based on soil moisture and germination
	! rate modifier (Step 1). See Eq. 4 of Hanbury-Brown et al., 2022

	! If SMP is below a pft-specific value, then no germination occurs
	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)
	else

	seedling_emerg_rate = 0.0_r8

	end if ! End soil-moisture based seedling emergence rate

	! Step 4. Calculate the amount of carbon germinating out of the seed bank
	litt%seed_germ_in(pft) = litt%seed(pft) * seedling_emerg_rate

	end if if_tfs_or_def

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Seed pool can become negative when using TRS with seedling dynamics, eventually crashing fates #1504

Describe the issue

Relevant log output

FATES tag

Host land model tag

Machine

Other supported machine name

Additional context

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Seed pool can become negative when using TRS with seedling dynamics, eventually crashing fates #1504

Description

Describe the issue

Relevant log output

FATES tag

Host land model tag

Machine

Other supported machine name

Additional context

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions