Medlyn stomatal conductance error check fails #1398
Comments
|
On a related note, I was cross checking the code and CLM5 technote and saw that term b in equation 9.25 of the technote seems to be missing a minus sign out front, which exists in the code. Removing the minus sign gives negative conductance errors so I assume its inclusion is correct. But it might be worth checking the entire implementation of that quadratic with the source/author of the code. |
|
Would we expect the operational stomatal conductance to match the Medlyn equation in this way? Whether this is the best approach is something I'm certainly interested in exploring, but I'm not sure this is a bug per se. I can edit the the technical note, if this is not documented well. |
|
Thanks for this note @djk2120. I had wondered how PHS was influencing this problem but hadn't gotten a chance to look into it yet. Given that the error check implemented by @heroldn is for maximum conductance and PHS downregulates this value, it is perhaps not surprising that the values don't match. Prior to PHS implementation, conductance was downregulated through the effect of btran on Vcmax (and also its effect on g0), so water stress was included more directly in the conductance equation. I see that it is mentioned in the PHS chapter of the tech note, but I think it might be a good idea for us to update the stomatal conductance chapter to mention that the Medlyn conductance is downregulated by PHS. |
I was logged into a different GitHub account when I wrote this. |
|
Thanks Daniel, that seems like a reasonable explanation. The stomatal conductance from the error check (the equation 2.9.1) always seems to be greater than the operational conductance in my single-point simulation so that seems to be consistent with your explanation. It might be good to confirm that they agree in a case with no water stress? |
|
Actually my values using eq. 9.1 seemed consistently lower than what CLM gave. And they still weren't matching even when PHS was off. But that was using the edits you made @djk2120 manually transferred into my own code, so it's possible I got something wrong when I combined. But still doesn't explain my lower values when using eq. 9.1. |
|
Digging a little further (with a healthy injection of print statements) I notice that the different gs values given by CLM (compared to eq. 9.1) are already present from the quadratic formula given by eq. 9.24 (and implemented in ci_func_PHS/ci_func). Should this make us suspect the quadratic formula? The technote says that the quadratic is arrived at by substituting eq. 9.23 for es into eq. 9.1, but I don't see where es factors into eq. 9.1, or how eq. 9.23 could be rearranged to substitute for something besides es. @djk2120, @danicalombardozzi , do you know who derived the numerical implementation? A colleague here also noted that the technote doesn't explicitly define a for the quadratic (implicitly it's "1"), and I note a typo or two in this section as well, so I wonder if something got missed and I'd be keen to confirm these equations with the author. Maybe it's nothing but still. |
|
I'm just copying in an offline conversation that happened over the last week or so, so that there is a record of it here. Hi Rosie and Gordon Thanks for this. And thanks Gordon your book looks like a great reference and it looks like it's resolved the issue too, but I'd appreciate all your eyes on this to confirm, and there are definitely some inconsistencies remaining (see below). Apologies for the long email. Basically, the analytical version of Medlyn Gs, according to Gordon's book eq. 12.21, should use the leaf surface VPD. The VPD in the source code and that I've used in my check is the leaf-to-air VPD. The technote could use revising in this regard because it uses the same symbol D to refer to both types of VPD (but I see it does say in paragraph form just under eq. 9.1 that that equation uses leaf surface VPD, but people like me just look at symbols!). Once I implemented leaf surface VPD in my check the values from the analytical equation matched much more closely with what comes from the quadratic form, often to within a few decimal places, but occasionally it's still extremely different too (attached is grepped output from a 5 day F case that I've been eye-balling). Most errors seem < 50 and thus represent a small % error, since Gs is usually on the order of tens of thousands. So, IF there was going to be an internal check in the code that compares the analytical and numerical solutions (as there is for BB), the error tolerance would have to be much higher than that used for BB (which is 0.1) and maybe % based. It would still fail occasionally as I see extremely large differences once in a while (factor of 2 or 3 different and Gs values can go into the millions, but these look rare). Not sure what the tolerance in any check should be but of course your intuitions are better than mine. Now, some small inconsistencies to note/clarify/revise in future tech notes: Thanks for everyone's input. If we can confirm this issue is resolved (am I correct that the occasional large discrepancy is OK?), then I can add the above details to the github report. Cheers Nick. Indeed, the quadratic solver predates the Medlyn implementation, and Gordon's book is probably your best resource on that... FATES isn't just for Paleo! We now have a 'satellite phenology' version just for looking into canopy physiology issues like this too ;) Keep us in the loop with what you find. -Rosie Le mer. 14 juil. 2021 à 18:24, Gordon Bonan [email protected] a écrit : Look at my book Climate Change and Terrestrial Ecosystem Modeling (Cambridge University Press, 2019). Chapter 12 on stomatal conductance gives the derivation and numerical solution needed for the Ball-Berry and Medlyn models (see specifically eqs. 12.14-12.23). There is accompanying matlab code (explained in the book). The code can be found at https://github.com/gbonan/bonanmodeling Gordon On Wed, Jul 14, 2021 at 6:00 AM Nicholas Herold [email protected] wrote: Thanks very much for your reply. I can understand if this is stretching your memory so really appreciate your thoughts on it. What you say about the technote consisting of old text from the BB method makes sense from my comparison of the different versions of that document. I actually have been looking at the code and there is actually a github issue logged for this here. I didn't include it initially to try and keep my email from being an information dump! If the es equation in the technote is redundant then I guess my question would be how does one arrive at the quadratic equation used in the code to get Gs (line 4159 of your link and described in eq. 9.24-9.26)? It's this quadratic that's giving different values compared to the Medlyn equation (eq. 9.1). You can see from the github comments Daniel suggested it could be because PHS downregulates Gs, but if I understand the code and methods correctly the Gs coming out of the quadratic should be prior to any downregulating and thus should match eq. 9.1 (I'm not a land modeller so have been learning many a new concept during all this!). I also tested this without PHS and still get inconsistent values, so I'm ruling out PHS at the moment, and that's led me to the quadratic getting used. Do you have any recollection of how this is derived by any chance? Thanks for the heads up about FATES too, I wasn't aware. In my paleo CESM days this would have been right up my alley (I was using the immensely simple BIOME4!). I'll definitely keep this in mind. Cheers Nick. -actually- cc-ing you guys! ---------- Forwarded message --------- Dear Nick, Thanks for this. On your first question (did I implement this?) I had to sit down for a long while to try and remember! I think I might have done, but my recollection appears to be a little foggy. A lot has happened since then, to put it mildly. Either way, looking at the tech note, what I think might have happened somewhat prosaically, is that we updated the stomatal conducts ce équation to Medlyn, but then didn't realize that the later text on es refera to it. My suspicion is that this (ew 23) was relevant to the previous explanation of BB. Of course, Medlyn is using VPD and so that's not how it works any more. My suggestion is thus probably to look directly at the code (e.g. working from line 3311 here) As you can see there, the stomatal model is really a one-time switch, (with a couple of upstream calculations of vod, etc) so if we did something wrong it shouldn't be too hard to identify, hopefully. Let me/us know if that throws up any questions. Also you might be interested to know (if you don't already) that we are also testing the implementation of Medlyn vs BB in the FATES model extension of CLM with Alister Rogers at BNL and his post doc Cherry. They are at the stage of writing up their work but I'm quite sure they have no plans (/funding) to look at the climate impacts. You might be interested in exchanging note with them at some point. Best wishes, Dr Rosie A. Fisher. On Wed, 14 Jul 2021, 01:24 Nicholas Herold, [email protected] wrote: My name's Nicholas Herold and I'm using CLM5 to study the climatic impact of different stomatal conductance schemes, namely Ball-Berry and Medlyn. I've been conversing with Keith, Danica and Daniel from NCAR on an issue I'm having reconciling Gs values I get from CLM with ones I calculate with the Medlyn equation. And I'm trying to trace where the source of this difference comes from. I was told you may have implemented the numerical solution for Medlyn in CLM5. Is that correct? If so I was wondering if you could help join some dots for me. The CLM5 technote says that eq. 9.23 is substituted into eq. 9.1 to get the quadratic equation described in equations 9.24 - 9.26. But I don't see how this substitution works (i.e. how es can substitute into eq. 9.1). I'm interested in this because the solution to the quadratic that CLM is coming up with is different to what I get from eq. 9.1, which, if I understand correctly, shouldn't be happening (but I'd be happy to be told otherwise). This happens whether the plant hydraulic stress module is active or not. Anyway, any insight would be appreciated. Thanks for your time. Cheers Nick. |
|
Ok, a bit more progress. Below I'm implementing the analytical equation after calculating es and vpd. There are a couple of modifying factors applied to VPD and g0 in the CLM source which I've included below (a minimum constraint for VPD, and g0 being multiplied by bsun). Removing the latter brings gs values into even closer agreement. But still some large exceptions occur. Removing the minimum constraint for VPD seems to improve things a little further, i.e. agreement is to the point that the 99.4th percentile of errors is 8.38190317153931E-09. But with the VPD constraint removed I see the maximum discrepancy goes up to 16,797,683,955 and is due to the analytical solution going bonkers when the numerical solution is close to zero - so I imagine that the VPD constraint is there for good reason and should be kept (comment in the code implies Jinyun Tang might know something about this constraint)! So, it seems removing the bsun multiplication is a genuine improvement (and I see this modification isn't applied in the Medlyn numerical solution so that makes sense) but maybe the VPD constraint should be kept. Can someone comment on this? Does the VPD minimum constraint need tweaking to improve treatment of edge cases? I can't see any other reason for the remaining large discrepancies. When only removing the bsun multiplication, the max discrepancy is 5,039,221 and this seems to be when the numerical solution gets big (6,841,128) compared to the analytical solution (1,801,906). And the 95th percentile of the errors is 91.35. So it seems at first glance that the VPD constraint improves the most extreme discrepancies (~16.7 billion down to ~5 million) but increases discrepancies elsewhere. |
|
I've tried a slightly different approach following some of Gordon's Matlab code. The following code for sunlit doesn't throw any errors in a 2 year run at the US-UMB tower site. No errors on the shaded side either. |
|
That makes sense to me, Keith.
Also, sorry if I missed something, but will this PR add a bsun/bsha factor
for g0 with Medlyn?
That seems more like a science change than a numerical bug fix.
My impression was that many supported the change with the Medlyn
implementation to not attentuate g0.
…On Tue, Jul 20, 2021 at 9:55 AM Keith Oleson ***@***.***> wrote:
I've tried a slightly different approach following some of Gordon's Matlab
code. The following code for sunlit doesn't throw any errors in a 2 year
run at the US-UMB tower site. No errors on the shaded side either.
According to Gordon's code, the error check is only done (valid) when
(esat_tv(p) - ceair) .gt. 50. I've removed the multiplication of g0 by bsun
except when an < 0 (where gsminsun is multiplied by bsun in the
photosynthesis code). hs and rh_leaf_sun already exist in the code because
of the Ball-Berry error check, but this could be simplified using es = hs *
esat_tv as you've done.
hs = (gb_mol(p)*ceair + gs_mol_sun(p,iv)*esat_tv(p)) / ((gb_mol(p)+gs_mol_sun(p,iv))*esat_tv(p))
rh_leaf_sun(p) = hs
if ((esat_tv(p) - ceair) .gt. 50) then
vpd = max(esat_tv(p) - rh_leaf_sun(p) * esat_tv(p), 0.1_r8) * 0.001_r8
if (an_sun(p,iv) .ge. 0._r8) then
gs_mol_err = gsminsun + 1.6_r8*(1._r8+(gs_slope_sun/sqrt(vpd))) * (max(an_sun(p,iv),0._r8)/(cs_sun/forc_pbot(c)))
else
gs_mol_err = max(bsun(p)*gsminsun, 1._r8) + 1.6_r8*(1._r8+(gs_slope_sun/sqrt(vpd))) * (max(an_sun(p,iv),0._r8)/(cs_sun/forc_pbot(c)))
end if
if (abs(gs_mol_sun(p,iv)-gs_mol_err) > 1.e-01_r8) then
write (iulog,*) 'Medlyn error check - SUNLIT stomatal conductance error:'
write (iulog,*)'gs_sun: ',gs_mol_sun(p,iv),'gs_err: ',gs_mol_err
end if
end if
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
<#1398 (comment)>, or
unsubscribe
<https://github.com/notifications/unsubscribe-auth/ADNHR5WY2FUMHV4UVN7P62LTYWTA5ANCNFSM46RKMPGQ>
.
|
|
The existing CLM5 release code, that has Medlyn with PHS, has attenuated g0, but only for an<0: if (an_sun(p,iv) < 0._r8) gs_mol_sun(p,iv) = max( bsun(p)*gsminsun, 1._r8 ) and so the error check has to accommodate that. But my impression for the Medlyn-SMS PR was that we don't want to attenuate g0. |
|
That seems like a bug to me. I think the intent is to have the Medlyn model
g0 unattenuated.
…On Tue, Jul 20, 2021 at 10:35 AM Keith Oleson ***@***.***> wrote:
The existing CLM5 release code, that has Medlyn with PHS, has attenuated
g0, but only for an<0:
if (an_sun(p,iv) < 0._r8) gs_mol_sun(p,iv) = max( bsun(p)*gsminsun, 1._r8 )
if (an_sha(p,iv) < 0._r8) gs_mol_sha(p,iv) = max( bsha(p)*gsminsha, 1._r8 )
and so the error check has to accommodate that.
But my impression for the Medlyn-SMS PR was that we don't want to
attenuate g0.
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
<#1398 (comment)>, or
unsubscribe
<https://github.com/notifications/unsubscribe-auth/ADNHR5XVQKIY4G65T4S7WVDTYWXVTANCNFSM46RKMPGQ>
.
|
|
I'm fine with changing that, it will make the error check simpler also. |
|
We'll just have to wrap a medlyn/BB if around those statements, because I
think we do still want the bsun/bsha attenuation for ball-berry
…On Tue, Jul 20, 2021 at 10:46 AM Keith Oleson ***@***.***> wrote:
I'm fine with changing that, it will make the error check simpler also.
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
<#1398 (comment)>, or
unsubscribe
<https://github.com/notifications/unsubscribe-auth/ADNHR5XS3IZFLNJS2TJKCFDTYWY7DANCNFSM46RKMPGQ>
.
|
|
Good point. I'm also going to do a short global simulation to see if the error check is robust. Maybe Nick can review my changes to the error check and verify that it seems correct and works in his own simulation. |
|
I have to confess, I'm only following this discussion superficially, but I wondered if it would be helpful to discuss at an upcoming CLM meeting? |
|
Great. Ok a 10 day global case doesn't show any errors on my end using @olyson's formulation. And it's fine with PHS on or off. For the record I'm testing in my own modified CLM code but don't see why this wouldn't be robust. To implement this check in the Photosynthesis routine (i.e. when PHS is off) I removed the if condition for an >= 0 which is what I understood from the above exchange. Combined BB and MED error check for PHS = off: |
|
And to confirm what you were suggesting @djk2120. Should the following two lines inside PhotosynthesisHydraulicStress: be replaced with something like: Is there any equivalent change required when PHS = off? |
|
@heroldn would you mind reviewing the photosynthesis chapter in the technical note? I've updated the document per your suggestions: Let me know if I missed anything. |
|
I've read the parts that I made suggestions on (certainly not qualified to review 2.9 in general).
|
|
Thanks for the review. I'm not quite sure what you mean by a formatting error with "a=1", unless you mean how these equations are aligned. Yes, I think your code sample is right and this needs to be reproduced when PHS is off. |
Updates to Tech Note Photosynthesis chapter ### Description of changes Updates to photosynthesis chapter in technical note per discussion in #1398. ### Specific notes Correct some errors and clarify derivation of stomatal conductance quadratic equation for Medlyn method. Contributors other than yourself, if any: @djk2120 @heraldn CTSM Issues Fixed (include github issue #): Documentation aspects of #1398 Are answers expected to change (and if so in what way)? No Any User Interface Changes (namelist or namelist defaults changes)? No Testing performed, if any: I did a local build on the documentation.
|
Ah yes, that equation is just right-aligned, it looks inconsistent with other equations. Ok thanks re g0 attenuation. And when PHS=off, I now see that the way variables are named ensures g0 attenuation is already only applied to BB (in my code anyway - I was taking edits to allow Medlyn from Daniel's files). Thanks for your time on all this. |
Brief summary of bug
Implementing an error check on medlyn stomatal conductance, similar to the Ball-Berry error check, indicates an error
General bug information
CTSM version you are using: release-clm5.0.34
Does this bug cause significantly incorrect results in the model's science? Maybe
Configurations affected: All using PhotosynthesisHydraulicStress (use_hydrstress=.true.)
Details of bug
A user tried to implement an error check for Medlyn stomatal conductance (sunlit and shaded) in PhotosynthesisHydraulicStress, similar to the error check that is done for Ball-Berry in Photosynthesis.
Specifically, the error check calculates the difference between stomatal conductance produced by the solution in PhotosynthesisHydraulicStress and the stomatal conductance produced by equation 2.9.1 in the technical note, and triggers an error statement if the difference is greater than 0.1.
Here is the added code for sunlit stomatal conductance:
Here, gs_mol_err represents the stomatal conductance as calculated by equation 2.9.1 in the technical note.
I've replicated this behavior in a single point simulation at US-UMB.
Here is some sample output:
Medlyn error check - SUNLIT stomatal conductance error:
30503.1558514552 30136.0322013949
It looks to me as if the code for gs_mol_err is consistent with the tech note equation.
rh_can is not actually relative humidity, it is the vapor pressure deficit at the leaf surface which is calculated correctly in the code.
I don't know if the equation is actually wrong or if there is something wrong in how gs_mol_sun(sha) is calculated in the PhotosynthesisHydraulicStress. Also not sure how to figure this out.
Important details of your setup / configuration so we can reproduce the bug
This case demonstrates the problem:
/glade/work/oleson/release-clm5.0.34/cime/scripts/spinclm50conr34sp_US-UMB_chkmedlyn
The text was updated successfully, but these errors were encountered: