Velocity model

We approach the problem of estimating using a Bayesian formulation that follows Reference TarantolaTarantola (2005). Inverse methods are, of course, well known, epitomized by the least-squares method and its regularized variants. Our motivation for using a Bayesian formulation is to apply a probabilistic approach to derive a generalized model of the desired quantity, in this case , and estimates of the uncertainties of the model parameters. There are three parameters that form the conceptual basis of this approach: the a posteriori conditional probability density function (PDF), or posterior, of the model, m, given the observed data; the a priori PDF, or likelihood, which relates the observations to the model; and the a priori PDF of the model, known as the prior, which incorporates prior expectations of all model parameters. The basic strategy of Bayesian inversion is then to represent the posterior, P(m|d), as a combination of the likelihood, P(d|m), and prior, P(m). Maximizing the posterior yields an expression for the model that is comparable with classical regularized least-squares formulations, but includes prior model estimates and a prior model covariance matrix. Though we use Bayesian inversion to infer a relatively simple physical model (3-D velocity fields) we note that the formulation of the posterior model is generalized and can be used to infer almost any linear or nonlinear model. Because Bayesian inverse methods are well developed (Reference TarantolaTarantola, 2005; Reference StuartStuart, 2010), the following derivation of the posterior model for is concise, meant only as an overview for readers unfamiliar with Bayesian methods.

Let m be a model for , such that , where is the set of all realizable models. From Bayes’ theorem, the posterior probability distribution is given as (Reference TarantolaTarantola, 2005)

Assuming a Gaussian model, the likelihood and prior are defined as (Reference TarantolaTarantola, 2005)

respectively, where m ₀ is the prior model estimate, and C _m and C _d are the prior model and data covariance matrices, respectively. Plugging Eqns (4) and (5) into Eqn (3) yields

where

The best-fit posterior model, , maximizes the posterior probability (i.e. minimizes β) such that (Reference TarantolaTarantola, 2005)

The first term in Eqn (8) is the posterior model covariance matrix:

which provides an estimate of the uncertainties in (Reference TarantolaTarantola, 2005). Higher amplitudes in the components of indicate higher uncertainty in . It is often desirable to encapsulate the error of the posterior model as

where tr is the trace operator.

Due to the size of most InSAR data, iterative approaches, which do not explicitly invert the parenthetical term on the right-hand side of Eqn (9), may be the only computationally tractable way to estimate . Therefore, when calculating it may be more practical to consider the case where , when Eqn (9) simplifies to

whose elements can be estimated for each independent pixel, assuming C _d is diagonal (see next subsection), instead of as an ensemble of interdependent pixels as required in Eqn (9). From Eqn (11), we can see that is a function of the viewing geometries (via G) and interferometric noise (via C _d).

Estimates of the uncertainty attributable to non-ideal viewing geometries are contained in the sensitivity matrix:

The diagonal terms of S are the sums of the squares of the components of the LOS vector, and the off-diagonal terms are the sums of the cross products of the LOS vector components. The off-diagonal components indicate the coupling between the respective inferred velocity field components that result from a non-ideal set of viewing geometries, while the diagonal components quantify how InSAR measurement errors propagate into the components of . An ideal set of viewing geometries can be generally described as having a constant, oblique incidence angle and full azimuthal coverage, with constant azimuthal spacing between flight lines. Differential incidence angles, inconsistent azimuthal spacing or incomplete azimuthal coverage in the viewing geometries leads to nonzero off-diagonal components and increased sensitivity to measurement noise. Sensitivity to measurement noise decreases with increasing amounts of data. To characterize the contribution of the LOS geometry to the model uncertainty, we define the variance term:

Readers familiar with GPS analysis might recognize Λ_g as the position dilution of precision (PDOP), the spatial component of the geometric dilution of precision (GDOP) (Reference Misra and EngeMisra and Enge, 2006).

