2.1. Ice-stream (response) model

Reference GudmundssonGudmundsson (2007) modelled the total forward surface velocity, u _s, of Rutford Ice Stream as

where r is the ratio between the mean sliding velocity and the mean forward deformational velocity, is the mean basal shear stress, K is a site-specific constant of proportionality between the tidally related hydrostatic pressure variation, ρ _w gh(t), and the resulting perturbation in basal shear stress, ρ _w is the density of sea water (1023 kg m⁻³), g is acceleration due to gravity (9.81 m s⁻²), C is the sliding coefficient and m is the power-law exponent. Reference GudmundssonGudmundsson (2007) fit model parameters as listed in Table 1 (solution A), but noted that this set of model parameters is not unique in providing a satisfactory fit to the data segment analysed. Reference GudmundssonGudmundsson (2007) provided a detailed discussion of physically plausible values of , K and m in Equation (1), which we do not repeat here.

Table 1.

Model parameters: A from Reference GudmundssonGudmundsson (2007), the others fitted to data from Reference Murray, Smith, King and WeedonMurray and others (2007); see text for details. Formal errors on the parameters determined using least squares (solutions C, D and E) are <1 × 10⁻⁴. Also shown are the mean velocity and the contribution of the tides to mean velocity

The time-variable component of the modelled motion is controlled purely by the tidal height, h(t), at a given time, t, at the grounding line. Rounding the value of m = 3.04 found by Reference GudmundssonGudmundsson (2007) to exactly 3 (consistent with, for example, Reference PatersonPaterson, 1994) allows us to investigate the frequency response given by the model. We show below that this value of m is also in agreement with a much longer data time series. We expand the relevant part of Equation (1) to obtain:

The time-variable response in this model is therefore governed by powers 1–3 of h(t). Consider now pairs (i and j) of tidal constituents such that:

where A, ω = 2π/T and φ are the constituent-specific amplitude, angular frequency and Greenwich phase, respectively, for tidal constituents with period T. Simplifying further to study the frequency content of the response only, we take A_i = A_j = 1 and φ_i = φ_j = 0. After manipulation:

and

In the form of Equations (4) and (5), the modelled frequencies of response, hidden in Equation (2), may be read directly. We examine the frequencies in Equations (4) and (5) more closely in section 3.

2.2. Tidal (forcing) model

We initially force Equation (1) with a recent regional Antarctic tide model, CATS2008a, an updated version of the CATS02.01 model described by Reference Padman, Fricker, Coleman, Howard and ErofeevaPadman and others (2002). CATS02.01 predictions were used by Reference GudmundssonGudmundsson (2007) to force Equation (1). At the time of writing, CATS2008a is the most accurate ocean-tide model for the circum-Antarctic seas, with per-constituent error generally <0.05 m. Included in CATS2008a are four major semi-diurnal (M₂, S₂, K₂, N₂) and four diurnal (K₁, O₁, P₁, Q₁) tidal constituents (see Table 2 for frequencies). From these we infer many other related tidal constituents in these bands in a standard way (Reference PughPugh, 1987). We do not include forcing from longer-period signals, such as fortnightly terms, since the ice stream does not respond elastically at these frequencies, while it does perfectly so in the semi-diurnal band and very nearly perfectly in the diurnal band (Reference GudmundssonGudmundsson, 2007). Adopting an ocean-tide model means that we can make tidal predictions very close to the grounding zone, but any model errors will propagate into the predicted ice-stream response. To complement this approach, we also force Equation (1) with tidal predictions based on GPS observations of vertical tidal motion ∼20 km downstream of the grounding line (Reference GudmundssonGudmundsson, 2007). To do this, we use the analysis of Reference King and PadmanKing and Padman (2005) who separated related constituents in this short record (P₁/K₁ and S₂/K₂) using phase and amplitude relationships in CATS02.01 and we make tidal predictions based on those constituents, inferring related constituents as above. There will be some error in not accounting for the difference in tides at this location compared with those in the grounding zone. There is also error in CATS2008a; at the ice-shelf GPS site, differences between the CATS02.01/CATS2008a and observed tidal constituents reach several centimetres (Reference King and PadmanKing and Padman, 2005).

Table 2.

Response periods (days) for the major semi-diurnal and diurnal constituent pairs with m = 3

