Introduction
Glaciers transport ice from areas of mass accumulation at high elevation to areas of mass loss at lower elevations through a combination of internal deformation and slip at the ice/bed interface. While the constitutive relation for viscous flow of glacier ice is relatively well known (Reference GlenGlen, 1955; Reference NyeNye, 1957; Reference MacAyealMacAyeal, 1989), understanding the mechanics of basal slip, which includes the sliding of ice relative to a stationary bed and deformation of the bed, remains an open problem (e.g. Reference Howat, Tulaczyk, Waddington and BjörnssonHowat and others, 2008; Reference SchoofSchoof, 2010; Reference Bartholomaus, Anderson and AndersonBartholomaus and others, 2011; Reference HewittHewitt, 2013; Reference Werder, Hewitt, Schoof and FlowersWerder and others, 2013). Slip at the glacier bed is an important component of velocity fields of many glaciers (e.g. Reference BoultonBoulton, 1979; Reference Engelhardt and KambEngelhardt and Kamb, 1998; Reference Tulaczyk, Kamb and EngelhardtTulaczyk and others, 2000; Reference Kamb, Alley and BindschadlerKamb, 2001), accounting for observed seasonal and diurnal velocity variations (e.g. Reference Rignot and KanagaratnamRignot and Kanagaratnam, 2006; Reference Shepherd, Hubbard, Nienow, McMillan and JoughinShepherd and others, 2009; Reference JoughinJoughin and others, 2012), and is imperative for erosion to occur (Reference BoultonBoulton, 1979; Reference HalletHallet, 1996; Reference IversonIverson, 2012). As a result, understanding subglacial mechanics is crucial for developing predictive models of the future states of glaciers, estimating the contribution of glacier melt to sea-level rise, and improving our knowledge of how glaciers shape the landscape.
Direct observations of glacier beds are often impractical, whereas surface velocity fields are relatively easy to observe and are useful for inferring subglacial mechanical and hydrological properties (e.g. Reference Iken and BindschadlerIken and Bindschadler, 1986; Reference KambKamb, 1987; Reference Tulaczyk, Kamb and EngelhardtTulaczyk and others, 2000; Reference Zwally, Abdalati, Herring, Larson, Saba and SteffenZwally and others, 2002; Reference Magnússon, Rott, Björnsson and PálssonMagnússon and others, 2007, Reference Magnússon, Björnsson, Rott and Pálsson2010, Reference Magnússon2011). Repeat-pass interferometric synthetic aperture radar (InSAR) can provide synoptic-scale observations of glacier surface velocities and has been used to map the velocity fields of many glaciers (e.g. Reference Joughin, Fahnestock, MacAyeal, Bamber and GogineniJoughin and others, 2001; Reference Rignot, Mouginot and ScheuchlRignot and others, 2011), often during winter. To gain a better understanding of the subglacial mechanics, and their interdependence with basal hydrology, it is desirable to observe glaciers during the melt season, in particular the early melt season (Reference SchoofSchoof, 2010), when surface meltwater flux can induce variations in basal slip on hourly-to-monthly timescales (e.g. Reference Zwally, Abdalati, Herring, Larson, Saba and SteffenZwally and others, 2002; Reference Shepherd, Hubbard, Nienow, McMillan and JoughinShepherd and others, 2009; Reference JoughinJoughin and others, 2012). However, the amplitudes of velocity fluctuations can be small relative to the mean background velocity (Reference Bartholomew, Nienow, Mair, Hubbard, King and SoleBartholomew and others, 2010), necessitating accurate and robust InSAR analysis techniques.
It is difficult to make useful InSAR measurements of glaciers during the early melt season because the surface often changes rapidly between SAR acquisitions, inducing high noise levels in InSAR data. Because most InSAR data are acquired by spaceborne systems, the time between repeated acquisitions is fixed and typically on the order of days to weeks, depending on the radar system. To overcome these limitations, we acquired repeat-pass InSAR data using NASA’s uninhabited aerial vehicle synthetic aperture radar (UAVSAR), an airborne, L-band (24 cm wavelength) SAR system that allows us to choose the repeat-pass time interval and to design custom flight lines (e.g. Reference Hensley, Zebker, Jones, Michel, Muellerschoen and ChapmanHensley and others, 2009a). In June 2012, we collected InSAR data over the ice caps Langjökull and Hofsjökull, central Iceland (Fig. 1), from multiple vantage points for six non-continuous days over a 12 day period. We designed the flight lines to provide complete spatial coverage of both ice caps from at least three different look directions with 24 hour repeat-pass times. Short repeat-pass times and L-band radar provided data with acceptable signal-to-noise ratios (SNR) everywhere on both ice caps.
Fig. 1. (a) Shaded relief map of central Iceland showing Langjökull and Hofsjökull. Inset map shows the location of the region of interest. (b, c) Horizontal transmit, horizontal receive (HH) amplitude images (expressed as normalized radar cross section, σ o) with ice divides (cyan lines) and major outlet glaciers of Langjökull (b) and Hofsjökull (c) (Reference BjörnssonBjörnsson, 1988). Contours indicate surface elevation in 150 m increments. The major outlet glaciers of Langjökull are: Þjófadalajökull (LJ), Leiðarjökull (LL), Kirkjujökull (LK), Norðurjökull (LN), Suðurjökull (LS), Eystri-Hagafellsjökull (LE), Vestari-Hagafellsjökull (LV), Svartárjökull (LT), Flosaskarðsjöklar (LF), Þrístapajökull (LÞ) and Baldjökull (LB). The major outlet glaciers of Hofsjökull are: Illviðrajökull (HI), Þjórsárjökull (HÞ), Múlajökull (HM), Blautukvislarjökull (HT), Blágnípujökull (HB), Blöndujökull (HL), Kvíslajökull (HK) and Sátujökull (HS).
Langjökull and Hofsjökull are natural laboratories that we can use to investigate basal mechanics more easily and in much greater detail than is practical in many other regions. Both ice caps are land-terminating and cover areas of ∼900 km2 with typical ice thicknesses >200 m (Reference Björnsson and PálssonBjörnsson and Pálsson, 2008). The beds of both ice caps have been near-completely mapped with ice-penetrating radar (Reference BjörnssonBjörnsson, 1986), and recently obtained surface digital elevation models (DEMs) are available, mostly from lidar surveys. Previous studies of the bedrock lithology show that porous lavas underlie southern Langjökull, whereas the remainder of Langjökull and all of Hofsjökull rest on impermeable bedrock (Reference Björnsson, Pálsson, Sigurðsson and FlowersBjörnsson and others, 2003). Both ice caps have been previously studied using InSAR data acquired with 3 day repeat-pass time by European Remote-sensing Satellite (ERS) in 1994 (Reference Palmer, Shepherd, Björnsson and PálssonPalmer and others, 2009; Reference Gourmelen, Kim, Shepherd, Park, Sundal and BjörnssonGourmelen and others, 2011).
In this study, we present a Bayesian approach to inferring three-dimensional (3-D) velocity fields from multiple InSAR acquisitions. Our approach incorporates a data correlation length to reduce large offsets in the velocity field, which often occur at InSAR scene boundaries. We use this Bayesian method to infer the horizontal velocity fields of Langjökull and Hofsjökull from UAVSAR data. We also present evidence that differential surface moisture content causes phase offsets that corrupt estimates of the vertical velocity component but do not cause significant errors in the inferred horizontal velocity components. We compare the inferred velocity field with a simple viscous flow model and co-located GPS data, discuss the observed characteristics of the outlet glaciers, and show that basal slip is likely to account for more than half of the observed surface velocities in many outlet glaciers.
Methodology
InSAR encompasses commonly used methods for measuring deformation or topography of an area that has been imaged at least twice by a SAR system. SAR data are complex-valued, providing information on both the amplitude and phase of the radar waves. InSAR processing takes two SAR scenes and aligns them such that sub-wavelength changes in the path distance between the radar antenna and a given target can be calculated using the difference in the phase of each image (e.g. Reference RosenRosen and others, 2000). In this study we consider repeat-pass InSAR, which requires two SAR images acquired over a given area at different times and from approximately the same sensor position. The final products, called interferograms, provide measurements of displacement along the radar line-of-sight (LOS) unit vector, 
, pointing from the sensor position, r′
i
, to a ground position, r. The InSAR phase per unit time in interferogram i is given as
where u i (r) is the displacement vector at r over repeat-pass time interval Δti .
We can write the InSAR measurements in matrix form as
where 
is the mean velocity vector and G(r, r′) is the design matrix. Rows in G consist of LOS unit vectors associated with the entries of the InSAR phase vector, d(r), which we take to be of the form in Eqn (1). If G contains three or more LOS vectors that are sufficiently different from one another, we can solve Eqn (2) for an estimate of the mean velocity vector, 
.
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 0 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 Ni , 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 sz 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 ∇2 m 0 ≈ 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
0 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
0, as would be the case if we assumed 
. In other words, we implicitly apply a correlation length to m
0 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
0, 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
0 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−3 m3 kg−1), ρ
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
0 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−3 (Reference Hallikainen, Ulaby and AbdelrazikHallikainen and others, 1986), f 0 = 1.25 GHz for L-band and f 0 = 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).)
Data
In June 2012, we collected GPS and InSAR data over Langjökull and Hofsjökull. NASA’s UAVSAR system collected InSAR data for 6 days beginning 3 June. Ten GPS stations deployed on Vestari-Hagafellsjökull collected 15 s dual-frequency data during the InSAR campaign. Two of these GPS stations, at elevations relative to mean sea level of ∼490 and 1100 m, were co-located with automatic weather stations that recorded air temperatures throughout the campaign.
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:
Averaging and decimating each interferogram using a 3 pixels × 12 pixels (range × azimuth) averaging window (a process commonly called ‘looking’), which yields 5 m × 7.2 m pixels in radar coordinates;
Filtering the interferogram using a 10 × 10 pixel equiweighted moving average, or low-pass, filter;
Unwrapping the filtered interferogram using the statistical cost, network flow algorithm for phase unwrapping (SNAPHU) (Reference Chen and ZebkerChen and Zebker, 2000, Reference Chen and Zebker2001, Reference Chen and Zebker2002). During unwrapping, we do not distinguish between pixels that image the glaciers and those that image the bare ground surrounding the glaciers, so InSAR phase is continuous when transitioning between rock and ice;
Geolocating the unwrapped interferograms;
Flattening the geolocated interferograms by fitting and removing a phase surface from the bare ground surrounding the glaciers. We assume the ground around the glaciers was stationary during the 24 hour period between data acquisitions and that this flattening process removes the majority of tropospheric delay (Reference Zebker, Rosen and HensleyZebker and others, 1997).
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.
Results
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−3 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−24 Pa−3 s−1, 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−5 (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.
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.
Discussion
Our derived horizontal velocity fields on both ice caps qualitatively agree with previously published results that used ERS data collected in February 1994. Over Hofsjökull, Reference Gourmelen, Kim, Shepherd, Park, Sundal and BjörnssonGourmelen and others (2011) show the same outlet glaciers and general flow pattern, but flow velocities on Illviðrajökull, Þjórsárjökull and Blautukvislarjökull are markedly higher in February 1994 relative to our results. Because the earlier data were collected in winter and show higher velocities, it is likely that these three glaciers were surging in February 1994. Þjórsárjökull is known to have surged in 1994 (Reference Björnsson, Pálsson, Sigurðsson and FlowersBjörnsson and others, 2003), but these observations mark the first recorded surges on Illviðrajökull and Blautukvislarjökull. Our results on Langjökull generally agree with Reference Gourmelen, Kim, Shepherd, Park, Sundal and BjörnssonGourmelen and others (2011) and Reference Palmer, Shepherd, Björnsson and PálssonPalmer and others (2009). Lower flow velocities on Suðurjökull in our study relative to 1994 lend credence to the postulate by Reference Palmer, Shepherd, Björnsson and PálssonPalmer and others (2009) that this glacier was surging in 1994. Unlike Reference Palmer, Shepherd, Björnsson and PálssonPalmer and others (2009), we do not observe elevated velocities at the edges of the ice cap.
Spatial patterns in the horizontal velocity fields of Hofsjökull and Langjökull are broadly consistent with the models of viscous ice flow, and the directionality of the measured ice flow is generally along the ice surface slope on both ice caps. Areas of low velocity not located near the edge of the glacier indicate ice divides, and areas of high velocity correspond to known outlet glaciers (Fig. 1; Reference BjörnssonBjörnsson, 1988; Reference Björnsson and PálssonBjörnsson and Pálsson, 2008; Reference Gourmelen, Kim, Shepherd, Park, Sundal and BjörnssonGourmelen and others, 2011). Significant differences between the magnitude of estimated viscous flow and measured velocity fields indicate that basal slip is likely to be an important component of the ice flow. Basal slip rates may be greater on Hofsjökull, which typically has higher driving stresses than Langjökull and comparable mechanical properties at the bed (Reference BjörnssonBjörnsson, 1986, Reference Björnsson1988; Reference Jóhannesson and SæmundssonJóhannesson and Sæmundsson, 1998). It remains an open question for future work to ascertain the extent to which surface meltwater flux may be influencing basal slip on each glacier.
On Vestari-Hagafellsjökull, southwest Langjökull, the InSAR-derived velocities agree with co-located GPS velocity measurements. Expected errors in the InSAR-derived velocity field near the upstream GPS stations, L04 and M02, are relatively low and the errors increase by more than a factor of five downstream. Stations L03, L04 and M03, whose locations have small Λm values, are in close agreement with the InSAR-derived horizontal velocities. High noise in the east component of the InSAR-derived velocity field causes a significant eastward shift in the inferred velocity vector at M05 but not the velocity magnitude, which matches the GPS-measured speed. GPS-derived horizontal speed is higher at station M02 than the InSAR-derived velocity field. Because M02 is located in an area with a steep velocity gradient, spatial filtering of the InSAR data is likely to account for the discrepancy between GPS and InSAR speeds. The InSAR-derived velocity near station L01 should be considered unreliable, due to high noise caused by suboptimal viewing geometries and low correlation. During fieldwork, we observed significant surface lowering in the vicinity of L02 and postulate that surface dynamics are the source of the disparate InSAR and GPS velocities. We secured the GPS antenna poles several meters into the ice, making GPS measurements largely independent of changes in the local free-surface height. By contrast, radar waves incident on the area around L02, where the surface was either wet snow or exposed ice throughout the data collection window, are likely to be influenced by free-surface dynamics because of high radar reflectivity at the surface. The generally good agreement between InSARderived and GPS velocities at most station locations suggests that the InSAR-derived velocity fields of both ice caps are reliable enough to allow for analysis of the ice flow.
A variety of characteristic features are evident in the outlet glaciers on both ice caps. Approximately one-third of the outlet glaciers, Múlajökull (HM), Blautukvislarjökull (HT), Blágnípujökull (HB), Kvíslajökull (HK), Vestari-Hagafellsjökull (LV), Þrístapajökull (LÞ) and Suðurjökull (LS), form concentrated ice streams that have distinct regions of considerably higher velocity than the surrounding ice. Most of these glaciers occur in areas with high driving stress, though the viscous flow model is a poor indicator of the location and magnitude of the fastest-moving ice, further supporting the idea that basal slip is an important contributor to the total ice flux. The remaining outlet glaciers tend to transport ice more diffusively across their widths, having lower shear rates at their margins relative to the more stream-like glaciers. The general characteristics of these sheet-like outlet glaciers are represented, for the most part, in the simple viscous flow models, though their inferred horizontal speeds can be more than a factor of two larger than the model predicts.
The fastest glaciers on Hofsjökull are Múlajökull, Blágnípujökull and Kvíslajökull, with Múlajökull having the highest transport of the three. These three fast glaciers, along with Sátujökull (HS) in the north, Þjórsárjökull in the east and Blöndujökull in the west, are known to surge (Reference Björnsson, Pálsson, Sigurðsson and FlowersBjörnsson and others, 2003). It is worth noting that only the high-velocity area of Kvíslajökull correlates with an expected high in viscous ice flow, based on the simple viscous model, indicating that a large amount of basal slip may be present beneath these glaciers.
Fast glaciers on Langjökull include Þrístapajökull, Norðurjökull, and Vestari-Hagafellsjökull, the deployment site for the GPS stations used in this study. The highest-velocity areas in these fast glaciers generally correlate with areas that have higher predicted rates of viscous flow, though the measured velocity is a factor of 1.5–2 times higher than the viscous deformation component, likely indicating significant basal slip in these areas. Vestari-Hagafellsjökull and its eastern neighbor, Eystri-Hagafellsjökull, are known surge-type glaciers (Reference Björnsson, Pálsson, Sigurðsson and FlowersBjörnsson and others, 2003).
We omit results for the vertical component of the InSAR-derived velocity field from this study because moisture-induced phase offsets make the true vertical component of the ice velocity inaccessible. For example, inferred vertical velocity components corresponding to the horizontal velocity fields given in Figure 4 have a median value of ∼20 m a−1 in the up direction, with higher values manifested in the southeast quadrant, which receives relatively high levels of solar radiation during the early melt season, and few areas that indicate downslope flow. None of the glaciers that are apparent in the horizontal field are manifested in the vertical velocity component. These results are not physically justifiable. We note that phase offsets occur along the radar LOS and are approximately constant on given days, meaning that the vertical component of the velocity field is particularly sensitive to moisture-induced phase offsets. Future work should focus on methods to decouple the moisture-induced phase offsets from ice motion, to allow for accurate estimates of the surface-normal component of the velocity field, an important quantity for ice-flow modeling and studies of basal slip.
Moisture-induced phase offsets in InSAR data have the potential to pose problems for high-precision InSAR applications in areas where liquid water is present at the glacier surface. These areas include southern Greenland, the Canadian Arctic Archipelago, the Antarctic Peninsula and peripheral zones, and many mountain glaciers. Here we have shown examples of moisture-induced phase offsets in repeat-pass interferograms, but we note that the dependence of phase on the radar incidence angle (Eqn (24)) suggests that InSAR-derived DEMs, which utilize SAR data collected from different antenna positions at approximately the same time, may also contain biases that can be linked to volumetric surface moisture content. As a result, InSARbased studies of glacial mass balance as well as kinematics will benefit from careful scrutiny of the data and consideration of the possible presence of moisture-induced phase offsets.
Conclusion
We have presented a new Bayesian approach for inferring 3-D velocity fields from multiple InSAR acquisitions, and use the new method to infer the horizontal velocity fields for two temperate ice caps, Langjökull and Hofsjökull, from airborne, L-band, repeat-pass InSAR scenes collected in June 2012. The flow directions of the horizontal fields on both ice caps closely agree with the free surface slope, and patterns in the horizontal speed are broadly consistent with estimates of viscous flow calculated from surface and bedrock DEMs using a simple model. InSAR-derived horizontal velocities correspond to co-located GPS velocities for the same time period on Vestari-Hagafellsjökull, except in areas where the ice motion is not well constrained by the InSAR data.
The InSAR-derived horizontal speeds differ markedly from the predicted viscous flow speeds, likely indicating the importance of basal slip to total ice flow on Langjökull and Hofsjökull. Both ice caps contain numerous outlet glaciers with various flow characteristics. These characteristics range from stream-like outlet glaciers, which flow significantly faster than the surrounding ice, to sheet-like glaciers, which transport ice over broad regions with little internal strain and low strain rates at the margins. Due to the variety of outlet glaciers, similar bed properties, and consistent climate forcing, these ice caps offer a valuable natural laboratory with which to study the mechanics of basal slip.
Differential surface moisture content on the glaciers prevented reliable estimates of the vertical component of the velocity fields. Because moisture-induced phase offsets are approximately constant across all interferograms and are small relative to the phase offsets over the flowing ice, the error is manifested primarily in the vertical velocity field and has little influence on horizontal velocities.
Acknowledgements
Funding for the UAVSAR campaign was provided by NASA’s Cryospheric Sciences Program. B.M. was supported by a NASA Earth and Space Sciences fellowship and an ARCS Foundation fellowship. We thank the UAVSAR team and the NASA-Dryden flight crew who collected and helped process the InSAR data, T. Jóhannesson for providing the Hofsjökull surface DEM, S. Owen for assistance with GIPSY-OASIS and F. Ortega for reviewing the manuscript.
