Fix land mask around Antarctic coast for 0.25deg topography

[micaeljtoliveira]

This issue was first discussed here. More details in this comment and comments following.

The proposed solution is:

We decided that the way forward is to change the land mask for the Antarctic coast (and nowhere else) to match GEBCO 2014, which will convert some ocean regions to land, but also convert some land to ocean, which will change the number of active tiles and may require tweaks to remove non-advective cells.

[micaeljtoliveira]

I've used the topog2mask.py script in topogtools to create a mask from the GEBCO 2014 topography. Here is the old (left) versus new masks (right):

The differences look reasonable (red means new land, blue means new water):

Considering that the grid extends beyond 80º S, I would say that there's no need to extend it any further, as they don't seem to be any changes to the mask close to the grid boundary.

@aekiss Does this look reasonable to you?

[aekiss]

Great that the grid doesn't need extending.

Currently the minimum x-size of the grid is 6km, occurring at the Antarctic margin. Can you tell from the grid data how much smaller dx will be if we move the coastline further south? We don't want it too small, as it may reduce the maximum timestep we can use while avoiding CFL instability (particularly "remap transport: bad departure points" errors in CICE).

Do you know if extending the ocean further south will require 1 more row of tiles, or will it need more?

The top 2 plots are partly obscured by the filled coastline data, making it hard to see the changes. Could they be replotted with only the model data displayed?

It looks like this removes some land at the tip of the Antarctic Peninsula and adds and removes some islands to the west of the Antarctic Peninsula? This would be clearer if the regions that are land in both old and new masks were plotted in a different colour.

It looks like this also fixes a glitch in the coastline at 0° longitude?

[ofa001]

