Questions about Cable Canopy

[mcuntz]

Dear all, @har917 @juergenknauer @Lo-hatch @mdekauwe @ManonSabot @gabsun

I create an issue here but if you think that the forum would be a better place, please tell me. I also added some people because I do not know who gets notifications for new issues.

We are working to include plant hydraulics in CABLE-POP (starting from the code of Martin de Kauwe and Manon Sabot). We came across a few points in cable_canopy that we do not understand.

Here the points that are in common between the CABLE main branch and the CABLE-POP_TRENDY branch:

1. The stomatal conductance models are supposed to depend on the humidity at the surface of the leaf. It is like this in Wang and Leuning (1998) and also in the CABLE documentation (Kowalczyk et al. 2006, Eq. 67).
   In the code, the variable is called dsx. I should be esat(T_leaf) - e_s with e_s the actual vapour pressure at the leaf surface. However, dsx is rather esat(T_leaf) - e_v with e_v the actual vapour pressure of the canopy air-space (met%qvair). The update of the variable in dryLeaf is: dsx(i) = met%dva(i) + air%dsatdk(i) * (tlfx(i)-met%tvair(i)) on line 548 in src/science/canopy/cbl_dryLeaf.F90 in MAIN and on line 2259 in core/biogeophys/cable_canopy.F90 in CABLE-POP_TRENDY. met%dva is (qstvair - met%qvair) * const.

   Q: Is this a deliberate change?
   Q: Do you know why?
   Q: When did this happen? I have code from 2015 that already has this change.
   Q: Do you think it is worth it reinstalling the dependence on surface vapour pressure.

2. There are x- and y-variables in the code, such as tlfx and tlfy for leaf temperature. The y-variables will be written to the output variables at the end of dryLeaf such as canopy%frday = 12.0 * SUM(rdy, 2) for leaf respiration. The y-variables are only updated in the iterations if the change in leaf temperature was smaller than any change of leaf temperature before in the iterations.

   There is also a line that does reduce the calculated change in temperature after 5 iterations (of 20) to avoid oscillations between states.

        tlfx(i) = ( 0.5 * ( MAX( 0, k-5 ) / ( k - 4.9999 ) ) ) *tlfxx(i) + &
             ( 1.0 - ( 0.5 * ( MAX( 0, k-5 ) / ( k - 4.9999 ) ) ) ) &
              * tlfx(i)
   

   So there should be oscillations in leaf temperature changes that are not suppressed by the regularisation of it to have a difference between x- and y-variables at the end. We ended up with oscillations in stomatal conductance, transpiration, etc. when including a new stomatal conductance model (Tuzet) because the regularisation after 5 iterations does only reduce the change in leaf temperature but not the other variables.

   Q: Why is only leaf temperature stabilised and not Cs and Ds as well?
   Q: If stabilisation worked, the difference between x- and y-variables would make no sense. What is your experience?
   Q: There are 4 outer iterations in define_canopy and 20 inner iterations in dryLeaf. Given that the y-variables are supposedly the "correct" output, should the inner iterations not start with the y-variables then (they start with the x-variables currently)?

3. The CO2 compensation point (Gamma not Gamma*) is set to 0: co2cp3 = 0.

   Q: Is there a reason? Is it supposed to be not necessary? Is it simply laziness to redo the quadratic equations?

There is one more remark for CABLE-POP:

1. Vanessa rewrote the Leuning model to reduce the correlation between the two parameters a1 and D0, which is large. This makes sense if you optimise parameters.
   However, the term a1 / (1 + D / D0)was rewritten to 1 / (1 / a1 + D / D0). The latter D0 should then be seen as a1*D0. But it is still around 1.5 kPa in the parameter file (rather than be around 13 kPa). So this should be changed somewhere.

   Q: We could adapt the parameter file. Or we could just undo the change and stay with the current parameters, which are easier to interpret.

Kind regards,
Matthias

[har917]

@mcuntz @juergenknauer @ManonSabot @mdekauwe @Lo-hatch @gabsun

Many of these questions pre-date my involvement with the code - what I can say is that this area of code was extensively rewritten by Yingping Wang somewhere in between CABLE1.4b and CABLE2.

On Q1. this part of the code was in the same/similar form back in cable1.4b. I suspect that this relates to the way goes about it's iterations - in that CABLE explicitly starts in isothermal conditions (where the leaf, canopy air and reference air are at the same temperature). That assumption makes the starting point correct - and it's the update that is questionable (but could be correct at a given level of linerisation of esat).

On Q2. - I personally hate the double iteration in CABLE - I think it's inefficient and I know it's ends up in an unconverged state more often than we realise. I have found that if you attempt to modify this part of the code (e.g. to correctly build the Penman-Monteith option for soil evaporation into within_canopy(), or build in alternate forms for the soil to canopy air resistances) the whole thing becomes highly twitchy and often settles on unrealistic values.

The stabilisation of temperature after 4 iterations looks to me like an (ugly) attempt to numerically solve an algorithm structure problem - with the hope that (since C and D are largely functions of leaf temperature) squashing the leaf temperature oscillations would also impact the other parts of the chain. It's been in the code since at least cable1.4b.

Remember that the two iteration loops are doing different things (the inner loop converges on leaf temperature given a u*, and the outer loop converges on static stability [MO length, u* etc.] given a stomatal conductance. The Raupach et al. (1997) [Soil-Canopy-Atmosphere-Model report - a reference that seems to have been lost from the collective CABLE memory] outer iteration assumes that you start each iteration in isothermal conditions, and evaluate a non-isothermal increment (which depends on the leaf and soil energy balance) - this logic means you shouldn't inherit the converged leaf temperature from the previous iteration (since it's no longer relevant given the update to u*).

Basically I think that this entire section of code should be rewritten - including an extension to allow proper sunlit-shaded leaf discrimination - but I have no idea who would do that/when that could happen.

On Q3 and Q4 - sorry I have no real idea. On the parameter value issue - could you state which parameter file you are referring to (there are lots going around)? This is possibly important for our recent ESM work @rml599gh