Spurious supercooling in ROMS-CICE simulation

[k.alexander]

My name is Kaitlin Alexander and I am a PhD student in Australia, working with Katrin Meissner, Ben Galton-Fenzi, and Matthew England. As part of my PhD I am using the new MetROMS (ROMS + CICE + MCT) system coupled by the Norwegian Meteorological Institute. I have adapted this to a circumpolar Antarctic domain (quarter-degree resolution, scaled with cos of latitude) and merged in Ben's ice shelf code.

Something I have been struggling with for months now is the issue of super thick sea ice. After 3 years of simulation the sea ice is around 4 times as thick as ICESat measurements, with some coastal areas as thick as 80 metres (!) This means I can't keep polynyas open, and the lack of strong localised dense water plumes causes all sorts of other problems.

After going through just about every part of the CICE code, I concluded the problem was most likely with ROMS, specifically the supercooling which ROMS calculates and passes to CICE to form new ice. This code is based on the POP code for coupling with CICE: the ocean integrates the net supercooling over the upper few layers, sets the supercooled temperatures back to the freezing point, and sends the energy to CICE in W/m^2. Since POP is a z-coordinate model but ROMS is terrain-following, rather than integrating over a specified number of layers, MetROMS integrates over a specified surface thickness (I'm using 5 metres) for which the number of terrain-following layers can vary.

I was checking the original POP code and found this warning:

Code: Select all

!  WARNING: unless a monotone advection scheme is in place, 
!  advective errors could lead to temps that are far below freezing
!  in some locations and this scheme will form lots of ice.
!  ice formation should be limited to the top layer (kmxice=1)
!  if the advection scheme is not monotone.

I don't want to limit ice formation to the upper layer, because this is very thin over the continental shelf (10-20 cm) but much thicker (>5 metres) offshore. I don't think it would be a good idea to have wildly different thicknesses for ice formation depending on the region.

So I started looking for evidence of spurious supercooling, and this is what I found. Here is a slice through the Ross Sea with no time-averaging or spatial interpolation/averaging. It shows the difference from the freezing point (as seen by the supercooling code, Tf = -0.054*salt, depth dependence shouldn't matter over the upper 5 metres):



There is definitely some spurious supercooling going on over the steeper areas of the continental shelf, far from the ice shelf cavity (which deals with supercooling separately). The spurious mixing is even more obvious with the absolute temperature:



and to some extent with salinity:



The velocities u and v look better, I think any weirdness there is largely just in response to the tracers, but I might be wrong:





If the spurious supercooling happens below 5 metres it won't directly form any ice (although I'm wondering if it would indirectly form ice by cooling the upper layers or just generally causing more instability?). But if it happens within 5 metres of the surface it won't show up in these plots because it's immediately set to the freezing point and sent to CICE to form new ice. I wasn't sure if this would even happen at the surface or if spurious mixing was largely a problem at depth. But then I compared plots of the positive freeze/melt potential (i.e. energy from supercooling sent to CICE):



to the bathymetry:



and if you line them up, the highest areas of supercooling are right near the continental shelf breaks of the Ross and Weddell Seas. In reality frazil formation should be concentrated right by the coast or ice shelf fronts, where there should be polynyas. The only reason supercooling would be strongest at the continental shelf break, as far as I can see, would be advective errors due to the steepness of the bathymetry there.

These are the CPP options I am using for advection and related things (eg mixing) right now:

Tracers: TS_C4HADVECTION, TS_C4VADVECTION, MIX_ISO_TS, TS_DIF2, DIFF_GRID
Momentum: UV_ADV, MIX_S_UV, UV_VIS4, VISC_GRID

and I have Vtransform=2, Vstretching=2.

By looking through other posts on this forum I think the culprit is TS_C4HADVECTION and TS_C4VADVECTION. I am writing to ask what would be better schemes to use. It looks like TS_MPDATA is the only monotonic advection scheme available, but Ben mentioned it might have other weaknesses. I have also seen posts suggesting TS_U3ADV_SPLIT or TS_A4HADVECTION/TS_A4VADVECTION are less prone to spurious maxima and minima, even though they aren't monotonic.

I can't test any of these until 1 October (when the CPU time allocation for next quarter becomes available) but I thought it would be wise to ask for expert advice in the meantime, rather than going down a path that is likely not to work. What do you think would be my best option(s)?

Many thanks,
Kaitlin

[kate]

I just got the metroms code running last week and I too have too much ice. I have a smaller test case of just the Bering Sea, starting with no ice in September and running for a year to watch the ice grow and retreat. It doesn't retreat anywhere near enough, unlike my ROMS-CICE which retreats a bit too quickly (as does the Budgell ice).

Seth Danielson and I have been talking about our Arctic sea-ice troubles and I just spent the day adding more diagnostics to the ice thermo (Budgell ice). We want to be conserving energy and the model wasn't giving us enough information to know what was going on. I'll be leaving town in a few days and may not get it working by then, but we'll eventually figure it out. Seth also has a student struggling with ice thermodynamics (Kofan Lu).

[k.alexander]

Hi Kate, that's interesting that you have the same problem with MetROMS but not your internally coupled ROMS/CICE. It's not the frazil (supercooling) code because one of the things I tried was using your frazil code instead of the (very slightly different) MetROMS code. Were you using the same advection schemes for both setups?

[kate]

Yes, just these:

TS_U3HADVECTION Third-order upstream horizontal advection of tracers.
TS_C4VADVECTION Fourth-order centered vertical advection of tracers.

No explicit horizontal diffusion. I've got more ice in the northern Bering Sea than I did a year ago with the bad old ROMS ice - there's go to be something more wrong than an advection scheme can account for.

[k.alexander]

The one thing I didn't try from your frazil code was ICE_LOG_LAYER. I don't suppose that would make a difference?

Regardless of whether advection is *the* problem, it clearly seems to be *a* problem for my setup at least...I tried TS_U3HADVECTION ages ago and it absolutely destroyed my sea ice everywhere except the Ross Sea, where there was way too much of it. Not sure what was going on there.

[kate]

I haven't tried ICE_LOG_LAYER since fixing several other things. It's more like what the Budgell model is doing, so shouldn't be that big a difference.

[msd]

Hi Kaitlin,

For what it's worth, I use the same tracer advection schemes as Kate (TS_U3HADVECTION, TS_C4VADVECTION) for my circum-Antarctic ROMS model with Budgell ice and I don't seem to have any issues with the advection scheme making the ice go away everywhere (except the Ross).

[k.alexander]

I've fixed a bunch of things since I last tried TS_U3HADVECTION. I'll try it (and a few other advection schemes) next week when I have CPU time again.

Kate, are your frzmlt values very different between the ROMS-CICE "fake coupling" and MetROMS? Or is it just the ice thickness?

[kate]

My frzmlt looks odd, so I'm rerunning after a tweak. Not sure when I'll be able to check on it though.

[kate]

OK, I think I fixed yet another bug. I am now melting a little less ice than I was before with my CICE coupling, still melting far more than with metroms. metroms gets more congelation ice growth in the spring.

I'm also fussing with the shortwave in the Budgell ice model and need to run a year of the Arctic domain with both models (again).

[k.alexander]

I ran 4 months with TS_U3HADVECTION (instead of TS_C4HADVECTION) and the thick ice is gone! No more signs of spurious supercooling, and no more crazy thick ice over the continental shelf break.

It remains to be seen whether this will fix my lack of coastal polynyas (and resulting erosion of Winter Water) or whether the sea ice will eventually disappear everywhere except the Ross (as it did last time I tried U3HADVECTION). I'll run a few years and see.

I also want to try TS_U3ADV_SPLIT, TS_MPDATA, and TS_A4HADVECTION/TS_A4VADVECTION. If TS_U3HADVECTION behaves I also want to try removing explicit mixing (right now using TS_DIF2) as several people have told me this isn't necessary with TS_U3HADVECTION. Would I remove explicit mixing of momentum as well (right now using UV_VIS4) or keep it on?

[k.alexander]

Okay I have run a whole year now, and the good news is that the crazy thick ice seems to be fixed. It's still a bit thicker than observations in some regions (especially the Ross Sea) but it at least passes the sanity check! I have more coastal polynya activity, much deeper mixed layers especially right near the coast, and much more winter water. It looks like what was happening before was that spurious supercooling at depth was causing ice to form, but then some of it was almost immediately melting, and that freshwater overstabilised the water column so polynyas couldn't stay open.

However, now that the water column is no longer overstabilised, a giant polynya has formed in the Eastern Weddell Sea. I know this is a common problem and it means I have too much vertical mixing, possibly related to not enough freshwater flux. Some ideas I have to fix it:

- experiment with vertical advection options and vertical mixing parameterisations (maybe just turning off TS_DIF2 would be enough?)
- decrease the atmosphere-ice drag to reduce turbulent mixing
- if I'm really stuck, surface salinity restoring or a flux correction in this region

[k.alexander]

I have run for 3 years now and the sea ice is getting really weird. Here is the concentration (top) and ice thickness (bottom) averaged over the last year:



Notice the giant Weddell Sea polynya opening up, as well as the lack of ice beyond the continental shelf in East Antarctica. In both of those regions there is CDW upwelling to the surface as you can see in this plot of mixed layer depth at the end of the simulation (scale capped at 1000 metres, but the mixed layer depths go all the way to the seafloor):



I compared to Paul Holland's 2014 Journal of Climate paper (which also used CICE with ERA-Interim forcing, just a different ocean model) and the IceSAT data shown there. My sea ice is too thick on the continental shelf, especially in the Ross Sea, and too thin further from the coast.

Here is my total sea ice volume:



Paul and the observations he cites have a seasonal cycle from 3-13 million km^3. So I have about the right annually averaged volume, just too small of a seasonal cycle and getting smaller.

Also, total sea ice area is falling off as the seasonal cycle progresses. It should be 2-17 million km^2 over the seasonal cycle, but I'm getting this as the CDW upwelling starts:



Kate, is this at all similar to what you are seeing in the Arctic with MetROMS? I know CDW upwelling won't be an issue there, but the very thick ice? I'm wondering if the sea ice is just forming so fast that it's too thick to easily blow around. I have been through the dynamics code over and over so I think the issue is most likely thermodynamic.

[kate]

Hi Kaitlin,

I'm currently on holiday in New Zealand, not ideal conditions for sorting this out, but what I noticed was a large difference in the "congel" field between the different CICE runs. Excuse the ncview look, but I should note that the color bar for the fake coupler goes to 0.17 cm/day while the metroms goes to 1.3 cm/day. I think the fake coupler has to be closer to reality.

[k.alexander]

Hi Kate, that's interesting about the congelation. I don't have an alternate coupling setup to compare my runs to, but for what it's worth my annually averaged congelation values are generally around 0.2-0.3, with the exception of right at the coast where they are around 1 (up to 1.4). I *think* that's ok for the Antarctic....We can chat about it in Hobart, but for now enjoy NZ! I'm impressed you found some internet there!

Just for the purposes of organising my brain, and in the hopes that someone with extensive knowledge of advection schemes sees this, I have some more information about my crazy sea ice. It looks like TS_U3HADVECTION/TS_C4VADVECTION is still producing spurious supercooling, although not as much as TS_C4HADVECTION/TS_C4VADVECTION or TS_A4HADVECTION/TS_A4VADVECTION. As a test I ran with TS_MPDATA, which is the only guaranteed monotone advection scheme as far as I understand, for 84 days before it blew up. Ben has also had issues with MPDATA's stability, but I reckoned it would at least provide a benchmark for the first few months.

And the difference was quite dramatic! Here is ice height at the end of 80 days with TS_MPDATA:



and here it is with TS_U3ADVECTION/TS_C4VADVECTION:



Note especially the differences in the Ross and Weddell Seas, where I was getting too thick ice with TS_U3HADVECTION/TS_C4VADVECTION as per my last post.

Here is ice concentration after 80 days with TS_MPDATA:



and with TS_U3HADVECTION/TS_C4VADVECTION:



The source of these differences is the the frazil ice formation. Since this can vary widely over timestep I have averaged it over the 80 days. Here it is with TS_MPDATA:



and here it is with TS_U3HADVECTION/TS_C4VADVECTION:



Lots of weird things going on there, with quite high frazil formation over the continental slope and at the edge of the initial conditions (1m ice everywhere south of 60S). Since the only difference between the two simulations is the advection scheme, the issue must be advective errors causing temperatures to go below freezing and form sea ice where it shouldn't.

I think this can also explain the CDW upwelling I'm getting in the Eastern Weddell and along East Antarctica. Here is frazil ice formation averaged over the second year, once things have settled down a bit, with TS_U3HADVECTION/TS_C4VADVECTION:



Note the super high rates of frazil formation just off the continental shelf break in the Eastern Weddell and in East Antarctica. There's no way this much frazil should be forming there. But I suspect all the resultant brine rejection is triggering the deep convection of CDW that I see in zonal slices of T and S (and in the ultimate destruction of sea ice in these regions).

Something else interesting is that the first few months of ice shelf melt rate are much higher for TS_MPDATA than for the other two schemes. However this might just be that the cavities don't flush out the initial conditions (which are quite warm) as fast with MPDATA, which might also be linked to sea ice formation and the resulting water masses.

While there might be further issues related to congelation (as Kate showed) I think this spurious supercooling caused by advective errors can explain the problems I've been having. The wildly different results between TS_C4HADVECTION/TS_C4VADVECTION, TS_U3HADVECTION/TS_C4VADVECTION, and TS_MPDATA - and the fact that the only realistic looking results (at least for the first 3 months) come from TS_MPDATA, which is the only monotone advection scheme - are very suspicious.

Here are my options as I see them:

1) Figure out why TS_MPDATA is blowing up (it's happening somewhere in the ice shelf cavities but I can't see the details, the restart file is all NaNs) and try to get it to be stable. Has anyone got it to work before? If so, what other mixing options did you use?
2) Maybe there is another tracer advection scheme that is both stable and doesn't suffer from these spurious extrema? Please let me know if you know of anything.
3) I'm sure advective errors are domain-dependent, so maybe there's something I can do to my domain (e.g. more smoothing, higher resolution although I can't really afford that) to fix it.
4) Rewrite the supercooling code (frazil_ice_prod_mod) to something less affected by advective errors. Any ideas? I am curious about Kate's log layer scheme like the Budgell sea ice model uses.

I would love to hear anyone's ideas, either here or in Hobart. I think this is a really interesting issue but there might not be a quick fix for it!

[kate]

There is an HSIMT advection scheme which has many of the properties of MPDATA. The COAWST code has it in their ROMS as an option.

[k.alexander]

Just read a bit about the HSIMT scheme and it sounds great (I particularly like the bit about being more efficient than MPDATA!) Who knows if it will blow up like MPDATA does, but it's worth a go. What's the best way to get access to the COAWST code so I can merge the option into MetROMS?

[kate]

Simply email John Warner at USGS - a quick search in here should turn up that info.

[kate]

Looking through more fields, the surface temperature in April in metroms is rather colder than in the fake coupler. I also saved monthly means of all the surface fluxes and the incoming longwave is uniformly 180 W/m2 in the metroms solution with zero incoming shortwave. By April, I have more incoming longwave in the fake coupler plus of course non-zero shortwave. Can you check those fields? The sensible and latent heat fluxes look OK, I think.

[k.alexander]

Hi Kate, I'm not getting that problem - in spring (October) shortwave is nonzero, longwave is non-constant, and both look reasonable. I assume it is CICE radiation fields that you mean? Maybe you have the wrong CICE radiation option turned on?

[kate]

In my code, ROMS passes the atmospheric fields to CICE - it has to read them anyway, so that seems most efficient. Since coupling.dat mentions those fields, I was expecting metroms to behave the same way. I guess I was wrong.

[k.alexander]

Ah that makes sense. Yes, ROMS and CICE read the forcing fields separately. A bit inefficient but since they aren't coupling every timestep I think it's the only way to do it properly (at least without rewriting a lot of code!)

FYI CICE wants specific humidity, whereas ROMS accepts either specific or relative humidity - that tripped me up for a while.

[Paul_Budgell]

With regard to MPDATA you were discussing a couple of days ago, unless the code has changed it won't work properly for subzero temperatures. It works at keeping things "positive". I have needed to, say, add 5 degrees to temperature, do the advection with MPDATA, then subtract 5 degrees. I don't know about HSIMP.

[k.alexander]

Ah ok interesting about MPDATA. I'm not sure how HSIMT handles negative temperatures, Kate do you know?

The MetROMS developers in Norway suggest TS_A4HADVECTION + TS_C4VADVECTION with quite a lot of diffusion (TS_DIF2) to smooth out the errors. I'm going to try that next. Feels a bit dirty to use diffusion errors to get rid of advection errors, but whatever works...

[kate]

No, I don't know about HSIMT and negative temperatures.

[stlaur]

Just a suggestion---sensitivity to the choice of the advection scheme may possibly indicate that your grid stretching is more than it should be.

Could be worth trying:
(1) decreasing the number of vertical levels,
(2) use conservative stretching parameters,
(3) increase the horizontal resolution.