We can glean some idea about the sensitivity of Λ_g to the number of independent data points and a constant incidence angle by considering an idealized set of p ≥ 3 viewing geometries. Let the members of a set of p LOS unit vectors be defined in a spherical coordinate system described by equispaced azimuth angles (ψ_k = 2πk/p, for k = 1, …, p) and a constant incidence angle, θ, measured relative to vertical. These simplifications yield

Therefore, Λ_g decreases as the square root of the number of observations for a given θ and is approximately constant, for a given p, over the range of incidence angles common in InSAR data.

Data covariance matrix

The data covariance matrix, C _d, can have contributions from atmospheric phase delay (e.g. Reference HanssenHanssen, 2001; Reference Emardson, Simons and WebbEmardson and others, 2003; Reference Lohman and SimonsLohman and Simons, 2005), interferometric decorrelation (e.g. Reference Rodriguez and MartinRodriguez and Martin, 1992; Reference Zebker and VillasenorZebker and Villasenor, 1992; Reference HanssenHanssen, 2001) and spatial dependences within the InSAR data. Because the spatial scales of our study areas are no more than a few kilometers, while the spatial wavelengths of atmospheric phase delays are ∼10 km (Reference Emardson, Simons and WebbEmardson and others, 2003) and we do not know the a priori spatial dependencies in the InSAR data, we assume that the data covariance matrix C _dis a function of only the interferometric SNR and is defined as

where is the phase variance for a given pixel in scene i. If any of the interferograms considered in the estimation share a common SAR scene (i.e. acquisition), the data covariance matrix will have off-diagonal components (Reference Emardson, Simons and WebbEmardson and others, 2003).

The phase variance can be estimated as a function of the interferometric correlation, γ_i , and N_i , the number of pixels in the incoherent averaging window (i.e. the number of independent looks) for scene i, using the Cramer–Rao bound (Reference Rodriguez and MartinRodriguez and Martin, 1992):

The interferometric correlation, γ_i , is defined as (e.g. Reference RosenRosen and others, 2000)

where s_z is the complex scattered signal in SAR image z, 〈·〉 indicates averaging over numerous realizations of the argument and * represents the complex conjugate. The phase in Eqn (17) is the interferometric phase, and the complement of the amplitude (ς_i = 1 − |γ_i |) is commonly called the interferometric decorrelation.

The correlation amplitude, sometimes called the coherency, provides an estimate of the interferometric noise in a given interferogram. Values near unity indicate a small amount of interferometric noise, whereas correlation values near zero mean the data are dominated by noise. The interferometric correlation is generally defined as a product, γ_i = [γ _n γ _b γ _z γ _t] _i , where the four independent components are due to noise, viewing geometry (perpendicular baseline), volumetric effects and temporal variations in the scatterers, respectively (e.g. Reference RosenRosen and others, 2000).

Prior model covariance matrix

We define the prior model covariance matrix for a given study area based on the physical processes being studied. Choosing reduces Eqn (8) to the classical weighted least-squares problem, the least computationally demanding approach. This choice implicitly assumes that points within the prior model are independent of all other points, which is not necessarily true for many geophysical problems.

In this study, we expect the surface velocity to vary smoothly in space, such that ∇² m ₀ ≈ 0. Therefore, we define the model prior covariance matrix as (Reference Ortega CulaciatiOrtega Culaciati, 2013)

where κ is a scalar weighting parameter, whose value can be chosen to reduce high-frequency variations in the posterior model, and

Because we assume m ₀ to be smoothly varying, it follows that . We note that applying Eqn (18) to Eqn (8) reduces Eqn (8) to a form similar to Tikhonov-regularized least-squares.

Applying this form of Ω_d results in a spatially varying damping factor and weights the elements of the Laplacian in terms of the interferometric phase variances, through , and the viewing geometry, via G. A key advantage in using the Laplacian form of is that we do not impose specific values for our prior model, m ₀, as would be the case if we assumed . In other words, we implicitly apply a correlation length to m ₀ that effectively makes velocity models with unphysical discrete jumps less likely in areas with non-ideal data coverage.

