Prescribe vegetation structure for regional/global simulations with FATES(ST3 Mode)?

[huitang-earth]

I am trying to run FATES with static stand structure (ST3) over Scandinavia. We want to prescribe vegetation structure according to our inventory based forest dataset over this region.

I notice that prescribing vegetation structure for a large-scale region in FATES has not been implemented yet. It is metioned in issue #166 that "possibly ST3 mode could be used at larger scales by initializing from a restart file or large-scale vegetation structure dataset (if/when such a thing exists)". I am wondering if there have been anyone doing some tests on this direction?

To do this quick, I am planning to modify a restart file to prescribe the vegetation structure, although it might not be the most decent way to do this. I foresee that there might be more problems than expected, and would be glad to know your comments and suggestions. Hopefully, our attempt could contribute to the development of FATES on this direction.

[pollybuotte]

Just yesterday I ran a regional simulation over a particular forest type in California that uses initialized stand structure. This advice from Charlie allowed me to get this run going...

"The first file that fates reads which lists all the css/pss files to read also has lats on lons associated with each one, so it uses the nearest neighbor site for each gridcell. One thing though is there is a hard-coded maximum distance here:
https://github.com/NGEET/fates/blob/master/main/FatesInventoryInitMod.F90#L88
which you may need to change to a larger value so that it lets you use more distant sites…"

I'm currently running this region with uniform stand structure, so one entry in the first stand initialization file and a max distance set to 10 degrees. I hope this helps you get going!

[rgknox]

hi @huitang-earth ,

We have some documentation on this:
https://github.com/NGEET/fates/wiki/Model-Initialization-Modes#Guide_for_a_Cold_Start_Simulation_with_Inventory_Case_3

The first file it mentions is the control file, where you can list out numerous pss/css files to associate with a coordinate. It will perform nearest neighbor when assigning an inventory couplet to the lat-lon location. Note, that multiple inventory couplets has not been tested in any detail though.

[huitang-earth]

@pollybuotte @rgknox Thanks very much for your useful information.
I've been thinking of using the existing "inventory" module (i.e., FatesInventoryInitMod.F90) in FATES. But I notice that in "FatesInventoryInitMod.F90", it states that "This procedure, if called from a massive grid, is both probably inappropriate, and probably time consuming".

I am planning to run FATES with quite high resolution (5x5km) for a large region. The sum of the grid cells are over 80000, and each grid cell should have its own unique stand structure. So, I am worried that "FatesInventoryInitMod.F90" will not work for my case, as it will have to deal with huge number of pss/css files. Do you think so?

Since our dataset is in netcdf format, I am wondering if making the module to read netcdf file instead of "pss/css" files will work better for the simulations with large number of grid cells?

[rgknox]

@huitang-earth , I think an alternative initialization mode that reads stand-structure data in a self-describing format like netcdf would be awesome and fully support this endeavor.

The only tricky part is that FATES does not haven any internal support for netcdf, as we rely on the ESMs to read this type of info, and then pass in via a boundary condition.

I personally feel good about supporting netcdf in FATES, as its pretty ubiquitous. But if this type of code goes back into NGEET/master, it would be good to get some consensus with the different teams using it. E3SM, NorESM, CESM all require netcdf support anyway, but maybe we should just check in with ATS to see how they feel about it. @xuchongang @ecoon @dasvyat

[huitang-earth]

@rgknox Thanks very much for the useful insights.
It seems to me that relying on ESMs (e.g., CLM5) to read the inventory data (in netcdf format) and pass them to FATES might be a good solution, and may be worth trying?

[ckoven]

hi @huitang-earth for a near-term solution, it may be easier for you to either (a) try using the text file method to read in your census data one time, and then write the resulting forest patch/cohort structure to a restart file, or (b) directly modify the patch/cohort structure in a restart file to correspond to your census data.

[rosiealice]

@huitang-earth I think I'd second @ckoven's suggestion b) in that you could write an arbitrary restart file and then overwrite it with the desired canopy structure. Maybe we should think about some jupyter notebook-type common tool for modifying restart files?

More broadly, we should discuss/think about what features of the model we would like to preserve in an ideal global coldstart. One could imagine, say, trying to restart with a particular PFT distribution and LAI, and then making some sort of educated guess r.e. the canopy structure underlying that LAI?

