Sensitivity of ice sheet surface velocity and elevation to variations in basal friction and topography in the full Stokes and shallow-shelf approximation frameworks using adjoint equations

Predictions of future mass loss from ice sheets are afflicted with uncertainty, caused, among others, by insufficient understanding of spatiotemporally variable processes at the inaccessible base of ice sheets for which few direct observations exist and of which basal friction is a prime example. Here, we present a general numerical framework for studying the relationship between bed and surface properties of ice sheets and glaciers. Specifically, we use an inverse modeling approach and the associated time-dependent adjoint equations, derived in the framework of a full Stokes model and a shallow-shelf/shelfy-stream approximation model, respectively, to determine the sensitivity of grounded ice sheet surface velocities and elevation to time-dependent perturbations in basal friction and basal topography. Analytical and numerical examples are presented showing the importance of including the time-dependent kinematic free surface equation for the elevation and its adjoint, in particular for observations of the elevation. A closed form of the analytical solutions to the adjoint equations is given for a two-dimensional vertical ice in steady state under the shallow-shelf approximation. There is a delay in time between a seasonal perturbation at the ice base and the observation of the change in elevation. A perturbation at the base in the topography has a direct effect in space at the surface above the perturbation, and a perturbation in the friction is propagated directly to the surface in time.


Introduction
Over the last decades, ice sheets and glaciers have experienced mass loss due to global warming, both in the polar regions and also outside of Greenland and Antarctica (Farinotti et al., 2015;Mouginot et al., 2019;Pörtner et al., 2019;. The most common benchmark date for which future ice sheet and glacier mass loss and associated global mean sea level rise is predicted is the year 2100 CE, but recently, even the year 2300 CE and beyond are considered (Pörtner et al., 2019;Steffen et al., 2018). Global mean sea level rise is predicted to continue well beyond 2100 CE in the warming scenarios typically referred to as RCPs (representative concentration pathways; see van Vuuren et al., 2011), but rates and ranges are afflicted with uncertainty, caused by, among others, insufficient understanding of spatiotemporally variable processes at the inaccessible base of ice sheets and glaciers (Pörtner et al., 2019;Ritz et al., 2015). These include the geothermal heat regime, subglacial and base-proximal englacial hydrology, and particularly the sliding of the ice sheet and glaciers across their base, for which only few direct observations exist (Fisher et al., 2015;Key and Siegfried, 2017;Maier et al., 2019;Pattyn and Morlighem, 2020).
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 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), which, once solved, render, e.g., ice velocities part of the solution. Studies of frictional models and resulting sliding laws for glacier and ice sheet flow emerged in the 1950ssee, e.g., Fowler (2011), Iken (1981), Lliboutry (1968), Nye (1969), Schoof (2005), and Weertman (1957) -have subsequently been implemented into numerical models of ice sheet and glacier behavior -e.g., in Brondex et al. (2017), Brondex et al. (2019), Gladstone et al. (2017), Tsai et al. (2015), Wilkens et al. (2015), and 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 few or no observational data are 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 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), the geothermal heat flux at the ice base (Zhu et al., 2016), or to estimate basal topography beneath an ice sheet (Monnier and des Boscs, 2017;van Pelt et al., 2013). The latter is not only difficult to separate from the determination of the 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 and 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 is computed by solving the so-called adjoint equations to the forward equations, where the latter often are slightly simplified, such as 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 (steadystate) 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 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). A related problem is to find the relation between perturbations of the friction parameters and the topography at the base and the resulting perturbations of the ice velocity and the elevation at the surface. The parameter sensitivity of the solutions is studied in Gudmundsson (2008) and Gudmundsson and Raymond (2008) by linearizing the governing equations. They show that basal perturbations of short wavelength are damped when they reach the surface, but the effect of long wavelengths can be observed at the top.
Here, we present an analysis of the sensitivity of the velocity field and the elevation of the surface of a dynamic, grounded ice sheet (modeled by both FS and SSA, briefly described in Sect. 2) to perturbations in the sliding parameters contained in Weertman's law (Weertman, 1957) and the topography at the ice base. The perturbations in a velocity component or the ice surface elevation at a certain location and time are determined by the solutions to the forward equations and the associated adjoint equations of a grounded ice sheet. A certain type of perturbation at the base may cause a very small perturbation at the top of the ice. Such a basal perturbation will be difficult to infer from surface observations in an inverse problem. High-frequency perturbations in space and time are examples when little is propagated to the surface. This is also the conclusion in Gudmundsson and Raymond (2008), derived from an analysis differing from the one presented here. The adjoint problem that is solved here to determine this sensitivity (Sect. 3) goes beyond similar earlier works by MacAyeal (1993) and Petra et al. (2012) because it includes the time-dependent advection equation for the kinematic free surface. The key concepts and steps introduced in Sect. 3 are supplemented by detailed derivations, collected in Appendix A. The same adjoint equations are applicable in the inverse problem to compute the gradient of the objective function and to quantify the uncertainty in the surface velocity and elevation due to uncertainties at the ice base. Examples of uncertainties are measurement errors in the basal topography and unknown variations in the parameters in the friction model. Analytical solutions in two dimensions of the stationary adjoint equations subjected to simplifying assumptions are presented, from which the dependence of the parameters, on, e.g., friction coefficients and bedrock topography, becomes obvious. The time-dependent adjoint equations are solved numerically, and the sensitivity to perturbations varying in time is investigated and illustrated.
The sensitivity of the surface velocity and elevation to perturbations in the friction and topography is quantified in extensive numerical computations in a companion paper (Cheng and Lötstedt, 2020). The adjoint equations derived and studied analytically in this paper are solved numerically for the FS and SSA models in Cheng and Lötstedt (2020) and compared to direct calculations of the surface perturbations with the forward equations. Discrete transfer functions are computed and analyzed as in Gudmundsson (2008) for the relation between surface perturbations and basal perturbations. While Gudmundsson's analysis is based on Fourier analysis, the analysis in Cheng and Lötstedt (2020) relies on analytical solutions of the SSA equations. The analysis in this paper and the numerical experiments in Cheng and Lötstedt (2020) confirm the conclusion in Gudmundsson (2008) that the perturbations with a long wavelength and low frequency will propagate to the surface while those of a short wavelength and high frequency are damped.