Limitations of the Bayesian method

The primary limitations of the Bayesian method discussed above are attributable to limitations of the chosen form of the prior model, m ₀, its associated prior covariance, C _m, and incomplete representations of errors described by the data covariance matrix, C _d. These limitations are consistent with any inverse problem, but we discuss them here in the context of our formulation of . For this study, we chose a posterior model that contains only an average velocity vector field, because the temporal sampling of our dataset negates the possibility of resolving other model components. This model implicitly assumes that there is no acceleration. However, our data were collected over a finite time during the early melt season, when diurnal variations in velocity driven by variable surface meltwater flux might be present. As a result, a time-invariant average velocity does not perfectly represent glacier motion and this shortcoming in the model introduces a prediction error (Reference Duputel, Agram, Simons, Minson and BeckDuputel and others, 2014) for which we currently cannot account. Though beyond the scope of this study, we note that, in general, a model could be composed of a mean velocity, a secular acceleration, a series of harmonic functions and any number of transient functions (Reference Hetland, Musé, Simons, Lin, Agram and DiCaprioHetland and others, 2012; Reference Riel, Simons, Agram and ZhanRiel and others, 2014).

The most appropriate form of the prior data covariance matrix must be based on a number of factors. As previously discussed, we assume that the data covariance matrix used in this study is diagonal and thereby neglect off-diagonal components. Our main motivations for this assumption are simplicity and computational tractability. We recognize that the resulting posterior model covariance matrix does not fully capture the total, or true, errors, but we expect it to capture most of the formal errors because the off-diagonal components of C _d are small relative to the diagonal components, for the reasons discussed above. It is important to note that formal errors are not necessarily a complete or accurate representation of total errors. In the model defined in Eqn (8), formal errors are described by the posterior model covariance matrix, a function of viewing geometry and InSAR correlation only. Any noise sources that do not impact either of these parameters are unaccounted for in the formal error. The classic example of such a noise source is tropospheric delay (e.g. Reference Lohman and SimonsLohman and Simons, 2005). Phase shifts caused by differences in the volumetric moisture content at or near the surface of the glacier between radar acquisitions will also introduce errors that are not fully manifest in our current error estimation.

Moisture-induced error

Just as propagation through water vapor in the troposphere can cause erroneous deformation signals in interferograms (e.g. Reference Zebker, Rosen and HensleyZebker and others, 1997; Reference HanssenHanssen, 2001; Reference Emardson, Simons and WebbEmardson and others, 2003; Reference Lohman and SimonsLohman and Simons, 2005), moisture contained in a volume through which the signal propagates can cause phase offsets. There is even evidence to suggest changes in the moisture content of natural media that are typically thought of as purely surface scatterers, such as bare agricultural fields, can cause phase shifts (Reference Nolan and FatlandNolan and Fatland, 2003; Reference Nolan, Fatland and HinzmanNolan and others, 2003; Reference Khankhoje, Van Zyl and CwikKhankhoje and others, 2012). Our interest is in accurate estimates of surface velocity during the early melt season, so it is important to consider the potential influence of moisture content on the radar signal in order to properly understand how InSAR phase estimates relate to actual glacier motion. Hereafter, when discussing moisture content or surface moisture in the context of InSAR measurements and as they relate to moisture-induced errors, we are referring to the volumetric concentration of liquid water in the uppermost region of the glacier that extends down to at least the penetration depth of the radar signal.

Whenever air temperatures exceed the melting temperature, snow at or near a glacier’s surface will be infused with liquid water and exhibit variations in moisture content over timescales shorter than repeat-pass time intervals. Changes in moisture content in the near-surface influence SAR phase values via the permittivity of the media. The real component of the permittivity of liquid water is ∼80, which differs markedly from the real permittivity for dry snow, which is ∼1 (Reference Ulaby, Moore and FungUlaby and others, 1986). As a result, even small changes in surface moisture content can significantly influence the electromagnetic properties of media observed with shallow-penetrating radar.

We use a simple empirical model to show the dependence of permittivity on the moisture content. Based on laboratory measurements of the scattered electric fields for snow samples with various moisture contents, Reference Hallikainen, Ulaby and AbdelrazikHallikainen and others (1986) modeled the relative permittivity of wet snow, ε, as

where f ₀ is the radar frequency in free-space, f _r ≈ 9 GHz is the relaxation frequency of liquid water in snow, c_ρ is a constant (1.83 × 10⁻³ m³ kg⁻¹), ρ _s is the density of dry snow and ν _m is the dimensionless volumetric concentration of liquid water. Eqns (20–23) show that, while both the real and imaginary permittivity components increase with the liquid water content, the real component is more sensitive to changes in moisture content over the range of standard SAR frequencies (Fig. 2a and b). Increasing permittivity can: (1) reduce the penetration depth, δ _p, of the radar signal , where k ₀ is the radar wavenumber in free space; Fig. 2c); (2) increase the wavenumber of the penetrating radar signal ; and (3) increase the power scattered by the free surface due to increased surface reflectivity, which scales with permittivity and local incidence angle. (Reference Ulaby, Moore and FungUlaby and others (1986) offer a more thorough treatment of the electromagnetic properties of various media, as it pertains to radar remote sensing.)

Fig. 2.

(a) Real and (b) imaginary components of the permittivity of snow as a function of liquid water content, ν _m, calculated from Eqns (20–23) (solid curves) using ρ _s = 500 kg mm⁻³ (Reference Hallikainen, Ulaby and AbdelrazikHallikainen and others, 1986), f ₀ = 1.25 GHz for L-band and f ₀ = 9.65 GHz for X-band. Dashed curves show the square root of the like-colored solid curve. (c) The penetration depth, δ _p, of a homogeneous medium with a constant permittivity profile dictated by the liquid-water content.

An idealized, heuristic model of InSAR phase shows the effect of changing permittivity between two SAR scenes. This conceptual model is concerned only with changes in moisture that occur over short (hourly-to-daily) timescales, because InSAR is generally ineffective at longer timescales, particularly in areas that experience surface melt. We adopt a basic two-layer model, one thin layer overlying a half-space, to describe moisture changes. Over the timescales of interest, we assume ν _m is a function of only depth z and varies only in the uppermost layer extending to depth h; ν _m is constant for z > h. If h ≪ δ _p, δp will be approximately constant between SAR acquisitions. If we neglect the surface scattered component of the received radar signal, a reasonable approximation when ν _m is small (Reference MätzlerMätzler, 1998; Reference Oveisgharan and ZebkerOveisgharan and Zebker, 2007), we can write a simple model for the InSAR phase, φ, of a unit-amplitude incident electric field scattered from a stationary volume as

where subscripts a and b indicate the two SAR scenes used to generate the interferogram, indicates the real component of the complex argument and , where θ_i is the local radar incidence angle. By comparing Eqns (20–23) and Eqn (24) and the permittivity values shown in Figure 2a, we see that even small changes in ν _m for h > 0 will influence the InSAR phase. (More in-depth descriptions of the salient physics are given by Reference IshimaruIshimaru (1978) and Reference Ulaby, Moore and FungUlaby and others (1986).)

UAVSAR acquisitions

The data used in this study were collected as part of a UAVSAR campaign in which we collected data for six non-continuous days over the course of 12 days. At the time of data acquisition, UAVSAR was being flown aboard a NASA Gulfstream III aircraft that cruises at ∼12.5 km altitude, providing an incidence angle range of 22–65°, which we trim to 40–65°. Data were collected along 15 unique flight lines that were designed to image Langjökull and Hofsjökull (Fig. 1) from at least three different LOS vectors during each data collection. The flights were scheduled such that the first 3 days of data were collected in the afternoon, a few hours after the expected maximum daily melt based on temperature, and the final 3 days of data were collected in the early morning, ∼12 hours offset from the afternoon data collection. The viewing geometries of the flight lines provide good constraints on the ice flow over most areas (Fig. 3).

Fig. 3.

Λ_g and (inset) number of SAR scenes collected each day by UAVSAR in June 2012 for (a) Langjökull and (b) Hofsjökull. Black contours denote UAVSAR scene boundaries.

UAVSAR is a fully polarimetric, L-band (1.25 GHz) SAR system, whose integrated autopilot system, inertial navigation unit (INU) and real-time GPS system are capable of piloting the aircraft through a 10 m diameter tube that encases the proposed flight line. This aerial precision facilitates repeat-pass interferometric observations whose temporal baseline and LOS vectors can be programmed. High bandwidth (80 MHz) and a large along-track antenna length give UAVSAR a raw spatial resolution of 1.9 m in range (cross-track) and 0.8 m in azimuth (along-track) (Reference Hensley, Zebker, Jones, Michel, Muellerschoen and ChapmanHensley and others, 2009a).

Random aircraft motions complicate the InSAR processing task for UAVSAR data relative to data acquired by satellite-based SAR systems. During processing, these random motions are largely accounted for using data from the INU and GPS system. Centimeter-scale motion between aircraft repeat passes (i.e. the residual interferometric baseline) that is not accounted for using the INU and GPS data is estimated by calculating the amplitude cross-correlation of the two scenes and considering only the range-correlated signals to be artifacts of aircraft motion (Reference HensleyHensley and others, 2009b). Small residuals can remain after this process, but because UAVSAR maintains very small perpendicular baselines (typically <2 m), the baseline correlation component, γ _b, is ∼1. Small perpendicular baselines and UAVSAR’s high SNR help ensure that the majority of decorrelation in the repeat-pass interferograms is due to temporal variations in the scatterers and volumetric decorrelation.

InSAR post-processing

After most of the random motion components were removed from the phase, along with estimates of Earth’s curvature and local topography provided by a DEM, we further processed the data to retrieve interferograms that are useful for inferring the velocity field. We employ a custom DEM that combines data from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER; version 1) with DEMs for Langjökull, derived from a GPS survey conducted in 1997 (Reference Palmer, Shepherd, Björnsson and PálssonPalmer and others, 2009), and Hofsjökull, derived from lidar surveys conducted in 2008 and 2010 (Reference JóhannessonJóhannesson and others, 2013). Post-processing includes:

Hereafter, we refer to the flattened, unwrapped, geolocated interferograms simply as interferograms.

We calculated the LOS unit vectors for each interferogram from the geometry of the flight track and the platform position. We reference the LOS vectors to a local east–north–up coordinate system, whose origin is coincident with the pixel location in the geolocated image and referenced to the WGS84 spheroid. The modeled velocity field is referenced to the same coordinate system.

GPS collection and processing

We deployed ten GPS stations on Vestari-Hagafellsjökull, an outlet glacier chosen for logistical simplicity and safety of the field crew. Data were continuously collected every 15 s for ∼2 weeks. The GPS collection window began 2 days prior to the first UAVSAR acquisition. Eight of the GPS stations operated throughout the 2 week deployment period while the other two stations lost power and stopped collecting data after ∼1 week. All GPS receivers were mounted on poles sunk several meters into the ice so that the GPS data capture the kinematics of the underlying ice, not the free surface.

We processed the raw GPS data using kinematic Precise Point Positioning (PPP) methods available as part of GNSS Inferred Positioning System and Orbit Analysis Simulation Software (GIPSY-OASIS) (e.g. Reference BertigerBertiger and others, 2010). PPP eliminates the need for a ground reference station by using precise clocks along with predetermined satellite orbits (Reference Zumberge, Heflin, Jefferson, Watkins and WebbZumberge and others, 1997), and kinematic processing allows for higher-frequency position updates during processing, which is most suitable for continuously deforming areas. We smoothed the processed positions over a 6 hour window and referenced all motions to the same reference frame as the InSAR-derived velocity fields.

Inferred velocity field

