fwsoil with Haverd2013 appears broken

[har917]

Current runs from MAIN with the fwsoil_switch = 'Haverd2013' appear in some (maybe all) cases to result in a constant value of canopy%fwsoil . Some - but not all - fluxnet site runs also exhibit poor water conservation.

A quick comparison between the MAIN and CABLE-POP-TRENDY branches shows that there is an additional (re)initialisation of canopy%fwsoil at line 194 in MAIN - but it's not clear whether this is the cause or not.

This issue is a blocker for the use of this switch in AM3/CM3 and ESM1.6 - we should not use this switch/codebase combination in the coupled models as they stand.

[har917]

Interestingly - the (re)initialisation line above also appear in the CM2 codebase on MOSRS - so this bug (if it is one) was in the codebase and published simulations.

[har917]

Further information: the following is a list of initial & final values of fwsoil in the output from benchcab runs for #545 alongside the mean value of the absolute change in fwsoil over the time step i.e. |fwsoil[t] - fwsoil[t-1]|

For the R0_S2 (i.e. MAIN, fwsoil = 'standard') case:

site	start	final	test
AU-ASM_2011-2017_OzFlux_Met_R0_S2_out.nc	0.25894716	0.6965682	0.00048
AU-Ctr_2010-2017_OzFlux_Met_R0_S2_out.nc	0.41043082	1.0	0.00018
AU-Cum_2013-2018_OzFlux_Met_R0_S2_out.nc	0.78071135	0.80972254	0.00074
AU-GWW_2013-2017_OzFlux_Met_R0_S2_out.nc	0.2510517	0.13480236	0.00063
AU-How_2003-2017_OzFlux_Met_R0_S2_out.nc	1.0	1.0	0.00046
AU-Stp_2010-2017_OzFlux_Met_R0_S2_out.nc	0.42416772	0.8791535	0.00045
AU-Tum_2002-2017_OzFlux_Met_R0_S2_out.nc	0.4048594	0.9022395	0.00109
BR-Sa3_2001-2003_FLUXNET2015_Met_R0_S2_out.nc	1.0	0.55046844	0.00044
CA-Qfo_2004-2010_FLUXNET2015_Met_R0_S2_out.nc	1.0	1.0	0.00002
CH-Dav_1997-2014_FLUXNET2015_Met_R0_S2_out.nc	1.0	1.0	0.00006
CN-Cha_2003-2005_FLUXNET2015_Met_R0_S2_out.nc	1.0	0.57827055	0.00045
CN-Din_2003-2005_FLUXNET2015_Met_R0_S2_out.nc	1.0	0.7241024	0.00029
DE-Geb_2001-2014_FLUXNET2015_Met_R0_S2_out.nc	1.0	1.0	0.00082
DE-Gri_2004-2014_FLUXNET2015_Met_R0_S2_out.nc	1.0	1.0	0.00013
DE-Hai_2000-2012_FLUXNET2015_Met_R0_S2_out.nc	1.0	1.0	0.00030
DE-Tha_1998-2014_FLUXNET2015_Met_R0_S2_out.nc	1.0	1.0	0.00026
DK-Sor_1997-2014_FLUXNET2015_Met_R0_S2_out.nc	1.0	1.0	0.00032
FI-Hyy_1996-2014_FLUXNET2015_Met_R0_S2_out.nc	1.0	1.0	0.00008
FR-Gri_2005-2013_FLUXNET2015_Met_R0_S2_out.nc	1.0	1.0	0.00071
FR-Pue_2000-2014_FLUXNET2015_Met_R0_S2_out.nc	1.0	1.0	0.00075
GF-Guy_2004-2014_FLUXNET2015_Met_R0_S2_out.nc	1.0	1.0	0.00034
IT-Lav_2005-2014_FLUXNET2015_Met_R0_S2_out.nc	1.0	1.0	0.00010
IT-MBo_2003-2012_FLUXNET2015_Met_R0_S2_out.nc	1.0	1.0	0.00030
IT-Noe_2004-2014_FLUXNET2015_Met_R0_S2_out.nc	1.0	1.0	0.00074
NL-Loo_1997-2013_FLUXNET2015_Met_R0_S2_out.nc	1.0	1.0	0.00043
RU-Fyo_2003-2014_FLUXNET2015_Met_R0_S2_out.nc	1.0	1.0	0.00006
US-Blo_2000-2006_FLUXNET2015_Met_R0_S2_out.nc	1.0	1.0	0.00034
US-FPe_2000-2006_LaThuile_Met_R0_S2_out.nc	0.08081804	0.5607219	0.00062
US-GLE_2009-2014_FLUXNET2015_Met_R0_S2_out.nc	0.9892275	1.0	0.00019
US-Ha1_1992-2012_FLUXNET2015_Met_R0_S2_out.nc	0.9512565	1.0	0.00014
US-MMS_1999-2014_FLUXNET2015_Met_R0_S2_out.nc	1.0	1.0	0.00078
US-Me2_2002-2014_FLUXNET2015_Met_R0_S2_out.nc	1.0	1.0	0.00056
US-Myb_2011-2014_FLUXNET2015_Met_R0_S2_out.nc	1.0	1.0	0.00059
US-NR1_1999-2014_FLUXNET2015_Met_R0_S2_out.nc	0.7753052	1.0	0.00073
US-PFa_1995-2014_FLUXNET2015_Met_R0_S2_out.nc	1.0	1.0	0.00066
US-SRG_2009-2014_FLUXNET2015_Met_R0_S2_out.nc	0.76381534	0.6342633	0.00076
US-SRM_2004-2014_FLUXNET2015_Met_R0_S2_out.nc	0.76381534	0.66914326	0.00070
US-Ton_2001-2014_FLUXNET2015_Met_R0_S2_out.nc	1.0	1.0	0.00047
US-UMB_2000-2014_FLUXNET2015_Met_R0_S2_out.nc	1.0	1.0	0.00140
US-Var_2001-2014_FLUXNET2015_Met_R0_S2_out.nc	1.0	1.0	0.00045
US-Whs_2008-2014_FLUXNET2015_Met_R0_S2_out.nc	0.76381534	0.37733287	0.00055
US-Wkg_2005-2014_FLUXNET2015_Met_R0_S2_out.nc	0.85242164	0.3563604	0.00059

None of the sites are invariant in time - though some start and end at the same (unrestricted saturation) value

For the R0_S0 case ( MAIN, 'Haverd2013')
case:

site	start	final	test
AU-ASM_2011-2017_OzFlux_Met_R0_S0_out.nc	0.06776762	0.06776762	0.00000
AU-Ctr_2010-2017_OzFlux_Met_R0_S0_out.nc	0.3176543	0.3176543	0.00000
AU-Cum_2013-2018_OzFlux_Met_R0_S0_out.nc	0.8118866	0.8118866	0.00000
AU-GWW_2013-2017_OzFlux_Met_R0_S0_out.nc	0.09268184	0.09268184	0.00000
AU-How_2003-2017_OzFlux_Met_R0_S0_out.nc	0.79319286	0.79319286	0.00000
AU-Stp_2010-2017_OzFlux_Met_R0_S0_out.nc	0.25389865	0.25389865	0.00000
AU-Tum_2002-2017_OzFlux_Met_R0_S0_out.nc	0.6679908	0.6679908	0.00000
BR-Sa3_2001-2003_FLUXNET2015_Met_R0_S0_out.nc	0.9151175	0.9151175	0.00000
CA-Qfo_2004-2010_FLUXNET2015_Met_R0_S0_out.nc	0.74438316	0.74438316	0.00000
CH-Dav_1997-2014_FLUXNET2015_Met_R0_S0_out.nc	1.0	0.9080068	0.00000
CN-Cha_2003-2005_FLUXNET2015_Met_R0_S0_out.nc	0.9318116	0.9318116	0.00000
CN-Din_2003-2005_FLUXNET2015_Met_R0_S0_out.nc	0.92355466	0.92355466	0.00000
DE-Geb_2001-2014_FLUXNET2015_Met_R0_S0_out.nc	0.9103051	0.9103051	0.00000
DE-Gri_2004-2014_FLUXNET2015_Met_R0_S0_out.nc	0.8603512	0.8603512	0.00000
DE-Hai_2000-2012_FLUXNET2015_Met_R0_S0_out.nc	0.9103051	0.9103051	0.00000
DE-Tha_1998-2014_FLUXNET2015_Met_R0_S0_out.nc	0.8603512	0.8603512	0.00000
DK-Sor_1997-2014_FLUXNET2015_Met_R0_S0_out.nc	0.9368293	0.9368293	0.00000
FI-Hyy_1996-2014_FLUXNET2015_Met_R0_S0_out.nc	0.9242546	0.9242546	0.00000
FR-Gri_2005-2013_FLUXNET2015_Met_R0_S0_out.nc	0.9062017	0.9062017	0.00000
FR-Pue_2000-2014_FLUXNET2015_Met_R0_S0_out.nc	0.89488924	0.89488924	0.00000
GF-Guy_2004-2014_FLUXNET2015_Met_R0_S0_out.nc	0.931098	0.931098	0.00000
IT-Lav_2005-2014_FLUXNET2015_Met_R0_S0_out.nc	0.86442274	0.86442274	0.00000
IT-MBo_2003-2012_FLUXNET2015_Met_R0_S0_out.nc	0.9117805	0.9117805	0.00000
IT-Noe_2004-2014_FLUXNET2015_Met_R0_S0_out.nc	0.86444867	0.86444867	0.00000
NL-Loo_1997-2013_FLUXNET2015_Met_R0_S0_out.nc	0.9382439	0.9382439	0.00000
RU-Fyo_2003-2014_FLUXNET2015_Met_R0_S0_out.nc	0.9657198	0.9657198	0.00000
US-Blo_2000-2006_FLUXNET2015_Met_R0_S0_out.nc	0.7982836	0.7982836	0.00000
US-FPe_2000-2006_LaThuile_Met_R0_S0_out.nc	0.0	0.0	0.00005
US-GLE_2009-2014_FLUXNET2015_Met_R0_S0_out.nc	0.0	1.0	0.00024
US-Ha1_1992-2012_FLUXNET2015_Met_R0_S0_out.nc	0.82209635	0.82209635	0.00000
US-MMS_1999-2014_FLUXNET2015_Met_R0_S0_out.nc	0.91186976	0.91186976	0.00000
US-Me2_2002-2014_FLUXNET2015_Met_R0_S0_out.nc	0.6708635	0.6708635	0.00000
US-Myb_2011-2014_FLUXNET2015_Met_R0_S0_out.nc	0.8319932	0.8319932	0.00000
US-NR1_1999-2014_FLUXNET2015_Met_R0_S0_out.nc	0.1764814	0.1764814	0.00000
US-PFa_1995-2014_FLUXNET2015_Met_R0_S0_out.nc	0.8470388	0.8470388	0.00000
US-SRG_2009-2014_FLUXNET2015_Met_R0_S0_out.nc	0.6186793	0.6186793	0.00000
US-SRM_2004-2014_FLUXNET2015_Met_R0_S0_out.nc	0.6186793	0.6186793	0.00000
US-Ton_2001-2014_FLUXNET2015_Met_R0_S0_out.nc	0.7982836	0.7982836	0.00000
US-UMB_2000-2014_FLUXNET2015_Met_R0_S0_out.nc	0.9387465	0.9387465	0.00000
US-Var_2001-2014_FLUXNET2015_Met_R0_S0_out.nc	0.7982836	0.7982836	0.00000
US-Whs_2008-2014_FLUXNET2015_Met_R0_S0_out.nc	0.6186793	0.6186793	0.00000
US-Wkg_2005-2014_FLUXNET2015_Met_R0_S0_out.nc	0.6573827	0.6573827	0.00000

Nearly all of the sites have an invariant fwsoil - this is certainly not the intent of the parametrisation.

We need to check the corresponding values from a run triggered off the CABLE-POP_TRENDY branch .

[rml599gh]

Is this the same issue as #484, which we found when I was trying to get the Haverd2013 switch working in ESM1.5 (#483)?
If so, I thought we'd fixed this - but perhaps it was only in my local code - or we assumed it would be picked up with AM3 code - or has it been done in the CABLE code in the UM7 repository but not in the CABLE repo??

[har917]

Is this the same issue as #484 ?

Possibly - but probably not. Tests associated with #494 are showing the same issue - and that's PR has the adjustments to %wbliq in cbl_soil_snow_main.

The results above - from #545 - are from a branch which doesn't have the ice fixes (or wbliq term) in.

Further investigations seem to be indicating that the water balance errors are focussed in short periods of time when the soil T > freezing, and soil moisture <= field capacity. fwsoil = 'standard' is reducing fwsoil and hence transpiration, fwsoil = 'Haverd2013' seems to be leaving the value of fwsoil alone (allowing for unrestricted transpiration) but restricting amount that can be taken from the soil.

[har917]

I now understand why we're getting somewhat inconsistent impacts from the 'Haverd2013' parameterisation.

Below are 4 plots showing time series (blue = 'standard' R0_S2, orange = 'Haverd2013', R0_S0) for fwsoil, latent heat Qle (basically transpiration), wbal (the integrated water imbalance) and soil moisture in each of the 6 layers - from the benchcab runs for UMB from #545.

As noted 'Haverd2013' is not varying the fwsoil, this allows unrestricted transpiration, which is driving the soil moisture below the standard case (and even below the wilting point swilt). However it is only if and when this drives the water content below one or other of the hard-wired limits (0 soil moisture in soil layer 4) that this manifests as a water imbalance. If you don't hit such a case then 'Haverd2013' is conserving well, albeit the dynamics are incorrect - hence the inconsistency in results across cases.

All this is pointing towards the need to understand what's happening with fwsoil with Haverd et al. as the route to explain this issue.

[ccarouge]

Going through the code, the likely culprit is this:

CABLE/src/science/canopy/cbl_dryLeaf.F90

Lines 744 to 758 in d744c2f

	IF (ANY(((rex*dt) > MAX((theta(:)-thetaw(:)),zero)*dx(:)) .AND. (Etrans > zero))) THEN
	fws = zero
	! distribute extraction according to available water
	!rex(:) = (theta(:)-0.01_r_2)*dx(:)
	rex(:) = MAX((theta(:)-thetaw(:))*dx(:),zero)
	trex = SUM(rex(:))
	IF (trex > zero) THEN
	rex(:) = rex(:)/trex
	ELSE
	rex(:) = zero
	ENDIF
	rex(:) = Etrans*rex(:)
	ELSE
	fws = MAXVAL(alpha_root(2:)*delta_root(2:))
	ENDIF

The ELSE condition will give a fwsoil value of about 0.94. So it looks like the condition is always false. Investigating more to see why that would be.

[ccarouge]

I obviously don't know anything about that science but this line is weird:

CABLE/src/science/canopy/cbl_dryLeaf.F90

Line 740 in d744c2f

rex(:) = Etrans*rex(:)

In the case, I looked at Etrans is very small (1e-11), so it makes rex very small.

Before that line we have:

CABLE/src/science/canopy/cbl_dryLeaf.F90

Lines 729 to 731 in d744c2f

	WHERE (((rex*dt) > (theta(:)-thetaw(:))*dx(:)) .AND. ((rex*dt) > zero))
	alpha_root = alpha_root*(theta(:)-thetaw(:))*dx(:)/(1.1_r_2*rex*dt)
	endwhere

This is pretty much the same condition as the one used to update fws but at line 729, rex is much larger than at the line 744 since it wasn't multiply by Etrans. So either we need to add Etrans in the second condition or rex shouldn't be multiply by it before calculating fws (note rex=Etran*rex in the code that updates fws as well, so it is done twice potentially).

[har917]

@ccarouge the CABLE-POP_TRENDY has additional comments in this section - which in effect are reminders from @mcuntz to ask Vanessa about the logic here.

The basic idea is that if there is any soil layer which cannot support the requested transpiration then you have to reshuffle where the water comes from, and if necessary rescale the amount requested. I agree with @mcuntz though that in this situation - only if rex(:) ends up at zero should fws be zero, and that's not what is happening.

Equivalently though - if there's been some reshuffling of rex(:) because of water limitations I'm not sure that fws should be =1.

So yes - this area is suspect.

@juergenknauer does this raise any memories?

Another thing that I don't like about all of this is that the way this is coded effectively makes fwsoil a state variable - fwsoil is initialised from canopy%fwsoil at the start of dryLeaf . This means it should go in the restart. It is in the CABLE-POP_TRENDY restarts but it's not in MAIN or the coupled models.

[har917]

@ccarouge Part of the problem here is that rex changes it's meaning part way through the routine. Also be careful in the test cases that you look at - many start at local midnight so, naturally, Etrans and rex will be very small.

rex starts off being the relative extraction - so the fraction of the total extraction that's to come from each layer.
Line 725 then converts the relative extraction to an actual extraction rate.

The conditions at 729 then check whether actual extraction rate over the time step is sufficient to force the soil moisture content below wilting point (thetaw in this routine) - and updates relative extraction if needed.
Line 740 then converts the (potentially) revised relative extraction to an actual extraction.

Line 742 then goes over the whole set up again trying to catch edge cases - using different logic to reshuffle water extraction (including - magically - allowing extraction from soil layers without any roots.)

I suspect that most cases will be caught by the section from 728-740 - but, yes, it's highly likely that 744-755 would cause water imbalance if the model were to get into that case, especially if rex ends up as zero but Etrans not.

[ccarouge]

After re-reading the code, it seems to me that rex=Etrans*rex should be done at the end of the subroutine only.

This is because line 748, tells us the relative extraction is comparable to theta-thetaw. But at other places in the code it compres the actual extraction to theta-thetaw. Both can't be correct.

[har917]

After re-reading the code, it seems to me that rex=Etrans*rex should be done at the end of the subroutine only.

Unfortunately it's not that simple:

• the WHERE at line 729 and 744 are testing actual extraction over the time step (rex * dt with rex being actual extraction) against available water ((theta-thetaw)*dx, dx is soil layer thickness), i.e. in kg water per m2.
• the lines at 748-754 are resetting the rex as relative extraction - and does based on available water amounts - the units don't match; then line 755 converts the new relative extraction to actual extraction.

In short - I think the logic around what rex is compared to at each point is sound, just confusing.

I (and Matthias) am more worried about the combination (in that 742-758 section) of fws=0, Etrans>0 and rex>0. Matthias noted that (because of the way the canopy%fwsoil is inherited from call to call within dryLeaf, call to call from define_canopy, and [CABLE-POP_TRENDY] from time step to time step) once fws=0 then it is very difficult to get out of that condition* **. Additionally I don't like it because the code can set fws=0, rex=0.0 but leave Etrans>0, or it can set fws=0 ,SUM(rex)=Etrans>0 but rex in a particular layer still violate theta>thetaw. Conservation implies that Etrans must change in those conditions and one of the other energy fluxes should take up the difference.

*this is actually why I think this conversation is a bit of a mis-direction to be honest. I want to check the output from #494 to see what is happening in those runs in more detail.
**Earlier I noted that MAIN has the extra canopy%fwsoil = 1.0 when it enters define_canopy each time step. This actually partially solves this possible never-ending fws=0 trap.

[mcuntz]

The idea behind Haverd2013 is that plants take up water as long there is water within the root zone. This is also called compensating roots in the literature. Plants take up water from deeper layers if the upper layers fall dry even when there is a smaller amounts of roots. This is why fwsoil should be 1 most of the time.
Then you also have to say where the water is taken from the soil. The uptake follows the amount of water in a soil layer and the root length distribution. The function of the dependence on soil water is the _root shut-down function" of Lai and Katul (2000) (their alpha_2).

That's the idea. But we in Nancy also have problems with Haverd2013.

rex is not really changing its meaning. It is root extraction. Only, we do alpha*root_length_density first, normalise it to 1, and then multiply it by Etrans to have rex. We could have defined another variable with a random name, which would not have helped a lot, so we used rex for alpha*root_length_density already and then update it.

Then we check the water balance. Say a = alpha, g = root_length_density, E = Etrans. delta in the code is only there to limit root water extraction to the root zone.
Then we have rex = a * g * E. We check if we take out too much water from a layer: rex*dt > (theta-theta_wilt)*dx. If so then we do rex = a * (theta-theta_wilt)*dx / (1.1 * rex * dt) * g * E, which is rex = a * (theta-theta_wilt)*dx / (1.1 * a * g * E * dt) * g * E = (theta-theta_wilt)*dx / (1.1 * dt) So in theory, we should not take too much water from a layer.
But I think that the normalisation step rex = rex / sum(rex), which is a*g in this moment, destroys this safeguard and we can still take out too much water from a layer.

