Evidence for non-tidal diurnal velocity variations of Helheim Glacier, East Greenland

Abstract

We have used tide-gauge data from near Helheim Glacier, East Greenland, and GPS data acquired on the glacier to investigate the spectra of tidal forcing and flow response. For both the tide-gauge and GPS time series, we calculated amplitudes and phases for a harmonic expansion using a limited set of harmonic constituents. We find that the semidiurnal constituents of the glacier flow are well modeled using a single admittance and lag with respect to the tide-gauge data. However, diurnal variations in the glacier flow cannot be simply described using this model. We find an additional signal in glacier position, in phase with the S₁ solar tide, with some modulations at other frequencies. These non-tidal variations account for a peak-to-peak variation in glacier flow speed at a site close to the terminus of ~0.7 m d^-1, compared with a mean flow rate at this location of ~22 m d^-1. The speed variations reach their daily maximum value ~6 hours after local noon. We hypothesize that these additional diurnal variations are associated with peaks in lubrication of the glacier bed due to surface melting driven by diurnal solar heating.

Introduction

Glaciers are subject to periodic or quasi-periodic external forcing from a number of sources. For marine-terminating glaciers, these sources include variations in back pressure on the face of the calving front exerted by ocean tide fluctuations (e.g. Meier and Post, 1987; Walters and Dunlap, 1987; Vaughan, 1995; Anandakrishnan and Alley, 1997; de Juan and others, 2010; Marsh and others, 2013). Another class of periodic behavior is stick–slip motion modulated by stress changes due to ocean tides (e.g. Bindschadler and others, 2003; Zoet and others, 2012).

Hydrological forcing driven by changes in water input is also a major source of periodic ice-flow variations. Glacier-flow responses to meltwater forcing have long been known to occur at alpine glaciers (e.g. Iken and Bindschadler, 1986; Bartholomaus and others, 2008), where diurnal and seasonal velocity variations are observed. In Greenland, seasonal signals observed on the ice sheet have been inferred to result from meltwater variations (Zwally and others, 2002; Van de Wal and others, 2008), and seasonal and diurnal signals have been observed at land-terminating glaciers (Shepherd and others, 2009; Bartholomew and others, 2010). At Greenland’s marine-terminating outlet glaciers, some seasonal response has been observed (Joughin and others, 2008; Howat and others, 2010), and day-to-day speed variations attributed to meltwater forcing have been described (Andersen and others, 2010, 2011). However, to our knowledge, diurnal variations in flow due to meltwater variation have not been documented at these marine-terminating glaciers.

Study of the response of tidewater glaciers to tidal forcing provides useful insight into their rheological properties, conditions at the glacier bed and the hydrology of these systems. However, separation of a diurnal velocity signal due to meltwater effects, tidally driven variations and short-term responses to calving has proven challenging (e.g. Walter and others, 2012). In part this is because all of these forcings will contain energy at or near the frequency of the solar constituent S₁, and harmonics thereof. The S₁ tidal harmonic is associated with the longitude of the Sun in the terrestrial coordinate system and is present as a component of the ocean tide, mainly due to diurnal atmospheric pressure loading (e.g. Ray and Egbert, 2004). We also expect a harmonic at the S₁ frequency to be a primary component in a melt signal driven, ultimately, by solar heating. Other spectral characteristics of tidal hydrological forcing also appear similar to forcing driven by melt. For instance, the tidal lines will be broadened, since the tidal spectrum is made up of a large number of individual species (e.g. Kudryavtsev, 2004), and, since meltwater production is primarily driven by temperature and insolation variations that have stochastic-like fluctuations, it will have a broadened line centered on S₁ and a reddish background continuum (Davis and others, 2012).

Given the importance of separating the tide- and melt-driven components of glacier flow, a number of approaches have been tried. When tide observations are available, comparison of the amplitudes of tidal forcing and glacier response has been fruitful. O’Neel and others (2001) pointed out that by assuming that the response of the glacier to tidal forcing was constant across the spectrum, the melt signal could be estimated. In other work, de Juan and others (2010) used a time-domain approach with a similar assumption, and found that a reasonable fit to a linear tidal admittance model could be achieved only through inclusion of parameters representing the daily acceleration of the glacier. Such a model assumes a linear response of glacier motion to tidal forcing, with admittance being the ratio of response to forcing.