We present examples of the inferred horizontal component of ice flow for Langjökull and Hofsjökull, along with estimates of the viscous component of ice flow (Fig. 4). Arrows in Figure 4a and c indicate the direction of horizontal flow, and the color map represents the horizontal speed, which we smoothed using a 200 m × 200 m low-pass filter. The InSAR data were collected in the early mornings of 13 and 14 June 2012 (Δt ≈ 24 hours), ∼14 hours after the maximum daily melt. Unusually warm weather over much of Iceland in early June 2012 caused the atmospheric temperature over central Iceland to peak above freezing in the late afternoon for several days prior to and during data collection.

Fig. 4.

(a) Horizontal velocity field for Hofsjökull, inferred from InSAR data collected on 13 and 14 June 2012 (Δt ≈ 24 hours). Arrows indicate the direction of the ice flow, and the color map indicates the horizontal speed. (b) Velocity estimated from a simple viscous flow model that does not account for basal slip (Eqn (25)). The color map indicates the speed of viscous flow, and arrows indicate the ice surface gradient. The difference between the estimated viscous flow speed and the measured speed in the outlet glaciers is indicative of slip at the glacier bed. Black contours are the same as in Figure 1, and tan-colored areas surrounding the glacier show ground elevation. (c, d) Same as (a, b), but for Langjökull.

Juxtaposed with the horizontal velocity fields are velocity fields estimated from a simple ice-flow model that accounts for only the internal viscous deformation in the ice, neglecting all sliding along the glacier bed (Fig. 4b and d). Assuming surface-parallel flow and a linear depth-dependent driving stress profile, it can be shown that the viscous component of the ice flow has a power-law relation to ice thickness and surface slope and is given as (Reference Cuffey and PatersonCuffey and Paterson, 2010)

where τ _b is basal shear stress and H is ice thickness. We assume τ _b = ρ _ice gHα, the gravitational driving stress, where α is the ice surface slope and ρ _ice = 900 kg m⁻³ is the depth-averaged density of glacier ice. Over broad areas, basal stress cannot exceed gravitational driving stress, meaning that our results for Eqn (25) approximate the maximum viscous deformation rate. The variables n and A arise from Glen’s flow law, the nonlinear constitutive relationship between the effective strain rate, , and the effective stress, τ _E, within the ice , and are taken to be 3 and 2.4 × 10⁻²⁴ Pa⁻³ s⁻¹, respectively (Reference Cuffey and PatersonCuffey and Paterson, 2010). We averaged slope and thickness over a window that is ∼10 times the average ice thickness in all directions. Arrows in Figure 4b and d, which are co-located with arrows in Figure 4a and c, respectively, indicate the free surface gradient, and the color map represents the speed calculated from Eqn (25). Because the surface slopes are small (<20°) and the viscous flow model is a simple model, we do not convert the surface-parallel values from Eqn (25) to true horizontal motion.

InSAR and GPS results

GPS data collected on Vestari-Hagafellsjökull allow us to validate a portion of the InSAR-derived velocity field (Fig. 5). Horizontal speeds calculated from GPS data represent the mean velocity over the same time window as the UAVSAR data. Black, white and red arrows in Figure 5a indicate the GPS-measured flow vector, the mean surface gradient vector and the InSAR-derived flow vector, respectively, and are not drawn to scale. Co-located GPS and InSAR-derived horizontal speeds are shown in Figure 5b, where the solid gray curve is the one-to-one regression line and the vertical error bars are derived from Λ_m. Horizontal speeds along a transect that runs from the black X in Figure 5a through GPS stations L04–L01 are shown in Figure 5c, along with Λ_m and Λ_g for the same transect.

Fig. 5.