So then there is the second check if any rex*dt > (theta-theta_wilt)*dx. Then indeed the whole soil column is taken and fwsoil = 0. This is because in the logic of the code it means that there is not enough water within the root zone for transpiration. Take it from the rest of the layer in this time step but fwsoil = 0 next time.

After looking at it now, I think

1. We should not normalise a * g a second time.
2. A first check would be if there is enough water in the total root zone to fulfil transpiration.
   a) If so, we should perhaps code a loop from the top to the bottom, checking each layer, and adding the excess root extraction to the next layer. If there is still excess extraction in the last rootzone soil layer, we would loop up again.
   b) If there is not enough water in the root zone, then we could do the first loop from top to bottom and take the excess water demand from the layer underneath to have a closed water balance.

All this was for water extraction. But what is fwsoil? At the moment it is maxval(a(2:)) or 0.
I would have taken a mean soil water sum(theta/dx)/sum(dx) over the root zone and calculate fwsoil = alpha with it. This follows abaout what is observed in forests (e.g. Granier et al. (2007)). The gamma parameter in Lai & Katul can then be different for eucalypt (gamma <<) and conifers (gamma >>). It is a parameter in def_veg_params.txt in CABLE-POP, ranging from 0.01 to 0.035.

[har917]

@mcuntz - I agree that the second normalization could lead to too much water being extracted from a layer, but not doing so prevents the plants from preferentially extracting water from a layer with less roots but more water and also breaks SUM(rex*dx)=Etrans .

The 'standard' code takes the top-bottom loop approach (on loop, applied in all circumstances) - in cbl_remove_trans - but that scheme has other conservation issues.

---

Just one question on this - I didn't quite follow your last point about how to evaluate fwsoil. Did you mean xx=sum(theta/dx)/sum(dx) or xx=sum(theta * dx)/sum(dx)? I thin kit should be the * option.

and then (to check) did you mean ?

lthetar_mean = LOG( MAX(xx,e3)/(SUM(thetaS*dx)/SUM(dx)) )  !with thetaS*dx or thetaS/dx in as per definition of xx
fwsoil = EXP( gamma/MAX(xx,e3) * lthetar_mean )

These suggested changes are likely for 'future consideration' - the immediate work is to try and figure out why #494 isn't conserving as well as MAIN when the 'Haverd2013' switch is on (since that indicates there may be a problem with the liquid/ice density changes code).

---

Is there a physiological reason/interpretation for the factors 1.1? This appears in the 'standard' code as well. Given that there is also the e3 safeguard this seems to be somewhat overly safe.

[mcuntz]

Sorry, should have been a multiplication for dx, of course.

The 1.1 comes from the fact that too much water was extracted from the layers if 1.0, most probably due to the normalisation. So the 1.1 makes it more likely that we stay above the wilting point after root extraction.

For fwsoil, I meant:

theta_mean = sum(theta(1:rlayer) * dx(1:rlayer)) / sum(dx((1:rlayer))
thetaw_mean = sum(thetaw(1:rlayer) * dx(1:rlayer)) / sum(dx((1:rlayer))
if ((theta_mean - thetaw_mean) > e3) then
    thetas_mean = sum(thetas(1:rlayer) * dx(1:rlayer)) / sum(dx((1:rlayer))
    ! this is: ((theta-thetaw)/(thetas)**(gamma/(theta-thetaw))
    lthetar_mean = LOG(MAX(theta_mean-thetaw_mean, e3) / thetas_mean)
    fwsoil = EXP(gamma / MAX(theta_mean - thetaw_mean, e3) * lthetar_mean)
else
    fwsoil = 0.0_r_2
endif

[mcuntz]

By the way, in CABLE main, the constants e3, one and zero are REAL(r_2) but the values are given as single precision, e.g. one = 1.0. This makes small errors that do are do not add up. Kills the isotopes, though ;-)
It should be 1.0_r_2, etc.