After meeting with Paul and Keith at BAS I tried again to work out what what was going on with the funny boundary layers for my crevasse setups.

Running

python geostrophic_balance.py

diverges after 1 day. This script is based on simplifying the crevasse setups to remove melting and crevasses.

I was expecting it to reach a steady geostrophic balance, where the body force in the y direction (equivalent to a fixed free surface tilt) is balanced by negative u flow.

The velocity field at 1 day is very noisy at the boundary, especially at x = 0

The log is here
02.06.24_nobcs_geostrophicspinup.txt

I also realised that if you replace the velocity pressure bcs with

vp_bcs = {1:{}, 2:{}, 3:{}, 4:{}, 5: {'un': no_normal_flow},
          6: {'un': no_normal_flow}}

i.e. add empty dictionary to each boundary the simulation runs for the alloted ten days and you actually get completely different flow field!

This is almost identical to the same pattern as I saw in all my crevasse runs previously. The reason there is a difference is that this line

I found this out by doing quite a bit of testing (for more details see here
Geosptrophic box.docx)
and ultimately turning off the momentum advection term proved that you could get geostrophic balance with open bcs. I then realised we are missing a surface integral after IBP in the momentum advection term and adding this in fixes the problem! But I will describe that more in the PR to fix this...!