fixing broken lake tiles in offline #335 and addresses #347

[har917]

CABLE

Description

This PR aims to achieve two tasks

• removing the hard-wired lake index from the offline code (including GW module)
• maintaining lakes at saturation in offline model (soil-snow only)

Fixes #341
Addresses #347

Type of change

• Bug fix
• Science advancement

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

Please add a reviewer when ready for review.

---

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

[har917]

commit 4aaccb5 produces bitwise equivalence on a case where the tumbarumba test case has been converted to a lake (iveg=16) with very small canopy height (hc=0.001)

[har917]

@ccarouge @JhanSrbinovsky @rml599gh @bibivking - in implementing this fix to the offline code I've noticed that the coupled model sets the lake tile water balance at field capacity (%sfc) and not at saturation (%ssat). Can anyone advise which would be correct? (I think it should be saturation but I'm not a hydrologist)

[har917]

commit b326199 enforces saturation at the start of the time step.

The original time series of soil temperature and moisture were

This is revised to

A comparison of original vs revised soil layer 1 moisture is

noting that this is generated from output at end of time step so soil layer 1 can drop below %ssat=0.42 in this output.

[har917]

The revision also impacts the surface energy balance


and diagnostic screen level temperature and humidity

The causal relationship for some of these differences is difficult to assign - there is the immediate impact of maintaining saturation (more latent, less sensible), an immediate impact through the lower surface temperatures (less latent, less sensible and less outgoing longwave), and a delayed response because of the temporal dynamics of soil moisture and temperature (in particular the timing of the transition from <> freezing changes)

[JhanSrbinovsky]

My only comment with this is that we've just pivoted in AM3 to using a namelist input to define lakes_cable. Something which will need to be ported to offline anyway - so the is a fair solution until that happens

[har917]

I am a bit concerned about the impact on the deepest soil layer moisture seen above - and the consequent dry down of layers 4 and 5. This is because of the lines in surfbv that remove water from layer 6 to fulfill the water requirements of the surface layer. For this particular case this removal accumulates through time faster than the water can be replenished by downwards transport through the soil column and deep soil layer moisture reaches (%sfc+%swilt)/2 for large parts of the year.

Note that this removal does not respect the frozen_limit nor any frozen fraction of soil moisture implying that the code may not conserve energy (and that the code is inefficient in that lines 112-115 is redundant given line 116-119)

Obvious tests would be

• remove 116-119 or 112-115 (and use the higher limit only)
• remove 112-119 (and not allow any water to be extracted from the deepest soil layer)

and/or

• extend the approach and spread the water increment across multiple layers (say 4, 5 and 6)

and/or

• modify these lines to acknowledge that only liquid water can be extracted and that frozen_limit must be abided by.

[har917]

Looking to separate concerns - will focus this issue on implementing the same physics as in coupled model. We will leave the investigation as to whether to use %sfc or %ssat in formulations until later - noting that this may need additional flexibility to choose between formulations of %wetfac.

[har917]

Following the coupled model approach and reverting to using %sfc in the code results in

!

The revised fix still extracts substantive water (over time) from the lowest layer but the gradient in soil moisture is reduce (which should have less impact on the dynamics).

[har917]

The corresponding (using %sfc not %ssat) impacts on the variables are

There is a much reduced impact on screen level temperature (which is closer to expectation) using this approach.

[har917]

Benchcab tests ran successfully - noting that this is a largely meaningless assessment given that none of the sites should encounter the revised code. These required the check%ranges patch -
config.yaml.txt

benchmark_cable_qsub.sh.o121637899.txt

me.org is still down so sending to review.

[har917]

In the current gw_hydro code there are some conditions on %iveg .GE. 16 which implies lakes or ice. The current implementation removes the hard-wired indexes (replacing 16 with lakes_cable) but retains the .GE. condition (i.e. the assumed tile ordering).

It should be possible to implement a minor variation to the changeset here to use (veg%iveg == lakes_cable) .OR. (veg%iveg == ice_cable) as the condition and resolve this.

[JhanSrbinovsky]

sounds like a solid fix

[ccarouge]

@har917 Sorry I'm getting to this a bit late. I'm a bit puzzled with the code in surfbv. You say some of it is redundant. But that's not how I read it.

I read this code for lake in this way:

       xxx = MAX(0.0_r_2, (ssnow%wb(:,ms) - REAL(soil%sfc(:),r_2))*soil%zse(ms)*1000.0)
       ssnow%sinfil  = MIN( REAL(xxx), ssnow%wb_lake )
       ssnow%wb(:,ms) = ssnow%wb(:,ms) - REAL(ssnow%sinfil / (soil%zse(ms)*1000.0), r_2)

will potentially remove some water from the 6th layer and then

       xxx = MAX(0.0_r_2, (ssnow%wb(:,ms) - 0.5*(soil%sfc + soil%swilt))*soil%zse(ms)*1000.0)
       ssnow%sinfil  = MIN( REAL(xxx), ssnow%wb_lake )
       ssnow%wb(:,ms) = ssnow%wb(:,ms) - ssnow%sinfil / (soil%zse(ms)*1000.0)

will potentially remove more water.

So I don't see these as being redundant. Maybe in reality, there will never be enough water to draw twice and one of this will be 0 but there is nothing in the code stopping it to happen.

Is the logic here trying to get wb down to sfc or 0.5*(sfc+swilt) depending on which is smaller and some restrictions according to wb_lake value?

Unit issue?
I'm also not sure about units here. If we take the first case and wb > sfc, then we have:

xxx = wb - sfc*zse*1000
sinfil = wb - sfc*zse*1000.    (assumption on the min here)
wb = wb - sinfil/(zse*1000) = wb - (wb -sfc*zse*1000)/(zse*1000) = wb - wb/(zse*1000)+sfc

This "wb-wb/(zse*1000)" is very weird., I'm not sure the division by zse*1000 should be in there when we update wb.

I'm not sure this should be part of this fix or come as a different issue or if I'm completely wrong, so looking for your input first.

[har917]

@ccarouge A couple of quick replies:

• I don't see why we need the first

       xxx = MAX(0.0_r_2, (ssnow%wb(:,ms) - REAL(soil%sfc(:),r_2))*soil%zse(ms)*1000.0)
       ssnow%sinfil  = MIN( REAL(xxx), ssnow%wb_lake )
       ssnow%wb(:,ms) = ssnow%wb(:,ms) - REAL(ssnow%sinfil / (soil%zse(ms)*1000.0), r_2)

if we then go and do the second extraction. Since swilt<=sfc always, 0.5*(sfc+swilt)<=sfc and the net effect with/without these first 4 lines will be the same. The code removes as much water from the 6th soil layer as required to satisfy %wb_lake without forcing %wb below 0.5*(sfc+swilt) either way; these four lines just force it to do it in 2 steps not 1.

• the factor 1000 here is a density_liq. We need to conserve water while noting that sinfil and wb_lake are in kg/m2 whereas wb is in m3 water/m3 soil (and sinfil in particular has other dependencies). You're correct that the conversions are meaningless if you're at the limit on xxx, but they're not if you're not at the limit.

• there is a problem here in that none of this respects any frozen water in the soil (i.e. lake) or the frozen limit - but that's a separate issue and highly unlikely to be encountered.

[ccarouge]

All good.

[har917]

@ccarouge one final comment on units: you said you were concerned about units in

xxx = wb - sfc*zse*1000
sinfil = wb - sfc*zse*1000.    (assumption on the min here)
wb = wb - sinfil/(zse*1000) = wb - (wb -sfc*zse*1000)/(zse*1000) = wb - wb/(zse*1000)+sfc

You missed a bracket - it should be

xxx = (wb - sfc)*zse*1000
sinfil = (wb - sfc)*zse*1000.    (assumption on the min here)
wb = wb - sinfil/(zse*1000) = wb - ((wb - sfc)*zse*1000)/(zse*1000) = sfc

so everything works out