soilsnow developments from AM3

[JhanSrbinovsky]

CABLE

Description

Developments from AM3 work. Crazy number centre around using an ice density in place of liquid density, or 0.9. * liquid density

iThis is why I left this til the end as there will be differences resulting from these mods

Fixes #494

Type of change

• model improvement

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--497.org.readthedocs.build/en/497/

[JhanSrbinovsky]

benchcab WILL reveal changes to output. In the process of downloading/uploading. My 3 site local test on gadi is pretty reasonable. Will post link to ME.org when done

[JhanSrbinovsky]

still sorting out the ME.org plots. My testing shows seasonal diurnal quantities at Tumba (as expected) negligible. One of the sites I test at is Hyttialia (because there is snow): Note seasonal labels appropriate for Sothern Hemisphere

[JhanSrbinovsky]

I think it looks as expected, I'll do a global one and see how it goes (I have to find the code first!)

[JhanSrbinovsky]

This is using GSWP2 forcing. Most of the plots that flashed in front of my eyes as the scripts kept turning over were unremarkable I thought. But the surface SoilTemp caught my eye as it is one that it make sense to do a relative difference as it never gets near zero etc. As with all the others it is pretty bland, except for this bright spot in Northern Russia. On closer inspection this bullseye represents a ~10% diff. But it goes away after a few months. Odd, I wonder what could be causing this? @har917 any ideas? Greenland is slightly cooler in Summer but not nearly as shocking. The Russian 10% difference appears to be an initialisation issue as it doesnt come back towards the end of the year. But still what is so dramatically different here. There was nothing here in AM3 was there?

[har917]

@JhanSrbinovsky could be a feature of the initialisation - could be something to do with the ice fixes.

I'd be cautious about looking at %change in soil temperature - there's little context of a 10% difference in temperature - depending on the units (C vs K) you can identify different regions for concern. e.g. C will be particularly susceptible to divide by near-zero regions, using K will results in even large (+3K) differences only showing up at the 1% level.

[JhanSrbinovsky]

@har917, its only really the ice fixes that are in here. but why should they apply to this bullseye in particular. I just assumed it was Kelvin as almost everything in the UM is K. If it is in C, Im not so surprised. I just assumed this 1 spot was experiencing a ~30degree cooling, so probably C. I'll check

[JhanSrbinovsky]

damn it. its in Kelvin

[JhanSrbinovsky]

@har917 This lazy screenshot shows that almost everything is within -1% to +1% diff. That might be 3K and I can't comment on whether or not that is acceptable (or even expected) BUT the bullseye in Russia is 10 times that. Is this something to worry about, or should we move on?

[har917]

unfortunately I think we need to understand the origin of a ~30K difference in soil temperature - the most likely cause is that something that is now being used is getting through uninitialised (or to a default value of 0).

[JhanSrbinovsky]

I dont have any ideas about what to look at. I might just look at the usual suspects and see what the moisture is doing at the same point. what else?
because it is isolated, and seems to go away by the next year I am hoping it is just a dodgy initialisation. Keep in mind as well that this is a no feedback, forced simulation.

[ccarouge]

Since Ian says there is some digging to do to understand the differences, I've removed the request for review. This way I'm not coming back in here constantly to see if it's ready. @JhanSrbinovsky just ask for reviews again when ready.

[JhanSrbinovsky]

Trouble spot zoomed in (rel. Diff ) T(K)



SoilMoist - Diff


SoilMoistIce - Diff



Iveg



Rather than show you the dozens of field which are unremarkablekable:
snowdepth

Curiously the snow water equivalent does not show this. A clue that I can’t make sense of. I can imagine an initiallized snow depth which is all of a sudden subjected to a different density, but this wouldn't be the only place

[JhanSrbinovsky]

This is interesting, the albedo:

[JhanSrbinovsky]

are we fortuitously hitting a limit?

[har917]

Jhan - can you check snow density? There is a link between snow density and albedo (and a limit on that impact) in surface_albedosn()

This is a region with a lot of (frozen) lakes - remember we changed how met%tk/ssnow%tggsn got used - so that could be a factor as well.

[JhanSrbinovsky]

sounds promising - I'll have to re-run both models - it isnt in my output. I imagine snow density output will be available through an easy switch

[JhanSrbinovsky]

bugger - it isnt!! it is available in a restart but this isnt much good to us. I'll come up with something a bit later.

[JhanSrbinovsky]

