diff --git a/images/challenge/FATES_complexity_modes.png b/images/challenge/FATES_complexity_modes.png new file mode 100644 index 000000000..aa21cbc92 Binary files /dev/null and b/images/challenge/FATES_complexity_modes.png differ diff --git a/images/challenge/FATES_schematic.png b/images/challenge/FATES_schematic.png new file mode 100644 index 000000000..ef4d2fc9e Binary files /dev/null and b/images/challenge/FATES_schematic.png differ diff --git a/images/diagnostics/clm_ctsm/fates_annual_gpp_diff.png b/images/diagnostics/clm_ctsm/fates_annual_gpp_diff.png new file mode 100644 index 000000000..22db28aad Binary files /dev/null and b/images/diagnostics/clm_ctsm/fates_annual_gpp_diff.png differ diff --git a/images/diagnostics/clm_ctsm/fates_annual_gpp_diff_pft.png b/images/diagnostics/clm_ctsm/fates_annual_gpp_diff_pft.png new file mode 100644 index 000000000..9c21cb22c Binary files /dev/null and b/images/diagnostics/clm_ctsm/fates_annual_gpp_diff_pft.png differ diff --git a/images/diagnostics/clm_ctsm/fates_annual_lh_diff.png b/images/diagnostics/clm_ctsm/fates_annual_lh_diff.png new file mode 100644 index 000000000..331d9db99 Binary files /dev/null and b/images/diagnostics/clm_ctsm/fates_annual_lh_diff.png differ diff --git a/images/diagnostics/clm_ctsm/fates_bart_canopy_structure.png b/images/diagnostics/clm_ctsm/fates_bart_canopy_structure.png new file mode 100644 index 000000000..fb6f5dd64 Binary files /dev/null and b/images/diagnostics/clm_ctsm/fates_bart_canopy_structure.png differ diff --git a/images/diagnostics/clm_ctsm/fates_bart_growth_mortality.png b/images/diagnostics/clm_ctsm/fates_bart_growth_mortality.png new file mode 100644 index 000000000..5db5ecaa2 Binary files /dev/null and b/images/diagnostics/clm_ctsm/fates_bart_growth_mortality.png differ diff --git a/images/diagnostics/clm_ctsm/fates_bart_layered_LAI_bg.png b/images/diagnostics/clm_ctsm/fates_bart_layered_LAI_bg.png new file mode 100644 index 000000000..a86504031 Binary files /dev/null and b/images/diagnostics/clm_ctsm/fates_bart_layered_LAI_bg.png differ diff --git a/images/diagnostics/clm_ctsm/fates_bart_layered_LAI_init.png b/images/diagnostics/clm_ctsm/fates_bart_layered_LAI_init.png new file mode 100644 index 000000000..523eb1d23 Binary files /dev/null and b/images/diagnostics/clm_ctsm/fates_bart_layered_LAI_init.png differ diff --git a/images/diagnostics/clm_ctsm/fates_bart_mortality.png b/images/diagnostics/clm_ctsm/fates_bart_mortality.png new file mode 100644 index 000000000..ccd984314 Binary files /dev/null and b/images/diagnostics/clm_ctsm/fates_bart_mortality.png differ diff --git a/images/diagnostics/clm_ctsm/fates_bart_mortality_by_size.png b/images/diagnostics/clm_ctsm/fates_bart_mortality_by_size.png new file mode 100644 index 000000000..8bf1e58c2 Binary files /dev/null and b/images/diagnostics/clm_ctsm/fates_bart_mortality_by_size.png differ diff --git a/images/diagnostics/clm_ctsm/fates_bart_mortality_snapshots.png b/images/diagnostics/clm_ctsm/fates_bart_mortality_snapshots.png new file mode 100644 index 000000000..5b75e1b65 Binary files /dev/null and b/images/diagnostics/clm_ctsm/fates_bart_mortality_snapshots.png differ diff --git a/images/diagnostics/clm_ctsm/fates_bart_patch_age.png b/images/diagnostics/clm_ctsm/fates_bart_patch_age.png new file mode 100644 index 000000000..628574373 Binary files /dev/null and b/images/diagnostics/clm_ctsm/fates_bart_patch_age.png differ diff --git a/images/diagnostics/clm_ctsm/fates_bart_promotion_demotion.png b/images/diagnostics/clm_ctsm/fates_bart_promotion_demotion.png new file mode 100644 index 000000000..746826c4c Binary files /dev/null and b/images/diagnostics/clm_ctsm/fates_bart_promotion_demotion.png differ diff --git a/images/diagnostics/clm_ctsm/fates_bart_stand_structure.png b/images/diagnostics/clm_ctsm/fates_bart_stand_structure.png new file mode 100644 index 000000000..b7d65e72a Binary files /dev/null and b/images/diagnostics/clm_ctsm/fates_bart_stand_structure.png differ diff --git a/images/diagnostics/clm_ctsm/fates_bart_vegc_composition_bg.png b/images/diagnostics/clm_ctsm/fates_bart_vegc_composition_bg.png new file mode 100644 index 000000000..65e3f0b3b Binary files /dev/null and b/images/diagnostics/clm_ctsm/fates_bart_vegc_composition_bg.png differ diff --git a/images/diagnostics/clm_ctsm/fates_bart_vegc_composition_init.png b/images/diagnostics/clm_ctsm/fates_bart_vegc_composition_init.png new file mode 100644 index 000000000..063a7777e Binary files /dev/null and b/images/diagnostics/clm_ctsm/fates_bart_vegc_composition_init.png differ diff --git a/images/diagnostics/clm_ctsm/fates_monthly_gpp_pft.png b/images/diagnostics/clm_ctsm/fates_monthly_gpp_pft.png new file mode 100644 index 000000000..481592c8a Binary files /dev/null and b/images/diagnostics/clm_ctsm/fates_monthly_gpp_pft.png differ diff --git a/images/diagnostics/clm_ctsm/fates_plot_monthly_gpp.png b/images/diagnostics/clm_ctsm/fates_plot_monthly_gpp.png new file mode 100644 index 000000000..f837e4e1c Binary files /dev/null and b/images/diagnostics/clm_ctsm/fates_plot_monthly_gpp.png differ diff --git a/notebooks/challenge/clm_ctsm/clm_exercise_3.ipynb b/notebooks/challenge/clm_ctsm/clm_exercise_3.ipynb index 21dc528d7..bfebd7b86 100644 --- a/notebooks/challenge/clm_ctsm/clm_exercise_3.ipynb +++ b/notebooks/challenge/clm_ctsm/clm_exercise_3.ipynb @@ -5,7 +5,7 @@ "id": "f406f992-92bd-4b17-9bd3-b99c5c8abaf3", "metadata": {}, "source": [ - "# 3: Modify input data" + "# 3: Run CLM with FATES and modify input data" ] }, { @@ -13,11 +13,59 @@ "id": "0037b73f-f174-48e7-8e4f-0744d7d23fe0", "metadata": {}, "source": [ - "We can modify the input to CLM by changing one of the plant functional type properties. We will then compare these results with the control experiment.\n", + "## FATES\n", "\n", - "Note that you will need to change a netcdf file for this exercise. Because netcdf are in binary format you will need a type of script or interperter to read the file and write it out again. (e.g. ferret, IDL, NCL, NCO, Perl, Python, Matlab, Yorick). Below in the solution we will show how to do this using NCO.\n", + "An important component of representing the land surface is accurately capturing vegetation dynamics and their impact on and interaction with the Earth system. Ecosystem demography models explicitly represent the size structure and successional state of vegetation, through direct simulation of plant growth, mortality, and regeneration. Thus, important vegetation characteristics such as vegetation canopy height, succession, and even potential biome shifts become emergent properties of the model rather than being prescribed.\n", "\n", - "NOTE: For any tasks other than setting up, building, submitting cases you should probably do these tasks on the Data Visualization Cluster - casper, and not on the derecho login nodes." + "**FATES** is the \"**F**unctionally **A**ssembled **T**errestrial **E**cosystem **S**imulator\". FATES is a cohort model of vegetation competition and co-existence, allowing a representation of the biosphere which accounts for the division of the land surface into successional stages, and for competition for light between height structured cohorts of representative trees of various plant functional types. Individual plants within FATES are grouped into \"cohorts\" of the same size and PFT, and these cohorts compete for light and resources on individual \"patches\" that represent different disturbance histories. This type of ecosystem heterogeneity is in contrast to the default vegetation model in CLM, which uses two (sunlit & shaded) \"big leaf\" canopies per PFT, each on their own patch, with no representation of within-canopy structural heterogeneity or disturbance history.\n", + "\n", + "When CLM is coupled to FATES (\"CLM-FATES\"), CLM provides site and soil conditions and atmospheric forcing, while FATES simulates plant physiological, vegetation demography, and biogeochemical processes.\n", + "\n", + "
\n", + "\"Conceptual\n", + "
\n", + "\n", + "*

Processes simulated in CLM-FATES by each model. Top: processes simulated by FATES when connected to CLM. Arrows in purple indicate conditions supplied to FATES by CLM. Arrows in green indicate conditions supplied to CLM by FATES. Bottom: Processes simulated by CLM when connected to FATES. Green starred variables are simulated and provided by FATES or in the case of aerodynamic resistance (ra) are influenced by the FATES-provided roughness length, displacement height, and leaf dimension. †: Only used in FATES hydraulics mode. From Foster et al. (202).

*" + ] + }, + { + "cell_type": "markdown", + "id": "d7001292-e327-4218-8089-867099549f32", + "metadata": {}, + "source": [ + "### FATES Complexity Modes\n", + "\n", + "Currently, FATES can be run in several different \"complexity modes\", where parts of the vegetation model are driven by input data rather than simulated. These modes can be used to facilitate calibration, test features, or run simulations more quickly. These modes are:\n", + "\n", + "1. **Satellite Phenology** (SP) mode: this mode is designed to run with leaf area index (LAI), stem area index (SAI), and canopy height (HTOP) as input to the model. As such, all processes that are normally used to calculate these values are turned off (e.g., mortality, allocation, etc.)\n", + "\n", + "2. **No-Competition Mode**: this mode runs with full complexity in terms of processes, but places each FATES PFT on its own patch. As such, PFTs do not compete with one another. The patch area of each PFT is determined from the input CLM surface dataset. However, please note that the PFTs in the FATES parameter file do not always map one-to-one with the CLM PFTs on the surface dataset. See the FATES parameter *fates_hlm_pft_map* on the FATES parameter file for the correct mapping of FATES to CLM PFTs.\n", + "\n", + "3. **Fixed Biogeography Mode**: this mode turns off prognostic spatial changes in the distribution of vegetation and instead, the model uses input data to determine which PFTs are present at any given gridcell. The PFT composition in each gridcell is derived from the input CLM surface dataset. \n", + "\n", + "4. **Full FATES Mode**: All processes are turned on an PFTs are allowed to grow anywhere.\n", + "\n", + "Note that there are different combinations of no-competition and fixed biogeography mode that will result in different model behaviors. See the FATES namelist documentation for these options.\n", + "\n", + "
\n", + "\"FATES\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "b322a528-9e23-45b9-aaed-b05dfe0292f5", + "metadata": {}, + "source": [ + "We will be running FATES in Satellite Phenology mode. This can be done by either choosing an SP compset (e.g. ) or manually setting this via the `user_nl_clm` file, with the parameters `use_fates_sp`, `use_fates_nocomp`, `use_fates_fixed_biogeog`. Check out the online documentation for more information on these parameters." + ] + }, + { + "cell_type": "markdown", + "id": "6afa3c46-9939-46bf-9d37-d5a8874abd6e", + "metadata": {}, + "source": [ + "## FATES Control Case" ] }, { @@ -26,14 +74,13 @@ "metadata": {}, "source": [ "
\n", - "Exercise: Run an experimental case

\n", + "Exercise: Run a FATES case

\n", " \n", - "Create a case called **i.day5.a_pft** using the compset `I2000Clm50Sp` at `f09_g17_gl4` resolution. \n", + "Create a case called **i.fates.year1** using the compset `I2000Clm60FatesSpCrujraRsGs` at `f19_f19_mt233` resolution. \n", "\n", - "Look at variable “rholvis” in the forcing file using ncview or ncdump –v rholvis. This is the visible leaf reflectance for every pft. Modify the rholvis parameter to .\n", - "`/glade/campaign/cesm/cesmdata/cseg/inputdata/lnd/clm2/paramdata/clm5_params.c171117.nc`\n", - " \n", - "Set the run length to **5 days**. \n", + "Set the run length to **1 year**.\n", + "\n", + "Set the `finidat` appropriately.\n", "\n", "Build and run the model.\n", " \n", @@ -42,132 +89,413 @@ }, { "cell_type": "markdown", - "id": "f639e182-f48a-431c-a594-9c34323417eb", + "id": "dbb791d8-e002-4f48-8bc7-eac232cee6b5", "metadata": {}, "source": [ + "
\n", + "
\n", "\n", + " Click here for hints \n", "\n", - "
\n", - "
\n", - "Click here for the solution
\n", - " \n", - "Create a clone from the control experiment i.day5.a_pft :\n", + "
\n", + "\n", + "**How do I compile?**\n", + "\n", + "You can compile with the command:\n", "```\n", - "cd /glade/u/home/$USER/code/my_cesm_code/cime/scripts\n", - "./create_clone --case ~/cases/i.day5.a_pft --clone ~/cases/i.day5.a\n", + "qcmd -- ./case.build\n", "```\n", + "\n", "
\n", "\n", - "Modify the rholvis parameter in the physiology file:\n", - "``` \n", - "cd /glade/derecho/scratch/$USER\n", - "cp /glade/campaign/cesm/cesmdata/cseg/inputdata/lnd/clm2/paramdata/clm5_params.c171117.nc .\n", - "chmod u+w clm5_params.c171117.nc\n", - "cp clm5_params.c171117.nc clm5_params.c171117.new.nc\n", - "ncap2 -A -v -s 'rholvis(4)=0.4' clm5_params.c171117.nc clm5_params.c171117.new.nc\n", - "```\n", + "**How do I change the run length to 1 year?**\n", + "\n", + "Use xml variables: ``STOP_OPTION`` and ``STOP_N``. \n", + "\n", "
\n", "\n", - "Check the new rholvis parameter to be sure the modification worked:\n", - "``` \n", - "ncdump -v rholvis clm5_params.c171117.new.nc\n", - "# and compare it to the original file\n", - "ncdiff clm5_params.c171117.nc clm5_params.c171117.new.nc ncdiff.nc\n", - "ncdump -v rholvis ncdiff.nc\n", + "**How do I check my solution?**\n", + "\n", + "When your run is completed, go to the archive directory and navigate to the subdirectory `lnd/hist`\n", + "\n", + "```\n", + "cd /glade/derecho/scratch/$USER/archive/i.fates.year1\n", + "cd lnd/hist\n", + "```\n", + "\n", + "(1) Check that your archive directory contains the files:\n", + "\n", + "```\n", + "i.fates.year1.clm2.h0a.2000-01.nc i.fates.year1.clm2.h0i.2000-01.nc\n", + "i.fates.year1.clm2.h0a.2000-02.nc i.fates.year1.clm2.h0i.2000-02.nc\n", + "i.fates.year1.clm2.h0a.2000-03.nc i.fates.year1.clm2.h0i.2000-03.nc\n", + "i.fates.year1.clm2.h0a.2000-04.nc i.fates.year1.clm2.h0i.2000-04.nc\n", + "i.fates.year1.clm2.h0a.2000-05.nc i.fates.year1.clm2.h0i.2000-05.nc\n", + "i.fates.year1.clm2.h0a.2000-06.nc i.fates.year1.clm2.h0i.2000-06.nc\n", + "i.fates.year1.clm2.h0a.2000-07.nc i.fates.year1.clm2.h0i.2000-07.nc\n", + "i.fates.year1.clm2.h0a.2000-08.nc i.fates.year1.clm2.h0i.2000-08.nc\n", + "i.fates.year1.clm2.h0a.2000-09.nc i.fates.year1.clm2.h0i.2000-09.nc\n", + "i.fates.year1.clm2.h0a.2000-10.nc i.fates.year1.clm2.h0i.2000-10.nc\n", + "i.fates.year1.clm2.h0a.2000-11.nc i.fates.year1.clm2.h0i.2000-11.nc\n", + "i.fates.year1.clm2.h0a.2000-12.nc i.fates.year1.clm2.h0i.2000-12.nc\n", + "```\n", + "\n", + "(2) Investigate the contents of one of the files with ``ncdump``.\n", + "\n", + "```\n", + "ncdump -h i.fates.year1.clm2.h0a.2000-01.nc\n", + "```\n", + "\n", + "(3) Check to make sure we see some FATES variables:\n", + "\n", + "```\n", + "ncdump -h i.fates.year1.clm2.h0a.2000-01.nc | grep \"FATES\"\n", + "```\n", + "\n", + "
\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "f639e182-f48a-431c-a594-9c34323417eb", + "metadata": {}, + "source": [ + "
\n", + "
\n", + "Click here for the solution
\n", + "\n", + "**Create a new case**\n", + "\n", + "Create a new case i.fates.year1.a with the command:\n", + "```\n", + "cd /glade/u/home/$USER/code/my_cesm_code/cime/scripts\n", + "./create_newcase --case ~/cases/i.fates.year1 --compset I2000Clm60FatesSpCrujraRsGs --res f19_f19_mt233 --run-unsupported\n", "```\n", "
\n", + "\n", + "**Case setup**\n", + " \n", + "Invoke case.setup with the command:\n", " \n", - "Case setup:\n", "``` \n", - "cd ~/cases/i.day5.a_pft\n", + "cd ~/cases/i.fates.year1\n", "./case.setup\n", "```\n", "
\n", "\n", - "Change the clm namelist using user_nl_clm to point at the modified file. Add the following line:\n", + "Check the namelist by running:\n", "``` \n", - "paramfile = '/glade/derecho/scratch/$USER/clm5_params.c171117.new.nc' \n", + "./preview_namelists\n", "```\n", - "
\n", + "
\n", + "\n", + "**Set run length**\n", " \n", - "Check the namelist by running:\n", + "Change the `run length`:\n", "``` \n", - "./preview_namelists\n", + "./xmlchange STOP_N=1,STOP_OPTION=nyears\n", "```\n", + "\n", "
\n", "\n", - "If needed, change job queue, account number, or wallclock time. \n", - "For instance:\n", + "**Change the job queue and account number**\n", + "\n", + "If needed, change `job queue` and `account number`.
\n", + "For instance, to run in the queue `tutorial` and the project number ``UESM0015`` (You should use the project number given for this tutorial), use the command:\n", "``` \n", - "./xmlchange JOB_QUEUE=tutorial,PROJECT=UESM0015 --force,JOB_WALLCLOCK_TIME=0:15:00\n", + "./xmlchange JOB_QUEUE=tutorial,PROJECT=UESM0015 --force\n", "```\n", + "\n", "
\n", "\n", - "Build case:\n", + "**Change the runtime length**\n", + "\n", + "Change wallclock time to use the maximum of the allowed wallclock on Derecho:\n", "```\n", - "qcmd -- ./case.build\n", + "./xmlchange --subgroup case.run JOB_WALLCLOCK_TIME=12:00:00\n", "```\n", "
\n", - " \n", - "Compare the namelists from the two experiments:\n", + "\n", + "**Build case:**\n", "```\n", - "diff CaseDocs/lnd_in ../i.day5.a/CaseDocs/lnd_in\n", + "qcmd -- ./case.build\n", "```\n", "
\n", " \n", - "Submit case:\n", + " \n", + "**Submit case:**\n", "```\n", "./case.submit\n", "```\n", "
\n", "\n", + "**Check the run:**\n", "When the run is completed, look into the archive directory for: \n", - "i.day5.a. \n", + "i.fates.year1.a. \n", " \n", - "(1) Check that your archive directory on derecho (The path will be different on other machines): \n", + "(1) Check that your archive directory on derecho (the path will be different on other machines): \n", "```\n", - "cd /glade/derecho/scratch/$USER/archive/i.day5.a_pft/lnd/hist\n", + "cd /glade/derecho/scratch/$USER/archive/i.fates.year1/lnd/hist\n", "\n", "ls \n", "```\n", "
\n", "\n", - "(2) Compare to control run:\n", + "(2) Check that your archive directory contains the files:\n", + "\n", "```\n", - "ncdiff i.day5.a_pft.clm2.XXX.nc /glade/derecho/scratch/$USER/archive/i.day5.a/lnd/hist/i.day5.a.clm2.XXX.nc i_diff.nc\n", + "i.fates.year1.clm2.h0a.2000-01.nc i.fates.year1.clm2.h0i.2000-01.nc\n", + "i.fates.year1.clm2.h0a.2000-02.nc i.fates.year1.clm2.h0i.2000-02.nc\n", + "i.fates.year1.clm2.h0a.2000-03.nc i.fates.year1.clm2.h0i.2000-03.nc\n", + "i.fates.year1.clm2.h0a.2000-04.nc i.fates.year1.clm2.h0i.2000-04.nc\n", + "i.fates.year1.clm2.h0a.2000-05.nc i.fates.year1.clm2.h0i.2000-05.nc\n", + "i.fates.year1.clm2.h0a.2000-06.nc i.fates.year1.clm2.h0i.2000-06.nc\n", + "i.fates.year1.clm2.h0a.2000-07.nc i.fates.year1.clm2.h0i.2000-07.nc\n", + "i.fates.year1.clm2.h0a.2000-08.nc i.fates.year1.clm2.h0i.2000-08.nc\n", + "i.fates.year1.clm2.h0a.2000-09.nc i.fates.year1.clm2.h0i.2000-09.nc\n", + "i.fates.year1.clm2.h0a.2000-10.nc i.fates.year1.clm2.h0i.2000-10.nc\n", + "i.fates.year1.clm2.h0a.2000-11.nc i.fates.year1.clm2.h0i.2000-11.nc\n", + "i.fates.year1.clm2.h0a.2000-12.nc i.fates.year1.clm2.h0i.2000-12.nc\n", + "```\n", + "\n", + "(2) Investigate the contents of one of the files with ``ncdump``.\n", "\n", - "ncview i_diff.nc\n", + "```\n", + "ncdump -h i.fates.year1.clm2.h0a.2000-01.nc\n", "```\n", "\n", + "(3) Check to make sure we see some FATES variables:\n", "\n", + "```\n", + "ncdump -h i.fates.year1.clm2.h0a.2000-01.nc | grep \"FATES\"\n", + "```\n", "
\n", "
" ] }, { "cell_type": "markdown", - "id": "b3ffd3cc-676e-4e7c-9ff4-cd65d4745397", - "metadata": { - "tags": [] - }, + "id": "58f9a16c-9845-40ed-be33-6c6b61b7449a", + "metadata": {}, + "source": [ + "## Modify FATES input data" + ] + }, + { + "cell_type": "markdown", + "id": "e4437454-21ab-4d49-922b-b2b98042c7cf", + "metadata": {}, "source": [ - "## Test your understanding" + "We can modify the input to CLM-FATES by changing one of the plant functional type properties. We will then compare these results with the control experiment.\n", + "\n", + "Instead of using a netcdf parameter file, FATES uses a json parameter file as its input. This makes it easy to modify with any text editing software." + ] + }, + { + "cell_type": "markdown", + "id": "3b160a06-3a78-4a37-9a2c-40618fde1ed8", + "metadata": {}, + "source": [ + "
\n", + "Exercise: Run an experimental case

\n", + " \n", + "Create a case called **i.fates.year1.vcmax** using the compset `I2000Clm60FatesSpCrujraRsGs` at `f19_f19_mt233` resolution.\n", + "\n", + "Look at variable “fates_leaf_vcmax25top” in the FATES parameter file (located in the FATES source directory: `/glade/u/home/$USER/code/my_cesm_code/CTSM/src/fates/parameter_files/fates_params_default.json`).\n", + "\n", + "This is the maximum carboxylation rate of Rubisco at 25ºC, at canopy top, and is indexed by FATES PFT and leaf age class. Modify the `fates_leaf_vcmax25top` parameter of the first PFT to 20.0.\n", + "\n", + "You can either modify the parameter file in place, or you can copy it (with `cp`), modify the copy, and tell the model to use the updated parameter file. This is done via the `user_nl_clm` file by setting `fates_paramfile=my/new/parameter_file`.\n", + "\n", + "Set the run length to **1 year**. \n", + "\n", + "Build and run the model.\n", + " \n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "c46470ce-f6ec-4715-a75d-3f54599ea448", + "metadata": {}, + "source": [ + "
\n", + "
\n", + "\n", + " Click here for hints \n", + "\n", + "
\n", + "\n", + "**How do I compile?**\n", + "\n", + "You can compile with the command:\n", + "```\n", + "qcmd -- ./case.build\n", + "```\n", + "\n", + "
\n", + "\n", + "**How do I modify the FATES parameter?**\n", + "\n", + "FATES automatically uses the parameter file in the FATES source code: `/glade/u/home/$USER/code/my_cesm_code/CTSM/src/fates/parameter_files/fates_params_default.json`. You can copy this parameter file to a new file:\n", + "\n", + "`cp /glade/u/home/$USER/code/my_cesm_code/CTSM/src/fates/parameter_files/fates_params_default.json ~/cases/fates_params_vcmax_update.json`\n", + "\n", + "You can open this file using any text editing software. Use Ctrl+F to find the `fates_leaf_vcmax25top` parameter.\n", + "
\n", + "\n", + "**How do I tell the model to use the new parameter file?**\n", + "\n", + "Set the `fates_paramfile` parameter in the `user_nl_clm` file to your new parameter file path.\n", + "
\n", + "\n", + "\n", + "**How do I change the run length to 1 year?**\n", + "\n", + "Use xml variables: ``STOP_OPTION`` and ``STOP_N``. \n", + "\n", + "
\n", + "\n", + "**How do I check my solution?**\n", + "\n", + "When your run is completed, go to the archive directory and navigate to the subdirectory `lnd/hist`\n", + "\n", + "```\n", + "cd /glade/derecho/scratch/$USER/archive/i.fates.year1.vcmax\n", + "cd lnd/hist\n", + "ls\n", + "```\n", + "
\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "ff9b1b70-4c64-4ba1-900c-e47575f47324", + "metadata": {}, + "source": [ + "
\n", + "
\n", + "Click here for the solution
\n", + "\n", + "**Clone a new case**\n", + "\n", + "Create a clone from the control experiment i.fates.year1.a_vcmax:\n", + "```\n", + "cd /glade/u/home/$USER/code/my_cesm_code/cime/scripts\n", + "./create_clone --case ~/cases/i.fates.year1.vcmax --clone ~/cases/i.fates.year1\n", + "```\n", + "
\n", + "\n", + "**Modify the FATES parameter file**\n", + "\n", + "Copy the default parameter file to a new parameter file.\n", + "``` \n", + "cd /glade/u/home/$USER/code/my_cesm_code/CTSM/src/fates/parameter_files/\n", + "cp fates_params_default.json ~/cases/fates_params_vcmax_update.json\n", + "```\n", + "\n", + "Modify the fates_leaf_vcmax25top parameter in the FATES parameter file.\n", + "\n", + "Edit the parameter file by opening it with any text editing software. Update the fates_leaf_vcmax25top parameter in the first index location to be 20.0:\n", + "``` \n", + "\"fates_leaf_vcmax25top\": {\n", + " \"dtype\": \"float\",\n", + " \"dims\": [\"fates_leafage_class\", \"fates_pft\"],\n", + " \"long_name\": \"maximum carboxylation rate of Rub. at 25C, canopy top\",\n", + " \"units\": \"umol CO2/m^2/s\",\n", + " \"data\": [[20.0, 62.0, 39.0, 61.0, 58.0, 58.0, 62.0, 54.0, 54.0, 38.0, 54.0, 86.0, 78.0, 78.0]]\n", + "},\n", + "```\n", + "\n", + "
\n", + "\n", + "**Case setup**\n", + "``` \n", + "cd ~/cases/i.fates.year1.vcmax\n", + "./case.setup\n", + "```\n", + "
\n", + "\n", + "**Update the user_nl_clm file to use the new parameter file**\n", + "\n", + "Modify the `user_nl_clm` for the case and add the following line (or wherever you put it):\n", + "\n", + "`fates_paramfile=/glade/home/$USER/cases/fates_params_vcmax_update.json`\n", + " \n", + "Check the namelist by running:\n", + "``` \n", + "./preview_namelists\n", + "```\n", + "
\n", + "\n", + "**Set run length**\n", + " \n", + "Change the `run length`:\n", + "``` \n", + "./xmlchange STOP_N=1,STOP_OPTION=nyears\n", + "```\n", + "\n", + "
\n", + "\n", + "**Change the job queue and account number**\n", + "\n", + "If needed, change `job queue` and `account number`.
\n", + "For instance, to run in the queue `tutorial` and the project number ``UESM0015`` (You should use the project number given for this tutorial), use the command:\n", + "``` \n", + "./xmlchange JOB_QUEUE=tutorial,PROJECT=UESM0015 --force\n", + "```\n", + "\n", + "
\n", + "\n", + "**Change the runtime length**\n", + "\n", + "Change wallclock time to use the maximum of the allowed wallclock on Derecho:\n", + "```\n", + "./xmlchange --subgroup case.run JOB_WALLCLOCK_TIME=12:00:00\n", + "```\n", + "
\n", + "\n", + "**Build case:**\n", + "```\n", + "qcmd -- ./case.build\n", + "```\n", + "
\n", + " \n", + "**Submit case:**\n", + "```\n", + "./case.submit\n", + "```\n", + "
\n", + "\n", + "\n", + "**Check the run:**\n", + "When the run is completed, look into the archive directory for: \n", + "i.fates.year1.a. \n", + " \n", + "Check that your archive directory on derecho (the path will be different on other machines): \n", + "```\n", + "cd /glade/derecho/scratch/$USER/archive/i.fates.year1.vcmax/lnd/hist\n", + "\n", + "ls \n", + "```\n", + "
\n", + "
" ] }, { "cell_type": "markdown", - "id": "03ac3664-5360-45d7-a3ad-0797a839a1d3", + "id": "ce042041-de2c-4de5-b384-5a2bbb22c954", "metadata": {}, "source": [ - "- How did rholvis change (increase/decrease)? Given this, what do you expect the model response to be?\n", - "- What changes do you see from the control case with the modified rholvis parameter?\n", - "- ... OTHERS? " + "We will investigate these runs in the next exercise." ] }, { "cell_type": "code", "execution_count": null, - "id": "27f332b2-5799-43a8-9060-50315ebdf6dc", + "id": "23b5c799-e502-4114-b29f-cfeeb421c893", "metadata": {}, "outputs": [], "source": [] @@ -175,9 +503,9 @@ ], "metadata": { "kernelspec": { - "display_name": "CMIP6 2019.10", + "display_name": "Python 3 (ipykernel)", "language": "python", - "name": "cmip6-201910" + "name": "python3" }, "language_info": { "codemirror_mode": { @@ -189,7 +517,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.7.10" + "version": "3.10.13" } }, "nbformat": 4, diff --git a/notebooks/challenge/clm_ctsm/clm_exercise_4.ipynb b/notebooks/challenge/clm_ctsm/clm_exercise_4.ipynb new file mode 100644 index 000000000..6c9abacbe --- /dev/null +++ b/notebooks/challenge/clm_ctsm/clm_exercise_4.ipynb @@ -0,0 +1,547 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# 4: Investigating Output" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**BEFORE BEGINNING THIS EXERCISE** - Check that your kernel (upper right corner, above) is `NPL 2023a`. This should be the default kernel, but if it is not, click on that button and select `NPL 2023a`." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "_______________\n", + "This activity was developed primarily by Adrianna Foster." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## CLM output\n", + "\n", + "While CLM-FATES and regular CLM share many output variables (e.g. `ASA`, `FSR`, `FSH`, etc.), most of the vegetation-related variables are now being simulated, calculated, and output by FATES. The variable names and their structure (as well as the PFTs!) are slightly different from a \"big leaf\" CLM run.\n", + "\n", + "Here we will show you how to look at some basic FATES output for the global cases you just ran." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "*As before, start by loading some packages*" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import glob\n", + "import xarray as xr\n", + "import functools\n", + "\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Helper functions\n", + "\n", + "These are a few helper functions that we will use in this notebook." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def preprocess(data_set: xr.Dataset, data_vars: list[str]) -> xr.Dataset:\n", + " \"\"\"Preprocesses an xarray Dataset by subsetting to specific variables - to be used on read-in\n", + "\n", + " Args:\n", + " data_set (xr.Dataset): input Dataset\n", + "\n", + " Returns:\n", + " xr.Dataset: output Dataset\n", + " \"\"\"\n", + "\n", + " return data_set[data_vars]\n", + "\n", + "\n", + "def annual_sum(raw_values: xr.DataArray, conversion_factor: float = 1.0) -> xr.DataArray:\n", + " \"\"\"Computes annual sum\n", + "\n", + " Args:\n", + " raw_values (xr.DataArray): input raw data\n", + " conversion_factor (float, optional): conversion factor. Defaults to 1.0\n", + "\n", + " Returns:\n", + " xr.DataArray: annual sum output\n", + " \"\"\"\n", + "\n", + " months = raw_values[\"time.daysinmonth\"]\n", + " return conversion_factor * (months * raw_values).groupby(\"time.year\").sum()\n", + "\n", + "\n", + "def annual_mean(raw_values: xr.DataArray, conversion_factor: float = 1.0) -> xr.DataArray:\n", + " \"\"\"Computes weighted annual mean using daysinmonth for missing-aware inputs.\n", + "\n", + " Args:\n", + " raw_values (xr.DataArray): input raw data\n", + " conversion_factor (float, optional): conversion factor. Defaults to 1.0\n", + "\n", + " Returns:\n", + " xr.DataArray: annual mean output\n", + " \"\"\"\n", + "\n", + " months = raw_values[\"time.daysinmonth\"]\n", + " \n", + " # multiply by number of days in month and conversion factor\n", + " weighted = (raw_values * conversion_factor) * months\n", + "\n", + " # compute number of valid days per year\n", + " valid_days = months.where(raw_values.notnull())\n", + "\n", + " # group and sum weighted data and valid days\n", + " ann_sum = weighted.groupby(\"time.year\").sum(dim=\"time\", skipna=True)\n", + " days_per_year = valid_days.groupby(\"time.year\").sum(dim=\"time\", skipna=True)\n", + "\n", + " return ann_sum / days_per_year.where(days_per_year > 0.0)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1. Reading and formatting data\n", + "\n", + "**Note**: the drop-down solutions, below, assume you used i.fates.year1.a and i.fates.year1.a_vcmax output for plotting for this section." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 1.1 Read in the data\n", + "The first step is to grab the history files from the runs you completed in the FATES challenge.\n", + "\n", + "For this example we will use:\n", + "- gross primary production (`FATES_GPP` and `FATES_GPP_PF`), \n", + "- latent heat flux (`EFLX_LH_TOT`), and\n", + "- sensible heat flux (`FSH`)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "NOTE: These are the raw history files that CTSM writes out. \n", + "\n", + "By default, they include grid cell averaged monthly means for different state and flux variables.\n", + "\n", + "
\n", + " TIP: If you want to look at other variables, the data_vars variable in the cell below is where you can modify what we're reading off of the CLM history files.\n", + "
\n", + "\n", + "#### Printing information about the dataset is helpful for understanding your data. \n", + "- *What dimensions do your data have?*\n", + "- *What are the coordinate variables?*\n", + "- *What variables are we looking at?*\n", + "- *Is there other helpful information, or are there attributes in the dataset we should be aware of?*\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# note: change the user here to your user name\n", + "user = 'afoster'\n", + "\n", + "# this should be the history directory for the control directory and the one with vcmax updated\n", + "control_hist = f\"/glade/derecho/scratch/{user}/archive/i.fates.year1.a/lnd/hist\"\n", + "vcmax_hist = f\"/glade/derecho/scratch/{user}/archive/i.fates.year1.a_vcmax/lnd/hist\"\n", + "\n", + "# create path to all files, including unix wild card for all dates\n", + "control_files = sorted(glob.glob(os.path.join(control_hist, \"i.fates.year1.a.clm2.h0a.*\")))\n", + "vcmax_files = sorted(glob.glob(os.path.join(vcmax_hist, \"i.fates.year1.a_vcmax.clm2.h0a.*\")))\n", + "\n", + "# read in files as xarray datasets:\n", + "data_vars = ['FATES_GPP', 'FATES_GPP_PF', 'EFLX_LH_TOT', 'FSH',\n", + " 'FATES_FRACTION', 'area', 'landfrac']\n", + "\n", + "ds_control = xr.open_mfdataset(\n", + " control_files,\n", + " preprocess=functools.partial(preprocess, data_vars=data_vars)\n", + ")\n", + "ds_vcmax = xr.open_mfdataset(\n", + " vcmax_files, \n", + " preprocess=functools.partial(preprocess, data_vars=data_vars)\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# print information about one of the datasets\n", + "ds_control" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "You can also print information about the variables in your dataset. The example below prints information about one of the data variables (`FATES_GPP_PF`) we read in. You can modify this cell to look at some of the other variables in the dataset.\n", + "\n", + "*What are the units, long name, and dimensions of your data?*" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ds_control.FATES_GPP_PF" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 1.2 Simple Calculations\n", + "To begin with, we need to modify all `FATES_` variables by multiplying them by `FATES_FRACTION`. This is because the FATES model doesn't know how much of each gridcell is occupied by natural vegetation, and so the grid-cell averaging assumings all of it is. We need to correct for this by multiplying by the variable `FATES_FRACTION`:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ds_control['GPP'] = ds_control.FATES_GPP*ds_control.FATES_FRACTION\n", + "ds_control['GPP'].attrs = ds_control.FATES_GPP.attrs\n", + "\n", + "ds_control['GPP_PF'] = ds_control.FATES_GPP_PF*ds_control.FATES_FRACTION\n", + "ds_control['GPP_PF'].attrs = ds_control.FATES_GPP_PF.attrs\n", + "\n", + "ds_vcmax['GPP'] = ds_vcmax.FATES_GPP*ds_vcmax.FATES_FRACTION\n", + "ds_vcmax['GPP'].attrs = ds_vcmax.FATES_GPP.attrs\n", + "\n", + "ds_vcmax['GPP_PF'] = ds_vcmax.FATES_GPP_PF*ds_vcmax.FATES_FRACTION\n", + "ds_vcmax['GPP_PF'].attrs = ds_vcmax.FATES_GPP_PF.attrs" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---\n", + "## 2. Basic plotting with global FATES\n", + "### 2.1 Easy plots using xarray\n", + "To get a first look at the data, we can plot a year of data from the simulation.\n", + "\n", + "
\n", + " NOTE: The plotting function only works with 1D or 2D data. Our data are 3D (time, lat, lon), so we need to specify a specific value `x`, `y`, and `col`.\n", + "
\n", + "\n", + "- We will plot GPP (variable = `GPP`). Note that we select the variable by specifying our dataset, `ds_control`, and the variable. \n", + "- This plotting function will plot `GPP` for each simulation in our dataset, and we have set the column to be month (`col=\"time\"`).\n", + "\n", + "*More plotting examples are on the [xarray web site](https://docs.xarray.dev/en/latest/user-guide/plotting.html)*" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ds_control.GPP.plot(x='lon', y='lat', col=\"time\", col_wrap=6);" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "
\n", + "Click here for the solution
\n", + "\n", + "![plot example](../../../images/diagnostics/clm_ctsm/fates_plot_monthly_gpp.png)\n", + "\n", + "*

Figure: Plotting solution for monthly GPP for our control simulation.

*\n", + " \n", + "
\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Question\n", + "\n", + "What do you notice in the Northern Hemisphere over the course of the year?" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 2.2 Calculating differences\n", + "\n", + "We can calculate the differences between our control run with FATES and the one in which we updated the vcmax parameter for PFT 1.\n", + "\n", + "The below code:\n", + "- Calculates annual GPP for both simulations\n", + "- Defines the difference as a new variable, `gpp_diff`" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "gpp_cf = 24*60*60 # kgC/m2/yr" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "gpp_ann_control = annual_sum(ds_control.GPP, gpp_cf)*ds_control.landfrac.isel(time=0)\n", + "gpp_ann_vcmax = annual_sum(ds_vcmax.GPP, gpp_cf)*ds_vcmax.landfrac.isel(time=0)\n", + "gpp_diff = gpp_ann_vcmax - gpp_ann_control\n", + "gpp_diff.attrs = {'units': 'kgC m-2 yr-1', 'long_name': 'annual gross primary production difference'}\n", + "gpp_diff.plot(cmap='RdBu_r', vmin=-1.5, vmax=1.5);" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "
\n", + "Click here for the solution
\n", + "\n", + "![plot example](../../../images/diagnostics/clm_ctsm/fates_annual_gpp_diff.png)\n", + "\n", + "*

Figure: Plotting solution for annual GPP for our updated solution minus the control simulation.

*\n", + " \n", + "
\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Questions\n", + "- How is GPP different for our updated simulation, relative to the control simulation?\n", + "- Where are the differences the largest?\n", + "- Where are they absent?\n", + "- What is causing these GPP changes in different regions and why is it only in some areas?" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Do the differences affect all PFTs?\n", + "To find out, we can plot the `GPP_PF` variable, which is indexed by PFT.\n", + "\n", + "We will also grab June (`.isel(time=5)`; 0-indexed):" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ds_control.GPP_PF.isel(time=5).plot(x='lon', y='lat', col=\"fates_levpft\", col_wrap=6);" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "
\n", + "Click here for the solution
\n", + "\n", + "![plot example](../../../images/diagnostics/clm_ctsm/fates_monthly_gpp_pft.png)\n", + "\n", + "*

Figure: Plotting solution for comparing GPP by PFT for our updated simulation from our control.

*\n", + " \n", + "
\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Do the differences\n", + "To find out, do the same calculation as above, but with the `GPP_PF` variable." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "gpp_ann_control_pf = annual_sum(ds_control.GPP_PF, gpp_cf)*ds_control.landfrac.isel(time=0)\n", + "gpp_ann_vcmax_pf = annual_sum(ds_vcmax.GPP_PF, gpp_cf)*ds_vcmax.landfrac.isel(time=0)\n", + "gpp_diff_pf = gpp_ann_vcmax_pf - gpp_ann_control_pf\n", + "gpp_diff_pf.attrs = {'units': 'kgC m-2 yr-1', 'long_name': 'annual gross primary production difference'}\n", + "gpp_diff_pf.plot(x='lon', y='lat', col=\"fates_levpft\", col_wrap=6)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "
\n", + "Click here for the solution
\n", + "\n", + "![plot example](../../../images/diagnostics/clm_ctsm/fates_annual_gpp_diff_pft.png)\n", + "\n", + "*

Figure: Plotting solution for comparing GPP by PFT for our updated simulation from our control.

*\n", + " \n", + "
\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Questions:\n", + "- Which PFT(s) show any difference? Why is this?\n", + "- Why don't other PFTs show any differences? Shouldn't there be some kind of compensating factor?\n", + "\n", + "*Hint: you can check the `fates_pftname` variable in `/glade/u/home/$USER/code/my_cesm_code/CTSM/src/fates/parameter_files/fates_params_default.json` to see which PFT corresponds to which index.*" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Do the differences affect other variables?\n", + "To find out, we can plot `EFLX_LH_TOT`, which is latent heat flux. How do you think this will be affected?\n", + "\n", + "`EFLX_LH_TOT` is a regular CLM variable (that is, it is calculated by CLM, not FATES, so we should not multiply it by `FATES_FRACTION`)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ds_control.EFLX_LH_TOT" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Do the differences\n", + "Let's calculate annual latent heat flux. This time we will use annual mean instead of an annual sum:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "lh_ann_control = annual_mean(ds_control.EFLX_LH_TOT)*ds_control.landfrac.isel(time=0)\n", + "lh_ann_vcmax = annual_mean(ds_vcmax.EFLX_LH_TOT)*ds_vcmax.landfrac.isel(time=0)\n", + "lh_diff = lh_ann_vcmax - lh_ann_control\n", + "lh_diff.attrs = {'units': 'W m-2', 'long_name': 'annual latent heat flux difference'}\n", + "lh_diff.plot(vmin=-20, vmax=20, cmap='RdBu_r')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "
\n", + "Click here for the solution
\n", + "\n", + "![plot example](../../../images/diagnostics/clm_ctsm/fates_annual_lh_diff.png)\n", + "\n", + "*

Figure: Plotting solution for comparing latent heat flux for our updated simulation from our control.

*\n", + " \n", + "
\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Questions:\n", + "- How is latent heat flux different?\n", + "- Try looking at sensible heat flux (`FSH`). What differences do you see? \n", + "\n", + "*Note that you might want to change the minimun (`vmin`) and maximum (`vmax`) colorbar values for the plot when you change variables*" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.13" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/notebooks/challenge/clm_ctsm/clm_exercise_5.ipynb b/notebooks/challenge/clm_ctsm/clm_exercise_5.ipynb new file mode 100644 index 000000000..600545c4a --- /dev/null +++ b/notebooks/challenge/clm_ctsm/clm_exercise_5.ipynb @@ -0,0 +1,1301 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# 5: Investigating vegetation demographics with FATES" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**BEFORE BEGINNING THIS EXERCISE** - Check that your kernel (upper right corner, above) is `NPL 2023a`. This should be the default kernel, but if it is not, click on that button and select `NPL 2023a`." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "_______________\n", + "This activity was developed primarily by Adrianna Foster." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## FATES output\n", + "\n", + "As we described above, when run in higher complexity modes (i.e., not SP mode) FATES simulates more than just GPP and other biophysical fluxes. It also simulates vegetation demography (e.g., growth, allocation, competition, mortality, etc.) over time. This means we can look at really interesting outputs like vegetation structure (sizes and PFTs of trees), canopy closure, and successional dynamics.\n", + "\n", + "We have already run a simulation at a NEON site: Bartlett Experimental Forest (**BART**). BART is a terrestrial NEON field site located within the Saco River Valley of the White Mountain National Forest in Carroll County, New Hampshire. Climate and glacial till soil make BART an ideal area for old-growth northern deciduous hardwoods consisting of beech (*Fagus sp.*) and sugar maple (*Acer saccarum*). White pines (*Pinus strobus*) are dispersed throughout the site but are primarily found in lower elevations. Softwood trees such as hemlock (*Tsuga canadensis*), balsam fir (*Abies balsamea*), and spruce are frequently found on cool steep slopes or in lower elevations with poor drainage.\n", + "\n", + "For more information about BART, see the NEON webpage for the site. \n", + "\n", + "For this simulation, we ran FATES in fixed biogeography mode, but ran two simulations:\n", + "\n", + "1) with FATES initialized with inventory data from NEON and run for 10 years\n", + "2) with FATES run from bare ground and run for 100 years\n", + "\n", + "For information on how to do inventory initialization with FATES, see the online documentation." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "*As before, start by loading some packages*" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import xarray as xr\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import matplotlib.colors as mcolors\n", + "import matplotlib.cm as cm\n", + "from matplotlib.ticker import FuncFormatter\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Helper functions\n", + "\n", + "These are a few helper functions that we will use in this notebook." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def annual_sum(raw_values: xr.DataArray, conversion_factor: float = 1.0) -> xr.DataArray:\n", + " \"\"\"Computes annual sum\n", + "\n", + " Args:\n", + " raw_values (xr.DataArray): input raw data\n", + " conversion_factor (float, optional): conversion factor. Defaults to 1.0\n", + "\n", + " Returns:\n", + " xr.DataArray: annual sum output\n", + " \"\"\"\n", + "\n", + " months = raw_values[\"time.daysinmonth\"]\n", + " return conversion_factor * (months * raw_values).groupby(\"time.year\").sum()\n", + "\n", + "\n", + "def annual_mean(raw_values: xr.DataArray, conversion_factor: float = 1.0) -> xr.DataArray:\n", + " \"\"\"Computes weighted annual mean using daysinmonth for missing-aware inputs.\n", + "\n", + " Args:\n", + " raw_values (xr.DataArray): input raw data\n", + " conversion_factor (float, optional): conversion factor. Defaults to 1.0\n", + "\n", + " Returns:\n", + " xr.DataArray: annual mean output\n", + " \"\"\"\n", + "\n", + " months = raw_values[\"time.daysinmonth\"]\n", + " \n", + " # multiply by number of days in month and conversion factor\n", + " weighted = (raw_values * conversion_factor) * months\n", + "\n", + " # compute number of valid days per year\n", + " valid_days = months.where(raw_values.notnull())\n", + "\n", + " # group and sum weighted data and valid days\n", + " ann_sum = weighted.groupby(\"time.year\").sum(dim=\"time\", skipna=True)\n", + " days_per_year = valid_days.groupby(\"time.year\").sum(dim=\"time\", skipna=True)\n", + "\n", + " return ann_sum / days_per_year.where(days_per_year > 0.0)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---\n", + "## 1. Reading and formatting data\n", + "\n", + "**Note**: all exercises below make use of pre-run simulations located in `/glade/campaign/cesm/tutorial/diagnostics_tutorial_archive/`." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 1.1. Reading in the data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "prestage_dir = \"/glade/campaign/cesm/tutorial/diagnostics_tutorial_archive/fates.BART\"\n", + "inv_file = os.path.join(prestage_dir, \"i.fates.BART.inv.nc\")\n", + "bg_file = os.path.join(prestage_dir, \"i.fates.BART.bareground.nc\")\n", + "\n", + "ds_inv = xr.open_dataset(inv_file).isel(lndgrid=0)\n", + "ds_bg = xr.open_dataset(bg_file).isel(lndgrid=0)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Because this is a single site, and not a regional or global run, we can just grab the one gridcell with `.sel(lndgrid=0)`, and we additionally don't need to multiply any of the FATES variables by `FATES_FRACTION`, because we just assume the gridcell is completely naturally vegetated." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# print information about the inventory-initialized dataset\n", + "ds_inv" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 1.2 FATES multi-plexed dimensions\n", + "\n", + "Before we plot the simulations, we will go over a FATES-specific feature of its history output.\n", + "\n", + "A single FATES variable can be resolved across several categorical axes at once (e.g., size class, plant functional type (PFT), canopy layer, patch age, etc.) FATES writes these out by multiplexing: when it writes the file, it flattens a pair (or more) of these category axes into one packed dimension.\n", + "\n", + "You can read which axes are packed into a variable using its name. One of the variables we'll start with, `FATES_NPLANT_SZPF` (number of plants), carries the `SZPF` suffix. It has been flattened across size class (`SZ`) and PFT (`PF`) onto a single dimension, `fates_levscpf`, whose length is `n_sizeclass * n_pft`. So instead of a tidy `(time, fates_levscls, fates_levpft)` array, what's actually on the file is `(time, fates_levscpf)`, with the two categories interleaved along one axis. This is not really useful in its raw form because we can't select, say, \"10–20 cm trees\" or \"broadleaf deciduous trees\" until we split the multiplexed axis back into its two actual dimensions.\n", + "\n", + "The helper functions below uses information about each individual dimension to unstack `fates_levscpf` back into separate size-class and PFT dimensions." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def deduplex(ds: xr.Dataset, da: xr.DataArray, dim1_short: str, dim2_short: str, preserve_order: bool = True) -> xr.DataArray:\n", + " \"\"\"Reshape a duplexed FATES dimension into its constituent dimensions\n", + "\n", + " Args:\n", + " ds (xarray Dataset): Dataset containing the variable with dimension to de-duplex\n", + " da (xr.DataArray): (Name of) variable with dimension to de-duplex\n", + " dim1_short (str): Short name of first duplexed dimension. E.g., when de-duplexing\n", + " fates_levagepft, dim1_short=age.\n", + " dim2_short (str): Short name of second duplexed dimension. E.g., when de-duplexing\n", + " fates_levagepft, dim2_short=pft.\n", + " preserve_order (bool, optional): Preserve order of dimensions of input DataArray? Defaults\n", + " to True. Might be faster if False.\n", + "\n", + " Returns:\n", + " xr.DataArray: the deduplexed output\n", + " \"\"\"\n", + " if dim1_short == dim2_short:\n", + " raise ValueError(\"deduplex() can't currently handle dim1_short==dim2_short\")\n", + "\n", + " dim_combined = _get_dim_combined(dim1_short, dim2_short)\n", + " if dim_combined not in da.dims:\n", + " raise NameError(\n", + " f\"Dimension {dim_combined} not present in DataArray with dims {da.dims}\"\n", + " )\n", + "\n", + " dim1 = _get_check_dim(dim1_short, ds)\n", + " dim2 = _get_check_dim(dim2_short, ds)\n", + "\n", + " n_dim1 = len(ds[dim1])\n", + " da_out = (\n", + " da.rolling({dim_combined: n_dim1}, center=False)\n", + " .construct(dim1)\n", + " .isel({dim_combined: slice(n_dim1 - 1, None, n_dim1)})\n", + " .rename({dim_combined: dim2})\n", + " .assign_coords({dim1: ds[dim1]})\n", + " .assign_coords({dim2: ds[dim2]})\n", + " )\n", + "\n", + " if preserve_order:\n", + " new_dim_order = []\n", + " for dim in da_out.dims:\n", + " if dim == dim2:\n", + " new_dim_order.append(dim1)\n", + " if dim != dim1:\n", + " new_dim_order.append(dim)\n", + " da_out = da_out.transpose(*new_dim_order)\n", + " return da_out\n", + "\n", + "\n", + "def _get_dim_combined(dim1_short: str, dim2_short) -> str:\n", + " \"\"\"Get duplexed dimension name, given two short names\n", + "\n", + " Args:\n", + " dim1_short (str): Short name of first duplexed dimension. E.g., when de-duplexing\n", + " fates_levscpf, dim1_short=scls.\n", + " dim2_short (str): Short name of second duplexed dimension. E.g., when de-duplexing\n", + " fates_levscpf, dim2_short=pft.\n", + "\n", + " Returns:\n", + " str: duplexed dimension name\n", + " \"\"\"\n", + " dim_combined = \"fates_lev\" + dim1_short + dim2_short\n", + "\n", + " # handle further-shortened dim names\n", + " if dim_combined == \"fates_levcanleaf\":\n", + " dim_combined = \"fates_levcnlf\"\n", + " elif dim_combined == \"fates_levcanpft\":\n", + " dim_combined = \"fates_levcapf\"\n", + " elif dim_combined == \"fates_levcdamscls\":\n", + " dim_combined = \"fates_levcdsc\"\n", + " elif dim_combined == \"fates_levsclsage\":\n", + " dim_combined = \"fates_levscag\"\n", + " elif dim_combined == \"fates_levsclspft\":\n", + " dim_combined = \"fates_levscpf\"\n", + "\n", + " return dim_combined\n", + "\n", + "\n", + "def _get_check_dim(dim_short: str, dataset: xr.Dataset) -> str:\n", + " \"\"\"Get dim name from short code and ensure it's on Dataset\n", + "\n", + " Probably only useful internally to this module; see deduplex().\n", + "\n", + " Args:\n", + " dim_short (str): The short name of the dimension. E.g., \"age\"\n", + " dataset (xr.Dataset): The Dataset we expect to include the dimension\n", + "\n", + " Raises:\n", + " NameError: Dimension not found on Dataset\n", + "\n", + " Returns:\n", + " str: The long name of the dimension. E.g., \"fates_levage\"\n", + " \"\"\"\n", + "\n", + " dim = \"fates_lev\" + dim_short\n", + " if dim not in dataset.dims:\n", + " raise NameError(\n", + " f\"Dimension {dim} not present in Dataset with dims {dataset.dims}\"\n", + " )\n", + " return dim" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---\n", + "## 2. What's growing at BART? PFT composition\n", + "\n", + "First, we will investigate BART's PFT composition. Remember that we ran these simulations in \"fixed biogeography mode.\"\n", + "\n", + "*Remember that fixed biogeography mode means that we allow PFTs to compete with each other, but restrict **which** PFTs can grow at each site*.\n", + "\n", + "For our runs, we only allowed two PFTs to grow: broadleaf deciduous (e.g., beech and sugar maple) and needleleaf evergreen (e.g., white pine, hemlock, fir, and spruce).\n", + "\n", + "Quantities resolved by PFT alone are written to history with a `_PF` suffix on the `fates_levpft` dimension, so they come out as a tidy `(time, pft)` array we can plot directly. We'll look at vegetation carbon (`FATES_VEGC_PF`) as our composition metric." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 2.1 Inventory composition" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ds_inv[\"FATES_VEGC_PF\"]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# these are the PFTs at this site\n", + "active_pfts = {\n", + " 2: (\"needleleaf evergreen\", \"darkslategray\"),\n", + " 6: (\"broadleaf deciduous\", \"forestgreen\"),\n", + "}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "vegc_inv = annual_mean(ds_inv[\"FATES_VEGC_PF\"], 10.0) # tC/ha\n", + "\n", + "fig, ax = plt.subplots(figsize=(8, 4.5))\n", + "ax.stackplot(vegc_inv.year,\n", + " *[vegc_inv.sel(fates_levpft=pft).values for pft in active_pfts],\n", + " labels=[v[0] for v in active_pfts.values()],\n", + " colors=[v[1] for v in active_pfts.values()],\n", + " alpha=0.85)\n", + "ax.set(xlabel=\"Simulation Year\", ylabel=\"Biomass (tC ha$^{-1}$)\",\n", + " title=\"FATES PFT composition: initialized run\")\n", + "ax.legend(loc=\"upper left\")\n", + "fig.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "
\n", + "Click here for the solution
\n", + "\n", + "![plot example](../../../images/diagnostics/clm_ctsm/fates_bart_vegc_composition_init.png)\n", + "\n", + "*

Figure: Biomass for our inventory-initialized case.

*\n", + " \n", + "
\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Question: \n", + "\n", + " - What do you notice about this plot? Which PFT is dominating? Does that match with the site description above?" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 2.2 Comparison to bare-ground run\n", + "\n", + "With the above plot, we were sort of \"cheating\" at getting the composition right because we initialize the model with actual inventory data. Let's also look at our bareground run and see how that compares:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "vegc_bg = annual_mean(ds_bg[\"FATES_VEGC_PF\"], 10.0) # tonnes C/ha\n", + "\n", + "fig, ax = plt.subplots(figsize=(8, 4.5))\n", + "ax.stackplot(np.arange(len(vegc_bg.year)),\n", + " *[vegc_bg.sel(fates_levpft=pft).values for pft in active_pfts],\n", + " labels=[v[0] for v in active_pfts.values()],\n", + " colors=[v[1] for v in active_pfts.values()],\n", + " alpha=0.85)\n", + "ax.set(xlabel=\"Stand Age\", ylabel=\"Biomass (tC ha$^{-1}$)\",\n", + " title=\"FATES PFT composition: bareground run\")\n", + "ax.legend(loc=\"upper left\")\n", + "fig.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "
\n", + "Click here for the solution
\n", + "\n", + "![plot example](../../../images/diagnostics/clm_ctsm/fates_bart_vegc_composition_bg.png)\n", + "\n", + "*

Figure: Biomass over time for the bare-ground case.

*\n", + " \n", + "
\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Question:\n", + " - What do you notice about the bareground run? How is it different from the inventory-initialized run?" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---\n", + "## 3. Anatomy of a FATES forest\n", + "\n", + "### 3.1. Size structure\n", + "\n", + "So the biomass plots above show what is growing and how much carbon is there, but remember that FATES keeps track of individual cohorts of specific size, age, and PFT. This architecture is what distinguishes FATES from a more \"big leaf\" structure like the default CLM vegetation model.\n", + "\n", + "So, let's look at size structure of our initialized stand, using the multiplexed variable `FATES_NPLANT_SZPF`.\n", + "\n", + "First we will \"deduplex\" it and then sum across the `fates_levpft` dimension to look at total size structure over time." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "nplant = deduplex(ds_inv, ds_inv[\"FATES_NPLANT_SZPF\"], 'scls', 'pft')\n", + "total = annual_mean(nplant.sum(\"fates_levpft\"), 1e4) # plants/ha\n", + "total = total.sel(fates_levscls=slice(1, None)) # not really interested in 0-1 cms" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "edges = total[\"fates_levscls\"].values\n", + "y_edges = np.append(edges, edges[-1] + (edges[-1] - edges[-2]))\n", + "years = total.year\n", + "x_edges = np.append(years, years[-1] + 1)\n", + "\n", + "vals = total.transpose(\"fates_levscls\", \"year\").values\n", + "vmin = np.nanmin(vals[vals > 0])\n", + "\n", + "fig, ax = plt.subplots(figsize=(9, 5))\n", + "cmap = plt.cm.viridis.copy()\n", + "cmap.set_bad(\"0.92\")\n", + "mesh = ax.pcolormesh(\n", + " x_edges, y_edges, np.ma.masked_less_equal(vals, 0),\n", + " norm=mcolors.LogNorm(vmin=vmin, vmax=np.nanmax(vals)),\n", + " cmap=cmap, shading=\"flat\",\n", + ")\n", + "ax.set(xlabel=\"Simulation Year\", ylabel=\"DBH Size Class (cm)\",\n", + " title=\"FATES stand size structure\")\n", + "ax.set_yticks(edges)\n", + "fig.colorbar(mesh, ax=ax, label=\"Stem Density (plants ha$^{-1}$)\")\n", + "fig.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "
\n", + "Click here for the solution
\n", + "\n", + "![plot example](../../../images/diagnostics/clm_ctsm/fates_bart_stand_structure.png)\n", + "\n", + "*

Figure: Stem density across DBH size classes (vertical) over time (horizontal), on a log color scale.

*\n", + " \n", + "
\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Questions:\n", + " - What do you notice about the size structure at this site? \n", + " - Why do we need to use a log scale to look at this variable?" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 3.2. Canopy structure\n", + "\n", + "The heatmap shows how many trees of each size there are, but it doesn't show where those trees sit in the light environment. FATES sorts every cohort into a canopy layer (crowns in full sun) or one of several understory layers (layers shaded beneath), following the Perfect Plasticity Approximation (see more [here](https://fates-users-guide.readthedocs.io/projects/tech-doc/en/stable/fates_tech_note.html#canopy-structure-and-the-perfect-plasticity-approximation)).\n", + "\n", + "Splitting the same size distribution by layer turns the heatmap into a canopy pyramid.\n", + "\n", + "We will first grab the variable `FATES_NPLANT_CANOPY_SZPF`, which is the number of plants in the top-most canopy layer, indexed by size and pft, as well as `FATES_NPLANT_USTORY_SZPF`, which is the same variable but for the understory layers.\n", + "\n", + "As before, we will \"deduplex\" these variables, sum by PFT, and take the annual mean, then remove the smallest size class.\n", + "\n", + "We will then plot the first year." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# canopy plants\n", + "canopy = annual_mean(deduplex(ds_inv, ds_inv[\"FATES_NPLANT_CANOPY_SZPF\"], 'scls', 'pft').sum(\"fates_levpft\"), 1e4)\n", + "canopy = canopy.sel(fates_levscls=slice(1, None))\n", + "\n", + "# understory plants\n", + "ustory = annual_mean(deduplex(ds_inv, ds_inv[\"FATES_NPLANT_USTORY_SZPF\"], 'scls', 'pft').sum(\"fates_levpft\"), 1e4)\n", + "ustory = ustory.sel(fates_levscls=slice(1, None))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# helper function for drawing a canopy/size class pyramid at a specific year\n", + "def draw_pyramid(ax, ustory, canopy, scls_labels, yr):\n", + " \"\"\"Mirrored size pyramid at year index yr: understory left, canopy right.\"\"\"\n", + " ax.barh(y, -ustory.isel(year=yr).values, height=0.85,\n", + " color=LAYER_COLORS[\"understory\"], label=\"understory\")\n", + " ax.barh(y, canopy.isel(year=yr).values, height=0.85,\n", + " color=LAYER_COLORS[\"canopy\"], label=\"canopy\")\n", + " ax.axvline(0, color=\"0.3\", lw=0.8)\n", + " ax.set_yticks(y)\n", + " ax.set_yticklabels(scls_labels, fontsize=8)\n", + "\n", + "\n", + "edges = ustory[\"fates_levscls\"].values\n", + "scls_labels = [f\"{int(edges[i])}–{int(edges[i+1])}\" if i+1 < len(edges)\n", + " else f\"{int(edges[i])}+\" for i in range(len(edges))]\n", + "\n", + "LAYER_COLORS = {\"understory\": \"#9CA86E\", \"canopy\": \"#2E6E5E\"}\n", + "y = np.arange(len(edges))\n", + "xmax = max(ustory.max(), canopy.max())\n", + "\n", + "fig, ax = plt.subplots(figsize=(7, 5))\n", + "draw_pyramid(ax, ustory, canopy, scls_labels, 0) # plotting .isel(year=0)\n", + "ax.set_xscale(\"symlog\", linthresh=10)\n", + "ax.set_xlim(-xmax * 1.05, xmax * 1.05)\n", + "ticks = ax.xaxis.get_major_locator().tick_values(0, xmax)\n", + "ax.set_xticks(np.concatenate([-ticks, ticks]))\n", + "ax.xaxis.set_major_formatter(FuncFormatter(lambda x, _: f\"{abs(x):g}\"))\n", + "ax.set(xlabel=\"Stem Density (plants ha$^{-1}$)\", ylabel=\"DBH size class (cm)\",\n", + " title=\"Canopy vs. understory size structure\")\n", + "ax.legend(loc=\"upper left\")\n", + "fig.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "
\n", + "Click here for the solution
\n", + "\n", + "![plot example](../../../images/diagnostics/clm_ctsm/fates_bart_canopy_structure.png)\n", + "\n", + "*

Figure: The size distribution split by canopy layer: understory stems to the left, canopy to the right.

*\n", + " \n", + "
\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Questions:\n", + " - How do the distributions differ between the sizes of the canopy and understory trees?\n", + " - Why are there some really small canopy trees?" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 3.3 Layered LAI\n", + "\n", + "The same canopy vs. understory split can also be seen in the forest's leaf area. We can separate LAI by its layer just as we did with stems.\n", + "\n", + "For this plot we will look at the `FATES_LAI_CANOPY_SZ` and `FATES_LAI_USTORY_SZ` variables, which are LAI indexed by size class for the canopy (top-most) and understory (lower) layers. \n", + "\n", + "We don't really care about the size class, so we will just sum across that index by doing `.sum(dim='fates_levscls')`.\n", + "\n", + "To look at peak LAI for the year we will use the `.groupby` method and take the maximum value:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "lai_can = ds_inv.FATES_LAI_CANOPY_SZ.sum(dim='fates_levscls')\n", + "lai_ust = ds_inv.FATES_LAI_USTORY_SZ.sum(dim='fates_levscls')\n", + "\n", + "peak_can = lai_can.groupby(\"time.year\").max()\n", + "peak_ust = lai_ust.groupby(\"time.year\").max()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Then we will just grab the first year's data to plot:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "peak_ust_year = peak_ust.values[0]\n", + "peak_can_year = peak_can.values[0]\n", + "year = [str(years[0].values)]\n", + "\n", + "fig, ax = plt.subplots(figsize=(2.5, 4.5))\n", + "ax.bar(0, peak_ust_year, color=LAYER_COLORS[\"understory\"], label=\"understory\")\n", + "ax.bar(0, peak_can_year, bottom=peak_ust_year,\n", + " color=LAYER_COLORS[\"canopy\"], label=\"canopy\")\n", + "ax.set_xticks([0])\n", + "ax.set_xticklabels(year)\n", + "ax.set_ylabel(\"Annual Maximum LAI (m$^2$ m$^{-2}$)\")\n", + "ax.set_title(\"Layered LAI\")\n", + "\n", + "handles, labels = ax.get_legend_handles_labels()\n", + "ax.legend(handles[::-1], labels[::-1], loc=\"center left\",\n", + " bbox_to_anchor=(1.02, 0.5));" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "
\n", + "Click here for the solution
\n", + "\n", + "![plot example](../../../images/diagnostics/clm_ctsm/fates_bart_layered_LAI_init.png)\n", + "\n", + "*

Figure: Leaf area index partitioned into canopy and understory layers.

*\n", + " \n", + "
\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### LAI over time\n", + "\n", + "One of the more interesting things about FATES is that you can see successional dynamics over time. In a forest landscape, we should see a difference between early-, mid-, and late-successional stands in terms of composition and structure. Let's take a look at how LAI and the canopy structure evolves for our bare-ground case.\n", + "\n", + "\n", + "This time we will plot monthly LAI over time (rather than peak annual) so we can take a look at the sub-annual dynamics." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "tyear = ds_bg[\"time\"].dt.year + (ds_bg[\"time\"].dt.month - 1) / 12\n", + "stand_age = tyear - tyear.min()\n", + "\n", + "lai_can_bg = ds_bg.FATES_LAI_CANOPY_SZ.sum(dim='fates_levscls')\n", + "lai_ust_bg = ds_bg.FATES_LAI_USTORY_SZ.sum(dim='fates_levscls')\n", + "\n", + "fig, axes = plt.subplots(1, 3, figsize=(14, 4.3), sharey=True)\n", + "windows = [(0, 9, \"Early stand\"), (20, 29, \"Mid-successional\"), (90, 99, \"Old growth\")]\n", + "\n", + "for ax, (y0, y1, title) in zip(axes, windows):\n", + " ax.stackplot(stand_age.values, lai_ust_bg.values, lai_can_bg.values,\n", + " labels=[\"understory\", \"canopy\"],\n", + " colors=[LAYER_COLORS[\"understory\"], LAYER_COLORS[\"canopy\"]], alpha=0.9)\n", + " ax.set_xlim(y0, y1)\n", + " ax.set_title(title, fontsize=10)\n", + " ax.set_xlabel(\"Stand Age (years)\")\n", + " ax.margins(x=0)\n", + "\n", + "axes[0].set_ylabel(\"LAI (m$^2$ m$^{-2}$)\")\n", + "\n", + "handles, labels = ax.get_legend_handles_labels()\n", + "ax.legend(handles[::-1], labels[::-1], loc=\"center left\",\n", + " bbox_to_anchor=(1.02, 0.5))\n", + "\n", + "fig.suptitle(\"Layered LAI over time\", y=1.03)\n", + "fig.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "
\n", + "Click here for the solution
\n", + "\n", + "![plot example](../../../images/diagnostics/clm_ctsm/fates_bart_layered_LAI_bg.png)\n", + "\n", + "*

Figure: Monthly leaf area index partitioned into canopy and understory layers over time at three different snapshots: early-, mid-, and late-successional.

*\n", + " \n", + "
\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Questions:\n", + "- What do you notice happening over time at this site?\n", + "- What are the major differences between the early-, mid-, and late-successional stands?\n", + "- How does the bare-ground case differ from our single-year plot using the inventory-initialized stand?" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---\n", + "## 4. Demographic Processes\n", + "\n", + "### 4.1. Canopy Layer Promotion and Demotion\n", + "\n", + "The processes that drive the LAI-over-time plot above include phenology, growth, mortality (plant death from a variety of factors), and recruitment (new trees establishing). As these play out, FATES continually re-sorts trees between canopy layers (promotion and demotion). Together these are the core of how FATES simulates succession, competition, and vegetation dynamics.\n", + "\n", + "To compare to our LAI plot above, let's first look at two useful variables: `FATES_PROMOTION_RATE_SZ` and `FATES_DEMOTION_RATE_SZ`. These are the *fluxes* (in plants m$^{-2}$ yr$^{-1}$) from the a lower to higher canopy layer (`PROMOTION`) and from a higher to a lower canopy layer (`DEMOTION`).\n", + "\n", + "Under PPA, trees fill the canopy from the top down by height, occupying space until total crown area equals the patch area (i.e., canopy closure). Beyond that point there's no room left in the layer, so the lowest-ranked trees are demoted into a lower, shaded understory layer. Conversely, when a tree dies and opens a gap, its layer is no longer full, and the tallest trees from below are promoted up to fill it. \n", + "\n", + "Let's look at these two variables, focusing again on those three windows of time in our bareground run." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "promo = annual_mean(ds_bg[\"FATES_PROMOTION_RATE_SZ\"].sum(\"fates_levscls\"), 1e4)\n", + "demo = annual_mean(ds_bg[\"FATES_DEMOTION_RATE_SZ\"].sum(\"fates_levscls\"), 1e4)\n", + "stand_age = promo.year - promo.year.min()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "FLUX_COLORS = {\"promotion\": \"#3A7CA5\",\n", + " \"demotion\": \"#D98A3D\"}\n", + "\n", + "fig, axes = plt.subplots(1, 3, figsize=(14, 4.3), sharey=True)\n", + "windows = [(0, 9, \"Early stand\"), (20, 29, \"Mid-successional\"), (90, 99, \"Old growth\")]\n", + "\n", + "for ax, (y0, y1, title) in zip(axes, windows):\n", + " ax.fill_between(stand_age.values, 0, promo.values,\n", + " color=FLUX_COLORS[\"promotion\"], alpha=0.85,\n", + " label=\"promotion\")\n", + " ax.fill_between(stand_age.values, 0, -demo.values,\n", + " color=FLUX_COLORS[\"demotion\"], alpha=0.85,\n", + " label=\"demotion\")\n", + " ax.axhline(0, color=\"0.3\", lw=1)\n", + " ax.set_xlim(y0, y1)\n", + " ax.set_title(title, fontsize=10)\n", + " ax.set_xlabel(\"Stand Age (years)\")\n", + " ax.margins(x=0)\n", + " ax.yaxis.set_major_formatter(FuncFormatter(lambda v, _: f\"{abs(v):g}\"))\n", + "\n", + "axes[0].set_ylabel(\"Layer Flux (plants ha$^{-1}$ yr$^{-1}$)\")\n", + "handles, labels = ax.get_legend_handles_labels()\n", + "ax.legend(handles[::-1], labels[::-1], loc=\"center left\",\n", + " bbox_to_anchor=(1.02, 0.5))\n", + "fig.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "
\n", + "Click here for the solution
\n", + "\n", + "![plot example](../../../images/diagnostics/clm_ctsm/fates_bart_promotion_demotion.png)\n", + "\n", + "*

Figure: Fluxes between canopy layers at three different snapshots: early-, mid-, and late-successional.

*\n", + " \n", + "
\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Questions:\n", + "- What do you notice happening over time with the promotion and demotion fluxes?\n", + "- Looking side-by-side with LAI over time, can you explain what is going on?" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 4.2. Growth and Mortality\n", + "\n", + "#### Diameter Growth\n", + "\n", + "So, as stated, the drivers of these promotion and demotion fluxes are growth and mortality of individual cohorts. We can investigate growth over time using the `FATES_DDBH_SZPF` variable, which shows the diameter increment growth of cohorts, weighted by the number of plants. We can un-normalize it using `FATES_NPLANT_SZPF`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ddbh = deduplex(ds_bg, ds_bg[\"FATES_DDBH_SZPF\"], 'scls', 'pft').sum(\"fates_levpft\")*100.0 # cm/m2/yr\n", + "nplant = deduplex(ds_bg, ds_bg[\"FATES_NPLANT_SZPF\"], 'scls', 'pft').sum(\"fates_levpft\")\n", + "\n", + "# un-normalize and take annual average\n", + "dbh_growth = annual_mean(xr.where(nplant > 0.0, ddbh / nplant, np.nan)).sel(fates_levscls=slice(1, None))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Mortality\n", + "\n", + "We can also investigate mortality in FATES. FATES outputs mortality by size class (plants m$^{-2}$ yr$^{-1}$) for all of the different causes of mortality:\n", + "\n", + "- **background**: background, \"random\" mortality; should be a static rate (`FATES_MORTALITY_BACKGROUND_SZ`)\n", + "- **freezing**: mortality due to cold intolerance (`FATES_MORTALITY_FREEZING_SZ`)\n", + "- **carbon starvation**: mortality due to low growth (`FATES_MORTALITY_CSTARV_SZ`)\n", + "- **hydraulic failure**: mortality from moisture stress and/or cavitation (`FATES_MORTALITY_HYDRAULIC_SZ`)\n", + "- **impact**: mortality from getting hit by a falling tree (`FATES_MORTALITY_IMPACT_SZ`)\n", + "- **termination**: forced mortality from when a cohort gets too small in number density (`FATES_MORTALITY_TERMINATION_SZ`)\n", + "\n", + "There is also a `FATES_MORTALITY_SENESCENCE_SZ` output, which is mortality due to age- or size-related \"senescence\". But for these runs we have turned this feature off." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "MORT_VARS = {\n", + " \"background\": \"FATES_MORTALITY_BACKGROUND_SZ\",\n", + " \"freezing\": \"FATES_MORTALITY_FREEZING_SZ\",\n", + " \"carbon starvation\": \"FATES_MORTALITY_CSTARV_SZ\",\n", + " \"hydraulic failure\": \"FATES_MORTALITY_HYDRAULIC_SZ\",\n", + " \"impact\": \"FATES_MORTALITY_IMPACT_SZ\",\n", + " \"termination\": \"FATES_MORTALITY_TERMINATION_SZ\",\n", + "}\n", + "mort = sum(annual_mean(ds_bg[v], 1e4) for v in MORT_VARS.values()).sel(fates_levscls=slice(1, None))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# helper function for plotting heatmap of size over time\n", + "def size_time_heat(ax, da, cmap, label, title, vmin=None, vmax=None, norm=None):\n", + "\n", + " stand_age = da.year - da.year.min()\n", + " e = da[\"fates_levscls\"].values\n", + " y_edges = np.append(e, e[-1] + (e[-1] - e[-2]))\n", + " x_edges = np.append(stand_age.values, stand_age.values[-1] + 1)\n", + " vals = da.transpose(\"fates_levscls\", \"year\").values\n", + " cm = plt.get_cmap(cmap).copy()\n", + " cm.set_bad(\"0.92\")\n", + " m = ax.pcolormesh(x_edges, y_edges, np.ma.masked_invalid(vals),\n", + " cmap=cm, shading=\"flat\", vmin=vmin, vmax=vmax,\n", + " norm=norm)\n", + " ax.set(xlabel=\"Stand Age\", title=title)\n", + " ax.set_yticks(e)\n", + " plt.colorbar(m, ax=ax, label=label)\n", + "\n", + "\n", + "# just grab DBH size classes that are occupied\n", + "keep = dbh_growth.notnull().any(\"year\")\n", + "dbh_growth_keep = dbh_growth.isel(fates_levscls=keep)\n", + "mort_keep = mort.isel(fates_levscls=keep)\n", + "\n", + "# plot growth and mortality side by side\n", + "fig, (a0, a1) = plt.subplots(1, 2, figsize=(13, 5), sharey=True)\n", + "size_time_heat(a0, dbh_growth_keep, \"magma\",\n", + " \"Diameter Increment (cm yr$^{-1}$)\", \"Growth\", vmin=0.0, vmax=0.5)\n", + "mort_keep = mort_keep.where(mort_keep > 0)\n", + "vals = mort_keep.values\n", + "pos = vals[np.isfinite(vals) & (vals > 0)]\n", + "norm = mcolors.LogNorm(vmin=vmin, vmax=np.nanmax(vals))\n", + "size_time_heat(a1, mort_keep, \"magma\", \"Mortality (plants ha$^{-1}$ yr$^{-1}$)\",\n", + " \"Mortality\", norm=norm)\n", + "a0.set_ylabel(\"DBH Size Class (cm)\")\n", + "fig.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "
\n", + "Click here for the solution
\n", + "\n", + "![plot example](../../../images/diagnostics/clm_ctsm/fates_bart_growth_mortality.png)\n", + "\n", + "*

Figure: Growth (left) and mortality (right) by size class over time.

*\n", + " \n", + "
\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Questions:\n", + "- How does diameter growth change with size class?\n", + "- Where (which size classes) and when do the largest diameter growths occur?\n", + "- How does mortality change over time and with size class?" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Mortality by cause\n", + "\n", + "We can also look at the causes of mortality in FATES:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "mort_by_cause = {label: ds_bg[v].sel(fates_levscls=slice(1, None)).sum(\"fates_levscls\")\n", + " for label, v in MORT_VARS.items()}\n", + "mort_annual = {k: annual_mean(v, 1e4) for k, v in mort_by_cause.items()}\n", + "years = next(iter(mort_annual.values()))[\"year\"].values\n", + "stand_age = years - years.min()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "CAUSE_COLORS = {\n", + " \"background\": \"#7f7f7f\",\n", + " \"freezing\": \"#4C72B0\",\n", + " \"carbon starvation\": \"#DD8452\",\n", + " \"hydraulic failure\": \"#C44E52\",\n", + " \"impact\": \"#8172B3\",\n", + " \"termination\": \"#CCB974\",\n", + "}\n", + "order = [k for k in CAUSE_COLORS if k in mort_annual]\n", + "\n", + "fig, ax = plt.subplots(figsize=(10, 5))\n", + "ax.stackplot(stand_age, *[mort_annual[k].values for k in order], labels=order,\n", + " colors=[CAUSE_COLORS[k] for k in order], alpha=0.9)\n", + "ax.set(xlabel=\"Stand Age\", ylabel=\"Mortality (plants ha$^{-1}$ yr$^{-1}$)\",\n", + " title=\"Mortality partitioned by cause\")\n", + "ax.legend(loc=\"upper right\", fontsize=8, ncol=2)\n", + "ax.margins(x=0)\n", + "fig.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "
\n", + "Click here for the solution
\n", + "\n", + "![plot example](../../../images/diagnostics/clm_ctsm/fates_bart_mortality.png)\n", + "\n", + "*

Figure: Mortality over time by cause.

*\n", + " \n", + "
\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Questions:\n", + "- What is the dominant cause of mortality at BART?\n", + "- Does this change over time?" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Snapshots of Mortality\n", + "\n", + "Let's look at those snapshots again, but instead of LAI or promotion/demotion fluxes, lets look at the same plot as above:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, axes = plt.subplots(1, 3, figsize=(14, 4.3))\n", + "windows = [(0, 9, \"Early stand\"), (20, 29, \"Mid-successional\"), (90, 99, \"Old growth\")]\n", + "\n", + "for ax, (y0, y1, title) in zip(axes, windows):\n", + " ax.stackplot(stand_age, *[mort_annual[k].values for k in order],\n", + " labels=order, colors=[CAUSE_COLORS[k] for k in order],\n", + " alpha=0.9)\n", + " ax.set_xlim(y0, y1)\n", + " ax.set_title(title, fontsize=10)\n", + " ax.set_xlabel(\"Stand Age (years)\")\n", + " ax.margins(x=0)\n", + "axes[0].set_ylabel(\"Mortality (plants ha$^{-1}$ yr$^{-1}$)\")\n", + "handles, labels = ax.get_legend_handles_labels()\n", + "ax.legend(handles[::-1], labels[::-1], loc=\"center left\",\n", + " bbox_to_anchor=(1.02, 0.5))\n", + "fig.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "
\n", + "Click here for the solution
\n", + "\n", + "![plot example](../../../images/diagnostics/clm_ctsm/fates_bart_mortality_snapshots.png)\n", + "\n", + "*

Figure: Mortality over time by cause: snapshots in time.

*\n", + " \n", + "
\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Questions:\n", + "- What is happening in each snapshot?\n", + "- Comparing to the LAI and promotion/demotion fluxes figures, can you formulate a story for what is going on?" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---\n", + "## 5. A forest landscape through succession\n", + "\n", + "So far we've looked at the stand as a whole. In this final section we will step back in two ways: first to see how the causes of mortality differentially impact cohorts of different sizes, and as a stand ages, and then to the landscape scale, where FATES represents the site not as one forest but as a mosaic of patches of different ages.\n", + "\n", + "### 5.1. Mortality causes over time\n", + "\n", + "Interestingly, FATES can also show us the causes of mortality across time *and* for different size structures. Here we plot three snapshots again (early-, mid-, and late-successional) and plot the fraction of mortality for each cause by size class." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "mort_by_size = {label: ds_bg[v].sel(fates_levscls=slice(1, None))\n", + " for label, v in MORT_VARS.items()}\n", + "mort_annual_bysize = {k: annual_mean(v, 1e4) for k, v in mort_by_size.items()}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# helper function for plotting mortality by size class\n", + "def window_stack(y1, y2, mort_annual):\n", + " y0 = int(mort_annual[order[0]][\"year\"].values[0])\n", + " M = np.vstack([mort_annual[k].sel(year=slice(y0+y1, y0+y2)).mean(\"year\").values for k in order])\n", + " tot = M.sum(axis=0)\n", + " frac = np.divide(M, tot, out=np.full_like(M, np.nan), where=tot > 0)\n", + " return frac\n", + "\n", + "\n", + "windows = {\"Early stand\": (0, 9), \"Mid-successional\": (20, 29), \"Old growth\": (90, 99)}\n", + "\n", + "order = [k for k in CAUSE_COLORS if k in mort_by_size]\n", + "edges = mort_by_size[order[0]][\"fates_levscls\"].values\n", + "labels = [f\"{int(edges[i])} – {int(edges[i+1])}\" if i+1 < len(edges)\n", + " else f\"{int(edges[i])}+\" for i in range(len(edges))]\n", + "\n", + "fig, axes = plt.subplots(1, 3, figsize=(14, 6), sharey=True)\n", + "for ax, (name, (y1, y2)) in zip(axes, windows.items()):\n", + " frac = window_stack(y1, y2, mort_annual_bysize)\n", + " left = np.zeros(len(edges))\n", + " for ci, k in enumerate(order):\n", + " ax.barh(y, np.nan_to_num(frac[ci]), left=left, height=0.8,\n", + " color=CAUSE_COLORS[k], label=k)\n", + " left = left + np.nan_to_num(frac[ci])\n", + " ax.set(xlabel=\"Fraction of Mortality\", title=name, xlim=(0, 1))\n", + "axes[0].set_yticks(y)\n", + "axes[0].set_yticklabels(labels, fontsize=8)\n", + "axes[0].set_ylabel(\"DBH Size Class (cm)\")\n", + "axes[-1].legend(fontsize=8, framealpha=0.95, loc=\"center left\",\n", + " bbox_to_anchor=(1.02, 0.5))\n", + "fig.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "
\n", + "Click here for the solution
\n", + "\n", + "![plot example](../../../images/diagnostics/clm_ctsm/fates_bart_mortality_by_size.png)\n", + "\n", + "*

Figure: Mortality at three different stand ages and by size class.

*\n", + " \n", + "
\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Questions:\n", + " - Within a given size class, which cause dominates mortality? Is it the same across most size classes?\n", + " - Do any causes appear only in the smallest size classes? Why?\n", + " - Comparing the three stand ages, does the cause composition change much with succession, or mostly stay the same?" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 5.2. The patch mosaic\n", + "\n", + "Every plot so far has treated BART as a single stand, but FATES actually represents each gridcell as a mosaic of patches, each at a different time since it was last disturbed. When a canopy tree dies it can open a gap and start a new, younger patch. Undisturbed patches age and their canopies close. The gridcell-level numbers we've been plotting are area-weighted averages over this shifting distribution of patch ages. The plot below shows how the site's area is partitioned among patch-age classes as the forest develops." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "patch_area = ds_bg[\"FATES_PATCHAREA_AP\"].sel(fates_levage=slice(0, 100)) # we only allowed the forest to run for 100 years\n", + "\n", + "edges = patch_area.fates_levage.values\n", + "na = len(edges)\n", + "labels = [f\"{int(edges[i])}–{int(edges[i+1])}\" if i+1 < na else f\"{int(edges[i])}+\"\n", + " for i in range(na)]\n", + "\n", + "colors = cm.YlGn(np.linspace(0.25, 0.95, na))\n", + "tyear = ds_bg[\"time\"].dt.year + (ds_bg[\"time\"].dt.month - 1) / 12\n", + "stand_age = tyear - tyear.min()\n", + "\n", + "fig, ax = plt.subplots(1, 1, figsize=(7, 5))\n", + "ax.stackplot(stand_age.values, *[patch_area.isel(fates_levage=i).values for i in range(na)],\n", + " colors=colors, labels=[f\"{l} yr\" for l in labels])\n", + "ax.set(xlabel=\"Stand Age\", ylabel=\"Patch Area (fraction of site)\")\n", + "ax.margins(x=0)\n", + "ax.set_ylim(0, 1.0)\n", + "\n", + "ax.legend(loc=\"center left\", fontsize=7, title=\"Patch Age\",\n", + " ncol=2, bbox_to_anchor=(1.02, 0.5))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "
\n", + "Click here for the solution
\n", + "\n", + "![plot example](../../../images/diagnostics/clm_ctsm/fates_bart_patch_age.png)\n", + "\n", + "*

Figure: Patch area over time and across different patch ages.

*\n", + " \n", + "
\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Questions:\n", + " - Early in the run, which patch-age class holds nearly all of the site's area? Why does that make sense for a bare-ground start?\n", + " - Does a fraction of young-patch area persist even late in the run? Since there's no fire or harvest in this simulation, where do new young patches keep coming from?" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---\n", + "## 6. Next Steps\n", + "\n", + "As you can see, FATES demography can result in a lot of complicated ways to visualize a forest stand over time. Feel free to explore some other history variables or other plotting methods. You can see what variables are available to look at by looking at the `ds_bg.data_vars` or `ds_inv.data_vars`.\n", + "\n", + "In most of our plots, we only looked at things organized by size class (`_SZ`), even though there were many variables that were indexed by size *and* PFT (`_SZPF`). As a challenge, can you explore these history variables, but comparing PFT and/or size *and* PFT?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "NPL 2023a", + "language": "python", + "name": "npl-2023a" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.16" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}