Ice models
In this section, the equations emerging from adopting a forward approach to describing ice dynamics are presented, together with relevant boundary conditions, for the FS (4) and SSA model (9), respectively. These, and the notation and terminology introduced here, provide the framework in which the adjoint equations are discussed in Sect. 3.
The flow of large bodies of ice is described with the help of the conservation laws of mass, momentum, and energy (Greve and Blatter, 2009), which together pose a system of nonlinear partial differential equations (PDEs) commonly referred to as the FS equations in glaciological applications. In the FS equations, nonlinearity is introduced through the viscosity in Glen's flow law, a constitutive relation between strain rates and stresses (Glen, 1955). Continental sized ice masses (ice sheets and, if applicable, their floating extensions known as ice shelves) are shallow in the sense that their vertical extension V is orders of magnitudes smaller than their horizontal extension L, such that the aspect ratio V /L is on the order 10 −2 to 10 −3 . The aspect ratio is used to introduce simplifications to the FS equations, resulting, e.g., in the shallow-ice (SIA) (Hutter, 1983), shallowshelf (Morland, 1987), and shelfy-stream (MacAyeal, 1989) approximations, parts of which can be coupled to FS using ISCAL (Ahlkrona et al., 2016), a method applicable to any FS framework and implemented in Elmer/Ice. They are all characterized by substantially reduced computational costs for numerical simulation, compared to using the FS model. A common simplification, also adopted in our analysis in the following, is the assumption of isothermal conditions, which implies that the balance of energy need not be considered.
The upper surface of the ice mass and also the ice-ocean interface constitute a moving boundary and satisfy an advection equation describing the evolution of its elevation and location (in response to mass gain, mass loss, or/and mass advection). For ice masses resting on bedrock or sediments, sliding needs to be parameterized at the interface. The interface between floating ice shelves and sea water in the ice shelf cavities is usually regarded as frictionless.

Full Stokes model
We adopt standard notation and denote vectors in bold italics and matrices in three-dimensional space in bold, and we denote derivatives with respect to the spatial coordinates and time by subscripts x, y, z, and t. The horizontal plane is spanned by the x and y coordinates, and z is the coordinate in the vertical direction (see Fig. 1a). 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, γ GL , z b = b, and downstream of γ GL we have z b > b (see Fig. 1). In two dimensions, γ GL consists of one point with x coordinate x GL .
The boundary enclosing the domain occupied by the ice has different parts (see Fig. 1b). The lower boundaries of 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 is denoted by s (ice surface) at h(x, y, t) in Fig. 1a. The footprint (or projection) of in the horizontal plane is denoted by ω and its boundary is γ .
The vertical lateral boundary (in the x-z plane, Fig. 1b) 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 . For a two-dimensional flow line geometry (Fig. 1a), this corresponds to x = (x, z) T , ω = [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 s , and gravitational acceleration by g. The values of these physical parameters are given in Table 1. 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 Tσ n = −βTu, where 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 The friction coefficient has a general form β(u, x, t) = C(x, t)f (u), where C(x, t) is independent of the velocity u and f (u) represents some linear or nonlinear function of u.
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 h 0 (x) > b(x, y) is the initial surface elevation and h γ (x, t) is a given height on the inflow boundary. In the floating ice, the equations and the boundary conditions are as above plus an equation for the free surface z b with z b ≥ b and a boundary condition on the wetted ice w : where a b is the basal mass balance and p w is the water pressure at w . With the sea level at z = 0 and the water density ρ w , p w = −ρ w gz b . The solution at the grounding line satisfies a nonlinear complementarity problem. Let d be the distance between the ice base and the bedrock such that On w , we have d(x, t) > 0 and n · σ n + p w = 0 and on b , we have d(x, t) = 0 and n · σ n + p w < 0. The complementarity relation at the ice base b ∪ w can be summarized as d(x, t) ≥ 0, n · σ n + p w ≤ 0, d(x, t)(n · σ n + p w ) = 0.
There are additional constraints on the solution. For example, the thickness of the ice is non-negative H ≥ 0, there is a lower bound on the velocity in Weertmann's friction law, and there is an upper bound on the viscosity η. These conditions have to be handled in a numerical solution of the equations but are not discussed further here.
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 ice velocity at the calving front is defined as u d to simplify the analysis. The vertical component of σ n vanishes on u .

Shallow-shelf approximation
The three-dimensional FS problem (4) in can be reduced 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 (MacAyeal, 1989). The SSA is often also used in regions of fast flow over lubricated bedrock (MacAyeal, 1989). 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 (7), the Frobenius inner product between two matrices A and C is used, defined by A : C = ij A ij C ij . The vertically integrated stress tensor ς (u) (cf. 3 for the FS model) is given by where H is the ice thickness (see Fig. 1). 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 β(u, x, t) = C(x, t)f (u) = C(x, t) u m−1 with a friction coefficient C(x, t) ≥ 0, just as in the FS model. In summary, the forward equations describing the evolution of the ice surface and ice velocities based on an SSA model (in which u is not divergence-free) read Above, t is the tangential vector on γ = γ u ∪ γ d such that n · t = 0 and a = a s − a b . 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 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 . When the equations are solved for the grounded part, H t = h t with a time-independent topography.
A flotation criterion determines the position of the grounding line (see, e.g., Seroussi et al. (2014)). The thickness H is compared to a flotation height H f given by If H > H f then the ice is grounded and C > 0 in (9). If H < H f then it is floating with C = 0. The grounding line is defined by H = H f .
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):