ok im back now. Does anyone know if it is possible to change the dump frequency in the offline model? This would be the easiest way I can think of to get a single point out of a global offline run?

[JhanSrbinovsky]

Hey @SeanBryan51, you've recently been looking at the _drivers for offline. Is writing the restart coded after the timestep loop?

[SeanBryan51]

Hey @SeanBryan51, you've recently been looking at the _drivers for offline. Is writing the restart coded after the timestep loop?

Apologies, I haven't looked into the restart IO in the drivers so I'm not too familiar on where the logic is in the code. This has not been touched yet so the restart behavior should be unchanged

[JhanSrbinovsky]

I've got a better handle on things myself, the ssatcurr limiting "fixed" things slightly. The rest of the globe has only been affected by a fraction of a percemt. Even here, the biggest change is only ~4%, when the cooling was of order 10%. At any rate the initial temps seem quite high in these regions. Whatever, that could be the case. The crux of the discrepancy between main:494 is that main: broadly continues this pattern, whereas the adjustment to tgg due to the latent heat of melting excess water/ice effectively "cools the hot spots". The region is more homogenous but IDK that this is necessarily what we want. By hot spots I mean these cells are initially at ~273K where as surrounding region is commonly 20-40 below. I've already put a snippet of the code above. The bit with dfactor. How do you want proceed @har917 ?

[JhanSrbinovsky]

The same region - NO ssat limiting Straight K -K difference

and global

[har917]

@JhanSrbinovsky I'm trying to unpick what you're showing and asking

By the 'ssat_curr' fix do you mean changing line 425(ish) in cbl_smoisturev.F90 from

ssatcurr(:,k) = soil%ssat - ssnow%wbice(:,k)

to

ssatcurr(:,k) = max(soil%ssat - ssnow%wbice(:,k),0.0)?

While this appears sensible I'm struggling to see how we get negative ssatcurr - i.e. there's something else going on which has allowed %wbice > frozen_limit* %ssat (equivalently %wb > %ssat)

I'm also struggling to understand what has been changed from the runs where you had a 30K difference to these latest runs where the differences appear to be more modest (up to 3K). could you restate?

[JhanSrbinovsky]

The two plots I just showed do NOT have the ssat range limiting, BUT yes it is like below:

ssatcurr(:,k) = max(soil%ssat - ssnow%wbice(:,k),0.0)?

While this appears sensible I'm struggling to see how we get negative ssatcurr - i.e. there's something else going on which has allowed %wbice > frozen_limit* %ssat (equivalently %wb > %ssat)

this I have no insight into, I clobbered it here and changed tac

I'm also struggling to understand what has been changed from the runs where you had a 30K difference to these latest runs where the differences appear to be more modest (up to 3K). could you restate?

ALL the AM3 developments in soilsnow are included EXCEPT the code which increments tgg by the LE required to melt water above frozen limit as in the code in smoisture :

! If redistribution of liquid has caused a layer to exceed frozen limit
! melt the excess and account for latent heat.
dfactor = 1.0 - Cdensity_ice / Cdensity_liq
DO k = 1, ms
   DO i = 1, mp
      IF (ssnow%wbice(i,k) > frozen_limit * ssnow%wb(i,k)) THEN
         sicemelt = (ssnow%wbice(i,k) - frozen_limit * ssnow%wb(i,k)) / &
                       (1.0 - frozen_limit*dfactor)
         ssnow%wbice(i,k) = MAX( 0.0_r_2, ssnow%wbice(i,k) - sicemelt )
         ssnow%wb(i,k) = ssnow%wb(i,k) - dfactor*sicemelt
         ssnow%tgg(i,k) = ssnow%tgg(i,k) -                                 &
                 sicemelt*soil%zse(k)*Cdensity_ice*CHLF / REAL(ssnow%gammzz(i,k) )
      END IF
   END DO
END DO

And then it is only a problem at a few points globally. And then also, it seems to go away after spinning up for a while. I will see if I can spot it in the ESM16(-ish) output.

I am having trouble reading it. I considered that it was just my misunderstanding or not knowing a missing piece of the puzzle, especially because it seems to work elsewhere in the globe. It seems to be kinda the reverse of what is in soilfreeze. It did occur to me though that the surrounding IF likely excludes almost everywhere anyway. I dont know how it gets in here. I haven't verified but almost certainly if these few lines alone are added to main, we'll get the same 30K discrepancy. Was it needed for the energy balance in this region in particular. It could be a PFT(=lakes) thing. I plotted iveg above but its not clear.

[JhanSrbinovsky]

Yeah - definitely Tundra.

[JhanSrbinovsky]

i.e. there's something else going on which has allowed %wbice > frozen_limit* %ssat (equivalently %wb > %ssat)

Given that this is expected to be the case, then shouldn't getting into the loop in the first place be excluded. I wonder how many times it gets into the loop and there isnt a problem.

Anyway, no sign of it in AM3. Curiously AM3 doesnt have these near 0 C Temps in the region. Its all 20-40 below.

[har917]

@JhanSrbinovsky to pick up on a couple of things in here

• we need the %tgg = %tgg - sicemelt*soil%zse(k)*Cdensity_ice*CHLF / REAL(ssnow%gammzz(i,k) ) line in there. We get into that loop on occasions where evaporation/drainage reduces the liquid part of the soil water (in surfbv) so that %wbice > frozen_fraction * %wb. It is is more likely to encounter that case now that we have different densities for liquid and frozen water.
• we shouldn't have been getting to %wbice > frozen_limit*%ssat
• I would not be surprised if there are a small number of grid cells where the initial conditions in gridinfo are incompatible with the updated science. Such cases will likely show up on the first time step.
• we've been a lot more careful around initialisation, especially around the application of frozen_limit, in the AM3 codebase - so again not that surprising that something has been triggered.
• I'd be happy if the new code ends up smoothing out any 'hotspots' or 'coldspots'

---

I'm not really sure about the way forward - it seems highly unlikely that we can create a codebase-gridinfo combination that preserves the previous offline performance while also not adding lots of adhocness to the codebase.

The most scientifically credible way forward is to either

1. include the %tgg line, take the initial temperature hit, run out for several years and then use that to create a new gridinfo file; or

2. do two phases: i) run a month or so with the %tgg line copied out (or wrapped in a if first_call condition), write to restart, then ii) add the %tgg line back in, rerun from restart for a while, then create a new gridinfo file.

3. is probably easier but does run the risk that the model doesn't come back to the same approximate climate. Either way we break comparability in a benchmarking sense. @ccarouge - any thoughts?

[JhanSrbinovsky]

@har917

• we need the %tgg = %tgg - sicemelt*soil%zse(k)*Cdensity_ice*CHLF / REAL(ssnow%gammzz(i,k) ) line in there.

This being the case then I am in favour of below:

1. include the %tgg line, take the initial temperature hit, run out for several years and then use that to create a new gridinfo file; or

As overwhelming and dramatic as it may seem to re-create the grid info file

Tthe grid info file is decades old anyway and no-one is really confident about where all the bits and pieces came from.

Presumably we could correct this.

We should take this opportunity to establish an offline grid-info file that is reflective of the ACCESS online runs. There would be reasons to keep it in ACCESS resolution, however it might be better to adopt a more generic 1*1 degree resolution.

It would be satisfying ti know what it is in the grid info file that is messing with us?

[JhanSrbinovsky]

So this has been run out for 30-50 years in AM3. Is there someone at NRI that can create a grid info file @ccarouge ?

[ccarouge]

