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

https://doi.org/10.5194/tc-15-715-2021

© Author(s) 2021. This work is distributed under the Creative Commons Attribution 4.0 License.

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

Gong Cheng ¹ , Nina Kirchner ^2,3 , and Per Lötstedt ¹

1 Department of Information Technology, Uppsala University, P.O. Box 337, 751 05 Uppsala, Sweden

2 Department of Physical Geography, Stockholm University, 106 91 Stockholm, Sweden

3 Bolin Centre for Climate Research, Stockholm University, 106 91 Stockholm, Sweden Correspondence: Gong Cheng (cheng.gong@uci.edu)

Received: 16 April 2020 – Discussion started: 13 May 2020

Revised: 15 December 2020 – Accepted: 15 December 2020 – Published: 15 February 2021

Abstract. Predictions of future mass loss from ice sheets are afflicted with uncertainty, caused, among others, by insuffi- cient understanding of spatiotemporally variable processes at the inaccessible base of ice sheets for which few direct obser- vations exist and of which basal friction is a prime example.

Here, we present a general numerical framework for study- ing 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 equa- tions, derived in the framework of a full Stokes model and a shallow-shelf/shelfy-stream approximation model, respec- tively, to determine the sensitivity of grounded ice sheet sur- face velocities and elevation to time-dependent perturbations in basal friction and basal topography. Analytical and numer- ical examples are presented showing the importance of in- cluding 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 verti- cal 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 pertur- bation in the friction is propagated directly to the surface in time.

1 Introduction

Over the last decades, ice sheets and glaciers have expe- rienced mass loss due to global warming, both in the po- lar regions and also outside of Greenland and Antarctica (Farinotti et al., 2015; Mouginot et al., 2019; Pörtner et al., 2019; Rignot et al., 2019). The most common benchmark date for which future ice sheet and glacier mass loss and as- sociated 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 be- yond 2100 CE in the warming scenarios typically referred to as RCPs (representative concentration pathways; see van Vu- uren et al., 2011), but rates and ranges are afflicted with un- certainty, caused by, among others, insufficient understand- ing of spatiotemporally variable processes at the inaccessi- ble base of ice sheets and glaciers (Pörtner et al., 2019; Ritz et al., 2015). These include the geothermal heat regime, sub- glacial and base-proximal englacial hydrology, and particu- larly 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 differ-

ent 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 slid- ing laws for glacier and ice sheet flow emerged in the 1950s – see, e.g., Fowler (2011), Iken (1981), Lliboutry (1968), Nye (1969), Schoof (2005), and Weertman (1957) – have subse- quently been implemented into numerical models of ice sheet and glacier behavior – e.g., in Brondex et al. (2017), Bron- dex 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 con- troversially (Minchew et al., 2019; Stearn and van der Veen, 2018).

Because few or no observational data are available to con- strain 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 computation- ally by solving an inverse problem provided that observa- tions 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 de- scribed for the case of the sliding law, can be used to deter- mine 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 circula- tion (Heimbach and Losch, 2012), the geothermal heat flux at the ice base (Zhu et al., 2016), or to estimate basal topog- raphy beneath an ice sheet (Monnier and des Boscs, 2017;

van Pelt et al., 2013). The latter is not only difficult to sepa- rate 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 (2003, 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 coun- terparts 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 ob- jective 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 thick- ness 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 (Gold- berg 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 advec- tion 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 parame- ters 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 ele- vation at the surface. The parameter sensitivity of the solu- tions 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 ve- locity 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 parame- ters 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 equa- tions 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 observa- tions 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) be- cause it includes the time-dependent advection equation for the kinematic free surface. The key concepts and steps in- troduced in Sect. 3 are supplemented by detailed derivations, collected in Appendix A. The same adjoint equations are ap- plicable 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 parame- ters in the friction model. Analytical solutions in two dimen- sions of the stationary adjoint equations subjected to simpli- fying 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 per- turbations 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 perturba- tions with the forward equations. Discrete transfer functions are computed and analyzed as in Gudmundsson (2008) for the relation between surface perturbations and basal pertur- bations. 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öt- stedt (2020) confirm the conclusion in Gudmundsson (2008) that the perturbations with a long wavelength and low fre- quency will propagate to the surface while those of a short wavelength and high frequency are damped.

2 Ice models

In this section, the equations emerging from adopting a for- ward approach to describing ice dynamics are presented, to- gether with relevant boundary conditions, for the FS (4) and SSA model (9), respectively. These, and the notation and ter- minology 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 re- ferred 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 exten- sions 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 ra- tio V /L is on the order 10 ⁻² to 10 ⁻³ . The aspect ratio is used to introduce simplifications to the FS equations, result- ing, e.g., in the shallow-ice (SIA) (Hutter, 1983), shallow- shelf (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 advec- tion 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 in- terface between floating ice shelves and sea water in the ice shelf cavities is usually regarded as frictionless.

2.1 Full Stokes model

We adopt standard notation and denote vectors in bold italics and matrices in three-dimensional space in bold, and we de- note 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 ₃ the velocity components of u = (u ₁ , u 2 , u 3 ) ^T in the x, y, and z direction, where x = (x, y, z) ^T is the posi- tion 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 down- stream of γ _GL we have z _b > b (see Fig. 1). In two dimen- sions, γ _GL consists of one point with x coordinate x _GL .

The boundary 0 enclosing the domain  occupied by the ice has different parts (see Fig. 1b). The lower boundaries of

 are denoted by 0 b (where the ice is grounded at bedrock) and 0 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 up- per boundary is denoted by 0 _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 0 _u in black and a down- stream part denoted by 0 _d in blue, where 0 = 0 _u ∪ 0 _d . Obvi- ously, if x ∈ 0 _u , then (x, y) ∈ γ _u , or if x ∈ 0 _d , then (x, y) ∈ γ _d , where γ = γ _u ∪ γ _d . Letting n be the outward-pointing nor- mal on 0 (or γ in two dimensions (x, y)), the nature of ice flow renders the conditions n · u ≤ 0 at 0 u and n · u > 0 at 0 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

 = {x|(x, y) ∈ ω, z _b (x, y, t ) ≤ z ≤ h(x, y, t )}, 0 s = {x|(x, y) ∈ ω, z = h(x, y, t)},

0 b = {x|(x, y) ∈ ω, z = z _b (x, y, t )

= b(x, y), x < x _GL (y)},

0 w = {x|(x, y) ∈ ω, z = z _b (x, y, t ), x > x GL (y)}, 0 _u = { x|(x, y) ∈ γ _u , z _b (x, y, t ) ≤ z ≤ h(x, y, t )},

0 d = {x|(x, y) ∈ γ _d , z b (x, y, t ) ≤ z ≤ h(x, y, t )}. (1)

Before the forward FS equations for the evolution of the

ice surface 0 _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 0 _s by a _s , and grav-

itational acceleration by g. The values of these physical pa-

rameters are given in Table 1. On 0 _s , h = (h _x , h _y , −1) ^T de-

scribes the spatial gradient of the ice surface (in two vertical

Figure 1. A schematic view of an ice sheet in the (a) x–z plane and (b) x–y plane.

dimensions h = (h _x , −1) ^T ). The strain rate D and the viscos- ity η are given by

D(u) = ¹ ₂ (∇u + ∇u ^T ), η(u) = ¹ ₂ A ⁻

¹ⁿ

(trD ² (u))

¹⁻ⁿ²ⁿ

, (2) where trD ² is the trace of D ² . 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

σ (u, p) = 2ηD(u) − Ip, (3)

where p is the isotropic pressure and I is the identity matrix.

Turning to the ice base, the basal stress on 0 _b is re- lated to the basal velocity using an empirical friction law Tσ n = −βTu, where a projection (Petra et al., 2012) on the tangential plane of 0 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

(a ⊗ c) _ij = a _i c _j , (A ⊗ C) _{ij kl} = A _ij C _kl .

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.

For instance, f (u) = kuk ^m−1 with the norm kuk = (u · u)

¹²

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 ¹ ₃ and 1.

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

h _t + h · u = a _s , on 0 _s , t ≥ 0,

h(x, 0) = h 0 (x), x ∈ ω, h(x, t) = h _γ (x, t), x ∈ γ _u ,

− ∇ · σ (u, p) = −∇ · (2η(u)D(u)) + ∇p = ρg,

∇ · u = 0, in (t), σ n = 0, on 0 s ,

Tσ n = −Cf (Tu)Tu, n · u = 0, on 0 b , (4) where h ₀ (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 0 _w :

z bt + uz bx + vz by − w = a b , on 0 w , t ≥ 0,

σ n = −p _w n, on 0 w , (5)

where a _b is the basal mass balance and p _w is the water pres- sure at 0 _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

d(x, t) = z _b (x, t) − b(x).

On 0 _w , we have d(x, t ) > 0 and n · σ n + p _w = 0 and on 0 _b , we have d(x, t ) = 0 and n · σ n + p _w < 0. The complemen- tarity relation at the ice base 0 _b ∪ 0 _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 exam- ple, 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 0 u and 0 d are of Dirichlet type such that

u| ₀

_u

= u _u , u| ₀

_d

= u _d , (6)

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 0 _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 0 _u .

2.2 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 consid- ered. 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 of- ten also used in regions of fast flow over lubricated bedrock (MacAyeal, 1989).

The simplifications associated with adopting the SSA im- ply that the viscosity (see 2 for the FS model) is now given by

η(u) = 1 2 A ⁻

¹ⁿ



u ² _1x + u ² _2y + 1

4 (u _1y + u _2x ) ² + u _1x u _2y



¹⁻ⁿ_2n

= 1 2 A ⁻

¹ⁿ

 1

2 B : D



¹⁻ⁿ_2n

, (7)

where B(u) = D(u) + ∇ · u I, with ∇ · u = tr D(u). This η differs from (2) because B 6= D due to the cryostatic ap- proximation of p in the SSA. In (7), the Frobenius inner product between two matrices A and C is used, defined by A : C = P

ij A _ij C _ij . The vertically integrated stress tensor ς (u) (cf. 3 for the FS model) is given by

ς (u) = 2H ηB(u) , (8)

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 horizon- tal velocity are neglected in SSA. Then, Weertman’s law is β(u, x, t ) = C(x, t )f (u) = C(x, t )kuk ^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 sur- face and ice velocities based on an SSA model (in which u is not divergence-free) read

H t + ∇ · (uH ) = a, t ≥ 0, x ∈ ω,

h(x, 0) = h 0 (x), x ∈ ω, h(x, t) = h _γ (x, t), x ∈ γ _u ,

∇ · ς − Cf (u)u = ρgH ∇h, x ∈ ω,

n · u(x, t) = u _u (x, t), x ∈ γ _u , n · u(x, t) = u _d (x, t), x ∈ γ _d ,

t · ς n = −C _γ f _γ (t · u)t · u, x ∈ γ _g ,

t · ς n = 0, x ∈ γ w . (9)

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 ve- locities u _u ≤ 0 and u _d > 0 are specified on γ _u and γ _d . The lat- eral side of the ice γ is split into γ _g and γ _w with γ = γ _g ∪ γ _w . There is friction in the tangential direction on γ _g which de- pends on the tangential velocity t · u with the friction coef- ficient 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 ground- ing line (see, e.g., Seroussi et al. (2014)). The thickness H is compared to a flotation height H _f given by

H f = − ρ w z b /ρ.

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 as- sumed to be at hydrostatic equilibrium. In such a case, a calv- ing 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):

ς (u) · n = 1 2 ρgH ²

 1 − ρ

ρ _w



n, x ∈ γ _d . (10)

2.2.1 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 x–

z plane (see Fig. 1). The forward SSA equations are de- rived from (9) by letting H and u ₁ be independent of y and setting u ₂ = 0. Since there is no lateral force, C _γ = 0.

The position of the grounding line is denoted by x _GL , and 0 _b = [ 0, x _GL ] , 0 _w = (x _GL , L]. 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

H _t + (uH ) _x = a, 0 ≤ t ≤ T , 0 ≤ x ≤ L, h(x, 0) = h 0 (x), h(0, t ) = h L (t ), (H ηu _x ) _x − Cf (u)u − ρgH h _x = 0,

u(0, t ) = u _u (t ), u(L, t ) = u d (t ), (11) 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

u _x (L, t ) = A  ρgH (L, t)

4 (1 − ρ

ρ _w )

 n

.

Assuming that u > 0 and u _x > 0, the viscosity becomes η = 2A ⁻

¹ⁿ

u

1−n

x

n

, and the friction term with a Weertman law turns into Cf (u)u = Cu ^m .

2.2.2 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 longi-

tudinal stress can be ignored in the steady-state solution (see

Table 1. The model parameters.

Parameter Quantity

a = 0.3 m yr ⁻¹ Surface mass balance A = 1.38 × 10 ⁻²⁴ s ⁻¹ Pa ⁻³ Rate factor of Glen’s flow law C = 7.624 × 10 ⁶ Pa m ^−1/3 s ^1/3 Basal friction coefficient g = 9.8 m s ⁻² Acceleration of gravity

m = 1/3 Friction law exponent

n = 3 Flow-law exponent

ρ = 900 kg m ⁻³ 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

(uH ) x = a, 0 ≤ x ≤ x GL , H (x GL ) = H GL ,

− Cu ^m − ρgH h x = 0, u(0) = 0.

(12)

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):

H (x) =



H _GL ^m+2 + m + 2 m + 1

Ca ^m

ρg (x _GL ^m+1 − x ^m+1 )



_m+2¹

, 0 ≤ x ≤ x _GL ,

H (x) = H GL , x GL < x < L, u(x) = ax

H , 0 ≤ x ≤ x _GL , u(x) = ax H _GL ,

x GL < x < L. (13)

The solution is normalized with the ice thickness H _GL = H (x GL ) at the grounding line. Both x GL and H GL are as- sumed 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 involv- ing H (0) is the solution of the equations. By including the longitudinal stress in the ice, an approximate, more compli- cated 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 ⁶ m on a single slope bed defined by b(x) = 720 −

778.5x

7.5×10

⁵

and lifts from it at the grounding line position x _GL = 1.035×10 ⁶ 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. Nu- merical 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).

3 Adjoint equations

In this section, the adjoint equations are discussed, as emerg- ing 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 coeffi- cient can be perturbed with δC 6= 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 deriva- tions 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 0 _s and over the time interval [0, T ], we consider the functional F :

F =

T

Z

0

Z

0

s

F (u, h) dx dt . (14)

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 func- tion in x. We distinguish two cases: u and h satisfy either the FS equation (4) or the SSA equation (9). Given F , the for- ward 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 La- grangian for the FS equations is

L(u, p, h; v, q, ψ; C) =

T

Z

0

Z

0

s

F (u, h) + ψ(h _t + h · u − a) dx dt

+

T

Z

0

Z

ω h

Z

b

− v · (∇ · σ (u, p)) − ρg · v

− q∇ · u dx dt, (15)

with the Lagrange multipliers v, q, and ψ corresponding to the forward equations for u, p, and h. The effect of perturba- tions δC and δb in C and b on F is given by the perturbation δL, viz.

δF = δL = L(u + δu, p + δp, h + δh; v + δv, q + δq,

ψ + δψ ; b + δb, C + δC) − L(u, p, h; v, q, ψ; b, C). (16)

Figure 2. The analytical solutions u(x) and H (x) in (13) for a grounded ice in [0, x _GL ] .

Examples of F (u, h) in (14) are F = ku − u _obs k ² , and F = |h − h _obs | ² 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 0 s (see also Gillet-Chaulet et al., 2016; Isaac et al., 2015; Morlighem et al., 2013; and Petra et al., 2012).

Another example is F (u, h) = _T ¹ u 1 (x, t)δ(x − x ∗ ) with the Dirac delta δ at x ∗ to measure the time-averaged horizontal velocity u ₁ at x ∗ on the ice surface 0 _s with

F =

T

Z

0

Z

0

s

F (u, h) dx dt = 1 T

T

Z

0

u ₁ (x ∗ , t )dt, (17)

where T is the duration of the observation at 0 _s . If the horizontal velocity is observed at (x ∗ , t ∗ ) then F (u, h) = u ₁ (x, t)δ(x − x ∗ )δ(t − t ∗ ) and

F =

T

Z

0

Z

0

_s

F (u, h) dx dt = u 1 (x ∗ , t ∗ ). (18)

The sensitivity in F and u ₁ 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 for- ward and adjoint equations are then solved once. In the in- verse problem with velocity observations, F is the objective function in a minimization procedure and F = ku − u obs k ² . 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 up- dated 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.

3.1 Adjoint equations based on the FS model

In this section, we introduce the adjoint equations and the perturbation of the Lagrangian function. The detailed deriva- tions 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:

ψ _t + ∇ · (uψ) − h · u _z ψ = F _h + F _u · u _z , on 0 _s , 0 ≤ t ≤ T ,

ψ (x, T ) = 0, ψ(x, t) = 0, on 0 d ,

− ∇ · ˜ σ (v, q) = −∇ · (2 ˜η(u) ? D(v)) + ∇q = 0,

∇ · v = 0, in (t),

σ (v, q)n = −(F ˜ _u + ψh), on 0 s ,

T ˜ σ (v, q)n = −Cf (Tu) (I + F b (Tu)) Tv, on 0 b ,

n · v = 0, on 0 b , (19)

where the derivatives of F with respect to u and h are F _u =  ∂F

∂u ₁ , ∂F

∂u ₂ , ∂F

∂u ₃

 T

, F _h = ∂F

∂h .

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 steady-

state adjoint equation for the FS equations discussed in Pe-

tra et al. (2012), an advection equation is added in (19) for

the Lagrange multiplier ψ (x, t ) on 0 _s with a right-hand side

depending on the observation function F and one term de-

pending on ψ in the boundary condition on 0 _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 6= 0 and F u = 0, then the adjoint elevation

equation must be solved together with the adjoint stress equa-

tion. 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 are η(u) ˜ = η(u) 

I + _nD(u):D(u) ¹⁻ⁿ D(u) ⊗ D(u)  ,

σ (v, q) ˜ = 2 ˜η(u) ? D(v) − qI; (20) 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 (A ? C) _ij = P

kl A _{ij kl} C _kl . In general, F _b (Tu) in (19) is a linearization of the friction law f (Tu) in (4) with respect to the vari- able Tu. For instance, with a Weertman-type friction law, f (Tu) = kTuk ^m−1 ,

F _b (Tu) = m − 1

Tu · Tu (Tu) ⊗ (Tu). (21)

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 veloci- ties, Tu and Tv at the ice base 0 b , and is given by

δF = δL =

T

Z

0

Z

0

_b

f (Tu)Tu · Tv δC dx dt. (22)

For this formula to be accurate, δC has to be small. Other- wise, nonlinear effects may be of importance.

3.1.1 Time-dependent perturbations

Let us now investigate the effect of time-dependent perturba- tions in the friction parameter on modeled ice velocities and ice surface elevation. Suppose that the velocity component u _1∗ = u ₁ (x ∗ , t ∗ ) is observed at (x ∗ , t ∗ ) at the ice surface as in (18) and that t ∗ < T , then

u _1∗ = F =

T

Z

0

Z

0

s

F (u) dx dt, (23)

with F (u) = u ₁ δ(x − x _∗ )δ(t − t ∗ ), F u

₁

= δ(x − x _∗ )δ(t − t ∗ ), F u

₂

= F u

₃

= 0, and F h = 0. 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 ⁱ denote the unit vector with 1 in the ith compo- nent, the boundary condition in (19) becomes ˜ σ (v, q)n =

− e ¹ δ(x − x ∗ )δ(t − t ∗ ) − ψh at t = t ∗ . For t < t ∗ , ˜ σ (v, q)n =

− ψh. Since ψ is small for t < t ∗ (see Sect. 3.1.4), the domi- nant part of the solution is v(x, t ) = v ₀ (x)δ(t − t ∗ ) for some v ₀ .

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 0 _s is, according to (22), given by

δu _1∗ = δL =

T

Z

0

Z

0

b

f (Tu)Tu · Tv δC dx dt

≈ Z

0

_b

f (Tu)Tu · Tv ₀ δC(x, t ∗ ) dx. (24)

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 fric- tion.

Further, to investigate the response of the ice surface eleva- tion, h ∗ at 0 _s , to perturbations in basal friction, one considers F (h) = h(x, t )δ(x − x _∗ )δ(t − t ∗ ), F h = δ(x − x _∗ )δ(t − t ∗ ), F u = 0.

The solution of the adjoint equation (19) with ˜ σ (v, q)n =

− ψh at 0 s for v(x, t ) is non-zero since ψ (x, t ) 6= 0 for t <

t ∗ .

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 ₀ (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 (1 year, from January to December), then Northern Hemi- sphere 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 con- stant 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

δh ∗ = δL = R T 0

R

0

_b

f (Tu)Tu · Tv δC(x, t ) dx dt

≈ J R t

∗

0 cos(2π t /τ ) dt = J _2π ^τ sin(2π t ∗ /τ ),

(25)

where the integral J is J =

Z

0

_b

f (Tu)Tu · Tv δC 0 dx. (26)

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 sum- mer), while δh ∗ 6= 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 ₁ in (24) where a pertur- bation 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 mea- sured in the summer with δh(x, t ∗ ) = 0, then the wrong con- clusion 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 ₀ , t 1 ] and 0 otherwise. Then with J in (26), δh ∗ in (25) is

