Ice sheet model output

This notebook demonstrates how to load, process and plot output from an ensemble of simulations of the Antarctic Ice Sheet from the last interglacial to present.

  1. 2020 by Torsten Albrecht (torsten.albrecht@pik-potsdam.de) | Potsdam Institute for Climate Impact Research (PIK)

data: Albrecht, Torsten (2019): PISM parameter ensemble analysis of Antarctic Ice Sheet glacial cycle simulations. PANGAEA, https://doi.pangaea.de/10.1594/PANGAEA.909728

More details on the ensemble can be found here.

[1]:

# load a package called intake, which is used to load the data.
import intake

1. Lazily load the model data

Selected outputs from the model ensemble are stored in two zarr directories in the google cloud. The data are cataloged in intake catalogs to simplify loading.

The zarr directory ‘present’ contains fields that describe the state of each ensemble member at the end of the simulation, i.e. the simulated ‘present-day’ state.

The zarr directory ‘mask_score_time_series’ contains time series of an integer mask (indicating ice-free/grounded/floating/ocean) for each ensemble member.

The cell below first loads the intake catalog from a github repo, then loads data from these two zarr directories in turn. However, it loads the data ‘lazily’, meaning that the data does not enter into memory untill we need them to, and ‘present’ and ‘mask_score_time_series’ are really just pointers to where the data are held, so that when we want to do a computation the machine will know where to look.

[2]:

# load the intake catalog
cat = intake.open_catalog('https://raw.githubusercontent.com/ldeo-glaciology/pangeo-pismpaleo/main/paleopism.yaml')

# load each of the zarr diretories contained in the intake catalog
present  = cat["present"].to_dask()
mask_score_time_series  = cat["mask_score_time_series"].to_dask()

2. Let’s take a look at the resulting xarrays

present

dimensions

present is a 7-dimensional xarray:

It has dimensions corresponding to each of the four ensemble parameters: par_esia, par_ppq, par_prec and par_visc

It has one time dimension, but this is only one element long so it can be mostly ignored.

It has two spatial dimensions, x and y.

Each dimensionn has a corresponding coordinate, which provides the values for the dimensions (i.e. the parameter values, the UPS X and Y). There is also an addition pair of coordinates called lat and lon, which provide the geographic positions of the grid cells.

variables

It contains variables corresponding to ice thickness (thk), ice surface elevation (usurf), ice speed (velsurf_mag), an integer mask (mask; indicating ice-free/grounded/floating/ocean), and bedrock uplift rate (dbdt).

It also contains a variable called ‘score’, which quantifies how well that ensemble member matched observationally constrain ice-sheet extent.

the xarray can be viewed by simply typing it name and executing the cell:

[3]:

present

[3]:

<xarray.Dataset>
Dimensions:      (par_esia: 4, par_ppq: 4, par_prec: 4, par_visc: 4, time: 1, x: 381, y: 381)
Coordinates:
    lat          (y, x) float64 dask.array<chunksize=(191, 191), meta=np.ndarray>
    lon          (y, x) float64 dask.array<chunksize=(191, 191), meta=np.ndarray>
  * par_esia     (par_esia) float64 1.0 2.0 4.0 7.0
  * par_ppq      (par_ppq) float64 0.25 0.5 0.75 1.0
  * par_prec     (par_prec) float64 0.02 0.05 0.07 0.1
  * par_visc     (par_visc) float64 1e+20 5e+20 2.5e+21 1e+22
  * time         (time) float64 50.0
  * x            (x) float64 -3.04e+06 -3.024e+06 ... 3.024e+06 3.04e+06
  * y            (y) float64 -3.04e+06 -3.024e+06 ... 3.024e+06 3.04e+06
Data variables:
    dbdt         (time, y, x, par_esia, par_ppq, par_prec, par_visc) float64 dask.array<chunksize=(1, 381, 381, 1, 4, 4, 4), meta=np.ndarray>
    index        (par_esia, par_ppq, par_prec, par_visc) int64 dask.array<chunksize=(4, 4, 4, 4), meta=np.ndarray>
    mask         (time, y, x, par_esia, par_ppq, par_prec, par_visc) int8 dask.array<chunksize=(1, 381, 381, 1, 4, 4, 4), meta=np.ndarray>
    score        (par_esia, par_ppq, par_prec, par_visc) float64 dask.array<chunksize=(4, 4, 4, 4), meta=np.ndarray>
    thk          (time, y, x, par_esia, par_ppq, par_prec, par_visc) float64 dask.array<chunksize=(1, 381, 381, 1, 4, 4, 4), meta=np.ndarray>
    topg         (time, y, x, par_esia, par_ppq, par_prec, par_visc) float64 dask.array<chunksize=(1, 381, 381, 1, 4, 4, 4), meta=np.ndarray>
    usurf        (time, y, x, par_esia, par_ppq, par_prec, par_visc) float32 dask.array<chunksize=(1, 381, 381, 1, 4, 4, 4), meta=np.ndarray>
    velsurf_mag  (time, y, x, par_esia, par_ppq, par_prec, par_visc) float32 dask.array<chunksize=(1, 381, 381, 1, 4, 4, 4), meta=np.ndarray>
