Temporal variations in the flow of a large Antarctic ice stream controlled by tidally induced changes in the subglacial water system

Abstract. Observations show that the flow of Rutford Ice Stream (RIS) is strongly modulated by the ocean tides, with the strongest tidal response at the 14.77-day tidal period (Msf). This is striking because this period is absent in the tidal forcing. A number of mechanisms have been proposed to account for this effect, yet previous modelling studies have struggled to match the observed large amplitude and decay length scale. We use a nonlinear 3-D viscoelastic full-Stokes model of ice-stream flow to investigate this open issue. We find that the long period Msf modulation of ice-stream velocity observed in data cannot be reproduced quantitatively without including a coupling between basal sliding and tidally induced subglacial water pressure variations, transmitted through a highly conductive drainage system at low effective pressure. Furthermore, the basal sliding law requires a water pressure exponent that is strongly nonlinear with q = 10 and a nonlinear basal shear exponent of m = 3. Coupled model results show that sub-ice shelf tides result in a ∼12 % increase in mean horizontal velocity of the adjoining ice stream. Observations of tidally induced variations in flow of ice streams provide stronger constraints on basal sliding processes than provided by any other set of measurements.


Introduction
The majority of ice streams in Antarctica are forced at their boundary by ocean tides, either directly or through the motion of an adjoining ice shelf.Measurements have shown the flow of ice streams to be greatly affected by ocean tides over large distances upstream from the grounding line (Anandakrishnan et al., 2003;Bindschadler et al., 2003a, b;Gudmunds-son, 2006;Murray et al., 2007;Marsh et al., 2013).On Rutford Ice Stream (RIS), West Antarctica, for example, flow velocities change by more than 10 % in response to tides over distances of 50 km upstream from the grounding line.Several different types of tidally induced perturbations in ice flow have been observed on Antarctic ice streams.These include stick-slip motion observed at Williams Ice Stream (Bindschadler et al., 2003a, b;Winberry et al., 2009Winberry et al., , 2011)), smooth diurnal variations observed on Kamb and Bindschadler Ice Streams (Anandakrishnan and Alley, 1997;Anandakrishnan et al., 2003), and long-periodic response found on RIS and on several other ice streams flowing into the Ronne Ice Shelf (Gudmundsson, 2006;Murray et al., 2007;Aðalgeirsdóttir et al., 2008;King et al., 2010;Marsh et al., 2013).
An interesting aspect of the tidal observations on RIS is the long period (> 1 day) modulation in ice-stream flow that clearly demonstrates a nonlinear response to the tidal forcing (Fig. 1).In the response to the ice stream, the dominant tidal amplitude is found at the M sf tidal frequency (14.77 days), despite this tidal component being statistically insignificant in the tidal forcing.Hence, the strongest response is found at a frequency absent in the forcing.The same pattern is seen in observations of the tidal response of other ice streams flowing into Ronne Ice-Shelf (unpublished), as well as on the Larsen C Ice-Shelf (King et al., 2011).Note that flow modulation at M sf frequency is not simply a harmonic beat of the two semidiurnal frequencies; in fact it is a property of spectral analysis that tidal amplitudes can never arise through linear superposition of other frequencies.
One of the key motivations for studying the impact of tides on ice-stream flow is that modelling work has shown the response to reflect mechanical conditions at the glacier bed.Hence, observing and modelling tidally induced modulations in ice-stream motion provides a window into the mechanisms that influence basal sliding.
As initially suggested by Gudmundsson (2006), a nonlinear sliding law offers a potential explanation for the RIS observations, and various flow-line and full 3-D full-Stokes models have now successfully reproduced the general aspects of the long-period modulation in ice-stream flow as arising from a nonlinear response to tidal forcing (Gudmundsson, 2007;King et al., 2010;Gudmundsson, 2011;Walker et al., 2012;Rosier et al., 2014b).These previous studies, however, have primarily focused on identifying a potential mechanism giving rise to the observed nonlinear tidal response on RIS by reproducing the observations qualitatively.So far, with the notable exception of the recent work by Thompson et al. (2014), no modelling work has attempted to replicate the RIS observations in any quantitative detail.The models presented so far have shown that the qualitative aspects of the long-period RIS response can arise through transmission of tidally induced stresses across the grounding line, provided the sliding law is sufficiently nonlinear.In these models the physical conditions upstream of the grounding line, as defined in these models through their sliding-law parameters, do not change with time in response to tides.
The motivation for this work are recent modelling studies that suggest that any models using time-invariant sliding-law parameters, while ignoring the effects of tidally induced subglacial pressure variations on sliding, will fail to reproduce the RIS observations in quantitative terms.Recent work by Thompson et al. (2014), which does not explicitly investigate long-period modulation but includes the effects of ice-stream margins, found that for realistic ice-stream geometries, the effect of tidal stress perturbation on flow is too small to account for observations.In addition to this, our own 3-D modelling study including side drag and capable of reproducing the long period modulation, produced M sf amplitudes much smaller than those observed (Rosier et al., 2014b).As a result of the discrepancies outlined above, the question as to what mechanism can lead to the observed fluctuations in surface ice velocity still remains an open one.
The first measurements of this effect made by Gudmundsson (2006), suggested M sf amplitudes of ∼ 0.3 m at the grounding line, decaying to ∼ 0.1 m at 40 km upstream and still present at 73 km upstream.The model described by Gudmundsson (2011), although correctly producing strongest tidal response at the M sf frequency, appears only to be capable of reproducing M sf amplitudes of ∼ 0.1 m at most.In a more recent fully 3-D study, that in contrast to Gudmundsson (2011) included lateral drag, this amplitude is decreased further to ∼ 0.05 m at the grounding line when forced with the same tidal regime as that of the RIS (Rosier et al., 2014b).Hence, the observed response at the M sf frequency in that model is an order of magnitude too small.Thompson et al. (2014) conclude that the observed effect is too strong to be produced by transmission of tidal stresses only and suggest that a tidally driven time-dependent variability in till strength through hydrological coupling could explain the observed M sf response.
Here we use a 3-D nonlinear visco-elastic model with a geometry closely matching that of RIS to investigate the causes for the observed tidal response.We couple our icemechanical model to a model describing the changes in basal water pressure due to ocean tides, by allowing basal velocity to change in response to changes in effective basal water pressure.
The paper is organised as follows.We first describe our nonlinear visco-elastic model and present the basic governing equations.We then perform a full-Stokes surface-to-bed inversion of medial line surface velocities to determine the time-averaged spatial distribution of basal slipperiness.We then establish in a thorough parameter study that the model of Rosier et al. (2014b) cannot reproduce the observed longperiod velocity fluctuations of sufficient amplitude to agree with observations.In particular, and in an agreement with Thompson et al. (2014), we find that the observations can not be replicated through the effects of mechanical transmission of stresses through the ice and the till alone, but that in addition the effects of subglacial water pressure variations on sliding must be included.Finally we simulate perturbations in effective basal pressures due to ocean tides, and allow those changes in subglacial pressure to impact sliding through a commonly-used parameterisation relating sliding velocity and effective basal water pressure.After a new model parameter optimisation, we are able to replicate the RIS observations in considerable detail, but only within a fairly strict range of parameters.