In this paper we describe harmonic analyses of both the relative sea-level record for Sermilik Fjord and position estimates from GPS observations on Helheim Glacier obtained during 21 days in 2007. The analyses are used to test the assumption that the semidiurnal flow variations can be described using only a single sea-level admittance, and to evaluate rigorously whether non-tidal signals that occur at or near tidal frequencies can be extracted robustly from the glacier-flow time series. We use the differences between the sea-level and position spectra to identify the non-tidal glacier-flow signal and to estimate its amplitude. We then use a Kalman filter to examine the temporal variability of the non-tidal signal.

Data

During each of the boreal summer seasons of 2006–08 we deployed an array of instrumentation on and near Helheim Glacier. In this study, we focus on data from 2007, and compare measurements of sea level and estimates of time-dependent glacier position inferred from analysis of GPS data. The collection of these datasets has been described elsewhere (Nettles and others, 2008; de Juan and others, 2010).

Tide-gauge data from Sermilik Fjord

We infer sea level from pressure-gauge observations acquired in a cove off Sermilik Fjord, at a distance of ~35 km from the terminus (Fig. 1) of Helheim Glacier. These data were recorded continuously at 1 min intervals during days 203–220 of year 2006 and 186–218 of 2007. The pressure gauge (hereafter referred to as a ‘tide gauge’) has a nominal resolution of ~1 mm.

Fig. 1. Landsat image of Helheim Glacier, southern Greenland, acquired 2001, showing the location of the GPS site IS22. Black dotted and dashed curves show the positions of the calving front at two times (4 July and 15 August) during the summer of 2007. Top left inset: Arrow marks the location of Helheim Glacier in southern Greenland, with white dots showing locations of glacial earthquakes (Tsai and Ekström, 2007). Top right inset: Map of the region showing Sermilik Fjord, and locations of the tide gauge and the tide prediction point. (Figure after de Juan and others, 2010.)

We focus on the sea-level time series for 2007 (Fig. 2a). A purely visual inspection indicates that the variations during this period are dominated by a semidiurnal component, modulated at ~15 days. There is also a smaller diurnal component that contributes to the pattern of semidiurnal maxima with alternating large and small amplitudes. This particular pattern persists even when the sinusoidal amplitude is small, due to the 15 day modulation (i.e. at neap tide). Shorter periods of high-frequency variability are observed during and following periods of large calving events associated with glacial earthquakes, as seen around days 189–190 (Nettles and others, 2008).

Fig. 2. Sea-level observations from Sermilik Fjord in 2007, along with a comparison with an open-ocean tide model. (a) Blue: observations. A mean value of 2.176 m has been subtracted. Red: AOTIM-5 ocean tide model (Padman and Erofeeva, 2004) calculated at a location outside Sermilik Fjord (Fig. 1) for days 185–220 in 2007. (b) Observed minus model differences. (Note the different scale.)

The sea-level observations are in close agreement with tidal predictions (red curve in Fig. 2a) from the Arctic tidal model AOTIM-5 (Padman and Erofeeva, 2004), evaluated for a location in the open ocean just outside Sermilik Fjord (Fig. 1). However, the differences between the observations and model (Fig. 2b) are not negligible. In particular, the differences show harmonic variations with diurnal and semidiurnal periods of maximum amplitude ~0.1 m. Differences between sea-level changes within the fjord and the tide model may be due to atmospheric pressure, winds and runoff within the fjord, all of which can have significant diurnal variability (e.g. Svendsen and Thompson, 1978; Mortensen and others, 2011).

Other differences between the sea-level observations and the tide model can also be observed in Figure 2. The short time spans of high-frequency variations are due to tsunamis produced by calving (Nettles and others, 2008). Longer-period variations can also be observed, particularly an apparent rate change around day 190. The longer-period differences are potentially problematic for such a short time series because they have to be dealt with in any spectral analysis. Although the time span of data collected before day 190 is short, there appears to be a change in the secular drift of the differences at this epoch. We speculate that this may be due to one or more of several possible causes: (1) the pressure-gauge electronics adjusting to environmental conditions; (2) a real long-period variation in sea level caused, for example, by winds and pressure variations; or (3) large wave arrivals associated with the glacial earthquakes on days 189–90 may have caused the instrumental drift to change. Since we are interested in a harmonic analysis rather than longer-period variations, we omit the tide-gauge data before day 190 and refer to the resulting dataset as the ‘truncated’ sea-level time series. The truncated sea-level observations have a slightly smaller root-mean-square difference than the tide model (92 compared to 96 mm).

