Generating diagnostics for photosynthesis

[serbinsh]

Hi @mpaiao and/or @rgknox

@Viskari and I are attempting to generate some standard photosynthesis "response curves" using some simple diagnostic ED2 runs, in some cases with made up meteorology, but in other cases for a real site but with a simple single cohort single patch runs. We am doing this to better understand the photosynthetic response of ED2 (not GPP per se but individual leaf/cohort for a "single leaf" and with increasing LAI, i.e. scaled within the canopy), in terms of leaf and canopy phys, to different meteorological drivers as well as VPD and soil moisture. This is for a manuscript as well as for our general interest in what the formulation within ED2 actually produce, once all the other ins/outs are taken into account. Of course we could simply implement the code outside of ED2 but I would rather obtain the appropriate output from ED2 so that we know it used the ED2 solver and so that we can examine and diagnose outputs for typical met conditions.

As such, I am attempting to figure out what outputs we would need for this. The met drivers are easy enough, however VPD and soil moisture will be a little tricky. There a several variables that could be of interest (e.g. A_open/close, RUBP, LIGHT, CO2) but most do not produce sensible graphs, e.g. photo versus light/PAR (i.e. light response curve). What should I be looking for that would represent the photosynthetic rate/assimilation rate for a typical "leaf" for a cohort? In other words, is there an output in the fast files that I can use for this?

In summary we want to generate

A vs PAR for a "leaf" and various LAIs A vs temperature A vs Ca cureves A vs VPD A vs SWC

If we know what output(s) to look at that we could use to represent photosynthetic rate than we can design some runs to get these results. Also, if this is possible we could consider creating these as diagnostics that would let us know if the response surfaces make sense or if there is an introduced bug.

[mpaiao]

Hi @serbinsh and @Viskari,

I'm not sure why the ARUBP/LIGHT/CO2 aren't working, last time I checked they looked fine. Which HDF5 variable and which output file type (-S-, -I-, -Q-) are you using? I think for you analysis you probably want to use the -I- files, and FMEAN[variable]_CO, though you could also use the history files as long as you only use instantaneous variables.

Unless there is some bug, these are the fmean variables I would recommend for this analysis: PAR – FMEAN_PAR_L_CO / LAI. This is the PAR absorbed by leaves, but you must divide it by LAI to be m2(leaf) Temperature – FMEAN_LEAF_TEMP_CO Ca – FMEAN_CAN_CO2_PA (same for all cohorts within each patch) VPD – FMEAN_LEAF_VPDEF_CO – VPD at the leaf surface SWC – Probably the best variable is FMEAN_WATER_SUPPLY_CO, though this is not exactly available water. FMEAN_AVAILABLE_WATER_PA may be an alternative, but this is shared amongst all cohorts in the patch.

I just saw that we don't have FMEAN_A_OPEN, FMEAN_A_CLOSED and FMEAN_A_NET in this version. It may be worth calculating these variables and including them in the output, if you track every occurrence of A_RUBP, you should be able to find all places that you must add the variables. As a quick solution, you could get FMEAN_A_NET_CO = (FMEAN_GPP_CO-FMEAN_LEAF_RESP_CO)*NPLANT/LAI (plus unit conversions, at some point we should probably standardise units in ED).

A few other variables that may be useful too: FMEAN_PSI_OPEN_CO and FMEAN_PSI_CLOSED_CO are transpiration rates. If you use the cohort level, I think they are already in m2(leaf), so you don't need to divide by LAI. FMEAN_LEAF_GSW_CO: stomatal conductance, already down-regulated by fs_open.

[serbinsh]

Hi @mpaiao

Thanks! This is very, very helpful. Also, I agree that it may be useful to output FMEAN_A_OPEN, FMEAN_A_CLOSED and FMEAN_A_NET.

We have been looking into the -S-, -I- and -T- files, but it is probably best to use the -I- outputs.

@Viskari I think this gives us what we need, provided we generate some test case met driver data for the runs.

[serbinsh]

Hi @mpaiao

I have been racking my head trying to get the values to come out to sensible numbers for FMEAN_A_NET. First, I can't run ED2 with _CO level outputs as there is currently a bug causes the runs to crash.

In the interim I have been running ED2 with a small (<2) number of cohorts for 2 days in mid July to try and generate sensible response curves. To get around the _CO issue I am using _PY outputs in the -I- files, running with 1800 timestep (30 minutes).

@Viskari and I backed out the conversion and are using this to try and convert GPP and leaf resp into Anet in umols/m2/s

FMEAN_A_NET_PY = (GPP-Leaf_Resp_PY / 0.38)*(NPLANT/LAI)

But the values do not come out to sensible numbers. I think I am missing something or am approaching this wrong. Basically, I want to calculate photosynthesis for a single leaf (or almost 0 LAI) and for increasing LAI in umols/m2/s. Any thoughts on the errors I might be making or is there another conversion I need to make for polygon level variables?

@mdietze Any thoughts? Any chance we could chat tomorrow in case you have some ideas? I get the right shapes but just not the right values.

[serbinsh]

@mpaiao

We are assuming GPP/Leaf Resp are in kgC / m2 / yr, even at the fast time steps. This is based on the units expressed in the code.

[mpaiao]

Hi @serbinsh Yes, it is in kgC/m2/yr, and at some point we should probably make all units consistent...

Did you divide both GPP and Leaf respiration by 0.38? From your previous post it looks like you were dividing leaf respiration only. Otherwise, the conversion to m2(leaf) looks fine.

Also, what is the error message that you get when you try to run -I- files with cohort-level output? I just tried to run a quick test with my version (very close to the master, but not with the most recent updates) and it seems to be working fine even with the strict debugging flags.

By the way, I just updated my branch, now it has FMEAN_A_OPEN, FMEAN_A_CLOSED and FMEAN_A_NET in the output, feel free to try it. I used the same unit convention for the other "A" variables: cohort-level variables in µmol/m2(leaf)/s, and polygon-level variables in µmol/m2(ground)/s. I haven't tested it much yet so I'm not submitting a pull request for the time being.

[serbinsh]

Hi @mpaiao

Yes, I divided both GPP and LR by 0.38. I just pulled down your version of ED2 and am going to give that a shot now. I will let you know what I find!

[serbinsh]

@mpaiao

Here is what I was able to come up with using a Vcmax of 60 at 25C (scaled to a Vm0 at 15) from leaf (LAI 0), LAI1,3,and 7. Does this seem as expected? Just trying to get a handle on what to expect from ED2.

@Viskari @mdietze

[serbinsh]

@mpaiao

I attempted to have a single dominate cohort in these runs. Also, I used the PY versus CO for the LAI1-7 plots. I think that should be correct, if I want to show what Anet looks like for an ED "canopy" of increasing LAI. Interestingly, the Anet_CO values (I had 2 cohorts, 1 dominate and 1 sub-dominate) decreased with increasing LAI. Also, I noted that there were two specific incident PAR vales (~652 and ~847 umols/m2/s) where the Anet values would jump well above the rest of the response curve. This was true from leaf through increasing LAI. This was also consistent around these two values of PAR. Maybe a region of inflection in the solver?

[mpaiao]

@serbinsh Leaf is the cohort variable and LAI are the polygon, right? If so, I think the leaf makes sense (the response per leaf should have the kink that represents the transition from light-limited to RuBP-saturated productivity). I ran a quick test in R with your Vm0 and the switch happens more or less at the same point (other things like temperature, D0, respiration, etc will affect the point).

Not sure if it is the case, but higher LAI may also result in more light being "wasted". Top cohorts may absorb more light than what they need to overcome light limitation, and the understory cohorts may not receive enough light. Also, if the temperature is high (thinking of 25°C or higher) then more light at the top may warm leaves too much, increase respiration and VPD, and reduce A_net. ED simulations in the tropics usually end up with more biomass and LAI in higher elevation and this is one of the reasons.

I've never noticed spikes in the assimilation rates, and they are definitely not expected. Do you know if the jumps are associated with the light or RuBP cases? The photosynthesis solver that Naomi and I came up should have only one root in the reasonable range. I don't know what could cause the high values, but two guesses:

1. The reasonable range is defined by two points when the function goes to -Infinity and +Infinity; we determine these bounds analytically before searching the root and make sure guesses never go outside the bounds, but maybe you found some exception.
2. This is a terrible case for root finding: the function is very flat near the answer so the tolerance has to be very strict, otherwise it can accept values that are quite far from the answer. I think this is what you meant by inflection, right? Not sure why it could fail at these two specific values, but maybe the tolerance is not sufficient.

[serbinsh]

Hi @mpaiao

I need to explore this further, eg. light or RuBP case, but here is an example of what I am talking about. Here is a run with a single cohort and an average LAI of 1. I am plotting Anet (umols/m2/s) versus PAR (umol/m2/s) with an average leaf temp of 25C and Vcmax of 60 @ 25C. Note the two points around 652 and 847 umols/m2/s

@Viskari @mdietze

[serbinsh]

@mpaiao Here is an example of a leaf level response:

@Viskari @mdietze

[serbinsh]

OK, after changing:

10

so that I would get a single cohort for all runs, here are the updated light response curves for 25C

@mpaiao @mdietze @Viskari

[mpaiao]

@serbinsh after looking at the scatter plots, I'm wondering whether the spikes are really caused by photosynthesis or something like instabilities on leaf VPD or leaf temperature solution. Perhaps just to isolate the problem, could you try a similar singular cohort test, but now going to the second lphysiol_full call in photosyn_driv.f90, and replacing all environmental controls (except leaf_par) by constants? Thanks.

            call lphysiol_full( & !
               100000.          & ! Canopy air pressure                         [       Pa]
             , 1.158067         & ! Ideal gas density (1000hPa,25C)             [    kg/m3]
             , 0.015            & ! Canopy air sp. humidity                     [    kg/kg]
             , 400.             & ! Canopy air CO2 mixing ratio                 [ µmol/mol]
             , ipft             & ! Plant functional type                       [      ---]
             , leaf_par         & ! Absorbed photos. active rad.                [ W/m²leaf]
             , 298.15           & ! Leaf temperature                            [        K]
             , 0.01995147       & ! Leaf intercellular spec. hum. (1000hPa,25C) [    kg/kg]
             , green_leaf_factor(ipft) & ! Greenness rel. to on-allometry       [      ---]
             , leaf_aging_factor(ipft) & ! Ageing parameter to scale VM         [      ---]
             , cpatch%llspan(ico)      & ! Leaf life span                       [       yr]
             , cpatch%vm_bar(ico)      & ! Average Vm function                  [µmol/m2/s]
             , 0.012                   & ! Aerodyn. condct. of water vapour     [  kg/m2/s]
             , ...

@mdietze @Viskari

[mpaiao]

@serbinsh I ran the suggestion above of fixed environmental conditions, using a single cohort (PFT=3, LAI =1, but I don't think it matters for the problem), and the response looks like we would expect.

BTW, if you set NL%IDETAILED=2, ED will write a simple text file for each cohort with all the environmental information used to calculate the assimilation rates for every photosynthesis call. This may be useful to identify what is causing the spikes. The subroutine that writes the ascii files is print_photo_details (ED/src/dynamics/photosyn_driver.f90).

@mdietze @Viskari

[serbinsh]

Thanks @mpaiao I will take a look at these suggestions, as well as setting NL%IDETAILED=2. The figure you generated definitely suggests it isn't the photo. I noticed the inflection point in your graph is rather low (~150 umols) but that is probably a result of the test conditions.

Also, I was plotting Anet (using your version of the model the outputs this variable) so I will give your same test a go looking at that as well. I will also look at the response of leaf resp.

Thanks and I will keep you informed of anything I may find out!

[serbinsh]

@mpaiao

Are you planning to check in your changes to the outputs to inlcude the A\ options, such as A_OPEN, RUBP, and ANET? I have found those very useful and it would be good to include moving forward.

[mpaiao]

@serbinsh Sure, I will submit a pull request, probably tomorrow. I just want to make sure that the other stuff I've been working (tropical allometry) doesn't go in the same pull request, since I'm still testing it.