Losing Dense Shelf Water in ACCESS-OM2-01

[AndyHoggANU]

Over recent months we have been investigating a problem with ACCESS-OM2-01 that caused Dense Shelf Water near the Antarctic to disappear. It would appear that we have diagnosed the source of the problem - this issue is designed to document the problem and to start a discussion on the best way to solve it.

The issue can be most easily diagnosed by averaging salinity in the SW corner of the Ross Sea at 350 m depth and over a 4°x5° box. A timeseries of this diagnostic is shown in the plot below. Here, we are comparing 3 runs, branching from our original simulation (_spinup6). These runs are

• _spinup6_newexe: using updated executables (from March 2019)
  • _spinup6_redux: Using the original executables (from 2018)
  • _spinup6_cice_d5b4d4d: Using the new executables, but reverting a change in ice_therm_itd.F90.

This plot makes it fairly clear that the change in ice_therm_itd.F90 is the culprit. To see this change, see here: https://github.com/CICE-Consortium/CICE-svn-trunk/commit/faa71893593ce88e6e96711d20cac5ff74cca2d8

The change we made in _spinup6_cice_d5b4d4d is here: https://github.com/COSIMA/cice5/commit/d5b4d4db399a6ecdd4707c6cda740c0152c707a0

I think the intent of this change is that when we use ktherm=2, we really should set update_ocn_f = false. However, note that this would still involve a code change from the original code.

Alternate possibilities are to:

• Go back to ktherm=1 (we would like to do this, but I think it is unstable at 0.1°)?
• Reinstate the original code, or rewrite the code.

Any suggestions?

[ofa001]

Hi Andy, we stayed with ktherm=1 for the coupled model, the Met office coupling of it multi layer thermodynamics is only set up for this option, though they originally wanted to extend to mushy thermodynamics but have now switched to the European ice model.

I think you originally had problems with ktherm-=1 in the fresh parts of the Arctic ( Baltic, Hudson Bay) rather than the Ross sea, didn't Russ also introduce another fix for salinity getting too fresh under the ice (its too far back to remember it exactly, I think I mentioned the change the met office used as well to you for the same stability problem ( can put you in contact with the person behind the correction if you like).

[russfiedler]

I think MOM really wants update_ocn_f=.true.. My reading of the code's intent is that POP automatically makes the change to salt and freshwater internally based on the freezing potential. Therefore, POP doesn't want to include the salt and water fluxes due to frazil in the ice model apart from a left over contribution in the mushy ice thermodynamics which gets sent back to the ocean.

Have a look here for what POP does. www.cesm.ucar.edu/models/cesm1.0/pop2/doc/sci/POPRefManual.pdf

pages 106-108. It makes it pretty clear that the salinity and freshwater changes are made in situ.

[AndyHoggANU]

Hi @russfiedler - are you saying that we want update_ocn_f=.true. but that we need to reinstate an option for what to do if ktherm=2?

Or are you saying that the original code was incorrectly double-counting salt fluxes?

[ofa001]

Russ that description in the POP manual is very similar to what Dave (Bi) copied over in to the 'auscom' set up, the names are even the same, I need to follow up on this Update_ocn_f = true option. The freezing potential is an intrinsic part of how CICE works.

[russfiedler]

Siobahn, yes, the frazil code in ACCESS-OM has been modified not to make changes to salinity (unless requested) and was based on code written before POP started altering the fresh water. If we alter the salt in the ocean then we would have to stop the salt flux from CICE but allow the freshwater. How we would deal with the left over bit from the mushy thermodynamics is something I don't want to think about...

[russfiedler]

@AndyHoggANU I'm pretty sure that the changes in https://github.com/COSIMA/cice5/commit/d5b4d4db399a6ecdd4707c6cda740c0152c707a0 are double counting in the case of ktherm=2.

[nichannah]

Do we have a run (or part) with ktherm=1 that we can run this diagnostic on?

[aekiss]

@nichannah these use ktherm=1: /g/data3/hh5/tmp/cosima/access-om2-01/01deg_jra55v13_ryf8485_spinup8/output* /g/data3/hh5/tmp/cosima/access-om2-01/01deg_jra55v13_ryf8485_spinup6/output00[0-2] there's not much with ktherm=1 at 0.1 deg since this tended to abort due to a failure to converge the thermodynamic iterations.

[AndyHoggANU]

Hi @nichannah, I can try to run a case with ktherm = 1 today if that's useful?

In the meantime, I thought the best I could was to look at the double counting issue more directly. The best way I know of to do this is to look at salt_global_ave. If this is steady, then we mustn't be double counting, right @russfiedler ??

So, two plots following. The first shows globally averaged salinity from the start of spinup6 (our original spinup). Alongside this I'm showing spinups from a series of recent attempts, all with the new executables:

I think this shows that the original scheme was doing the right thing, and that the changes induced by the new bit of CICE code was breaking our salinity budget. I am assuming this was because of the fragment of code we are discussing. In support of this, here is the global salt balance from the most recent runs. I'n showing the original spinup6 as well as the _redux case. The _newexe run has the latest CICE code, while the _cice_d5b4d4d code reverses the change in ice_therm_itd.F90:

I think these runs demonstrate that the COSIMA/cice5@d5b4d4d changes are doing the right thing. (I, for one, am a bit relieved by this, but if you think I am wrong please weigh in.)

I did add one new case to our repertoire overnight. This is labelled _spinup6_update_ocn_f in the plot above, and it keeps ktherm=2 but sets update_ocn_f=.false.. You will see from the code that this adds the vi0tmp correction to dfresh. You will see that this partly corrects our salt imbalance, but salt still seems to decrease a little.

Comments welcome. In particular:

• Any more diagnostics you would like to see?
• Any more runs I can do to tease this out?
• Should I run a test with ktherm=1? (@aekiss - I think I may have to start from zero time to do this, right??)
• Any suggestions for what we should put in the code to fix this?

[russfiedler]

@AndyHoggANU I have a feeling that we're missing something subtle with the ktherm=2 case. The note in the code says that the freshwater and salt fluxes for the case update_ocn_f are added somewhere else. The thing is, where? Your plots make it pretty clear that either salt isn't being conserved. I wonder if there is some new code that wasn't included or if it's handled by preprocessor directives so that we're missing it.

[AndyHoggANU]

Good point @russfiedler. But is it possible that the new CICE code just accounts for the use cases that they know about, and that we are using a combination of parameters which they did not consider?

One possibility here, suggested by @kialstewart, is that we run test cases at 1°. We could do this with both ktherm values, while using less computer time, etc. It seems this is a worthwhile exercise, so will try to do this unless anyone objects.

[russfiedler]

@AndyHoggANU Yes, I was coming to a similar conclusion myself. I can't find anywhere in the code (cice5 or cice6) where the update_ocn_f is used. It's like they were planning to do something but didn't get around to it and the only case that works with mushy physics currently is with POP. I think the correct version of that block (for the moment) is just to remove the conditional on ktherm <= 1 and run with update_ocn_f=.true..

  if (update_ocn_f) then
        dfresh = -rhoi*vi0new/dt
        dfsalt = ice_ref_salinity*p001*dfresh
        fresh  = fresh  + dfresh
        fsalt  = fsalt  + dfsalt
  else ! update_ocn_f = false
     if (ktherm == 2) then ! return mushy-layer frazil to ocean (POP)
        vi0tmp = fnew*dt / (rhoi*Lfresh)
        dfresh = -rhoi*(vi0new - vi0tmp)/dt
        dfsalt = ice_ref_salinity*p001*dfresh
        fresh  = fresh + dfresh
        fsalt  = fsalt + dfsalt
        frazil_diag = frazil - vi0tmp
     ! elseif ktherm==1 do nothing
     endif
  endif

``

[AndyHoggANU]

Hi @russfiedler - I would be happy with that. Does anyone else (@nichannah, @aekiss, ?) want to weigh in on this before we make the call?

Also, should we be pushing this change back up the main CICE trunk?

[russfiedler]

I think we need to find out what added elsewhere means!

[aekiss]

Nice work pinning this down. I'm very relieved our previous runs were OK.

@russfiedler's suggested code change looks good to me. I agree theupdate_ocn_f stuff seems half-finished. There's nothing else relevant in the commit in which it first appears: https://github.com/COSIMA/cice5/commit/5b3d0b6bb74a00364bbeaa1140d44b294ed48f9a

[nichannah]

Can we just ask Elizabeth what "! elseif (ktherm == 2) the fluxes are added elsewhere" means?

Rewording to try to make things clearer.

update_ocn_f = .true. means that the ocean model does not make an assumption about the fresh and salt content of frazil and therefore needs to be told about this via the coupling fields when there is frazil growth.

update_ocn_f = .false. is the POP way. Fluxes from frazil growth are not included in fresh and salt. It will need to calculate these itself based on frazil growth (I guess).

ktherm = 2 the fresh and salt content of frazil is accounted for in the calculated frazil growth. ktherm = 0, 1 reverse of the above. i.e. calculated frazil growth does not include fresh and salt.

If update_ocn_f = .false. and ktherm = 2 (POP) the ktherm calculation uses it's own assumptions about salt and fresh amounts in frazil, this may be inconsistent with the ocean model assumption therefore the fluxes need to be split and passed separately.

If update_ocn_f = .true. and ktherm = 0,1 no fresh and salt included in frazil growth so tell ocean model what this should be.

It's a long shot but the new code (maybe) makes sense in the context of just calculating the change in frazil, fresh and salt related to ITD (ice thickness distribution) - regardless of which ktherm option is used. With ktherm=2 the fresh and salt fluxes are not changed by ice_therm_itd - only frazil growth is. It is then the responsibility of some code outside this calculate fresh and salt fluxes based on frazil growth.

A question about the old (now re-instated code). Is this corrrect anyway for ktherm = 2? It is not taking into account the actual salinity of the frazil ice which ktherm = 2 is meant to do - instead it uses ice_ref_salinity.

[nichannah]

I've decided that I agree that this change is a oversight and we should re-instate the old code. However we should make sure that for kthem=2 the salt is balanced between frazil formation and melt given that this option allows for variable salinity frazil.

[ofa001]

I have now spotted that you pulled this section of code in from CICE6 (that's what puzzled me earlier Andy). I think Nick is partly right. Chasing pack the meaning of the (vi0new-vi0tmp){rhoi/dt} that is being returned to the ocean. vi0tmp=-fnewdt/rhoiLfresh whilst vi0new (when ktherm=2)=-fnewdt/enthalpy_mush( Ti, SiOnew).

So the variable salinity is passing indirectly into these terms, of the mushy ice layer is passing into the frazil, that the ocean is forming. The ocean should be forming frazil based on its ambient salinity and only when it joins on to the ice pack will it start taking on the profile of the existing mushy ice salinity profile (and drains its brine), if it adds to the ice bottom or more usually frazil forms in leads. so adds to the thinnest ice thickness category. How that fits into this mushy ice formula, I am not sure, I am no expert on mushy ice. And as Russ said above its complicated, but the salinity issues is also one for the ocean model frazil scheme.

[aekiss]

Is there anything left to follow up on this or can it be closed? (has anyone checked that salt is balanced between frazil formation and melt with ktherm=2 as @nichannah suggested?)