Geodetic data from Helheim Glacier

We acquired continuous GPS observations (dual-frequency phase and pseudo-range data) at sites near the calving front of Helheim Glacier and rock-based sites near the glacier margins (Nettles and others, 2008). The on-glacier network spanned a distance of up to 37 km behind the calving front. The sites were generally distributed along the glacier’s center flowline, although several sites were established off this line to measure cross-glacier flow variation. All systems employed geodetic-quality, dual-frequency Trimble NetRS or Trimble R7 receivers and Trimble Zephyr antennas. We recorded phase, pseudo-range and other data at sampling intervals of 1–5 s. Our on-glacier networks consisted of 23 sites on the glacier surface in 2006, 16 in 2007 and 23 in 2008.

We used the TRACK software package (Chen, 1998) to perform a fully kinematic analysis of the GPS data. This analysis yielded time-dependent estimates of position every 15 s for the GPS antennas on the glacier relative to the static bedrock antennas. Geocentric Cartesian coordinate time series from each day were transformed into glacier-based coordinate systems by a rotation about the local vertical axis to define ‘along-flow’ and ‘cross-flow’ components. Typical formal uncertainties for the horizontal coordinates are 5–10 mm. We deleted position estimates for epochs at which any ambiguities were unresolved.

In 2007, we acquired data from various Helheim Glacier sites from day 185 (4 July) to 237 (25 August) (Nettles and others, 2008). Site IS22, located a distance of 1.3–2.5 km from Helheim Glacier’s terminus near the center of the glacier (Fig. 1) provides a nearly complete time series for 21 days (186–206) in 2007, a period of significant overlap with the tide-gauge record. The glacier was grounded at or near the calving front during the period of observation, based on the lack of a vertical response at this site.

The time series of estimated along-flow position for IS22 from the GPS analysis, after removal of a mean flow of 22.3 m d^-1, is shown in Figure 3. Days 189 and 190 have been omitted from this plot and from the analysis below, since these time series are complicated by a series of glacial earthquakes (Nettles and others, 2008; de Juan and others, 2010). Day boundaries are shown to emphasize the fact that the GPS data were analyzed independently in 1 day batches (Nettles and others, 2008; de Juan and others, 2010), as is typical for such data.

Fig. 3. Estimated relative along-flow position for GPS site IS22 (Fig. 1). (a) Along-flow position relative to mean flow of 22.3 md⁻¹. (b) Along-flow position relative to a model for long-term variability of position, revealing harmonic variability. Data during the days of glacial earthquakes (189–190) have been omitted.

At the scale of Figure 3a, non-periodic variations in velocity are dominant, and periodic variations are barely visible. To view these, we first remove an approximate model for the long-period variability. The resulting variations in flow are shown in Figure 3b. At this scale, the glacier flow can be seen to contain large diurnal and semidiurnal variations of time-dependent amplitude.

Methods

The goal of our study is to understand the source of the observed harmonic variations of glacier flow for Helheim Glacier, a fast-flowing tidewater glacier. As we discussed in the Introduction, previous studies (e.g. Meier and Post, 1987) have demonstrated that such glaciers respond to tidal forcing. A small number of studies (e.g. Andersen and others, 2010, 2011) have demonstrated that large tidewater glaciers can respond to atmospheric forcing on multi-day timescales. However, it has not been demonstrated that Greenland’s tidewater glaciers respond to diurnal atmospheric forcing. Our approach will therefore be first to explore whether we can fully explain the harmonic variations of Helheim Glacier using a tidal admittance (e.g. O’Neel and others, 2001). We will then use the internal consistency of the results, reflected in the variability of the admittance across the tidal spectrum, to search for the presence of non-tidal components of glacier flow.

Because of the short lengths of the fjord sea-level and glacier-flow time series, it is not possible to determine the complete spectrum of sea-level variations, or to separate closely spaced harmonics uniquely. However, the goal of this study is not to provide such a complete model, but rather to understand the amplitude of variations of the main constituents within the diurnal and semidiurnal bands. Our analysis thus allows for some differentiation among constituents within each band. In addition, we wish to perform the same analysis on both the tide-gauge and position time series, where the latter is dominated by large non-periodic flow-speed variations. Although biases between daily analyses of GPS data of position appear to be small in Figure 3b, we allow for such a possibility in order not to misinterpret such potential biases in terms of harmonic variability. We therefore seek an alternative to the traditional harmonic analysis approach, that will enable simultaneous estimation of the variable flow signal and harmonic components and that will minimize the impact of these biases.

