Sensitivity of ice sheet surface velocity and elevation to variations in basal friction and topography in the Full Stokes and Shallow Shelf Approximation frameworks

basal friction and topography in the Full Stokes and Shallow Shelf Approximation frameworks Gong Cheng1, Nina Kirchner2,3, and Per Lötstedt1 1Department of Information Technology, Uppsala University, P. O. Box 337, SE-751 05 Uppsala, Sweden 2Department of Physical Geography, Stockholm University, SE-106 91 Stockholm, Sweden 3Bolin Centre for Climate Research, Stockholm University, SE-106 91 Stockholm, Sweden Correspondence: Gong Cheng (cheng.gong@it.uu.se)

In computational models of ice dynamics, the description of sliding processes, including their parametrization, plays a central role, and can be treated in two fundamentally different ways, viz. using a so-called forward approach on the one hand, or an inverse approach on the other hand. In a forward approach, an equation referred to as a sliding law is derived from a conceptual 5 friction model, and provides a boundary condition to the equations describing the dynamics of ice flow (in glaciology often referred to as the Full Stokes (FS) model) and which, once solved, render e.g. ice velocities as part of the solution. Studies of frictional models and resulting sliding laws for glacier and ice sheet flow emerged in the 1950s, see e.g. Fowler (2011); Iken (1981); Lliboutry (1968); Nye (1969); Schoof (2005); Weertman (1957), and have subsequently been implemented into numerical models of ice sheet and glacier behavior e.g. in Brondex et al. (2017Brondex et al. ( , 2019; Gladstone et al. (2017); Tsai et al. (2015); 10 Wilkens et al. (2015); Yu et al. (2018), and continue to be discussed (Zoet and Iverson, 2020), occasionally controversially (Minchew et al., 2019;Stearn and van der Veen, 2018).
Because little or no observational data is available to constrain the parameters in such sliding laws (Minchew et al., 2016;Sergienko and Hindmarsh, 2013), actual values of the former, and their variation over time (Jay-Allemand et al., 2011;Schoof, 2010;Sole et al., 2011;Vallot et al., 2017), often remain elusive. Yet, they can be obtained computationally by solving an 15 inverse problem provided that observations of e.g. ice velocities at the ice surface, and elevation of the ice surface, are available (Gillet-Chaulet et al., 2016;Isaac et al., 2015). Note that the same approach, here described for the case of the sliding law, can be used to determine other "inaccessibles", such as optimal initial conditions for ice sheet modeling (Perego et al., 2014), the sensitivity of melt rates beneath ice shelves in response to ocean circulation (Heimbach and Losch, 2012), or to estimate basal topography beneath an ice sheet (van Pelt et al., 2013). The latter is not only difficult to separate from the determination of the 20 sliding properties (Kyrke-Smith et al., 2018;Thorsteinsson et al., 2003), but also has limitations related to the spatial resolution of surface data and/or measurement errors, see Gudmundsson (2003Gudmundsson ( , 2008; Gudmundsson and Raymond (2008).
Adopting an inverse approach, the strategy is to minimize an objective function describing the deviation of observed target quantities (such as the ice velocity) from their counterparts as predicted following a forward approach when a selected parameter in the forward model (such as the friction parameter in the sliding law) is varied. The gradient of the objective function 25 is computed by solving the so-called adjoint equations to the forward equations, where the latter often are slightly simplified, such as e.g. by assuming a constant ice thickness or a constant viscosity (MacAyeal, 1993;Petra et al., 2012). However, when inferring friction parameter(s) in a sliding law using an inverse approach, recent work (Goldberg et al., 2015;Jay-Allemand et al., 2011) has shown that it is not sufficient to consider the time-independent (steady state) adjoint to the momentum balance in the FS model. Rather, it is necessary to include the time-dependent advection equation for the ice surface elevation in the 30 inversion. Likewise, but perhaps more intuitively understandable, the choice of the underlying glaciological model (FS model, vs. e.g. Shallow shelf/stream approximation (SSA) model, see Sect. 2), has an impact on the values of the friction parameters obtained from the solution of the corresponding inverse problem (Gudmundsson, 2008;Schannwell et al., 2019).
Here, we present an analysis of the sensitivity of the velocity field and the elevation of the surface of a dynamic ice sheet (modelled by both FS and SSA, respectively, briefly described in Sect. 2) to perturbations in the sliding parameters contained y coordinates, and z is the coordinate in the vertical direction, see Fig. 1(a). Specifically, we denote by u 1 , u 2 , and u 3 the velocity components of u = (u 1 , u 2 , u 3 ) T in the x, y and z direction, where x = (x, y, z) T is the position vector and T denotes the transpose. Further, the elevation of the upper ice surface is denoted by h(x, y, t), the elevation of the bedrock and the location of the base of the ice are b(x, y) and z b (x, y, t), and the ice thickness is H = h − z b . Upstream of the grounding line, The boundary Γ enclosing the domain Ω occupied by the ice has different parts, see Fig. 1(b): The lower boundaries in the ω plane are denoted by Γ b (where the ice is grounded at bedrock), and Γ w (where the ice has lifted from the bedrock and is floating on the ocean). These two regions are separated by the grounding line γ GL , defined by (x GL (y), y) based on the assumption that ice flow is mainly along the x-axis. The upper boundary in the ω plane is denoted by Γ s (ice surface) at 10 h(x, y, t) in Fig. 1(a). The boundary of ω itself is denoted by γ.
The vertical lateral boundary (in the x − z plane, Fig. 1(b)) has an upstream part denoted by Γ u in black and a downstream part denoted by Γ d in blue, where Γ = Γ u ∪ Γ d . Obviously, if x ∈ Γ u , then (x, y) ∈ γ u or if x ∈ Γ d , then (x, y) ∈ γ d where γ = γ u ∪ γ d . Letting n be the outward pointing normal on Γ (or γ in two dimensions (x, y)), the nature of ice flow renders the conditions n·u ≤ 0 at Γ u and n·u > 0 at Γ d . In a two dimensional vertical ice ( Fig. 1(a)), this corresponds to x = (x, z) T , ω = 15 [0, L], γ u = 0, and γ d = L where L is the horizontal length of the domain. In summary, the domains are defined as: Before the forward FS equations for the evolution of the ice surface Γ s and the ice velocity in Ω can be given, further notation needs to be introduced: ice density is denoted by ρ, accumulation and/or ablation rate on Γ s by a, and gravitational acceleration by g. On Γ s , h = (h x , h y , −1) T describes the spatial gradient of the ice surface (in two vertical dimensions h = (h x , −1) T ).
The strain rate D and the viscosity η are given by where trD 2 is the trace of D 2 . The rate factor A in (2) depends on the temperature and Glen's flow law determines n > 0, here taken to be n = 3. The stress tensor is where p is the isotropic pressure and I is the identity matrix.
Turning to the ice base, the basal stress on Γ b is related to the basal velocity using an empirical friction law. The friction 10 coefficient has a general form β(u, is independent of the velocity u and f (u) represents some linear or nonlinear function of u. For instance, f (u) = u m−1 with the norm u = (u · u) 1 2 introduces a Weertman type friction law (Weertman, 1957) on ω with a Weertman friction coefficient C(x, t) > 0 and an exponent parameter m > 0.
Common choices of m are 1 3 and 1. Finally, a projection (Petra et al., 2012) on the tangential plane of Γ b is denoted by T = I − n ⊗ n where the Kronecker outer product between two vectors a and c or two matrices A and C is defined by With these prerequisites at hand, the forward FS equations and the advection equation for the ice sheet's elevation and velocity for incompressible ice flow are where p w is the water pressure at Γ w , h 0 (x) is the initial surface elevation, and h γ (x, t) is a given height on the inflow 20 boundary. The boundary conditions for the velocity on Γ u and Γ d are of Dirichlet type such that where u u and u d are known. These are general settings of the inflow and outflow boundaries which keep the formulation of the adjoint equations as simple as in Petra et al. (2014). Should Γ u be at the ice divide, the horizontal velocity is set to u| Γu = 0.
The vertical component of σn vanishes on Γ u .

Shallow shelf approximation.
The three dimensional FS problem (4) in Ω can be simplified to a two dimensional, horizontal problem with x = (x, y) ∈ ω, by adopting the SSA, in which only u = (u 1 , u 2 ) T is considered. This is because the basal shear stress is negligibly small at the base of the floating part of the ice mass, viz. the ice shelf, rendering the horizontal velocity components almost constant in the z direction (Greve and Blatter, 2009;MacAyeal, 1989;Schoof, 2007). The SSA is often also used in regions of fast-flow 5 over lubricated bedrock (MacAyeal, 1989;Pattyn et al., 2012).
The simplifications associated with adopting the SSA imply that the viscosity (see (2) for the FS model) is now given by where B(u) = D(u)+∇·u I with ∇·u = tr D(u). This η differs from (2) because B = D due to the cryostatic approximation of p in the SSA. In (6), the Frobenius inner product between two matrices A and C is used, defined by A : The vertically integrated stress tensor ς(u) (cf. (3) for the FS model) is given by The friction law in the SSA model is defined as in the FS case. Note that basal velocity is replaced by the horizontal velocity. This is possible because vertical variations in the horizontal velocity are neglected in SSA. Then, Weertman's law is Above, t is the tangential vector on γ = γ u ∪ γ d such that n · t = 0. The inflow and outflow normal velocities u u ≤ 0 and u d > 0 are specified on γ u and γ d . The lateral side of the ice γ is split into γ g and γ w with γ = γ g ∪ γ w . There is friction in the 20 tangential direction on γ g which depends on the tangential velocity t · u with the friction coefficient C γ and friction function f γ . There is no friction on the wet boundary γ w .
For ice sheets that develop an ice shelf, the latter is assumed to be at hydrostatic equilibrium. In such a case, a calving front boundary condition (Schoof, 2007;van der Veen, 1996) is applied at γ d , in the form of the depth integrated stress balance (ρ w is the density of seawater) With (9), a calving rate u c can be determined at the ice front, see Sect. 3.2.1.

Adjoint equations
In this section, the adjoint equations are discussed, as emerging in a FS framework (Sect. 3.1) and in a SSA framework (Sect. 3.2), respectively. The adjoint equations follow from the Lagrangian of the forward equations after partial integration.
Lengthy derivations have been moved to Appendix A.

5
On the ice surface Γ s and over the time interval [0, T ], we consider the functional F We wish to determine its sensitivity to perturbations in both the friction coefficient C(x, t) at the base of the ice, and the topography b(x) of the base itself. We distinguish two cases: either u and h satisfy the FS equations (4), or the SSA equations (8). Given F, its forward solution (u, p, h) to (4) or (u, h) to (8), and its adjoint solution (v, q, ψ) or (v, ψ) to the adjoint FS and 10 adjoint SSA equations (both derived in the following and in Appendix A), we introduce a Lagrangian L(u, p, h; v, q, ψ; b, C).
Examples of F (u, h) in (10) are F = u−u obs 2 , and F = |h−h obs | 2 in an inverse problem, in which the task is to find b and C such that they match observations u obs and h obs at the surface Γ s , see also Gillet-Chaulet et al. (2016); Isaac et al. (2015); 15 Morlighem et al. (2013); Petra et al. (2012). Another example is F (u, h) = 1 T u 1 (x, t)δ(x − x * ) with the Dirac delta δ at x * to measure the time averaged horizontal velocity u 1 at x * on the ice surface Γ s with where T is the duration of the observation at Γ s . If the horizontal velocity is observed at ( The sensitivity in F and u 1 in (12) to perturbations in C and b is then given by (11) with the forward and adjoint solutions.
The same forward and adjoint equations are solved both for the inverse problem and the sensitivity problem but with different forcing function F . In the inverse problem, the change δF in F is of interest when δC and δb are changed during the solution procedure. In order to minimize F, δC and δb are chosen such that δF < 0 and F decreases. This procedure is repeated 25 iteratively until δF → 0 and F has reached a minimum.

Adjoint equations based on the FS model
In this section, we introduce the adjoint equations and the perturbation of the Lagrangian function. The detailed derivations of (13) and (16) below are given in Appendix A, starting from the weak form of the FS equations (4) on Ω, and by using integration by parts, and applying boundary conditions as in Martin and Monnier (2014); Petra et al. (2012).
The definition of the Lagrangian L for the FS equations is given in (A15) in Appendix A where (v, q, ψ) are the Lagrange 5 multipliers corresponding to the forward equations for (u, p, h). To determine (v, q, ψ), the following adjoint problem is solved: where the derivatives of F with respect to u and h are

10
Note that (13) consists of equations for the adjoint elevation ψ, the adjoint velocity v, and the adjoint pressure q.
cf. also Petra et al. (2012). For the rank four-tensor I, I ijkl = 1 only when i = j = k = l, otherwise I ijkl = 0. The operation 20 in (14) between a rank four-tensor A and a rank two-tensor (viz., a matrix) C is defined by (A C) ij = kl A ijkl C kl . In general, F b (Tu) in (13) is a linearization of the friction law f (Tu) in (4) with respect to the variable Tu. For instance, with a Weertman type friction law, f (Tu) = Tu m−1 , The perturbation of the Lagrangian function with respect to a perturbation δC in the slip coefficient C(x, t) involves the tangential components of the forward and adjoint velocities, Tu and Tv at the ice base Γ b , and is given by:

Time-dependent perturbations
Let us now investigate the effect of time-dependent perturbations in the friction parameter on modelled ice velocities and ice 5 surface elevation. A numerical example, based on the Marine Ice Sheet Model Intercomparison Project (MISMIP) (Pattyn et al., 2012) used also in Cheng and Lötstedt (2020) illustrates the findings presented.
Suppose that the velocity component u 1 * = u 1 (x * , t * ) is observed at (x * , t * ) at the ice surface as in (12) and that t * < T , then Above, we have introduced the sim-10 plifying notation that a variable with subscript * is a short-hand for it being evaluated at (x * , t * ), or, if it is time independent, at x * .
The procedure to determine the sensitivity is as follows. First, the forward equation (4) (13) becomes We start by investigating the response of ice velocities to perturbations in friction at the base: When the slip coefficient at the 20 ice base is changed by δC, then the change in u 1 * at Γ s is, according to (16), given by This implies that the perturbation δu 1 * mainly depends on δC at time t * and that contributions from previous δC(x, t), t < t * , are small. If we observe the horizontal velocity, then it responds instantaneously in time to the change in basal friction.
Further, to investigate the response of the ice surface elevation, h * at Γ s , to perturbations in basal friction, one considers The solution of the adjoint equation (13) In applied scenarios, friction at the base of an ice sheet is expected to exhibit seasonal variations. These can be expressed by δC(x, t) = δC 0 (x) cos(2πt/τ ), viz. a time dependent perturbation added to a stationary time average C(x), with 0 < δC 0 ≤ C. If, for illustrational purposes, τ = 1 (so, one year, from January to December), then Northern hemisphere cold and warm seasons can in a simplified manner be associated with nτ, n = 0, 1, 2, . . . (winter) and (n + 1/2)τ , n = 0, 1, 2, .... (summer).
Assume further that f (Tu)Tu · Tv is approximately constant in time. This is the case if u varies slowly. Then ψ ≈ const and 5 v ≈ const for t < t * . The change in ice surface elevation, δh, due to time-dependent variations in basal friction varies as Obviously, from the properties of the cosine function, the friction perturbation δC is large at t * = 0, τ /2, τ . . . , and vanishes at t * = τ /4, 3τ /4, . . .. Yet, (18) shows that δh * = 0 at t * = 0, τ . . . (so, during maximal friction in the winter) and at t * = τ /2, 3τ /2 . . . (so, during minimal friction in the summer), while δh * = 0 when δC = 0 at t * = τ /4, 3τ /4, . . . in the spring and 10 the fall. The response in h by changing C is delayed in phase by π/2 or in time by τ /4 = 0.25 yr. This is in contrast to the observation of u 1 in (17) where a perturbation in C is directly visible.
Particularly in an inverse problem where the phase shift between δC and δh in (18) is not accounted for, if h * is measured in the summer with δh(x, t * ) = 0, then the wrong conclusion would be drawn that there is no change in C.
To illustrate this, a two-dimensional numerical example is shown in Fig. 2, where the time scale and friction coefficient are 15 chosen as follows: τ = 1 yr, δC(x, t) = 0.01C cos(2πt) in x ∈ [0.9, 1.0] × 10 6 m. The ice sheet flows from x = 0 to L = 1.6 × 10 6 m on a single slope bed defined by b(x) = 720 − 778.5x 7.5×10 5 . u 1 and h are observed at x ∈ [0.85, 1.02] × 10 6 m. The steady state solution with the GL located at x GL = 1.035 × 10 6 m is perturbed by δu and δh when C is perturbed by δC. After perturbation, the GL position will oscillate. For additional details of the setup, we refer to the MISMIP (Pattyn et al., 2012) test case EXP 1 (step 1) used already in Cheng and Lötstedt (2020). The ice sheet is simulated by FS with Elmer/Ice (Gagliardini 20 et al., 2013) for 10 years. Fig. 2 shows that the perturbations δu 1 and δh in the grounded part of the ice sheet, specifically at x * = 0.85, 0.9, 0.95, 1.0, and 1.02×10 6 m for which individual panels are shown, oscillate regularly with a period of 1 year. The perturbations are small outside the interval [0.9, 1.0] × 10 6 . The initial condition at t = 0 is the steady state solution of the MISMIP problem and the FS solution with a variable C is essentially that steady state solution plus a small oscillatory perturbation, as in Fig. 2. 25 The weight f (Tu)Tu · Tv 0 in (17) is negative and an increase in the friction, δC > 0, leads to a decrease in the velocity and δC < 0 increases the velocity in all panels of Fig. 2. The velocities δu 1 and the surface elevations δh are separated by a phase shift in time, ∆φ = π/2, as predicted by (17) and (18). The weight in (18) for δC 0 in the integral over x changes sign when the observation point is passing from x * = 0.9 × 10 6 to 1.0 × 10 6 explaining why the shift changes sign in the red dashed lines shown in the two lower panels of Fig. 2.

The sensitivity problem and the inverse problem.
From a theoretical point of view, it is interesting to note that there is a relation between the sensitivity problem where the effect of perturbed parameters in the forward model is estimated and the inverse problem used to infer 'unobservable' parameters such as basal friction from observable data, e.g. ice velocity at the ice sheet surface. The same adjoint equations (13) are solved in both problems but with different driving functions defined by F (u, h) in (10).
. . , d, be the solution to (13) when u i is observed atx and F u = e i δ(x −x). The adjoint solution and the variation δF of the inverse problem can be expressed in It is shown in Appendix B that is a solution of (13) with arbitrary weights w i (x), i = 1, . . . , d, when When C is perturbed, the first variation of the functional in (16) is In the inverse problem in Petra et al. (2012), The weights in (20) whereṽ is a weighted sum of the solutions of the sensitivity problem v i over the whole domain ω The same construction of the solution is possible when h obs is given. Then d = 1, F (h) = 1 2 (h−h obs ) 2 , and F h = w = h−h obs . If we are interested in solving the inverse problem and determine δF in (21) in order to iteratively compute the optimal solution with a gradient method, then we solve (13) directly with F u = u − u obs or F h = h − h obs to obtainṽ without computing v i .

Steady state solution to the adjoint elevation equation in two dimensions.
A further theoretical consideration shows that the solution ψ to the adjoint elevation equation need not be computed to estimate 15 perturbations in the velocity for a two-dimensional vertical ice sheet at steady state. We show with the analytical solution in the FS model that the influence of ψ is negligible. It is sufficient to solve the adjoint stress equation for v to estimate the perturbation in the velocity.
The adjoint steady state equation in a two dimensional vertical ice in (13) is

20
The velocity from the forward equation is u(x, z) = (u 1 , u 3 ) T and the adjoint elevation ψ satisfies the right boundary condition The analytical solution ψ to (24) So, this solution has a jump −g( With a small h · u z (y) ≈ 0 in (25), an approximate solution is ψ( (25) and ψh ≈ 0 in (13). This is the case in the SSA of the FS model where u 1z (x) = 0 and in the SIA of the FS equations where u 1z (x, h) = 0 (Greve and Blatter, 2009;Hutter, 1983). When these approximations are accurate then u 1z will be small. Consequently, when u 1 is observed, the effect on v in the adjoint stress equation of the solution ψ of the adjoint advection equation in (13) 2012) yields an adequate answer. The situation is different when h is observed and ψ = 0 in (25).

Shallow shelf approximation
Starting from (8), a Lagrangian L of the SSA equations is defined, using the technique described and applied to the FS equations in Petra et al. (2012). By evaluating L at the forward solution (u, h) and at the adjoint solution (v, ψ), the adjoint SSA equations are obtained. Then, the effect of perturbed data at the ice base manifests itself at the ice surface as a perturbation δL; for details, see Appendix A. The adjoint SSA equations read: where the adjoint viscosityη and adjoint stressς are (cf. (14) for the case of FS) From (26) it is seen that the adjoint SSA equations have the same structure as the adjoint FS equations (13). There is one stress equation for the adjoint velocity v, and one equation for the Lagrange multiplier ψ corresponding to the surface elevation 20 equation in (8). However, the advection equation for ψ in (26) depends on v, implying a fully coupled system for v and ψ.
Equations (26) are solved backward in time with a final condition on ψ at t = T . As in (8), there is no time derivative in the stress equation. With a Weertman friction law, viz. f (u) = u m−1 and f γ (t · u) = |t · u| m−1 (cf. also Appendix A1), If the friction coefficient C at the ice base (both where it is grounded on bedrock (C > 0) and floating (C = 0)) is changed by δC, if the bottom topography is changed by δb, and if the lateral friction coefficient C γ is changed by δC γ , then it follows from Appendix A2 that the Lagrangian L is changed by (note that the weight in front of δC in (28) is actually the same as in (16)) The same perturbations in δC, δb, and δC γ could be allowed for the FS equations in (16) but because the FS equations are 5 more complicated than the SSA equations, the complexity of the derivation in the appendix and the expression for δL would increase considerably, which is why we refrain from considering them here.
Suppose that only h is observed with F h = 0 and F u = 0 in (26) where u u is the speed of the ice flux at x = 0 and u d is the calving speed at x = L. If x = 0 is at the ice divide, then u u = 0.
By the stress balance (9), the calving front satisfies Assuming that u > 0 and u x > 0, the viscosity becomes η = 2A − 1 n u 1−n n x , and the friction term with a Weertman law turns into Cf (u)u = Cu m . The adjoint equations for v and ψ follow either from simplifying (26), or from (29) and read as follows: Note that the viscosity above is multiplied by a factor 1/n, n > 0 which represents an extension of the adjoint SSA in MacAyeal (1993) where n = 1. The effect on the Lagrangian of perturbations δb and δC is obtained from (28) The weights or sensitivity functions w b and w C multiplying δb and δC in the integral are defined by 3.2.2 The two-dimensional forward steady state solution.
We turn now to a discussion of steady-state solutions to the system (30). Except from letting all time derivatives vanish, even the 10 longitudinal stress can be ignored in the steady state solution, see Schoof (2007). Moreover, the viscosity terms in (30) and the adjoint equations (30) reduce, under the assumption that the basal topography is characterized by a small spatial The solution to the forward equation (33) is presented for the case when a and C are constant, (for details, see (D3) and (D4) in Appendix D): The solution is calibrated with the ice thickness H GL = H(x GL ) at the fixed grounding line.

3.2.3
The two-dimensional adjoint steady state solution with F u = 0.
In this section, the analytical solution to the adjoint equation (34) is discussed. The derivation of the solution is detailed in 15 Appendix E to Appendix F. It is here sufficient to recall that the below given solution is derived under the assumptions that b x H x , and that a and C are constants.
For observations of u at x * ,   or, after scaling with u * : The weights w b and w C in (37) multiplying δb and δC, defined in the same manner as in (31) and (32), are shown in Fig. 6(a) and Fig. 7(a) with the solutions in Fig. 4(a) and Fig. 5(a). The Dirac term is plotted as a vertical line at x * in Fig. 7(a).
All perturbations in C between x * and x GL will result in a perturbation of the opposite sign in u * at the surface because 5 w C < 0 in (x * , x GL ) in Fig. 6(a) and (37). The same conclusion holds true for perturbations in b because w b < 0 in (x * , x GL ) in Fig. 7(a) but an additional contribution is added from δb at x * by the Dirac delta in w b . A perturbation is less visible in u the farther away from x GL the observation point is since the amplitude of both w C and w b decays when x * decreases.
The relation in (38) between the relative perturbations δb/H, δC/C and δu/u can also be interpreted as a way to quantify the uncertainty in u. An uncertainty δC in C and an uncertainty δb in b is propagated to an uncertainty δu * in u at x * by (38).

10
The following conclusions can be drawn from (37) and (38) and Figs. 6 and 7: (i). The closer perturbations in basal friction are located to the grounding line, the larger perturbations of velocity will be observed at the surface. This is because the weight in front of δC increases when x * → x GL , see Fig. 6, which in turn is an effect of the increasing velocity u * and the decreasing thickness H * , as the grounding line is approached, see Fig. 3.
Or, compactly expressed, δC with support in [x * , x GL ] will cause larger perturbations at the surface the closer x * is to 15 x GL , and the closer δC(x) is to x GL .
(ii). Variations in the observed velocity δu * at the surface at observation point x * will include contributions from changes in the frictional parameter, δC, between x * and the grounding line x GL , and from changes in basal topography, δb, but it is impossible to disentangle their individual contributions to δu * . (iii). When the variation in ice thickness is small compared to the overall ice thickness, H x H, a small perturbation in basal topography δb is directly visible in the surface velocity. This is because in such a case, δu * ≈ u * δb * /H * and the main effect on u * from the perturbation δb is localized at each x * , see (37).
(iv). For an unperturbed basal topography, two different perturbations of the friction coefficient will result in the same perturbation of the velocity. In other words: the perturbation δC cannot be uniquely determined by one observation of δu. This The same δu * is observed with a constant perturbation in 2 ).

10
(v). A rapidly varying friction coefficient at the base of the ice sheet will be difficult to identify by observing the velocity at the ice surface. In contrast, a smoothly varying friction coefficient at the base will be easily observable at the ice sheet surface. This is seen as follows: Perturb C by δC = cos(kx/x GL ) in (37) for some wave number k which determines the smoothness of the friction at the bedrock and amplitude and let δb = 0 and m = 1. The wavelength is λ = 2πx GL /k.
When k is small then the wavelength is long and the variation of C + δC is smooth. When k is large then the friction 15 coefficient varies rapidly in x with a short λ. The perturbation in the velocity is For a thin ice with a small H * , a perturbation in C is easier to observe at the surface than for a thick ice. When k grows at the ice base, the amplitude of the perturbation at the ice surface decays as 1/k. Thus, the effect of high wave number perturbations of C will be difficult to observe at the top of the ice but smooth perturbations at the base will propagate to 20 the surface. If k is large and the surface velocity is of interest in a numerical simulation, then there is no reason to use a fine mesh at the base to resolve the fast variation in C because it will not be visible at the top of the ice.
(vi). A perturbation in the topography with long wavelength is easier to detect at the surface than a perturbation with short wavelength. If δC = 0 and b is perturbed by δb = cos(kx/x GL ), then any perturbation at x * is propagated to the surface by u * δb * H * , which is the first term on the right hand side of (37). The effect is larger if the ice is thin and moving fast. 25 The integral term will behave in the same way as in (40) with mainly perturbations with small wave numbers and long wavelengths visible at the surface.
Finally, let us comment on other approaches to investigate the sensitivity of surface data to changes in b and C, e.g. using three linear models as in Gudmundsson (2008) and along a flow line at steady state in Gudmundsson and Raymond (2008) with a linearized FS model with n = 1 and m = 1. In these papers, transfer functions for the perturbations from base to surface corresponding to our formulas (37) and (38) are derived by Fourier and Laplace analysis. The perturbations with long wavelength λ and small wave number k are propagated to the surface but short wavelengths are effectively damped in Gudmundsson (2008). The transfer functions are utilized in Gudmundsson and Raymond (2008) to estimate how well basal data can be retrieved from surface data. Retrieval of basal slipperiness C is possible for perturbations δC of long wavelength 5 and if the errors in the basal topography δb is small. Short wavelength perturbations δb can be determined from surface data.
The same conclusions as in Gudmundsson (2008) and Gudmundsson and Raymond (2008) can be drawn from our explicit expressions for the dependence of δu * and δh * on δC and δb. For example, it follows from (40) that only δC with a long wavelength is visible at the surface and that δb also with a short wavelength affects δu * in (38). If δb is small or zero in (38), then it is easier to determine the δC that causes a certain δu * . In the case when h is observed at x * and F u = 0 and F h = δ(x − x * ), the expressions for ψ and v satisfying (34) are The corresponding formulas when u is observed are found in (36). There is a discontinuity at the observation point x * in v(x) illustrated in Fig. 5(b), but ψ(x) is continuous in the solution of (34) and in Fig. 4(b). 15 The second derivative term ( 1 n ηHv x ) x is neglected in the simplified equation (34) but is of importance at x * . A correctionψ of ψ at x * in (41) is therefore introduced to satisfy 1 n ηHv x x −Hψ Fig. 4(b) withψ as a vertical line representing the negative Dirac delta.
As in (37), (42) is rewritten with the weights w b and w C in (32) These weights are shown in Fig. 6(b) and Fig. 7(b). The negative derivative of the Dirac delta is depicted in Fig. 7(b) as a vertical line in the negative direction immediately followed by one in the positive direction.
The contribution from the integrals in (38) and (42) is identical except for the sign (compare w C in Fig. 6(a) and Fig. 6(b) 5 and w b in Fig. 7(a) and Fig. 7(b)). The first term in (38) (38) and (42) (37) and (38) for δu * are valid also for δh * in (42). Finally, the time dependent adjoint equation (30)  year, 1 yr, where the beginning of the year is associated with winter) in the forward equation (29): C(x, t) = C 0 (1 + κ cos(2πt)), 0 < κ < 1. 15 Apparently, C has its highest value at t = n, n = 0, 1, 2, . . ., i.e. the winter, and its lowest value at t = n + 1/2, i.e. the summer, as in Fig. 2. The amplitude of the perturbation is set to κ = 0.5 and the forward equation (29) is solved for 11 years. The topography b is kept constant in time. Observations of u and h are taken at x * = 9 × 10 5 m for 0.1 a in the four seasons starting from the summer of the tenth year, e.g., in the summer (t * = 9.5), the fall (t * = 9.75), the winter (t * = 10), and the spring (t * = 10.25). The forward equations (29) are solved numerically from t = 0 with the steady state solution as initial data to 20 the observation points t = t * and the adjoint equations (30) are solved from t = t * backward in time to t = 0. According to a convergence test, the time step is chosen to be 0.01 yr and the spatial resolution is 10 3 m. A visual inspection of the computed solutions after halving the step sizes indicates that a sufficiently converged solution has been reached. Most of the weights in space and time are negligible implying that perturbations in those domains are not visible at (x * , t * ).
Only δC and δb in a narrow interval around x * for t in [0, t * ] have an influence on δu * and δh * . Therefore, we take a snapshot of the x axis (or abscissa) with the width of 10 5 m in space around x * in Fig. 8. The weights oscillate in time because of the 30 seasonal variation in the basal conditions. A perturbation at the base is propagated to the x * position on the surface but with a possible delay in time. The earlier a perturbation in C or b takes place in the interval [0, t * ), the smaller the effect of it is at t * .
After five years a perturbation can hardly be detected at the surface.
The temporal variations of the adjoint weights at x * in Fig. 8 are shown in Fig. 9 for the four seasons with four different colors.
As expected, the weights vanish when t > t * . In Fig. 9(a) and (b), the perturbations δC * and δb * have a direct effect on δu * at t * , where both w C (x * , t * ) and w b (x * , t * ) are negative. The same direct effect of δC is found for δu 1 * solving the FS equations 5 (17) in Sect. 3.1.1. A change in δC * at the base is observed immediately as a change in u at the surface. The effect of δC on δu * for t < t * is weak in Fig. 9(a), i.e. the memory of old perturbations is short. The largest effect of δC on δu * and δh * appears with t * in the summer when C is small in (44) (the blue lines in Fig. 9(a) and (b)).
However, when h is observed, the effects of δC * and δb * are not visible directly because w C * ≈ 0 and w b * ≈ 0 in Fig. 9(c) and Fig. 9(d). An intuitive explanation is that there is an immediate effect on the velocity but there is a delay in h since it 10 is integrated in time from the velocity field. Additionally, the effects of δC and δb are difficult to separate, since the weight w b (x * , t) has a shape similar to w C (x * , t). The largest effect on δh * is from δC in the summer due to the peaks in w C in Fig. 9(c). For the same δC, the largest δh * is observed in the fall (orange), then the second largest δh * is in the winter (green) followed by the spring observation (red). If δh * is observed in the fall and the time dependency is ignored, then the wrong conclusion is drawn that δC in the fall has the strongest effect (but it is the summer perturbation). There is a delay in time Suppose that the temporal perturbation is oscillatory δC 0 cos(2πf t) with frequency f . A low frequency f with f 1 corresponds decennial or centennial variations and a high frequency with f 1 corresponds to diurnal or weekly variations. Then the perturbation in h at t = t * is 25 cf. (40). With a high frequency, f 1, then δh * ∝ 1/f and high frequency perturbations are damped efficiently. If the frequency is low, f 1, then δh * ∝ τ and the change in h * is insensitive to the frequency. The same conclusions hold true for δb where decennial perturbations seem more realistic.

Conclusions
The perturbation of the basal condition at x * has the strongest impact at x * on the surface, possibly with a time delay. Such a time delay occurs when a perturbation at the ice base is visible at the surface in h, but in u it is observed immediately (Fig. 9). The effect of a perturbation disappears more quickly, the older the perturbation is.

10
Perturbations in the friction coefficient at the base observed in the surface velocity determined by SSA are damped inversely proportional to the wave number and the frequency of the perturbations in (40) and (45), thus making very oscillatory perturbations in space and time difficult to register at the ice sheet surface. In such a case, there is no need to have a fine mesh and a small time-step in a numerical solution to resolve the rapid oscillations in C at the base.
Appendix A: Derivation of the adjoint equations 15

A1 Adjoint viscosity and friction in SSA
The adjoint viscosityη(u) in SSA in (14) is derived as follows. The SSA viscosity for u and u + δu is Determine B(u) such that Then use the operator to define B 1−n 2nη

Thus, let
or in tensor form Replacing B in (A2) by D we obtain the adjoint FS viscosity in (14).
(A6) 15 Terms of order two or more in δL are neglected. Then the first term in δL satisfies Using partial integration, Gauss' formula, and the initial and boundary conditions on u and H and ψ(x, T ) = 0, x ∈ ω, and ψ(x, t) = 0, x ∈ γ w , in the second integral we have The first integral after the second equality vanishes since h is a weak solution and I 2 is Using the weak solution of (8), the adjoint viscosity (27), (A2), the friction coefficient (A3), Gauss' formula, the boundary conditions, and neglecting the second order terms, the third and fourth integrals in (A5) are I 3 = I 31 + I 32 , where Collecting all the terms in (A7), (A9), and (A10), the first variation of L is The forward solution (u * , p * , h * ) and adjoint solution (v * , q * , ψ * ) satisfying (8) and (26) are inserted into (A4) resulting in Then (A12) yields the variation in L in (A13) with respect to perturbations δb, δC γ , and δC in b, C γ , and C

A3 Adjoint equations in FS
The FS Lagrangian is In the same manner as in (A5), the perturbed FS Lagrangian is 10 L(u + δu, p + δp; v + δv, q + δq; h + δh, ψ + δψ; C + δC) = L(u, p, h; v, q, ψ; C) + I 1 + I 2 + I 3 . (A16) Terms of order two or more in δu, δv, δh are neglected. The first integral I 1 in (A16) is Partial integration, the conditions ψ(x, T ) = 0 and ψ(x, t) = 0 at Γ s , and the fact that h is a weak solution simplify the second integral Define Ξ, ξ, and Υ to be Then a weak solution, (u, p), for any (v, q) satisfying the boundary conditions, fulfills The third integral in (A16) is I 3 = I 31 + I 32 , The integral I 31 is expanded as in (A10) and (A11) or Petra et al. (2012) using the weak solution, Gauss' formula, and the 10 definitions of the adjoint viscosity and adjoint friction coefficient in Appendix A1. When b < z < h we have Υ(u, p; v, q) = 0.
If Υ is extended smoothly in the positive z-direction from z = h, then with z ∈ [h, h + δh] for some constant c > 0 we have |Υ| ≤ cδh. Therefore, where |ω| is the area of ω. This term is a second variation in δh which is neglected and I 3 = I 31 .
Consider a target functional F for the steady state solution with a weight vector w(x) with components w i (x) multiplying δu i in the first variation of F. Using (16), δF is Appendix C: Steady state solution of the adjoint height equation in the FS model In a two dimensional vertical ice with u(x, z) = (u 1 , u 3 ) T , the stationary equation for ψ in (13) is When x > x * , where F h = 0 and F u = 0, we have ψ(x) = 0 since the right boundary condition is ψ(L) = 0.
If u 1 is observed at Γ s then F (u, h) = u 1 (x)χ(x) and F u = (χ(x), 0) T and F h = 0. The weight χ on u 1 may be a Dirac delta, a Gaussian, or a constant in a limited interval. On the other hand, if F (u, h) = h(x)χ(x) then F h = χ(x) and F u = 0.

Appendix D: Simplified SSA equations
The forward and adjoint SSA equations in (33) and (34) The equation for H m+2 for x ≤ x GL is integrated from x to x GL such that By including the viscosity term in (29) and assuming that H(x) is linear in x, a more accurate formula is obtained for u(x) on the floating ice in (6.77) of Greve and Blatter (2009).
Since H is continuous and u and v are bounded, when x − * → x + * , then A similar relation for ψ can be derived With F u = 0 and F h = 0 for x < x * and v(0) = ψ x (0) = 0, we find that If F (u, h) = uδ(x − x * ), then by (E5) and (E6) Let H(x − x * ) = x−x * −∞ δ(s) ds be the Heaviside step function at x * . Then To satisfy the jump condition in (E8) and (E9), the constant C v is 15 Combine (F1) with the relation ψ x = (F u − Cmu m−1 v)/H and integrate from x to x GL to obtain With the jump condition in (E8) and (E9), ψ(x) at 0 ≤ x < x * is