(a) GPS receiver locations, flow direction (black arrows) and horizontal speed (colored circles) overlaying the contemporaneous InSAR-derived horizontal speed (colored surface) and co-located flow direction (red arrows) calculated using κ = 10⁻⁵ (Eqn (18)) – the same velocity field as in Figure 4c. Circles are colored on the same color scale as the InSAR-derived field. White arrows indicate the horizontal component of the co-located mean ice surface gradient. Arrow lengths are not to scale. Inset map indicates the extent and location of the main map. Black contours are the same as in Figure 1. (b) InSAR versus GPS-derived horizontal speed. InSAR values are taken from the velocity field shown in (a). Vertical error bars are derived from Λ_m. GPS errors are too small to represent on this scale. (c) InSAR (blue curve) and GPS-derived horizontal speeds (blue circles), Λ_m (solid gray curve) and Λ_g (dashed gray curve) along the transect that begins at the black X in (a) and follows GPS stations L04–L01.

Posterior model covariance

East and north variances from provide estimates of the reliability of the respective inferred velocity component (Fig. 6a and b for Hofsjökull and Fig. 6d and e for Langjökull), while Λ_m indicates the total error in the posterior model (Fig. 6c and f). Hofsjökull and, to a lesser extent, Langjökull have low expected errors over most of their surface. High variance values occur in areas where the ice motion is poorly constrained by InSAR observations (Fig. 3). High-frequency features in Figure 6 are attributable to variations in interferometric correlation.

Fig. 6.

East (left column) and north (center column) variance and Λ_m (right column) for the InSAR-derived velocity fields in Figure 4a and c for Hofsjökull (a–c) and Langjökull (d–f). Black contours are the same as in Figure 1.

Moisture-induced error

The interferograms used to derive the velocity fields contain phase offsets that we attribute to differential volumetric moisture content in the glacier near-surface. Examples of these offsets are shown, along with estimates of ambient air temperature, in Figure 7. We present interferograms collected from two different flight lines that image Hofsjökull on different days (Fig. 7b–e), along with double-differenced interferograms, i.e. the differences of the respective interferograms (Fig. 7f–h). Because Hofsjökull is roughly dome-shaped and we designed the flight lines to look up the surface slope as much as possible, these lines are representative of all flight lines. Data used to derive the velocity field in Figure 4a are shown in Figure 7e. Interferograms and double-differenced interferograms corresponding to data collected in the afternoon (Fig. 7b–c and f), when surface moisture content should be highest, show a more distinct phase offset relative to data collected in the morning at higher elevations. The differential phase sign change in Figure 7f occurs at an elevation where temperatures are ∼0°C on 5 June.

Fig. 7.

Evidence of moisture-induced phase offsets. (a) Estimates of air temperature at elevations indicated in the legend and by like-colored contours in (b, h). Vertical black and gray dashed lines indicate the time of the first and last UAVSAR acquisitions, respectively. (b–e) Examples of flattened, unwrapped interferograms of Hofsjökull for various daily pairs. The respective LOS directions are shown by arrows in (b), and the parenthetical ‘am’ and ‘pm’ indicate whether the data were collected in the morning or afternoon. Note that color maps are saturated to elucidate phase variations in regions of slow-moving ice. (f–h) Double-differenced interferograms for the two flight lines shown in (b–e) and an additional flight line covering the area in between. LOS direction in the central flight line is the same as the southern flight line. The striped pattern evident outside the glacier boundary in the northern flight line in (f, g) is residual aircraft motion, and the region southwest of Hofsjökull was covered with snow throughout the UAVSAR campaign. Contour elevations are the same as in Figure 1.

Differential surface moisture content is the only plausible explanation for the residual phase offsets shown in Figure 7. All InSAR data used in this study have perpendicular baselines, B _⊥ < 10 m, and most data have B _⊥ < 5 m. Because topographic sensitivity scales with B _⊥, the InSAR data shown here have virtually no sensitivity to DEM errors. Furthermore, the residual phase increases from near-zero over the bare ground to >1 rad over ∼10 m near the edge of the ice. The troposphere is not capable of supporting moisture gradients steep enough to account for such steep phase gradients. Instead, atmospheric phase offsets should smoothly vary from the ice to bare ground over length scales that are at least an order of magnitude larger than observed in Figure 7.