In performing this analysis for the tide-gauge data, we start with components in the harmonic spectrum of AOTIM-5 for the location shown in Figure 1. However, the qualitative analysis based on visual inspection of the differences between the time series and the tide model in the previous section indicated significant differences between sea-level variations in the fjord and those predicted for tides in the open ocean. Given the short time series, we cannot resolve the exact frequency of the differences. The tide model includes components at two nearly diurnal tidal frequencies, i.e. frequencies, f, close to 1 cycle per (solar) day (cpd), but it does not include the S₁ tide (f = 1 cpd) because this tide is very small in the open ocean (e.g. Ray and Egbert, 2004). The P₁ tide has a frequency of 1–1/365: 243 cpd and the K₁ tide has a frequency of 1 + 1/365: 243 cpd However, due to the possibility of environmental forcings in the fjord that are not present, or are much smaller, in the open ocean (e.g. Svendsen and Thompson, 1978; Mortensen and others, 2011), we allow for the possibility of an S₁ component in the fjord sea level.

To estimate components at these close frequencies in the sea-level data, we constrained the ratio of the P₁ and K₁ amplitudes to their AOTIM-5 values. This is similar to the technique of inference (Godin, 1972), except only the real amplitudes were constrained. Due to the short length of the time series, we made no nodal modulation corrections (e.g. Godin, 1972). Long-period tides could not be determined, and we simply remove a best-fit straight line prior to harmonic analysis.

Given the approach described in the first paragraph of this section, we wish to decompose the glacier-flow time series into harmonic series based on tidal constituents, but we do not wish to constrain the amplitudes based on the estimated model for the sea-level observations. A challenge in the analysis of the glacier-flow data arises because, in comparison with the tide-gauge data, the glacier-flow uncertainties are larger relative to the periodic signals (Fig. 3), and there are also data gaps. The problem of separation of various constituents is thus more acute, especially for the P₁, S₁ and K₁ components. We nevertheless must choose specific frequencies to use for estimation. Therefore, we model any combined diurnal signals that cannot be resolved from one another using a frequency of 1 cpd (i.e. S₁).

We used the method described in the Appendix to analyze the sea-level and position time series consistently. The initial step of this approach involves estimating, and removing, amplitudes of signals with nominal frequencies within each band for daily batches of data. (We choose 1 and 2 cpd for simplicity.) This procedure has the effect of shifting downward the frequencies of the components making up the original time series. For example, a signal with the M₂ frequency when shifted by 2 cpd has a frequency of –0.0677 cpd (period of –14: 764 days), which can be easily distinguished from a signal with the shifted N₂ frequency of –0.104 cpd (period 9.614 days). The S₂ component itself shows up in the shifted signal as a zero-frequency offset. The amplitude of each identified harmonic will potentially represent the summation over dozens or hundreds of tidal harmonics that are indistinguishable over so short a time span. We compared results for our code against T_TIDE (Foreman, 1977; Pawlowicz and others, 2002) using simulated time series, and found good agreement.

Results

Power spectral densities

The estimated power spectral density (PSD) for each observed time series is shown in Figure 4. For the sea-level data, a straightforward Fourier-transform method was used on the truncated time series. For the position time series, we first limited ourselves to the nearly linear span of data from days 191–200. Even during this period there are daily speed variations of 0.1–0.3 md⁻¹, so we first removed a model consisting of best-fit daily piecewise continuous straight lines. To account for unequal data spacing (due to data deletions) and unequal standard deviations, we used the Lomb–Scargle method (Lomb, 1976; Scargle, 1982) to calculate the PSD.

Fig. 4. Power spectral densities of the observed time series, with initial identifications of the spectral peaks. (a) Truncated sea-level series. (b) Power spectral density of along-flow position from GPS, after subtraction of a best-fit model consisting of daily piecewise-continuous straight lines. Diurnal (D) and semidiurnal (S) peaks are labelled. The PSDs are normalized so that integration across the band shown in the plot yields the sample variance.