@JhanSrbinovsky for the gridinfo file the only things we are working on are creating the file from UM output/restarts and creating the file using ANTS (that's far from being ready). Neither of these options is what you need here. So no, we don't have anything for what you want to do.

[JhanSrbinovsky]

Didnt @Whyborn not long ago help @rml599gh with extracting a grid info file from ESM1.5?

[rml599gh]

A description of the gridinfo file I built based on ESM1.5 is here: https://forum.access-hive.org.au/t/cable-site-runs-with-access-forcing/2108/8
I've only ever tested this in a single site set-up so have no idea if it would run a global case.

[Whyborn]

I have built such a gridinfo file, based mostly off of Rachel's method in the forum post she linked. I ran a few spin-up stages of CABLE offline (on the CABLE-POP_TRENDY branch) using an "alpha" version of the gridinfo file, and it seemed to work fine. You can find the conversion script at convert_from_fields.py, which should just work out of the box in the hh5 conda environment, assuming you have access to rp23, fs38 and vk83.

Whether this gridinfo is something we want to use longer term, I can't say.

[JhanSrbinovsky]

OK thanks @Whyborn , @rml599gh

@Whyborn - I might try to use your grid info file as a test in this GSWP2 fiasco here. Can you please point me to it?

[Whyborn]

The gridinfo I've run with is at /g/data/rp23/experiments/2024-03-12_CABLE4-dev/lw5085/CABLE-as-ACCESS/gridinfo_files/ACCESS-ESM1p5-1p875x1p25-gridinfo-CABLE_v0.1.nc. This was on the CABLE-POP_TRENDY branch- I don't think there are major differences in the required gridinfo between main and CABLE-POP_TRENDY?

[har917]

@JhanSrbinovsky @Whyborn

I don't think there are major differences in the required gridinfo between main and CABLE-POP_TRENDY?

Probably not - however I think CABLE-POP (so TRENDY and BIOS applications) routinely uses a 4 tile gridinfo whereas as main will expect a 17 tile version. I don;t know the offline i/o code well enough to advise whether this is an issue or not.

[JhanSrbinovsky]

Unfortunately it doesnt like it.

I'll see if I can use your python script. a couple of questions @Whyborn ,

1. (out of) Where did you get the fields mentioned
2. Was you CABLE run global?

Nan in evap flux, 5219 65.50000 -17.50000
fe nan 5219 2 2.8106968E-03 0.0000000E+00 0.0000000E+00
220.3700 0.0000000E+00 0.0000000E+00 269.2953 7.291198
NaN 946.4033 NaN NaN NaN
NaN NaN NaN NaN

[JhanSrbinovsky]

reading the script it seems like an ESM1.5 dump is where you were grabbing things from, which is unfortunate as this is where I was thinking anyway, but it doesnt like something :(

[Whyborn]

It might be to do with iveg? I'm a bit confused- do we treat ice as bare ground in CABLE offline? The default CABLE gridinfo appears to set most of Greenland and high arctic regions to be bare ground, while ACCESS-ESM sets those areas to the actual ice type (or at least appears to, based on the ancillaries). Do we have valid handling for the ice tile in offline?

[JhanSrbinovsky]

Really, bare ground. Im not familiar with offline, I assumed it was treated as ICE (i.e. iveg=17). Im going to check GSWP, Im sure it was ICE but could be wrong. Maybe is just a TRENDY thing.

[JhanSrbinovsky]

yeah GSWP treats Greenland as ICE and dodges the poles completely

[rml599gh]

Not sure about Greenland but I don't think offline typically runs for Antarctica so the 1x1 gridinfo doesn't go that far south.
When I was making my access gridinfo, I did treat ice points differently. Even though there are values for the soil parameters in the restart file I didn't use these but overwrote them with the values for isoil=9 from def_soil_params.txt. I think this is what happens in ACCESS-ESM1.5 (another quirk that it would be good to clean up). Hence my code:

           if (iveg(i,j).eq.17) then
                   isoil(i,j) = 9
                   rhosoil(i,j) = 910.
                   css(i,j) = 2100.
                   bch(i,j) = 7.1
                   clay(i,j) = 0.3
                   css(i,j) = 2100.
                   hyds(i,j) = 0.000001
                   sand(i,j) = 0.37
                   sfc(i,j) = 0.301
                   silt(i,j) = 0.33
                   ssat(i,j) = 0.479
                   sucs(i,j) = -0.153
                   swilt(i,j) = 0.216
           else
                   isoil(i,j) = 2
                   rhosoil(i,j) = 1600.
                   css(i,j) = 850.
           endif

[har917]

This issue of ensuring consistency between iveg and isoil and how to fill in the soil/pft values in the offline appliciation has been a focus of some recent work (linked to GW model implementation) - see #338.

I think we're okay offline but the initialisation in ESM1.6 may need to be looked at again.

[Whyborn]

@JhanSrbinovsky I've made some fixes to the gridinfo generated from the ACCESS-ESM1.5, and now it seems to run for a few select points without producing NaN. Had a quick look at GPP and Ebal and the results seem reasonable? The new version of the gridinfo is at /g/data/rp23/experiments/2024-03-12_CABLE4-dev/lw5085/CABLE-as-ACCESS/gridinfo_files/ACCESS-ESM1p5-1p875x1p25-gridinfo-CABLE_v0.3.nc.

[JhanSrbinovsky]

Thanks @Whyborn , I'll give it a shot

@Whyborn Unfortunately I got the same error. :(

[har917]

these are new terms from AM3 development. they are the counterpoint to the melt terms to %dtmlt and come about by the freezing of liquid precip onto frozen ground or snow pack.

good that they've been picked up - but not that surprising given they sit in a different part of the code base to the rest of the ice fixes