MATLAB Examples

bedhead documentation

The bedhead function returns static pressure head at an ice sheet base in units of meters freshwater equivalent. This function calculates only a simple first-order approximation and does not consider effects such as velocity or strain or really anything fancy.



This function requires no special toolboxes and can be used for any glacier or ice sheet DEM, but if you do not supply your own DEM, the bedhead function will automatically use the Bedmap2 DEM. If you do not supply a DEM of your own you'll need the Bedmap2 Toolbox for Matlab and Antarctic Mapping Tools, both of which are available on the Mathworks File Exchange site.


psi = bedhead
psi = bedhead(lat,lon)
psi = bedhead(x,y)
psi = bedhead(...,extrakm)
psi = bedhead(...,'surface',SurfaceDEM)
psi = bedhead(...,'bed',BedDEM)
psi = bedhead(...,'thickness',ThicknessDEM)
psi = bedhead(...,'rhoi',IceDensity)
psi = bedhead(...,'rhow',WaterDensity)
[X,Y,psi] = bedhead(...)


psi = bedhead returns a 6667x6667 1 km gridded static pressure head for Antarctica from Bedmap2.

psi = bedhead(lat,lon) returns 1 km gridded static pressure head from Bedmap2 bounded by a polar stereographic quadrangle just large enough to encompass all the points in lat,lon. This syntax follows the bedmap2_data function.

psi = bedhead(x,y) returns 1 km gridded static pressure head from Bedmap2 bounded by a polar stereographic quadrangle just large enough to encompass all the polar stereographic (ps71 meters) coordinates x,y. This syntax follows the bedmap2_data function.

psi = bedhead(...,extrakm) adds a specified extra number of kilometers around x,y or lat,lon if the syntax above is used. This adds a buffer around your data, or if you would like a 500-km wide grid centered on a single point (75°S,100°E) let extrakm equal 250, i.e., bedhead(-75,100,250),

psi = bedhead(...,'surface',SurfaceDEM) specifies a gridded surface elevation DEM. If you specify a surface DEM you must also specify a BedDEM or a ThicknessDEM.

psi = bedhead(...,'bed',BedDEM) specifies a gridded bed elevation DEM. If you specify a bed elevation DEM you must also specify a SurfaceDEM or a ThicknessDEM.

psi = bedhead(...,'thickness',ThicknessDEM) specifies a gridded ice thickness DEM. If you specify a thickness DEM you must also specify a BedDEM or a SurfaceDEM.

psi = bedhead(...,'rhoi',IceDensity) specifies ice density. If you specify your own DEM you may also define IceDensity as a spatially-variable grid, say if you want to fully model variations in column- averaged density due to firn air content. Default value is 917 kg/m^3.

psi = bedhead(...,'rhow',WaterDensity) specifies density of fresh water. Default value is 1000 kg/m^3.

[X,Y,psi] = bedhead(...) returns gridded polar stereographic meters X,Y corresponding to the psi grid.

Example 1: Hydrostatic potential under Kamb Ice Stream:

Get a 600-km wide map of hydrostatic potential under Kamb Ice Stream:

[kamblat,kamblon] = scarloc('kamb ice stream');
[X,Y,psi] = bedhead(kamblat,kamblon,300);

And plot with pcolor:

shading flat
axis equal
scarlabel('Kamb Ice Stream','fontangle','italic')
cb = colorbar;
ylabel(cb,'hydraulic head (m)')

Example 2: Define your own DEM:

Suppose you have some arbitrary 1 km grid, 1000 km wide, centered on Byrd Camp. For this example I'll use psgrid to create such a grid:

[lat,lon] = psgrid('byrd camp',1000,1);

To make this example easy I'm using Bedmap2 surface and bed elevations as a starting point:

sfz = bedmap2_interp(lat,lon,'surface');
bed = bedmap2_interp(lat,lon,'bed');

But Bedmap2 has its uncertainties. Let's say it's 10 m for the surface and the bed uncertainty is given by:

bedunc = bedmap2_interp(lat,lon,'beduncertainty');

% Create surface and surface and bed DEMs with random values added to simulate
% uncertainties (crude approach, I know):

sfz2 = sfz + 10*randn(size(sfz));
bed2 = bed + bedunc.*randn(size(bed));

% Pressure head for the DEM with random uncertainties added is then given by:

psi2 = bedhead('surface',sfz2,'bed',bed2);
mapzoom('byrd camp',1000,'inset','nw')


There is only one simple mathematical expression in this function, which I adapted from Section 3.3 of

Fricker, H.A. and T. Scambos. "Connected subglacial lake activity on lower Mercer and Whillans ice streams, West Antarctica, 2003?2008." Journal of Glaciology 55.190 (2009): 303-315,

Fricker and Scambos and just about everyone else who has done work in this area cite Shreve, who seems to be the originator of the model:

Shreve, R. L. "Movement of water in glaciers." Journal of Glaciology 11 (1972): 205-214.

Author Info

This function and supporting documentation were written by Chad A. Greene of the University of Texas at Austin's Institute for Geophysics (UTIG), January 2016. I did not come up with the model for subglacial hydraulic head. If you use this function please cite the Shreve 1972 paper listed above.