For each spectrum in Figure 4, an initial identification of spectral peaks is shown. For the sea-level data, the PSD reveals four main peaks, which we nominally identify using Darwin notation. The S₁, K₁ and P₁ components cannot be separated visually in this plot, and the diurnal peak has been labelled ‘K₁/P₁/S₁’. Peaks at 3, 4 and 5 cpd are also visible, which may be present in the open-ocean tide or produced by internal reflections within the fjord (e.g. Valle-Levinson and others, 2007). Higher-frequency power associated with calving events and, presumably, data noise is also present in the PSD.

For the glacier-flow time series, peaks at diurnal, semidiurnal and other diurnal harmonics are visible. Due to the relatively poor resolution of the flow PSD compared with the sea-level PSD (due to larger errors and data gaps), we have labelled the peaks simply ‘D’ for diurnal and ‘S’ for semidiurnal. In contrast to the tide-gauge spectrum, the glacier-flow spectrum has more energy in the diurnal frequency band than in the semidiurnal.

Harmonic constituents

The results of the harmonic analysis of the sea-level data (Table 1, columns labelled ‘Sea level’) are consistent with the visual inspections of the time series and its PSD, in indicating that sea-level variation is dominated by semi-diurnal tides. The primary harmonic is M₂, with an amplitude of ~1.0 m. The next largest component is S₂, with an amplitude only 34% of the M₂. The root-sum-square (rss) of all the diurnal amplitudes (labelled ‘D’ in Table 1) is 18% of that for the semidiurnal (SD) amplitudes, a result consistent with the visual observation of the primarily semidiurnal nature of the ocean tide (Fig. 2).

Table 1. Results of the analysis of the truncated sea-level and glacier-flow time series. Rows labelled ‘D’ (diurnal) and ‘SD’ (semidiurnal) give the root-sum-square of the amplitudes for the respective band

The S₁ frequency we found for the sea level is too large to be reasonable for the S₁ tide in the open ocean (e.g. Ray and Egbert, 2004), but we cannot say with certainty whether this is a real sea-level variation that occurs within the fjord or whether this represents entanglement with K₁ and P₁. The S₁ ocean tide can reach the centimeter level near Helheim Glacier (Ponchaut and others, 2001). However, we are not so much concerned with the details of the tidal model in the fjord as the comparison of energy in the diurnal and semidiurnal bands.

The analysis of the glacier-flow time series (Table 1, columns labelled ‘Glacier flow’) produced fewer components, because of the larger uncertainties, relative to the signal, of these data. Two relatively large diurnal components were estimated, although we lack the ability to resolve the closely spaced nearly diurnal signals (see Methods section). In contrast to the sea-level spectrum, the diurnal amplitude for the flow spectrum is larger than the semidiurnal by a factor of ~1.5. We therefore suspect that these components of glacier flow may be primarily non-tidal in origin.

The residuals for both the sea-level and glacier-flow time series, relative to their respective models composed of the harmonic components of Table 1 (along with the long-term trend variations), are shown in Figure 5. The rms residual for the sea-level data is 42 mm, indicating, as might be expected for an analysis based on a short record, a much closer agreement of this model with the sea-level observations than achieved by the open-ocean tide model. The residuals do not appear to contain any harmonic variability that is coherent over the time span of the data. We did not estimate any higher-order harmonics, so those would still be present in the residuals. However, there is no obvious pattern in the residuals that would modify our conclusion that the sea-level variability in the fjord is dominated by the semidiurnal components of the ocean tide. We note that, due to the short time span of the tide-gauge time series, these results do not represent a complete tidal model in the sense of, for example, AOTIM.

Fig. 5. Residuals relative to the best-fit sea-level and glacier-flow models. (a) Tide-gauge residuals. The rms residual is 42 mm. (b) Glacier-flow residuals. The rms residual is 20 mm.

The rms residual for the glacier flow (Fig. 5b) is 20 mm. While this fit is better than that for the sea-level data in absolute terms, relative to the amplitude of the signal it is worse. The signal-to-noise ratio is ~25 for the sea-level data, ~2.8 for the glacier-flow data. These results are not surprising given the much larger uncertainties of the position estimates compared with the nominal tide-gauge uncertainties. However, the residuals from both time series show systematic variations. For sea level, ocean-dynamic signals play an important role. For the glacier-flow data, the systematic variations at short periods are more pronounced, and the residuals (Fig. 5) are on the same order as the estimated amplitudes (Table 1). The model fit is especially poor around day 203.

Analysis

One way to compare the spectra of the ocean tide and the periodic glacier flow is to model the glacier-flow variations for the k th harmonic component, x_k (t), as a lagged linear response to sea-level forcing, F_k :

(1)

where A_k is the admittance and τ_k the time delay for the k th harmonic component (Walters and Dunlap, 1987; O’Neel and others, 2001). Our estimates of the admittances and delays are shown in the last two columns of Table 1. We calculated the admittance by dividing the along-track flow amplitude by the ocean-tide amplitude. A negative sign for the admittance indicates that a positive sea-level change yields a negative flow displacement. The lag, τ_k , for the k th harmonic component is calculated from the phase (°) using

(2)

where is the tide-gauge phase, is the along-track flow phase, T_k is the period (hours) and the quantity in parentheses is taken modulo 360°. The 180° phase advance takes into account the fact that, for sea-level back-pressure forcing, the along-track position minima occur near forcing maxima (e.g. de Juan and others, 2010).

Can we assume a frequency-independent response of the glacier back-pressure forcing? In many models, glacier speed is proportional to the driving stress, τⁿ , with n ≃ 3 (e.g. Walters, 1989; Paterson, 1994; Howat and others, 2005). For ice streams with small mean driving stresses and large tidally imparted hydrostatic stresses, this nonlinearity yields tidal flow responses that are unequal across tidal frequencies (e.g. Gudmundsson, 2007; King and others, 2010). However, in glaciers for which the tidally imparted stress is small compared with the mean driving stress (e.g. Helheim Glacier), these models predict a linear response to the tide. Using the model described by Gudmundsson (2007), we find that the term that has a linear dependence with tidal height is at least two orders of magnitude larger than higher-order, nonlinear terms.

The mean semidiurnal admittance (Table 1) is –0. 028 ± 0. 002. Applying this admittance to the diurnal components, it is not surprising that no position spectral components for P₁, K₁, K₂, ν ₂ and O₁ appear, since their amplitudes in the tide-gauge spectrum are too small. The admittances for the M₂, S₂ and N₂ components are consistent at better than the 2σ level, whereas the admittance for the S₁ component is a factor of ~27 larger, and that for the Q₁ component is ~36 times larger. The phases of the position spectra reveal a time lag of position relative to the tide-gauge observations. Again, the M₂, S₂ and N₂ lags are consistent, and within 2σ of the weighted mean of 1.9 hours. However, the S₁ lag is 21.7 hours and the Q₁ lag is 11 hours.

Figure 6 demonstrates that the three primary observed semidiurnal components (M₂, S₂ and N₂) of the along-track glacier motion at site IS22 can be accurately predicted from the tide-gauge observations, using a single admittance and response lag. The two primary diurnal components (S₁ and Q₁), however, cannot be predicted from these same admittance and lag values; in fact they cannot be predicted from any single pair of admittance and lag. As we argue above, the response of Helheim Glacier to sea-level back-pressure forcing should vary little, if at all, between diurnal and semidiurnal bands, in which case all the admittances and lags should be equal, whereas our results require a large difference in the responses between the diurnal and semidiurnal bands.

Fig. 6. Phasor diagram of the IS22 position spectral components (blue), along with predictions (red) based on the tide-gauge spectra using the mean semidiurnal admittance and lag from Table 1. The 95% confidence ellipses are shown. The P₁, K₁, K₂, O₁ and ν ₂ components are not shown, since no corresponding component for the glacier position was estimated. The Q₁ tide-gauge estimate is not shown since it is nearly zero. The predicted M₂, N₂ and S₂ components lie within the 95% confidence ellipses for their respective observed values, whereas those for the S₁ and Q₁ components do not. Phase angle is measured clockwise from the horizontal x –axis (not shown).

Our analyses of the sea-level forcing and the glacier response lead us to conclude that the diurnal-band glacier motion is largely not a response to back pressure caused by sea-level variation. Based on previous studies of other types of glacier (e.g. Walters and Dunlap, 1987; Kamb and others, 1994; O’Neel and others, 2001) we hypothesize that the diurnal variations in glacier flow are associated with variations in lubrication of the glacier bed, due to surface melting driven by diurnal solar heating. In this scenario, the estimated Q₁ harmonic is an artifact that, for our data span, absorbs some of the (possibly weather-dependent) temporal variability of the diurnal signal.