Ice flow model
Our numerical ice flow model solves the field equations for conservation of mass, linear momentum (equilibrium equations) and angular momentum: where D/Dt is the material time derivative, ρ is density, ν 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.We use the comma to donate partial derivatives and the summation convention, in line with notation commonly used in continuum mechanics.None of the terms in the equilibrium equations are omitted.
In glaciology such models are commonly referred to as full-Stokes models.
We use an upper-convected Maxwell rheological model that relates deviatoric stresses τ ij and deviatoric strains e ij with ėij = where A is the rate factor, the superscript denotes the upper-convected time derivative, n is the constant in Glen's flow law (a nonlinear relation with n = 3 is used throughout), G is the shear modulus ν is the Poisson's ratio and E is the Young's Modulus.The upper-convected Maxwell model allows for calculation of large strain under rotation which, although not essential for the strains present in our model, we have chosen to use for completeness.More details of this rheological model can be found in Gudmundsson (2011).The deviatoric stresses are defined as and the deviatoric strains as where σ ij and ij are the stresses and strains, respectively.This rheological model approximates the visco-elastic behaviour of ice at tidal timescales, and can be thought of as a spring and dashpot in series such that the resulting strain is the sum of the elastic and viscous components and the stresses are equal.These equations are solved using the commercial finiteelement software package MSC.Marc (MARC, 2013).The ice stream and the underlaying till are treated as two separate deformable bodies.In a previous study we have calculated the migration of the grounding line in response to ocean tides, and accounted for the resulting effect on ice flow upstream from the grounding line in a flow-line setting (Rosier et al., 2014b).Due to computational considerations we have here, however, not allowed the grounding line to migrate over tidal cycles.
Basal velocity is given by a commonly used empirical form that includes effects of hydrology (e.g.Budd and Keage,, 1979;Bindschadler,, 1983): where τ b is the tangential component of the basal traction, N is the effective pressure (kept constant for the initial parameter study and subsequently perturbed due to the tide), c is basal slipperiness and both m and q are exponents.The slipperiness, tangential basal traction and effective pressure are all spatially variable.The effective subglacial water pressure N at the ice-till interface is defined as the difference between the normal component of the basal traction (σ nn , with a positive stress acting upwards) and the subglacial water pressure (p w ), i.e.N = −σ nn − p w , where a positive value for N indicates grounded ice where the downwards pressure of ice exceeds water pressure (as is the case everywhere upstream of the grounding line in this model).

