How to estimate subglacial flow accumulation

Here's a quick and easy way to make maps of subglacial water accumulation using TopoToolbox. This example uses Bedmap2 surface and bed elevations for for Thwaites Glacier.



This example requires

Step 1: Create a DEM of hydraulic head

You may know the nominal geo coordinates of Thwaites Glacier, but I sure don't, so I use scarloc.

[thwaiteslat,thwaiteslon] = scarloc('thwaites glacier');

The equation for static head is pretty simple, but I wrote it up in a function called bedhead so all the data importing and whatnot only requires one line. If you enter a single lat,lon location, a Bedmap2-derived static head DEM is centered on that location with a buffer of specified width in kilometers around that that location. A 600 km wide buffer around thwaiteslat,thwaiteslon creates a 1200 km wide map centered on Thwaites:

[X,Y,head] = bedhead(thwaiteslat,thwaiteslon,600);

Step 2: Turn the bedhead DEM into a GRID object

TopoToolbox uses GRID objects and FLOW objects. Turning X, Y, and head into a GRID object is a piece of cake with GRIDobj:

DEM = GRIDobj(X,Y,head);

For water routing to work, all sinks must be filled. In some cases, sinks are locations where subglacial lakes form. In other places, sinks simply represent errors in Bedmap2:

DEMf = fillsinks(DEM);

Step 3: Calculate flow accumulation

The TopoToolbox flowacc function calculates the total number of pixels contributing to each pixel in the DEM. You must convert the filled DEMf to a FLOW object, then call flowacc:

FD = FLOWobj(DEMf);
Accum = flowacc(FD);

You've now created a simple flow accumulation map. Show it like this:


You don't see much in that map because the grid is 1200 pixels wide whereas the figure is only ~500 pixels wide. All we see are a few aliased pixels where single-pixel-wide streams exist. Thus, it helps to dilate the image. TopoToolbox has a dilate function which is similar to Matlab's imdilate, but handles NaNs better. Here we dilate to a 5x5 neighborhood:

cb = colorbar;
ylabel(cb,'total contributing pixels')

Spatially-variable basal melt

The flowacc function was designed to allow spatially-variable precipitation weight, and this option works just as well for spatially-varying subglacial basal melt. By entering a map of basal melt, flowacc weights the value of each contributing pixel relative to the maximum and minimum values in the basal melt rate map.

Here I'm using geothermal flux as a crude stand-in for a spatial distribution of basal melt.

HeatFlux = heatflux_interp(X,Y);

% Set heat flux to zero except where there is grounded ice:
HeatFlux(bedmap2_interp(X,Y,'icemask','nearest')~=0) = 0;

% Create GRID object for basal melt:
BasalMelt = GRIDobj(X,Y,HeatFlux);

% Calculate flow accumulation
Accumw = flowacc(FD,BasalMelt);

Prettier displays

TopoToolbox has a nice version of imagesc which applies hillshading. But you may have noticed the parula colormap is not intuitive for water, so if you want something a bit more recognizable as ice and water, try <h cmocean function. I'm setting the first value in the colormap equal to the highest value so the ocean will look like a proper dark blue sink.

cm = cmocean('-ice');
cm(1,:) = cm(end,:);

The imageschs function applies hill shading to the hydrostatic head map by default, but I'd rather apply hillshading to the underlying bed topography and let color represent flow accumulation. The imageschs function wants bed topography as a GRID object, so start by turning Bedmap2 bed topography into a GRID object:

DEMbed = GRIDobj(X,Y,bedmap2_interp(X,Y,'bed'));

To illuminate some of the small-scale flow features consider taking the log or cube root of accumulation when you display it:

colorbar off

s = scarlabel({'Thwaites Glacier','Pine Island Glacier','Dotson Ice Shelf'},...

Turning grids into streams

Perhaps you'd rather have water routes as line objects rather than gridded flow accumulation values at each pixel. Turning flow accumulation grids into lines is easy, but you'll need to determine a threshold value that says, "wherever flow accumulation exceeds a certain value, consider that a river and ignore everything else." Let's use a threshold of 200 to capture many streams (use a higher threshold to only include larger streams) and plot the streams as orange lines:

thresh = 200;
S = STREAMobj(FD,Accumw>thresh);

hold on
p = plot(S,'-','linewidth',1,'color',rgb('orange'));

Just in case it isn't clear in that image, those are orange line objects plotted on top of the flow accumulation grid. The advantage here is you can plot the line objects on some other map, see:

axis equal tight

And for fun, plot flowlines from seed locations:

% Create a grid of seed locations:
[latseed,lonseed] = psgrid('thwaites glacier',1200,40);

% Get axis limits now because calling flowline might adjust them:
ax = axis;

% Plot flowlines:
f = flowline(latseed,lonseed,'plotxy','color',rgb('dark green'),'linewidth',.1)

% Make the plot pretty:
axis(ax); % Sets the axis limits back to the way they were
uistack(s,'top') % moves the text labels to the top of the stack
f =

Author info and citing these tools

This example was written by Chad A. Greene of the University of Texas at Austin's Institute for Geophysics (UTIG), January 2016, updated August 2016.

This page was written as a tutorial to get you started with some of the tools in AMT and TopoToolbox. If you use TopoToolbox, it'd be nice to cite:

Schwanghart, W., and D. Scherler. "Short Communication: TopoToolbox 2-MATLAB-based software for topographic analysis and modeling in Earth surface sciences." Earth Surface Dynamics 2.1 (2014): 1.

If you use Antarctic Mapping Tools, the nicest hat tip you can give me is a citation:

Greene, Chad A., David E. Gwyther, and Donald D. Blankenship. "Antarctic mapping tools for Matlab." Computers & Geosciences (2016).

And of course, if you use any publicly-available dataset, please the dataset(s) accordingly.