Although separation of tidal and non-tidal components of forcing has previously been challenging (e.g. Walter and others, 2012), we can use the admittances for this purpose. We first formed a residual diurnal position spectrum (as described by O’Neel and others, 2001). We predict the diurnal spectrum of along-flow position using the diurnal sea-level spectrum and mean admittance and lag from Table 1, for the set of diurnal tides estimated for position, also given in Table 1. We then subtract the predicted diurnal spectrum from the observed along-flow diurnal spectrum to obtain the non-tidal variations. One-day model time series for position and velocity formed from the residual non-tidal diurnal components are shown in Figure 7. This signal consists of two components: S₁ (with an exactly 1 day period) and Q₁, with which the mixing produces the 9.4 day beating period of Figure 7. The ~60 mm position amplitude of this signal produces, in velocity, a speed change of ~0.7 m d^-1 amplitude peak-to-peak, or a peak-to-peak variation of ~3% compared with the glacier mean flow speed of ~22 m d^-1 at IS22 during this period. From Figure 7, the maximum speed occurs ~6 hours after local noon during this period and varies slightly over a longer time period because of the Q₁ modulation. (The insets of Fig. 7 show longer time spans.)

Fig. 7. Model for the non-tidal diurnal signal in IS22 along-flow position and speed predicted from the residual glacier-flow spectrum, for day 195. The solid gray line indicates local noon based on the longitude of Helheim Glacier. Insets: model time series for days 190–210.

In their investigation of speed variations at Columbia Glacier, Alaska, Walters and Dunlap (1987) found that a component of the variations was correlated with solar insolation, which they used as a proxy for melt runoff. They found a lag of ~8 hours between the glacier speed response and the insolation. We find that the maximum glacier speed from the non-tidal signal occurs ~6 hours past local noon; a good agreement with Walters and Dunlap (1987), if we assume that solar altitude should also be a good proxy for melt. However, Walters and Dunlap (1987) found their melt signal to be temporally variable, a reasonable finding given that insolation at the surface is impacted by weather. The analysis presented above, which was focused on the tidal spectrum, estimated a single S₁ amplitude and phase for the entire position time series. Visual inspection of the postfit residuals (Fig. 5), however, leads us to infer that a temporally variable diurnal signal might be more appropriate. The postfit residuals have some periodicity. In the earlier half of the time series, there is a clear residual diurnal variation, and the model is unable to reproduce the variability during some days (e.g. days 201–203).

To investigate the variability of the diurnal component further, we designed a Kalman filter (e.g. Ravishanker and Dey, 2002) to estimate the parameters of an admittance model with a stochastic (random-walk) along-flow speed, tidal admittance, tidal-response lag and amplitude of the diurnal sinusoidal variability. The model for the along-flow position, x (t), is given by

(3)

Here x˳ = x(t ˳)is the position at initial epoch t˳ v (t) is the non-periodic along-flow speed, A (t) is the tidal admittance, F (t) is the tide-gauge observation, τ (t) is the lag, (t) is assumed to be zero-mean Gaussian random error and the diurnal signal, x _D(t), is given by

(4)

This form was used to insure continuity across day boundaries, at which epochs the contribution of x _D(t) is identically zero. The stochastic amplitudes, a _c(t) and a _s(t), were modeled as slowly varying random-walk stochastic processes with variance rates of 10^-4 m² d^-1. For this analysis, only data from the period when the GPS observations overlapped with the tide-gauge observations were used, and observations during the glacial earthquakes were retained to provide continuity.

The results of the Kalman-filter solution are shown in Figure 8. Figure 8a shows the diurnal contribution of Eqn (4). During this time span, the amplitude of the diurnal signal varies by a factor of ~2. The most rapid change is seen in days 189–190, during the period of glacial earthquake activity. De Juan and others (2010) reported an abrupt change in the tidal admittance during that time, that we also find. The postfit residuals shown in Figure 8b are significantly smaller than those calculated using the model with a constant diurnal amplitude (Fig. 3b), having an rms value less than half that of the previous solution.