Subglacial hydrology model
Our approach to including subglacial hydrology within the finite element model framework described in Rosier et al. (2014b) is to reduce the problem to the simplest possible set of equations.Rather than attempt to model a complex system of connected channels and distributed flow, we treat the drainage system as a homogenous porous medium with a characteristic 'conductivity' that, once coupled to the ice-flow model, can be tuned so that the velocity response matches observations.This approach to modelling subglacial hydrology has been used successfully in previous coupled studies eg.de Fleurian et al. (2014).
As a starting point we must lay out how the tide perturbs the subglacial water pressure.We write the subglacial water pressure (p w ) at any location upstream from the grounding line as where S is the mean ocean surface elevation, and S(t) the ocean tide.We incorporate the effects of the tides on subglacial water pressure through the grounding-line boundary condition for the perturbation in the hydrological head h.We assume that at the grounding line the subglacial water system is in direct contact with the ocean, and the subglacial water pressure at that location is therefore equal to the ocean pressure, or at x = x gl , and hence The tidally induced perturbation in hydrological head is then modelled as a diffusion process, i.e.
where K is the hydraulic conductivity.In the context of Darcy groundwater flow, K can be expressed as where k is the permeability, µ the viscosity of water, and S s the specific storage capacity.In reality this parameter combination is poorly constrained and here treated as an unknown.Thus, our approach is to solve for tidal perturbations in hydraulic head (rather than water pressure) which is known at the grounding line and transmitted upstream through a simple diffusion process controlled by the conductivity K.When modelling the spatial and the temporal variations of the subglacial drainage system water, we only attempt to describe the perturbations in effective pressure due to tides.This avoids the complications of calculating the temporally averaged pressure field, which is unnecessary as the effects of the mean pressure on basal flow are already accounted for in the temporally averaged value of the basal slipperiness which we derive in our inversion (see below).This is coupled to our ice-stream model through the sliding law (Eq.8) which we expand to consider perturbations in N: where and N is mean effective pressure such that N = N + N.
Re-arranging this gives where c = c N−q and ξ = N/ N.This now puts slipperiness and mean effective pressure into a new c term which is a function of x but not a function of t.In this way the baseline effective pressure and slipperiness conditions that affect the mean velocity of the glacier are separated from the perturbed terms.The c term is what is inverted for, as described later, to match observed medial line flow.Re-arranging the equation in this way means that N only affects the relative size of the non-dimensionalised perturbation ξ and not the mean flow which is constrained by observations.The hydrological coupling leads to six constants: N, K, ν, E, m and q which are treated as unknowns.The rheological parameters E (Young's modulus) and ν (Poisson's ratio) are constrained to some extent from previous visco-elastic modelling efforts on tidally induced motion, with values of E expected to be ∼ 4.8 GPa and ν of ∼ 0.41 (Reeh et al., 2003;Gudmundsson, 2011).The sliding law exponents m and q are treated here as tunable parameters.Note that once c has been determined, through the inversion procedure outlined below, K and N only affect modeled flow through their combined effect on ξ .Sensitivity of the model to the choice of these parameters is presented later.

Model geometry
Our model geometry is based on the RIS, however, we have not attempted to reproduce its geometry exactly and our thickness distribution in along-flow direction corresponds to the mean ice thickness across the ice stream.The 3-D model domain (Fig. 2) has zero bed slope, a surface slope of 0.0036 and ice thickness at the grounding line of 2040 m.This simple geometry is derived from average bed and surface profiles along the RIS medial line from BEDMAP2 data (Fretwell et al., 2013).While using constant slopes is a simplification and in reality the bed undulates considerably over the 100 km length being considered, there is no obvious overall shallowing or deepening, and the surface slope is relatively uniform.The width and length of the model domain are 16 and 120 km, respectively.The model width does not vary alongflow and the value chosen is an approximate average width for the region of interest.The hydrological component of the model extends a further 100 km upstream.

Boundary conditions
A no-slip condition is applied along one of the lateral boundaries and a free-slip condition along the other.The latter represents the ice stream medial line, giving an overall width of 32 km for the symmetrical problem approximately matching that of the RIS .Along the upper in-flow boundary, a surface traction is prescribed based on the analytical solution for the flow of a uniformly-inclined slab of ice. .The subglacial drainage system extends a further 100 km upstream from the ice-stream boundary.Note that since the problem is symmetrical, the medial line is considered to be the plane z = 0 and the ice stream being modelled is therefore 32 km wide.The term clamp is used to denote a node that cannot move in one or several degrees of freedom as indicated by the direction of the arrow.
where ρ i and ρ w are ice and water density respectively (910 kg m −3 and 1030 kg m −3 ), H is the ice thickness, s is the ice surface, z is the depth, g is gravitational acceleration and p b is a back pressure term that we treat as unknown and include in the inversion.In the case of the RIS, a non-zero value of p b could, for example, be expected to result from lateral resistance to ice-shelf flow.Two boundary conditions are necessary to solve for the diffusion of hydraulic head upstream from the grounding line.As mentioned earlier, at x = x gl subglacial water pressure and ocean pressure are assumed to be equal, leading to the boundary condition given in Eq. 11.At the upstream boundary the condition h → 0 as x → ∞ is strictly correct for this form of diffusion equation.Since this is not possible to implement in our model we use h = 0 at x = 200 km, assuming that h is very small at the upstream boundary.This can be justified analytically by solving Eq. 12 to give a decay length scale, for some periodic change in hydraulic head, of √ 2K/ω where ω is the tidal angular frequency being considered.For the range of conductivity values and tidal frequencies considered here, the model domain of 200 km is far larger than this length scale, thus this boundary condition can be safely applied without influencing the model results.
Ocean pressure is applied to the base of the floating ice shelf as a spring foundation (a more detailed description can be found in Rosier et al., 2014b) and the tidal forcing is introduced into the model as a perturbation in mean sea level.The tidal forcing is taken from the CATS2008 tidal model output (Padman et al., 2008), using the largest six tidal constituents at the RIS grounding line (M 2 , S 2 , O 1 , K 1 , K 2 and N 2 ).This model performs particularly well in this region since it is constrained by previous GPS measurements in this area and comparison with the vertical GPS record of Gudmundsson (2006) shows very close agreement.Tidal currents beneath the ice shelf are not included in the model since the effect on basal drag is negligible (Brunt, 2008;Makinson et al., 2012) and effects on basal melt are too slow to affect velocities at daily timescales.A schematic showing the various tidal processes, including some not included in the model, is shown in Fig. 3.

Model initialization
Preliminary experiments were conducted in which the stress exponent of the flow law (m) was changed to examine the effect on M sf response.Changing this parameter alters the mean flow in a non-trivial way that cannot be simply accounted for by altering slipperiness over the entire domain.Since the M sf response is sensitive to mean velocity it is important when comparing results to keep the mean velocity as close to observations as possible.To reproduce the general pattern of observed surface velocities on RIS, and in particular the general increase in velocities towards the grounding line, we invert for slipperiness (c ) using the medial line velocities obtained from the MEaSUREs InSAR velocity data set (Rignot, 2011) (note the term slipperiness here encompasses bed slipperiness and mean effective pressure).Although these InSAR-derived velocities are potentially flawed in regions with long period tidal modulation in flow (Gudmundsson, 2006), we address this by increasing the a priori error estimate (discussed later) to be larger than the errors provided in the data set.In general a comparison of the In-SAR velocities with in situ GPS measurements does show some differences but the only large discrepancy is on the ice shelf where we are not concerned with matching the velocities.
A Bayesian inversion approach was used to empirically calculate the i × j sensitivity matrix K describing the sensitivity of surface velocities to basal slipperiness.The method and equations are broadly similar to those presented in Raymond and Gudmundsson (2009) except that, rather than using analytical expressions for the sensitivity matrix, it is computed as the partial derivative of the forward model with respect to the state vector.The sensitivity matrix is given by where p and q are nodal numbers along the upper and lower surfaces of the finite element mesh.Here the measurement vector u has i elements and is the surface velocity, and the state vector c has j elements and is the slipperiness at the bed.Thus we calculate, for each element of the state vector, the change in measurement vector, giving one entire column of K.This is repeated for every element of the state vector to build up a complete sensitivity matrix.Since the model response to a change in slipperiness is nonlinear, the inversion will not converge to an optimum so- lution in a single iteration and so a Newton-Gauss iterative approach is used of the form where is the Fisher information matrix, S e is the covariance of measurement errors, S a is the covariance of a priori errors and F(c ) is the forward model (Rodgers, 2000).Measurement errors (σ e ) are assumed to be uncorrelated and have a normal distribution, such that the measurement error covariance matrix is proportional to the identity matrix, in the form S e = σ 2 e I m .We choose a large value of 0.2 m d −1 for σ e to account for errors arising from undersampling of tidal effects in this area.
Our treatment of the prior covariance matrix is the same as Gudmundsson and Raymond (2008), based on the assumption that basal slipperiness is spatially correlated, whereby each prior estimate of c at location i is related to a neighbouring location i − 1 by where c has variance σ 2 .The elements of S −1 a can then be given by where λ is a decay length scale, related to φ c by λ = −1/ ln φ c and the variance is This results in a covariance matrix which has σ 2 a along the diagonal and non-zero off-diagonal elements.
We reduce the number of calculations needed by only taking into account along-flow variations in slipperiness.This simplification is justified due to the simple geometry and because we only seek to match the medial line ice-stream velocity.Buttressing (p b , which is particularly relevant for flow velocities near the grounding line) is inverted for by adding a single non-dimensionalised element to the end of the state vector.This is treated in the same way as the other state vector elements apart from having its own (uncorrelated) prior error estimate.
Although this brute force approach to inverting for basal slipperiness is computationally more expensive than others such as the adjoint method, there are a number of advantages to this method such as giving an explicit estimate of the inversion error.Furthermore, because each element of the K matrix is independent of all the others, it is possible to easily parallelize its calculation, meaning that run times need not be orders of magnitude greater if sufficient computing resources are available.The sensitivity matrix need not be calculated for each iteration and in fact it is advantageous to iterate a number of times using the same matrix before re-calculating it.The iteration was continued until it converged on the maximum a posteriori solution, in contrast to many other similar studies which stop iterating once the misfit between model output and observations is below a given threshold.

Results
As discussed above, to date no model has been presented that can reproduce the tidally induced horizontal velocity variation observed on the RIS.Admittedly, most models have focused on trying to identify the mechanism responsible for the rather striking observation that the response of the ice stream is concentrated at tidal frequencies absent in the forcing.However, it would be expected that if the mechanism has been correctly identified, and is the primary cause for the velocity fluctuations, modeled amplitudes would be close to those measured.In fact modelling work presented so far has always produced too small a response at the M sf tidal period, and too strong at both diurnal and semidiurnal periods.

Modelling the tidal response of RIS assuming no temporal changes in water pressure
To address the open question of whether RIS observations can be replicated through stress transmission alone, our first modelling aim is to establish an upper bound on the possible M sf amplitude in the absence of any temporal changes in bed conditions, i.e. ōther than those resulting from direct stress transmission through the ice due to the flexing of the ice in response to tides.In the context of our modelling methodology described above this is equal to setting the stress exponent, q, of the effective water pressure in the sliding law (see Eq. 8) to zero.In effect we repeat the fully 3-D simulations conducted in Rosier et al. (2014b) but with a broader range of parameters, an ice-stream geometry closer to that of RIS and a basal slipperiness distribution (c(x)) determined through a formal inversion of surface velocities.Our tunable model parameters with no subglacial hydrological coupling are the Poisson's ratio (ν), the Young's modulus (E) and the stress exponent (m).We set the stress exponent (n) in Glen's flow law to n = 3, and determine the rate factor A from a static temperature distribution defined in the model using the commonly used temperature relation given by Cuffey and Patterson (2010).
We performed an extensive parameter study, with the stress exponent m of 1, 3, 5 and 10, and the Young's modulus of 3, 4.8 and 6 GPa.The Poisson's ratio was varied between 0.3 to 0.45, but was found to have almost no effect on the modelled tidal response and we do not discuss those results further.For every value of the basal sliding-law stress exponent m, we first determined the maximum a posteriori distribution of basal slipperiness (c ) using our inversion approach.In the surface-to-bed inversion the long-term average flow in the absence of tidal forcing was matched to the observed velocity, and a (purely) viscous flow model was therefore used in the forward step.We then forced our visco-elastic timedependent model by tides.For each given value of m and the associated basal slipperiness distribution, tidal response was calculated for a range of elastic rheology parameters.From modelled horizontal displacements curves, we then calculated tidal amplitudes and phases as a function of distance along the medial line.By fitting an exponential curve to the spatial variation in tidal amplitudes, we then determined decay length scales for each tidal component, as well as phase velocities.The decay length scale we refer to here is defined as the e-folding length scale, or the distance for a given signal (in this case the horizontal tidal signal) to decay by factor e.
The results of the parameter study are summarized in Fig. 4. In Fig. 4a the amplitude of the M sf frequency 10 km upstream from the grounding line is shown.The modelled M sf amplitudes are never larger than a few centimetres.The largest values are found for high m and high E values.Although somewhat higher M sf amplitudes could be obtained by increasing m even further, the modelled results show that this increase is sub-linear as a function of m.Furthermore, for m > 10 other model outputs that must match observations such as phase velocity, decay length scale and notably M 2 amplitude, would also increase beyond the range of desired values.The model is, thus, not able to reproduce the observed magnitude of the M sf tidal amplitude.
Both the decay length scale (Fig 4c) and phase velocity (Fig 4d) increase with increasing m, in agreement with the analytical solution derived in previous work (Rosier et al., 2014b).
The amount of buttressing needed to match observed velocities increases as m is increased and varied from 650 KPa to 850 KPa for m = 1 to 10.Note that the inversion procedure, in minimising the cost function, tries to find a solution that does not vary significantly from the a priori estimates of slipperiness and buttressing, and therefore this buttressing value may be to some extent artificial if the a priori buttressing estimate and error are poorly chosen.For this reason a large value (1000 kPa) is chosen for the error estimate of buttressing used in the inversion.Decay of the M sf amplitude upstream of the grounding line for all parameter study simulations is plotted in Fig. 5 (blue lines) and compared with the observed amplitudes (crosses).This clearly shows the disparity between desired amplitude and the range of possible amplitudes using the mechanism described above.The conclusion from this parameter study, in agreement with Thompson et al. (2014), is that stress transmission alone cannot explain the large amplitude of M sf modulation, with maximum amplitudes 10 km upstream approaching ∼ 0.05 m, considerably smaller than the desired 0.3 m.Clearly an additional nonlinear effect is needed to match observations.Although stress-transmission can reproduce the qualitative aspects of the data, in particular the generation of M sf response, the effects are (at the most) about an order of magnitude smaller than revealed by measurements.

Modelling the tidal response of RIS assuming temporal changes in water pressure
We now couple our hydrological model (Sect.2.2) to the 3-D full-Stokes model by using values of q > 0 in order to see whether this can explain measurements made on the RIS.Coupled model results obtained through optimization of hydrological parameters are shown in Fig. 6.This provides a much better agreement with GPS measurements than any previous combination of parameters for the model with no subglacial water pressure coupling.Notably, the M sf amplitude and decay length scale are both large and match very closely with data (Fig. 5).The hydrological model used a mean effective pressure ( N) of 105 kPa, pressure exponent (q) of 10 and conductivity (K) of 7 × 10 9 m 2 d −1 .Other The Cryosphere, 9, 1649-1661, 2015 model parameters were E = 4.8 GPa, ν = 0.41 (both in accordance with the optimum Maxwell rheology given by Gudmundsson, 2011) and m = 3.
The only feature of these results that is arguably not in agreement with observations is the amplitude of the semidiurnal tidal constituent detrended displacements.Comparison between Figs. 1 and 6 shows an M 2 amplitude (visible in both figures as the higher frequency modulation overlain on the long period M sf signal) that is approximately twice as large at the grounding line as the amplitude determined by tidal analysis of the data.Possible explanations for this are that the M 2 amplitude may be too small to be sufficiently resolved by the GPS receivers that originally made the measurements or limitations of the simple Maxwell rheology.Errors in the GPS measurements are of the order of centimetres; more details of the original data set can be found in Gudmundsson (2006) and a description of similar processing in Dach et al. (2009).
We perform a sensitivity analysis to determine whether the M sf response is robust or highly sensitive to certain parameters.Figure 7 shows change in M sf amplitude (panel a), M 2 amplitude (panel b), M sf decay length scale (panel c) and M sf phase velocity (panel d) compared to the optimized model for a ±10 % change in each parameter.
Comparison in Fig. 7a and b suggests that the calculated M sf and M 2 amplitudes are closely correlated and thus, for the parameters tested here, there is no clear modification of the model that would decrease the semidiurnal (M 2 ) amplitude without also reducing the M sf response.Softening the ice by reducing E may be one possible route, since this appears to increase M sf amplitude more than M 2 amplitude, however this parameter is more tightly constrained than others since the rheology of ice is not entirely unknown and the sensitivity is too small to solve the issue.M sf amplitude is most sensitive to normalized changes in N and q, as might be expected since it is the nonlinearity here that drives the majority of the long period modulation in flow.
A reduction in m increases the nonlinear response of the modeled ice stream, the reverse of the response with no hydrological coupling, but increases the M sf length scale and phase velocity.Overall, all parameters are most sensitive to the choice in N.This is not surprising, since N is small and as N gets large the dimensionless number ξ will drop out and that source of nonlinearity disappears.
The large difference in M sf amplitude between the parameter study simulations and those that include tidally induced subglacial pressure variation poses an important question; is a nonlinear sliding law where m > 1 required at all, given that the M sf modulation appears to be largely generated by water pressure changes.Results from the sensitivity analysis suggest that the stress exponent m remains a crucial parameter in altering characteristics of the M sf response.To look at this in more detail, the model was rerun with varying exponents q and m, with the aim of examining the characteristics of the M sf response given changes in the dominance of the two mechanisms.The four characteristics of the model's tidal response are plotted against exponents q and m, each varying between 1 and 10, in Fig. 8.These results show that reducing m leads to an increase in amplitude of both tidal frequencies investigated, but a decrease in the length scale and phase velocity.An M sf decay length scale of ∼ 50 km is observed on the RIS but panel c shows that for m = 1 the length scale is smaller up to q = 10 and in fact appears to have reached an asymptote.Increasing m for any given value of q however leads to a large increase in the length scale.The mechanism by which increasing m reduces M sf amplitude but increases length scale is discussed later but suggests that a flow low with m > 1 is still required to reproduce the RIS tidal response.

Discussion
We find that stress transmission alone cannot fully explain the observed M sf modulation of surface velocities on the RIS.An additional mechanism whereby a tidally induced pressure wave travels up a subglacial drainage system, altering the effective pressure at the base of the ice stream, is required to produce a sufficiently large M sf amplitude.The drainage system must be highly conductive and sufficiently nonlinear, such that a small change in basal water pressure leads to a large change in surface velocity.
This nonlinearity arises largely in two of the parameters: N and q.The model does not take into account feedback between ice flexure and water pressure.water pressure and it has been suggested that this mechanism could pump brackish water upstream (Walker et al., 2013;Sayag and Worster, 2013).This flexure may have the additional effect of opening crevasses beneath the ice or dilating the subglacial till, leading to changes in local water storage and thereby altering the distribution of water.Our justification in ignoring these additional processes is that ice flexure is limited to within several ice thicknesses of the grounding line and the M sf modulation is observed to travel much further upstream.Spatial variations in N are accounted for in the inverted c and cannot be separated from spatial variability in c.In reality if N varied spatially this would affect the nonlinearity in ξ .Ultimately we ignore this additional complication and the decay in ξ is only a function of the spatially uniform conductivity, K.In doing so several processes are combined to provide a more general picture of the subglacial drainage characteristics.A fit to observations could to some extent still be obtained if N was altered by compensating with a change in q since the two parameters are correlated.In general though, a relatively low value of effective pressure with no large gradient going upstream from the grounding line is needed, since a gradient would cause the nonlinearity to be rapidly reduced in the upstream direction.
In order to understand the interaction between the hydrology and stress transmission mechanisms it is important to consider the relative timing with which they act on the ice stream.As explained previously, an exponent m > 1 causes an increase in ice velocity during low tide and decrease at high tide.Conversely, at high tide near the grounding line the water pressure within the subglacial drainage system will be at its highest, lowering the effective pressure and increasing ice velocity.The two effects are therefore opposite in phase at the grounding line (although in both cases the peak velocities are still during the spring tide, so there is no phase shift in the M sf frequency that they generate at this point).Since the subglacial pressure effect is larger it dominates at the grounding line and the reduction in M sf amplitude at this point for m > 1 is a result of the stress transmission effect being 180 • out of phase with the subglacial pressure variations, thereby dampening the velocity modulation.
Results from Fig. 8 suggest that, while it may be possible to reproduce the observed M sf response of the RIS for m = 1, this would necessitate an almost infinite conductivity in order to transmit the signal far enough upstream.With the set of model parameters presented, the effect of subglacial pressure variations dominates at the grounding line and can produce very large M sf amplitudes, but what is much more difficult is to reproduce the long decay length scale of this frequency.
The key parameter then becomes m, which can substantially increase the decay length scale given values m > 1.Any reduction in the M sf amplitude from using a high value of m can be compensated for by increasing the nonlinearity of the drainage system (reducing N or increasing q).
None of the other parameters within the model had such a large effect on the length scale and the implication is that a nonlinear sliding law is required in addition to any nonlinear response to subglacial pressure variations.Matching the observed long period modulation of ice-stream flow requires a balance between large M sf amplitude and decay length scale.A choice of m that is too small means the M sf signal will decay too rapidly upstream of the grounding line, but too large and the generation of the signal due to subglacial hydrology becomes hindered.
An explanation for this increase in length scale with m > 1 can be thought of intuitively as follows.Consider the propagation of nonlinear M sf period up the RIS as two waves, generated by the upper and lower terms on the right of Eq. 8.These two waves clearly have the same frequency but since they propagate up the ice stream by different mechanisms it is reasonable to assume they have different phase velocities.At x gl they are 180 • out of phase but with different phase speeds this destructive interference becomes constructive interference as you move away from the source.As a consequence the M sf amplitude is reduced at the grounding line but its decay may be slowed as a result of constructive interference upstream.
The requirement of high conductivity in order to transmit the tidal signal far enough upstream to match observations suggests that there must be a channelized drainage system beneath the RIS.This could consist of a few large channels that transmit the tidal pressure wave far upstream which then permeates through the till on either side of the channel, leading to changes in effective pressure over large portions of the ice-stream base.Gudmundsson (2011) demonstrated that the nonlinearity described above leads to an increase in the RIS mean velocity of ∼ 5 % due to the presence of the tides.A simulation with identical model setup to that used in Fig. 6 but with tidal amplitude set to zero everywhere was done to examine this process with the larger M sf amplitudes presented in this work.The result with this new model, that successfully replicates the amplitude of long period modulation, is that mean surface velocity is increased by ∼ 12 % due to the presence of the tides.This is a considerable increase on the previous value which is expected since the M sf amplitudes in that model were smaller.It demonstrates that tidal forcing can not necessarily be ignored over longer timescales.Future changes in ice-shelf thickness and extent could lead to interesting feedbacks between tidal amplitudes and ice-stream velocities (Arbic et al., 2008;Rosier et al., 2014a).

