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... :disappointed:

[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).

[russfiedler]

@aekiss Only on topography. It infers the mask from the topography file (topog<=0.0). Number of levels is calculated via the vgrid.nc file. Make sure that all double precision numbers are exactly representable as singles or there could be unforseen results. There were some bad 75 level vgrid.nc files lurking about at various stages.

[aekiss]

Thanks @russfiedler and apologies @micaeljtoliveira for giving you a bad suggestion.

[micaeljtoliveira]

@russfiedler My topography files are here: /g/data/v45/mo1833/topography/make_025deg_topo. This folder is a clone of https://github.com/COSIMA/make_025deg_topo with a single modification: I'm calling a script to modify the land mask before it is applied to the topography. The topography after applying the mask is the topog_new_edited_deseas_partialcell_mindepth_masked.nc file and topog_new_edited_deseas_partialcell_mindepth_masked_fixnonadvective.nc is the topography after calling fix_nonadvective_mosaic.f90.

I didn't check this in detail, but after thinking a bit more about what is going on, my guess is that the tool is creating some new land to fix the nonadvective cells.

@aekiss No worries, I'm still learning quite a lot with all these issues ;)

[russfiedler]

@micaeljtoliveira Yes, your suspicions are correct. There's definitely land being created to eliminate non advective cells. The regions I looked at seemed ok to me.

[micaeljtoliveira]

@russfiedler Thanks for looking into this. After discussing further with @aekiss we decided that the best would be to fix most of these non-advective cells at the level of the mask. I now have a script to do this and will fix a couple of remaining cases by hand.

While looking into this, I think I found a bug in the current version of fix_nonadvective_mosaic.f90. I've created a pull request (https://github.com/COSIMA/topogtools/pull/5). Would be great if you could have a look!

[micaeljtoliveira]

Reopening this, as it needs changes to https://github.com/COSIMA/025deg_jra55_ryf and https://github.com/COSIMA/025deg_jra55_iaf, plus some further testing.

[micaeljtoliveira]

Some updates regarding this issue. After successfully generating the updated topography and fixing the non-advective cells, I've regenerated the regriding weights used to couple the atmosphere with the ocean/ice models and updated the mask table. All the new files can be found here: /g/data/ik11/inputs/access-om2/input_20220919_025deg_topog.

There is now a new branch named update_antarctic_topo in both https://github.com/COSIMA/025deg_jra55_ryf and https://github.com/COSIMA/025deg_jra55_iaf using the updated inputs. I've done a quick check by running the 025deg_jra55_iaf configuration for one month with daily output. As expected, when comparing with a similar run with the old topography, changes in the outputs are only observed in the antarctic region, and they become larger with time. All the differences look reasonable to me, but it would be good if people could test the new files and report any feedback here.

[aekiss]

@ofa001, @davebi, @martindix we'd be interested in hearing whether this new topography /g/data/ik11/inputs/access-om2/input_20220919_025deg_topog works better for ACCESS-CM2-025 than the previous one (TOPO5 /g/data/ik11/inputs/access-om2/input_20211125_025deg_topog/mom_025deg/topog.nc; see https://github.com/COSIMA/access-om2/issues/158)

[aekiss]

Note that there are changes to the land mask and processor layout.

[ofa001]

Hi @aekiss noting that @MartinDix is starting a new 0.25 coupled run to try out new GM values this might be a good time to test the topography as well though he was going to start part way through the existing run, is that correct @MartinDiX @AndyHoggANU

[MartinDix]

The UM uses fractional land and a land mask that includes all points where the land fraction is non-zero. Process is to average the MOM land mask to the UM grid using the ESMF weights and set any grid boxes with land fraction < 0.01 to all ocean. This was an arbitrary limit and the Met Office don't use a limit at all so have some tiny land fractions. This map shows change in land fraction with the new MOM mask, red more land, blue less land.

The change in land sea mask shows some extra islands that I'd never heard of which have land fraction just > 0.01 on the N96 grid. Peter I island at approx 70S 90W and the Balleny Islands at 66S, 164E. Google maps doesn't even know the name of these! https://www.google.com/maps/@-66.9249036,163.6697566,7.37z. Bing maps does name them.

[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 https://github.com/COSIMA/access-om2/issues/265#issuecomment-1231159693 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.