Bug in SFMainMod? Don't ever remove trunks from `sum_fuel`?

[adrifoster]

In refactoring the SPITFIRE code, it seems like we don't actually ever remove trunk litter from patch%sum_fuel. We do remove trunk characteristics from the weighted averages of bulk density, SAV, etc., but from my read of the code we do not ever correct patch%sum_fuel from it's original calculation (here):

       currentPatch%sum_fuel =  sum(litt_c%leaf_fines(:)) + &
                                sum(litt_c%ag_cwd(:)) + &
                                currentPatch%livegrass

Note that litt_c%ag_cwd(:) is defined here:

   integer, public, parameter :: ncwd  = 4    ! number of coarse woody debris pools 
                                              ! (twig,s branch,l branch, trunk)
   real(r8)                   :: ag_cwd(ncwd) ! above ground coarse wood debris (cwd) [kg/m2]

It is then used in the rate of spread subroutine (here):

ir = reaction_v_opt*(currentPatch%sum_fuel/0.45_r8)*SF_val_fuel_energy*moist_damp*SF_val_miner_damp 

Which, according to the original Rothermel equations and the SPITFIRE paper , should not include the 1000-hr (i.e. trunk) fuels.

If my read is correct, updating this will certainly change answers...

Is this an update that happened? or was it always like this?

tagging relevant folks... @ckoven @rgknox @glemieux @samsrabin

[adrifoster]

I'll note that the patch-level variable sum_fuel (here) is noted as:

real(r8)              :: sum_fuel                ! total ground fuel related to ROS (omits 1000 hr fuels) [kgC/m2]

So at some point we either intended to or actually did omit the trunks from this variable

[rgknox]

@adrifoster , my interpretation is that sum fuel should interpret was is actually there, and evaluate all litter that can be considered fuel. If that portion of the litter (trunk wood) is not there or has recently been changed, then it should be reflected in litt_c%ag_cwd as that is the prognostic variable. If there is no change in the amount of litter in existence, and this is more about what litter is included in the sum, then we should be changing how we calculate sum_fuel. Is this kind of what you are getting at?

[adrifoster]

If there is no change in the amount of litter in existence, and this is more about what litter is included in the sum, then we should be changing how we calculate sum_fuel. Is this kind of what you are getting at?

Yes, so the 1000-hr fuels (i.e. trunks) are not supposed to be used to calculate fire rate of spread, intensity, etc. So the sum_fuel should either be modified to remove trunks if we are going to use it in those equations, or we need to use a different variable.

[samsrabin]

I think you're right, Adrianna; I don't see anywhere that sum_fuel gets modified.

In addition to the comment at the definition of sum_fuel you mentioned, there's this bit further indicating an intention to exclude trunks from the ROS calculation:

! Correct averaging for the fact that we are not using the trunks pool for fire ROS and intensity (5)
! Consumption of fuel in trunk pool does not influence fire ROS or intensity (Pyne 1996)
currentPatch%fuel_bulkd     = currentPatch%fuel_bulkd     * (1.0_r8/(1.0_r8-currentPatch%fuel_frac(tr_sf)))
currentPatch%fuel_sav       = currentPatch%fuel_sav       * (1.0_r8/(1.0_r8-currentPatch%fuel_frac(tr_sf)))
currentPatch%fuel_mef       = currentPatch%fuel_mef       * (1.0_r8/(1.0_r8-currentPatch%fuel_frac(tr_sf)))
currentPatch%fuel_eff_moist = currentPatch%fuel_eff_moist * (1.0_r8/(1.0_r8-currentPatch%fuel_frac(tr_sf)))

[adrifoster]

Okay thanks @samsrabin. So I guess we need to update this... Should this come in as a hot-fix? I'd like to change it first so I can ensure my refactors aren't doing something wonky.

[rgknox]

I'm happy to help fast track this and get it in before we update to API 35. I'd recommend getting the fix based on main though, and not wrapping it into the fire refactor

[adrifoster]

I'm happy to help fast track this and get it in before we update to API 35. I'd recommend getting the fix based on main though, and not wrapping it into the fire refactor

thanks! working on sending a quick fix PR now.

[adamhb]

Oh wow, thanks @adrifoster for finding this. I've had a lot of issues with fire intensity being too high for my simulations in Californian mixed conifer forests, and I guess this would explain why. So if I understand correctly what you found, the model would essentially treat logs on the landscape as if they had the properties of finer fuels (SAV etc.) and then include all that extra mass in the calculations of intensity / spread?

[samsrabin]

@adamhb From what I can tell, this should only influence rate of spread. The code I reference here should mean that trunks aren't included in fuel_bulkd, fuel_sav, fuel_mef, and fuel_eff_moist.

[ckoven]

@adamhb From what I can tell, this should only influence rate of spread. The code I reference here should mean that trunks aren't included in fuel_bulkd, fuel_sav, fuel_mef, and fuel_eff_moist.

@samsrabin but if this affects rate of spread, then it should also affect intensity and area burned too, since those are both dependent on rate of spread, correct?

[samsrabin]

Oh yeah, it would definitely affect burned area. I didn't think that ROS entered into the intensity (kW/m) calculations, but if that's the case then yes that will be affected too.

[samsrabin]

And yes, you're right: https://github.com/NGEET/fates/blob/b054cd00ee3740f828b57b57768f8ec0bd1bc874/fire/SFMainMod.F90#L857-L859

So yes, @adamhb, this will affect intensity!

[adamhb]

Ok thanks! Most simulations Ive run have very large 1000 fuel loads, so this will def change results. Especially because the 1000 hr fuel isnt actually consumed...

[adrifoster]

Hi @adamhb yes, this bug fix should decrease rate of spread, intensity, mortality, etc. like @samsrabin said.

But also yeah as @samsrabin said, it's only erroneously included in the sum_loading variable part of ROS, NOT in the calculations for SAV, bulk density, etc.

I'm still running my "very long" (i.e. 10 years) simulations to test potential differences so I will post those here, but please feel free to add any if you have a simulation you want to kick off.

[adrifoster]

Ok thanks! Most simulations Ive run have very large 1000 fuel loads, so this will def change results. Especially because the 1000 hr fuel isnt actually consumed...

@adamhb is that true? I think it should be consumed... just not used to calculate ROS or FI...

[adamhb]

@adrifoster these lines of code suggest that trunk fuel is not consumed right? !---calculate overall fuel consumed by spreading fire --- ! ignore 1000hr fuels. Just interested in fuels affecting ROS currentPatch%TFC_ROS = sum(FC_ground)-FC_ground(tr_sf)

"currentPatch%TFC_ROS" is defined here as "total fuel consumed by flaming front (kgC/m2 of burned area)"

I think the theory behind this is that logs generally are consumed by smoldering, which is not really a capability of spitfire, which is more geared towards the dynamics of ROS and flaming front.

[adamhb]

@adrifoster However, I think it would be great to having smoldering as a capability to avoid unrealistic build up of trunk fuel.

[adrifoster]

@adamhb true, it's not included in TFC_ROS but here it's being consumed:

       do c = 1,ncwd

          ! Transfer above ground CWD
          donatable_mass     = curr_litt%ag_cwd(c) * patch_site_areadis * &
                               (1._r8 - currentPatch%burnt_frac_litter(c))

          burned_mass        = curr_litt%ag_cwd(c) * patch_site_areadis * &
                               currentPatch%burnt_frac_litter(c)

          new_litt%ag_cwd(c) = new_litt%ag_cwd(c) + donatable_mass*donate_m2
          curr_litt%ag_cwd(c) = curr_litt%ag_cwd(c) + donatable_mass*retain_m2

          site_mass%burn_flux_to_atm = site_mass%burn_flux_to_atm + burned_mass

          ! Transfer below ground CWD (none burns)

          do sl = 1,currentSite%nlevsoil
             donatable_mass         = curr_litt%bg_cwd(c,sl) * patch_site_areadis
             new_litt%bg_cwd(c,sl)  = new_litt%bg_cwd(c,sl) + donatable_mass*donate_m2
             curr_litt%bg_cwd(c,sl) = curr_litt%bg_cwd(c,sl) + donatable_mass*retain_m2
          end do

[adamhb]

@adrifoster ok yes I see. I suspect that the brunt_frac_litter value for trunk fuel is typically near zero as calculated here (but that's a separate issue and I haven't confirmed in).

[adrifoster]

@adamhb okay good to know. Happy to discuss with you if you want to go over it, otherwise I can be sure to add some specific unit testing for it during my refactoring work.

Incidentally, I actually updated this set of equations (for a different model) based on some observational data I had (specific to boreal region...) but we could think about investigating

        real, dimension(FL_LEVS), parameter :: PARA = [-1.36, -1.36, -1.32,    &
            -0.83, -0.17, 0.06, -1.72, -1.72, 0.06, -1.72]
        real, dimension(FL_LEVS), parameter :: PARB = [0.54, 0.54, 0.41,       &
            -0.22, -0.82, -0.87, 0.57, 0.57, -0.87, 0.57]
        real, dimension(FL_LEVS), parameter :: PARC = [0.95, 0.95, 0.96, 1.02, &
            0.98, 0.55, 0.98, 0.98, 0.82, 0.98]

       do ic = 1, FL_LEVS
            if (soil%fuels(ic) > epsilon(1.0)) then

                lmoist = soil%litter_moist(ic)/soil%litter_mef(ic)

                if (lmoist <= MIN_MOISTURE(ic)) then
                    ! Dry litter
                    soil%frac_burnt(ic) = 1.0
                    burned(ic) = soil%fuels(ic)

                else if (lmoist > MIN_MOISTURE(ic) .and. lmoist <= 1.0) then

                    soil%frac_burnt(ic) = max(0.0, min(1.0, PARA(ic)*lmoist**2 +   &
                        PARB(ic)*lmoist + PARC(ic)))
                    burned(ic) = soil%frac_burnt(ic)*soil%fuels(ic)

                else if (lmoist > 1.0) then
                    ! High moisture
                    soil%frac_burnt(ic) = 0.0
                    burned(ic) = 0.0
                end if
            end if
        end do

[adamhb]

@adrifoster Thanks for sharing! I think a unit test would be a good idea. Just looking back at some old output, the trunk fuel moisture didn't really change going into the dry season (stayed very high), so not much of it burned. But, this might be what should happen in a model that focuses on the dynamics at the flaming front...

[glemieux]

@adrifoster can we close this per #1180? Or is this issue broader than just that fix?

[adrifoster]

@glemieux no we can close

[ckoven]

just noticed that this was fixed by #1180, but still open. closing.