Attributes:
    NCO:              4.6.8
    command:           /p/tmp/albrecht/pism18/pismOut/pism_paleo/pism1.0_pale...
    history:          Thu Dec  5 15:34:39 2019: ncatted -O -a history,global,...
    parameter_space:  {'visc': [1e+20, 5e+20, 2.5e+21, 1e+22], 'sia_e': [1.0,...

mask_score_time_series is also 7-dimensional xarray:

It has the same dimensions as present (above), but in len(mask_score_time_series.time) returns 125, indicating that this dataset has multiple time slices.

[4]:

mask_score_time_series

[4]:

<xarray.Dataset>
Dimensions:   (par_esia: 4, par_ppq: 4, par_prec: 4, par_visc: 4, time: 125, x: 381, y: 381)
Coordinates:
  * par_esia  (par_esia) float64 1.0 2.0 4.0 7.0
  * par_ppq   (par_ppq) float64 0.25 0.5 0.75 1.0
  * par_prec  (par_prec) float64 0.02 0.05 0.07 0.1
  * par_visc  (par_visc) float64 1e+20 5e+20 2.5e+21 1e+22
  * time      (time) float64 -1.24e+05 -1.23e+05 -1.22e+05 ... -2e+03 -1e+03 0.0
  * x         (x) float64 -3.04e+06 -3.024e+06 -3.008e+06 ... 3.024e+06 3.04e+06
  * y         (y) float64 -3.04e+06 -3.024e+06 -3.008e+06 ... 3.024e+06 3.04e+06
Data variables:
    index     (par_esia, par_ppq, par_prec, par_visc) int64 dask.array<chunksize=(4, 4, 4, 4), meta=np.ndarray>
    mask      (time, y, x, par_esia, par_ppq, par_prec, par_visc) int8 dask.array<chunksize=(125, 381, 381, 1, 4, 4, 4), meta=np.ndarray>
    score     (time, par_esia, par_ppq, par_prec, par_visc) float64 dask.array<chunksize=(125, 4, 4, 4, 4), meta=np.ndarray>
Attributes:
    NCO:              4.6.8
    command:           /p/tmp/albrecht/pism18/pismOut/pism_paleo/pism1.0_pale...
    history:          Thu Dec  5 16:44:54 2019: ncatted -O -a history,global,...
    parameter_space:  {'visc': [1e+20, 5e+20, 2.5e+21, 1e+22], 'sia_e': [1.0,...

3. Let’s plot the ice velocity field of one ensemble members

To select one of the ensemble members you need to choose which of the different choices for each paramater you want to plot. This example picks this randomly.

[5]:

present.usurf.isel(par_esia=3,par_ppq = 2, par_prec= 1,par_visc = 2).plot(size = 10)

[5]:

<matplotlib.collections.QuadMesh at 0x7f5964493520>
../../../_images/repos_ldeo-glaciology_pangeo-glaciology-examples_04_paleo_PISM_9_1.png

4. Do some calculations using all ensemble members

xarrays are great for dealing with high-dimensional data like this ensemble.

4a. calculate the standard deviation of thickness across the ensemble for each location

[6]:

mean_std = present.thk.std({'par_esia','par_visc','par_ppq','par_prec'},keep_attrs=True)
mean_std.plot(size = 10)

[6]:

<matplotlib.collections.QuadMesh at 0x7f59643717f0>
../../../_images/repos_ldeo-glaciology_pangeo-glaciology-examples_04_paleo_PISM_12_1.png

4b. Compute anomalies

Compute the difference between several ensemble members and the ensemble-wide mean computed above, then plot the results as an array of plots.

[7]:

mean_thk = present.thk.mean({'par_esia','par_visc','par_ppq','par_prec'},keep_attrs=True)
thickness_anomaly = (present.thk.isel(par_ppq=3,par_prec=2)-mean_thk)
thickness_anomaly.attrs['long_name'] = 'thickness anomaly'
thickness_anomaly.attrs['units'] = 'm'
thickness_anomaly.plot(x='x',y='y',col='par_esia',row='par_visc');

[7]:

<xarray.plot.facetgrid.FacetGrid at 0x7f59645f4700>
../../../_images/repos_ldeo-glaciology_pangeo-glaciology-examples_04_paleo_PISM_14_1.png

5. Plot the time slice data.

5a. plot a time series of ice shelf area

[8]:

cellArea = mask_score_time_series.x.attrs['spacing_meters'] * mask_score_time_series.y.attrs['spacing_meters']

[9]:

import xarray as xr
ice_shelf_area = (xr.where(mask_score_time_series.mask.isel(par_esia=1, par_ppq=2, par_prec=2, par_visc=3) == 3 ,1 ,0 ).sum(['x','y'])*cellArea)
ice_shelf_area.attrs['units'] = 'm^2'
ice_shelf_area.attrs['long_name'] = 'ice shelf area'
ice_shelf_area.time.attrs['long_name'] = 'simulation time'
ice_shelf_area.time.attrs['units'] = 'years BP'
p = ice_shelf_area.plot(size=8)
../../../_images/repos_ldeo-glaciology_pangeo-glaciology-examples_04_paleo_PISM_18_0.png

5b. Choose one of the ensemble members and plot masks for several time slices.

[10]:

mask_score_time_series.mask.isel(par_esia=1, par_ppq=2, par_prec=2, par_visc=3,
                       time=slice(0,124,15)).plot(x='x',y='y',col='time',col_wrap=3,size = 5)

[10]:

<xarray.plot.facetgrid.FacetGrid at 0x7f59545ce310>
../../../_images/repos_ldeo-glaciology_pangeo-glaciology-examples_04_paleo_PISM_20_1.png

6. Plot a map of the most deglaciated each location becomes across the whole ensemble and all time.

[11]:

mask_score_time_series.mask.isel(par_visc=2).max({"par_esia","par_ppq","par_prec","time"}).plot(size = 10)

[11]:

<matplotlib.collections.QuadMesh at 0x7f59643be2e0>
../../../_images/repos_ldeo-glaciology_pangeo-glaciology-examples_04_paleo_PISM_22_1.png 
[ ]: