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/challenge/clm_fates_lh_diff.png b/images/challenge/clm_fates_lh_diff.png new file mode 100644 index 000000000..ce8e858a5 Binary files /dev/null and b/images/challenge/clm_fates_lh_diff.png differ diff --git a/images/challenge/clm_sp_bgc_diff.png b/images/challenge/clm_sp_bgc_diff.png new file mode 100644 index 000000000..f8d7dac0b Binary files /dev/null and b/images/challenge/clm_sp_bgc_diff.png differ diff --git a/images/challenge/clm_sp_plot_monthly_gpp.png b/images/challenge/clm_sp_plot_monthly_gpp.png new file mode 100644 index 000000000..04100f8f0 Binary files /dev/null and b/images/challenge/clm_sp_plot_monthly_gpp.png differ diff --git a/images/challenge/fates_bart_canopy_structure.png b/images/challenge/fates_bart_canopy_structure.png new file mode 100644 index 000000000..fb6f5dd64 Binary files /dev/null and b/images/challenge/fates_bart_canopy_structure.png differ diff --git a/images/challenge/fates_bart_gpp_lai.png b/images/challenge/fates_bart_gpp_lai.png new file mode 100644 index 000000000..55dde855a Binary files /dev/null and b/images/challenge/fates_bart_gpp_lai.png differ diff --git a/images/challenge/fates_bart_growth_mortality.png b/images/challenge/fates_bart_growth_mortality.png new file mode 100644 index 000000000..5db5ecaa2 Binary files /dev/null and b/images/challenge/fates_bart_growth_mortality.png differ diff --git a/images/challenge/fates_bart_layered_LAI_bg.png b/images/challenge/fates_bart_layered_LAI_bg.png new file mode 100644 index 000000000..a86504031 Binary files /dev/null and b/images/challenge/fates_bart_layered_LAI_bg.png differ diff --git a/images/challenge/fates_bart_layered_LAI_init.png b/images/challenge/fates_bart_layered_LAI_init.png new file mode 100644 index 000000000..523eb1d23 Binary files /dev/null and b/images/challenge/fates_bart_layered_LAI_init.png differ diff --git a/images/challenge/fates_bart_mortality.png b/images/challenge/fates_bart_mortality.png new file mode 100644 index 000000000..ccd984314 Binary files /dev/null and b/images/challenge/fates_bart_mortality.png differ diff --git a/images/challenge/fates_bart_mortality_by_size.png b/images/challenge/fates_bart_mortality_by_size.png new file mode 100644 index 000000000..8bf1e58c2 Binary files /dev/null and b/images/challenge/fates_bart_mortality_by_size.png differ diff --git a/images/challenge/fates_bart_mortality_snapshots.png b/images/challenge/fates_bart_mortality_snapshots.png new file mode 100644 index 000000000..5b75e1b65 Binary files /dev/null and b/images/challenge/fates_bart_mortality_snapshots.png differ diff --git a/images/challenge/fates_bart_patch_age.png b/images/challenge/fates_bart_patch_age.png new file mode 100644 index 000000000..628574373 Binary files /dev/null and b/images/challenge/fates_bart_patch_age.png differ diff --git a/images/challenge/fates_bart_promotion_demotion.png b/images/challenge/fates_bart_promotion_demotion.png new file mode 100644 index 000000000..746826c4c Binary files /dev/null and b/images/challenge/fates_bart_promotion_demotion.png differ diff --git a/images/challenge/fates_bart_stand_structure.png b/images/challenge/fates_bart_stand_structure.png new file mode 100644 index 000000000..b7d65e72a Binary files /dev/null and b/images/challenge/fates_bart_stand_structure.png differ diff --git a/images/challenge/fates_bart_vegc_composition_bg.png b/images/challenge/fates_bart_vegc_composition_bg.png new file mode 100644 index 000000000..65e3f0b3b Binary files /dev/null and b/images/challenge/fates_bart_vegc_composition_bg.png differ diff --git a/images/challenge/fates_bart_vegc_composition_init.png b/images/challenge/fates_bart_vegc_composition_init.png new file mode 100644 index 000000000..063a7777e Binary files /dev/null and b/images/challenge/fates_bart_vegc_composition_init.png differ diff --git a/images/challenge/fates_clm_gpp_diff.png b/images/challenge/fates_clm_gpp_diff.png new file mode 100644 index 000000000..fefc56d26 Binary files /dev/null and b/images/challenge/fates_clm_gpp_diff.png differ diff --git a/images/challenge/fates_monthly_gpp_pft.png b/images/challenge/fates_monthly_gpp_pft.png new file mode 100644 index 000000000..481592c8a Binary files /dev/null and b/images/challenge/fates_monthly_gpp_pft.png differ diff --git a/images/challenge/fates_plot_monthly_gpp.png b/images/challenge/fates_plot_monthly_gpp.png new file mode 100644 index 000000000..f837e4e1c Binary files /dev/null and b/images/challenge/fates_plot_monthly_gpp.png differ diff --git a/images/challenge/fates_vcmax_diff.png b/images/challenge/fates_vcmax_diff.png new file mode 100644 index 000000000..4a998ec24 Binary files /dev/null and b/images/challenge/fates_vcmax_diff.png differ diff --git a/images/challenge/fates_vcmax_lh_diff.png b/images/challenge/fates_vcmax_lh_diff.png new file mode 100644 index 000000000..7b2b921f8 Binary files /dev/null and b/images/challenge/fates_vcmax_lh_diff.png differ diff --git a/notebooks/challenge/clm_ctsm/clm_ctsm.ipynb b/notebooks/challenge/clm_ctsm/clm_ctsm.ipynb index fd254f236..0cef3f3a7 100644 --- a/notebooks/challenge/clm_ctsm/clm_ctsm.ipynb +++ b/notebooks/challenge/clm_ctsm/clm_ctsm.ipynb @@ -11,7 +11,7 @@ "\n", "It can be useful for people interested in land science to run simulations with only an active land components and atmospheric forcing. In this exercise, you will learn how to run one of these land-only simulations.\n", "\n", - "This exercise was created by Peter Lawrence, Erik Kluzek and Alice DuVivier." + "These exercises were created by Adrianna Foster and Erik Kluzek." ] }, { @@ -27,9 +27,10 @@ "id": "346cbd7b-3b8e-41f0-b120-b369ab20f6cc", "metadata": {}, "source": [ - "- Student will learn what a I compset is, the types of forcing available to run one, and how to run one.\n", - "- Student will learn how to run a case with satellite phenology and one with prognostic crops and compare the two experiments.\n", - "- Learn what inputs CLM needs and what they look like.\n" + "- Students will learn what an \"I\" compset is, the types of forcing available to run one, and how to run one.\n", + "- Students will learn how to run a default CLM (\"big leaf\") case with satellite phenology and one with prognostic BGC and crops and compare the two experiments.\n", + "- Students will learn how to run FATES and compare to big leaf CLM output.\n", + "- Students will learn how to modify an input parameter file and investigate its effect on output." ] }, { @@ -46,10 +47,10 @@ "metadata": {}, "source": [ "- This exercise uses the same code base as the rest of the tutorial. \n", - "- You will be using the I2000Clm50Sp and IHistClm50BgcCrop compsets at the f09_g17_gl4 resolution.\n", - "- You will run a control simulation and two experimental simulations. \n", - "- You will modify netcdf input files. \n", - "- You will use simple, command line netcdf tools to evaluate how the experiments differ from the control simulation." + "- You will be using the I2000Clm60Sp, IHistClm60BgcCrop, and I2000Clm60FatesSpCrujraRsGs compsets at the f19_f19_mt233 resolution.\n", + "- You will run a CLM SP simulation, a CLM BGC simulation, and two FATES SP simulations.\n", + "- You will modify a json input file.\n", + "- You will use the `xarray` and `matplotlib` libraries to evaluate how the simulations differ." ] }, { @@ -57,7 +58,7 @@ "id": "fd2ff959-6ec0-4534-8713-f5dfdb13b955", "metadata": {}, "source": [ - "## Useful CLM references" + "## Useful CLM and FATES references" ] }, { @@ -67,7 +68,7 @@ "source": [ "
\n", "\n", - "[CESM CLM/CTSM User's Guide](https://escomp.github.io/ctsm-docs/versions/master/html/users_guide/index.html)\n", + "[CESM CLM/CTSM User's Guide](https://escomp.github.io/CTSM/users_guide/index.html)\n", "\n", "
" ] @@ -86,6 +87,18 @@ "" ] }, + { + "cell_type": "markdown", + "id": "e9024cab-f818-4d8a-99bb-720a69bc19ce", + "metadata": {}, + "source": [ + "
\n", + "\n", + "[FATES Users Guide](https://fates-users-guide.readthedocs.io/en/latest/)\n", + "\n", + "
" + ] + }, { "cell_type": "markdown", "id": "cb8a0e10-a920-41dd-8bc4-5c03395d231e", @@ -99,11 +112,12 @@ "id": "69988cae-fd26-4fda-8fee-6ab52e79f42d", "metadata": {}, "source": [ - "The I compset has active clm with a data atmosphere. The sea ice, ocean, atmosphere, and wave models are not active. There are two types of options for I compsets:\n", - "- SP: Satellite Phenology\n", - "- BGC: Biogeochemistry \n", + "\"I\" compsets have active CLM with a data atmosphere (DATM). The sea ice, ocean, atmosphere, and wave models are not active. There are three types of options for I compsets:\n", + "- SP: Satellite Phenology (with CLM big leaf)\n", + "- BGC: Biogeochemistry (with CLM big leaf)\n", + "- FATES: The FATES cohort vegetation model instead of big leaf, either SP or BGC\n", "\n", - "We will use the Satellite Phenology option for our control case." + "We will use the Satellite Phenology option for our first simulation." ] }, { @@ -113,7 +127,7 @@ "source": [ "![icase](../../../images/challenge/i2000.png)\n", "\n", - "*

Figure: I2000 compset definition.

*" + "*

Figure: I2000 compset definition.

*" ] }, { @@ -135,9 +149,7 @@ "source": [ "## Post processing and viewing your output\n", "\n", - "You will use [ncview](https://ncar.github.io/CESM-Tutorial/notebooks/tools/netcdf_tools.html#ncview) and [NCO operator](https://ncar.github.io/CESM-Tutorial/notebooks/tools/netcdf_tools.html#netcdf-operators-nco) tools to evaluate how the experiments differ from the control simulation.\n", - "\n", - "These tools will need to be loaded into your environment using the NCAR HPC [modules](https://ncar.github.io/CESM-Tutorial/notebooks/tools/unix/modules.html)." + "You will use python and the `xarray` and `matplotlib` libraries to evaluate how the experiments differ." ] }, { @@ -145,22 +157,20 @@ "id": "828e538d-3558-4747-9b61-3b29da04e294", "metadata": {}, "source": [ - "1) You can create an annual average of the first year's data for each simulationg using the `ncra` (netCDF averager) command from the netCDF operators package \\([NCO](https://nco.sourceforge.net/)\\). \n", + "1) You can read in multiple history files using the [xarray](https://docs.xarray.dev/en/stable/) command:\n", "```\n", - "ncra $OUTPUT_DIR/*.clm.*nc $CASENAME.cln.h.0001.nc\n", + "ds = open_mfdataset(files)\n", "```\n", "\n", - "2) Create a file that contains differences between each of the experiments and the control simulation\n", + "2) Look at differences between the same variable on two different xarray datasets:\n", "```\n", - "ncdiff i.day5.b.clm2.XXX.nc /glade/derecho/scratch/$user/archive/i.day5.a/lnd/hist/i.day5.a.clm2.XXX.nc i_diff.nc\n", + "gpp_diff = ds_1.GPP - ds_2.GPP\n", "```\n", "\n", - "3) Examine variables within each annual mean and the difference files using `ncview`\n", - "```\n", - "ncview i_diff.nc\n", + "3) Examine variables using the `.plot()` function\n", "```\n", - "\n", - "4) You can also look at other monthly-mean outputs or component log files." + "ds.GPP.plot()\n", + "```" ] }, { @@ -174,9 +184,9 @@ ], "metadata": { "kernelspec": { - "display_name": "CMIP6 2019.10", + "display_name": "Python 3 (ipykernel)", "language": "python", - "name": "cmip6-201910" + "name": "python3" }, "language_info": { "codemirror_mode": { @@ -188,7 +198,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_1.ipynb b/notebooks/challenge/clm_ctsm/clm_exercise_1.ipynb index c3fa9a716..e50482abf 100644 --- a/notebooks/challenge/clm_ctsm/clm_exercise_1.ipynb +++ b/notebooks/challenge/clm_ctsm/clm_exercise_1.ipynb @@ -5,7 +5,9 @@ "id": "f406f992-92bd-4b17-9bd3-b99c5c8abaf3", "metadata": {}, "source": [ - "# 1: Control case" + "# 1: CLM SP\n", + "\n", + "We will first run CLM in satellite phenology (SP) mode. In SP mode, CLM prescribes vegetation state (i.e., leaf area, height, and stem area) from satellite observations rather than growing it. The model only computes the biophysical exchange of energy, water, and carbon over that fixed canopy, which makes SP runs cheap and a good baseline when you want realistic surface fluxes without a prognostic carbon cycle." ] }, { @@ -14,45 +16,105 @@ "metadata": {}, "source": [ "
\n", - "Exercise: Run a control case

\n", + "Exercise: Run a CLM SP case

\n", " \n", - "Create a case called **i.day5.a** using the compset `I2000Clm50Sp` at `f09_g17_gl4` resolution. \n", + "Create a case called **i.clm.sp.year1** using the compset `I2000Clm60Sp` at `f19_f19_mt233` resolution. \n", " \n", - "Set the run length to **5 days**. \n", + "Set the run length to **1 year**.\n", "\n", - "Build and run the model. Since this is a control case, we want to build it \"out of the box\" without any modifications. \n", + "Build and run the model.\n", "\n", "
\n" ] }, { "cell_type": "markdown", - "id": "f639e182-f48a-431c-a594-9c34323417eb", + "id": "656858f2-ff60-44cc-b9fb-a48cd4ed3506", "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 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.clm.sp.year1\n", + "cd lnd/hist\n", + "```\n", + "\n", + "Check that your archive directory contains the files:\n", + "\n", + "```\n", + "i.clm.sp.year1.clm2.h0a.2000-01.nc i.clm.sp.year1.clm2.h0i.2000-01.nc\n", + "i.clm.sp.year1.clm2.h0a.2000-02.nc i.clm.sp.year1.clm2.h0i.2000-02.nc\n", + "i.clm.sp.year1.clm2.h0a.2000-03.nc i.clm.sp.year1.clm2.h0i.2000-03.nc\n", + "i.clm.sp.year1.clm2.h0a.2000-04.nc i.clm.sp.year1.clm2.h0i.2000-04.nc\n", + "i.clm.sp.year1.clm2.h0a.2000-05.nc i.clm.sp.year1.clm2.h0i.2000-05.nc\n", + "i.clm.sp.year1.clm2.h0a.2000-06.nc i.clm.sp.year1.clm2.h0i.2000-06.nc\n", + "i.clm.sp.year1.clm2.h0a.2000-07.nc i.clm.sp.year1.clm2.h0i.2000-07.nc\n", + "i.clm.sp.year1.clm2.h0a.2000-08.nc i.clm.sp.year1.clm2.h0i.2000-08.nc\n", + "i.clm.sp.year1.clm2.h0a.2000-09.nc i.clm.sp.year1.clm2.h0i.2000-09.nc\n", + "i.clm.sp.year1.clm2.h0a.2000-10.nc i.clm.sp.year1.clm2.h0i.2000-10.nc\n", + "i.clm.sp.year1.clm2.h0a.2000-11.nc i.clm.sp.year1.clm2.h0i.2000-11.nc\n", + "i.clm.sp.year1.clm2.h0a.2000-12.nc i.clm.sp.year1.clm2.h0i.2000-12.nc\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 i.day5.a with the command:\n", + "Click here for the solution
\n", + "\n", + "**Create a new case**\n", + "\n", + "Create a new case i.clm.sp.year1 with the command:\n", "```\n", "cd /glade/u/home/$USER/code/my_cesm_code/cime/scripts\n", - "./create_newcase --case ~/cases/i.day5.a --compset I2000Clm50Sp --res f09_g17_gl4 --run-unsupported\n", + "./create_newcase --case ~/cases/i.clm.sp.year1 --compset I2000Clm60Sp --res f19_f19_mt233 --run-unsupported\n", "```\n", "
\n", "\n", - "Case setup:\n", + "**Case setup**\n", + "\n", "``` \n", - "cd ~/cases/i.day5.a \n", + "cd ~/cases/i.clm.sp.year1\n", "./case.setup\n", "```\n", + "\n", "
\n", "\n", - "Change the clm namelist using user_nl_clm by adding the following lines:\n", + "**Set run length**\n", + " \n", + "Change the `run length`:\n", "``` \n", - "hist_nhtfrq = -24\n", - "hist_mfilt = 6\n", + "./xmlchange STOP_N=1,STOP_OPTION=nyears\n", "```\n", "
\n", " \n", @@ -62,38 +124,61 @@ "```\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", + "\n", + "For instance, to run in the queue `tutorial` and the project number `UESM0016`:\n", "``` \n", - "./xmlchange JOB_QUEUE=tutorial,PROJECT=UESM0016 --force,JOB_WALLCLOCK_TIME=0:15:00\n", + "./xmlchange JOB_QUEUE=tutorial,PROJECT=UESM0016 --force\n", "```\n", "
\n", "\n", - "Build and submit:\n", + "**Change the wallclock time**\n", + "\n", + "```\n", + "./xmlchange --subgroup case.run JOB_WALLCLOCK_TIME=01:00:00\n", + "```\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", + "**Check the run:**\n", "When the run is completed, look into the archive directory for: \n", - "i.day5.a. \n", + "i.clm.sp.year1. \n", " \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/lnd/hist\n", + "cd /glade/derecho/scratch/$USER/archive/i.clm.sp.year1/lnd/hist\n", "\n", "ls \n", "```\n", - "
\n", - "\n", - "(2) Look at the output using ncview\n" + "
" + ] + }, + { + "cell_type": "markdown", + "id": "289c5251-a0bb-4a9d-9b40-d320953c621a", + "metadata": {}, + "source": [ + "We will investigate the output in a subsequent exercise." ] }, { "cell_type": "code", "execution_count": null, - "id": "17cf7e19-1211-45f2-97cd-fe2badc69dac", + "id": "de84fb8e-aa52-4bf6-b7c7-583af74b1328", "metadata": {}, "outputs": [], "source": [] @@ -101,9 +186,9 @@ ], "metadata": { "kernelspec": { - "display_name": "CMIP6 2019.10", + "display_name": "Python 3 (ipykernel)", "language": "python", - "name": "cmip6-201910" + "name": "python3" }, "language_info": { "codemirror_mode": { @@ -115,7 +200,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_2.ipynb b/notebooks/challenge/clm_ctsm/clm_exercise_2.ipynb index 1398f0dd6..7cbee1ef5 100644 --- a/notebooks/challenge/clm_ctsm/clm_exercise_2.ipynb +++ b/notebooks/challenge/clm_ctsm/clm_exercise_2.ipynb @@ -5,7 +5,7 @@ "id": "f406f992-92bd-4b17-9bd3-b99c5c8abaf3", "metadata": {}, "source": [ - "# 2: Use the BGC model" + "# 2: CLM BGC" ] }, { @@ -13,7 +13,11 @@ "id": "0037b73f-f174-48e7-8e4f-0744d7d23fe0", "metadata": {}, "source": [ - "We can use a different I compset: IHistClm50BgcCrop. This experiment is a 20th century transient run using GSWP3v1 and the biogeochemistry model including crops." + "We can use a different I compset: IHistClm60BgcCrop. This experiment is a 20th century transient run using CRUJRA and the biogeochemistry model including crops.\n", + "\n", + "In biogeochemistry (BGC) mode, CLM grows vegetation rather than prescribing it, computing photosynthesis, allocating carbon to leaves, stems, and roots, and cycling carbon and nitrogen through vegetation, litter, and soil. Leaf area becomes a prognostic result of that carbon balance rather than a satellite input, which makes BGC far better for studying how vegetation responds to change, but costlier, since the carbon pools must be spun up to equilibrium first.\n", + "\n", + "By default, IHIST compsets have a default `finidat` (initialization file) that they use to initialize the model. We will leave this at default." ] }, { @@ -23,7 +27,7 @@ "source": [ "![icase](../../../images/challenge/ihist.png)\n", "\n", - "*

Figure: IHIST compset definition.

*" + "*

Figure: IHIST compset definition.

*" ] }, { @@ -34,15 +38,71 @@ "
\n", "Exercise: Run an experimental case with prognostic BGC

\n", " \n", - "Create a case called **i.day5.b** using the compset `IHistClm50BgcCrop` at `f09_g17_gl4` resolution. \n", + "Create a case called **i.clm.bgc.month3** using the compset `IHistClm60BgcCrop` at `f19_f19_mt233` resolution. \n", " \n", - "Set the run length to **5 days**. \n", + "Set the run length to **3 months**. \n", "\n", "Build and run the model.\n", " \n", "
" ] }, + { + "cell_type": "markdown", + "id": "846e0789-67bd-4186-b02f-1867793d1659", + "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 change the run length to 3 months?**\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.clm.bgc.month3\n", + "cd lnd/hist\n", + "```\n", + "\n", + "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", + "
" + ] + }, { "cell_type": "markdown", "id": "f639e182-f48a-431c-a594-9c34323417eb", @@ -52,32 +112,36 @@ "\n", "
\n", "
\n", - "Click here for the solution
\n", - " \n", - "Create a new case i.day5.b :\n", + "Click here for the solution
\n", + "\n", + "**Create a new case**\n", + "\n", + "Create a new case i.clm.bgc.month3 :\n", "```\n", "cd /glade/u/home/$USER/code/my_cesm_code/cime/scripts\n", - "./create_newcase --case ~/cases/i.day5.b --compset IHistClm50BgcCrop --res f09_g17_gl4 --run-unsupported\n", + "./create_newcase --case ~/cases/i.clm.bgc.month3 --compset IHistClm60BgcCrop --res f19_f19_mt233 --run-unsupported\n", "```\n", "
\n", "\n", - "Case setup:\n", + "**Case setup**\n", + "\n", "``` \n", - "cd ~/cases/i.day5.b\n", + "cd ~/cases/i.clm.bgc.month3\n", "./case.setup\n", "```\n", "
\n", "\n", - "Note differences between this case and the control case:\n", + "**Set run length**\n", + " \n", + "Change the `run length`:\n", "``` \n", - "diff env_run.xml ../i.day5.a/env_run.xml\n", + "./xmlchange STOP_N=3,STOP_OPTION=nmonths\n", "```\n", "
\n", "\n", - "Change the clm namelist using user_nl_clm by adding the following lines:\n", + "Note differences between this case and the SP case:\n", "``` \n", - "hist_nhtfrq = -24\n", - "hist_mfilt = 6\n", + "diff env_run.xml ../i.clm.sp.year1/env_run.xml\n", "```\n", "
\n", " \n", @@ -87,14 +151,26 @@ "```\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", + "\n", + "For instance, to run in the queue `tutorial` and the project number `UESM0016`:\n", "``` \n", - "./xmlchange JOB_QUEUE=tutorial,PROJECT=UESM0016 --force,JOB_WALLCLOCK_TIME=0:15:00\n", + "./xmlchange JOB_QUEUE=tutorial,PROJECT=UESM0016 --force\n", "```\n", "
\n", "\n", - "Build case:\n", + "**Change the wallclock time**\n", + "\n", + "```\n", + "./xmlchange --subgroup case.run JOB_WALLCLOCK_TIME=01:00:00\n", + "```\n", + "
\n", + "\n", + "\n", + "**Build case**\n", + "\n", "```\n", "qcmd -- ./case.build\n", "```\n", @@ -102,63 +178,44 @@ " \n", "Compare the namelists from the two experiments:\n", "```\n", - "diff CaseDocs/lnd_in ../i.day5.a/CaseDocs/lnd_in\n", + "diff CaseDocs/lnd_in ../i.clm.sp.year1/CaseDocs/lnd_in\n", "```\n", "
\n", " \n", - "Submit case:\n", + "**Submit case**\n", + "\n", "```\n", "./case.submit\n", "```\n", "
\n", "\n", - "When the run is completed, look into the archive directory for: \n", - "i.day5.b. \n", - " \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.b/lnd/hist\n", + "**Check the run:**\n", "\n", - "ls \n", - "```\n", - "
\n", + "When the run is completed, look into the archive directory for: \n", + "i.clm.bgc.month3. \n", "\n", - " \n", - "(2) Compare to control run:\n", + "Check that your archive directory on derecho (The path will be different on other machines): \n", "```\n", - "ncdiff -v TLAI i.day5.b.clm2.XXX.nc /glade/derecho/scratch/$USER/archive/i.day5.a/lnd/hist/i.day5.a.clm2.XXX.nc i_diff.nc\n", + "cd /glade/derecho/scratch/$USER/archive/i.clm.bgc.month3/lnd/hist\n", "\n", - "ncview i_diff.nc\n", + "ls \n", "```\n", - "\n", - "\n", "
\n", - "
\n" - ] - }, - { - "cell_type": "markdown", - "id": "d69c456a-fdc6-4625-bbcc-ed32ab6ae8e8", - "metadata": { - "tags": [] - }, - "source": [ - "## Test your understanding" + "
" ] }, { "cell_type": "markdown", - "id": "b3f28ecf-02b8-4cc0-bdc4-c36c8fc9e7aa", + "id": "d6066cbd-09bc-4dfe-945e-d0f5506e8b36", "metadata": {}, "source": [ - "- What changes do you see from the control case with the prognostic BGC?\n", - "- ... OTHERS?" + "We will investigate the output in a subsequent exercise." ] }, { "cell_type": "code", "execution_count": null, - "id": "5b5a8cee-0ca9-4076-a731-e6ec200b70d4", + "id": "4e3a0a21-7f9c-4621-b4fe-7b696eb3a7f9", "metadata": {}, "outputs": [], "source": [] @@ -166,9 +223,9 @@ ], "metadata": { "kernelspec": { - "display_name": "CMIP6 2019.10", + "display_name": "Python 3 (ipykernel)", "language": "python", - "name": "cmip6-201910" + "name": "python3" }, "language_info": { "codemirror_mode": { @@ -180,7 +237,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_3.ipynb b/notebooks/challenge/clm_ctsm/clm_exercise_3.ipynb index 8df149c34..daa52ed4f 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,15 @@ "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` in `user_nl_clm`:\n", + "\n", + "`finidat = '/glade/campaign/cesm/tutorial/diagnostics_tutorial_archive/i.fates.sp.spinup.update.clm2.r.2060-01-01-00000.nc'`\n", "\n", "Build and run the model.\n", " \n", @@ -42,132 +91,427 @@ }, { "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 set the finidat?**\n", + "\n", + "After running `./case.setup`, open the `user_nl_clm` file using any text editing software and add:\n", + "\n", + "`finidat = '/glade/campaign/cesm/tutorial/diagnostics_tutorial_archive/i.fates.sp.spinup.update.clm2.r.2060-01-01-00000.nc'`\n", + "\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 `h0a` 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 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", + "**Set finidat**\n", + "\n", + "After running `./case.setup`, open the `user_nl_clm` file using any text editing software and add:\n", + "\n", + "`finidat = '/glade/campaign/cesm/tutorial/diagnostics_tutorial_archive/i.fates.sp.spinup.update.clm2.r.2060-01-01-00000.nc'`\n", + "\n", + "
\n", + "\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 ``UESM0016`` (You should use the project number given for this tutorial), use the command:\n", "``` \n", - "./xmlchange JOB_QUEUE=tutorial,PROJECT=UESM0016 --force,JOB_WALLCLOCK_TIME=0:15:00\n", + "./xmlchange JOB_QUEUE=tutorial,PROJECT=UESM0016 --force\n", "```\n", + "\n", "
\n", "\n", - "Build case:\n", + "**Change the wallclock time**\n", + "\n", "```\n", - "qcmd -- ./case.build\n", + "./xmlchange --subgroup case.run JOB_WALLCLOCK_TIME=02: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. \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", + "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", - "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", "\n", - "ncview i_diff.nc\n", + "(2) Investigate the contents of one of the `h0a` 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", "
" ] }, { "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": [ + "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": [ - "## Test your understanding" + "
\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.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/u/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 ``UESM0016`` (You should use the project number given for this tutorial), use the command:\n", + "``` \n", + "./xmlchange JOB_QUEUE=tutorial,PROJECT=UESM0016 --force\n", + "```\n", + "\n", + "
\n", + "\n", + "**Change the wallclock time**\n", + "\n", + "```\n", + "./xmlchange --subgroup case.run JOB_WALLCLOCK_TIME=02: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. \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 +519,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 +533,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..ac63eeb13 --- /dev/null +++ b/notebooks/challenge/clm_ctsm/clm_exercise_4.ipynb @@ -0,0 +1,640 @@ +{ + "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": [ + "## CLM and CLM-FATES output\n", + "\n", + "You've run four cases spanning two land model structures and two vegetation configurations:\n", + "\n", + "* CLM-SP: standard \"big-leaf\" CLM with satellite-prescribed phenology\n", + "* CLM-BGC: big-leaf CLM with prognostic biogeochemistry (carbon/nitrogen cycling, prognostic LAI)\n", + "* FATES-SP: FATES with satellite-prescribed phenology\n", + "* FATES-SP (modified parameters): the same FATES-SP configuration with a perturbed parameter file\n", + "\n", + "Together these let you make three comparisons, each isolating a different piece:\n", + "\n", + "1. CLM-SP vs. FATES-SP: same meteorology and prescribed phenology, different vegetation model.\n", + "2. FATES-SP vs. FATES-SP (modified): how FATES output responds to a parameter change.\n", + "3. CLM-SP vs. CLM-BGC: prescribed vs. prognostic vegetation within big-leaf CLM.\n", + "\n", + "While CLM-FATES and big-leaf CLM share many output variables (mostly CLM quantities like `ASA`, `FSR`, and `FSH` that don't depend on the vegetation model) most vegetation-related quantities in a FATES run are simulated, calculated, and output by FATES itself. FATES variables carry their own names (typically prefixed `FATES_`), and their dimensions differ from big-leaf CLM: instead of the CLM PFT axis, FATES output is organized by FATES PFTs, size classes, and patch/cohort structure. Even the PFT definitions themselves differ between the two.\n", + "\n", + "This notebook walks through some basic CLM and CLM-FATES output for the global cases you just ran, and shows how to do some basic comparisons." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "*We will 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", + "import matplotlib.pyplot as plt\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 exercises below assume you use `i.clm.sp`, `i.clm.bgc`, and `i.fates.year1` and `i.fates.year1.vcmax` output from the previous challenges 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 challenge exercises. For these examples we will use gross primary production, latent heat flux, and sensible heat flux.\n", + "\n", + "Two of these are CLM variables that share a name across all four cases:\n", + "- latent heat flux (`EFLX_LH_TOT`)\n", + "- sensible heat flux (`FSH`)\n", + "\n", + "GPP is where the model structural differences show up. The variable name differs across all three configurations:\n", + "\n", + "- **CLM-SP**: `FPSN` - photosynthesis rate (SP mode doesn't run the full carbon cycle, so there is no prognostic \"GPP\")\n", + "- **CLM-BGC**: `GPP` - gross primary production from the active biogeochemistry\n", + "- **FATES-SP**, **FATES-SP (modified)**: `FATES_GPP` - GPP calculated output by FATES\n", + "\n", + "Note that the units for these variables are also all different! We will deal with this later..." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# note: change the user here to your user name\n", + "user = 'afoster'\n", + "\n", + "# per-case config: archive case name + the GPP-like variable\n", + "cases = {\n", + " 'clm_sp': {'case': 'i.clm.sp.year1', 'gpp': 'FPSN'}, # SP big-leaf: photosynthesis rate\n", + " 'clm_bgc': {'case': 'i.clm.bgc.month3', 'gpp': 'GPP'}, # prognostic BGC: true GPP\n", + " 'fates_sp': {'case': 'i.fates.year1', 'gpp': 'FATES_GPP'}, # FATES control\n", + " 'fates_vcmax': {'case': 'i.fates.year1.vcmax', 'gpp': 'FATES_GPP'}, # FATES, vcmax modified\n", + "}\n", + "\n", + "# variables shared by every case\n", + "shared_vars = ['EFLX_LH_TOT', 'FSH', 'area', 'landfrac']\n", + "\n", + "# FATES-only: gridcell fraction occupied by FATES\n", + "fates_vars = ['FATES_FRACTION']\n", + "\n", + "# read in all datasets, supplying the correct `data_vars`\n", + "datasets = {}\n", + "ref_coords = None\n", + "for name, cfg in cases.items():\n", + " hist = f\"/glade/derecho/scratch/{user}/archive/{cfg['case']}/lnd/hist\"\n", + " files = sorted(glob.glob(os.path.join(hist, f\"{cfg['case']}.clm2.h0a.*\")))\n", + "\n", + " data_vars = shared_vars + [cfg['gpp']]\n", + " if cfg['gpp'].startswith('FATES'):\n", + " data_vars += fates_vars\n", + " \n", + " ds = xr.open_mfdataset(\n", + " files,\n", + " preprocess=functools.partial(preprocess, data_vars=data_vars),\n", + " )\n", + " datasets[name] = ds" + ] + }, + { + "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 shared_vars list in the cell above (for variables common to every case), the fates_vars list for variables only for FATES, or the cases config if you want to change the per-case GPP variable. Each case's read list is assembled from those.\n", + "
\n", + "\n", + "#### Printing information about the dataset is helpful for understanding your data.\n", + "\n", + "Let's start with one of the FATES cases:\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?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# print information about the FATES control dataset\n", + "datasets['fates_sp']" + ] + }, + { + "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`) 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": [ + "datasets['fates_sp'].FATES_GPP" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 1.2 Simple Calculations\n", + "\n", + "Before we can compare GPP across the four runs, we need to get them into the correct units and format. Two corrections are needed:\n", + "\n", + "1. **`FATES_FRACTION` correction (FATES cases only)**: FATES doesn't know how much of each gridcell is occupied by natural vegetation, so the gridcell average implicitly assumes all of it is. We correct by multiplying the FATES history area-weighted values by `FATES_FRACTION`.\n", + "\n", + "2. **Unit conversion (all cases)**: The three GPP-like variables come out in different units (`FPSN` as a photosynthesis rate in μmol CO2 m-2 s-1, `GPP` in gC m-2 s-1, and `FATES_GPP` in kgC m-2 s-1) We convert all three to a common gC m-2 s-1 so the runs are actually comparable.\n", + "\n", + "We'll store the result under a common variable name (`GPP`) in every dataset." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "GPP_CONVERSION = {\n", + " 'FPSN': 1e-6 * 12.011, # umol CO2 m-2 s-1 -> gC m-2 s-1 (1e-6 mol/umol * 12.011 gC/mol)\n", + " 'GPP': 1.0, # gC m-2 s-1 -> gC m-2 s-1 (already there)\n", + " 'FATES_GPP': 1000.0, # kg C m-2 s-1 -> gC m-2 s-1 (1000 g/kg)\n", + "}\n", + "\n", + "for name, cfg in cases.items():\n", + " ds = datasets[name]\n", + " gpp_orig = cfg['gpp'] # native GPP-like variable for this case\n", + " da = ds[gpp_orig]\n", + "\n", + " # FATES only: correct for FATES fraction\n", + " if gpp_orig.startswith('FATES'):\n", + " da = da * ds['FATES_FRACTION']\n", + "\n", + " # unify units -> gC m-2 s-1\n", + " gpp = da * GPP_CONVERSION[gpp_orig]\n", + " gpp.attrs = dict(ds[gpp_orig].attrs) # copy source attrs\n", + " gpp.attrs['units'] = 'gC m-2 s-1' # then fix the ones we changed\n", + " gpp.attrs['long_name'] = 'gross primary production'\n", + " ds['GPP'] = gpp # common name across all four cases" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "After this, every dataset has a `GPP` in identical units, and we can now plot with this variable instead of using a different one for each dataset." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---\n", + "## 2. Basic plotting\n", + "\n", + "### 2.1. Easy plots using xarray\n", + "\n", + "To get a first look at the data, we can plot a year of data from one of the simulations. We'll start with **CLM-SP**, the first case you ran.\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'll plot GPP (variable = `GPP`). Note that we select the variable by specifying our dataset (here `datasets['clm_sp']`) and then the variable. \n", + "- Because we standardized GPP in the previous step, this same code works for any of the four cases. Just swap `datasets['clm_sp']` for another key (e.g. `datasets['fates_sp']`).\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": [ + "datasets['fates_sp']['GPP'].plot(x=\"lon\", y=\"lat\", col=\"time\", col_wrap=4);" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "
\n", + "Click here for the solution
\n", + "\n", + "![plot example](../../../images/challenge/clm_sp_plot_monthly_gpp.png)\n", + "\n", + "*

Figure: Plotting solution for monthly GPP for our CLM SP 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", + "#### FATES-SP vs. FATES-SP (modified)\n", + "\n", + "We can calculate the difference between the FATES control run and the one where we updated the vcmax parameter for PFT 1. The code below:\n", + "\n", + " - Calculates annual GPP for both simulations (converting the standardized gC m-2 s-1 fluxes to an annual kgC m-2 yr-1 total)\n", + " - Weights by land fraction so we're only looking at land\n", + " - Defines the difference as a new variable, `gpp_diff`" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "gpp_cf = 24 * 60 * 60 * 1e-3\n", + "\n", + "fates_sp = datasets['fates_sp']\n", + "fates_vcmax = datasets['fates_vcmax']\n", + "\n", + "gpp_ann_control = annual_sum(fates_sp.GPP, gpp_cf) * fates_sp.landfrac.isel(time=0)\n", + "gpp_ann_vcmax = annual_sum(fates_vcmax.GPP, gpp_cf) * fates_vcmax.landfrac.isel(time=0)\n", + "\n", + "gpp_diff = gpp_ann_vcmax - gpp_ann_control\n", + "gpp_diff.attrs = {'units': 'kgC m-2 yr-1',\n", + " 'long_name': 'annual gross primary production difference (vcmax − control)'}\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/challenge/fates_vcmax_diff.png\")\n", + "\n", + "*

Figure: Plotting solution comparing the two FATES simulations.

*\n", + " \n", + "
\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Questions\n", + "- How is GPP different for FATES simulation with the updated parameter value, relative to the control FATES simulation?\n", + "- Where are the differences the largest?\n", + "- Where are they absent?\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": [ + "#### FATES-SP vs. CLM-SP \n", + "\n", + "Using the same methods, we can also compare our FATES-SP and CLM-SP cases:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "clm_sp = datasets['clm_sp']\n", + "\n", + "gpp_ann_clm = annual_sum(clm_sp.GPP, gpp_cf) * clm_sp.landfrac.isel(time=0)\n", + "gpp_ann_fates = annual_sum(fates_sp.GPP, gpp_cf) * fates_sp.landfrac.isel(time=0)\n", + "\n", + "gpp_diff_model = gpp_ann_fates - gpp_ann_clm\n", + "gpp_diff_model.attrs = {'units': 'kgC m-2 yr-1',\n", + " 'long_name': 'annual GPP difference (FATES-SP − CLM-SP)'}\n", + "gpp_diff_model.plot(cmap='RdBu_r');" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "
\n", + "Click here for the solution
\n", + "\n", + "![plot example](../../../images/challenge/fates_clm_gpp_diff.png)\n", + "\n", + "*

Figure: Plotting solution for comparing GPP between CLM-SP and FATES-SP.

*\n", + " \n", + "
\n", + "
" + ] + }, + { + "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", + "Unlike GPP, `EFLX_LH_TOT` is computed and gridcell-averaged by the CLM for every case, so it needs no `FATES_FRACTION` correction and no unit conversion. We take an annual mean here rather than a sum, since it's a flux density." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# latent heat flux difference: FATES control vs FATES vcmax \n", + "\n", + "lh_ann_control = annual_mean(fates_sp.EFLX_LH_TOT) * fates_sp.landfrac.isel(time=0)\n", + "lh_ann_vcmax = annual_mean(fates_vcmax.EFLX_LH_TOT) * fates_vcmax.landfrac.isel(time=0)\n", + "\n", + "lh_diff_param = lh_ann_vcmax - lh_ann_control\n", + "lh_diff_param.attrs = {'units': 'W m-2',\n", + " 'long_name': 'annual mean latent heat flux difference (vcmax − control)'}\n", + "lh_diff_param.plot(cmap='RdBu_r', vmin=-20, vmax=20);" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "
\n", + "Click here for the solution
\n", + "\n", + "![plot example](../../../images/challenge/fates_vcmax_lh_diff.png)\n", + "\n", + "*

Figure: Plotting solution for comparing latent heat flux between our two FATES cases.

*\n", + " \n", + "
\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# latent heat flux difference: CLM-SP vs FATES-SP\n", + "\n", + "lh_ann_clm = annual_mean(clm_sp.EFLX_LH_TOT) * clm_sp.landfrac.isel(time=0)\n", + "lh_ann_fates = fates_sp.EFLX_LH_TOT.mean('time') * fates_sp.landfrac.isel(time=0)\n", + "\n", + "lh_diff_model = lh_ann_fates - lh_ann_clm\n", + "lh_diff_model.attrs = {'units': 'W m-2',\n", + " 'long_name': 'annual mean latent heat flux difference (FATES-SP − CLM-SP)'}\n", + "lh_diff_model.plot(cmap='RdBu_r', vmin=-30, vmax=30);" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "
\n", + "Click here for the solution
\n", + "\n", + "![plot example](../../../images/challenge/clm_fates_lh_diff.png)\n", + "\n", + "*

Figure: Plotting solution for comparing latent heat flux between CLM-SP and FATES-SP.

*\n", + " \n", + "
\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Questions:\n", + "- How is latent heat flux different for the two FATES cases?\n", + "- What about FATES vs. CLM?\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": "markdown", + "metadata": {}, + "source": [ + "## 3. Comparing CLM-SP and CLM-BGC\n", + "\n", + "For the last comparison we will look at just CLM big leaf for prescribed vs prognostic vegetation. In CLM-SP, LAI and vegetation height are prescribed; in CLM-BGC, the carbon–nitrogen cycle runs prognostically and vegetation state evolves with the simulation.\n", + "\n", + "Additionally, the climate is different between the two runs. `IHIST` compsets by default starts in 1850 whereas the `I2000` compsets start in year 2000.\n", + "\n", + "The BGC case only ran for a few months, so we can't take an annual total the way we did for the SP runs. Instead we will compare the two runs directly over the months they share." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "clm_bgc = datasets['clm_bgc']\n", + "\n", + "# months available in the shorter BGC run\n", + "n_months = clm_bgc.time.size\n", + "\n", + "clm_sp_short = clm_sp.isel(time=slice(0, n_months))\n", + "\n", + "sec_per_day = 24 * 60 * 60\n", + "gpp_sp_mean = (clm_sp_short.GPP.groupby('time.month').mean() * sec_per_day) * clm_sp_short.landfrac.isel(time=0)\n", + "gpp_bgc_mean = (clm_bgc.GPP.groupby('time.month').mean() * sec_per_day) * clm_bgc.landfrac.isel(time=0)\n", + "\n", + "gpp_diff_phen = gpp_bgc_mean - gpp_sp_mean\n", + "gpp_diff_phen.attrs = {'units': 'gC m-2 day-1',\n", + " 'long_name': 'mean GPP difference over shared months (BGC − SP)'}\n", + "gpp_diff_phen.plot(cmap='RdBu_r', col='month', vmin=-5, vmax=5);" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "
\n", + "Click here for the solution
\n", + "\n", + "![plot example](../../../images/challenge/clm_sp_bgc_diff.png)\n", + "\n", + "*

Figure: Plotting solution for comparing GPP between CLM-SP and CLM-BGC

*\n", + " \n", + "
\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Questions:\n", + "- What is the primary difference in GPP between the IHIST BGC case and the I2000 SP case?\n", + "- Why might this be occuring? " + ] + }, + { + "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 +} 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..c94ad6e76 --- /dev/null +++ b/notebooks/challenge/clm_ctsm/clm_exercise_5.ipynb @@ -0,0 +1,1363 @@ +{ + "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": [ + "## 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. Starting familiar: canopy fluxes and leaf area\n", + "\n", + "Before we use any of FATES's structural detail, let's look at output that any land model (including CLM's default big-leaf vegetation) would produce: the seasonal cycles of GPP and LAI. Averaged across all years of the inventory-initialized run, these show the forest as one canopy, with deciduous leaves flushing in spring, carbon uptake peaking in summer, and both falling off through autumn." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# seasonal climatology: average each calendar month across all years\n", + "gpp_clim = ds_inv[\"FATES_GPP\"].groupby(\"time.month\").mean()*86400*1e3 # gC/m2/day\n", + "lai_clim = ds_inv[\"FATES_LAI\"].groupby(\"time.month\").mean()\n", + "\n", + "months = gpp_clim[\"month\"].values" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax1 = plt.subplots(figsize=(8, 4.5))\n", + "l1 = ax1.plot(months, gpp_clim.values, color=\"#2E6E5E\", marker=\"o\", label=\"GPP\")\n", + "ax1.set(xlabel=\"Month\", ylabel=\"GPP (gC m$^{-2}$ day$^{-1}$)\")\n", + "\n", + "ax2 = ax1.twinx()\n", + "l2 = ax2.plot(months, lai_clim.values, color=\"#9CA86E\", marker=\"s\", label=\"LAI\")\n", + "ax2.set_ylabel(\"LAI (m$^2$ m$^{-2}$)\")\n", + "\n", + "ax1.set_xticks(range(1, 13))\n", + "lines = l1 + l2\n", + "ax1.legend(lines, [line.get_label() for line in lines], loc=\"upper left\")\n", + "ax1.set_title(\"Seasonal cycle of GPP and LAI\")\n", + "fig.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "
\n", + "Click here for the solution
\n", + "\n", + "![plot example](../../../images/challenge/fates_bart_gpp_lai.png)\n", + "\n", + "*

Figure: GPP and LAI climatology for our inventory-initialized case.

*\n", + " \n", + "
\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Questions: \n", + "\n", + " - How do GPP and LAI track each other through the year? Which leads the other in spring, and does that make physical sense?\n", + " - Does LAI fall all the way to zero in winter? What does the residual tell you about the mix of plant functional types at this site? (think back to which PFTs grow at BART.)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---\n", + "## 3. What's growing at BART? PFT composition\n", + "\n", + "Next, let's 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": [ + "### 3.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/challenge/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": [ + "### 3.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/challenge/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", + "## 4. Anatomy of a FATES forest\n", + "\n", + "### 4.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/challenge/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": [ + "### 4.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/challenge/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": [ + "### 4.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/challenge/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/challenge/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", + "## 5. Demographic Processes\n", + "\n", + "### 5.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/challenge/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": [ + "### 5.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/challenge/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/challenge/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/challenge/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", + "## 6. 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", + "### 6.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/challenge/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": [ + "### 6.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/challenge/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", + "## 7. 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 +}