Owing to the non-linear nature of the model, the low-frequency terms are particularly sensitive to changes in the predicted amplitude and phase of the vertical semi-diurnal and diurnal tides. To examine the sensitivity to tidal errors, we modified the CATS2008a predicted M₂ tide by 0.1 m in amplitude and evaluated Equation (1) (using the Reference GudmundssonGudmundsson (2007) values) and compared the output to results using an unmodified CATS2008a prediction. This level of vertical tidal error produced only an 8 mm difference in modelled horizontal displacement amplitude at 14.76 days and only a 3 mm difference in the semi-diurnal band. The same error applied to K2 produced a 50 mm difference at 183 days and a 1 mm difference in the semi-diurnal band. This suggests the modelled displacements are largely invariant to reasonable ocean-tide prediction errors, except at the longest periods where tide-related errors may be substantial.

2.3. Ice-stream position data

Our comparison of model and data is based on the ∼2 year dataset of Reference Murray, Smith, King and WeedonMurray and others (2007). However, there is a large data gap during the first ∼0.5 years (Reference Murray, Smith, King and WeedonMurray and others, 2007), so we excluded the entire earlier part of the time series. This yielded a highly periodic ice-stream response as shown in Figure 1 (after detrending the along-flow position) of ∼1.5 years in duration (although with some significant gaps). This is substantially longer than the similarly located records of Reference GudmundssonGudmundsson (2006) and Reference AðalgeirsdóttirAðalgeirsdóttir and others (2008), both of which were ∼7 weeks in duration. The record of Reference AðalgeirsdóttirAðalgeirsdóttir and others (2008) forms the initial portion of the record of Reference Murray, Smith, King and WeedonMurray and others (2007) and hence is not included in our analysis. In using this long record, we sacrifice some time-series precision and accuracy due to necessarily different GPS data analysis techniques (the shorter records were derived relative to a local base station, with some GPS systematic errors differencing, whereas the longer record we use was derived using Precise Point Positioning where they do not difference). Most notably, the record of Reference Murray, Smith, King and WeedonMurray and others (2007) has larger GPS-related systematic bias at some tidal periods (e.g. K₁, K₂ and S₂; Reference King, Watson, Penna and ClarkeKing and others, 2008). Here we typically present detrended along-flow position in figures, but the displacement due to mean flow is also included in the model tuning and analysis.

Fig. 1.

Observed and modelled response (after detrending) using fits to the displacement data. Data are shown in terms of (a) velocity and (b) residual to mean along-flow displacement. Velocities with periods ≤1 day have been smoothed (a). Date format: Oct05 = October 2005.

3.1. Fitting to displacement

For the displacement analysis we followed the procedure of Reference GudmundssonGudmundsson (2007) by first forward-integrating the modelled velocities to produce modelled displacements, meaning that the model tuning is dominated by fortnightly and semi-annual signals, plus the mean velocity. We tuned the model by minimizing the RMS of the observed minus-modelled detrended displacement signal on a parameter-byparameter basis. We tested a wide range of values around the values of Reference GudmundssonGudmundsson (2007), modifying one parameter while holding the others fixed. However, we note that r and C are mutually dependent in this model and hence cannot be determined independently. The determined model parameters are given in Table 1 (solution B). We found the model parameters to be well defined at the level of precision given in Table 1.

The model predictions of velocity based on these parameter values are shown in Figure 1, along with detrended displacement computed using forward integration. The velocity signal is shown after filtering out signals with periods ≤1 day. These are shown alongside the predictions using the parameter values from Reference GudmundssonGudmundsson (2007). The values of Reference GudmundssonGudmundsson (2007) clearly are not appropriate for our site, as they predict a higher mean velocity; otherwise the output is very similar. Regardless of model parameters, the amplitude of the semi-diurnal and diurnal signals is clearly overestimated. As we expect, given that mean flow, semi-annual and fortnightly signals dominate the displacement time series, the fit for these frequencies is better than for the semi-diurnal and diurnal frequencies.

3.2. Fitting to velocity

