AUS2200 instability analysis

Hi All (specifically @Scott and @MartinDix)

Following on from our chat in the AUS2200 meeting the other day, I’ve made this fun (though scientifically useless) visualisation tracking the maximum vertical wind across the AUS2200 domain for both the 30-second timestep and 75-second timestep model runs here: multiplot.mp4. Apologies in advance for my lack of design sense. Anyway, the green dot/line represents the 75-second model and the red dot/line is the 30-second. Your guess was spot on, it appears as if the ‘instability’ forms right over PNG on the northern boundary of the region. There is no green dot/line after 2:16:15 as the 75-second model crashes at that point.

For a bit of background, we (CLEX) encountered some model instability when running the AUS2200 2.2km regional UM model during the northern NSW floods of Feb-Mar 2022. We didn’t have a lot of time to get the data ready, so we shortened the timestep and hoped (which worked, fortunately). Now that the pressing deadline of the original WACI hackathon has passed, we have some time to drill down and figure out what went wrong. I’m no UM expert, so I bought it up at the last AUS2200 meeting hoping to learn what to look for and what can be done about model instability. “Vertical winds over PNG” was the first hint, and that does seem to be the case. The next hint was that the orography may need to be smoothed out, so I’m hoping to learn what that would entail.

Dale

For the smoothing -

In the regional ancillary suite, make sure rgNN_rsNN_orog_srtm=True is set (use shuttle orography)

Then a further option will appear rgNN_rsNN_srtm_smooth. Select “121” for this, click on the parameter name to get the help text on how to increase the smoothing.

Once it’s smoothed enough to start running you could combine the dataset with the original orography if you like, using the smoothed orography only over the maritime continent and the original orography over australia - you’d need to use mule for this.

Thanks @Scott

It looks like the ancil suite is already producing a bunch of orographies with different smoothing factors applied, so that part is done. I don’t know how its picked which one to use though, it has linked qrparm.orog.mn to qrparm.orog.srtm.121x1, but SMOOTH_FACTORS="5 10 20 35 50 75 100". Does that mean it has done a single sweep of 121 smoothing on all points, or just over 1500m? In order to test this out, do I need to run the ancil suite again, or can I just link qrparm.orog.nm to a different file and run the forecast suite?

Dale

My recollections a bit fuzzy there - try linking one of the other files and see how it runs

I think I need a crash course in how orography works in the UM because I’m struggling to determine exactly what orography the model is loading in. It must be doing some kind of transform on the input data, because even simple things like “the highest and lowest points” don’t match between the ancil, the model dump created by reconfiguration. Most notably, the lowest point on the input dump is apparently around -105m(?). There is also a 200m discrepancy between the highest points in the model dump vs the ancil file. Any help would be appreciated.

Edit: Here is a picture, this is the northern-most part of the region. The actual location of the -105m point (the arrow points to it) is 5 grid points (11km) in from the northern boundary. The eastern-most edge of the image is about 1100km from the eastern boundary of the domain. Maybe an interpolation/regridding gone wrong?

More edit: It took me far too long to realise that the instability occurs more or less at that point as well. And as if to confirm it, I ran the model with 100 sweeps of 121 smoothing and it still crashed.

I’ve finally had some time to revisit this, and with ideas from @Scott and @MartinDix, I think I’ve narrowed the issue down to orography blending. Scott was right, in that the ERA5Land-generated orography is replaced by the ancil orography during the reconfiguration. The exception to this is on the very edge of the model, where the ancil orogoraphy is blended with the ERA5Land orography. I assume that this is to avoid discontinuities on the boundary. In AUS2200, the blending is configured to take the first 9 boundary rows/columns from the ERA5Land orography, then weight it against the ancil orography by 0.8, 0.6, 0.4 and 0.2 for the next 4 rows, then use the ancil orography for the rest of the domain. This appears to be some sort of default, as there is not much in the way of optional configuration that changes this. As per my last post, the instability occurs in 5 rows in from the boundary, a region that is entirely ERA5Land orography.

The main issue with that is it appears to have put Mt Bosavi in the wrong spot. Its the only topographical feature I can see that comes even close to what the orography is doing in that area. Its a collapsed extinct volcano about 2.5km high, with a 1km deep crater. The crater is very exaggerated he ERA5Land-derived orography, extending 100m below sea-level. I’ve made a thing that shows the issue here: AUS2200instability.pptx

Does anyone have any advice on how to proceed?

We should have the following orography sources:

  1. Base ERA5 orography, computed from ERA5 geopotential
  2. Orography used when converting from GRIB to UM formats, generated by the ancillary suite
  3. Orography of the target Aus2200 domain, generated by the ancillary suite

The ERA5 GRIB data gets interpolated to the boundary domain using the reconfiguration. As part of this the geopotential field (1) gets replaced with the boundary domain orography ancillary (2) so we don’t need to think about the raw ERA5 data. This is then processed by a second stage of reconfiguration to the model domain, which would be where the blending comes in I’d think.

The orography in the outer band of your plots should be coming from (2), so you should be able to manipulate it either through the ancillary suite or manually. Did you end up using extra smoothing for the model domain orography (3)? Consider doing the same for the boundary orography (2) - I would expect slide 2 to be showing the same data for both the model and boundary domains (except maybe the boundary domain is at 0.1 degree rather than 0.0198 degree grid spacing).

Yes, I did end up smoothing the model domain orography, however, that had no effect, as it turned out that it just wasn’t used. I did intentionally cut off the boundary domain in 2, but yes, it is continuous across the boundary. I guess I’m more concerned that there is this significant topological feature that just shouldn’t be there. No amount of smoothing will fix that.

What I really would like to understand is how the geopotential on slide 6 (that I’ve mis-labeled as orography and now realise I’ve set unuseful colour mapping on) becomes the orography in the second plot. The reconfiguration can’t be generating features that fine out of one or two pixels, so it must be attempting to match this geopotential data to known topological features. It just seems to have put a mountain that’s rather important to us in this case about 0.2 degrees too far south. I’m sure smoothing the boundary domain will bring make the crater a bit more sensible, but it still shouldn’t be in the model domain at all.

Are the initial conditions and ancillary files available somewhere that’s visible to me? (e.g. ly62)

Yep, ancils are in /scratch/ly62/dr4292/cylc-run/u-cp145, the initial condition are in /scratch/ly62/dr4292/cylc-run/u-cs142. Please let me know if you need me to chmod or something.

Here’s the SRTM reference orography data at the wonky location. It’s picking up the volcano correctly but something odd happened in the ridgeline to the SE


Summary
plt.figure(figsize=(15,10))
ax = plt.axes(projection=ccrs.PlateCarree())

srtm = iris.load_cube('/g/data/access/projects/access/data/SRTM/pp/S07.0E142.5.pp')
srtm_x = xarray.DataArray.from_iris(srtm)
srtm_x.plot(ax=ax, vmin=0, vmax=3000, cmap='Reds', add_colorbar=True)

srtm = iris.load_cube('/g/data/access/projects/access/data/SRTM/pp/S07.0E143.0.pp')
srtm_x = xarray.DataArray.from_iris(srtm)
srtm_x.plot(ax=ax, vmin=0, vmax=3000, cmap='Reds', add_colorbar=False)

srtm = iris.load_cube('/g/data/access/projects/access/data/SRTM/pp/S07.0E143.5.pp')
srtm_x = xarray.DataArray.from_iris(srtm)
srtm_x.plot(ax=ax, vmin=0, vmax=3000, cmap='Reds', add_colorbar=False)

srtm = iris.load_cube('/g/data/access/projects/access/data/SRTM/pp/S07.5E142.5.pp')
srtm_x = xarray.DataArray.from_iris(srtm)
srtm_x.plot(ax=ax, vmin=0, vmax=3000, cmap='Reds', add_colorbar=False)

srtm = iris.load_cube('/g/data/access/projects/access/data/SRTM/pp/S07.5E143.0.pp')
srtm_x = xarray.DataArray.from_iris(srtm)
srtm_x.plot(ax=ax, vmin=0, vmax=3000, cmap='Reds', add_colorbar=False)

srtm = iris.load_cube('/g/data/access/projects/access/data/SRTM/pp/S07.5E143.5.pp')
srtm_x = xarray.DataArray.from_iris(srtm)
srtm_x.plot(ax=ax, vmin=0, vmax=3000, cmap='Reds', add_colorbar=False)

plt.savefig('srtm.png')

/scratch/ly62/dr4292/cylc-run/u-cp145/share/data/ancils/aus2200/era5/orog_srtm/qrparm.orog.srtm is showing more reasonable orography

I guess the coarse default orography is being selected for the era5 to UM conversion rather than the SRTM orography. Looking at the logs this seems to be the case

$ grep 'orog' /scratch/ly62/dr4292/cylc-run/u-cs142/work/20220219T1200Z/ec_um_recon_000/pe_output/umgla.fort6.pe0
Special treatment for orography
Ancil filename     : /home/563/dr4292/cylc-run/u-cp145/share/data/ancils/aus2200/era5/qrparm.orog
OPEN: File /home/563/dr4292/cylc-run/u-cp145/share/data/ancils/aus2200/era5/qrparm.orog to be Opened on Unit 13 Exists
IO: Open: /home/563/dr4292/cylc-run/u-cp145/share/data/ancils/aus2200/era5/qrparm.orog on unit  13
CLOSE: File /home/563/dr4292/cylc-run/u-cp145/share/data/ancils/aus2200/era5/qrparm.orog Closed on Unit 13

(qrparm.orog should be instead the srtm version)

I’d recommend modifying /scratch/ly62/dr4292/cylc-run/u-cp145/share/data/ancils/aus2200/era5/ancil_filenames to point the orography to the SRTM version of the file and see if that improves the result. Presumably there’s a better way to do this in the suite that’s gone wrong somewhere

Yes, I see, the ERA5Land geopotential has nothing to do with it. According to the ancil suite, the orography is derived from /g/data/access/TIDS/UM/ancil/atmos/master/orography/GLOBE30/v1/qrparm.orog, which has the same weird ridgeline going on. Looks like it should be using SRTM from the start. Perhaps the ERA5->AUS2200 reconfiguration step should use SRTM orography instead?
OriginalOrog

Yeah using SRTM (100m gridlength) is a good idea at these resolutions

1 Like

@Scott that seems to have sorted it. I switched to the SRTM orography at the ERA5 → UM conversion, and the model did not fall over at the 2 hour mark. Fingers crossed we can revert back to the 75s timestep for future AUS2200 runs. Thanks everyone for your help!

1 Like

The Met Office GAL9 ancillary suite uses the newer GMTED orography rather than GLOBE30, ~access/umdir/ancil/atmos/master/orography/GMTED_RAMP2/v2/qrparm.orog. Still 30" resolution.

This is largely based on SRTM, https://www.usgs.gov/coastal-changes-and-impacts/gmted2010 and fixes the PNG problem.

GLOBE30 was pretty awful over Australia as well and we used a version patched with the GA data.

Thanks @MartinDix , that’s really helpful. The STRESS2023 AUS2200 runs are well underway now and its a bit late to change the orography from SRTM. Good news is that its running at a 75 second timestep and there is no sign of model instability so far.

1 Like