1. Introduction
Greenland outlet glaciers have been delivering highly variable and increasing quantities of ice to the adjacent seas from the vast inland ice-sheet reservoir, leading to increased rates of global sea-level contribution. This insight is due largely to satellite remote sensing that has revealed glacier surface velocity (e.g. Reference Joughin, Abdalati and FahnestockJoughin and others, 2004; Reference Luckman, Murray, de Lange and HannaLuckman and others, 2006), calving-front position (e.g. Reference Sohn, Jezek and van der VeenSohn and others, 1998; Reference Csatho, Schenk, van der Veen and KrabillCsatho and others, 2008), thinning/thickening rates (e.g. Reference Howat, Joughin and ScambosHowat and others, 2007) and ice discharge (Reference Rignot and KanagaratnamRignot and Kanagaratnam, 2006). Aerial photography (Reference BauerBauer, 1968; Reference Carbonnell and BauerCarbonnell and Bauer, 1968; Reference FountainFountain, 1983; Reference Fastook, Brecher and HughesFastook and others, 1995; Reference Motyka, Fahnestock, Howat, Truffer, Brecher and LuethiMotyka and others, 2008), terrestrial photogrammetry (Reference Brecher and ThompsonBrecher and Thompson, 1993; Reference Amundson, Truffer, Lüthi, Fahnestock, West and MotykaAmundson and others, 2008) and repeated airborne laser altimeter surveys (Reference KrabillKrabill and others, 2000; Reference AbdalatiAbdalati and others, 2001) have also proven effective in quantifying glacier dynamics.
Analysis of glacier-change time series on annual to decadal scales has been conducted by, for example, Reference Sohn, Jezek and van der VeenSohn and others (1998), Reference Joughin, Abdalati and FahnestockJoughin and others (2004), Reference Luckman, Murray, de Lange and HannaLuckman and others (2006) and Reference Csatho, Schenk, van der Veen and KrabillCsatho and others (2008). However, observing ice-dynamical quantities at a frequency shorter than ∼1 week has remained at the technological frontier. Approaches employing in situ GPS (e.g. Reference Hamilton and WhillansHamilton, and Whillans, 2000), auto-theodolite reflectors (e.g. Reference Harper, Humphrey, Pfeffer and LazarHarper and others, 2007), terrestrial laser scanning (Reference Schwalbe, Maas, Dietrich and EwertSchwalbe and others, 2008) and radar interferometry (Reference Chapuis, Rolstad and NorlandChapuis and others, 2010) have begun to fill this observational gap.
Velocities determined from time-lapse photogrammetry have the advantage of providing high spatial and temporal resolution sampling. In situ GPS or theodolite reflectors, moving with the glacier, provide discrete velocity points along a trajectory in the Lagrangian frame. Time-lapse photogrammetry offers an alternative observational frame, enabling velocity to be derived through fixed locations in the Eulerian frame.
Time-lapse photogrammetry in glaciology is not new. Reference Harrison, Echelmeyer, Cosgrove and RaymondHarrison and others (1992) described the difficulties of feature tracking in the calculation of absolute speed. Recent advances in commercially available digital camera technology, solid-state data storage and GPS have greatly facilitated terrestrial time-lapse photogrammetry. Reference EvansEvans (2000) summarized image-processing stages for motion tracking in sequential images by image-matching and -filtering techniques. Reference Motyka, Hunter, Echelmeyer and ConnorMotyka and others (2003) and Reference O’Neel, Echelmeyer and MotykaO’Neel and others (2003) used time-lapse photogrammetry to determine the horizontal terminus-area changes at LeConte Glacier, Alaska, USA. Applying photogrammetric techniques, Reference Maas, Dietrich, Schwalbe, Bässler and WestfeldMaas and others (2006) and Reference DietrichDietrich and others (2007) observed ice velocity and tidal effects on the lower reach of Jakobshavn Isbræ, West Greenland.
This study focuses on the development of techniques for calculating ice displacement and velocity variability from single monoscopic camera stations installed on land beside glacier termini in West Greenland. We first discuss the various problems encountered. We then present techniques for image stabilization. Next, solutions in image-feature matching, displacement calculation and the conversion of image-space displacements to metric velocity are described. Error budgets and temporal/spatial resolution analyses are made. Finally, glacier dynamics insights are provided, with a more comprehensive application of this technique in glacier–climate sensitivity studies deferred to future work. This work is to serve as the methodological foundation for subsequent quantitative work of the Extreme Ice Survey (EIS).
2. Time-Lapse Cameras
2.1. Data
Digital single-lens reflex cameras with fixed focal length lenses were installed on land beside the termini of marine-terminating glaciers (Rink, Store, Umiamako and Jakobshavn glaciers) in West Greenland during May–June 2007 (Fig. 1; Table 1). Table 2 lists site information. For the sunlit period of the yeaar, photos were acquired each hour. The interval timer reduced the number of photos, depending on the time of year, such that during winter the image acquisition rate decreased to a minimum of one or two frames per day.
Fig. 1. (a) Time-lapse camera locations; (b) station set-up at Umiamako B.
Table 1. EIS time-lapse station hardware specification
Table 2. EIS time-lapse camera site information
2.2. Problems encountered
The effects of falling snowflakes and rain droplets, snow and rain on the ground, lens internal reflections and enclosure window optical distortion all contributed to variable amounts of image clarity. Clouds blocked some glacier observations, leading to an image loss rate of 8–12% (Table 3). West Greenland fjord cloudiness is minimized in the lower atmosphere by the force and dryness of katabatic winds.
Table 3. Daily image loss rates (%)
Pixels affected by extremes in illumination, i.e. under- and overexposure, are not always completely compensated by automatic camera exposure. The shadows and ‘hot spots’ lack brightness variation and these can confound feature identification. Camera malfunctions have occurred due to excessive vibration presumably caused by winds. Temporary sticking of the camera shutter has occurred at <−30°C. The camera timer has failed in more than one instance, leading to data gaps of more than 10 days. Wild animals have scratched enclosure windows and elsewhere have left tooth-marks. Thick plastic conduit successfully protects power cables. Manual image inspection is performed to remove images rendered useless by these problems. The above problems account for a loss rate of 6–10% (Table 3). Including image loss owing to fog, the total loss rate is 15–20%.
2.3. Spatial resolution
Image spatial resolution varies across the glacier scene, especially given the oblique 75–85°-from-vertical view of these installations. We simulate the spatial resolution pattern by assuming a horizontal plane. The pixel ground coverage (Fig. 2) is 0.62, 0.93 and 1.24 m with distances 2, 3 and 4 km, respectively, from 20 mm nominal focal length and a 6.2 μm camera detector array element size. Because of the oblique view, what would be squares in satellite imagery are trapeziums (Fig. 2). Displacement magnitudes must increase to remain resolvable as distance from the camera increases.
Fig. 2. The effects of perspective: (a) Landsat velocity determined from image pairs spanning 7–26 July 2002; (b) time-lapse photogrammetry derived velocities for the interval 22–23 July 2007. Arrows are (a) velocity and (b) pixels d−1 image displacements.
2.4. Landsat imagery
Glacier surface velocities calculated using the offset of visible surface features between repeat, co-registered Enhanced Thematic Mapper Plus (ETM+) images (Reference Ahn and HowatAhn and Howat, 2009; Reference Howat, Box, Ahn, Herrington and McFaddenHowat and others, 2010) are available to this study. We use four image pairs for Rink Isbræ, one pair for Store Gletscher, five pairs for Umiamako and five pairs for Jakobshavn Isbræ.
3. Photogrammetric Methods and Results
3.1. Area-based image matching
Image area-based ‘correlation matching’ is implemented to identify feature motion between successive images. Glacier motion is evident as image pixel column and row offsets, referred to here as ‘displacements’. We use only the midday photos to provide the most continuous coverage of the year and with the highest sun illumination angles.
Square image subsets or ‘chips’ or ‘windows’ are extracted from a ‘reference image’ and a ‘search image’. A reference chip is coordinated within a search chip larger than the reference chip. As the reference chip is nudged in the image column and row directions within the search chip, a correlation coefficient is calculated from all the combinations of reference chip placement, producing a correlation grid. The maximum in the correlation map corresponds with a ‘conjugate point’ that is assumed to be the same point in the reference and search images. The pixel offset in the image column and row direction thus yields the feature displacement. A quadratic surface fitting around the correlation peak is performed to obtain sub-pixel precision. Further details of this technique are provided by Reference SchenkSchenk (1999), and for a glaciological application see Reference Scambos, Dutkiewicz, Wilson and BindschadlerScambos and others (1992).
3.2. Registration
Camera motion caused by wind gusts, camera platform settling and other factors (e.g. thermal expansion/contraction of the camera platform) introduces rotation and image column and row shifts between image pairs spanning time. Image registration compensates for camera motion, aligning images so that image comparison through time can accurately and precisely represent glacier displacements. Registration is a first processing step (flow chart, Fig. 3) and involves matching techniques described later. We refer to motion-compensated images as having been ‘registered’ (Image registration in flow chart, Fig. 3). In extreme camera-motion compensation, registration offsets are up to 38 pixels. Typical offsets are, however, in the 0–8 pixel range (Table 4). Nonzero average rectification offsets indicate camera ‘exterior orientation’ drift.
Fig. 3. The work flow of mono time-lapse camera processing to derive glacier velocity. AOI: area of interest; TL: time lapse; EOP: exterior orientation parameter.
Table 4. Pixel offset statistics from camera motion compensation
An approach for compensating camera motion involves a registration technique that transforms one image to an equivalent perspective in the other, i.e. from a chosen source image to a destination image (Reference BrownBrown, 1992). Registration, as such, requires conjugate points on stationary features, in our case, sharp rock edges and crack features. Feature identification is problematic when the land-surface properties change over time owing to changing surface wetness reflectivity, snow-cover accumulation/ablation, permafrost creep and vegetation. Whether the solar illumination is direct or diffuse, i.e. causes shadows or not, is important. Direct beam illumination introduces a time-varying orientation and length of shadows. Diffuse illumination is preferred, as it illuminates what would otherwise be obscured by shadows. Direct illumination in one image and diffuse illumination in another creates an undesirable image pair.
We solve the registration problem by choosing a distribution of control points (CPs) on land both on the near and far sides of the glacier. A sample of 25 neighboring points around each CP is generated. The correlation-matching procedures then identify conjugate points. The maximum correlation point among the 25 CP neighbor points is chosen for the registration via a transformation model. When foreground and ‘farground’ CPs are available, as in the case of the Rink, Umiamako and Store glacier sites, a ‘projective’ transformation model (Reference Slama, Theurer and HenriksenSlama and others, 1980; Reference Zitová and FlusserZitová and Flusser, 2003) is used. A simpler linear conformal transformation model (Reference Slama, Theurer and HenriksenSlama and others, 1980) is used for the Jakobshavn Isbræ site, where only foreground CPs are available.
For projective transformation, the outlier-insensitive random sample consensus (RanSaC) correlation procedure (Reference Fischler and BollesFischler and Bolles, 1981) is applied. The image is divided into quadrants, with CPs selected from each quadrant. Then all possible combinations of sampling with at least 1 CP from each quadrant are used to ensure a wide distribution of CPs and thus a more reliable transformation solution.
3.3. Multi-Image/Multi-Chip (MIMC) matching
We use ‘area-based’ image matching to (1) compensate camera motion by identifying motionless features on land and (2) determine ice displacements from multiple images spanning time. We use a set of 15 image pairs to obtain robust matching solutions. Here we discuss this MIMC approach (MIMC feature tracking in flow chart, Fig. 3).
The first MIMC procedure involves multiple chip sizes. Optimal chip size varies with feature displacement rate, viewing geometry, illumination, etc., and is determined empirically by balancing the requirements of feature uniqueness and the geometric distortion that grows with chip size. We use square chips of 30, 40 and 50 pixels on a side. The use of multiple chip sizes increases the matching success rate by enabling disqualification of infrequently occurring spurious correlation maxima, i.e. occurring only once among a larger set of valid results.
The second MIMC procedure is the use of multiple images (Image enhancement in flow chart, Fig. 3). In this case, before the image-matching procedures are followed, image enhancements are applied. First, a principal-component intensity image is extracted from the red, green and blue (RGB) bands (Reference Gonzales and WoodsGonzales and Woods, 2002). Then a 0° and 90° directional edge enhancement applied to the intensity image produces two more images. The 0° and 90° enhanced images are high-pass filtered using a 3 x 3 kernel (−1 −1 −1; −1 8 −1; −1 −1 −1) to extract image information corresponding to sharp features such as rifts and crevasse edges. Ultimately, five enhanced images are generated from the original image (Fig. 4).
Fig. 4. (a) Original image, (b) intensity image and (c) 0° edge-enhanced, (d) 90° edge-enhanced, (e) 0° edge + high-frequency enhanced and (f) 90° edge + high-frequency enhanced images of a heavily crevassed portion of an image of Rink Isbræ.
Given that we apply three search-chip sizes to five enhanced images, MIMC identifies up to 15 matches per chip, in 30 images. Among these, the displacement results vary. Furthermore, the highest correlation is not always the most accurate solution. To identify the most accurate solution, a filtering scheme identifies a cluster of displacement results using a drop-in-the-bucket ‘voting cell’, also akin to a cluster. The displacement within the cell with the maximum count is chosen as the most accurate matching result. The width of the voting cell is 2 pixels, and the dimension of the region where voting cells are counted is defined by the spatial range of the 15 resultant displacement vectors. Another set of voting cells with a 1-pixel shift in image column and row directions is sampled to prevent a cluster from being divided at a cell boundary. See MIMC application to satellite imagery (Reference Ahn and HowatAhn and Howat, 2009).
3.4. MIMC versus conventional matching
Conventional correlation matching results for three bands (RGB) and three chip sizes (30 × 30, 40 × 40 and 50 × 50) are compared with MIMC results at image ‘posts’ having a 50-pixel image interval (Fig. 5). The conventional correlation result is illustrated for the red band with a 30 × 30 reference chip size in Figure 5a and b. MIMC results are displayed in Figure 5c and d. We find that without any filtering (e.g. using a correlation threshold), MIMC results in fewer outliers than conventional matching. Enlarged Figure 5b and d show how erroneous displacements are excluded in the MIMC cluster filtering. It is impossible to exclude spurious vectors in conventional matching when only one displacement vector is determined.
Fig. 5. Illustration of 30 × 30 reference chip red-band displacement results from conventional matching (a, b) and from MIMC (c, d). Dashed arrows indicate outliers removed by the voting cell approach. Dashed arrows in (d) indicate cases where conventional matching fails.
Using MIMC, we can also identify what type of image enhancement – i.e. intensity, edge, high-pass filtered and what chip size (30, 40 and 50 pixels on a side) – contributes most successfully (highest correlations) to the final displacement results. Intensity image pairs in all chip sizes consistently result in the highest correlations and therefore the highest contribution to the chosen displacement (Table 5). Edge and high-pass filtered image pairs contribute only 2–5% to the selected displacement retrievals. The directional or high-pass filtered images seem to outperform the principal-component images in hazy cases.
Table 5. Contributions (%) from different image types to final matching selections
We applied additional filtering to MIMC results to eliminate outliers not removed by the voting cell approach, that is, (1) an empirically determined maximum possible displacement threshold of 40 pixels d−1, and (2) only accepting displacements that satisfy (Di − D mean) < D std, where Di is the ith displacement and D mean and D std are respectively the mean and standard deviation of a 5 × 5 window from a′displacement grid’ centered at the ith point. The displacement grid refers to the 50 × 50 pixel grid of posts (Fig. 5).
To verify the robustness of MIMC results, we first compare the MIMC principal-component intensity-image correlation values with that from conventional matching of RGB intensity images for the three chip sizes. Table 6 lists the proportion of posts where the correlation is higher for our approach compared with conventional matching. A correlation ratio of 0.5 indicates no improvement. When matching using principal-component intensity-image pairs, we find an average increase in correlation ratio of 0.25 more than conventional matching. Second, we compare the displacements between the nine conventional matching results and the MIMC filtered results. When the difference in image row and column displacement exceeds 2 pixels, we categorize the conventional matching results as ‘fail’ (Table 6). The failure rate is 3% or smaller.
Table 6. Comparison between MIMC intensity and conventional matching
There may be some consistency among the matching results on an image band basis. For example, red outperforms the other bands at Umiamako. The blue band performs best at Rink and Store glaciers (Table 6). A more comprehensive study of performance by bands may yield conclusive color-based matching improvements.
3.5. Daily displacement calculations
Once all images are co-registered, the matching procedures are applied to the ice-covered image areas in rectified images spanning 1–4 days. The image area that provides optimum results is usually near the glacier terminus, where the largest displacements are and distance to the camera is lowest. Pixel assortments nearer the camera contain more distinct feature information, i.e. better resolution. To develop a time series of glacier motion, we keep the target always on ice by choosing an image position represented by a relatively small (20 × 20 pixel) area close to the terminus yet far enough upstream that any front retreat does not leave the target off the glacier.
Daily displacement calculations are based on photos selected from within 3 hours of the same hour of day to minimize uncertainty caused by changing illumination (Pair selection in flow chart, Fig. 3). Changes that occur in images at times t 0 and t 1 are considered the average displacement at (t 0 + t 1)/2.
3.6. Exterior orientation parameter (EOP) calculation
EOPs describe the camera position in three-dimensional space, i.e. three X, Y and Z position coordinates and three rotation angles ω, φ and κ. As such, an object in three dimensions can be projected to image space, and the object in the image can also be projected to a surface or a plane (Reference Slama, Theurer and HenriksenSlama and others, 1980) (EOP available in flow chart, Fig. 3).
The ‘single photo resection’ that calculates the EOPs requires more than three ground control points (GCPs) with known positions and corresponding image points. The GCPs are taken from Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) imagery and an ASTER-derived digital elevation model (DEM).
The single photo resection is possible because the camera position is known by GPS and by identifying the conjugate features (e.g. mountain tops and valley intersections) in time-lapse images and in the ASTER image draped over the ASTER terrain model. We calculate that the uncertainties of ASTER imagery and elevation data are sources of error resulting in angular uncertainty of ±0.3°, corresponding to a ±30 m error at a 3 km distance, assuming the GPS position is error-free. The image space and ground space can be linked so that points in image space can be transferred to the ground space. We choose points in ground space and transfer them to the camera image space so that point intervals have a regular 120 m ground spacing (Fig. 2).
3.7. Image-space displacement conversion to velocity
We develop two approaches for converting pixel displacement to velocity (Displacement to velocity in flow chart, Fig. 3). The first is to calculate the ‘scale factor’. In this case, the scale factor is a function of ground footprint size and ice-flow direction. The footprint size increases with distance to the target, and changes with the changing position of displacements. In the case of the 20 mm focal length lens, we calculate a footprint size per pixel of 6.2 m at 1 km distance at the image center. The footprint size increases to 9.4 m at the image edges for the same distance. The distance between the camera station and the glacier can be approximated using satellite imagery. Because the EOPs and a DEM are available, we calculate the 3-D position where a ray from the camera intersects the DEM surface. First, the ray from the camera intersects an initial height, h 0, where the surface height, h 1, is then calculated. If (h 1−h 0) is greater than a threshold of 5 m, then the first surface height becomes the initial height for iteration, which stops when dh is 5 m or less. If the procedure does not converge after 50 iterations, the average height of the surface from the digital terrain model is assumed. The effect of a ±15 m ASTER DEM uncertainty on scale assignment is 5% of the scale factor, if the usual target is 2 km away and the camera elevation is 300 m above the DEM.
The flow direction ‘into’ the frame partially compensates a ‘foreshortening’ effect that results in shorter than actual displacement when the camera captures motions that are not parallel to the camera plane. If the camera is installed perpendicular to the flow direction, as in the case of Store B and Umiamako B, this effect is less noticeable.
Another approach to assign scale is to carry out the projection of pixel displacements to the orthographic satellite-derived velocity plane. With the availability of Landsat-derived velocities and EOPs, we set a vertical plane that contains each Landsat velocity vector. Projecting time-lapse displacements to the vertical plane enables comparison of velocities at equivalent posts, shown later.
3.8. Error budget
The two major sources of uncertainty in displacement calculations stem from image matching and projective transformation model-fitting errors. Reference DietrichDietrich and others (2007) analyzed the uncertainty in matching based on the steepness of parabola fitting in correlation, and concluded that the uncertainty is <0.1 pixel for a foreground reference target and 0.1–0.3 pixels for the glacier feature. Reference Skvarca, Raup and De AngelisSkvarca and others (2003) more conservatively estimated a 0.5-pixel uncertainty for correlation matching and a 1.5-pixel uncertainty in registration. We compute a 0.23 ± 0.08 average registration root-mean-square error (RMSE) from projective transformation (Table 7). Since our area of interest is located between registration control points, the uncertainty in registration procedure can be estimated to be less. Correlation-matching uncertainty is difficult to quantify, since quadratic surface fitting uses a correlation value around a 3 × 3 or 5 × 5 pixel peak. Here a 0.3-pixel uncertainty for matching can be considered conservative. When uncertainties are square-summed for each step, the daily displacement uncertainties for our sites (Table 7) are 0.52, 0.46, 0.47, 0.45, 0.49 and 1.04 pixels, corresponding to ‘on the ground’ velocity uncertainties of 0.18 m d−1 at 1000 m distance at the image center and with ice motion perpendicular to the camera.
Table 7. Registration error statistics
3.9. Spatial and temporal resolution
The minimum resolvable displacement and its related time interval are functions of camera focal length, pixel size projected onto the ground, target distance, ice velocity and registration precision. A first necessary condition for motion to be resolved is that the displacement exceed the registration and matching uncertainty of 0.3 pixels. In our case, the registration error is roughly twice the matching error.
Pixel footprint size, derived from the target distance, ice velocity and time-lapse interval are related parameters. For relatively slow glaciers, a longer time interval or closer features have a similar effect. The resolvable lapsing time, given displacement Δp (pixel units), time interval Δd (days), pixel footprint F (length units), the uncertainty e (pixels) and the error budget ∊ b expressed as a percentage, is Δd = (∊/(Δp∊b)) × 100%. In the case of Rink Isbræ, which has 14–18 pixels d−1 displacement, with a calculated absolute uncertainty of 0.52 pixels, the displacement uncertainty is 20%. The resolvable lapsing interval is thus ∼4.5 hours. Halving the distance doubles Δp and thus halves the minimum time interval of resolvable displacement to ∼2.5 hours. A relatively slow glacier, flowing 2 m d−1, having 3.2-pixel displacement at a distance of 1000 m, with 0.62 m pixel size, requires >1 day for a displacement to exceed the uncertainty. The distance to target and the glacier speed have a similar effect on resolution.
The temporal resolution depends on how the image pairs are selected. Pairs for sub-daily displacement are applicable by searching before and after the specific hour. Even hourly intervals (4 hours, 6 hours, … 24 hours) produce displacement results with even (integer) time coordinates. Combining with odd hour intervals (5 hours, 7 hours, … 23 hours) can yield effective sub-hourly resolution.
3.10. Validation with Landsat-derived quantities
The camera-derived velocities are averaged over ±3 days and binned to be equivalent to the relatively sparse Landsatderived quantities. Landsat-derived velocities indicate velocity differences less than 1.2 m d−1, or 8%, compared with camera-derived velocities (Fig. 6). Landsat-derived velocities are relatively sparse in time, having variable time separation, i.e. 9–30 days.
Fig. 6. Daily velocity variability at a point near the front of, respectively, Rink Isbræ, Store Gletscher A, Store Gletscher B, Umiamako A, Umiamako B and Jakobshavn Isbræ. Vertical lines indicate major calving events.
3.11. Validation with camera pairs
Two cameras were placed at the Umiamako and Store sites (Table 2). The Umiamako (Store) cameras are separated by 2.9 km (1.7 km). The ‘B’ sites are further up-glacier. The average velocities from these two perspectives are equivalent on average between the sites, i.e. within the budgeted error (Table 8). However, the velocity variations from these differing perspectives do not have equivalent correlations: Umiamako R 2 = 0.77 and Store R 2 = 0.64. We attribute the differing correlations to the different perspective and flow direction with respect to the cameras. The perspective difference strongly influences illumination. The time series of running correlation indicates a fluctuation with many values having correlation greater than 0.9 and periods of correlation dropping out.
Table 8. West Greenland glacier velocity statistics from time-lapse photogrammetry
4. Glaciological Results
Figure 6 illustrates the much higher temporal frequency of camera-derived velocities than is available from multi-day satellite image pairs. The time interval in sequential camera-derived pairs is 1–4 days, translating to effective 12 hour temporal resolution (section 3.5). Daily observation frequency reveals coherent multi-day velocity trends and relatively abrupt changes. The vertical solid lines in Figure 6 indicate major calving events identified by visual inspection in sequential hourly time-lapse images. At Rink Isbræ on 19 August and 12 September 2007, the calving of ‘ice islands’ >1 km in length led to ∼10 day progressive accelerations of 18% and 15% (Fig. 6a and b), respectively. A 22 June 2007 calving event at Umiamako led to a progressive ∼6 day 17% acceleration (Fig. 6c). Jakobshavn Isbræ, on 3 and 26 July 2007, exhibited acceleration of 6% and 7%, respectively, over 7 days after calving (Fig. 6d and e). These accelerations need to be analyzed in terms of the area of calving, the event location with respect to the medial flowline, the glacier front shape, the abruptness of the calving given that some events last days before complete detachment, the existence of sea ice, whether the tongue is grounded or floating, glacier thickness changes, etc.
The 17–38% velocity range during the mid-June to late-September period considered in this study (Table 8) is larger than has previously been reported. In this section, we discuss this finding in the context of sub-annual variability on a glacier-by-glacier basis.
4.1. Jakobshavn Isbræ
Seasonality in flow speed for Jakobshavn Isbræ appears to have emerged with the 1999–2003 progressive disintegration of its ice shelf. Between 1985 and 1986, Reference Echelmeyer and HarrisonEchelmeyer and Harrison (1990) found no measurable seasonally averaged speed variation. For 2004–07, Reference Joughin, Das, King, Smith, Howat and MoonJoughin and others (2008) attributed a maximum of 15% of speed fluctuations, at a point roughly 4 km from the terminus, to terminus-position changes, with maximum speed occurring when the front had retreated most. Our daily observations in 2007 indicate a gradual acceleration of 38% as the perennial ice shelf and glacier terminus disintegrated to a new retreated minimum. Reference Amundson, Fahnestock, Truffer, Brown, Lüthi and MotykaAmundson and others (2010) found that a multi-km long ice melange packed in the fjord slows the calving rate and regulates seasonal advances of the ice front, and therefore its speed and ice discharge.
4.2. Umiamako
According to Reference Howat, Box, Ahn, Herrington and McFaddenHowat and others (2010), in 1999 the Umiamako terminus position exhibited seasonal oscillations of a few hundred meters, advancing in the spring and then retreating with the clearing of the calved glacier ice and landfast sea-ice melange. Using approximately weekly Landsat imagery, in 2001 they observed a ∼50% speed increase, from 1.2 m d−1 in April to 2.5 m d−1 in mid-June, followed by a return to 1.3 m d−1 flow through mid-September. In 2003, Umiamako began a progressive 4 km retreat lasting at least until 2009. Ice-flow speed increased from a melt-season peak of 2.3 m d−1 before 2003, to 6.5 m d−1 in 2009, roughly a factor of three increase. We observe velocities consistent with the Umiamako post-2003 retreat acceleration, i.e. 4.8–7.4 m d−1 (Fig. 6; Table 8). In terms of seasonality, we compute a ∼36% range in June–September daily speeds. Thus, like Jakobshavn Isbræ, Umiamako has gained seasonality and speed during a period of retreat. We thus confirm that velocity variations seem strongly controlled by flow resistance imposed by glacier and sea-ice jamming at the front.
4.3. Store Gletscher
Similar to Jakobshavn and Umiamako, an acceleration (∼30%) is evident at Store Gletscher from the earliest observations of the year when an ice melange imposed flow resistance at the front (Fig. 6). Our observations ended before an end-of-melt deceleration could be confirmed.
4.4. Rink Isbræ
The ice-discharge rate of Rink Isbræ has been estimated to be roughly equivalent to that of Store Gletscher (Reference Weidick and BennikeWeidick and Bennike, 2007). We find a 3% higher mean flow velocity and a 17% higher minimum speed for Rink Isbræ than for Store Gletscher. Rink Isbræ calves larger bergs, some wider (>1 km) than they are deep (∼0.6 km), such that they do not immediately turn over. Rink Isbræ exhibits the most distinct multi-day speed variations among our sites (Fig. 6). We speculate that unlike Store Gletscher, Rink Isbræ calves ice islands because it is nearer flotation.
5. Conclusions
Daily time-lapse camera images can be used to precisely quantify glacier motion in image space, i.e. with pixels d−1 units. By projecting the camera view geometry onto a digital terrain model, it is possible to assign horizontal scale to converting displacements into velocities (m d−1). The ASTER image geolocation and elevation uncertainty results in angular uncertainty of ±0.3°, corresponding to planimetric ±30 m errors at 3 km distance. Comparison between terrestrial photogrammetric velocities and velocities derived from 15 m Landsat 7 imagery confirms the feasibility of retrieving velocities in orthometric coordinates.
Before pixel displacements and ice velocities are calculated, a series of procedures are applied, including manual image filtering of images confounded by weather (e.g. by fog or rain droplets on the camera enclosure window). Fortunately, fog is uncommon in fjords continually flushed by relatively warm and dry katabatic winds. Hardware failures have resulted in data gaps of >1 week, especially in the first year of operation before including a more robust timer. Camera motion is compensated using matching of motionless features on land on the near and far sides of the glacier that flows through the middle of the image. In the case of Jakobshavn Isbræ, with no motionless land features in the upper part of the image, a less robust ‘affine transformation’ approach is used to compensate camera motion, producing uncertainties roughly twice as large as the glacier situations with land in the foreground and farground. Fortunately, Jakobshavn Isbræ has the highest velocity in our survey (maximum 47.6 m d−1, average 37.6 m d−1), resulting in the lowest speed uncertainty (2%, or 3.9 m d−1) among the sites. The slowest glacier in our survey, Umiamako, flowing 5.9 m d−1 on average, has an associated mean speed uncertainty of 9%, or 0.5 m d−1.
Image enhancement includes principal-component intensity and 90° and 0° edge-detection filters. The use of multiple image enhancements and multiple image chip sizes in matching procedures produces distinctly better matching performance than matching single image pairs. Intensity image pairs produce 25% higher matching correlation values than red, green or blue greyscale image pairs. 90° and 0° edge and high–pass filtered image pairs contribute only a few percent to the displacement retrievals. The directional or high-pass filtered images seem to outperform the principal-component images only in hazy cases. A clustering technique, not possible with only one displacement vector, as in conventional matching, facilitates rejection of high-correlation false matches.
The red image ‘band’ outperforms the other (green and blue) bands at the Umiamako site. The blue band performs best (highest correlation) at Rink and Store glaciers. A more comprehensive study of performance by bands may yield conclusive image-band-based matching improvements.
At distances of ∼3 km and given the range of observed displacements, the uncertainty in matching becomes as large as the computed displacements, rendering the technique applicable to an area limited to the lowest few km of the glacier, the fastest areas of the glacier, or the areas nearest the camera. Higher camera placements are possible but are confounded by increased cloudiness and are harder to reach on foot.
Resolvable displacement is a function of (1) pixel footprint derived from target distance, (2) ice velocity and (3) time-lapse interval. For relatively slow glaciers, a longer time interval or closer features have a similar effect. In the case of Rink Isbræ, which has 14–18 pixels d−1 displacements, with a calculated absolute uncertainty of 0.52 pixels, the displacement uncertainty is 20%. The resolvable lapsing interval becomes ∼4.5 hours. Halving the distance halves the minimum interval of resolvable displacement to 2.5 hours. A relatively slow glacier, flowing 2 m d−1, having 3.2-pixel displacement at a distance of 1000 m, with 0.62 m pixel size, requires a >1 day time interval for a displacement to exceed the uncertainty.
The minimum resolvable displacement is a function of camera focal length, pixel size on the ground, distance to target, ice velocity, and registration precision. A first necessary condition for motion to be resolved is that the displacement exceeds the registration and matching uncertainty, here 0.3 pixels. In our case, the registration error is roughly twice as large as the matching error. Thus, the success of this approach depends most on compensating camera motions that result from wind, the settling of the ground, interference by wild animals, human visitors, etc.
A preliminary analysis of daily velocities suggests that major calving events are responsible for the largest day-to-day velocity changes, which are progressive accelerations requiring up to 10 days for maximum speed to be reached. After the calving-induced acceleration, the velocity gradually declines, presumably as resistance to flow rebuilds. A greater seasonality is evident, not just given the ability to observe more frequently than the at-best weekly satellite perspective, but as the fronts of marine-terminating glaciers retreat.
Terrestrial photogrammetry opens new possibilities for investigating the sensitivity of glacier dynamics to climate variability. Daily displacements retrieved from time-lapse camera images yield higher spatial/temporal resolution near the calving front than is typical of conventional remote-sensing techniques. Continuous GPS or auto-theodolite reflectors can provide higher temporal resolution, yet provide more complex Lagrangian point measurements, while the time-lapse imaging technique can yield Eulerian data with broader coverage, especially at the hazardous calving front where the survivability of in situ sensors is limited.
Acknowledgements
Major support for this research was provided by NASA grant NNX06AE31G managed by Ming-Ying Wei and The Ohio State University’s Climate, Water and Carbon (CWC) program managed by D. Alsdorf. This work has benefited from camera system development by J. Balog, W.T. Pfeffer and A. Lewinter. A. Behar and N. Humphrey provided electrical engineering assistance. Guidance from H. Brecher, D. Merchant and A. Schenk provided an important conceptual foundation for this work. We thank J. Orlowski and O.J. Hammeken for assistance with fieldwork. Major logistical support was provided by J. Buenning, M. Begnauld, S. Zager and R. Abbott of CH2MHillPolar Field Services, and E. Vaughn of the US Air National Guard. I. Howat helped obtain Landsat and ASTER imagery and contributed constructive ideas and discussions on the MIMC approach. Constructive comments from T. James and an anonymous reviewer improved the clarity of the manuscript. We thank Nikon and The North Face for EIS sponsorship. This is Byrd Polar Research Center contribution No. 1397.