Are there any existing data products that could help us with this (e.g.
https://www.ilamb.org/CMIP6/historical/EcosystemandCarbonCycle/Biomass/GEOCARBON/GEOCARBON.html) ?

[huitang-earth]

@rosiealice ,@ckoven, @rgknox I have tried to modify a restart file and prescribed patch/cohort structures according to our inventory data. Techincally, it seem to work fine now. But there are some issues I encountered during the processes, and would appreciate for your comments to make the experiments more scientifically meaningful:

(1) How do you deal with LAI and plant phenology in ST3 mode? I would expect an option in FATES to prescribe LAI similar to the specified phenology mode in CLM. But so far, when ST3 is used, LAI will still be calculated based on some allometric relationship. When I prescribed the forest structure according to our inventory data, the LAI calculated by FATES is quite different from our observed LAI data. In addition, I notice that "phenology" is included in ST3 mode in the older model version, but not included in ST3 mode in the newest model version, I don't know why? If I want to prescribe LAI in FATES, do you have any suggestions how this should be done reasonably? Perhaps, tune some allometric parameter to fit the observed LAI?

(2) How the height of cohort from inventory data is used in FATES? I notice that in "FatesInventoryInitMod.F90", although "height of cohort" is a required variable, it is actually not used to create cohorts. The height of cohort is calculated based on allometric relation (h_allom). I am wondering why? Is it still useful to use our own data for the "height of cohort".

(3) I have been struggling with my ST3 mode results that the prescribed patch/cohort structure is actually changed during the simulation. Then, I notice that it is subroutine "ed_update_site" that cause the problem. After reading a restart file and each time step of "ed_ecosystem_dynamics", "ed_update_site" will be called. In this subroutine, "canopy_spread" is called to adjust "canopy spread" with a small value. This can give rise to a overcrowded canopy that some of the cohorts have to be demoted to understorey. I am wondering why "canopy_spread" is also used for ST3 mode? I would expect that in ST3 mode, the cohort/patch structure should be static.

(4) Which variables in restart file should be modify?

Here are the list of variables I have modified, with which I have done a successful simulation. Are there any variables missing in the list or any variables on the list that is not necessary to modify? Note that I use a relatively old FATES version so far and will test this on the newest FATES version.

fates_PatchesPerSite(:), # based on our own data
fates_CohortsPerPatch(:), # based on our own data
fates_age(:), # based on our own data
fates_area(:), # based on our own data
fates_pft(:), # based on our own data
fates_dbh(:), # based on our own data
fates_height(:), # based on our own data
fates_nplant(:), # based on our own data
fates_canopy_layer(:), # based on our own data, all the tree cohorts are considered as upper canopy
fates_canopy_layer_yesterday(:), # based on our own data
fates_isnew(:), # 0
fates_size_class_lasttimestep(:), # according to size class bin used in FATES
fates_spread(:), # 0.5
fates_canopy_trim(:), # 1.0
fates_status_coh(:), # according to phenology of each PFT
fates_laimemory(:), # according to "bleaf" used in FATES
leaf_c_val_001(:), # according allometric function used in FATES with our own dbh. This can produce overcrowded canopy. If this is the case, I adjust our tree density data to keep the canopy area smaller.
fnrt_c_val_001(:), # according allometric function used in FATES
sapw_c_val_001(:) # according allometric function used in FATES
store_c_val_001(:) # according allometric function used in FATES
struct_c_val_001(:) # according allometric function used in FATES
fates_solar_zenith_flag_pa(:) #copy from original restart file for each patch
fates_solar_zenith_angle_pa(:) #copy from original restart file for each patch
fates_gnd_alb_dir(:) #copy from original restart file for each patch
fates_gnd_alb_dif(:) #copy from original restart file for each patch

(5) I notice that when using FATES, the namelist parameter in CLM "create_crop_landunit" has to be set to FALSE. The namelist building will fail if I use "create_crop_landunit = .true." I don't know why this has to be the case. I test "create_crop_landunit = .true." with FATES, which will be comparable with our CLM5 simulation, and I don't notice any problem.

@rosiealice I am now using NCL script to modify the restart file. I am wondering which programming language are recommended for sharing within FATES community. If python is the common language, I can adapt my script to python. I will be happy to share my script on github and hopefully it can be constantly updated with the development of FATES before it can be hard coded into FATES.

[rosiealice]

Hi @huitang-earth. Happy New Year, and sorry for the delay. Great that you made some progress with this.

On 1), I know @ckoven had making a proscribed phenology version as an aspiration within the ST3 mode, but it's not implemented yet, perhaps because we are mostly using it for tropical systems without much in the way of phenology ;) He might have more to add about it. I think your stop-gap measure of adjusting allocation should work...

@rgknox would know the answer to #2. I guess one can't use both DBH -and- height to determine tree size when using a given allometric relationship, so we need to ignore one of them. Is that what you're getting at.

on #3 and #4 I am not sure of the answer. We'd need to look carefully at those. I think you're right that this shouldn't happen and maybe it didn't cause an issue in systems we already tested the code in?

On languages, I think we are converging to using python for these sorts of tools. I've been trying to get into using Jupyter notebooks (wherein @ckoven was ahead of the curve) along with many of the CLM folk too...

[ckoven]

Hi @huitang-earth just to say that for (3) if you set the min and max values for the crown area allometry (fates_allom_d2ca_coefficient_min and fates_allom_d2ca_coefficient_max in the parameter file) to be the same, then canopy spread changes will have no effect on the model dynamics. So that should solve that issue without having to rewrite any code?

For (1), yes this has been on the to-do list for a while but not yet been solved. For (4) and (5) I'm not sure. For (6) I've not actually modified restart files, only read from them, but in general I'd advocate for python (and possibly jupyter if its not meant to just be a command-line thing).

[rgknox]

Hi @huitang-earth ,

Regarding question 2, there are indeed several redundant and unused variables in that file format. That file format is specifically designed to be compatible with the ED2 model, that is why it seems so strange. It would be great to have another format that streamlined and specific to FATES.

[ckoven]

Closing this as an old discussion. @huitang-earth if you feel like this is still thing to do, feel free to open, but my sense I that the new reduced complexity configurations such as SP mode handle a lot of the science use cases that I think you were asking about.