Fig. 8. Results from the Kalman-filter solution with stochastic S₁ amplitude. (a) Estimated S₁ contribution (x _D(t) from Eqn (4)). (b) Postfit residual. The rms residual is 9.7 mm. Error bars (1) are shown in gray.

The Kalman-filter solution explicitly includes only those harmonic components inherent in the tide-gauge data (and having the same relative amplitude as in the tide-gauge data) plus an additional time-variable diurnal component. The time variability of a diurnal component will express itself (in power-spectrum terms) as a broadening of the peak centered on the diurnal frequency. The ability of the stochastic model to fit this variability reinforces the supposition, discussed above, that the large Q₁ component was simply an artifact of that variability over this short time series.

In future studies, we plan to use meteorological data collected on Helheim Glacier to improve our understanding of the influence of melt on glacier flow. We have previously been able to model the surface energy balance with some accuracy (Andersen and others, 2010, 2011), and comparison with a more complete analysis of speed variations throughout the network should be fruitful. Although such a study is beyond the scope of this paper, we can perform a simple comparison of the model for speed variations inferred from the Kalman-filter model for the diurnal contribution with measurements of temperature obtained at the nearby weather station at Tasiilaq (Fig. 9). The speed and temperature variations are clearly correlated, with the speed generally lagging the temperature by ~3.8 hours. The temporal variability of the amplitude response shows a more complex relation, as expected from previous work at Columbia Glacier (e.g. Kamb and others, 1994).

Fig. 9. Comparison of estimated diurnal speed variations for IS22 from the Kalman-filter (blue) and surface air temperature anomalies from the nearby weather station at Tasiilaq (red) relative to a Gaussian smoothed value calculated with a 1.5 hour window. Data from Cappelen (2014).

We have focused on GPS site IS22 for this study because it provides a complete, high-quality record of flow variations during a period in which tide-gauge observations were also available, and the goal of this study was to identify and quantify the non-tidal component of glacier flow. Interpretation of our previous analyses (de Juan and others, 2010), in light of the results presented here, leads us to expect that future study of a broader range of glacier sites and of longer time periods may prove illuminating.

Conclusions

We compared the spectrum of sea-level time series acquired in Sermilik Fjord near Helheim Glacier with the spectrum of along-flow position time series acquired from a GPS site on the glacier. We found that the position spectrum was poorly explained by the assumption of forcing by the sea level alone, and that there were additional components to the motion in the diurnal band, centered around the S₁ frequency of 1 cpd The non-tidal variations accounted for a peak-to-peak variation in glacier flow speed of ~0.7 m d^-1 at a location where the mean flow rate was ~22 m d^-1. The phase of the non-tidal variations was such that the maximum flow speed was reached ~6 hours after local noon, suggesting that the source of these variations was surface runoff from ice melt associated with solar heating. Further investigation of the non-tidal diurnal signal using a stochastic filter revealed that it could be better modeled using stochastic amplitudes (sine and cosine) that varied by a factor of ~2 over the ~21 day dataset.

Our analysis focused on rigorously establishing the presence of non-tidal diurnal variations in glacier speed. It is reasonable to conclude that the non-tidal signal is associated with melt runoff, both from the point of view of the above discussion and what is known regarding the influence of runoff on glacier flow. At periods exceeding 1 day, data from Helheim Glacier have already been used to investigate the relationship between flow-speed variations and melt (Andersen and others, 2010, 2011). However, unambiguous identification of a non-tidal melt component at timescales of ≤ 1 day is not straightforward. Any melt component associated with solar heating will be represented mainly at a frequency of 1 cpd and harmonics thereof, which are also valid tidal frequencies. Moreover, the melt component will maintain phase with these tides, since the Sun is the main contributor to both forcings. In our case, we were able to compare the spectrum of glacier-flow determinations from GPS to the tide spectrum from tide-gauge observations. Our approach also depends on having highly accurate continuous time series, such as those afforded by geodesy with GPS.

Our comparison of the spectra of the tidal forcing and glacier flow yielded tidal admittances in good agreement with the admittances calculated by de Juan and others (2010) using a purely time-domain modeling approach. This agreement is due both to the large semidiurnal component of the tidal forcing and also to the inclusion by de Juan and others (2010) of an additional acceleration parameter that absorbed much of the variable diurnal signal reported here. Our work suggests that a comparison of the spectra of the various forcings and the glacier flow are important to understanding their individual contributions.