We then attempted to fit the model to the velocity data using least squares. In contrast to the displacement fit, fitting to velocity will emphasize the higher-frequency signals. Since the highest-frequency terms are subject to unmodelled physics as well as GPS systematic error at some frequencies, we filtered these out by computing velocities based on the GPS positions over 2.5 days and also filtering the modelled velocities. To enable linearization of Equation (1), we took the logarithm of both sides, meaning that we actually fit to log₁₀(u_s). In the least-squares solution we fixed r = 107 (value determined by Reference GudmundssonGudmundsson, 2007) to overcome the mutuality with C but estimated the other parameters together. Using r = 110 instead (Table 1, solution B) will result in a slightly lower value of C. To help stabilize the solutions we added non-negativity constraints to , C and m. We found that allowing m to be completely freely varying produced unstable results, so we constrained m = 3 to within ±0.05, but noted that it adjusted by <0.001, suggesting that this constraint did not distort the solution. The solution converged quickly to the values shown in Table 1 (solution C). We repeated the solution on a shorter (∼7 week) section of data and found almost identical values, suggesting that the solution is not dependent on the frequency of the signal, at least at long periods. The solution showed no preference for the sign of K and, when negative initial values of K were specified, it converged to K = −0.18, instead of K = +0.18, with the other parameters being identical.

Comparing the determined model parameters in Table 1, we note that fitting to the velocity gives a ∼20% lower value for and a 100% increase in the value of C, when compared with the solution tuned to displacement. The value of K changes only marginally. Comparing Figures 1 and 3 shows that solution C provides a substantially better fit to the observations. The 14.76 day signal is particularly well replicated, with only very small error in mean velocity (Table 1).

To test the influence of tide model prediction errors in CATS2008a, we forced Equation (1) with the tidal prediction based on the downstream GPS observations rather than CATS2008a. The differences in best-fit parameters are very small (Table 1 solutions C and D). The modelled velocity and displacement for solutions C and D are shown in Figure 3, with the prediction based on CATS2008a giving too large an amplitude for the fortnightly and semi-annual signals (Table 2).

Fig. 3.

Same as Figure 1 but using the fits to the velocity data.

We then repeated the solution after first constraining m = 2. This solution did not converge. Table 1 shows the parameter values after 20 iterations (solution E), and the modelled velocity and displacements are shown in Figure 3. It is evident that the solution with m = 2 cannot replicate the observed long-period signal amplitude (Fig. 2) at the same time as replicating the mean flow rate (Table 1). However, the m = 3 solution does provide a good fit to both. Furthermore, Figure 2 shows that the higher-frequency signal seen in the observations at 6 cycles d⁻¹ (0.167 days; Table 2) is not evident with m = 2. We note signal in the observations also at 5, 7 and 8 cycles d⁻¹, but these are consistent with GPS observation error (Reference King, Watson, Penna and ClarkeKing and others, 2008). Following this analysis, we focus on the m = 3 solutions forced by the GPSderived tidal predictions for the remainder of this paper.

To investigate more closely the content of the observation-model fit across the entire time series, we computed Lomb amplitude spectra which are suited to gappy data (Reference ScargleScargle, 1982; Reference Press, Teukolsky, Vetterling and FlanneryPress and others, 1992). Figure 4 shows the amplitude spectra for the observed displacements and for those that remain after removing the model prediction. We note that the dominant fortnightly and semi-annual terms are well fit by the model. At the highest frequencies the model amplitude becomes obviously too large. The phases of the semi-diurnal and diurnal terms are particularly sensitive to model parameters (not shown). The amplitude over-prediction at semi-diurnal and diurnal amplitudes can be seen also in Reference GudmundssonGudmundsson (2007), using model parameters heavily influenced by the fortnightly term (and unaffected by the semi-annual terms). Reasons for the poor quality of the model fit at the higher frequencies are discussed below (see section 4).

Fig. 4.

Lomb amplitude spectra for the observations together with the residual after removing the model output (solution D). (a) Low-frequency signal and (b) high-frequency signal.

The semi-annual signal observed by Reference Murray, Smith, King and WeedonMurray and others (2007), and now shown to be reproducible by the model of Reference GudmundssonGudmundsson (2007), is not the longest-period signal that may be predicted by Equation (1). The lunar-related constituents (all considered here, except S2) undergo modulation throughout the 18.6 year lunar node tide cycle (Reference PughPugh, 1987). To examine the effect of this decadal change in the sub-daily term, we forced the model over 1.5 cycles (28 years from 1990.0). In addition to the signals evident in Figure 1, a signal with period 18.6 years and displacement amplitude ∼0.5 m is predicted by the model. Therefore, the effect on velocity is negligible and is unlikely to be of consequence.

4.1. Model fit at diurnal and semi-diurnal periods

