Calculation of `wbliq` when considering different ice and liquid densities

[ccarouge]

Recently, some effort has been made in CABLE to differentiate the ice and liquid water density in soil. Until now, they have been assumed equal.

The groundwater work has done some effort in this direction as well. One difference is that it introduces:

wbliq = wb - ( density_ice * wbice) / density_liquid

Is this a change we would like to make for CABLE throughout the code? That is, we would change all calculations of wbliq in CABLE to introduce the densities.

To note that currently soil_snow has code like:

    ssnow%wbliq = ssnow%wb - ssnow%wbice

    ssnow%wbtot = 0.0
    DO k = 1, ms
      ! tot moisture this timestep 
      ssnow%wbtot(:) = ssnow%wbtot(:) +                                         &
        (ssnow%wbliq(:,k)*density_liq + ssnow%wbice(:,k)*density_ice) * soil%zse(k)
    END DO

which seems to introduce inconsistencies between wb and wbtot.

@har917 Is that something we want? I realise it can have big implications for coupled simulations. If we don't want to change this entirely in CABLE then we will have to introduce a different calculation of wbliq for the GW model which does not seem to make much sense...

[har917]

okay - a couple of points to make here - but basically this comes down to definitions and units (and I think the GW code is wrong or at least different).

In soil_snow %wb, %wbliq and %wbice are volumetric measures of water - their units are m3 (of water in any form) per m3 of soil. wb = wbliq + wbice by definition. The Richard's equations that determine the evolution of the soil moisture are written in terms of these volumetric measures of water.

The ground water equation above appears to indicate that these variables have different definitions - possibly liquid water equivalent volumes (but that doesn't work as it should be density_liq/density_ice), possibly mass of liquid water equivalent per m3 soil

In contrast %wbtot is the total mass of water in the soil column (i.e. in kg water/m2 of land area) - so the sum over soil layers of the density*wb. This is the primary measure of conservation - the change in wbtot over the time step should match inputs (infiltration) - outputs (subsurface runoff+ET) - see the wbal code in cable_checks.

[ccarouge]

@har917 thanks for the input. I confirm that I got the equation as it is in the code so density_ice/density_liq. Also, the GW code calculates:

ssnow%wb    = ssnow%wbliq + den_rat * ssnow%wbice

den_rat = density_ice/density_liq. So it is internally consistent. Although it still has wbtot = wbliq * denstity_liq + wbice * density_ice

From what you are saying Ian, we need to differentiate the calculations of wb and wbliq between default and groundwater with an IF statement if they occur outside the GW code.

[har917]

@ccarouge I'm basically not really happy to comment further without confirmation from the GW people to check/document the units they are using and (implicitly) check that their code is as the intent. I realise that's possibly not possible.

ssnow%wb = ssnow%wbliq + den_rat * ssnow%wbice with den_rat = density_ice/density_liq and wbtot = wbliq * denstity_liq + wbice * density_ice indicates that %wb is being viewed as a volumetric measure of liquid water equivalent (and that's not what should be compared to %swilt, %sfc and %ssat) whereas %wbliq and %wbice` are volumetric measures of water in that form.

---

Other things that come into this:

• offline output appears to be indicating that ssnow%wb and ssnow%wbice are in m3/m3 - see lines 611
• there are lines in the GW code where the interpretation of the terms appears counter to the above (e.g. 816 the pore space available for liquid water is %ssat minus volume occupied by ice - so 816 indicates wbice is in kg/m3 not m3/m3 as the equation for wbtot but the equation is multiplying %wbice by rat_dens=0.921 )
• a lot of the rat_dens terms are accompanied by !MMY

Edit: On reflection I'm wondering whether there are some different assumptions about how to interpret %ssat etc. - so instead of %ssat being the maximum amount of water that can be carried by the soil, maybe it's being interpreted as the maximum amount of liquid water equivalent that can be carried by the soil.

[har917]

Also - a key point - if there has been a change in interpretation e.g. of units then the code coupling GW into AM3 will need to be thoroughly reviewed.

[ccarouge]

not really happy to comment further without confirmation from the GW people to check/document the units they are using

Yep, understood. Your input was already useful.

[ccarouge]

After consulting with Mengyuan (@bibivking), it appears ssnow%wb is the liquid-water equivalent soil moisture in m3/m3 in the groundwater code. This means wb is a different quantity for groundwater and for standard soil snow (I haven't checked SLI). Although ssnow%wb is a liquid-water equivalent, the conversion of ice in soil to liquid-water equivalent was originally forgotten in various places in the code and created unphysical amounts of snow. This has been fixed by Mengyuan, hence the comments in her name. Mengyuan has only worked in making the code internally consistent, ensuring the ice content is correctly converted, but has not verified whether it is correct to use the liquid-water equivalent volumetric ratio.

From this and previous comments in this issue, it definitely appears the calculations and usage of ssnow%wb, ssnow%wbliq and ssnow%wbice will need to be differentiated within CABLE especially in parts of the code common to the ground water and the default soil-snow schemes.

[har917]

@ccarouge @bibivking This is sort of what I figured out and is a key point of difference between the parameterisation schemes. (SLI uses %wb as per soil_snow so m3/m3 of water in the relevant form I think).

We can then note that (unfortunately)

1. the equations for %wbtot and the code for %wbal in cable_checks() is not consistent with the GW interpretation*
2. the coupling equations within the UM (so %wbice and %wb to/from UM variables smcl and sthf) will need to be modified accordingly.

*probably only an issue for runtime diagnostics at the moment, but would be a problem for the rivers within the coupled model.