Comparing Model Resolution

Load and explore model solutions at a range of horizontal resolutions (1°, 0.25°, and 0.1°). We use snapshots of model output to explore how the velocity fields differ between the simulations. We also look at eddy kinetic energy, which highlights that coarser resolution simulations have less variability in their velocity fields.

import cosima_cookbook as cc
import matplotlib.pyplot as plt
import xarray as xr
import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cft
import IPython.display

Given that the high resolution simulations contain quite a lot of data, let’s fire up a Dask cluster to parallelise the computation.

from dask.distributed import Client

client = Client('127.0.0.1:38079')
client

Client

Cluster

  • Workers: 4
  • Cores: 16
  • Memory: 68.72 GB

In preparation for the plots, let’s load the land for cartopy.

land_50m = cft.NaturalEarthFeature('physical', 'land', '50m',
                                   edgecolor='black', facecolor='gray', linewidth=0.5)

Now instantiate a database session.

session = cc.database.create_session()

Pick three experiments to analyse. If you want to look at a list of all the experiments available, run the following snippet.

print(cc.querying.get_experiments(session)['experiment'].values)

expts = ['01deg_jra55v13_iaf', '025deg_jra55v13_iaf_gmredi6', '1deg_jra55v13_iaf_spinup1_B1_lastcycle']
titles = ['0.1° horizontal resolution', '0.25° horizontal resolution', '1° horizontal resolution']
sims = {}
for exp in expts:
    u = cc.querying.getvar(exp,'u',session, ncfile='ocean.nc')
    v = cc.querying.getvar(exp,'v',session, ncfile='ocean.nc')
    speed = cc.querying.getvar(exp,'speed',session, ncfile='ocean.nc')
    eke = (u-u.mean(dim='time'))**2 + (v-v.mean(dim='time'))**2

    sims[exp] = xr.merge([u, v, speed, eke.rename('eke')])

Let’s have a look at the data we just loaded into this dictionary.

sims
{'01deg_jra55v13_iaf': <xarray.Dataset>
 Dimensions:   (st_ocean: 75, time: 396, xu_ocean: 3600, yu_ocean: 2700)
 Coordinates:
   * st_ocean  (st_ocean) float64 0.5413 1.681 2.94 ... 5.511e+03 5.709e+03
   * xu_ocean  (xu_ocean) float64 -279.9 -279.8 -279.7 -279.6 ... 79.8 79.9 80.0
   * yu_ocean  (yu_ocean) float64 -81.09 -81.05 -81.0 -80.96 ... 89.92 89.96 90.0
   * time      (time) object 1985-01-14 12:00:00 ... 2017-12-14 12:00:00
 Data variables:
     u         (time, st_ocean, yu_ocean, xu_ocean) float32 dask.array<chunksize=(1, 7, 300, 400), meta=np.ndarray>
     v         (time, st_ocean, yu_ocean, xu_ocean) float32 dask.array<chunksize=(1, 7, 300, 400), meta=np.ndarray>
     speed     (time, st_ocean, yu_ocean, xu_ocean) float32 dask.array<chunksize=(1, 7, 300, 400), meta=np.ndarray>
     eke       (time, st_ocean, yu_ocean, xu_ocean) float32 dask.array<chunksize=(1, 7, 300, 400), meta=np.ndarray>,
 '025deg_jra55v13_iaf_gmredi6': <xarray.Dataset>
 Dimensions:   (st_ocean: 50, time: 300, xu_ocean: 1440, yu_ocean: 1080)
 Coordinates:
   * st_ocean  (st_ocean) float64 1.152 3.649 6.565 ... 5.034e+03 5.254e+03
   * xu_ocean  (xu_ocean) float64 -279.8 -279.5 -279.2 -279.0 ... 79.5 79.75 80.0
   * yu_ocean  (yu_ocean) float64 -81.02 -80.92 -80.81 ... 89.79 89.89 90.0
   * time      (time) object 1958-06-30 12:00:00 ... 2257-06-30 12:00:00
 Data variables:
     u         (time, st_ocean, yu_ocean, xu_ocean) float32 dask.array<chunksize=(1, 10, 216, 288), meta=np.ndarray>
     v         (time, st_ocean, yu_ocean, xu_ocean) float32 dask.array<chunksize=(1, 10, 216, 288), meta=np.ndarray>
     speed     (time, st_ocean, yu_ocean, xu_ocean) float32 dask.array<chunksize=(1, 10, 216, 288), meta=np.ndarray>
     eke       (time, st_ocean, yu_ocean, xu_ocean) float32 dask.array<chunksize=(1, 10, 216, 288), meta=np.ndarray>,
 '1deg_jra55v13_iaf_spinup1_B1_lastcycle': <xarray.Dataset>
 Dimensions:   (st_ocean: 50, time: 60, xu_ocean: 360, yu_ocean: 300)
 Coordinates:
   * st_ocean  (st_ocean) float64 1.152 3.649 6.565 ... 5.034e+03 5.254e+03
   * xu_ocean  (xu_ocean) float64 -279.0 -278.0 -277.0 -276.0 ... 78.0 79.0 80.0
   * yu_ocean  (yu_ocean) float64 -77.75 -77.51 -77.26 -77.01 ... 89.1 89.55 90.0
   * time      (time) object 2198-07-02 12:00:00 ... 2257-07-02 12:00:00
 Data variables:
     u         (time, st_ocean, yu_ocean, xu_ocean) float32 dask.array<chunksize=(1, 25, 150, 180), meta=np.ndarray>
     v         (time, st_ocean, yu_ocean, xu_ocean) float32 dask.array<chunksize=(1, 25, 150, 180), meta=np.ndarray>
     speed     (time, st_ocean, yu_ocean, xu_ocean) float32 dask.array<chunksize=(1, 25, 150, 180), meta=np.ndarray>
     eke       (time, st_ocean, yu_ocean, xu_ocean) float32 dask.array<chunksize=(1, 25, 150, 180), meta=np.ndarray>}

