MATLAB Examples

iceflex_interp documentation

The iceflex_interp function performs spatial interpolation to find local "coefficients" of ice flexure using the model presented by David Vaughan's 1995 JGR paper, Tidal flexure at ice shelf margins. It's worth noting that Vaughan's paper revived the model from Holdsworth's 1969 Journal of Glaciology paper, Flexure of a floating ice tongue.

The model described by Vaughan and Holdsworth is a 1D beam model, but here I'm implementing it in a semi-two-dimensional fashion. The 1D beam model considers distance from grounding line, ice thickness, and bulk material properties of the ice. The semi-two-dimensional model used in this function does not attempt to model two-dimensional stresses within the ice but instead calculates 1D flexure as a function of distance from any interpolation point to the nearest grounded ice. This approach has not been validated by either observations or comparison to other models, but in my tests it seems to perform quite well as a first-order approximation. In most cases, errors in ice thickness and/or exact location of the grounding line are likely larger contributors to errors than the semi-two-dimensional approximation.



C = iceflex_interp(lati,loni)
C = iceflex_interp(xi,yi)
C = iceflex_interp(...,'E',YoungsModulus)
C = iceflex_interp(...,'method',InterpolationMethod)
C = iceflex_interp(...,'rhosw',SeawaterDensity)
C = iceflex_interp(...,'nu',PoissonsRatio)
C = iceflex_interp(...,'grid',GridStruct)


C = iceflex_interp(lati,loni) returns ice flexure coefficient(s) at georeferenced location(s) given by lati,loni, which can be a point, line, or grid.

C = iceflex_interp(xi,yi) solves ice flexure for polar stereographic coordinates xi,yi. Coordinates are automatically determined by the islatlon function.

C = iceflex_interp(...,'E',YoungsModulus) specifies an elastic modulus in pascals. Default value is 0.88e9 from Vaughan's 1995 Ice flexure at ice shelf margins paper but it should be noted that Marsh et al. found a better fit with 1.4e9 via InSAR on Beardmore Glacier. Choice of Young's modulus has a relatively small impact on vertical displacements.

C = iceflex_interp(...,'method',InterpolationMethod) specifies a spatial interpolation method and can be any of the methods accepted by interp2. Default method is 'cubic'. Choice of interpolation method has a nearly negligible impact on vertical displacements.

C = iceflex_interp(...,'rhosw',SeawaterDensity) specifies seawater density in kg/m^3. Default value is 1028.

C = iceflex_interp(...,'nu',PoissonsRatio) specifies Poisson's ratio. Default value is 0.3.

C = iceflex_interp(...,'grid',GridStruct) lets you define a grounded ice mask and ice thickness grid of your own. GridStruct is a structure and must have these fields:

GridStruct.x         (a grid containing polar stereographic x in meters)
GridStruct.y         (a grid containing polar stereographic y in meters)
GridStruct.ground    (a logical grid where ones denote grounded ice)
GridStruct.thickness (an ice thickness grid with thickness values in meters)

Alternatively, GridStruct may be defined using geo coordinates:       (a grid of latitudes)
GridStruct.lon       (a grid of longitudes)
GridStruct.ground    (a logical grid where ones denote grounded ice)
GridStruct.thickness (an ice thickness grid with thickness values in meters)

The name of the GridStruct variable can be anything (e.g., GridStruct, G, bananas), but the fields within GridStruct are case sensitive and must be exactly .x, .y, .ground, and .thickness, or .lat, .lon, .ground and .thickness. Default GridStruct is generated automatically from Bedmap2.

Example 1: Profile along a track

Helen Fricker and Laurie Padman wrote a real nice GRL paper on ICESat-observed flexure of the Institute Ice Stream region of southern Ronne Ice Shelf. Here's a quick recreation of their Figure 3a and 3e using Bedmap2 data instead of ICESat observations. First get reference track number 218 at 100 m spacing using reftrack clip it to the latitude range 81.1°S to 80.3°S:

[lat,lon] = reftrack(218,'spacing',100);

% clip to latitude range of interest:
ind = lat>-81.1 & lat<-80.3 & lon<0;
lat = lat(ind);
lon = lon(ind);

If Bedmap2 grounding line location and ice thickness are correct here, relative ice flexure should be about this much:

C = iceflex_interp(lat,lon);

Quickly eyeballing the modeled tide anomalies for the laser campaigns shown in Figure 3e I get these values:

tide.l3c = 2.3;
tide.l3a = 0.35;
tide.l2a = 0;
tide.l3b = -2.6;

Now we can use bedmap2_interp to get nominal elevations and recreate Figure 3a:

z = bedmap2_interp(lat,lon,'surfw');

hold on
plot(lat,z+tide.l3b*C,'.','color',[0.98 0.45 0.02],'markersize',8);

axis([-81.1 -80.3 100 200])
text(-80.5,180,'Laser 2a','color','b')
text(-80.5,170,'Laser 3a','color','c')
text(-80.5,160,'Laser 3b','color',[0.98 0.45 0.02])
text(-80.5,150,'Laser 3c','color','m')
ylabel 'ICESat elevation (m)'
xlabel 'Latitude'

And Figure 3e should look about like this:

hold on
plot(lat,tide.l3b*C,'.','color',[0.98 0.45 0.02],'markersize',8);

axis([-81.1 -80.3 -3 3])
text(-80.5,180,'Laser 2a','color','b')
text(-80.5,170,'Laser 3a','color','c')
text(-80.5,160,'Laser 3b','color',[0.98 0.45 0.02])
text(-80.5,150,'Laser 3c','color','m')
ylabel 'Elevation anomaly (m)'
xlabel 'Latitude'

Example 2: Map the hydrostatic line

The iceflex_interp function can handle points, profiles, or grids. To make a map of an estimated hydrostatic line, start by defining a 120 km wide grid of 250 m resolution with psgrid and then get a 2D coefficient of ice flexure C:

[X,Y] = psgrid(-81.05,-67.05,120,0.25,'xy');
C = iceflex_interp(X,Y);

The 2D coefficient of ice flexure will look like this in polar stereographic coordinates:

shading interp
hold on
xlabel 'eastings (m)'
ylabel 'northings (m)'
cb = colorbar;
ylabel(cb,'coefficient of tidal flexure');

For added context we can add the Bedmap2 grounding line. Currently the bedmap2 function plots the whole continent's worth of grounding line and will adjust axis limits so we'll get the axis limits before plotting the grounding line, plot the grounding line, then set the axes to the original limits:

ax = axis;

The hydrostatic line may be defined where the coefficient of ice flexure is unity, and thus Matlab's contour function can be used to get the coordinates of everywhere C is equal to one. Included in AMT is a C2xyz function which will make the contour matrix a little more user friendly:

C1 = contour(X,Y,C,[1 1]);
[C1x,C1y] = C2xyz(C1);

If you look closely at C, you'll notice that on the ice shelf values tend to wobble slightly around unity, so C1 may contain multiple possible hydrostatic lines. We only want the hydrostatic lines closest to the grounding line, so it takes a little manual picking to see we only want cells 3 and 13 for this hydrostatic line. We'll plot the hydrostatic line in red:

for k = [3 13]

Now how well does our modeled hydrostatic line compare to ICESat observations? Let's recreate Figure 4a of Brunt et al.'s 2011 Journal of Glaciology paper. It looks like they used MOA 2004 and this work matches ICESat grounding zones:


Now here are our hydrostatic lines estimated from Bedmap2. Note we have to transform from polar stereographic to geo coordinates:

for k = [3 13]
    [C1lat,C1lon] = ps2ll(C1x{k},C1y{k});

Example 3: Beardmore Glacier

Here's a quick example of what flexure might look like at the grounding line of Beardmore Glacier if Bedmap2 were accurate. Marsh et al. showed that Bedmap2 thickness data may be incorrect on the order of hundreds of meters here. Let's take a look, and we'll use their best-fit Young's modulus value 1.4 GPa. Below we use psgrid to center a 40 km-wide grid of 100 meter spatial resolution on the outlet of Beardmore Glacier:

[lat,lon] = psgrid(-83.58,172.34,40,0.1);
C = iceflex_interp(lat,lon,'E',1.4e9);

Multiply C by a 55 cm tidal displacement and map it:

p = pcolorm(lat,lon,C*0.55);
caxis([0 0.55])
cb = colorbar;
ylabel(cb,'estimated tidal flexure (m)')
bedmap2 gl