At present, the relatively poor fit between the model and the diurnal and semi-diurnal flow variations is not explained fully. We can, however, speculate on three possible limitations of the model which could be significant (see also the discussion of Reference GudmundssonGudmundsson, 2007). Further data would be required to address these possibilities.

4.2. Model predictions for other ice streams

Reference GudmundssonGudmundsson (2007) considered the implications of the ice-stream response model for other ice streams, both those where tidal heights at the grounding line are similar to Rutford Ice Stream and those where they are significantly different. He concluded that different dynamic responses were mostly due to variations in tidal forcing, rather than any fundamental differences between the ice streams. Following our success applying the same model to our ∼1.5 year time series from Rutford Ice Stream, we now use it to make predictions for two other ice streams. This allows us to predict dynamic behaviour, which can be tested with field observations, and provide a further assessment of the applicability of the model. Figure 5 shows the predicted response using model parameters from Reference GudmundssonGudmundsson (2007). We used local tidal predictions from CATS2008a. The Rutford Ice Stream values are shown for comparison. We chose these two ice streams because we consider it possible that field observations to compare with the model predictions could be obtained in the near future.

Fig. 5.

Predictions of response by period (a, c) and residual to mean along-flow displacement (b, d) for three ice streams forced using the semi-diurnal constituents only (a, b) and both the semi-diurnal and diurnal constituents (c, d).

The semi-diurnal tides near Bindschadler Ice Stream and Lambert Glacier grounding lines are much smaller than those near Rutford Ice Stream and hence the predicted response is 2–10% of the Rutford response at semi-annual periods, as can be seen in Figure 5a and b. The response due to the diurnal terms is, however, similar. This is in disagreement with Reference GudmundssonGudmundsson (2007) who, after examining only semi-diurnal forcing of Bindschadler Ice Stream, suggested that a diurnal signal would dominate any long-period response (ignoring any frequency dependence in the elastic response of the ice stream). Instead, Figure 5c and d show that, assuming identical model parameters apply, tidal modulation of flow at semi-annual and fortnightly periods may also be observable at these other locations, with the amplitude of variations at along-flow positions greatest at lower frequencies (the converse will be true in terms of along-flow velocity). Since the signal-to-noise ratio of the long-period response is greatly increased relative to that for diurnal and semi-diurnal modulations, this could assist attempts to understand ice-stream/till interactions. Table 2 gives the periods of possible modulation that need to be considered. This will, of course, be subject to the model and its parameters being appropriate to these regions also. Tests reveal that amplitudes of variations in flow rate scale with the magnitude of the basal shear stress. This could be tested easily with appropriate GPS datasets, although long records of precise ice motion have not yet been published for either location.

4.3. Tidal component of mean ice flow

A necessary result of the non-linear response to tidally varying basal shear stress is that some component of the mean ice-stream velocity is present only because of the ocean tides (Reference GudmundssonGudmundsson, 2007). The contribution of tides to mean flow is given in Table 1 for the various solutions, showing a substantial increase (from 3.6% to 6.2% of mean flow) using our revised parameter values. This time-averaged component of the mean flow is almost entirely contributed by the (h(t))² term in Equation (2), with a near-negligible component provided by the fourth term. The (h(t))² term contributes to mean flow, according to this model, since it always possesses a positive mean, regardless of the sign of h(t) or K. The tidal component of mean flow in this model will therefore scale with the square of the tidal amplitude and is not sensitive to the sign of K.

Interestingly, this suggests that if the tidal range was to change for some reason, then the ice stream would also undergo a change in its mean velocity. Such a change in ocean-tide range may have occurred around Antarctica during the Holocene, for example (Reference Arbic, MacAyeal, Mitrovica and MilneArbic and others, 2004, Reference Arbic, Mitrovica, MacAyeal and Milne2008; Reference Griffiths and PeltierGriffiths and Peltier, 2009). The non-linear response of the model means that larger tidal range will give much faster mean ice-stream flow and smaller tidal range will give slightly slower flow (zero tidal range gives ∼6% slower flow in the case of the present-day Rutford Ice Stream). One implication of this is the potential influence on ice flow of sub-ice cavity shape near the grounding line. A restricted cavity normally amplifies the tidal range, although it is possible for damping to occur under some conditions. Falling sea level, for example, at a location where this reduced the cavity such that tidal amplification was significantly increased, would theoretically result in an increase in ice flow.