Bug at tripole seam

[aekiss]

There's a bug in CICE or MOM at the tripole seam that crosses the north pole along the 100W/80E meridian north of 65N. Apologies for the long post, I wanted to get it clear in my head.

From the grid you can see that on the Atlantic side of the seam positive u velocity runs from Canada to Siberia; on the Pacific side it is the reverse orientation.
v velocity also changes sign at the seam. On both sides of the seam positive v velocity has a northward component. There is also a zonal component of v velocity that changes sign at the seam (positive v is partly eastward in western Siberia/Canada and partly westward in eastern Siberia/Canada).

Here are some views off Siberia. Date is in the heading.

Starting with u velocity (in pairs, first ice, then ocean).

Most of the time the u velocity has a sign change at the seam as it should (positive u is down on the left of 80W and up on the right in this view):

But sometimes the sign of u does not change across the seam i.e. the physical tangential velocity abruptly changes direction across the seam.

It starts small


and grows

This strong shear at the seam leads to flow instability within a few days

The same sequence of days in v velocity shows no such problems (positive v is rightward on the left of 80W and leftward on the right of 80W in this view)

The wind stress also shows no problems at the seam.

Movies of the ice shear show this problem comes and goes in various locations along this seam.

Clearly there's something wrong with at least one term in the u tendency at the seam in MOM or CICE. The fact that the problem is not apparent most of the time suggests the usually-dominant terms are correct, and the problem appears only on occasions when the buggy term is unusually large.

I expect it's a term with a stencil that crosses the seam but doesn't take into account the sign change. The fact that it is asymmetrical (the positive u on the west spreads eastwards) and that the v velocity is eastward suggests the problem is with the y advection of x momentum in the u tendency in MOM or CICE, ie vdu/dx (remembering that x,u are along-seam and y,v across-seam).

It is unclear whether the problem arises first in the ocean or ice (it seems pretty simultaneous in the figures) but I'm trying runs with zero ocean-ice drag to try to identify the culprit. The fact that v is more strongly eastward in the ice than the ocean when the problem starts on 0017-08-03 makes me suspect the bug is in CICE rather than MOM.

[aekiss]

We are using incremental remapping (advection = 'remap') in CICE - see https://cice-consortium.github.io/CICE/cice_2_science_guide.html#horizontal-transport

[aekiss]

...but as far as I can see that only applies to tracer transport and there's no ice momentum advection in CICE: https://cice-consortium.github.io/CICE/cice_2_science_guide.html#dynamics

• so does that point the finger at MOM?

[nichannah]

It would be interesting to look at the ice-ocean coupling fields. One way to do this is to change the the words 'EXPORTED' and 'IGNORED' (these are equivalent) to EXPOUT in the namcouple file. Best to choose only the necessary fields because it writes out a v. large amount of data.

[aekiss]

In MOM we're using velocity_advect_centered=.true., ie advection_scheme = VEL_ADVECT_2ND_ORDER.
This is done by subroutine horz_advection_centered in ocean_velocity_advect.F90.
I'm probably missing something, but I can't immediately see how this handles the sign change of u across the seam: https://github.com/mom-ocean/MOM5/blob/master/src/mom5/ocean_core/ocean_velocity_advect.F90#L368

[aekiss]

thanks @nicjhan that's a good idea

[ofa001]

no you only advect tracers not velocities. One issue is that I thought of is that the u v position in cice are not the same as in MOM this led to some noise in ocean stress at coast in strong currents It could be this Apologies I am using phone not laptop to send this _______ From: Andrew Kiss [notifications@github.com] Sent: Wednesday, 28 March 2018 5:13 PM To: OceansAus/access-om2 Cc: Subscribed Subject: Re: [OceansAus/access-om2] Bug at tripole seam (#86) ...but as far as I can see that only applies to tracer transport and there's no ice momentum advection in CICE: https://cice-consortium.github.io/CICE/cice_2_science_guide.html#dynamics * so does that point the finger at MOM? — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub<#86 (comment)>, or mute the thread<https://github.com/notifications/unsubscribe-auth/ANrRggJMCqTjeO0XXeKOPqCqx-1WklmKks5tiyn8gaJpZM4S-Csh>.

[ofa001]

You could look just at the a few time steps Dave has aflag to do this in the coupled set up I will check tomorrow if its in the stand alone I would focus on the ice icean stress based on my earlier comment on uv mismatch

…

________________________________________ From: Nic Hannah [notifications@github.com] Sent: Wednesday, 28 March 2018 5:20 PM To: OceansAus/access-om2 Cc: Subscribed Subject: Re: [OceansAus/access-om2] Bug at tripole seam (#86) It would be interesting to look at the ice-ocean coupling fields. One way to do this is to change the the words 'EXPORTED' and 'IGNORED' (these are equivalent) to EXPOUT in the namcouple file. Best to choose only the necessary fields because it writes out a v. large amount of data. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub<#86 (comment)>, or mute the thread<https://github.com/notifications/unsubscribe-auth/ANrRgpYrm3Dr0j8-RKybrovKOzuaCeSdks5tiyuXgaJpZM4S-Csh>.

[ofa001]

You are on the right tracki in MOM in CICE its done differently so the mismatch may be the issue

…

________________________________________ From: Andrew Kiss [notifications@github.com] Sent: Wednesday, 28 March 2018 5:32 PM To: OceansAus/access-om2 Cc: Subscribed Subject: Re: [OceansAus/access-om2] Bug at tripole seam (#86) In MOM we're using velocity_advect_centered=.true., ie advection_scheme = VEL_ADVECT_2ND_ORDER. This is done by subroutine horz_advection_centered in ocean_velocity_advect.F90. I'm probably missing something, but I can't immediately see how this handles the sign change of u across the seam: https://github.com/mom-ocean/MOM5/blob/master/src/mom5/ocean_core/ocean_velocity_advect.F90#L385 — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub<#86 (comment)>, or mute the thread<https://github.com/notifications/unsubscribe-auth/ANrRgjv2QhqdHLQ-Pjt0js99u95aAasQks5tiy6NgaJpZM4S-Csh>.

[aidanheerdegen]

Forensic work! Nice pics as well ...

[StephenGriffies]

Do these issues arise in the MOM-SIS simulations? If not, then I suspect CICE.

[aidanheerdegen]

Also might be worth checking if the same thing happens in the 1 deg simulations. They're a lot coarser resolution, so maybe not? If it does, it would be a more tractable to debug with a cheaper simulation. Particularly if you're dumping coupling fields.

[nichannah]

Just some thoughts about the coupling that may be relevant here:

The coupling exchange between ice and ocean is based on grid cell indices not on lat, lon coordinates. This is because they should have identical grids.

Given the above the coupling could introduce some weirdness if the models have a different idea of where the seam is in a grid-cell sense, i.e. a given u-point does not have the same i, j index in the grid definition.

[aekiss]

The sea ice cover also disappears on the eastern side of 80E

[russfiedler]

@aekiss The flipping of velocities across the tripolar seam is done in mpp_update domains.

[aekiss]

Thanks @ofa001 but I think this might be a different problem. This animation shows the spurious shear at the seam comes and goes, so it seems unlikely to be a grid mismatch. (See the yellow line; ignore the dotted white grid line.)

[ofa001]

I will check the daily data from the CMIP 5 data I did just archive them off the system but can recover them. we did see an occasional polynya starting at the tripole in summer when there was no convective generation from the ocean to cause it. and also some wave noise, on one occasion but I didn't forensically examine those files as I was seriously ill the year they were generated 2011-2012, and was more interested in the Antarctic daily files. You may be right it could be in the advective term, but there was a grid issue at the coast again intermittently, during strong currents.
Definitely worth looking at the MOM-SIS runs, Love the shear zones in the movie.

[aekiss]

Here is a repeat of the ice, ocean u velocity pairs for the same days as above but for a run with ice-ocean drag coefficient dragio reduced by an order of magnitude (from 0.00536 to 0.0005), and timestep reduced from 240s to 24s and barotropic_split reduced from 80 to 8.

The ice is moving faster due to reduced drag. It also seems clear that u velocity of the same sign at the seam first appears in the ocean.

[aekiss]

These figures make my last post clearer. The first two are ice, the second two are ocean. The first of each pair is the original case, the second of each pair is with reduced dragio. You can see that reduced drag decreases the problem in ice and increases it in the ocean, suggesting the ocean is the source.

[aekiss]

Mind you, there's still a strong discontinuity at the seam in the ice (u on the west of the seam has much larger magnitude than on the east)...

[aekiss]

It doesn't look like MOM's horizontal velocity advection is to blame.

This is u vel from a run starting from the same nonlinear initial condition on 1 Aug as my first post but with zero_velocity_advect_horz=.true. so there's no horizontal velocity advection in MOM (the eddying structures are decaying fossils from the initial condition). All other parameters are the same as my first post. We still get u of the same sign on either side of the seam by day 4.

I'm wondering now whether there could be a bug in the ocean-ice stress coupling. I'll output strocnx, strocny etc in a new run to check.

[russfiedler]

It looks to me as if the problem is in the coupling of the stresses to MOM from CICE..

In cpl_forcing_handler.f90

tiostrsu(:,:,:) = strairx_ocn(:,:,:) * (1. - aice(:,:,:)) - strocnxT(:,:,:) * aice(:,:,:)
tiostrsv(:,:,:) = strairy_ocn(:,:,:) * (1. - aice(:,:,:)) - strocnyT(:,:,:) * aice(:,:,:)

But strairx_ocn and stroncxT are both on the T grid centre whereas MOM wants them on the NE corners. This means that the assumed symmetry in the stresses is broken and you can set up a shear.

There needs to be a remapping to the U grid so that everything matches along the fold. In MOM-SIS. this is hidden away inside SIS.

We need calls of t2ugrid_vector( tiostrs[uv]) somewhere before the coupling to the ocean.Alternatively do the rgeridding in MOM

[russfiedler]

Having thought about this for a while I believe that we should try to be consistent with MOM-SIS and only interpolate the air-sea stress to the U grid and use the ice-ocean stress strocn[xy] as directly calculated on the U grid in CICE. This means using aiceu for the weightings.

[russfiedler]

Ok. iostrs[u]v] are regridded in CICE_RunMod.F90. Need to look elsewhere.

[aekiss]

I've confirmed that this bug still exists in 01deg_jra55v13_ryf9091, despite using more recent executables, including a lot of updates to cice:

/short/public/access-om2/bin/yatm_b6caeab.exe 
/short/public/access-om2/bin/fms_ACCESS-OM_da2a93f_libaccessom2_b6caeab.x
/short/public/access-om2/bin/cice_auscom_3600x2700_722p_47650cc_libaccessom2_b6caeab.exe

(MOM commit da2a93f is on this branch)

Also, the script for the figures above is here.

[russfiedler]

I wonder if this behaviour is related to the inconsistent background viscosities that I identified occuring along the tripolar fold here mom-ocean/MOM5#282 Essentially you're solving two different equations at the same point. If you plot the friction terms along j=2700 you can see they don't align.

[aekiss]

I was wondering the same thing, but ncar_boundary_scaling=.false. for these experiments, so they shouldn't be affected by that bug, right?

[aekiss]

Just noting that the CICE grid misalignment issue #190 does not explain the shear on the seam, as the misalignment was fixed in /g/data/hh5/tmp/cosima/access-om2-01/01deg_jra55v13_ryf9091 but the problem still remains, as noted above.
For example, /g/data/v45/aek156/figures/access-om2-01/Ice_Crash_ACCESS-OM2-01/01deg_jra55v13_ryf9091_Arctic_closeup_shear_1934-08.png

[aekiss]

Sec 4.2.2 of the CICE5 manual says

• CICE5 supports having the tripole fold (seam) at either U points or T points. (We use U points, i.e. ns_boundary_type=‘tripole’, to match MOM).
• In the U-fold tripole grid, the poles have U-index nx_global/2 and nx_global on the top U-row of the physical grid, and points with U-index i and nx_global − i are coincident.
• within rows along the fold, coincident points must always have the same value. This is achieved by averaging them in pairs.

Taken literally, "averaging them in pairs" seems incorrect for velocity components because they have opposite signs in the coincident U-points. Perhaps there's a term in CICE where the velocity sign change isn't properly taken into account in this averaging?

[apcraig]

I'm looking at the tripole halo update in CICE at the moment, basically writing a unit tester for it, so am trying to understand the implementation. In CICE, you can set field_type_scalar or field_type_vector for a halo update. The vector implementation will switch signs across the seam, the scalar will not. Is it possible that a haloupdate for a vector field is done with "scalar" when it should be done with "vector"? I'm working with the latest version of CICE, https://github.com/CICE-Consortium/CICE, but if this is CICE5, maybe there's a bug there (or maybe there's a bug in the current version too).

I found this issue trying to get more info on the tripoleT implementation. I'm fairly confident the tripole halo update is working correctly in CICE6, but maybe the halo update call is incorrect.

[aekiss]

Thanks for this suggestion @apcraig and sorry for my very slow reply. This was a very helpful hint - it led me to finding an unrelated bug.

I'm not sure I've found anything to explain this current issue, but it was unexpected to me that the stress tensor components are updated as scalars - I guess that just reflects my lack of understanding?

[aekiss]

I think the code is correct in transforming the stress tensor as a scalar.

The stress tensor $\sigma_{ij}$ is contravariant. Coordinate transformations change the tensor components according to $\sigma_{ij}'=a_{im}a_{jn}\sigma_{mn}$, where $a_{ij}$ are the components of a rotation matrix and the Einstein summation convention is used.

The coordinate transformation across the tripole seam is a 180° rotation, so

$$a_{ij}=\left[{\begin{matrix} -1 &0 \\\ 0 & -1 \\ \end{matrix}}\right]$$

Therefore the coordinate transformation leaves the tensor components unchanged, i.e. $\sigma_{ij}'=\sigma_{ij}$, which is the same as what happens for a scalar.

This seems to make intuitive sense, since the 180° coordinate rotation changes the signs of the components in both surface normal vectors and traction vectors.