Cheers,
Pierre

[mitya]

Kaitlin, did you switch to github version of metroms? I noticed that you did some commits to source.uit.no version of it after Nils migrated it to github metno repositories.

Dmitry

[kate]

Kaitlin - great talk of the many ways to get into trouble in this modeling business. I was wondering, have you taken a really close look at your periodic join? If your ROMS grid has Lm, Mm of say 1024, 480, then the CICE grid should be dimensioned the same, except METROMS seems to want it to be 1026 by 482. Not a big deal except for periodic grids, where you will end up with two extra points at the join.

[k.alexander]

Just got access to HSIMT, and confirmation that it is okay with negative temperatures. I will try that today or tomorrow. However, I don't think my issues are due to advection anymore. That would show up in the frazil fields which are very low in the Ross and Weddell. Instead the ice there is getting very thick due to congelation. I keep going back and forth as to whether I think it's a thermodynamics issue (ice forming too fast) or a dynamics issue (ice not being exported fast enough). So HSIMT will hopefully clear up the last few advection errors I'm getting (I see a few in the Amundsen Sea, and out the near the northernmost extent of the sea ice) but I don't expect it to drastically change my solution.

With regards to GitLab, I have my own branch including ice shelf thermodynamics, but I am keeping an eye on the GitHub updates to make sure I don't miss anything important.

