Including C-based harvest in FATES

[sshu88]

Motivation: A common approach in land models to estimate the C emission from wood logging and corresponding LUC activities is achieved by reading the forest cover loss, i.e., the harvested area data. However, more accurate data based on national statistics is in the unit of harvested biomass/carbon amount. For example, in LUH2, quantities of harvested wood in biomass units are converted into harvested areas based on the biomass density estimation from the MIAMI-LU ecosystem model. Thus, C-based harvest, which reads in harvest rate in the unit of harvested biomass/carbon amount, can guarantee a much more precise estimation of the direct LUC emission from logging.
Another reason is that the FATES model can simulate a more realistic scenario of C-based harvest. Since FATES can determine the dominant plant community at each site, area-based harvest only transfers a certain amount of area into secondary forest without considering the actual loss of forest biomass amount, which sometimes can cause “blank harvest” if no woody biomass available for logging in certain sites.

Note: Introducing the ability for harvesting based on mass (carbon) will not replace the area-based harvesting, but just be included as an additional option.

Design:
First introduce some definitions here to avoid inconsistency:
Harvestable carbon: All available woody carbon for harvest in a site, which excluded the carbon from immature woody PFTs (forest and shrub), defined by all the plants with a dbh lower than a certain threshold (dbhmin).
Available carbon for harvest: The total amount of woody carbon available for harvest.
The basic logic of the C-based harvest is to convert the carbon-based harvest rate to the area fraction-based harvest rate, shown in the following flow chart:

Two subroutines are added in EDLoggingMortalityMod.F90:

1. get_harvestable_carbon(): This subroutine calculates the amount of available harvestable carbon of three harvest categories (primary, secondary mature, secondary young) for each site. These values are required by the following subroutine get_harvest_rate_carbon(), and need to be calculated before entering the cohort loop.

Below is the current workflow and description of what happens in: get_harvestable_patch_carbon(), but comments and suggestions are welcomed.

1. Enter the patch loop and check which category (primary, secondary mature and secondary young) the current patch belongs to.
2. Enter the cohort loop and examine if current cohort matches the requirement for wood harvest (i.e., pft = woody and dbh > dbhmin, which is 50cm by default and can be modified in FATES parameter file)
3. Once the criteria is met, retrieve the aboveground carbon of sapwood and structural wood of the current cohort, and also consider the loss of carbon into litter when harvest to calculate the harvestable carbon.
4. If the current cohort does not match the criteria, do the same calculation but add to the available carbon for harvest instead of harvestable carbon.
5. Aggregate the carbon from all cohorts in the current patch
6. Once finishing the patch loop, output the harvestable carbon of all three different categories under both circumstances (harvestable carbon and available carbon for harvest) in the unit of KgC.

2. get_harvest_rate_carbon(): first read in the C-based harvest rate from data, then determine if harvestable carbon can fulfill the harvest rate or not and then transfer actual harvested carbon into area fraction for each cohort. This subroutine is called within a cohort loop, thus any new cohort and patch level loops are not allowed in this subroutine, thus a new subroutine outside the loop is required.

Below is the current workflow and description of what happens in: get_harvest_rate_carbon()

1. Get the harvest rate (KgC per site) of three categories from the ELM-FATES interface. We call this harvest rate the harvest demand in the following text.
2. Check which category the current cohort belongs to and divide the corresponding harvest rate by the harvestable carbon amount and obtain the area fraction.
3. Check if biomass demand from the logging event can be fulfilled by harvestable carbon. If so (the fraction is less than 1), the calculated area fraction will be the final output.
4. If biomass demand cannot be fulfilled by the total harvestable amount, but can be fulfilled by potential carbon, we either release the restriction for harvestable carbon (dbhmin) and account for all potential carbon when the added tag “bypass_harvest_criteria” is set to true, or assign zero to the harvest rate if the tag is set to false. The harvest rate (in area) equals the harvest rate in carbon divided by the unharvestable carbon.
5. If biomass demand cannot be fulfilled by both harvestable and unharvestable potential carbon, we assign zero to the harvest rate and put all of the harvest demand as “wood harvest debt”.
6. Output the harvest rate (in area fraction) for the current cohort.

Subroutine get_harvestable_carbon() is called at the beginning of calculating the disturbance rate. Subroutine get_harvest_rate_carbon() is an alternative for the similar subroutine get_harvest_rate_area() and will be called to replace get_harvest_rate_area() in FATES when reading C-based harvest rate.

Site level evaluation:

1. Brazil site level evaluation:
   Performed the same site level simulation (100 year spinup and 50 years simulation 1951 - 2000).

1. Checked the consistency between the harvested carbon from the dataset and the actual harvested carbon calculated.
   

2. Compare the calculated harvested area fraction under the C-based harvest method to the area-based harvest rate from LUH2.
   

Seems like the calculated secondary forest fraction is less than the harvested area fraction by LUH2, which is due to a higher biomass on the current site compared to the biomass data used by LUH2.