Having loaded the data, let’s plot some output from each simulation and look at the differences. We use xarray’s ability to index using dimension names to select the model level closest to the surface, st_ocean = 0, and the final timestamp. The resulting 2-D array is then plotted.

The following plots are snapshots of zonal velocity. All three simulations show roughly the same results, with westward and eastward flows occuring in approximately the same locations. However, the 0.1° and 0.25° simulations show mesoscale variability that is entirely lacking from the 1° simulation.

plt.figure(figsize=(10,12))

for counter, exp in enumerate(expts):
    ax = plt.subplot(3, 1, counter+1, projection=ccrs.Robinson(), label='{0}'.format(counter))

    ax.coastlines(resolution='50m')
    ax.add_feature(land_50m)
    gl = ax.gridlines(draw_labels=False)
    p1 = ax.pcolormesh(sims[exp].u.xu_ocean, sims[exp].u.yu_ocean,
                       sims[exp].u.sel(st_ocean=0, method='nearest').sel(time=sims[exp].time[-1]),
                       transform=ccrs.PlateCarree(),
                       cmap='RdBu_r', vmin=-0.8, vmax=0.8)
    ax.set_title(titles[counter])

ax_cb = plt.axes([0.85, 0.15, 0.03, 0.7])
cb = plt.colorbar(p1, cax=ax_cb, orientation='vertical')
cb.ax.set_ylabel('Zonal velocity (m/s)');
../_images/Model_Resolution_Comparison_0.png

The mesoscale variability hinted at in the above global figures is more easily seen by zooming into a regional scale. The plots below show the Atlantic sector of the Southern Ocean, including the Agulhas Current, and Drake Passage. The figures below also show that higher horizontal resolution support much higher flow speeds.

projection=ccrs.Mercator(central_longitude=0.0, min_latitude=-70.0, max_latitude=-20.0)

plt.figure(figsize=(12,12))

for counter, exp in enumerate(expts):
    ax = plt.subplot(3, 1, counter+1, projection=projection, label='{0}'.format(counter))

    ax.set_extent([-100, 50, -70, -20], crs=ccrs.PlateCarree())
    ax.coastlines(resolution='50m')
    ax.add_feature(land_50m)
    gl = ax.gridlines(draw_labels=True)
    gl.xlabels_top = False
    p1 = ax.pcolormesh(sims[exp].speed.xu_ocean, sims[exp].speed.yu_ocean,
                       sims[exp].speed.sel(st_ocean=0, method='nearest').sel(time=sims[exp].time[-1]),
                       transform=ccrs.PlateCarree(),
                       cmap='inferno', vmin=0, vmax=0.8)
    ax.set_title(titles[counter])

ax_cb = plt.axes([0.85, 0.15, 0.03, 0.7])
cb = plt.colorbar(p1,cax=ax_cb, orientation='vertical')
cb.ax.set_ylabel('Speed (m/s)');
../_images/Model_Resolution_Comparison_1.png

The higher resolution simulations exhibit much more variability in their velocity fields. This can be easily seen by looking at eddy kinetic energy, defined as

\[(u - \overline{u})^{2} + (v-\overline{v})^{2}\]

which picks out variability in the velocity fields. Once again we select the uppermost level, and the final timestamp.

This plot takes quite some time to produce, since it requires averaging the entire time series of u and v for each simulation.

projection=ccrs.Mercator(central_longitude=0.0, min_latitude=-70.0, max_latitude=-20.0)

plt.figure(figsize=(12,12))

for counter, exp in enumerate(expts):
    ax = plt.subplot(3, 1, counter+1, projection=projection, label='{0}'.format(counter))

    ax.set_extent([-100, 50, -70, -20], crs=ccrs.PlateCarree())
    ax.coastlines(resolution='50m')
    ax.add_feature(land_50m)
    gl = ax.gridlines(draw_labels=True)
    gl.xlabels_top = False
    p1 = ax.pcolormesh(sims[exp].eke.xu_ocean, sims[exp].eke.yu_ocean,
                       xr.ufuncs.log10(sims[exp].eke.sel(st_ocean=0, method='nearest').sel(time=sims[exp].time[-1])),
                       transform=ccrs.PlateCarree(),
                       cmap='inferno', vmin=-4.5, vmax=-0.5)
    ax.set_title(titles[counter])

ax_cb = plt.axes([0.85, 0.15, 0.03, 0.7])
cb = plt.colorbar(p1,cax=ax_cb, orientation='vertical')
cb.ax.set_ylabel('Log(EKE)');
../_images/Model_Resolution_Comparison_2.png

Download python script: Model_Resolution_Comparison.py

Download Jupyter notebook: Model_Resolution_Comparison.ipynb