And with regards to the periodic boundary/domain sizes, that caused some problems earlier (off by one errors in the coupling code wherever there is a periodic boundary) but they are all fixed now. I have the same number of points in my ROMS rho-grid as in the CICE grid. My CICE input grid file is for the CICE u-grid, which corresponds to the ROMS psi-grid, but you have to extrapolate one extra row and column to make it the right size. The CICE code calculates the t-grid at runtime which is the same as the ROMS rho-grid. I'm happy to share my code if you're interested, Kate.

[kate]

I heard Hernan say that HSIMT is missing something in fully 3-D applications. Did you hear that too? What about Gary's advection scheme limiter?

Could you share your code with Rajesh? He's working on this stuff while I'm funded to do other things right now. I'm also intrigued by Hernan's NUOPC path if that's as imminent as he says.
http://journals.ametsoc.org/doi/10.1175 ... 14-00164.1

[k.alexander]

Yes I missed exactly what Hernan said about HSIMT, I think it was something about topography? However the COAWST version is apparently a bit different (or so Hui Wu tells me) so maybe that is fixed. I'm just testing it for now, if I end up keeping it I will do more investigation. I am also going to test Gary Brassington's TS_U3HADVECTION limiters (lots of people were excited about that weren't they!), he is just packaging up the code for me.

I'll email the python grid scripts for MetROMS to Rajesh and cc you.

[hwusklec]

Hi this is Hui Wu (HSIMT). Kaitlin, it seems you have some new results with HSIMT, right?

I am not familiar with ROMS code, Tarandeep in John's group helped to rewrite the HSIMT code for COWAST. It is very hard to evaluate an advection schemes in realistic simulations, since there are too many unknowns in model, such as turbulent mixing. But MPDATA doesn't work for zero temperature is expected. More precisely speaking, MPDATA works only when tracer value stay positive or negative. In the algorithm of MPDATA there is a step with 0.5(a(i-1)+a(i)) in the denominator. Hence if a(i-1) and a(i) have the opposite values, which is very like around the zero point, MPDATA will crash. I think the sea ice simulation is a wonderful testcase for advection schemes. If the temperature dropped below -2 the sea water will be converted to ice, hence the overshooting of seawater temperature (if below -2) will be "recorded".

Sigma-like coordinate,including S coordinate, often brings additional challenges for advection schemes. High order schemes such as HSIMT (or others) uses at least three grid points to calculate the cell-interface tracer flux. If there is a very steep topography, say a sea mount, the three neighbouring grid points on one sigma coordinate are actually at very different vertical location, hence very likely they have distinct temperature. Using them to calculate the advection will make a big mess. This is why sometimes HSIMT works not very good around rough topography. Low-order schemes, on the contrary, may work better in this case, since it uses only two grid points, thus the "sigma titling" effect is not that severe. A good practice is to design the stretching coordinate carefully so that the titling S layers (due to topography) are well below the thermocline.

[k.alexander]

Thanks Hui for your message. Yes, as we discussed over email, I was getting more advection errors with A4 + strong diffusion than I realised. Switching to HSIMT removed the weird "corner" of sea ice in the Ross Sea, and there is thinner sea ice in the outer Weddell Sea now too. I am currently running standardised experiments with the different advection schemes to try to quantify this difference. But I think I will stick with HSIMT.

Thanks for the explanation of HSMIT around rough topography. If there were errors occurring due to this phenomenon, what sort of warning signs should I look for? I assume there would be no spurious supercooling since HSIMT is monotone.

Since there is (correctly) less sea ice with HSIMT, the Cosmonaut Polynya keeps swallowing up everything after a few years, and even wrapping around East Antarctica. I think there is a dynamics issue preventing sea ice from exporting away from the coast, so it doesn't melt far enough away from the coast and this missing freshwater is leading to destabilisation of the water column, upwelling of CDW, the usual story. I'm working on it...but fixing the advection errors is a good first step.

[hwusklec]

It is not easy to say what kind of signal will occur around rough topography. If your S-coordinate intersects with pycnocline with big angle and if a strong cross-isobathic flow is present, artificial "diapycnal mixing" may occur around the slope. A general impression of this numerical phenomenon is "diffusive". It seems to be a problem of advection scheme, but in fact it is a problem of vertical coordinate. But it should be OK if your horizontal resolution is high and the topography is smoothed (which is often required in ROMS).