Ice-stream response to ocean tides and the form of the basal sliding law

The response of ice streams to ocean tides is investigated. Numerical modelling experiments are conducted using a two-dimensional flow-line model of coupled ice-stream and ice-shelf flow. The model includes all components of the equilibrium equations, and uses a non-linear viscoelastic constitutive equation for ice. Basal sliding is simulated with a Weertman type sliding law where basal sliding is proportional to some power of the basal shear stress. The response of ice-streams to tidal forcing is found to be profoundly affected by mechanical conditions at the bed. For a non-linear sliding law, a non-linear interaction between the two main semi-diurnal tidal constituents (M2 and S2) can give rise to a significant perturbation in ice-stream flow at the lunisolar synodic fortnightly (MSf) tidal period of 14.76 days. For a linear sliding law, in contrast, no such modulation in flow at the MSf frequency is found. For vertical ocean tides of the type observed on Filchner-Ronne Ice Shelf (FRIS), the amplitude of the horizontal modulation in ice-stream flow at the MSf frequency resulting from a non-linear interaction between the S2 and M2 tidal constitutes can be larger than the direct response at the S2 and the M2 frequencies. In comparison the non-linear interaction between K1 and O1 tidal components is weak. As a consequence, modelled icestream response to mixed oceanic tides of the type found on FRIS is stronger at the MSf period of 14.76 days than at both the semi-diurnal and diurnal frequencies, while at the same time almost absent at the similar Mf period of 13.66 days. The model results compare favourably with measurements of tidally induced flow variations on Rutford Ice Stream (RIS), West Antarctica. On RIS a strong tidal response is found at the MSf frequency with a smaller response at the semidiurnal and diurnal frequencies, and almost no response at the Mf frequency. A non-linear viscous sliding law appears to have the potential to fully explain these observations. Correspondence to: G. H. Gudmundsson (ghg@bas.ac.uk)


Introduction
Observations on a number of ice streams have shown their flow to react sensitively to ocean tides (Anandakrishnan et al., 2003;Bindschadler et al., 2003a,b;Gudmundsson, 2006;Murray et al., 2007;Wiens et al., 2008). In some instances, tidally induced variations in horizontal speed have been detected tens of kilometres upstream from the grounding line causing a periodic variation in flow speeds of up to 20%, depending on location (Gudmundsson, 2006). These variations are therefore neither small nor limited to the zone of elastic flexure around the grounding line. Significant variations in horizontal speed due to tidal action have also been observed downstream of the grounding line on floating ice shelves (Doake et al., 2002;Brunt et al., 2010;King et al., 2011a).
These observations are interesting and intriguing for a number of reasons. They demonstrate that changes in stresses downstream of the grounding line can have a significant and immediate effect on the large-scale flow regime of active ice streams. They also challenge our ability to theoretically characterise ice-stream dynamics and have the potential to provide an insight into the mechanical interaction between ocean, ice shelves, and ice streams. Previous modelling work of tide-induced lateral movement on ice streams have invoked specific assumptions about till rheology (e.g. Bindschadler et al., 2003a;Gudmundsson, 2007;Winberry et al., 2009;Sergienko et al., 2009;King et al., 2010) and ice-stream response to tides can be thought of as a natural experiment providing insight into the mechanics of ice-stream flow.
Rutford Ice Stream (RIS), West Antarctica, is an example of an ice stream where tides are known to significantly affect flow speeds. A puzzling aspect of the tidal response of RIS is the fact that the largest tidal modulation takes place over long tidal periods. (Long period tides are defined as having periods longer than those of any diurnal tides, i.e. periods significantly larger than one day.) The presence of a long-period tidally modulated ice-stream flow is intriguing because long-periodic ocean tidal amplitudes are small in comparison with the main semi-diurnal (S2 and M2) and the main diurnal components (K1 and O1). For example, a tidal analysis of a 55-day GPS record obtained about 20 km downstream of the RIS grounding line shows vertical amplitudes of MSf and Mf to be statistically insignificant, and no larger than few cm at the most, and the amplitudes of the S2, M2, K1, and O1 all to be on the order of a meter. Nevertheless, the response in ice-stream flow is stronger at the MSf period of 14.76 days than at any of the semi-diurnal and diurnal periods (Gudmundsson, 2006). A linear system, when forced over a given range of frequencies, will only produce a response at those same frequencies. RIS responds strongly at frequencies absent in the forcing, a clear evidence for some sort of a non-linear system response.
Vertical motion of a floating ice shelf is generally a faithful representation of the response of the ocean free-surface height to tide. The only exception to this is the region in the vicinity of the grounding line where stresses within the ice can support some of the weight of the ice column. Here the main focus is not on this vertical aspect of the ice motion, but on the horizontal flow response to tide of the grounded ice upstream from the grounding line.
This study is an extension of a previous modelling effort (Gudmundsson, 2007) suggesting that the observed longperiod tidally induced variations in flow on RIS are indicative of non-linear basal processes. The modelling work in Gudmundsson (2007) was done using a simple conceptual model of the interaction between ocean tides and ice stream flow. On the basis of that modelling work, it was concluded that a non-linear basal boundary condition of the type commonly used in glaciological modelling work, has the potential to produce the type of non-linear response observed on RIS. This study extends and complements earlier modelling efforts by including a number of processes not included previously. The ice is modelled as a non-linear visco-elastic medium and the effects of all the components of the equilibrium equations are included in the numerical model. In contrast, in Gudmundsson (2007), the contribution of ice deformation to ice-stream flow was ignored, and the basalstress perturbation was not calculated directly but rather parametrised in terms of the ocean tidal amplitude. Here the basal-stress perturbation is calculated from first principles, i.e. by solving the field equations describing the conservation of mass and momentum for given rheological models of ice and subglacial till.