The flow line model of SSA
In this section, the SSA equations are presented for the case of an idealized, two-dimensional vertical sheet in the xz plane (see Fig. 1). The forward SSA equations are derived from (9) by letting H and u 1 be independent of y and setting u 2 = 0. Since there is no lateral force, C γ = 0. The position of the grounding line is denoted by x GL , and . Basal friction C is positive and constant where the ice sheet is grounded on bedrock, while C = 0 at the floating ice shelf's lower boundary. To simplify notation, we let u = u 1 . The forward equations thus become where u u is the speed of the ice flux at x = 0 and u d is the speed at the calving front at x = L. If x = 0 is at the ice divide, then u u = 0. By the stress balance (10), 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 two-dimensional forward steady-state solution
We now discuss the steady-state solutions to the system (11). Except for letting all time derivatives vanish, even the longitudinal stress can be ignored in the steady-state solution (see Table 1. The model parameters.
Parameter Quantity Rate factor of Glen's flow law C = 7.624 × 10 6 Pa m −1/3 s 1/3 Basal friction coefficient g = 9.8 m s −2 Acceleration of gravity m = 1/3 Friction law exponent n = 3 Flow-law exponent ρ = 900 kg m −3 Ice density Schoof (2007)). With a sliding law in the form f (u) = u m−1 and the thickness given at x GL , (11) thus reduces to The solution to the forward equation (12) is derived for the case when a and C are constant (for details, see D3 and D4 in Appendix D): The solution is normalized with the ice thickness H GL = H (x GL ) at the grounding line. Both x GL and H GL are assumed to be known in the formula. Similar equations to those in (12) are derived in Nye (1959) using the properties of large ice sheets. A formula for H (x) resembling (13) and involving H (0) is the solution of the equations. By including the longitudinal stress in the ice, an approximate, more complicated expression for H (x) is obtained in Weertman (1961). Figure 2 displays solutions from (13) obtained with data from the Marine Ice Sheet Model Intercomparison Project (MISMIP) (Pattyn et al., 2012) test case EXP 1 chosen in Cheng and Lötstedt (2020). The modeling parameters in (13) are given in Table 1. 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 and lifts from it at the grounding line position x GL = 1.035×10 6 m. As x approaches x GL , H decreases to become H GL at x GL in Fig. 2b.
The larger the friction coefficient C and accumulation rate a are, the steeper the decrease in H is in (13). The numerator in u increases and the denominator decreases when x → x GL , resulting in a rapid increase in u. The MISMIP example is such that the SSA solution is close to the FS solution. Numerical experiments in Cheng and Lötstedt (2020) show that an accurate solution compared to the FS and SSA solutions is obtained with u and H in (13) solving (12).
Finally, it is noted that an alternative solution to (11) valid for the floating ice shelf, x > x GL , but under the restraining assumption of H (x) being linear in x, is found in Greve and Blatter (2009).

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. We derive the adjoint equations for the grounded part of the ice sheet. There the friction coefficient can be perturbed with δC = 0, and a perturbation δb of the topography has a direct influence on the flow of ice. The adjoint equations follow from the Lagrangian based on the forward equations after partial integration. Lengthy derivations have been moved to Appendix A. A numerical example based on the MISMIP (Pattyn et al., 2012) used also in Cheng and Lötstedt (2020) illustrates the findings presented.
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 of the base itself b(x), which is a smooth function in x. We distinguish two cases: u and h satisfy either the FS equation (4) or the SSA equation (9). Given F, the forward solution (u, p, h) to (4) or (u, h) to (9), and the adjoint solution (v, q, ψ) or (v, ψ) to the adjoint FS and adjoint SSA equations (both derived in the following and in Appendix A), we introduce a Lagrangian L(u, p, h; v, q, ψ; b, C). The Lagrangian for the FS equations is with the Lagrange multipliers v, q, and ψ corresponding to the forward equations for u, p, and h. The effect of perturbations δC and δb in C and b on F is given by the perturbation δL, viz. Examples of F (u, h) in (14) 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;Morlighem et al., 2013;and Petra et al., 2012).
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 (17) or (18) to perturbations in C and b is then given by (16) 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 (14). If we are interested in how u 1 changes at the surface when b and C are changed at the base by given δb and δC, then F is as in (18). The forward and adjoint equations are then solved once. In the inverse problem with velocity observations, F is the objective function in a minimization procedure and F = u − u obs 2 . 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 with C + δC and b + δb and u is closer to u obs . Precisely how δC and δb are chosen depends on the optimization algorithm. This procedure is repeated iteratively, and b and C are updated by b + δb and C + δC until δF → 0 and F has reached a minimum. The forward and adjoint equations have to be solved in each iteration in the inverse problem.

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 (19) and (22) below are given in Appendix A, starting from the weak form of the FS equation (4) on and by using integration by parts, and applying boundary conditions as in Martin and Monnier (2014) and Petra et al. (2012).
The definition of the Lagrangian L for the FS equations is given in (15) and (A15) in Appendix A. To determine (v, q, ψ) in (15), the following adjoint problem is solved: where the derivatives of F with respect to u and h are Note that (19) consists of equations for the adjoint elevation ψ, the adjoint velocity v, and the adjoint pressure q. The equations are the same as when the derivatives are computed in the inverse problem except for the terms depending on F , which is the misfit between the numerical solution and the observation in the inverse problem. Compared to the steadystate adjoint equation for the FS equations discussed in Petra et al. (2012), an advection equation is added in (19) for the Lagrange multiplier ψ(x, t) on s with a right-hand side depending on the observation function F and one term depending on ψ in the boundary condition on s . The adjoint elevation equation for ψ can be solved independently of the adjoint stress equation since it is independent of v. If h is observed and F h = 0 and F u = 0, then the adjoint elevation equation must be solved together with the adjoint stress equation. Otherwise, the term ψh is ignored in the right-hand side of the boundary condition of the adjoint stress equation and the solution is v = 0 with δF = 0 in (22); see below.
The adjoint viscosity and adjoint stress arẽ cf. also Petra et al. (2012). For the rank four-tensor I, I ij kl = 1 only when i = j = k = l; otherwise I ij kl = 0. The operation in (20) between a rank four-tensor A and a rank two-tensor (viz. a matrix) C is defined by (19) 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, 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 For this formula to be accurate, δC has to be small. Otherwise, nonlinear effects may be of importance.

