494 redux local ice fixes

[JhanSrbinovsky]

CABLE

Thank you for submitting a pull request to the CABLE Project.

Description

Please include a brief summary of the change and list the issues that are fixed.
Please also include relevant motivation and context.

You can link issues by using a supported keyword in the pull request's description or in a commit message:

Fixes #(issue)

Type of change

Please delete options that are not relevant.

• Bug fix
• New or updated documentation

Checklist

• The new content is accessible and located in the appropriate section
• I have checked that links are valid and point to the intended content
• I have checked my code/text and corrected any misspellings

Testing

• Are the changes bitwise-compatible with the main branch? If working on an optional feature, are the results bitwise-compatible when this feature is off? If yes, copy benchcab output showing successful completion of the bitwise compatibility tests or equivalent results below this line.

• Are the changes non bitwise-compatible with the main branch because of a bug fix or a feature being newly implemented or improved? If yes, add the link to the modelevaluation.org analysis versus the main branch or equivalent results below this line.

Please add a reviewer when ready for review.

---

📚 Documentation preview 📚: https://cable--527.org.readthedocs.build/en/527/

[JhanSrbinovsky]

@ccarouge I have several issues with testing required (apart from having limited time):

• The issue here didnt even show up until global run (is this even possible in benchcab/me.org). Now that that is otherwise "resolved" probably single -site tests are ok
• Anything short of uploading from gadi is untenable. My failed attempt and request for help (https://forum.access-hive.org.au/t/benchcab-me-org/4056) is still hanging

[JhanSrbinovsky]

with end of #494 comments/modifiations

[har917]

A couple of new things noted - now as new issues.

Please check i) how frozen_limit and Cdensity_ice are getting into spec_soilsnow_init and ii) why we need a REAL() in one of the SUBROUTINE argument calls - if they're okay then this PR is approved.

[JhanSrbinovsky]

So @ccarouge - what are we going to do about testing here?

[har917]

okay - we have a problem - lines 502-514 in cable_checks

I think the easiest way around this is to replace the local re-evaluation of the total water (liquid++frozen). and hence change in water content, with that evaluated in ssnow%wbtot. This is available from soil_snow and soil_snow_gw. SLI already has it's own wbal version in cable_checks which correctly incorporates frozen/liquid densities.

In effect the mass_balance routine is changed to (variable definition):

! Local variables
REAL(r_2), DIMENSION(mp,2), SAVE :: bwb         ! LOCAL total (liquid+ice) soil moisture
REAL(r_2), DIMENSION(mp)              :: delwb      ! change in soil moisture

and (first time step and main calculation)

IF(ktau==1) THEN
       ! initial value of soil moisture content - liquid/ice densities included
       bwb(:,1)=ssnow%wbtot
       bwb(:,2)=ssnow%wbtot
       delwb(:) = 0.0

ELSE
       ! Calculate change in soil moisture b/w timesteps:
       IF(MOD(REAL(ktau),2.0)==1.0) THEN            ! if odd timestep
             bwb(:,1)=ssnow%wbtot
             delwb(:) = bwb(:,1) - bwb(:,2)

       ELSE IF(MOD(REAL(ktau),2.0)==0.0) THEN    ! if even timestep
             bwb(:,2)=ssnow%wbtot
             delwb(:) = bwb(:,2) - bwb(:,1)
          
       END IF
END IF

There is an interesting quirk around how this water balance works with the lakes (where we keep adding water anyway)

[har917]

@ccarouge Looking at one of the problematic cases with the Haverd2023 fwsoil switch indicates that this is not varying (at least in the output). This seems to be common across all sites with this switch at true. Different cases seem to land on different values - but always constant. Depending on what value it takes we could easily end up in a position where different parts of the model want to transpire/evaporate/extract different amounts of water.

Assuming that the output is working - this seems to indicate that cbl_dryLeaf is not calling getrex_1d properly - given that's where the value of fwsoil gets set with that switch - or that there something else gone awry around the sequence to find fwsoil

Looking at the CABLE-POP_TRENDY branch and MAIN in a bit more detail indicates that a possible problem is that canopy%fwsoil is re-initialised (to 1.0) each entry into define_canopy (line 194) from MAIN but there is no equivalent line in CABLE-POP_TRENDY. Since cbl_dryLeaf with Haverd2013 uses the value of canopy%fwsoil to set the local value fwsoil at first then this difference could be introducing a consequential difference.

I don't know when that difference was implemented - likely around the refactorisation of the code.

[ccarouge]

okay - we have a problem - lines 502-514 in cable_checks

I think the easiest way around this is to replace the local re-evaluation of the total water (liquid++frozen). and hence change in water content, with that evaluated in ssnow%wbtot. This is available from soil_snow and soil_snow_gw. SLI already has it's own wbal version in cable_checks which correctly incorporates frozen/liquid densities.

If we are going to change how the wbal stuff is calculated, can we rewrite it in this way:

    REAL(r_2), DIMENSION(:),POINTER, SAVE :: owb         ! volumetric soil moisture
    REAL(r_2), DIMENSION(mp)                  :: delwb       ! change in soilmoisture

    IF(ktau==1) THEN
       ALLOCATE( owb(mp) )
       ! initial value of soil moisture
       owb=ssnow%wbtot
    END IF
           
    ! Calculate change in soil moisture b/w timesteps:
    delwb = ssnow%wbtot - owb
    
    ! Update old soil moisture to this timestep value.
    owb = ssnow%wbtot

It does the same thing and is faster to read (well, at least for me).

We could add this to this PR if we want to.

For the fwsoil value, it is worth looking into but let's not do that within this PR.

[har917]

Also - interestingly @mcuntz added canopy%fwsoil to the CABLE-POP_TRENDY restart in 2022-03-03 - indicating that this could be key.

... but I can't see how %fwsoil in the output would end up at a constant value /=0 or 1. It's being reset to 1.0 each entry to cable_canopy (on MAIN) but is then being modified by cbl_dryLeaf:getrex_1d - so how does it land on the constant value given soil moisture is varying through time?

[har917]

We could add this to this PR if we want to.

I think we have to do this in this PR - it's a related part of the whole liquid/ice densities aren't equal topic and could create confusion later on if we don't do it now.

Happy to go with your preferred implementation.

[har917]

opened new issue about the fwsoil quirk with 'Haverd2013' at #546

[ccarouge]

I'm changing the water balance calculation then.

[ccarouge]

New water balances after the fix to the densities:

station	main + Haverd2013	main + standard	494 + Haverd2013	494 + standard
FI-Hyy	-1.1594189E-02	-1.0896737E-02	-6.402957	-6.238629
US-UMB	-607.1661	-2.1411697E-03	-1.596815	-42.26907
CA-Qfo	-1.0523711E-02	-1.0954248E-02	-1.574845	-0.8182139

The fix seems to improve the cumulative water balances for this branch in general, although they are still worse than the main branch. Except for something weird at US-UMB for fwsoil_switch='standard' but considering US-UMB results on the main branch for Haverd2013, it looks like something site-specific.

[har917]

@ccarouge @JhanSrbinovsky, @rml599gh has just pointed me to #484

Just to note that the current branch includes the changes necessary to capture that issue.

[har917]

@ccarouge @JhanSrbinovsky @rml599gh 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 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.

• dotted lines mark the wilting point and field capacity for this site.

[JhanSrbinovsky]

Sorry - In and out today - what's the 4000-6000?

[JhanSrbinovsky]

fwsoil stuck at .95 is suspicious isnt it ?

[har917]

what's the 4000-6000?

time steps in the output - i.e. 0 on the x-axis relates to time step 4000 in the output. Sorry could have been cleaner.

[har917]

@JhanSrbinovsky @ccarouge @rml599gh - see #546

okay so it appears that @rml599gh was at least partially correct. I took MAIN and only edited cbl_soil_snow_main to make the LOCAL wbliq into ssnow%wbliq - this allows %fwsoil to vary through time. For the default Tumbarumba test case the water imbalance drops from -9.485mm to -6.07e-4 mm when fwsoil_switch = 'Haverd2013'

For reference, the 'standard' switch on this case gives a water imbalance of -4.42E-4mm .

However this doesn't explain what's happening with the #494 branch since I think that branch incorporates those edits to cbl_soil_snow_main.

[har917]

@ccarouge sorry this is incorrect. This updates owb before using it to evaluate delwb - i.e. delwb will always = 0

[ccarouge]

It's only on ktau == 1. so first time step initialisation. So at first time step we will have 0 but then it will work fine. Same behaviour as before.

[ccarouge]

If it helps on the review page on GitHub, it's possible to display the differences unified or split. Use the gear icon to choose (let me know if not clear on the screenshot).

[JhanSrbinovsky]

However this doesn't explain what's happening with the #494 branch since I think that branch incorporates those edits to cbl_soil_snow_main.

It does,. What's the wbal for Tumbarumba with #494 @ccarouge ?

[JhanSrbinovsky]

also what's the idea behind the wbal tracking "main" until timestep~5000 (in the above plots)?

[har917]

@ccarouge I think we need to look at the timeseries of wbal here - it could be possible that all of this imbalance is occurring in the first few time steps (which could reflect that the initial conditions may not be consistent with the new science).

Are these data still on gadi somewhere?

[JhanSrbinovsky]

I just re-read it, so the soil moisture just t takes that long to dry out I guess

[har917]

also what's the idea behind the wbal tracking "main" until timestep~5000 (in the above plots)?

In effect the Haverd2013 code is working at conserving water until the overly high value of fwsoil allows so much transpiration to occur that it drives the soil water content (layer 4) negative. At this point (timestep ~5000) the code forces soil moisture back to (near) zero values (while not decreasing transpiration) - so water conservation is broken.

Once there's some rain and the soil moisture content can naturally climb above zero again (time step ~5250) it reverts to conserving nicely at the wbal plot levels off.

With a sensibly behaving fwsoil (as provided by the changes to cbl_soil_snow_main) CABLE shuts off transpiration before this occurs and wbal is much better behaved (for this run).

[har917]

@JhanSrbinovsky In this #494 branch - is %wbliq initialised before the first call to cable_canopy and dryLeaf? If not we're passing 'garbage' into the first calls to both dryLeaf and fwsoil_calc_std if GW active.

This was in cable_parameters but that file has been removed. I suspect that we need a line in cbl_soilsnow_init_special

[JhanSrbinovsky]

ok - thanks, that was my next question - so the recovery had to wait for the rain.

wbliq init in AM3 happens. IDK about offline. It used to be there. wbliq was messed up, but the initialization was still there. Did someone remove it?

[har917]

@ccarouge @JhanSrbinovsky This is worse than it looks from the headline water balance numbers.

Looking at the time series of wbal from data on modelevaluation.org for #527 - it would appear that the end numbers picked out by Claire for these sites are actually end of year numbers related to something much worse. There is a significant annual cycle to the water imbalance at these 3 sites (both fwsoil switches) - e.g. ranging over an order of magnitude bigger than the end numbers (both signs).

BUT I'm not sure which set of numbers I've picked up from me.org - the time stamp indicates 12/3/2025 whereas I thought we'd (well @ccarouge) redone the analysis on the 14th. So it could be that this comment is out of date given the changes to cable_checks .

Also - my comment around cable_parameters was misleading - I got things confused. cable_parameters is still in the branch/repo but cable_params_mod has been deleted.

[ccarouge]

@har917 @JhanSrbinovsky Trying to look into this more systematically, I've created very rough plots for all the sites for:

• main
• main at the start of this branch
• the head of this branch
• this branch without the last fix for the water balance

The plots have timeseries of Wbal, Fwsoil and SoilMoist for all soil levels. There is a file for Haverd2013:
stations_S0.pdf

And a file for standard fwsoil_switch (only difference between the runs):
stations_S2.pdf

What I see from these is:

• this branch is fixing a problem with fwsoil_switch=Haverd2013. fwsoil used to be constant but it isn't anymore.
• the last fix on the Wbal calculation is essential to the water balance

So it looks like we have enough there to document the changes. I'm not too sure why the fwsoil_switch=Haverd2013 seems to work better with this branch changes but I'm fine with that. So we are good to go for me.

[ccarouge]

Good to go