Nutrient controls on leaf level photosynthetic processes #768
Conversation
| ! kn = 0.11. Here, derive kn from vcmax25 as in Lloyd et al | ||
| ! (2010) Biogeosciences, 7, 1833-1859 | ||
|
|
||
| kn = decay_coeff_kn(ft,currentCohort%vcmax25top) |
There was a problem hiding this comment.
or, alternative to my above comment, should we remove these and pass nscaler_prt into LeafLayerMaintenanceRespiration() ?
There was a problem hiding this comment.
I'd be a little cautious about creating these extra _prt parameters that are conceptually the same trait, just calculated differently. We had some similar discussion wrt gs traits 616 and decided to not proliferate variables when they are the same. We stuck with a single g0 as it's conceptually the same in BB and Medlyn models but different g1 because the two models mean that g1 has a different effect and so is conceptually slightly different.
I would suggest doing that for all variables labelled _prt, but especially these nscalar and kn given that they are essentially calculated in the same way. Couldn't we have a case statement for the calculation of kn (would be unnecessary if both Vcmax's were labeled the same way) and then just a single calculation of nscalar ?
| vcmax25top_prt = min(max(vcmax25top_prt, 10.0_r8), 150.0_r8) | ||
| jmax25top_prt = min(max(jmax25top_prt, 10.0_r8), 250.0_r8) | ||
|
|
||
| kn_prt = decay_coeff_kn(ft,vcmax25top_prt) |
There was a problem hiding this comment.
are the kn_prt and nscaler_prt variables used @jenniferholm ? Should we remove these?
|
@jenniferholm I noticed that there the leaf maintenance respiration calculations use the old or existing nscalar value, which is based on the stochiometry parameter constant: But the effective biophysical rates, vcmax, jmax etc, take arguments for both the old (nscalar) and new (nscalar_prt) values, which gets confusing. I'm wondering if it would be worth it to determine the lnc, kn and nscalar values before both subroutines (which you do), but only arrive at a singluar lnc, nscalar and kn value (not both nscalar and nscalar_prt, and kn/kn_prt, etc), and pass in those just the single value to the routines, consistently. Likewise, I think we need to add a control flag from the parameter file that allows us to override the prognostic N and P concentrations in favor of the values derived from the stoichiometry parameter constants. Sorry if I just don't get it, maybe I just dont understand why a different nscalar value would be used for maintenance respiration and the vcmax. |
There was a problem hiding this comment.
Looks good. I see no issues that will cause problems and all of the new calculations look good. See my general comments about using the _prt suffix for variables that are conceptually the same. This code could be significantly simplified imo if the suffix was removed and the newly calculated vcmax etc values added to the cohort lists / data structure.
|
|
||
| fates_eca_vcmax_np4 = 0.282, 0.282, 0.282, 0.282, 0.282, 0.282, 0.282, 0.282, | ||
| 0.282, 0.282, 0.282, 0.282 ; | ||
|
|
There was a problem hiding this comment.
How was variability in these parameters calculated? Also worth noting based on analysis by Luo et al 2021 that the values of Vcmax that we get are v likely to be low-biased based on the parameter values from my 2014 paper.
This will also low-bias Jmax, though if we can get better estimates of Vcmax then I don't think we'll need to edit the Jmax~Vcmax relationship. Am working on getting better estimates of the Vcmax parameters.
| ! kn = 0.11. Here, derive kn from vcmax25 as in Lloyd et al | ||
| ! (2010) Biogeosciences, 7, 1833-1859 | ||
|
|
||
| kn = decay_coeff_kn(ft,currentCohort%vcmax25top) |
There was a problem hiding this comment.
I'd be a little cautious about creating these extra _prt parameters that are conceptually the same trait, just calculated differently. We had some similar discussion wrt gs traits 616 and decided to not proliferate variables when they are the same. We stuck with a single g0 as it's conceptually the same in BB and Medlyn models but different g1 because the two models mean that g1 has a different effect and so is conceptually slightly different.
I would suggest doing that for all variables labelled _prt, but especially these nscalar and kn given that they are essentially calculated in the same way. Couldn't we have a case statement for the calculation of kn (would be unnecessary if both Vcmax's were labeled the same way) and then just a single calculation of nscalar ?
|
|
||
| end select | ||
|
|
||
|
|
There was a problem hiding this comment.
Similar comment to above here. This could all be simplified if the _prt suffix was removed.
| vcmax_np3 => EDPftvarcon_inst%vcmax_np3, & !vcmax~np relationship coefficient | ||
| vcmax_np4 => EDPftvarcon_inst%vcmax_np4 ) !vcmax~np relationship coefficient | ||
|
|
||
| ! vcmax25 parameters from Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593 |
There was a problem hiding this comment.
It will be cleaner if fnr, act25, jmax_np parameters are all moved to parameter file, just like vcmax_np parameters. It will also benefit parameter tuning later.
| vcmax25top_prt = min(max(vcmax25top_prt, 10.0_r8), 150.0_r8) | ||
| jmax25top_prt = min(max(jmax25top_prt, 10.0_r8), 250.0_r8) | ||
|
|
||
| kn_prt = decay_coeff_kn(ft,vcmax25top_prt) |
There was a problem hiding this comment.
Do you consider multi-layer canopy here? If so, the lnc_top, lpc_top calculation needs to consider the vertical distribution of leaf N by applying nscaler_prt term. If not, I think the nscaler_prt should be always one. @jenniferholm
…esis, fixed typo for vcmax_np1
There was a problem hiding this comment.
Thanks @jenniferholm this looks great and will be a really useful option to be able to use. I made a bunch of comments where there are things like inline numbers that don't have super clear provenance or units. Also it looks like a bunch of unrelated changes are in the parameter file, I think once you update relative to the main branch and resolve those conflicts, most of those will go away anyway.
| fnr = 7.16_r8 | ||
| act25 = 3.6_r8 !umol/mgRubisco/min | ||
| ! Convert rubisco activity units from umol/mgRubisco/min -> umol/gRubisco/s | ||
| act25 = act25 * 1000.0_r8 / 60.0_r8 | ||
|
|
||
| ! jmax~np relationship coefficients needed in the parteh calculation of | ||
| ! variable jmax25top_prt | ||
| jmax_np1 = 1.246_r8 | ||
| jmax_np2 = 0.886_r8 | ||
| jmax_np3 = 0.089_r8 |
There was a problem hiding this comment.
should all these numbers be in a parameter file rather than inline fortran? or at least with a clearer provenance to be able to track, as well as units etc.
| vcmax25top_prt = min(max(vcmax25top_prt, 10.0_r8), 150.0_r8) | ||
| jmax25top_prt = min(max(jmax25top_prt, 10.0_r8), 250.0_r8) |
There was a problem hiding this comment.
what is the source of these numbers?
| ! vcmax25top_prt_ft from above. | ||
| ! Approach taken from ELM photosynthesis code (Q. Zhu) | ||
| ! dayl_factor is a photoperiod acclimation correction term from Bauerle et al. 2012 | ||
| jmax25top_prt_c = (2.59_r8 - 0.035_r8*min(max((temp_a10-tfrz),11._r8),35._r8)) * & |
There was a problem hiding this comment.
where are these numbers from?
| real(r8) :: jmax_np1 ! jmax~np relationship coefficient | ||
| real(r8) :: jmax_np2 ! jmax~np relationship coefficient | ||
| real(r8) :: jmax_np3 ! jmax~np relationship coefficient |
There was a problem hiding this comment.
can we add units for each of these?
| char fates_hydr_organname_node(fates_hydr_organs, fates_string_length) ; | ||
| fates_hydr_organname_node:units = "unitless - string" ; | ||
| fates_hydr_organname_node:long_name = "Name of plant hydraulics organs (DONT CHANGE, order matches media list in FatesHydraulicsMemMod.F90)" ; | ||
| char fates_litterclass_name(fates_litterclass, fates_string_length) ; | ||
| fates_litterclass_name:units = "unitless - string" ; | ||
| fates_litterclass_name:long_name = "Name of the litter classes, for variables associated with dimension fates_litterclass" ; |
There was a problem hiding this comment.
should we pull these lines out of this PR as they pertain to a different one?
| char fates_hydr_organname_node(fates_hydr_organs, fates_string_length) ; | |
| fates_hydr_organname_node:units = "unitless - string" ; | |
| fates_hydr_organname_node:long_name = "Name of plant hydraulics organs (DONT CHANGE, order matches media list in FatesHydraulicsMemMod.F90)" ; | |
| char fates_litterclass_name(fates_litterclass, fates_string_length) ; | |
| fates_litterclass_name:units = "unitless - string" ; | |
| fates_litterclass_name:long_name = "Name of the litter classes, for variables associated with dimension fates_litterclass" ; |
There was a problem hiding this comment.
Lets see what happens when we merge up to master @ckoven , either way, we will preserve the code on the master side where needed. I see that there is a conflict on this file anyway, so these issues will probably get blow away at that time.
| double fates_hydr_htftype_node(fates_hydr_organs) ; | ||
| fates_hydr_htftype_node:units = "unitless" ; | ||
| fates_hydr_htftype_node:long_name = "Switch that defines the hydraulic transfer functions for each organ." ; | ||
| fates_hydr_htftype_node:possible_values = "1: Christofferson et al. 2016 (TFS); 2: Van Genuchten 1980" ; | ||
| double fates_prt_organ_id(fates_prt_organs) ; | ||
| fates_prt_organ_id:units = "index, unitless" ; | ||
| fates_prt_organ_id:long_name = "This is the global index the organ in this file is associated with in PRTGenericMod.F90" ; |
There was a problem hiding this comment.
should we pull these lines out of this PR as they pertain to a different one?
| double fates_hydr_htftype_node(fates_hydr_organs) ; | |
| fates_hydr_htftype_node:units = "unitless" ; | |
| fates_hydr_htftype_node:long_name = "Switch that defines the hydraulic transfer functions for each organ." ; | |
| fates_hydr_htftype_node:possible_values = "1: Christofferson et al. 2016 (TFS); 2: Van Genuchten 1980" ; | |
| double fates_prt_organ_id(fates_prt_organs) ; | |
| fates_prt_organ_id:units = "index, unitless" ; | |
| fates_prt_organ_id:long_name = "This is the global index the organ in this file is associated with in PRTGenericMod.F90" ; |
| @@ -11,6 +11,7 @@ dimensions: | |||
| fates_pft = 12 ; | |||
| fates_prt_organs = 4 ; | |||
| fates_string_length = 60 ; | |||
| fates_hlm_pftno = 14 ; | |||
There was a problem hiding this comment.
should we pull these lines out of this PR as they pertain to a different one?
| @@ -11,6 +11,7 @@ dimensions: | |||
| fates_pft = 12 ; | |||
| fates_prt_organs = 4 ; | |||
| fates_string_length = 60 ; | |||
| fates_hlm_pftno = 14 ; | |||
There was a problem hiding this comment.
| fates_hlm_pftno = 14 ; |
rgknox
added
the
PR status: Not Ready
This code takes into account leaf-level nutrient effects on photosynthesis when PARTEH is used, particularly updating how the vcmax25 and jmax25 terms are calculated. These new relationships between leaf N, leaf P, and regression coefficients estimated from a photosynthesis data set will allow for a dynamic approach to represent nutrient constraints on photosynthesis.
Connected to the issue discussion at #597 , #443 , and could later be important for the acclimation work being done by others.
Description:
Following the work that was implemented in ELMv1 by Qing Zhu 2019 (https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2018MS001571), the code has been updated so that "photosynthesis is a function of leaf N (lnc: gN/m2 leaf area) and leaf P (lpc: gP/m2 leaf area) contents, which are commonly measured leaf properties, and their fundamental controls on plant photosynthesis rates are widely supported by observations (Kattge et al., 2009; Reich et al., 2009)."
The equations for vcmax25top and jmax25top have now been updated to be a function of leaf N and leaf P, based on the results from (Walker et al., 2014), instead of the static values provided in the parameter file.
These changes only occur when PARTEH is turned in (in c-only mode and CNP mode in PARTEH). The default, non-PARTEH mode, still uses the vcmax25top from the parameter file, and jmax25top being derived from the parameter file vcmax25top within "FatesParameterDerivedMod".
Now that there is a varying vcmax25top term, additional changes need to be made where other terms are based on vcmax25top. For example, the "kn" term (the decay coefficient) is now based on the newly updated vcmax25top, which in turn produces a new "nscaler" term (leaf nitrogen scaling coefficient). Two new variables were introduced to keep track of these new terms: kn_prt, and nscaler_prt (standing for PARTEH). For now a newly updated "kn_prt" term is only being implemented in the photosynthesis code when using PARTEH, not in other parts of the code where "kn" is used, such as in FatesAllometryMod when using kn to calculate things like LAI an SLA. The default, parameter value of vcmax25top is still being used.
The "kp25top" term (initial slope of CO2 response curve) has also been updated to account for the varying vcmax25top, but only when using PARTEH. Otherwise kp25top is still being derived in "FatesParameterDerivedMod" using this%kp25top(ft,iage) = 20000._r8 * vcmax25top(ft,iage).
Note - I would appreciate any feedback on this change, and also make sure I'm not double charging kp25top * 20000 twice. In the LEAFN figure below the initial jump in leaf N in the "w/ photosyn. CO2 response" test is due to this sudden change in how kp25top is calculated, which wasn't used during the spinup procedure.
Fixes #597
Collaborators:
@rgknox , @walkeranthonyp , @ckoven , @Qing Zhu, @rosiealice , @alistairrogers
Expectation of Answer Changes:
These code changes will definitely change answers. Below are some output examples from different scenarios testing the changes to different forest states and fluxes.
All simulations began after running an AD spinup case (with supplemental N and P), a 200-year regular spinup (with AD turned off, and no supplemental N and P, until all carbon pools were stable in equilibrium). Then code changes to photosynthesis were turned on.
Legend key for figures below:
"Regular Spinup" - continuing the regular spinup mode, with no leaf level photosynthesis changes (as a pseudo default case)
"w/ photosyn." = all leaf level photosynthesis code changes, including adding the multiplier of "day length factor", which allows for the consideration of light limitation. Testing at a boreal forest site, so large impacts.
"w/ photosyn. no day-length" = removing the day length multiplying factor to test the impacts of only the leaf level nutrient additions.
"w/ photosyn. CO2 response" = all leaf level photosynthesis code changes, including adding the multiplier of day length factor, and slope of CO2 response curve based on updated vcmax25top.
"C-only PARTEH" = test of if case is (prt_carbon_allom_hyp), and using the new approach of vcmax25top_prt = lnc_top * flnr(ft) * fnr * act25.
Checklist:
Yes, the code changes requires these checklist items and they still need to be done.
Test Results:
CTSM (or) E3SM (specify which) test hash-tag:
CTSM (or) E3SM (specify which) baseline hash-tag:
E3SM = Ryan's PARTEH branch = "rgknox/lnd/fates-cnp cdbeea0 Made modifications to supplementation of NP during FATES-ECA."
FATES baseline hash-tag:
Test Output:
Single site test, in Alaska boreal forest -