Example 4: Effect of elasticity

Using the grid from the example above, investigate the effect of elasticity by comparing the 1.4 GPa Young's modulus to Vaughan's (default) 0.88 GPa. Below I'm using Stephen Cobeldick's brewermap function to get the nice red-white-blue colormap:

C088 = iceflex_interp(lat,lon);

cb = colorbar;
ylabel(cb,'\Delta z_{1.4 GPa} - \Delta z_{0.88 GPa} (m)')
caxis([-0.1 0.1])
bedmap2 gl

I'm not sure if the figure above makes it clear, but choice of elastic modulus does not have a large impact on vertical displacements in the grounding zone. The much more important and often uncertain factors are ice thickness and grounding line location.

Example 5: DIY DEM, or build your own grid

A recent Nature Geoscience paper by Greenbaum et al. describes a floating section of Totten Glacier which was previously thought to be grounded. Here's the area:

[x,y] = psgrid(-67.1731,116.6206,100,0.1,'xy');
C = iceflex_interp(x,y);

p = pcolor(x,y,C);
shading interp
hold on
ax = axis;
caxis([0 1])

That section in the middle is actually connected by ocean. Using [cx,cy] = coord to click around on the map I get these coordinates as approximate lateral bounds on the trough:

cx = [2244581;2244581;2244581;2243996;2244581;2247505;2248675;2247505;2247797;2247505];
cy = [-1121432;-1123771;-1126695;-1129326;-1131081;-1130788;-1129034;-1125525;-1122601;-1120847];

h = plot(cx,cy,'ro-');

To build your own grid, you'll need gridded coordinates, which we already have in the form of x and y. We'll also need a grounded ice mask and thickness data. Begin with the ice mask, which we'll create by modifying Bedmap2 grounded ice data. This is a logical dataset, so we'll use nearest-neighbor interpolation:

mask_bedmap2 = bedmap2_interp(x,y,'icemask','nearest');

Grounded ice in Bedmap2 is denoted by zeros, so that's how we start defining our mask:

ground = mask_bedmap2 == 0;

To define the trough as ungrounded, set the grid values inside the cx,cy polygon to zero:

ground(inpolygon(x,y,cx,cy)) = 0;

Now our grounded ice mask is defined. It looks like this, where white denotes grounded ice:

axis xy

See, the trough in the middle is now defined as not grounded.


We'll also need thickness data, which you can get from your favorite DEM. For convenience, here we'll use Bedmap2 for the thickness data:

thickness = bedmap2_interp(x,y,'thickness');

Now combine grid data into a structure which we'll call G:

G.x = x;
G.y = y;
G.thickness = thickness;
G.ground = ground;

Re-solve using our custom grid and replace the previous pcolor data with our new custom flexure:

C_custom = iceflex_interp(x,y,'grid',G);
caxis([0 1])

See, the ice in that trough should now rise and fall with the tides.


If you're interested in ice flexure and grounding lines, these fine articles are certainly worth a read:

Brunt, Kelly M., Helen A. Fricker, and Laurie Padman. "Analysis of ice plains of the Filchner?Ronne Ice Shelf, Antarctica, using ICESat laser altimetry." Journal of Glaciology 57.205 (2011): 965-975.

Fricker, Helen Amanda, and Laurie Padman. "Ice shelf grounding zone structure from ICESat laser altimetry." Geophysical Research Letters 33.15 (2006).

Greenbaum, J. S., et al. "Ocean access to a cavity beneath Totten Glacier in East Antarctica." Nature Geoscience 8.4 (2015): 294-298.

Holdsworth, Gerald. "Flexure of a floating ice tongue." Journal of Glaciology 8 (1969): 385-397.

Marsh, Oliver J., et al. "Grounding-zone ice thickness from InSAR: inverse modelling of tidal elastic bending." Journal of Glaciology 60.221 (2014): 526-536.

Vaughan, David G. "Tidal flexure at ice shelf margins." Journal of Geophysical Research 100.B4 (1995): 6213-6224.

Author Info

The iceflex_interp function and supporting documentation were written by Chad A. Greene of the University of Texas at Austin's Institute for Geophysics (UTIG), 2015.