δh ∗ ≈ J

t

∗

Z

0

s(t ) dt =







0, t ∗ ≤ t ₀ , J (t ∗ − t ₀ ), t ₀ < t ∗ < t ₁ , J (t ₁ − t ₀ ), t ∗ > t ₁ .

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.

3.1.2 Example with seasonal variation

To illustrate the phase delay in an oscillatory perturba- tion, a two-dimensional numerical example is shown in Fig. 3, where the timescale and friction coefficient are cho- sen as follows: τ = 1 year, δC(x, t ) = 0.01C cos(2π t ) in x ∈ [0.9, 1.0]×10 ⁶ 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 ₁ and h are observed at x ∈ [0.85, 1.02] × 10 ⁶ m.

The steady-state solution of the forward equation with the GL located at x _GL = 1.035 × 10 ⁶ m is perturbed by δu ₁ and δh when C is perturbed by δC as expressed in formulas δu ₁ = u ₁ (C +δC)−u ₁ (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 ⁶ 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 ⁶ . 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 oscil- latory perturbation, as in Fig. 3.

The weight f (Tu)Tu · Tv ₀ in (24) is negative, and an in- crease 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 ₁ and the surface elevations δh are separated by a phase shift in time, 1φ = π/2, as predicted by (24) and (25).

The weight in (25) for δC ₀ in the integral over x changes sign when the observation point is passing from x ∗ = 0.9 × 10 ⁶ to 1.0 × 10 ⁶ , explaining why the shift changes sign in the red dashed lines shown in the two lower panels of Fig. 3.

3.1.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 esti- mated 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 equa- tion (19) are solved in both problems but with different driv- ing functions defined by F (u, h) in (14).

Let (v ⁱ , q ⁱ , ψ ⁱ ), i = 1, . . ., d, be the steady-state solution to (19) when u _i is observed at ¯ x and F u = e ⁱ δ(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 ⁱ , q ⁱ , ψ ⁱ ). The perturbation δC is cho- sen 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



 Z

ω d

X

i=1

w _i ( ¯ x)v ⁱ dx, Z

ω d

X

i=1

w _i (x)q ⁱ dx, Z

ω d

X

i=1

w _i (x)ψ ⁱ dx



 is a solution of (19) with arbitrary weights w _i (x), i = 1, . . ., d, when

F u = Z

ω d

X

i=1

w _i (x)e ⁱ δ(x − x) dx =

d

X

i=1

w _i (x)e ⁱ . (27)

When C is perturbed, the first variation of the functional in (22) is

δF = Z

0

b

f (Tu)Tu · T



 Z

ω d

X

i=1

w _i (x)v ⁱ dx



 δC dx. (28)

In the inverse problem in Petra et al. (2012), F = 1

2 Z

ω

k u(x) − u _obs (x)k ² dx, F u = u(x) − u _obs (x). (29)

The weights in (27) for the inverse problem are w _i (x) = u _i (x) − u _obs,i (x). Let ˜v denote a weighted sum of the so- lutions of the sensitivity problem v ⁱ over the whole domain ω:

v(x) = ˜ Z

ω d

X

i=1

(u _i (x) − u _{obs, i} (x))v ⁱ dx. (30)

Figure 3. Observations at x ∗ = 0.85, 0.9, 0.95, 1.0, 1.02 × 10 ⁶ m with FS in time t ∈ [0, 10] of δu ₁ (solid blue) and δh (dashed red) with perturbation δC(t ) = 0.01C cos(2π t ) for x ∈ [0.9, 1.0] × 10 ⁶ m. Notice the different scales on the y axes.

Then the effect of δC on F in the inverse problem is by (28) δF =

Z

0

_b

f (Tu)Tu · T ˜v(x) δC dx. (31)

The same construction of the solution is possible when h _obs is given. Then d = 1, F (h) = ¹ ₂ (h − h _obs ) ² , and F _h = w = h − h _obs .

We have investigated the relation between the sensitiv- ity problem and the inverse problem. By solving d sensi- tivity problems with F u = e ⁱ δ(x − x), i = 1, . . ., d, to obtain their adjoint solutions (v ⁱ , q ⁱ , ψ ⁱ ) 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 gra- dient method, then we solve (19) directly with F _u = u − u _obs or F _h = h − h _obs to obtain ˜v without computing d vectors v ⁱ . Taking δC = −f (Tu)Tu·T ˜v in (31) guarantees that δF < 0.

3.1.4 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 ana- lytical 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

(u ₁ ψ ) _x = F _h + (hψ + F u ) · u _z , z = h, 0 ≤ x ≤ L. (32) The velocity from the forward equation is u(x, z) = (u ₁ , u ₃ ) ^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 ₁ is observed and let g(x) = 1 if h is observed. Then the adjoint solution is

ψ (x) =



 



 



− g(x ∗ ) u 1 (x) exp



− Z x

∗

x

h · u _z (y) u 1 (y) dy

 , 0 ≤ x ≤ x ∗ ,

0, x ∗ < x ≤ L.

(33)

So, this solution has a jump −g(x ∗ )/u ₁ (x ∗ ) at x ∗ .

With a small h · u _z (y) ≈ 0 in (33), an approximate solu-

tion is ψ (x) ≈ −g(x ∗ )/u ₁ (x). If u ₁ 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 ₁ is observed then ψ (x) is negligible. The situation is different when h is ob- served and ψ ≈ 1/u ₁ (x ∗ ) in (33).

3.2 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 ef- fect 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 ψ t + u · ∇ψ + 2ηB(u) : D(v) − ρgH ∇ · v

+ ρgv · ∇b = F _h , in ω, 0 ≤ t ≤ T , ψ (x, T ) = 0, in ω, ψ(x, t) = 0, on γ w ,

∇ · ˜ ς (v) − Cf (u)(I + F ω (u))v − H ∇ψ

= − F u , in ω,

t · ˜ ς (v)n = −C _γ f _γ (t · u)(1 + F γ (t · u))t · v, on γ _g ,

t · ˜ ς (v)n = 0, on γ w ,

n · v = 0, on γ , (34)

where the adjoint viscosity ˜η and adjoint stress ˜ ς are (cf. 20 for the case of FS)

η(u) = η(u) ˜



I + 1 − n

nB(u) : D(u) B(u) ⊗ D(u)

 ,

ς (v) = 2H ˜η(u) ? B(v). ˜ (35)

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) = kuk ^m−1 and f _γ (t · u) = |t · u| ^m−1 (cf. also Appendix A1),