Conclusions
Observations of surface motion of the RIS show a strong, nonlinear response that propagates a long way upstream from the grounding line.The nonlinear response of this ice stream and others in the region is striking both in its amplitude and extent and matching observations is not possible through stress transmission considerations alone.Coupling with a hydrological model that sends tidally induced subglacial pressure variations far upstream is required to explain these observations.Furthermore, three other requirements must be met; low effective pressure across the entire ice-stream bed, a highly conductive subglacial drainage system and a nonlinear sliding law such that m > 1.
Hydrological and basal sliding model parameters that produced a best fit to observations were m = 3, q = 10, K = 7 × 10 9 m 2 d −1 and N = 105 kPa.Although a complete exploration of the parameter space is not currently possible due to prohibitive computational expense, we are confident that the set of parameters outlined above is robust for our simplified 3-D model.Future models, incorporating detailed RIS topography, could further constrain these parameters.We know of no other approach that can provide these insights into the controls on basal motion.Our conclusion from attempting to match the observed nonlinear response of the RIS is that a channelized and highly efficient drainage system must exist at the bed in order to reproduce an M sf response of sufficient amplitude and extent.

Figure 1 .
Figure 1.Linearly detrended horizontal displacements on the RIS reproduced by a tidal fit to the original measured data.Measurements are shown from five GPS stations at 20 km downstream of the grounding line (R − 20 km), at the grounding line (R + 00 km) and distances of 20, 40 and 73 km upstream from the grounding line (R + 20 km, R + 40 km and R + 73 km, respectively).Data at 10 km upstream is not included for the sake of clarity.
where h(x, t) is the tidally induced perturbation in the hydrological head, ρ w is the ocean density, g the gravitational www.the-cryosphere.net/9etal.: Tidal forcing of Rutford ice stream acceleration, b(x) the bed elevation, and the ocean surface elevation S(t) is given by

Figure 2
Figure 2. 3-D model domain, showing the boundary forces (black arrows) and flow constraints (red arrows).The subglacial drainage system extends a further 100 km upstream from the ice-stream boundary.Note that since the problem is symmetrical, the medial line is considered to be the plane z = 0 and the ice stream being modelled is therefore 32 km wide.The term clamp is used to denote a node that cannot move in one or several degrees of freedom as indicated by the direction of the arrow.

Figure 3 .
Figure 3. Schematic showing the various mechanisms by which tides can influence ice-stream flow.Note that grounding line migration, crevassing and tidal currents are not included in the model.

Figure 4 .
Figure 4. Modelled M sf and M 2 tidal amplitudes 10 km upstream from the grounding line (a and b, respectively), and M sf decay length scales and phase velocities (c and d, respectively) as a function of the basal sliding law stress exponent m and the elastic Young's modulus (E) of ice.Here the potential effects of subglacial water pressure variations in response to tides on sliding were not included, i.e. in the sliding law (Eq.8), q = 0. Crosses indicate model simulations.The contour plot is based on interpolation of model results.M sf and M 2 amplitudes were taken at 10 km upstream from the grounding line.

Figure 5 .
Figure5.M sf amplitude as a function of distance upstream for parameter study simulations with no hydrology (yellow), compared with GPS measurements (crosses), all results from the coupled model sensitivity study (blue) and the coupled model best fit to GPS data (red).

Figure 7 .
Figure 7. Sensitivity analysis of model parameters (N , K, q, E and m), showing change in M sf and M 2 amplitudes (a and b), M sf decay length scale (c) and M sf phase velocity (d) for +10 % (white bar) and −10 % (grey bar) changes in each parameter.Model outputs were compared to the simulation presented in Fig. 6 and all other parameters were kept at the values defined in that plot.

Figure 8 .
Figure 8. Response of M sf amplitude (a), M 2 amplitude (b), M sf decay length scale (c) and M sf phase velocity (d) to choice of stress exponent (m) and hydrological exponent (q).Crosses indicate model simulations.M sf and M 2 amplitudes were taken at 10 km upstream from the grounding line.