Modified crown allometric spread logic #279
Conversation
|
Adding some comparison figures to show the dynamics of this code. In each, black line is control case, red line is code from this PR. They're all from a 150 year run of 1x1amazon in prescribed physiology mode. first is the mean canopy area divided by patch area, conditional on patch age distribution, for years 50-150: second is a biomass time series for the 150 years: third is a time series of the total tree coverage fraction across the site: I looked at some other variables (size distributions, etc) but nothing stood out in particular. To me these results make sense. The change achieves the goal of getting rid of the stepwise canopy closure pattern that you can see in the black line of the first figure. The overall growth is slower, which makes sense given that canopies are smaller. It also makes sense that more bare ground is evident since it now takes longer for the canopy to close in disturbed patches. I'm not totally sure what's going on in the initial oscillation in the canopy closure time series, but suspect it has to do with patch dynamics during the first few years. |
|
Ran tests, results below. This code actually passes regression tests, which underscores that our current testing regime isn't strict enough, but that will be fixed shortly by https://github.com/NGEET/fates-clm/pull/23 @rgknox, could you run this through your science test suite? |
|
I forgot to add biomass to the test list, but here is a comparison with master: |
biogeochem/EDGrowthFunctionsMod.F90
Outdated
| spreadterm = cohort_in%siteptr%spread * EDPftvarcon_inst%allom_d2ca_coefficient_max(cohort_in%pft) + & | ||
| (1._r8 - cohort_in%siteptr%spread) * EDPftvarcon_inst%allom_d2ca_coefficient_min(cohort_in%pft) | ||
| ! | ||
| c_area = 3.142_r8 * cohort_in%n * (spreadterm * dbh)**crown_area_to_dbh_exponent |
There was a problem hiding this comment.
its kinda like pi but only has 4 significant figures. happy to switch to actual pi now that i've confirmed that it's bit for bit as it currently stands.
There was a problem hiding this comment.
interesting, that's the circle thing
There was a problem hiding this comment.
which is itself slightly weird. I feel like that made more sense physically when the crown_area_to_dbh_exponent was 2 (as in the original PPA papers) but at this point its sort of an artifact. that term could easily be pulled into the allom_d2ca_coefficient but for simplicity I left as-is since doing so would change the dependency between allom_d2ca_coefficient and crown_area_to_dbh_exponent.
There was a problem hiding this comment.
ok, so then would people prefer I subsume it into the default value of the allom_d2ca_coefficient given whatever default crown_area_to_dbh_exponent we are currently using? or leave as-is and make it a named pi?
There was a problem hiding this comment.
well, honestly my comment started off as a joke, but now I see where you are going. There is no math reflecting a circular form anymore, so the pi is somewhat misleading. imho, the precision on that pi specifically isn't important and everyone knows what 3.14 is. If you want to remove it and subsume in the parameter, up to you, I would see it as a bonus prize.
There was a problem hiding this comment.
actually, while doing this, i realize that its weird that this particular power law is written as (ax)^b whereas all the others are just ax^b. does anyone mind if I pull the a out of the parentheses to conform to what the other pieces of code are doing (and adjust the value accordingly)?
There was a problem hiding this comment.
that comment got mangled by formatting, look at edited version
There was a problem hiding this comment.
thanks. that makes it easier to understand the crown plasticity as just the ratio of the max:min terms rather than that ratio raised to some power.
| @@ -2624,6 +2631,12 @@ subroutine define_history_vars(this, initialize_variables) | |||
| avgflag='A', vtype=site_age_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & | |||
| ivar=ivar, initialize=initialize_variables, index = ih_zstar_si_age ) | |||
|
|
|||
| call this%set_history_var(vname='BIOMASS_BY_AGE', units='m', & | |||
There was a problem hiding this comment.
This diagnostic is a nice addition, thanks for including this @ckoven
|
late to this party. 1) yes, it's meant to be PI. Not sure where the 4sf thing came from. 2) I guess if the exponent isn't 2 then the pi r^2 thing doesn't really apply anymore, and so maybe it is now arbitrary. The change from 2 was made to prevent the tree_lai constantly changing as the trees got older given that the dbh-> leaf exponent is <2.... |
biogeochem/EDCanopyStructureMod.F90
Outdated
| if(EDPftvarcon_inst%woody(currentCohort%pft) == 1)then | ||
| arealayer(currentCohort%canopy_layer) = arealayer(currentCohort%canopy_layer) + currentCohort%c_area | ||
| if(EDPftvarcon_inst%woody(currentCohort%pft) .eq. 1 .and. currentCohort%canopy_layer .eq. 1 ) then | ||
| arealayer = arealayer + currentCohort%c_area | ||
| endif | ||
| currentCohort => currentCohort%shorter | ||
| enddo | ||
|
|
There was a problem hiding this comment.
So this is not really arealayer, but area of all upper canopy layers in the site.
There was a problem hiding this comment.
right, i guess i ought to rename it
|
I guess the former, since if we ever use some real data to calibrate this
it won't necessarily conform to the term with PI in...
…On 2 October 2017 at 15:22, Charlie Koven ***@***.***> wrote:
***@***.**** commented on this pull request.
------------------------------
In biogeochem/EDGrowthFunctionsMod.F90
<#279 (comment)>:
> @@ -260,15 +258,12 @@ real(r8) function c_area( cohort_in )
! So we allow layer index exceedence here and force it down to max.
! (rgk/cdk 05/2017)
! ----------------------------------------------------------------------------------
-
- can_layer_index = min(cohort_in%canopy_layer,nclmax)
-
- if(EDPftvarcon_inst%woody(cohort_in%pft) == 1)then
- c_area = 3.142_r8 * cohort_in%n * &
- (cohort_in%patchptr%spread(can_layer_index)*dbh)**crown_area_to_dbh_exponent
- else
- c_area = 3.142_r8 * cohort_in%n * (ED_val_grass_spread*dbh)**crown_area_to_dbh_exponent
- end if
+
+ ! apply site-level spread elasticity to the cohort crown allometry term
+ spreadterm = cohort_in%siteptr%spread * EDPftvarcon_inst%allom_d2ca_coefficient_max(cohort_in%pft) + &
+ (1._r8 - cohort_in%siteptr%spread) * EDPftvarcon_inst%allom_d2ca_coefficient_min(cohort_in%pft)
+ !
+ c_area = 3.142_r8 * cohort_in%n * (spreadterm * dbh)**crown_area_to_dbh_exponent
ok, so then would people prefer I subsume it into the default value of the
allom_d2ca_coefficient given whatever default crown_area_to_dbh_exponent we
are currently using? or leave as-is and make it a named pi?
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
<#279 (comment)>, or mute
the thread
<https://github.com/notifications/unsubscribe-auth/AMWsQ4hTw-afYzWI7t8tCZ1DDpw12J-2ks5soVQUgaJpZM4PmDxX>
.
--
-----------------------------------------------------------------
Dr Rosie A. Fisher
Staff Scientist
Terrestrial Sciences Section
Climate and Global Dynamics
National Center for Atmospheric Research
1850 Table Mesa Drive
Boulder, Colorado, 80305
USA.
+1 303-497-1706
http://www.cgd.ucar.edu/staff/rfisher/
|
|
ok, will do. |
… it into the allom_d2ca_coefficient parameters
|
ok, for future reference, our default value of the fortran variable dbh2bl_b (which in the parameter file is called fates_allom_d2bl2) is 1.3; therefore the number that I'm multiplying the allom_d2ca_coefficient_max and allom_d2ca_coefficient_min terms by is 3.142 ^ 0.769230769230769 = 2.412494282999975. |
|
@ckoven , I'm curious, what does the site level spread parameter projection look like vis-a-vis the other plots you show? As far as a technical review goes, my impression is that the code seems to be doing what the intention is. I'm still grappling with the conceptual change though. Do you think it is likely that the first group of cohorts to exist in a site will still go through the de-spreading phase when they first start to fill out their layer(shown by the hump in your plots, for the old method)? And after that since they have already established a crowded site, that is why we don't see the saddle-point? Also, would you agree that it is expected that a closed-canopy forest should find an equlibrium around spread = 0? I think the thing I am grappling with is that I'm just not sure what this site level spread means anymore. At some point, I think we would expect spread to reach equlibrium and then I wonder why this spread isn't subsumed in the allometric parameters. My other thought, is that spreading is about individuals reacting (somehow) to the space that is afforded them. And now we have individuals that are reacting, but not really based on their own environment. I like the simplicity we gain from this improvement in spread, but now I wonder if its at the site level, why have it at all? I suppose if someone wanted to turn it off, the idea would be to just set the threshold extremely high and no de-spreading would ever trigger. |
|
@rgknox that's right that it starts off at one and drops quickly to zero once the canopy starts to close. figure shown below. To the rest of your questions, I think the point of this is that it should be subsumed into the allometric parameters, except for the special case (which is what I believe the problem that @rosiealice was intending to solve was) where we may want there to be plasticity with respect to the environment in thinking about forest to savanna transition zones. the within-disturbance transition zones ought to be part of the normal trajectory of an individual's growth rather than the stepwise logic that switches from one allometric trajectory to another (and possibly back and forth) depending on the patch-level light environment. so for closed-canopy forests, I think we want to work under the assumption that spread = 0 and only the fates_allom_d2ca_coefficient_min term is relevant. that way there is just the one trajectory which can be diagnosed from observations.
|
|
Thoughts on this
1) If we started off with a closed canopy, consisting of equal parts of
each PFT, rather than bare ground, then we wouldn't necessarily see the
'squashing' logic when the canopy closed. That, however, would create
something of a epic tragedy of the commons in semi-arid places, since all
the members of the identical large cohort(s) would run out of carbon at the
same time. (I'm seeing something like this in the grass simulations I'm
running at the moment. Which I was on the verge of reporting on and now
cheyenne is down for the day). I think in general it's better to start off
with a smaller init_d, for the sake of getting greater diversity in height
by the time canopy closure happens, therefore making dominance by a single
PFT/cohort less likely.
2) Remind me why it goes to zero? Is that a scalar from the maximum to the
minimum spread limit?
3) In the sense that we always do model spin-up, and we expect weird things
to happen while that is going on, do we need to worry if we see the 'step'
in the first ten years of the simulation when ultimately we will throw that
away?
…On 2 October 2017 at 16:00, Charlie Koven ***@***.***> wrote:
@rgknox <https://github.com/rgknox> that's right that it starts off at
one and drops quickly to zero once the canopy starts to close. figure shown
below.
To the rest of your questions, I think the point of this is that it
*should* be subsumed into the allometric parameters, except for the
special case (which is what I believe the problem that @rosiealice
<https://github.com/rosiealice> was intending to solve was) where we may
want there to be plasticity with respect to the environment in thinking
about forest to savanna transition zones. the within-disturbance transition
zones ought to be part of the normal trajectory of an individual's growth
rather than the stepwise logic that switches from one allometric trajectory
to another (and possibly back and forth) depending on the patch-level light
environment. so for closed-canopy forests, I think we want to work under
the assumption that spread = 0 and only the fates_allom_d2ca_coefficient_min
term is relevant. that way there is just the one trajectory which can be
diagnosed from observations.
[image: screen shot 2017-10-02 at 2 54 01 pm]
<https://user-images.githubusercontent.com/10852790/31101227-986f5f22-a781-11e7-9142-01428eb4b297.png>
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
<#279 (comment)>, or mute
the thread
<https://github.com/notifications/unsubscribe-auth/AMWsQ_8OTkR4nhHjlObxE5k5SZITU-bTks5soV0agaJpZM4PmDxX>
.
--
-----------------------------------------------------------------
Dr Rosie A. Fisher
Staff Scientist
Terrestrial Sciences Section
Climate and Global Dynamics
National Center for Atmospheric Research
1850 Table Mesa Drive
Boulder, Colorado, 80305
USA.
+1 303-497-1706
http://www.cgd.ucar.edu/staff/rfisher/
|
|
Hi Rosie- (1) I agree with the idea that it's better to have relatively low initd and a heterogeneous size structure than have a massive initial cohort. (2) yes, the scalar now goes between 0 and 1, with 0 being the PFT-specific minimum crown allometry and 1 being the maximum. (3) agreed too that we should expect various sorts of weirdness during spinup, and that the initial period of crown allometric narrowing is ok here since, unlike before, it doesn't affect the long-term steady-state dynamics. |
…changed the name of arealayer var in canopy_spread subroutine
|
No, I don't mind, if we are in the process of changing the meaning of 'a'
anyways.
…On 4 October 2017 at 12:27, Charlie Koven ***@***.***> wrote:
***@***.**** commented on this pull request.
------------------------------
In biogeochem/EDGrowthFunctionsMod.F90
<#279 (comment)>:
> @@ -260,15 +258,12 @@ real(r8) function c_area( cohort_in )
! So we allow layer index exceedence here and force it down to max.
! (rgk/cdk 05/2017)
! ----------------------------------------------------------------------------------
-
- can_layer_index = min(cohort_in%canopy_layer,nclmax)
-
- if(EDPftvarcon_inst%woody(cohort_in%pft) == 1)then
- c_area = 3.142_r8 * cohort_in%n * &
- (cohort_in%patchptr%spread(can_layer_index)*dbh)**crown_area_to_dbh_exponent
- else
- c_area = 3.142_r8 * cohort_in%n * (ED_val_grass_spread*dbh)**crown_area_to_dbh_exponent
- end if
+
+ ! apply site-level spread elasticity to the cohort crown allometry term
+ spreadterm = cohort_in%siteptr%spread * EDPftvarcon_inst%allom_d2ca_coefficient_max(cohort_in%pft) + &
+ (1._r8 - cohort_in%siteptr%spread) * EDPftvarcon_inst%allom_d2ca_coefficient_min(cohort_in%pft)
+ !
+ c_area = 3.142_r8 * cohort_in%n * (spreadterm * dbh)**crown_area_to_dbh_exponent
that comment got mangled by formatting, look at edited version
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
<#279 (comment)>, or mute
the thread
<https://github.com/notifications/unsubscribe-auth/AMWsQwJkvNBB2GhygMFg9oBs4GUMiFU-ks5so84igaJpZM4PmDxX>
.
--
-----------------------------------------------------------------
Dr Rosie A. Fisher
Staff Scientist
Terrestrial Sciences Section
Climate and Global Dynamics
National Center for Atmospheric Research
1850 Table Mesa Drive
Boulder, Colorado, 80305
USA.
+1 303-497-1706
http://www.cgd.ucar.edu/staff/rfisher/
|
|
As for signing off on the conceptual approach and overall goals of this PR, I'm ok with it. I think there are some minor syntax changes requested above, but after that, I will approve. |
|
Updated code to make the canopy coefficient parameter include the pi term and also the exponent. ran test suite: it now fails the baseline, as expected because raising the spread terms (which were 0.18 and 0.3) to the 1.3 power results in a sufficiently-long decimal representation that it is no longer roundoff. results below. I looked at bunch of variables from a long run and everything stays within a percent or so, and without a systematic shift, relative to the earlier versions of this PR (though they differ by slightly more than that from the master, and with a systematic shift towards less-full canopies at early patch ages, as discussed above). if this looks ok to @rosiealice and @rgknox now, please feel free to pull. |
|
@ckoven , can you confirm that tests were ran with fates-clm-045be03 and this branch against master? |
|
yep, they were fates-clm 045be03 and fates 06e3c88, and run against master. fails on baseline and namelist comparisons only, as expected. echo.... |
|
approved |
This PR addresses two issues related to the canopy spread logic: it moves the calculation from the patch to the site level (#274) and it also makes the spread parameters PFT-level rather than global parameters (#199). The upshot of this is that there is still within-indvidual-lifetime-plasticity in the crown allometry, but it is only allowed to express itself over large-scale spatial gradients (e.g. forest to savanna) and, at some point when we have multiple fates sites (columns) per grid cell, along gradients of heterogeneity in that dimension. It does not allow crown allometric plasticity along fine-scale disturbance gradients. However, it now opens the door for trait-filtered allometric variability along any of these gradients.
Changes to user interface are a few new parameters in the parameter file, and a few old ones are deleted:
3 global parameters: fates_grass_spread, fates_maxspread, and fates_minspread have been deleted and replaced with 2 PFT parameters: fates_allom_d2ca_coefficient_max and fates_allom_d2ca_coefficient_min. If these are the same for a given pft (e.g. grasses) then there is no crown allometric plasticity in that PFT.
1 new global parameter is added: fates_canopy_closure_thresh; this sets the boundary between where plants will have wide or narrow crown allometries. This was previously an inline hard-coded parameter set to 0.9 (at the patch level), it is now a specified parameter with (for now anyway) a default value of 0.8 since we expect less canopy closure at the site level than at the patch level.
Another change, consistent with the desire for this to only apply to savanna-to-forest gradients rather than disturbance gradients is that the spread term is applied only in the uppermost canopy level; the understory trees are assumed to have the same spread as their canopy neighbors. This avoids the logical issue of trees instantly adjusting their shape when promoted or demoted, but creates the practical issue of understory trees having small light-gathering footprints. I haven't worked out how big a deal this is yet, and so invite discussion on whether this seems like a good change; in any case one proposed solution is to make the crown allometry vary between successional PFTs so that in practice the understory has a different allometry than the canopy, but via a trait-filtering rather than individual-plasticity mechanism.
This may have an effect on the grass issues discussed in #277, at least in that the crown spread specification for grasses happens via the PFT rather than the global parameter.
For now, hold off on merging, as I have only run qualitative tests rather than the full test suite. I'll post figures of changes as they are done. This will be an answer-changing branch once incorporated. Mainly I wanted to throw this up as a PR to invite people (particularly @rosiealice and @rgknox) to look at the code and discuss before I get too far with it.