F _ω (u) = m − 1

u · u u ⊗ u, F _γ = m − 1.

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)

δL =

T

Z

0

Z

ω

(2ηB(u) : D(v) + ρgv · ∇h + ∇ψ · u) δb

− f (u)u · v δC dx dt −

T

Z

0

Z

γ

_g

f _γ (t · u)t · u t

· v δC _γ ds dt. (36)

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 6= 0 and F _u = 0 in (34). Then the adjoint elevation equation must be solved for ψ 6= 0 to have a v 6= 0 in the adjoint stress equation and a perturbation in the Lagrangian in (36). The same result fol- lows from the adjoint FS equations. If F _h 6= 0 and F _u = 0 in (19), then ψ 6= 0. Consequently, v 6= 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.

In a broader context, it is worth emphasizing that the ad- joint equation derived in MacAyeal (1993) is identical to the stress equation in (34), if H is constant, F ω = 0 (e.g., m = 1), and ˜η(u) = η(u).

3.2.1 The two-dimensional adjoint solution

The adjoint SSA equations in two vertical dimensions are de- rived from (34) in the same manner as (11), by letting ψ and v 1 be independent of y and setting v ₂ = 0 and C _γ = 0. To simplify the notation, we also let v = v ₁ . The adjoint equa- tions for v and ψ follow either from simplifying (34) or from the Lagrangian and (11) and read as follows:

ψ _t + uψ _x + (ηu _x − ρgH )v _x + ρgb _x v = F _h , 0 ≤ t ≤ T , 0 ≤ x ≤ L,

ψ (x, T ) = 0, ψ (L, t ) = 0, ( 1

n ηH v _x ) _x − Cmf (u)v − H ψ _x = − F _u ,

v(0, t ) = 0, v(L, t ) = 0. (37)

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):

δL =

T

Z

0 L

Z

0

(ψ x u + v x ηu x + vρgh x ) δb

− vf (u)u δC dx dt. (38)

The weights or sensitivity functions w b and w C multiplying δb and δC in the integral are defined by

w b (x, t ) = ψ x u + v x ηu x + vρgh x ,

w _C (x, t ) = −vf (u)u. (39)

The steady-state solutions to the system (37) can be ana- lyzed as in the forward equations in Sect. 2.2.2 after simpli- fications. The viscosity terms in (37) are often small and can hence be neglected, and we assume that the basal topogra- phy is characterized by a small spatial gradient b _x . The ad- vantage 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 an- alytical approximations are less accurate close to the ice di- vide where some of the above assumptions are not valid. The adjoint equation (37) reduce to

uψ x − ρgH v x = F h , 0 ≤ x ≤ x GL , ψ x (0) = 0, ψ (x _GL ) = 0,

− Cmu ^m−1 v − H ψ x = − F u ,

v(0) = 0. (40)

3.2.2 The two-dimensional adjoint steady-state solution with velocity observation

In this section, the analytical solution to the adjoint equa- tion (40) is discussed. The derivation of the solution is de- tailed in Appendix E to Appendix F. It is here sufficient to recall that the solution given below is derived under the as- sumptions that b _x  H _x , and that a and C are constants.

For observations of u at x ∗ ,

F =

L

Z

0

u(x)δ(x − x ∗ ) dx = u ∗ , F _u = δ(x − x ∗ ), F _h = 0,

and the adjoint solutions are ψ (x) = Ca ^m x ∗

ρgH ∗ ^m+3

x ^m _GL − x ^m
, x ∗ < x ≤ x _GL ,

ψ (x) = − 1 H ∗

+ Ca ^m x ∗

ρgH ∗ ^m+3

x _GL ^m − x _∗ ^m
, 0 ≤ x < x ∗ , v(x) = ax ∗

ρgH ∗ ^m+3

H ^m , x ∗ < x ≤ x GL ,

v(x) = 0, 0 ≤ x < x ∗ , (41)

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):

δu ∗ = δL =

x

_GL

Z

0

w _b δb + w _C δC dx

=

x

_GL

Z

0

(ψ _x u + v _x ηu _x + vρgh _x ) δb − vu ^m δC dx

=

x

_GL

Z

x

⁻∗

ax ∗ H ^m H ∗ ^m+3

[(m + 1)H _x H(x − x ∗ )

+ H δ(x − x ∗ )] δb − ax ∗ (ax) ^m ρgH ∗ ^m+3

δC dx

= u ∗ δb ∗

H ∗

− u ∗

ρgH ∗ ^m+2 x

GL

Z

x

∗

C(ax) ^m



(m + 1) δb H + δC

C



dx, (42)

or, after scaling with u ∗ , δu ∗

u ∗

= δb ∗

H ∗

− 1

ρgH ∗ ^m+2 x

GL

Z

x

∗

C(ax) ^m



(m + 1) δb H + δC

C



dx. (43)

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 distribu- tion 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

δu 1i = X

j

W Cij δC j .

The elements W _Cij of the transfer matrix correspond to

quadrature coefficients in the discretization of the first inte-

gral 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 ₁ at the surface. The same analysis is performed for the other combinations of δb, δC and δu ₁ , δh.

Finally, let us comment on other approaches to investigate the sensitivity of surface data to changes in b and C, e.g., us- ing 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 perturba- tions with long wavelength λ and small wave number k are propagated to the surface, but short wavelengths are effec- tively damped in Gudmundsson (2008). The transfer func- tions 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 pertur- bations δ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 de- pendence of δu ∗ and δh ∗ on δC and δb. For example, it fol- lows 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 dif- ferent 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 veloc- ity 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 veloc- ity 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 pertur- bations 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 pertur- bation 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 perturba- tion δb is localized at each x ∗ (see Eq. 42).

iv For an unperturbed basal topography, two different per- turbations 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 perturba- tion of the friction coefficient be a constant δC ₀ 6= 0 in [ x ₀ , x ₁ ] ∈ [ x ∗ , x _GL ] and evaluate the integral in (42) to obtain

δu ∗ = − u ∗

ρgH ∗ ^m+2 x

1

Z

x

0

(ax) ^m δC 0 dx

= − a ^m u ∗

(m + 1)ρgH ∗ ^m+2

(x ₁ ^m+1 − x ₀ ^m+1 )δC ₀ . (44)

The same δu ∗ is observed with a constant perturbation in [x ₂ , x 3 ] ∈ [ x ∗ , x GL ] with the amplitude δC ₀ (x ₁ ^m+1 − x ₀ ^m+1 )/(x ₃ ^m+1 − x ₂ ^m+1 ).

v A rapidly varying friction coefficient at the base of the ice sheet will be difficult to identify by observing the ve- locity 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 num- ber 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 λ.