Time-dependent perturbations
Let us now investigate the effect of time-dependent perturbations in the friction parameter on modeled ice velocities and ice surface elevation. Suppose that the velocity component at the ice surface as in (18) and that t * < T , then with Above, we have introduced the simplifying notation that a variable with subscript * is a shorthand for it being evaluated at (x * , t * ) or, if it is time independent, at x * . Here we have chosen to consider the perturbation at a certain point in space and time (x * , t * ), which is sufficient because other types of sensitivity over a certain period of time and space as in (17) are the linear combination of point-wise sensitivities. The procedure to determine the sensitivity is as follows. First, the forward equation (4) is solved for u(x, t) and h(x, t) from t = 0 to t = T . Then, the adjoint equation (19) is solved backward in time (from t = T to t = 0) with ψ(x, T ) = 0 as the corresponding final condition. Obviously, the solution for t * < t ≤ T is ψ(x, t) = 0 and v(x, t) = 0. Letting e i denote the unit vector with 1 in the ith component, the boundary condition in (19) We start by investigating the response of ice velocities to perturbations in friction at the base: when the slip coefficient at the ice base is changed by δC, then the change in u 1 * at s is, according to (22), 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 (19) 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 timedependent perturbation added to a stationary time average C(x), with 0 < δC 0 ≤ C. If, for illustrational purposes, τ = 1 (1 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 in time. Then ψ ≈ constant and v ≈ constant for t < t * . The change in ice surface elevation, δh, due to time-dependent variations in basal friction varies as where the integral J is 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, (25) shows that δh * = 0 at t * = 0, τ, . . . (during maximal friction in the winter) and at t * = τ/2, 3τ/2, . . . (during minimal friction in the summer), while δh * = 0 when δC = 0 at t * = τ/4, 3τ/4, . . . in the spring and the fall. The response in h by changing C is delayed in phase by π/2 or in time by τ/4 = 0.25 years. This is in contrast to the observation of u 1 in (24) where a perturbation in C is directly visible. Particularly in an inverse problem where the phase shift between δC and δh in (25) is not accounted for, if h * is measured in the summer with δh(x, t * ) = 0, then the wrong conclusion would be drawn such that there is no change in C.
In another example, suppose that there is an interval with a step change of C with δC(x, t) = δC 0 (x)s(t), where s(t) = 1 in the time interval [t 0 , t 1 ] and 0 otherwise. Then with J in (26), δh * in (25) is The effect of the basal perturbation successively increases in the elevation when t * > t 0 and stays at a higher level for t * > t 1 when δC has vanished.

Example with seasonal variation
To illustrate the phase delay in an oscillatory perturbation, a two-dimensional numerical example is shown in Fig. 3, where the timescale and friction coefficient are chosen as follows: τ = 1 year, δC(x, t) = 0.01C cos(2πt) in x ∈ [0.9, 1.0]×10 6 m. We reuse the MISMIP (Pattyn et al., 2012) test case EXP 1 as in Sect. 2.2.2. The parameters of the setup are the same as in Fig. 2 and are given in Table 1. The variables u 1 and h are observed at x ∈ [0.85, 1.02] × 10 6 m. The steady-state solution of the forward equation with the GL located at x GL = 1.035 × 10 6 m is perturbed by δu 1 and δh when C is perturbed by δC as expressed in formulas δu 1 = u 1 (C+δC)−u 1 (C) and δh = h(C+δC)−h(C). After perturbation, the GL position will oscillate in time. The ice sheet is simulated by FS with Elmer/Ice (Gagliardini et al., 2013) for 10 years using the method implemented there for the position of the grounding line. Figure 3 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. 3.
The weight f (Tu)Tu · Tv 0 in (24) is negative, and an increase in the friction, δC > 0, leads to a decrease in the veloc-ity, and δC < 0 increases the velocity in all panels of Fig. 3. The velocities δu 1 and the surface elevations δh are separated by a phase shift in time, φ = π/2, as predicted by (24) and (25).
The weight in (25) 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. 3.

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 equation (19) are solved in both problems but with different driving functions defined by F (u, h) in (14).
Let (v i , q i , ψ i ), i = 1, . . ., d, be the steady-state solution to (19) when u i is observed atx and F u = e i δ(x − x). These are solutions to the sensitivity problem. We will show that the adjoint solution and the variation δF of the inverse problem can be expressed in (v i , q i , ψ i ). The perturbation δC is chosen such that δF < 0 in each step in the iterative solution of the inverse problem. Then the objective function F decreases stepwise toward the minimum.
It is shown in Appendix B that is a solution of (19) with arbitrary weights w i (x), i = 1, . . ., d, when When C is perturbed, the first variation of the functional in (22) is In the inverse problem in Petra et al. (2012), The weights in (27) for the inverse problem are w i (x) = u i (x) − u obs,i (x). Letṽ denote a weighted sum of the solutions of the sensitivity problem v i over the whole domain ω: Then the effect of δC on F in the inverse problem is by (28) The same construction of the solution is possible when h obs is given.
We have investigated the relation between the sensitivity problem and the inverse problem. By solving d sensitivity problems with F u = e i δ(x − x), i = 1, . . ., d, to obtain their adjoint solutions (v i , q i , ψ i ) and combine them with the weights w i from F u in (29) for the inverse problem, the adjoint solution to the inverse problem is (30). This solution can then be inserted into (28) to evaluate the effect in F of a change in C as in (31). In practice, if we are interested in solving the inverse problem and determine δF in (28) in order to iteratively compute the optimal solution with a gradient method, then we solve (19) directly with F u = u − u obs or F h = h − h obs to obtainṽ without computing d vectors v i . Taking δC = −f (Tu)Tu·Tṽ in (31) guarantees that δF < 0.

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 perturbations in the velocity for a two-dimensional vertical ice sheet at a 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 (19) is The velocity from the forward equation is u(x, z) = (u 1 , u 3 ) T , and the adjoint elevation ψ satisfies the right boundary condition ψ(L) = 0. The analytical solution ψ to (32) is derived in Appendix C. Let g(x) = u 1z (x) if u 1 is observed and let g(x) = 1 if h is observed. Then the adjoint solution is So, this solution has a jump −g(x * )/u 1 (x * ) at x * . With a small h · u z (y) ≈ 0 in (33), an approximate solution is ψ(x) ≈ −g(x * )/u 1 (x). If u 1 is observed and g(x) = u 1z ≈ 0, then ψ(x) ≈ 0 in (33) and ψh ≈ 0 in (19). 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 (19) is small. Solving only the adjoint stress equation for v as in Gillet-Chaulet et al. (2016), Isaac et al. (2015), and Petra et al. (2012) yields an adequate answer. Numerical solution in Cheng and Lötstedt (2020) of the adjoint FS equation (19) in two dimensions confirms that when u 1 is observed then ψ(x) is negligible. The situation is different when h is observed and ψ ≈ 1/u 1 (x * ) in (33).

Adjoint equations based on SSA
Starting from (9), a Lagrangian L of the SSA equations is defined, using the technique described and applied to the FS equations in Petra et al. (2012). The SSA Lagrangian in (A4) in Appendix A is similar to the FS Lagrangian in (15). By partial integration in L and evaluation at the forward solution (u, h), 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. 20 for the case of FS) From (34) it is seen that the adjoint SSA equations have the same structure as the adjoint FS equation (19). There is one stress equation for the adjoint velocity v and one equation for the Lagrange multiplier ψ corresponding to the surface elevation equation in (9). However, the advection equation for ψ in (34) depends on v, implying a fully coupled system for v and ψ. Equations (34) are solved backward in time with a final condition on ψ at t = T . As in (9), 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 Eq. 36 is actually the same as in Eq. 22) The same perturbations in δC, δb, and δC γ could be allowed for the FS equations in (22), but because the FS equations are 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 (34). Then the adjoint elevation equation must be solved for ψ = 0 to have a v = 0 in the adjoint stress equation and a perturbation in the Lagrangian in (36). The same result follows from the adjoint FS equations. If F h = 0 and F u = 0 in (19), then ψ = 0. Consequently, v = 0 and a perturbation δC will cause a perturbation δL in (22). The conclusion that the adjoint elevation equation must be solved if the surface elevation is observed is independent of the two ice models.

The two-dimensional adjoint solution
The adjoint SSA equations in two vertical dimensions are derived from (34) in the same manner as (11), by letting ψ and v 1 be independent of y and setting v 2 = 0 and C γ = 0. To simplify the notation, we also let v = v 1 . The adjoint equations for v and ψ follow either from simplifying (34) or from the Lagrangian and (11) 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 implicitly. The effect on the Lagrangian of perturbations δb and δC is obtained from (36): The weights or sensitivity functions w b and w C multiplying δb and δC in the integral are defined by The steady-state solutions to the system (37) can be analyzed as in the forward equations in Sect. 2.2.2 after simplifications. The viscosity terms in (37) are often small and can hence be neglected, and we assume that the basal topography is characterized by a small spatial gradient b x . The advantage resulting from these simplifications is that both the forward and adjoint equations can be solved analytically on a reduced computational domain where x ∈ [0, x GL ]. The analytical approximations are less accurate close to the ice divide where some of the above assumptions are not valid. The adjoint equation (37) reduce to

The two-dimensional adjoint steady-state solution with velocity observation
In this section, the analytical solution to the adjoint equation (40) is discussed. The derivation of the solution is detailed in Appendix E to Appendix F. It is here sufficient to recall that the solution given below is derived under the assumptions that b x H x , and that a and C are constants. For observations of u at x * , and the adjoint solutions are where ψ(x) and v(x) have discontinuities at the observation point x * . The perturbation of the Lagrangian (38) is with the Heaviside step function H(x) and the Dirac delta δ(x) (cf. Appendix F): or, after scaling with u * , The relation in (43) between the relative perturbations δb/H, δC/C, and δu/u can also be interpreted as a way to quantify the uncertainty in u. The uncertainty may be due to measurement errors in the topography b. For example, it is known that the true surface is in an interval [b − δb, b + δb] around b where, e.g., δb = 1 m or δb has a normal distribution with zero mean and some variance. Such an uncertainty δb in b or similarly an uncertainty δC in C is propagated to an uncertainty δu * in u at x * by (43); see Smith (2014). As an example, suppose that a 5 % error in the surface velocity is acceptable, then the tolerable error in the topography could be 20-30 m with a 1000 m thick ice.
The perturbations δu 1i at discrete points x * ,i , i = 1, 2, . . . due to perturbations δC j at discrete points x j , j = 1, 2, . . . are connected by a transfer matrix W C in Cheng and Lötstedt (2020). The relation between δu 1i and δC j is for all i and j such that The elements W Cij of the transfer matrix correspond to quadrature coefficients in the discretization of the first integral in (42) with δb = 0. The properties of W C are examined numerically in Cheng and Lötstedt (2020). We conclude that certain perturbations of C (not only highly oscillatory) are difficult to observe in u 1 at the surface. The same analysis is performed for the other combinations of δb, δC and δu 1 , δh.
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 a 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 (42) and (43) 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 and if the error 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 (45) that only δC with a long wavelength is visible at the surface and that δb also with a short wavelength affects δu * in (43). If δb is small or zero in (43), then it is easier to determine the δC that causes a certain δu * .
The analytical adjoint solutions ψ(x) and v(x) in (41) of the MISMIP case in Fig. 2 with parameters in Table 1 at different x * positions are shown in Figs. 4a and 5a.
The weights w b and w C in (42) multiplying δb and δC, defined in the same manner as in (38) and (39), are shown in Figs. 6a and 7a with the solutions ψ and v in Figs. 4a and 5a. The Dirac term is plotted as a vertical line at x * in Fig. 7a.
All perturbations in C between x * and x GL will result in a perturbation of the opposite sign in u * at the surface because w C < 0 in (x * , x GL ) in Fig. 6a and (42). The same conclusion holds true for perturbations in b because w b < 0 in (x * , x GL ) in Fig. 7a, 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 following conclusions can be drawn from (42) and (43) 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. 2). Or, compactly expressed, δC with support in [x * , x GL ] will cause larger perturbations at the surface the closer x * is to x GL and the closer δC(x) is to x GL . The same conclusion is drawn in Cheng and Lötstedt (2020) with numerically computed SSA adjoint solutions.
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 , as well as 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 Eq. 42).
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 follows if we let the perturbation of the friction coefficient be a constant δC 0 = 0 in [x 0 , x 1 ] ∈ [x * , x GL ] and evaluate the integral in (42) to obtain The same δu * is observed with a constant perturbation .
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 (42) 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 of the perturbation 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 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 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. How the damping de-pends on λ in the FS equations is computed in Cheng and Lötstedt (2020).
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 (42). The effect is larger if the ice is thin and moving fast. The integral term will behave in the same way as in (45), with mainly perturbations with small wave numbers and long wavelengths visible at the surface. Figure 6. The analytical solution of the weights w C = −vu m on δC in (38) for (a) u and (b) h observed at x * = 0.25 × 10 6 , 0.5 × 10 6 , 0.7 × 10 6 , and 0.9 × 10 6 m. Figure 7. The analytical solution of weights w b = ψ x u+v x ηu x +vρgh x on δb in (38) for (a) u and (b) h observed at x * = 0.25×10 6 , 0.5× 10 6 , 0.7 × 10 6 , and 0.9 × 10 6 m.

The two-dimensional adjoint steady-state solution with elevation observation
In the case when h is observed at x * and F u = 0 and F h = δ(x − x * ), the expressions for ψ and v satisfying (40) are The corresponding formulas when u is observed are found in (41). There is a discontinuity at the observation point x * in v(x) illustrated in Fig. 5b, but ψ(x) is continuous in the solution of (40) and in Fig. 4b. The second derivative term ( 1 n ηH v x ) x is neglected in the simplified Eq. (40) but is of importance at x * . A correctionψ of ψ at x * in (46) is therefore introduced to satisfy 1 n ηH v x x −Hψ x = 0. With v x (x * ) = −δ(x −x * )/(ρgH * ), the correction isψ(x) = −δ(x − x * )η * /(nρgH * ). The solution ψ is updated at each x * in Fig. 4b, withψ as a vertical line representing the negative Dirac delta.
The perturbation in h is as in (42) with ψ and v in (46) and the additional termψ: where a(xδb/H ) x (x * ) = (uδb) x (x * ) represents the xderivative of uδb evaluated at x * . When δb = 0 then δu * in (43) and δh * = δH * in (47) satisfy δu * H * = −δH * u * as in the integrated form of the advection equation in (12) and in (D1) in the Appendix. As in (42), (47) is rewritten with the weights w b and w C in (39): These weights are shown in Figs. 6b and 7b. The negative derivative of the Dirac delta is depicted in Fig. 7b as a vertical line in the negative direction immediately followed by one in the positive direction. The contribution from the integrals in (43) and (47) is identical except for the sign (compare w C in Fig. 6a and b and w b in Fig. 7a and b). The first term in (43) depends on δb/H , and the first term in (47) depends on the derivative of axδb/H = uδb. The derivative of uδb at x * directly affects the perturbation of h at x * . A perturbation of b at the base is directly visible locally in u at the surface, while the effect of δC is non-local in the integral in (47). Because of the similarities between (43) and (47) (42) and (43) for δu * are valid also for δh * in (47).

The two-dimensional time-dependent adjoint solution
Finally, the time-dependent adjoint equation (37) is investigated. Equation (37) is solved numerically for the same MIS-MIP test case as Fig. 2 in Sect. 2.2.2 with the parameters in Table 1. As in Sect. 3.1.2, the friction coefficient C has a seasonal variation (period τ 1 year, where the beginning of the year is associated with winter) in the forward equation (11): 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. 3. The amplitude of the variation in C is set to κ = 0.5, and the forward equation (11) is solved for 11 years. The GL is determined by the conditions in Sect. 2.2 and will move in time because of the variation in C. The topography b is kept constant in time. Observations of u and h are made at x * = 9 × 10 5 m for 0.1 year in the four seasons starting from the summer of the 10th 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 equation (11) is solved numerically from t = 0 with the steady-state solution as initial data to the observation points t = t * , and the adjoint equation (37) is solved from t = t * backward in time to t = 0.
According to a convergence test, the time step is chosen to be 0.01 year 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 numerical solution has been reached. Figure 8 shows the results for the adjoint weights w C (x, t) and w b (x, t) multiplying the perturbations δC and δb, as defined in (38), for the observations of u and h at x * = 9×10 5 m in all four seasons, where each column represents one season. The friction coefficient C follows the seasonal variation in (49). Each row is one of the combinations of the weights w C and w b for the observations of u and h. The time axis (or ordinate) in the figure follows the time direction in the forward problem (11). 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 seasonal variation in the basal conditions in (49). 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 5 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. 9a 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 equation (24) 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. 9a; 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 (49) (the blue lines in Fig. 9a 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. 9c and d. An intuitive explanation is that there is an immediate effect on the velocity, but there is a delay in h since it 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. 9c. For the same δC, the largest δh * is observed in the fall (orange), and 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 between the perturbation and the observation of the effect in the surface elevation. The same shift in time is what we found in Sect. 3.1.1, (25), and Fig. 3 for the FS equations.  A reference adjoint solution at x * observed during the fall season (t * = 9.75) with time-independent C and b, κ = 0 in (49), as parameters in the forward equations is shown in black dashed lines in all the four panels of Fig. 9. The weight w b at x * for a constant b is well approximated by w b * exp(−(T −t)/τ ) in time with τ = 1.4 years for some w b * for the observation of both u and h in Fig. 9b and d. For the weight w C , the same exponential function holds with weight w C * , but the time constant τ = 1.8 years for the observation of h * in Fig. 9c and τ = 2.2 years for the u * case.
Suppose that the temporal perturbation is oscillatory with frequency f and located in space at x * with δC(x, t) = δC 0 cos(2πf t)δ(x − x * ).
A low-frequency f with f 1 corresponds to decennial or centennial variations and a high-frequency f with f 1 corresponds to diurnal or weekly variations. Then the perturbation in h at t = t * is cf. (45). If the frequency is high, f 1, then δh * ∝ 1/f and high-frequency perturbations are damped efficiently. At certain times of observation t * when sin(2πf t * ) = 0, the damping is even stronger with δh * ∝ 1/f 2 . 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 adjoint equations are derived in the FS and the SSA frameworks including time and the surface elevation equation. Time-dependent perturbations δC and δb in basal friction coefficient C and basal topography b are introduced, and their effect on observations of the velocity u and the surface elevation h at the top surface of the grounded ice is studied. With the solution of the adjoint equations, we can determine the perturbation at a given point in space and time on the surface due to all basal perturbations. By solving the forward equations twice with C and C + δC or b and b + δb, we can compute the perturbation in all points in space and time on the surface, e.g., δu 1 = u 1 (C + δC) − u 1 (C) in the first velocity component, for a given δC or δb.
The perturbations in the observations are determined numerically in Cheng and Lötstedt (2020) either using the adjoint equations and their solutions in Sect. 3 or by solving the forward equations with unperturbed and perturbed parameters to obtain δu * and δh * . The numerical solutions are compared to each other and to analytical solutions for SSA. The agreement is good in the comparisons.
In Sect. 3.1.3, a relation is established between the inverse problem (aiming to infer parameters from data) and the sensitivity problem (aiming to quantify the effect of variations in parameters): the same adjoint equations are solved. However, the forcing functions differ and are specific to the inverse problem and the sensitivity problem, respectively. Common to both problems is that the adjoint equations tell how perturbations in the parameters at the ice base are propagated to perturbations in the velocity and the elevation of the surface.
For steady-state problems, and in an FS setting where u is observed, we find (cf. Sect. 3.1.4) that the contribution of the solution of the adjoint elevation equation (32) is small and that it therefore suffices to solve only the adjoint stress equations -see, e.g., Gillet-Chaulet et al. (2016), Isaac et al. (2015), and Petra et al. (2012) -in order to be able to draw conclusions regarding perturbations of u. For steady-state problems in a two-dimensional SSA setting for a vertical ice, (43), (47), and Figs. 6 and 7 show that the sensitivity of the velocity and elevation increases (because the velocity increases and the ice thickness decreases) as the observation point x * approaches the grounding line.
In this setting, a non-local effect of a perturbation in C has been observed, in the sense that δC(x) affects both u(x * ) and h(x * ) even if x = x * , but a perturbation δb in b has a strong local effect concentrated at x * . Nevertheless, the shapes of the two sensitivity functions (or weights) for δb and δC are very similar except for the neighborhood of x * , which makes it difficult to separate their respective contribution in an observation. Different combinations of the perturbations in the basal friction and bedrock elevation can produce the same effect on the velocity and surface elevation changes at one observation point.
In the inverse problems based on time-dependent simulations of FS and SSA, it is necessary to include the adjoint elevation equation. If the perturbations in the basal conditions are time dependent and h is observed (see Figs. 3, 9c, and d), then time cannot be ignored in the inversion. If time dependence is ignored, wrong conclusions concerning the conditions at the ice base may be drawn from observations of h, in both the FS and the SSA model. In the time-dependent solution of SSA, a 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.
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 (45) and (50), 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.  (20) 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: Thus, let η ij kl (u) = η(u)(I ij kl + B ij kl (u)), or in tensor form Replacing B in (A2) by D, we obtain the adjoint FS viscosity in (20). The adjoint friction in SSA in ω and at γ g in (34) with a Weertman law is derived as in the adjoint FS equations (19) and (20). Then in ω with ξ = u, ζ = v, c = C, and F = F ω and at γ g with ξ = t · u, ζ = t · v, c = C γ , f = f γ , and F = F γ , we arrive at the adjoint friction term cf (ξ ) (I + F(ξ )) ζ , where

A2 Adjoint equations in SSA
The Lagrangian for the SSA equations is with the adjoint variables ψ, v, q after partial integration and using the boundary conditions. The perturbed SSA Lagrangian is split into the unperturbed Lagrangian and three integrals: L(u + δu, h + δh; v + δv, ψ + δψ; b + δb, C γ + δC γ , The perturbation in L is Terms of order 2 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 (9), the adjoint viscosity (35), (A2), the friction coefficient (A3), Gauss' formula, and the boundary conditions, as well as 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 (9) and (34) are inserted into (A4) resulting in L(u * , p * ; v * , q * ; h * , ψ * ; b, C γ , C) Then (A12) yields the variation in L in (A13) with respect to perturbations δb, δC γ , and δC in b, C γ , and C:
Terms of order 2 or more in δu, δv, and δ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: T 0 s δψ(h t + h · u − a s ) + ψ(δh t + u · δh + u z · hδh + h · δu) dx dt = + h · u z ψ)δh + h · δuψ dx dt.
The first variation of L is then δL =I 1 + I 2 + I 3 = T 0 s (F u + ψh) · δu dx dt + T 0 s (F h + F u u z − (ψ t + ∇ · (uψ) − h · uψ))δh dx dt With the forward solution (u * , p * , h * ) and the adjoint solution (v * , q * , ψ * ) satisfying (4) and (19), the first variation with respect to perturbations δC in C is (cf. A14) δL = or observation of h with d = 1 Introduce the weight functions w i (x), i = 1, . . ., d. It follows from (19) that (w i (x)v i (x), w i (x)q i (x), w i (x)ψ i (x)) is a solution with F u = w i (x)e i δ(x − x) or F h = w(x)δ(x − x). Therefore, also is a solution with F u = ω w i (x)e i δ(x − x) dx = w i (x)e i or F h = ω w(x)δ(x−x) dx = w(x). A sum over i, i = 1, . . ., d, of each integral in (B3) is also a solution.
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 (22), δ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 (19) 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), 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.
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) v To satisfy the jump condition in (E8) and (E9), the constant C v is Combine (F1) with the relation ψ x = (F u −Cmu m−1 v)/H and integrate from x to x GL to obtain ψ(x) = C v a m−1 C x m GL − x m , x * < x ≤ x GL .
With the jump condition in (E8) and (E9), ψ(x) at 0 ≤ x < x * is The weight for δC in the functional δL in (38) is non-zero for x * < x ≤ x GL : Use (F1) and (40) in (38) to determine the weight for δb in δL: