1. Introduction
This paper addresses the accuracy of three-dimensional ice velocities derived from a combination of interferometric synthetic aperture radar (InSAR) measurements, glacier flow modeling, and ice-sounding radar measurements of ice thickness. The approach to deriving three-dimensional ice-sheet velocities, is described in Reference Reeh, Madsen and MohrReeh and others (1999) and applied to Storstrømmen glacier, northeast Greenland in Reference Reeh, Mohr, Madsen, Oerter and GundestrupReeh and others (2003). The surface-parallel flow assumption approach described in Reference Joughin, Kwok and FahnestockJoughin and others (1996) is extended with a submergence/emergence term derived from the principle of mass conservation. For the ablation zone of a grounded glacier (ice of constant density from bed to surface), the equations for deriving the three-dimensional velocity become
where v
a and v
d are the projections of the unknown glacier velocity 
on the line-of-sight unit vectors 
for ascending and descending orbits, respectively. The surface normal, 
, is calculated by (−∂S/∂e,−∂S/∂n, 1), where S is the surface and (e, n) are Cartesian horizontal coordinates (Reference Reeh, Mohr, Madsen, Oerter and GundestrupReeh and others, 2003). The flux-divergence (FD) term on the righthand side of Equation (1c) is the submergence/emergence velocity, h(e, n) denotes ice thickness, 
is the horizontal surface velocity, and F is the column mean value divided by the surface value of 
.
One purpose of the present paper is to provide tools to establish a reasonable error budget for specific applications of InSAR to glaciological measurements. The second purpose of the paper is to present an error analysis for the dataset used in the companion paper (Reference Reeh, Mohr, Madsen, Oerter and GundestrupReeh and others, 2003), which presents three-dimensional flow results from Storstrømmen, based on ERS-1/-2 tandem data and airborne ice-sounding radar measurements. The paper is organized in three main parts:
-
1. Method descriptions (sections 2–5), providing the basic equations for the error analysis.
-
2. An assessment/application example (sections 6–8), describing how numbers are put into the basic equations.
-
3. Results (section 9), where the error analysis is summarized.
2. Method overview
The error analysis is carried out in three steps described in the following three sections. The error analysis assumes that error sources are so small that the output error is a linear superposition of the various error contributions. This assumption allows us to deal with the error sources independently.
Firstly, interferometric errors arising directly from v
a and v
d are assessed, assuming that the surface normal, , and the FD term in Equation (1c) are known.
Secondly, the errors in the underlying surface-parallel flow inversion are addressed. The coupling between Equations (1a–1b) and Equation (1c) is ignored, corresponding to assuming a level surface and ignoring the submergence/emergence velocity. This is acceptable for the error analysis.
Thirdly, the errors relating to the FD term in Equation (1c) are addressed. This term depends both on the interferometrically derived velocities and on the ice thickness derived from ice-sounding radar measurements.
After assessment of the direct error sources, the coupling of Equations (1a–1c) is addressed in “Results”, which also summarize our method. An overview of the different error sources considered and the corresponding sections is provided in Figure 1.
Fig. 1. Overview of the different error sources and mechanisms included in the error analysis. Numbers in parentheses refer to equations. Other numbers refer to sections in the method part of the paper and to the results section.
3. Method for interferometric errors
We assume that a double-difference method is used to obtain line-of-sight velocities and terrain elevations (Reference Joughin, Kwok and FahnestockKwok and Fahnestock,1996). In our approach we unwrap and calibrate each interferogram individually, by using ground-control points (GCPs). After calibration, the line-of-sight velocities and terrain elevations are found, assuming a stationary flow.
3.1. Geometric calibration
Geometric calibration (geocoding accuracy) is not addressed in the present paper. In Reference Mohr and MadsenMohr and Madsen (2001) it is shown that the European Remote-sensing Satellite (ERS-1/-2) mapping system can achieve a geometric accuracy corresponding to 10 m rms horizontally on the ground. This is not the final geocoding accuracy, as phase errors (static or noise) decrease the accuracy of the interferometrically derived elevations, which again cause a horizontal shift of the image in the cross-track direction. The horizontal error is 1/tan θ times the elevation error, where θ is the angle of incidence. Thus, a, typical elevation error of 10 m rms corresponds to a horizontal error of 25 m rms for ERS-1/-2 where θ ≈ 23°. This will degrade the co-registration of ascending and descending data. This is, however, negligible considering that ice velocities on glaciers in Greenland and Antarctica typically vary slowly spatially.
3.2. Interferometric path-length distortions
Both the troposphere and the ionosphere cause path-length distortions equivalent to path-length changes from surface displacements. Perturbations with magnitudes and spatial appearance similar to those of atmospheric path-length errors might additionally be caused by a dry uneven snow layer, which could potentially change the interferometric phase if superimposed between the two acquisitions (Reference Mohr, Reeh and MadsenMohr and others, 1998). Recently, Reference Guneriussen, Høgda, Johnsen and LauknesGuneriussen and others (2001) investigated this phenomenon, which is not described in detail in the literature. However, for our Storstrømmen glacier (northeast Greenland) example we disregard the effect (see section 6.3.2).
In Reference Mohr and MadsenMohr and Madsen (1996), a simplified interferometric measurement model is described and analyzed. In the following the resulting key equations for error analysis are summarized.
For convenience we use interferometric path lengh, φ′ for the error analysis, not interferometric phase, φ We define φ′ = λφ/(4π)+ B ||, where λ denotes radar wavelength and B || is the parallel baseline (corresponding to the flat-Earth phase). The direct effect of path-length distortions on the measured line-of-sight velocity, v, and elevation, h, depends on the orthogonal baselines B ⊥,1 and B ⊥,2 and temporal baselines T 1 and T 2 for the two interferograms denoted 1 and 2. For interferogram 1, the sensitivities are found to be
where R denotes slant range and θ the angle of incidence. Similar equations can be derived for interferogram 2. We here assume that the errors in the two interferograms acquired at different epochs are independent.
The sensitivity of the equation system (1a–1b) to errors in the line-of-sight velocities depends on the line-of-sight unit vectors. The line-of-sight unit vectors are parameterized by their angle of incidence (with respect to a level surface), θ, and angle of the ground swath track, ψ, with respect to north. An approximate expression for ψ, as function of the geographic latitude of the point of interest and its ground-track distance, is given in Appendix A. To simplify the error analysis, we assume that the point of interest is mapped at the same θ in both the ascending and descending pass, which again assures that ψ are numerically equal in both passes. Solving Equations (1a) and (1b) for v e, v n and v u yields
and
This shows that the east velocity, v e, can be found unambiguously from interferometric radar data alone. Thus, the surface-parallel flow assumption, i.e. Equation (1c), possibly modified with the FD term or another constraint, is required only to separate the north, v n, and up, v u, components. For the further error analysis we assume no vertical motion, which leads to
To summarize, the velocity errors caused by a path-length distortion can be assessed, by using Equation (3) or a similar equation for each of the two pairs of interferograms, followed by an application of Equations (4) and (6).
3.3. Interferometric phase noise
The Cramer–Rao bound for the interferogram phase standard deviation is
where N is the number of looks and γ the correlation coefficient (Reference Rodriguez and MartinRodriguez and Martin, 1992). The phase noise can be converted to an equivalent path-length noise by a division by 4π/λ, which subsequently can be used to assess velocity errors through the equations presented in section 3.2. The phase noise is not correlated from interferogram to interferogram. It is also uncorrelated from resolution cell to resolution cell (i.e. approximately from pixel to pixel), as opposed to the effects of the atmosphere and dry snow, which both are expected to exhibit large spatial correlations within one interferogram. A snow event between the two image acquisitions, though, would also increase phase noise since it inevitably causes decorrelation.
3.4. Other phase distortions
Phase unwrapping errors are not considered. Errors are usually detectable by inspection, as they appear as linearly fragmented features. In specific applications where unwrapping errors are expected, the error analysis could be carried out similar to that for atmospheric disturbances. Also, a possible clock drift (Reference Massonnet and VadonMassonnet and Vadon, 1995) in ERS-1/-2 is ignored. Again, the error analysis could be carried out similar to that for atmospheric disturbances. For small areas (e.g. one standard frame) clock drift is not expected to be significant.
3.5. Baseline calibration
In order to measure elevation and displacements using repeat-track interferometric data, the baselines have to be calibrated using GCPs. In the following, GCPs are assumed to be located on stationary terrain. The following sources of errors are considered:
Phase errors due to the atmosphere and dry snow.
Phase noise (described by the decorrelation).
GCP elevation errors due to digital elevation model (DEM) inaccuracies and misregistration between DEM and SAR data.
When a double-difference approach is used, and the same GCPs (with known line-of-sight displacement (e.g. zero)) are used for both interferograms, the GCP elevation errors will be identical in the two interferograms. Using equations (1) and (3) from Reference Mohr and MadsenMohr and Madsen (1996) it can be shown that such a GCP elevation error common to interferograms 1 and 2 will cancel out when the velocity is calculated. The derived elevations will, of course, be affected by the error in the GCP elevation. This decoupling of height references and velocity measurements is fundamentally different from the approach where the topographic phase is removed from one interferogram by using an independent DEM. In that case, the velocity error depends on the spatial and temporal baselines and the quality of the applied DEM.
In the following, a double-difference approach is assumed, so the velocity error from GCP elevation errors is ignored, leaving the phase noise and propagation delays as the error sources. It can be shown that a linear baseline variation in azimuth is approximately equivalent to a phase variation of the form
where x is the along-track coordinate, y the across-track coordinate, and a, b, c and d the constants to be determined by the calibration. The coefficients may be found by a least-squares estimation in a general linear model,
where 
is the phase errors at the GCP pixels, 
the phase distortions at the GCP pixels originating from phase noise, atmosphere, dry snow, etc., and n is the number of GCPs used. If the covariance matrix for the phase distortions at the GCP pixels is denoted 
, the covariance matrix for the parameter estimate 
, is 
, where superscript T indicates transposition. Thus, the standard deviation of the estimate of the phase correction value at a point (x; y) with 
is
assuming there is no correlation between the phase distortion at any GCP pixel and the phase distortion at the point of interest (x, y).This is the case for phase noise, but is often also a reasonable assumption for phase distortions caused by the troposphere and dry snow, provided that the GCPs (on stationary terrain) used to calibrate the interferograms are located at some distance (e.g. at least 10 km) from the areas of interest on the glacier. If the assumption of independence is violated, our method will overestimate the errors since both the direct effect of a phase distortion and the indirect effect through the baseline calibration are taken into account.
After calculation of the standard deviations of the phase errors caused by inaccurate baseline calibrations, the impact can be calculated by using a method similar to that used for path-length changes. It is further noted that the baseline calibration errors described by Equation (10) are highly correlated spatially. This is also the case for the combined effect of baseline calibration errors, which is a linear combination of the error in each of the four interferograms used (ignoring the stretch of the interferograms introduced during ground range transformation).
3.6. Non-stationary flow — same orbit family
The double-difference approach relies on an assumption of identical velocities during the observation periods for the two interferograms. For the purpose of analyzing the impact of a flow change, the velocity projections onto slant range, v 1 and v 2, corresponding to each of the two interferograms, are expressed in terms of an average value, v r, and a difference, vε , thus v 1 = v r − vε/2 and v 2 = v r + vε /2. Note that the velocity corresponding to an interferogram is the average velocity between the acquisition dates of the images composing the interferogram. If it is assumed, incorrectly, that vε is zero, the estimated velocity projection becomes
In the ERS-1/-2 tandem case, where T 1 = T 2 = 1 day, it is desirable to choose spatial baselines with opposite signs, as in that case |(B ⊥,1)+(B ⊥,2)/(B ⊥,1)+(B ⊥,2)| < 1. This implies that v m will be in the range [v 1 ; v 2]. In other words, the measured velocity projection is between the velocities at the two intervals of acquisitions. If spatial baselines are of equal sign, the ratio is larger than 1 and the combined estimate of the line-of-sight velocity thus more sensitive to flow changes. Optimal baselines are further discussed in Reference Mohr and MadsenMohr and Madsen (1996). The elevation estimate is similarly
Again, spatial baselines of opposite sign minimize the sensitivity to flow variations.
3.7. Non-stationary flow — different orbit family
A non-stationary flow between the ascending- and descending-orbit acquisitions also affects the direction and magnitude of the derived flow. To analyze the effect, the flow vectors corresponding to each look direction are expressed in terms of the average velocity 
and a difference 
, i.e. 
and 
. The measured line-of-sight velocities thus become 
and 
for the descending and ascending datasets, respectively. Assuming a level flow, the error velocity 
is written as 
, where η is the direction of the change. The resulting difference between the derived velocity 
and the average value 
becomes
(see Mohr and Madsen,1996). For glaciers with a seasonally varying velocity, the direction of flow is often constant, and only the magnitude changes. From Equation (13) it is seen that even given a constant flow direction, the measured flow will not necessarily be aligned with the flow. This is illustrated in Figure 2. It is also noted that the magnitude of the measurement error vector is never zero but always between 
and 
.
Fig. 2. Magnitude (a) and directional (b) deviation from the average flow with a ±5% flow magnitude change between ascending- and descending-orbit acquisitions. Track angles of 18°, 28° and 38° correspond to approximately 67°, 76° and 80° N for ERS-1/-2 at mid-swath.
4. Method for surface-parallel flow errors
The error in the surface-parallel flow term of the derived vertical velocity, v u, i.e. the v e(∂S/∂e) + v n(∂S/∂n) term in Equation (1c), is
Here the coupling between Equations (1c) and (5) is ignored.Since the velocity errors, 
and 
, in the east and north directions are generally different, the error generally depends on the direction of the terrain slope. The third and fourth contributions in Equation (14) are the surface slope errors multiplied by the actual horizontal velocities. With proper averaging, the dominating source of surface slope errors is usually not noise. With a good distribution of GCPs, baseline calibration errors are also of minor importance, due to large spatial correlation of the elevation errors caused by baseline errors. This leaves the atmospheric path-length distortions as the primary mechanism for creating slope errors. If atmospheric distortions, and thus slope errors, are assumed independent of direction, Equation (14) can be simplified to
where σ (∂S/∂x) denotes slope-error rms in one specific direction.
5. Method for flux-divergence errors
In this section, the errors in the FD, i.e. the righthand side of Equation (1c), are discussed. First, we describe how we calculate the FD. Then, errors arising from Fh uncertainties and 
uncertainties are discussed in two separate subsections.
5.1. Flux-divergence calculation
The FD, i.e. the submergence/emergence velocity, V, is calculated numerically, as
where E = Fhv e, N = Fhv n, and Δe and Δn denote the grid spacing in the east and north directions, respectively. It proved necessary to smooth the FD field (Reference Reeh, Mohr, Madsen, Oerter and GundestrupReeh and others, 2003) to obtain convergence of our solution. An analysis of the criteria for stability of the solution is beyond the scope of the present analysis. We use a (2m + 1) by (2m + 1) box filter, i.e.
For convenience we skip the (k, l) indices from 
, E and N, and assume a common grid spacing Δ
x
. We then find
Thus, the FD at the centre of the box filter only depends on the E and N values at the boundary lines of the filter andthe values at the adjacent points immediately outside the boundary.
5.2. Flux divergence — (Fh) errors
This subsection concerns errors in the derived vertical velocity, v
u, originating from the Fh term in the FD, i.e. 
in Equation (1c). For the error analysis we again ignore the coupling between Equations (1c) and (5). The error sources considered are:
Noise in h measurement.
Bias in h measurement.
Grid interpolation errors in h.
Variations of F associated with flow over an undulating base.
Bias in F due to changing basal sliding conditions.
Assumptions and notation for h and F errors, respectively, are described below, followed by an analysis of the effect of the different errors after the applied grid averaging (Equation (18)).
5.2.1. Ice-thickness errors
Ice-thickness noise is denoted εh and is assumed independent from gridpoint to gridpoint.
Ice-thickness bias is denoted βh . In order to simplify the error analysis, the bias is assumed constant for all gridpoints along one border of the filter (Equation (18)). The biases on different borders of the filter are assumed uncorrelated. A bias could, for example, be caused by spatial variations in surface and bottom echo amplitudes which again bias the echo location algorithm and thereby the derived ice thicknesses.
Ice thickness is assumed measured with profiling ice-sounding radar. This leaves gaps between the profiles which are filled by interpolation. Inspection of the measured thickness profiles shows an undulating behavior. For simplicity, we model these undulations as a sinusoidal with an amplitude d, and unknown angular frequency, ω = 2π/L, and unknown phase, φ
5.2.2. F-value variations
Flow over basal irregularities causes variations of the F factor, from large values over basal highs to smaller values over basal lows. In the following it will be shown that this reduces the impact of an unknown basal relief.
For the purpose of error analysis, we consider plane flow over an undulating base. The ice flux is assumed not to vary in the direction of the flow, i.e. flux changes due to surface mass balance are neglected. For ice-sheet conditions, which apply to Storstrømmen, the surface undulations are out of phase by π/2 and an order of magnitude smaller than the bottom undulations (Reference PatersonPaterson, 1994, p. 265). In that case, the ice thickness, h, can be assumed to vary as
where h 0 is the mean ice thickness, d is the ratio between the amplitude of the basal undulations and the mean ice thickness, ω = 2π/L, where L is the wavelength of the basal undulations, and s is the distance in the direction of flow. Using Reference ReehReeh (1987) it can be shown that in the plane-flow case the Fh product is approximately
where F 0 is the F factor of the mean flow and c a dimensionless parameter in the range [0; 1], depending on the h 0 /L ratio, ice properties, etc. A typical value is c = 0.5. Equation (20) shows that the variation over the undulation of Fh is reduced by a factor c as compared to the h variation.
If the ice-thickness distribution was known in a dense grid, a glacier flow model could in principle be used to calculate c and thus an accurate Fh value for each gridpoint. For the less detailed ice-thickness distribution of our Storstrømmen example, the latter term of the righthand side of Equation (20) is not known at the gridpoints, and therefore must be treated as an error term.
Bias in the F value is denoted βF . In order to simplify the error analysis, the bias is assumed constant for all gridpoints within the filter.
5.2.3. Effect of grid averaging
For the error analysis, the E and N terms in Equation (18) are treated separately but similarly. According to the discussion in section 5.2.2, we write F and h as
The ice-thickness noise, εh , is assumed independent from gridpoint to gridpoint. The phase of the unknown bottom undulations, φ, and ice-thickness bias, βh , are assumed fully correlated for all gridpoints along one border of the filter, but values on the different sides of the filter are assumed uncorrelated. The bias in F, βF , is assumed constant for all gridpoints.
Errors of v e are neglected for the moment, but will be treated in section 5.3. Ignoring second- and higher-order terms, we may express an Em,j value at a gridpoint (m, j) as
For clarity, indices of the error terms, i.e. φ, εh , βh and βF , are omitted. After some calculations (see Appendix B), we find
for the rms of the error in the vertical velocity, σ
V, caused by errors in the (Fh) term. In Equation (23), 
, 
and 
are the rms values of εh
, βh
and βF
, respectively. To simplify Equation (23), the center-point (or average) values F
0, h
0 and v
H are introduced. The terms 
and 
are the emergence/submergence velocities originating from east–west and north–south FD, respectively. In Appendix B, alternative expressions for 
and 
are given.
5.3. Flux divergence — interferometric errors
This subsection addresses errors in the derived vertical velocity, v
u, originating from the interferometric velocity term, 
, in the flux divergence, i.e. 
in Equation (1c). For the error analysis we ignore again the coupling between Equations (1c) and (5). We consider:
Noise; independent from gridpoint to gridpoint.
Undulations; with some point-to-point correlation.
Biases; i.e. long scale variations.
Noise on the horizontal ice velocity is denoted εv and εv for the east and north components, respectively.
Velocity undulations, Δv e and Δv n, of the east and north velocity components are here modelled as
where 
and 
are the rms values of the distortions, ω
e and ω
n the unknown spatial frequencies, and φ
e and φ
n the unknown phases. Undulations could, for example, be caused by atmospheric distortions or a dry snow layer.
Velocity biases are denoted 
and 
for the east and north components, respectively. In order to simplify the error analysis, the biases are assumed constant for all gridpoints. A bias could, for example, be caused by a baseline calibration error.
5.3.1. Effect of grid averaging
Similar to section 5.2.3, we write
for an E term. For clarity, indices of the error terms, i.e. 
, φ
e and 
, are omitted. The noise on the velocity, 
, is assumed independent from gridpoint to gridpoint. The phase of the unknown velocity undulations, φe
, is assumed fully correlated for all gridpoints along one border of the filter, but values on the different sides of the filter are assumed uncorrelated. The velocity bias 
is assumed constant for all gridpoints. Carrying out the analysis similar to that for the Fh term yields
where 
is the rms value of the horizontal velocity noise, 
the rms of horizontal velocity undulations, and 
and 
the rms of the velocity bias in the east and north directions, respectively.
6. Assessment of interferometric errors
In this section, errors in the interferometric measurements of the horizontal velocities and surface elevations are quantified. The aim is to compile a table with typical values for the impact of different error sources. ERS-1/-2 data from Storstrømmen, northeast Greenland (Reference Reeh, Mohr, Madsen, Oerter and GundestrupReeh and others, 2003), are used as an example. The characteristics of the data are summarized inTable 1.
Table 1. Characteristics of the ERS-1/-2 C-band (λ = 0.0566 m) example data from Storstrømmen (Reference Reeh, Mohr, Madsen, Oerter and GundestrupReeh and others, 2003)
The area covered by both satellite and ice-sounding radar measurements is approximately 30 km by 80 km.The glacier is building up from a surge. The vertical emergence velocities range from up to 6 m a−1 in the northwestern upper part to around 0 m a−1 in the southern lower part. The corresponding horizontal velocities are between 260 and 0 m a−1.
It is emphasized that it is not possible to compile a generally applicable error budget. The magnitude of the errors has to be calculated on a case-by-case basis based on acquisition geometry (baselines and latitude), processing approach (interferometric configuration, number of GCPs, etc.) and the characteristics of the glacier (penetration depth, flow direction, stationarity of flow, etc.).
6.1. Summary of interferometric error budget
The error budget for the direct effects of interferometric errors in the example data used is presented in Table 2. Each of the contributions is described in separate subsections below.
Table 2. Effect of different error sources using the example data fromTable 1
For reasons also described below, we consider the atmosphere and baseline errors (closely related to the atmospheric phase delays) the major error sources. With a 0.3 cm rms atmospheric path-length error, the combined effect of atmospheric and baseline errors in the area of interest is typically 6.2 m a−1 rms for the horizontal velocity and 16 m rms for the elevation. This is in line with the empirical results from Reference Mohr, Reeh and MadsenMohr and others (1998) and Reference Reeh, Mohr, Madsen, Oerter and GundestrupReeh and others (2003). It is noted that for larger areas of interest and for more humid climates, significantly larger atmospheric path-length distortions, and thus also baseline calibration errors, could be expected.
6.2. Geometric calibration
The effect of misregistration between ascending- and descending-orbit interferograms is ignored as glacier velocities in Greenland and Antarctica typically vary slowly. However, it should be noted that the processing system used in Reference Reeh, Mohr, Madsen, Oerter and GundestrupReeh and others (2003) was not fully calibrated geometrically, so the co-registration error is likely to be slightly larger than estimated here.
6.3. Interferometric path-length distortions
In this subsection, path-length distortions caused by the atmosphere and dry snow are assessed.
6.3.1. Atmosphere
Atmospheric path-length distortions may originate from both the ionosphere and the troposphere (Reference Tarayre and MassonnetTarayre and Massonnet,1996).
Ionospheric irregularities typically have greater spatial correlation lengths than tropospheric irregularities (Reference Hanssen, Weckwerth, Zebker and KleesHanssen and others, 1999). Over the small scales relevant in the present study, ionospheric distortions are only significant during unusual conditions, as the interferometric product is calibrated with GCPs spaced at most a few tens of kilometers from the study area (Reference Gray, Mattar and SofkoGray and others, 2000). Thus, in the present error analysis, ionospheric distortions are ignored as a separate error source, but are lumped together with the tropospheric path-length distortions.
The spatial characteristics of tropospheric path-length distortions are significantly different from day to day (Reference HanssenHanssen, 1998). Typical values for the spatially non-linear phase distortions (tropospheric and ionospheric) over a 21 km by 21 km arctic area are found to be ±0.5 cm peak to peak, corresponding to ±1 rad (Reference Mattar, Gray, Geudtner and VachonMattar and others, 1999). This is consistent with other reports. In the following analysis, we use 0.3 cm rms for the atmospheric path-length variation.
Given 0.3 cm path-length distortions in each of the four example interferograms, the effect on the decomposed line-of-sight velocities and elevations is calculated from Equations (2) and (3). The results are given in Table 3. The combined rms error is, using Equations (4) and (6) and thus assuming surface-parallel flow, found to be 2.1 m a−1 in v e and 3.9 m a−1 in v n. Thus, the total horizontal velocity error is 4.4 m a−1 rms, as listed inTable 2. In our approach, the surface elevation is determined from the ascending-orbit interferograms only, and thus the total elevation error is 9 m rms as listed inTable 2.
Table 3. Effect of a 0.3 cm path-length error (i.e. 0.67 rad) applied successively to each of the Storstrømmen interferograms inTable 1
The small spatial baselines of the descending-orbit data of course reduce the accuracy of the elevations (Equation (2)). However, since the baseline ratio is close to zero, the accuracy of the velocities is not severely affected (Equation (3)). The geocoding accuracy, which is important when ascending- and descending-orbit data are combined, is degraded, but this could in principle be circumvented by using the ascending-orbit DEM as topography reference for the velocities derived from descending-orbit data.
6.3.2. Dry snow
The effective path-length delay for a pixel is determined by the average value of the snow-layer thickness difference between the individual SAR acquisitions. Thickness variations within an interferogram pixel, which typically has a size of 50 m by 50 m, or more, will cause phase decorrelation. Since the interferograms are calibrated with GCPs, it is the spatial variations of the average snow layer from pixel to pixel which cause path-length distortions, not the average value.
The standard formula for the relative dielectric constant, ε, as a function of the snow density, ρ s, for dry non-dense snow is ε = 1 + 0.0019ρ s (Reference Ulaby, Moore and FungUlaby and others,1986, p.2063), which is consistent with Reference MätzlerMätzler (1996). A wind-compacted snow layer with a density of 250 kg m−3 is hence found to have a relative dielectric constant of approximately 1.5 corresponding to an index of refraction of 1.2. The index of refraction decreases to 1.1 for newly fallen snow with a density of 100 kg m−3. We here proceed using a 10 cm snow layer with an index of refraction 1.2, which thus causes a 2 cm delay difference.
At a glacier like Storstrømmen, it is expected that a dry snow layer will be distributed unevenly. Assuming a rectangular distribution between 0 and 10 cm for the averagesnow layer from pixel to pixel, the path-length variation is thus 
. The effect of such a differential snow layer simultaneously but statistically independent in each of the four interferograms is shown inTable 2.
However, it seems reasonable to ignore the dry-snow contribution as a significant error source for the example data. We did not observe significant decorrelation in any of the interferograms, which would indicate a precipitation event. This is considered the normal situation, since the snowfall on Storstrømmen is limited (most likely <0.5 m snow per year) and since the precipitation events only affect the interferometric measurements if coincident with one or more of the four 1 day periods of SAR acquisition.
6.4. Interferometric phase noise
For the Storstrømmen example, the number of looks is approximately 20, large enough for Equation (7) to be used to represent phase noise. Observed correlation values for rock and ice are provided in Table 1 for each of the four interferograms. Converting the correlation values to phase noise, combined with the results presented in Table 3, yields the estimated phase-noise error terms presented in Table 2. It should be noted that these values are applicable for the comparison with pole measurements presented in Reference Mohr, Reeh and MadsenMohr and others (1998), but during the derivation of the vertical emergence/submergence velocity presented in Reference Reeh, Mohr, Madsen, Oerter and GundestrupReeh and others (2003) a further averaging was applied, which reduced the phase variance with a factor on the order of 10.
6.5. Baseline calibration
For the example data, the number of GCPs is 6. Although two GCPs are located within 8 km of each other and two within 11 km, the phase distortions are assumed uncorrelated between the pixels where GCP calibration is performed, i.e. 
, where 
is the unity matrix of order n. Note that for a denser GCP network, it is important to include the co-variance structure of the (atmospheric) phase noise in Σ. This structure is discussed, for example, in Reference Williams, Bock and FangWilliams and others (1998). In Figure 3, the spatial variation of 
is shown assuming σ = 1 rad at the six GCPs in the example data. It is observed that in the area enclosed by GCPs the standard deviation of the calibration residual is typically 1 rad. In line with Reference Joughin, Kwok and FahnestockJoughin and others (1996), it is also seen that outside the area enclosed by GCPs the errors increase rapidly.
Fig. 3. Storstrømmen glacier, northeast Greenland. Phase (baseline) calibration errors in radians in one example interferogram assuming a 1 rad uncorrelated rms noise at the GCPs. All four example interferograms are calibrated using the six GCPs indicated by the grey “+”symbols. The location of the area is indicated by the insert. The light-gray area is covered by SAR interferometry. The dark gray area is additionally covered by ice-sounding radar. Arrows indicate line-of-sight directions from ascending and descending orbits at the center of the area.
For the evaluation of the baseline contribution to the velocity errors, only phase noise and atmospheric disturbances are included. Disturbances due to a dry snow layer are unlikely, as argued above. From Figure 3, it is seen that within the area in the example data enclosed by GCPs, the baseline error is typically equal to the error at the GCPs. A typical value for the rms velocity errors is thus estimated to be equal to the velocity errors originating from the atmosphere and phase noise (Table 2). In most cases, though, we use a spatially varying velocity error (of baseline calibration origin), based on the factor 
times the combined error originating from the atmosphere and phase noise.
With respect to elevation, the elevation errors on the GCPs also have to be included. Assuming a 10 m rms error on the GCPs, the corresponding phase error terms are, using the baselines from Table 1, found to be 0.92 rad and 0.13 rad for A1 and A2, respectively. Again ignoring the dry-snow contributions, adding the phase noise and atmospheric disturbances and using the results fromTable 3 yields an elevation rms error of 13 m originating from improper calibration of baselines (Table 2).
6.6. Non-stationary flow — same orbit family
For the example data, the baseline values for the ascending-orbit data have different signs (Table 1). This implies that in the case of a non-stationary flow, the derived line-of-sight velocity (displacement rate) will be a weighted average of the velocities at the two interferogram observation periods (Equation (11)). Thus, the derived velocity has been present at least at one instant in the period from the acquisition of the first to that of the last SAR image. In other words, the line-of-sight velocities derived from the ascending-orbit interferograms are correct for at least one specific but unknown epoch in the period between the acquisitions. The same argument applies to the descending-orbit data. The baseline ratios of 1/7 and 1/19 for ascending- and descending-orbit data, respectively, though, give large weights to the velocities at the acquisitions with short spatial baselines.
A non-stationary flow will, on the other hand, change the measured elevations. The change of line-of-sight velocity being known, the impact on the measured elevation can be calculated from Equation (12). For the Storstrømmen example, it is unknown whether or not the velocity changes during the winter, since there are no global positioning system (GPS) measurements from the winter. It is known that the summer velocities over much of the glacier are larger than the annual means.We believe that the increased summer velocity is caused by meltwater penetrating to the bottom and thereby enhancing sliding. During winter (September–May) there is no meltwater, so we assume a negligible seasonal velocity variation during the observation period. GPS measurements of surface velocities during the observation period (1992–95) show a general decrease on the order of 10 m a−1 per year (Reference Reeh, Mohr, Madsen, Oerter and GundestrupReeh and others, 2003). If assumed to be uniform, this is a decrease of 2 m a−1 during the observation period.
Using Equation (12), and the interferogram parameters from Table 1, it is found that a 1 m a−1 change in the across-track component of the horizontal velocity causes a −2.3 m change in the measured elevation. Changes in the along-track direction have no effect, as these do not change the line-of-sight velocity. If the flow change is averaged over all directions, the rms error is reduced to 
. Thus, a 2 m a−1 rms horizontal across-track velocity change corresponds to 3.2 m rms (Table 2).
6.7. Non-stationary flow — different orbit family
From Equation (13) it is seen that the ratio of the flow measurement error and the flow change 
is between tan ψ and 1 / tan ψ. The direction of the flow measurement error is not necessarily aligned with the direction of the flow change. Thus, for our study area, where ψ ≈ 28° (see Equation (27)), a 1 m a−1 flow change between the ascending- and descending-orbit acquisitions will at most result in a 0.9 m a−1 deviation from the average value.
For Storstrømmen, flow changes are expected to be primarily in the direction of the flow. This is supported by repeated measurements of positions of poles drilled into the ice. The pole measurements also suggest that the potential flow changes are rather independent of flow magnitude, except for the areas with very small or zero velocities. Since the flow directions are varying over a large azimuth interval, the correct approach would be to compile a map with the effect of a potential flow change. However, as shown below, the effect is expected to be less important than atmospheric distortions and baseline calibration errors. Therefore, we simplify the error analysis to provide one number for the effect instead of a map. This is done by estimating a value for the standard deviation of the potential flow change and assuming an unknown direction. For other glaciers, with a more unidirectional flow, it might be more appropriate to assume a fixed direction of the flow change.
If a flow change, 
, is assumed having a magnitude of 1 m a−1 and averaging over all directions, the east and north errors are found to be tan 
and cot 
, respectively,for the Storstrømmen example (Equation (13)). As described above, no direct information is available on the flow changes during winter. Again, a 10 m a−1 change per year is assumed, corresponding to 2.8 m a−1 in the period between the ascending- and descending-orbit interferograms with the shortest spatial baselines. Thus, the errors due to flow changes are estimated to be 0.5 m a−1 in the east direction and 1.8 m a−1 in the north direction, i.e. 1.9 m a−1 horizontally (Table 2).
Since the non-stationary flow contribution is small compared to the atmosphere and baseline errors and since we have no clear evidence for the phenomena, we do not consider the non-stationary flow a major error source for our example data.
7. Assessment of surface-parallel flow errors
Equation (15) describes the error in the vertical velocity originating from the surface-parallel flow term, i.e. the lefthand side of Equation (1c). It depends on the magnitude and error of the horizontal velocity and the error in the surface normal.
For the example data, the surface slope is generally <0.01, but it reaches 0.06 in a few areas.This implies that the 6.2 m a−1 rms on the horizontal velocity scales down to <0.36 m a−1 rms on v u, and in most areas even less than 0.06 m a−1.
The surface slope errors are more difficult to quantify, even when the phase noise and baseline calibration errors are ignored, so that distortions of the surface elevation measurements caused by the atmospheric path-length changes are the primary error source. Assuming a sinusoidal distortion with a wavelength L, 
amplitude, and averaging over all directions of the slope error, the slope error is found to be 2πA/L rms in one specific direction. For A = 9 m and L = 10 km, the slope error is thus 0.0057 rms. For horizontal velocities of 0–260 m a−1 in the study area, the resulting error in v
u is thus 0–1.5 m a−1 rms.
The combined effect of east velocity errors, north velocity errors and slope errors is shown in Figure 4a. Note that in this figure the spatial variation of the baseline error (see Equation (10)) is taken into account. The values are 0–1.4 m a−1 in the area of interest.
Fig. 4. Vertical velocity errors from the surface-parallel flow and FD terms, for the Storstrømmen example. (a) Rms of vertical velocity errors (see Equation (15)), originating from the surface-parallel flow term.We use 
(atmosphere), 
(baseline), 
(atmosphere), 
(baseline) and σ∂
S/∂x = 0.0057. (b) Rms of vertical velocity errors (see Equation (23)), originating from the Fh term in the flux divergence. In the calculation of flux divergence we use F0 = 0.95 (constant), a grid with Δ
x
= 500 m, and a box filter with m = 10.The errors assumed are 
(normal distribution), 
(rectangular distribution), 
(rectangular distribution), c = 0.5 (typical value) and d = 0.1 (i.e. ±10%). (c) Rms of vertical velocity errors (see Equation (26)), originating from the interferometric velocity term in the flux divergence.With the same inversion parameters as in (b) we use 
(phase noise), 
(atmosphere), 
(baseline), and 
(baseline) (Table 2). All maps are in Universal Transverse Mercator (UTM), zone 27, with the lower left corner at (435 000, 8 515 000), and 10 km between grid ticks.World Geodetic System 1984 (WGS84) datum. Dots and numbers indicate positions of stakes used in the companion paper (Reference Reeh, Mohr, Madsen, Oerter and GundestrupReeh and others, 2003).
8. Assessment of flux-divergence errors
8.1. Flux-divergence error assessment — (Fh) errors
In section 5.2 above, five types of errors are considered: noise, bias and grid interpolation errors in ice thickness h, errors due to undulations in F, and bias in F. The resulting error in the vertical velocity can be calculated by using Equation (23). In this subsection, the magnitude of each contribution, required for an application of Equation (23) is assessed, and the results for our Storstrømmen data shown in Figure 4b.
8.1.1. Noise in ice-thickness measurements
The top and bottom echoes from ice-sounding radars are located using the shape of the received (and averaged) radar pulses. The accuracy depends on the thermal noise. The varying signal to clutter level of the bottom echo might also disturb the location algorithm, in effect adding noise. For our ice-radar measurements we assume a radar noise term of 6 m rms (Reference Christensen, Reeh, Forsberg, Jørgensen, Skou and WoeldersChristensen and others, 2000).
In addition, there are uncertainties from short-scale variations of h that are smoothed out in the interpolation procedure we use to form a 500 m by 500 m grid of elevations from the ice-thickness profiles. Comparing the measured ice thickness along the radar flight-lines with the corresponding smoothed ice thickness of the gridded model shows a standard deviation of 12.5 m. The 12.5 m includes both real variations in the ice thickness and noise in the ice radar measurements. Since the gridding procedure includes averaging, the error at a gridpoint is primarily due to ice-thickness variations, and thus smaller than 12.5 m, since a part of the 12.5 m is caused by radar noise.
We will, however, use 
and assume it to be uncorrelated from gridpoint to gridpoint.
8.1.2. Bias in ice-thickness measurements
There are several effects which potentially cause biases in the ice-radar measurements. A firn/snow layer with unknown thickness and dielectric properties, and thus an unknown wave propagation speed, could be present, but this is not the case for the example dataset, which is in the ablation zone. Water on the surface or water percolating into the ice could also change the dielectric constant. During the ice-radar measurements, lakes were present at the surface, but corresponding disturbances in the radar echoes were not detected. Varying surface and bottom echo amplitudes may also bias the echo location algorithm. Finally, bottom undulations cause a bias as the radar detects the closest echo, not the nadir echo. In total a bias of up to 10 m is (somewhat arbitrarily) assumed. Assuming a rectangular distribution, the rms becomes 
.
8.1.3. Grid interpolation error in ice thickness
Inspection of the measured thickness profiles in the upper region of the glacier, where the ice thickness is 600–800 m, shows that sub-ice relief is characterized by amplitudes of approximately 100 m with wavelengths of approximately 10 km. Nearer to the glacier terminus, both the ice thickness and the amplitude of the basal relief decrease. As the distance between the ice-radar flight-lines is 5–10 km, thickness variations in between the flight-lines on these length scales will not be picked up by the gridding procedure, suggesting that the relative error of a grid value of h could be 10% or more. On average we assume a magnitude of the relative thickness undulations of d = 0.1.
8.1.4. Short-scale undulations of F
Flow over basal irregularities causes variations of the F factor, which reduce the effect of unknown ice-thickness variations by a factor c (Equation (20)). We assume a constant value of c = 0.5.
8.1.5. Long-scale variations of F
At Storstrømmen, variations of F with spatial correlation lengths of tens of kilometers are expected due to changes in the basal sliding conditions. The observed increase of the ratio of summer velocities to mean annual velocity along the glacier (Reference Reeh, Mohr, Madsen, Oerter and GundestrupReeh and others, 2003) indicates that, when approaching the glacial terminus, basal sliding constitutes an increasing fraction of the forward motion of the glacier. Accordingly, the F factor will increase along the glacier during summer months. A variation from F = 0.9 (no basal sliding) to F = 1 (fully developed basal sliding) can be expected. The InSAR velocities are measured in the winter, when observations indicate that basal sliding is less developed (Reference Reeh, Mohr, Madsen, Oerter and GundestrupReeh and others, 2003). Nevertheless, for the calculation of the FD term, we use a constant F = 0.95, and assume abias of ±0.05 with a rectangular distribution, i.e. F
0 = 0.95 and 
in Equation (23).
8.2. Flux-divergence error assessment — interferometric errors
In section 5.3 above, three types of errors are considered: noise, undulations, and bias of the horizontal velocity, 
. The resulting error in the vertical velocity can be calculated by using Equation (26). The magnitude of each contribution required for an application of Equation (26) is assessed in this subsection, and the results for our Storstrømmen data shown in Figure 4c.
8.2.1. Noise
Phase noise is assumed to have a normal distribution. The effect of the noise is summarized in Table 2, i.e. 
.
8.2.2. Undulations
Of the error sources considered in section 6 above, the atmosphere and dry snow potentially cause velocity undulations. For the reason explained in section 6.3 above, the effect of dry snow is ignored. This leaves the atmosphere as the primary error source, i.e. Av H = 4.4 m a−1.
8.2.3. Biases
The effect of a baseline calibration error is independent of the actual glacier velocities.The baseline error (see Equation (10)) is slowly varying within the FD filter area, but this variation is neglected. The effect of baseline calibration errors is summarized in Table 2, but here we take into account the spatial variation described by Equation (10). We use 
and 
.
9. Results
This section describes how the final error budget for the three velocity components, v e, v n, v u, is compiled.
The first step is to compile a table similar to Table 2 using the approach described in section 6 above. From this analysis the error in v e can be found. For our example dataset, the major error source is the atmospheric disturbances which also indirectly affect the accuracy through the spatially varying baseline error (Equation (10)). The rms value of the v e error is shown in Figure 5a.
Fig. 5. Final (ve, vn, vu) error estimates for the Storstrømmen example data. (a) Rms of the east velocity (ve) errors found by combining a constant error of 2.1 m a−1 (atmosphere) and a spatially varying error of 
(baseline). (b) Rms of the north velocity (vn) errors found by combining a constant error of 3.9 m a−1 (atmosphere) and a spatially varying error of 
(baseline) with the vertical velocity error from (c) scaled by a factor cot 23°/sin 28° ≈ 5.0. (c) Rms of the c vertical velocity component (vu) errors found by combining the results from Figure 4a–c. Map projection as in Figure 4.
The error in the v n and v u components requires further analysis since they cannot be found from the interferometric measurements alone (Equation (5)). Our approach to separate the v n and v u components includes a surface-parallel flow assumption. The corresponding error in the vertical velocity can be estimated using Equation (15), in conjunction with results fromTable 2 and Equation (10). Dependent on the origin of the surface DEM, the analysis may differ slightly from the one described in section 7. The result is shown in Figure 4a.
If the FD assumption is also applied, we estimate the corresponding errors in the vertical velocity by using Equations (23) and (26), in conjunction with results from Table 2 and Equation (10). The results are shown in Figure 4b and c. The error in v u can now be calculated as the sum of the three vertical velocity error terms from Figure 4. The result is shown in Figure 5c.
With the above approach, the coupling between v n and v u is ignored. In the error analysis, we account for the coupling by assuming that v u errors are independent of the direct interferometric errors in v n. In the error assessment for v n, the coupling between v u and v n is then accounted for by adding the v u errors scaled by a factor cot θ / sin ψ (see Equation (5)) to the direct effects of errors in the interferometric measurements.
For the direct effects of interferometric errors in our example data, we use a constant of 3.9 m a−1 for the v
n error(atmosphere) and a spatially varying v
n error of 
(baseline). The scale factor for v
u errors is cot 23° / sin 28° ≈ 5.0, resulting in a v
n error as shown in Figure 5b.
10. Conclusions
In this paper, a model for a differential interferometric radar system has been developed. It is shown that by using observations from two look directions (e.g. ascending- and descending-orbit data) one velocity component, typically the east component v e, can be found unambiguously. In order to separate the two other components, typically v n and v u, additional information or assumptions are required (e.g. an assumption of surface-parallel flow). With the combination of interferometric measurements from different epochs, a stationary-flow assumption is furthermore needed. It is shown that if the spatial baselines are of different sign, non-stationary flow between acquisitions from the same orbit family (e.g. ascending) results in a measured displacement which has actually occurred at a specific but unknown epoch between the two observations. The error arising from non-stationary flow between different orbit family acquisitions is also quantified.
Relations between interferometric phase errors and displacement and elevation measurements have been developed. These relations include an equation which relates the phase noise and the position of GCPs to corresponding phase errors after baseline calibration. This is used to compile an error budget for our example data from Storstrømmen (Reference Reeh, Mohr, Madsen, Oerter and GundestrupReeh and others, 2003), including the effect of the atmosphere and phase noise.
A relation between surface slope and velocity errors and the corresponding surface-parallel-flow vertical velocity has been developed. Also, a method for quantifying the vertical velocity error corresponding to the FD term was established. This includes both errors from the Fh term and errors from the horizontal velocity, v h.
For our Storstrømmen example, the rms error of the FD vertical velocity component relative to the magnitude of this component is typically < 30% in areas with FD components >1.0 m a−1. The magnitude of the FD component is between −2.3 and 6.2 m a−1, with typical values between −0.5 and 2.5 m a−1.
Acknowledgements
E. L. Christensen (Technical University of Denmark) is thanked for information on the ice radar system, and N. Gundestrup (Department of Geophysics, University of Copenhagen, Denmark) is thanked for information regarding the ice radar campaign on Storstrømmen. We also thank the two anonymous reviewers for their helpful review of the paper, and Scientific Editor H. Rott for all his help during the preparation of the manuscript.
Appendix A
The sensitivity of the three-dimensional interferometric measurements depends on the angle of incidence, θ, and the angle of the ground swath track, ψ, for the ascending- and descending-orbit data which are combined.
The ground swath track angle, ψ, measured clockwise relative to north, depends on the target latitude, φ, the orbit inclination, i, the angle, y′ = y/R e, between the ground swath track and the nadir track, where y is the ground swath track distance (positive to the right), the Earth radius, R e, and the Earth rotation, ω E. If Earth rotation is ignored, the Earth assumed spherical, and the orbit circular, the track angle, ψ s, can be calculated from
where “+”applies to ascending and “−” to descending orbit arcs. Note, ψ is in the interval −90° to +90°, so it will change sign at the northern- and southernmost points of the orbit. The Earth rotation contributes an additional latitude-dependent term, with a numerical value of approximately 3.5° at the Equator, decreasing gradually to zero at the northern-/southernmost points of the orbit. The resulting track angle, ψ, can be approximated by
where v E = ω E R E denotes the Earth rotation velocity at the Equator, v S is the satellite ground velocity, and “−” and “+” apply to ascending and descending orbit arcs, respectively. The accuracy of Equation (A1) is sufficient for error analysis, but not for actual three-dimensional velocity decomposition, which should be based on the actual state vectors.
Appendix B
In this appendix an expression for the error in the vertical velocity caused by uncertainties in Fh in the FD term of Equation (1c) is derived.
Because of the correlation of some of the error sources, it is convenient first to consider the variance of the sum of E values along lines parallel to the filter boundaries. According to the discussion of the Fh errors in section 5, we may write the sum of Em,j , given by Equation (22), along one of the filter boundary lines as
where the relevant values for p are {− (m + 1), −m, m, m + 1}.The bottom undulation term is described by the constants c, d, ω = 2π/L, and Δ
n
and an unknown phase, φp
, assumed equal for adjacent p values and uncorrelated from one side of the filter to the other. The ice-thickness noise, 
, is assumed independent from gridpoint to gridpoint. The ice-thickness bias, 
, is assumed equal for adjacent filter boundary lines but uncorrelated from one side to the other. The F bias, βF
, is assumed constant over the entire filter.
First the cosine error term is considered. Replacing Fp,j
, hp,j
and 
with their mean values along the boundary line Fp
,0, hp
,0 and 
approximating the summation with an integral, and introducing the filter length relative to the wavelength of the basal undulation, L′ = 2mΔ
n/L, yields
where the worst-case value of 1 for sin(x)/x is assumed. With an unknown phase of the basal undulation the variance of the term is 
.
For the remaining error terms it is also acceptable to replace Fp,j
, hp,j
and 
with their mean values. Taking the correlation into account, the total variance, 
of a ∑Ep,j
term becomes
where 
and 
are the rms of the noise and bias for an individual gridpoint, respectively, and 
the rms of the F bias. In order to simplify the notation in Equation (B3), we have omitted the subscript 0 in all quantities.
Analogous expressions can be derived for the variance of the other ∑Ep,j
terms in Equation (18). The contribution to the variance, 
,of the smoothed FD terms from all four ∑Ep,j
terms is found by combining the variances with due consideration of the degree of correlation between the different sub-contributions. For calculation of the variance contribution from the bias term, βF
, which has a correlation coefficient of 1 within the entire filter, we use the approximation F
+ = Fm
= Fm
+1 and F
− = F
−m
= F
−m−1, and similar approximations for hp
and 
. For the ice-thickness bias, βh
, and the ice-thickness noise, which are uncorrelated from one side of the filter to the other, we further approximate Fp
, hp
and 
by the central values F
0, h
0 and 
. We find
The variance of the sum of the four N terms of Equation (18) is found similarly with Δ
x
= Δ
e
= Δ
n
and 
the total Fh contribution to the variance becomes
In our case, where we assume a constant F within the entire area, the last two terms can be modified by multiplying with 
in both the numerator and the denominator. Using the constant F
0 approximation, the east term becomes
where 
is the flux divergence (i.e. the submergence/emergence velocity) averaged from one end of the filter to the other. Likewise for the north term.