Yes @aekiss and @micaeljtoliveira a little worried at what is occurring with the blue around the tip of peninsula are a lot of islands being removed. I guess I said for the large ice shelves (Ross Ronnie/Filchner they probably had originally very straight edges when in reality they are more complex like in the new topography.

[micaeljtoliveira]

Currently the minimum x-size of the grid is 6km, occurring at the Antarctic margin. Can you tell from the grid data how much smaller dx will be if we move the coastline further south? We don't want it too small, as it may reduce the maximum timestep we can use while avoiding CFL instability (particularly "remap transport: bad departure points" errors in CICE).

I'll calculate this and let you know.

Do you know if extending the ocean further south will require 1 more row of tiles, or will it need more?

No yet. I first need to combined the new mask with the old one, as we only want the changes from the new mask in the Antarctic region.

The top 2 plots are partly obscured by the filled coastline data, making it hard to see the changes. Could they be replotted with only the model data displayed?

Here they are:

It looks like this removes some land at the tip of the Antarctic Peninsula and adds and removes some islands to the west of the Antarctic Peninsula? This would be clearer if the regions that are land in both old and new masks were plotted in a different colour.

Yes, it seems there are quite some changes around the Antarctic Peninsula. Here is the new plot:

It looks like this also fixes a glitch in the coastline at 0° longitude?

The glitch seems to be coming from the old mask. I would say the new mask looks better in that area.

[aekiss]

Thanks, those plots are great.
Might be worth (visually) comparing to the 0.1deg land mask to see what choices @russfiedler made, particularly at the end of the peninsula.

[micaeljtoliveira]

@aekiss The new value for dx is approx. 5.5 km, corresponding to a new southernmost latitude of -78.65.

I've also plotted the 0.1deg land mask (right figure bellow). I'm not sure what to make out of it.

[aekiss]

Hmm, the old 0.25 mask looks more like the 0.1 mask at the end of the peninsula than the new 0.25 mask does. ie the old mask looks "better" there.

I'm guessing there's a mix of deep water and low-altitude ice shelves/land on the GEBCO grid, which averages to something below sea level when regridded to the 0.25 grid. In this situation it may be more valid to treat a 0.25 grid cell as land if it contains any land on the GEBCO grid.

Probably the simplest thing to do would be retain the old 0.25 mask everywhere north of 65S.

[adele-morrison]

Keeping the old mask north of 65S seems a bit arbitrary.

I think testing something more reproducible would be better (e.g. the suggestion 'to treat a 0.25 grid cell as land if it contains any land on the GEBCO grid').

[aekiss]

That could be done by basing the land mask on the frac (rather than depth) field in topg_new.nc (as output by gen_topo.f90), i.e. setting all points with frac < 1 to land.

[micaeljtoliveira]

That could be done by basing the land mask on the frac (rather than depth) field in topg_new.nc (as output by gen_topo.f90), i.e. setting all points with frac < 1 to land.

That should be straightforward to implement. Let me try it and see how it looks.

[aekiss]

[aekiss]

maybe too much land now

[aekiss]

here's frac < 0.8

[aekiss]

it starts to get a bit arbitrary

[micaeljtoliveira]

Yes, that's a bit arbitrary. Plus, it's not clear how it will impact the mask in other areas.

[aekiss]

Here's frac <0.5. I think this looks closest to the 0.1deg mask. There will be quite a few non-advective points to remove.

[aekiss]

here's the whole coastline with frac<0.5 (click to enlarge)

[aekiss]

setting frac<0.5 seems somehow less arbitrary than any other value, and seems to do a good enough job - it's be interesting to see what an updated difference plot would look like with this mask.

[adele-morrison]

Great, let's use 0.5!

[aekiss]

Actually, to be certain we won't have any ocean points with depth above sea level, we need to set land to (frac < 0.5) AND (depth <= 0.0) in https://github.com/COSIMA/topogtools/blob/master/topog2mask.py#L31

[micaeljtoliveira]

@aekiss Shouldn't the condition be (frac < 0.5) OR (depth <= 0.0)?

If I understood the issue correctly, the problem is that some points have depth > 0 and frac < 0.5, but are being marked as ocean, while we would like them to be marked as land. To do this, we must mark as land any points that have depth <= 0 OR frac < 0.5, no?

Anyway, the condition (frac < 0.5) AND (depth <= 0.0) leaves the new mask unchanged, while the OR condition yields the following new mask (left), which looks quite similar to the one from the 0.1deg topography (right):

And here is the plot showing the changes wrt the old mask:

I would say we are on the right track!

[adele-morrison]

Looks great!

[aekiss]

Oops, yes you're absolutely right re. using OR rather than AND. The new mask looks excellent.

Next step is to fix non-advective points in the mask.
They can be found with non-advective.ipynb and fixed with editTopo.py.
If I recall correctly it might also be possible to do some of this automatically with fix_nonadvective_mosaic.f90. It might not pick up all of them but could clean up most of them more easily.

[micaeljtoliveira]

Looks like there are quite a lot of new nonadvective cells (>200!). I tried using fix_nonadvective_mosaic.f90, but it doesn't seem to work properly, as the python notebook still identifies the same nonadvective before and after "fixing" those cells.

[micaeljtoliveira]

fix_nonadvective_mosaic.f90 also seems to have a strange side effect when used on top of the new mask. After "fixing" the advective cells with the tool, some ocean cells now have a depth of zero, which should definitely not happen (took me most of the day to realize it was the tool that was doing this, not my new mask... 😞

[russfiedler]

@micaeljtoliveira I have vague memories of a (logic) bug in old versions of this code. I'll have a hunt around and see where/when it was fixed. Can you point me to the location of the topogrhy files in question or put them in a public place and I'll see what is going on?

[aekiss]

@russfiedler is fix_nonadvective_mosaic.f90 designed to work on mask files, or only topography?

[aekiss]

@micaeljtoliveira non-advective.ipynb is overzealous, so you probably won't need to remove all the non-advective points it identifies. Only the points with no advective connection to the ocean need to be removed (or widened).

[ofa001]

Thanks @MartinDix Balleny Islands I have heard of Islands but not sure I know of Peter Is. Glad we know have matching grid set ups and can do a test with the new ocean topography.

[access-hive-bot]

This issue has been mentioned on ACCESS Hive Community Forum. There might be relevant details there:

https://forum.access-hive.org.au/t/scope-of-access-nri-access-om2-release/172/3

[access-hive-bot]

This issue has been mentioned on ACCESS Hive Community Forum. There might be relevant details there:

https://forum.access-hive.org.au/t/bathymetry-for-ocean-model-at-any-resolution/462/8

[aekiss]

@MartinDix, @ofa001, @DaveBi - just wondering if you've done a test run with this corrected 0.25deg bathymetry? I'd like to make this the default for ACCESS-OM2-025 but I want to know if it's suitable for ACCESS-CM2-0.25 first.

[MartinDix]

Just started a new run with this bathymetry.

[aekiss]

Great! I'm interested to see how that goes.

[MartinDix]

Some global means statistics from run in progress, https://accessdev.nci.org.au/p66/mrd599/access-cm2/control025/

Original run had a crash around Antarctica about every 20 years on average. No problems in the new run yet.

[aekiss]

Thanks @martin, glad to hear it's more stable now.

Those plots mostly look good, but these things jumped out for me:

The Baltic is still heading very low (as with the previous topo) - was the initial condition based on obs (eg WOA)?

cf. ACCESS-OM2-025 test runs

Min SSS is taking a weird plunge - hope it levels off!

Sept SH SIA has dropped slightly, might be a bit below obs now; SIV also has the same dip

[MartinDix]

The mask has some isolated ocean points near the Antarctic Peninsula and one of these is where the salinity goes to zero in year 24.

Image above at #265 (comment) shows these as connected but they're isolated in /g/data/ik11/inputs/access-om2/input_20220919_025deg_topog/mom_025deg/ocean_mask.nc.

Did something happen when removing the non-advective points?

[ofa001]

That area you have circled in the West Antarctic @MartinDix looks concerning, its probably topography under an ice shelf rather the coastline itself. We proably need to eliminate it. The runs not going to be happy if salinity reaches zero.

[micaeljtoliveira]

Did something happen when removing the non-advective points?

I think I didn't realize that the fixing of non-advective cells could create this kind of isolated points. In the workflow I used, isolated cells are only removed before fixing the non-advective cells.

[ofa001]

Hi @micaeljtoliveira, that means it might need a second pass to remove any new isolated point.

[micaeljtoliveira]

I'll try to get an updated topography and masks ASAP.

[aekiss]

thanks @micaeljtoliveira

[micaeljtoliveira]

@MartinDix I've finally updated the inputs with a new version of the topography that should have those isolated points removed. The new inputs are in /g/data/ik11/inputs/access-om2/input_20230515_025deg_topog. Would be great to have this tested at some point.

(see the following PR for a plot of the new topography vs old: COSIMA/make_025deg_topo#2)

[access-hive-bot]

This issue has been mentioned on ACCESS Hive Community Forum. There might be relevant details there:

https://forum.access-hive.org.au/t/access-cm2-025-storage-request/1649/1

[aekiss]

@wghuneke has a test run been done with this new topography and land mask? If it works well it would be good to make this the default in the standard configs.

[wghuneke]

Yes, the new bathymetry works well for me, no crashes anymore. I'd support making this the new default product.

[ofa001]

Thats good news Willma, we can start using the updated bathymetry for new 0.25 runs at CSIRO

[micaeljtoliveira]

I think this can be closed now 🎉

[aekiss]

Closed - 0.25° JRA55 configurations (RYF and IAF, master and master+bgc) now all use /g/data/ik11/inputs/access-om2/input_20230515_025deg_topog.