Data
Although the main focus of this study is on investigating the general role of basal control on ice-stream response to oceanic tides, and not on reproducing the exact response curves from any one particular ice stream, data collected  -20), at the grounding line (R+00), and at distances of 20 and 40 km upstream from the grounding line (R+20 and R+40, respectively). The displacement curves shown are based on a tidal fit to the original measured data. on RIS is most pertinent to this modelling study. Therefore, some of the observational data from RIS are shown in Figs. 1 and 2. A more detailed discussion of these data, and other similar data sets from the same ice stream, are found in Gudmundsson (2006Gudmundsson ( , 2007; Murray et al. (2007); Adalgeirsdóttir et al. (2008); Dach et al. (2009);King et al. (2010). Figure 1 shows linearly detrended displacements curves from RIS. The displacements are along the mean flow direction at each site. In Fig. 2 calculated long-periodic (longer than than one day) tidal modulations in flow speeds are depicted for the same sites as in Figure 1. As seen the in Fig. 1 there is a prominent long-periodic modulation found in the displacement curves from RIS.
As is evident from a simple inspection of the data shown in Fig. 1, and as quantified in a more detailed tidal analysis (Gudmundsson, 2006), the long-period amplitudes at all sites are larger than any of the semi-diurnal and the diurnal tidal amplitudes. As an example, a tidal analysis of a further 71 day long GPS record collected 73 km upstream from the grounding line in the period from December 2005 to mid February 2006 shows the MSf tidal amplitude to be several times larger than any of the other tidal constituents. In that record, the MSf amplitude was estimated to be 4.7 times larger than that of the similarly period Mf tide (The period of the MSf tide is 14.76 days, and that of the Mf tide 13.66 days). Other GPS record from RIS give similar results King et al., 2010). At all sites shown in Fig. 1 The long-period tidal signal arrives at slightly different times at different locations (see Fig. 2). In Gudmundsson (2006) an order-of-magnitude estimate of the phase velocity of 1 m s −1 is given. A more detailed analysis of the same data set indicates that although this estimate of phase velocity is correct within an order of magnitude, the actual propagation speed could be anywhere between 0.2 to 1 m/s. The considerable errors in this estimate are mainly due to the limited length of available temporal coinciding records of only about 50 days. Estimating the phase shifts is confounded by the fact that the signal is dispersive and phase shifts not independent of the tidal period. Consequently, an estimate of propagation speed based on cross correlation of displacement curves from different sites does not give the same result as a calculation derived from relative phase shift estimates for individual tidal components at different sites. The estimate of 0.2 to 1 m s −1 is based on an analysis of the phase relationship of the MSf tide at four different sites along the medial line spaced 10 km apart. Adalgeirsdóttir et al. (2008) give an estimate of 10 ± 4 m s −1 for the phase velocity based on data collected on RIS at sites approximately 3 km apart. The difference between the estimate by Adalgeirsdóttir et al. (2008) of 10 ± 4 m s −1 and the 0.2 to 1 m s −1 given above, appears too large to be due to methodological differences only, and the source for the discrepancy between these estimates is unclear.
Apart from these tidally induced variations in flow, there appears to be no significant temporal changes in the flow of RIS (Gudmundsson and Jenkins, 2009). There also appears to be no clear temporal pattern in basal seismicity related to either tidal forcing or to the long-period variation in flow (H. Pritchard, personal communication, 2011). However, there are conflicting reports on the relationship between tides and seismicity on RIS. Adalgeirsdóttir et al. (2008) Fig. 3. Schematic showing an ice stream flowing into the ocean and forming an ice shelf. The pressure from the ocean acting on the ice is depicted as red arrows. Tides cause temporal changes in the oceanic forces and can lead to changes in ice-stream flow. Note that the figure is only a schematic and that in model runs the actual detailed model geometry is not as shown. Information on model geometry for each numerical experiment is given in the text.
concluded that "there is no simple relationship between the ocean tide and the velocity and basal seismicity of the ice stream". On the other hand Murray et al. (2007) writes that "The ice stream has two-weekly cycles in downstream flow and basal seismicity". Somewhat confusingly, Murray et al. (2007) gives the source of the two-weekly cycle in seismicity as Adalgeirsdóttir et al. (2008).

Model
The model setup is shown schematically in Fig. 3. The numerical model is a two-dimensional full Stokes flow-line model of ice-stream/ice-shelf flow. The numerical calculations were performed with the commercial finite element analysis software MSC.Marc (MARC, 2010).
The field equations are Dρ Dt + ρv q,q = 0 (1) representing the conservation of mass, linear momentum, and angular momentum, respectively, for a slowly moving medium. In the above listed equations, D/Dt denotes the material time derivative, v i are the components of the velocity vector, σ ij are the components of the Cauchy stress tensor, and f i are the components of the gravity force per volume. All terms of the equations listed above are included in the numerical model.

Ice rheology
Over tidal time scales ice behaves as a visco-elastic medium (Jellinek and Brill, 1956;Morland and Spring, 1981). Linear elastic models have been used to describe ice rheology over tidal periods, but the limitations of this approach have been pointed out by Reeh (2003). In his study of tidal flexure, Reeh (2003) used a linear visco-elastic four-element Burgers model to describe the rheology of ice. As a part of this study both a non-linear four-element Burgers model, an extension of the Reeh (2003) model to non-linear viscous rheology, and a non-linear two-element Maxwell model were used. As explained in more detail in Appendix A, it was found that the parameters of the two-element Maxwell model could be selected in such a way as to closely mimic the rheological behaviour of the four-element Burgers model over all times scales of interest in this study. Using parameter values suggested by Reeh (2003), the only significant differences between these two models are for loading periods shorter than about 100 s. As there is no appreciable tidal loading at such short periods, all the modelling work presented here is based on the simpler two-element Maxwell model. In comparison to the use of the Burgers model, the Maxwell model allowed the selection of longer time steps and resulted in shorter computational times.
For the non-linear Maxwell rheological model the deviatoric stresses, τ ij , and the deviatoric strains, e ij , are related througḣ where G is the bulk modulus of the Maxwell model, A is a rate factor, and n the stress exponent (Christensen, 1982). The deviatoric strains and deviatoric stresses are defined as and respectively, where ij are the strains and σ ij the stresses. The effective stress τ is defined as i.e. as the square root of the negative of the second invariant of the deviatoric stress tensor. The superscript ∇ denotes the upper convected time derivative, i.e.
is the upper-convected time derivative of the deviatoric stress tensor τ , where v is the velocity. Equation (4) can also be written on the form and where λ is the effective relaxation time, and η the effective viscosity. These are referred to as "effective" quantities because for n = 1 both λ and η are not material proprieties but depend on the state of stress. For loading periods long in comparison to the relaxation time λ, the constitutive relations is the usual Glen-Steinemann constitutive law (Steinemann, 1954(Steinemann, , 1958Glen, 1955) commonly used in large-scale modelling of ice masses (see Eq. 4).
As is common in the treatment of viscoelastic materials (Shames and Cozzarelli, 1997) ice is considered elastic under hydrostatic pressure, i.e.
where K is the shear modulus of the Maxwell model. In a number of glaciological studies of tidally induced deformation, ice rheology has been approximated using linear elastic constitutive equations where τ ij = 2Ge ij and σ ii = 3K ii (e.g. Holdsworth, 1969;Lingle et al., 1981;Stephenson, 1984;Vaughan, 1995;Sykes et al., 2009). For a viscoelastic material such as the upper convected Maxwell model given by Eqs. (4) and (11), one can define an effective shear modulus G e through τ ij = 2G e e ij , and an effective bulk modulus K e , through σ ii = 3K e ii . However, these effective parameters will, in general, be dependent on time. Under oscillating loading, for example, the effective shear and bulk modulus of a viscoelastic material are functions of the loading period. Furthermore, for any viscoelastic material that responds purely elastically to hydrostatic pressure, the corresponding effective Poisson's ratio is also time dependent and approaches 0.5 for slowly varying loads.
Because ice is viscoelastic over tidal periods, studies using linear elastic models describing tidal deformation of ice use effective parameters that are not independent of the loading period. It is therefore somewhat difficult to use values derived from such studies to constrain a visco-elastic model of rheology. For that reason the rheological values used here are primarily based on Reeh (2003), which appears to be the most in-depth modelling study of visco-elastic behaviour of large ice masses done to date.
As shown in Appendix A, a Maxwell model with a Young modulus E = 4.8GPa and a Poisson's ratio µ = 0.41 gives the same response to tidal loading periods as the Burgers model used by Reeh (2003). In this study, values for Young modulus ranging from 1 to 5 GPa, and Poisson's ratios between 0.4 to 0.5 were used. None of the results presented depend critically on the particular numerical values used for these rheological parameters. Note that the Young's modulus of ice at loading frequencies of hours and days is considerably smaller than the dynamical Young's modulus of about 10 GPa derived from the propagation of sound waves (Schulson and Duval, 2009

Basal boundary condition
Upstream from the grounding line, and along the ice-bed interface, a power-law type sliding law of the form is used, where t b is the basal traction withn being a unit normal vector to the bed pointing into the ice, and v b is the basal sliding velocity The sliding law (12) is commonly used in glaciology (Cuffey and Paterson, 2010) and often referred to as Weertman sliding law. It has two adjustable parameters c and m. The parameter c is referred to as the basal slipperiness. The basal slipperiness can, in general, be expected to be function of various other quantities such was basal water pressure and small-scale basal topography, etc., and therefore to be a function of location. In most flow modelling work to date the basal slipperiness is tuned, sometimes using formal inverse methods, to mach measurements of velocity and geometry. The other free parameter of the sliding law is the stress exponent m. Despite the value of the stress exponent m demonstrably having a decisive effect on the results of transient modelling work on large ice masses (e.g. Joughin et al., 2010), and despite decades of intense efforts at putting some constraints on its possible range, no consensus has emerged on either realistic values for m or on the general applicability of Weertman sliding law in the context of large-scale ice-flow modelling work. Values ranging from 1 to infinity are commonly used in flow modelling (Cuffey and Paterson, 2010). Basal motion was simulated in the model by introducing a deformable layer of till. The till was modelled as viscous medium using a flow law of the same form as Glen-Steinemann constitutive law. This approach of introducing basal motion has been used in numerous numerical studies, and a recent example for this approach with detailed description can be found in Leysinger Vieli and Gudmundsson (2010).

Boundary condition along the ice-ocean interface
Downstream of the grounding line the ice is subjected to oceanic pressure (p w ) acting normal to the surface given by where ρ w is the specific density of ocean water, g is the gravity acceleration. The variable S stands for the vertical position of the ocean surface, which, because of tidal action, is a function of time. A coordinate system with the z axis pointing vertically upwards is used. The mechanical boundary condition along the ice-ocean interface is where σ is the stress tensor within the ice. The condition (16) was implemented as a linear elastic spring, where the pressure (p) acting normal to the ice is given by where k is the spring constant, z 0 the spring offset, and z the vertical position of the ice-ocean interface. Setting k = −ρ w g and z 0 = −S(t) gives Eq. (15). Using this approach, the ocean pressure acting on the ice is not specified directly as a boundary condition. Only the dependency of the ocean pressure on the geometry, as given by (17), is specified. Both the pressure p and the vertical position of the ice-ocean interface are solved for as a part of the solution procedure. The ocean pressure always acts normal to the ice-ocean interface.
Perturbations in stresses at the grounding line due to ocean tides are broadly caused by two different mechanism, (a) bending stresses (Holdsworth, 1969), and (b) an overall change in horizontal stress as the height of the ocean water column changes (e.g. Thomas, 2007). The first mechanism, i.e. flexure, only acts if a glacier has a floating tongue. Vertical deformation around the grounding line is primarily due to bending stresses, and beam theory, which ignores overall changes in horizontal stress, has successfully been used to analyse measurements of tidal flexure (e.g. Reeh, 2003). In the model used here the ocean pressure is at each location along the ice-ocean interface given as a function of water depth, and these two different processes are not separated in the treatment of the boundary, but both are included.

Finite element discretisation
In all calculations an eight-node isoparametric quadrilateral plane-strain element was used. nodal positions moved with the flow of the medium. Automated remeshing was used to limit element distortion. However, in most model runs calculations ended before the need of remeshing.

Results
The numerical model was used to calculate the tidal response of an idealised ice stream to ocean tides. Model calculations were performed using ice-stream geometries based on that of RIS and for tidal amplitudes typical for that region of Ronne Ice Shelf. In line with the generic character of the modelling exercise, the geometry of RIS along the medial line was not replicated in exact detail. However, average thickness and slope were based on the RIS geometry (details given below). Fig. 4 shows modelled ice stream response at a distance 11 km upstream from the grounding line. The figure illustrates the effect that changing the value of m from 1 to 3 has on modelled tidal response (red and blue curves shown in Fig. 4). The only differences between these two runs are the values of the stress exponent m and the mean basal slipperiness c. The tidal forcing, model geometry, and ice rheology are in both cases identical. The basal slipperiness was changed as m was changed to ensure that the surface velocity was similar in both cases, or about 1 m/d.
The domain of the finite-element model used in producing the data shown in Fig. 4  ical behaviour of ice was described by a non-linear Maxwell model (see Eq. 11). The model has four adjustable parameters: The Young modulus E, the Poisson's ratio ν, a rate factor A and a stress exponent n. The rate factor of the ice was set at 4.5×10 −12 a −1 kPa −3 which corresponds to a temperature of about −20 degrees Celsius. The Young modulus E and the Poisson's ratio ν were 1GPa and 0.45, respectively. The thickness of the till layer was set at 250 m and the rate factor of the till was tuned to give a surface velocity of about 1md −1 upstream from the grounding line. Ice and ocean densities were ρ = 917kgm −3 and ρ w = 1030kgm −3 , respectively. As Fig. 4 shows, no long-period tidal modulation in ice stream flow is generated for m = 1. On the other hand, for m = 3 not only is there a long-period variation in flow, the amplitude of the long-period tidal modulation upstream from the grounding line is several times larger than the modulation at the semi-diurnal and diurnal periods. The amplitude of the fortnightly period in Fig. 4 is about 30 cm, or of the same order of magnitude as measured variability in ice-flow at that period on RIS (see Fig. 1 Figure 5 shows modelled ice-stream response to ocean tidal forcing when forced with the two main semi-diurnal tidal components S2 and M2 only. Hence, in this run the model was not subjected to any long-periodic forcing. Response is shown for m = 3 at three different sites at distances of 11, 21, and 31 km upstream from the grounding line, respectively. The model parameters are slightly different from those used in 4 and are listed in the figure caption of Fig. 5.
As seen Fig. 5, despite no forcing at long tidal periods, the strongest response, as measured by the amplitude of the detrended horizontal surface displacement, is at the MSf frequency. The MSf frequency is the difference between M2 and S2 frequencies, and the response at the MSf frequency is a nonlinear contribution of the forcing by the M2 and S2 tidal constituents.
Close inspection of the displacement curves in Fig. 5 reveals that they are phase shifted with respect to each other. The phase is a consequence of the visco-elastic rheology of ice. Calculating the cross correlation between the displacement curves to determine the phase speed gives a phase speed of 0.25 m s −1 which is comparable to the observed phase speeds on RIS of 0.2 to 1.0 m s −1 . The modelled phase speed is expected to depend on the parameter values of the rheological model and a detailed sensitivity study has not been performed.
Forcing the same model with the diurnal tidal components K1 and O1 only, results in a rather complicated looking response were the long period tidal components are mostly absent (see Fig. 6). The difference between the O1 and the K1 frequency is the Mf frequency. The results shown in Fig. 6 indicate that the strength of the non-linear interaction between K1 and O1 is not sufficient to produce a sizable Mf tidal modulation in flow. Note that despite the large differences between the results shown in Figs. 5 and Fig. 6, the model setup, i.e. rheological parameters and geometry, are in both cases identical. The difference in response is entirely due to the different oceanic forcing applied.
A feature of the non-linearity of the tidal response is to cause a mean shift in surface velocities. An example of this effect is given in Fig. 7 showing the horizontal component of the surface velocity vector as a function of time for both m = 1 and m = 3. For m = 1 the perturbation in velocity is symmetrical around the mean velocity. For m = 3 the perturbation is, on the other hand, asymmetrical. As a consequence, for m = 1 the mean velocity is shifted, and there is a net contribution to forward motion through the tidal action. In the particular modelling experiment shown in Fig. 7 tidal forcing causes about 5% increase in mean speed for m = 3. Figure 8 shows the dependency of the tidal response to the value of the stress exponent m. Shown in the figure are detrended horizontal displacement curves 30 km upstream from the grounding line as a function of time for different values of m. In each run, the basal slipperiness was adjusted to ensure that the surface velocities was at this site was 1 m/d irrespectively of the value of m. Clearly visible in the figure is how the amplitude of the fortnightly MSf horizontal tide increases with increasing m. As discussed above, for m = 1 there is no corresponding response at the MSf frequency. The amplitude of the semidiurnal tide, which in Fig. 8 can be seen superimposed on the longer period MSf tide, also increases markedly with m. vertically. This effect can be seen for both the long-period MSf tide, where the maximum in detrended displacement is reached about two days earlier for m = 10 than for m = 2, and the in the semidurnal variation (see Fig. 8).

Discussion
Modelled tidal modulation in ice stream flow is strongly sensitive to the parameters of the sliding law. For example, simply changing the value of the stress exponent from m = 1 to m = 3 gives rise to a long-period response in ice stream flow that is absent for m = 1 (see Fig. 4). Strong tidal modulation in ice-stream flow at long-tidal periods can be generated from the action of the semi-diurnal and diurnal oceanic tidal components alone, provided the relationship between basal stress and basal motion is non-linear (Fig. 5). The long-period response requires a non-linear mechanism. In the model the source of this non-linearity is the basal sliding law (for m = 1).
The qualitative difference in model response when forced with the S2/M2 tides (see Fig. 5), as compared to the response when forced with the K1/O1 tides (see Fig. 6), indicates that the model response to semi-diurnal loading periods is different from the response to diurnal loading periods. The viscoelastic rheology model introduces an additional timescale, i.e. the Maxwell time scale which is the ratio of the Young modulus and the (effective) viscosity (see Eq. 9). The presence of this time scale allows for a different types of model response to these two different types of tidal forcing. If no such additional timescale were involved, the model would not distinguish between forcing at semidiurnal and diurnal time scales, and the non-linear interaction between O1 and K1 would be similar to that between S2 and M2.
One of the consequences of the different strengths of the non-linear interaction between the S2/M2 pair and the O1/K1 pair, is that the long-periodic response is primarily concentrated at the MSf frequency. In comparison, the response at the Mf frequency is small. In qualitative terms this is in good agreement with observations from RIS where the horizontal MSf tidal amplitude is several times larger than the Mf tidal amplitude .
Although an exact comparison with data from RIS is not justified for a two-dimensional flow-line model, the modelled tidal response using a viscous sliding law with a exponent of m = 3 shares all the main characteristics of observed tidal modulation of RIS. A full parameter study has still to be performed. However, a value of m = 1 can be excluded.
Note that the modelled temporal flow variability is not due to any corresponding temporal changes in the model parameters. All modelling parameters are kept constant and do not change in space or time. In particular, although basal stresses and basal motion varies in space and time, the parameters of the sliding law do not. In nature one can expect the basal slipperiness to vary across the bed, and possibly also in relation to tides. In modelling terms, prescribing such a variation poses no difficulties. However, although such variability in basal slipperiness can be expected to modify the modelled response, and introducing such a variability might be useful as a part of a model-optimisation study, no such variability is required to reproduce the general characteristics of the tidal response observed on RIS.
The model produces no clear long-term tidal response when forced with diurnal tides only. This raises the possibility that the difference between observed tidal response on RIS and on some of the Siple Coast ice streams, such as the Bindschadler Ice Stream (e.g. Anandakrishnan et al., 2003), is primarily related to differences in forcing rather than differences in basal conditions. Ross Sea tides are predominantly diurnal with K1 and O1 amplitudes of about 0.5 m or less, and semi-diurnal amplitudes of less than 0.1 m (MacAyeal, 1984). The tides on Ronne-Filchner Ice Shelf are mixed diurnal and semi-diurnal tides with M2, S2, K1 and O1 amplitudes of around 1 m (Robertson et al., 1998;Fricker and Padman, 2002;King et al., 2011b). Hence, in comparison to RIS, the Siple Coast ice streams are subjected to at least ten times smaller semi-diurnal forcing. For these different types of ocean forcing, the modelled response would be different (compare Fig. 5 and Fig. 6) and yet in both cases similar to the observations, i.e. largest response is concentrated at long-periodic tides on RIS with little or no response at long-periods for Siple Coast ice streams.
The numerical model used here supports conclusions based on the previous modelling approach of Gudmundsson (2007), and gives added confidence in the applicability of that model to quantify effects of ocean tides on ice-stream flow as done by King et al. (2010). Modelled MSf amplitudes in ice stream flow are, for example, for both models almost identical, and both models produce a similar shift in mean surface velocity. However, there are a number of important differences between the visco-elastic full Stokes model presented here and the simple conceptual model of Gudmundsson (2007). Here, the basal perturbation in stress is calculated and is a model output, whereas in Gudmundsson (2007) the basal perturbation at each measurement site was an unknown model parameter that was estimated from the observed temporal variation in flow. Due to its simplicity, the model of Gudmundsson (2007) can be described as "an educated guess" of the effects of ocean tides on ice-stream flow. In the model presented here, the full set of the momentum equations are solved for a non-linear visco-elastic rheology without resorting to any simplifying modelling assumptions regarding the stress state.
For both the conceptual model presented in Gudmundsson (2007), and the model presented here, predicted ice-stream response to ocean tides does not scale linearly with the amplitude of the tides, unless for m = 1 and n = 1. These are basic characteristics of any non-linear model. Murray et al. (2007) discuss the fact that the velocity on RIS appears not to be "simply" related to tidal height, and state that the this observation invalidates the Gudmundsson (2007) model. Although the exact meaning of the word "simply" as used in this context by Murray et al. (2007) is not fully clear, both the argument and the conclusion are incorrect. In fact, a further analysis of the data presented in Murray et al. (2007) done by King et al. (2010) showed that the model of Gudmundsson (2007) could be used to reproduce that data set using "strikingly similar" parameter values to those of Gudmundsson (2007).

Conclusions and outlook
Applying ocean tides to a model of a coupled ice-stream/iceshelf flow can give rise to tide-induced horizontal movement of the ice-stream that is strongest at frequencies not presented in the forcing, provided a non-linear function is used to describe the relationship between basal motion and tangential basal stresses. Forcing the model with semi-diurnal tides only, can cause a strong fortnightly response in horizontal displacement upstream of the grounding line. This model response bears strong similarities to observations made on Rutford Ice Stream (RIS).
Using a Weertman type basal sliding law with moderately large stress exponent (i.e. within the range from 2 to 10) the numerical model is found to replicate all main qualitative features of the observed tidal motion of RIS, such as the genesis of long-period tidal modulation in flow in response to diurnal and semi-diurnal tides. The model also explains why the long-period response is concentrated at the MSf pe-riod of 14.76 days and almost absent at the similar Mf period of 13.66 days. Furthermore, when compared with data from RIS the model gives realistic order-of-magnitudes for amplitudes and phases of tidal constituents. Forcing the model with the O1/K1 diurnal tides only results in a fairly complicated tidal response that is mostly devoid of long-period components.
The phase relationship between the tidal response and the ocean tide is dependent on the parameters of the rheological model and on distance from the grounding line. Lowest speed is not observed at high tide and highest speed not at low tide, as might be expected if the change in horizontal pressure due to the varying height of the water column was the primary variable affecting the flow. The phase relationship is complicated, and in the model runs presented above highest forward speed occurs at different times depending on location. Sufficient distance downstream of the grounding line highest horizontal speed occurs approximately at highest rate of rising ocean tide.
There are a number of important issues not fully resolved here that warrant further modelling efforts. Of interest is the prospect of conducting a fully three-dimensional study of tidal modulation in flow. Such a study would constitute a much stronger test on the validity of the mechanism proposed here for the generation of tidal motion on ice streams, and deliver firmer constrains on the basal boundary condition.

Visco-elastic rheology models
The constitutive equations of linear viscoelastic materials under multiaxial stress can be written as and where τ ij and d ij are the deviatoric stresses and strains, respectively, and σ ii and ii (summation implied) are the volumetric stresses and strains (Shames and Cozzarelli, 1997). P d , Q d , P v , and Q v are differential time operators specific to a particular rheological model. The Maxwell model is a two-element model where a spring element and a viscous dashpot element are connected in series. A two-element model where a spring element and viscous dashpot element are connected in parallel is referred to as the Kelvin model. The Burgers model is a four-element model where a Maxwell and a Kelvin model are connected in series.
The Burgers model (Shames and Cozzarelli, 1997) is defined by and where Here η K and η M are the viscosities of the Kelvin and the Maxwell parts of the Burgers model, respectively, while G K and G M are the corresponding shear moduli. The parameter K is the bulk modulus. The volumetric deformation is assumed to be elastic. For oscillatory deviatoric strain d ij = d • ij e iwt , and deviatoric stresses τ ij = τ * ij e iwt , the complex shear modulus of the Burgers model, defined asĜ B = τ * ij /(2d • ij ), is given by The complex bulk modulus, K , of the model is independent of frequency, i.e. K = K.
As is typically the situation for visco-elastic bodies (Findley et al., 1976), the Poisson's ratio (defined as the negative of the ratio between lateral and axial strain under uniaxial stressing), for the Burgers model is a function of the loading frequency, The Poisson's ratio (ν) is, therefore, not a material parameter and, in general, time dependent. For w → +∞, ν → (3K − R)/(6K + R), and for w = 0 we find ν = 1/2 corresponding to an incompressible material. Similarly, the Young modulus (E) of a visco-elastic body is also not a material parameter, but depends on the loading period. The Burgers model is one of the simplest rheological models possible that can represent instantaneous elastic strain, delayed elastic response (primary creep), and steady-state viscous deformation (secondary creep). The instantaneous elastic response of the model is determined by the bulk modulus K and the shear modulus G M . The model parameters G K and η K determine the delayed elastic response, and η M the steady-state viscous deformation.
Values for the instantaneous shear and bulk moduli are listed by Röthlisberger (1972). The temperature dependency of that data was analysed by Hutter (1983). Following Reeh (2003), here the values G M = 3.5GPa and K = 8.9GPa are adapted. These values correspond to an instantaneous Young module E = 9.3GPa and an instantaneous Poisson's ratio ν = 0.33. Using results by Brill and Camp (1961) from studies of primary creep, and further following Reeh (2003), gives G K = 3.3GPa and η K = 600GPas. These values imply a retardation time of Kelvin element of the four-element fluid of only a few minutes, and suggest that the simpler Maxwell model may well be an equally good approximation to ice rheology over tidal time periods of hours and longer.
The complex shear modulus (Ĝ) of the Maxwell model is where η and G are the material parameters of the model (not to be confused with η M and G M which are the material parameters of the Maxwell part of the Burgers model).
We want Maxwell to reproduce Burgers over periods of interest. This is done selecting an effective G of the Maxwell model such that the complex shear modulus of the Maxwell model is a good approximation to the complex shear modulus of the Burgers model for tidal frequencies. One way of achieving this goal is by ignoring the delayed elastic response of the Burgers model and setting the parameters of the Maxwell model to and η = η M .
Equations (A12) and (A13) give the relationship between the values of the Burgers model and that of a Maxwell model that reproduces the Burgers model for loading periods large in comparison to the retardation time of the Burgers model. As shown in Fig. A1, using these values, the Maxwell model closely reproduces the behaviour of the more complex Burgers model over all frequencies larger than about 0.05 per day. The shear modulus of the (effective) Maxwell model is 1.7GPa and the instantaneous Poisson's ratio ν = 0.41 implying a instantaneous Young modulus of 4.8GPa.