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.

1 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; Rignot 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 1950s – see, 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 (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 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 (steady-state) adjoint to the momentum balance in the FS model. Rather, it is necessary to include the time-dependent advection equation for the ice surface elevation in the 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.

2 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⁻² to 10⁻³. The aspect ratio is used to introduce simplifications to the FS equations, resulting, 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 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.

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 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_{\mathrm{1}},u_{\mathrm{2}},$ and u₃ the velocity components of $\mathit{u}=(u_{\mathrm{1}},u_{\mathrm{2}},u_{\mathrm{3}}{)}^{\mathrm{T}}$ in the x,y, and z direction, where $\mathit{x}=(x,y,z{)}^{\mathrm{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.

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

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 $\mathrm{\Gamma }={\mathrm{\Gamma }}_{\mathrm{u}}\cup {\mathrm{\Gamma }}_{\mathrm{d}}$. Obviously, if x∈Γ_u, then $(x,y)\in {\mathit{\gamma }}_{\mathrm{u}}$, or if x∈Γ_d, then $(x,y)\in {\mathit{\gamma }}_{\mathrm{d}}$, where $\mathit{\gamma }={\mathit{\gamma }}_{\mathrm{u}}\cup {\mathit{\gamma }}_{\mathrm{d}}$. Letting n be the outward-pointing normal on Γ (or γ in two dimensions (x,y)), the nature of ice flow renders the conditions $\mathit{n}\cdot \mathit{u}\le \mathrm{0}$ at Γ_u and $\mathit{n}\cdot \mathit{u}>\mathrm{0}$ at Γ_d. For a two-dimensional flow line geometry (Fig. 1a), this corresponds to $\mathit{x}=(x,z{)}^{\mathrm{T}},\mathit{\omega }=[\mathrm{0},L],{\mathit{\gamma }}_{\mathrm{u}}=\mathrm{0},$ and γ_d=L, where L is the horizontal length of the domain. In summary, the domains are defined as

$$\begin{array}{}\text{(1)}& \begin{array}{rl}\mathrm{\Omega }& =\mathit{\left\{}\mathit{x}\right|(x,y)\in \mathit{\omega },z_b(x,y,t)\le z\le h(x,y,t)\mathit{\}},\\ {\mathrm{\Gamma }}_{\mathrm{s}}& =\mathit{\left\{}\mathit{x}\right|(x,y)\in \mathit{\omega },z=h(x,y,t)\mathit{\}},\\ {\mathrm{\Gamma }}_{\mathrm{b}}& =\mathit{\left\{}\mathit{x}\right|(x,y)\in \mathit{\omega },z=z_b(x,y,t)\\ & =b(x,y),x<x_{\mathrm{GL}}\left(y\right)\mathit{\}},\\ {\mathrm{\Gamma }}_{\mathrm{w}}& =\mathit{\left\{}\mathit{x}\right|(x,y)\in \mathit{\omega },z=z_b(x,y,t),x>x_{\mathrm{GL}}\left(y\right)\mathit{\}},\\ {\mathrm{\Gamma }}_{\mathrm{u}}& =\mathit{\left\{}\mathit{x}\right|(x,y)\in {\mathit{\gamma }}_{\mathrm{u}},z_b(x,y,t)\le z\le h(x,y,t)\mathit{\}},\\ {\mathrm{\Gamma }}_{\mathrm{d}}& =\mathit{\left\{}\mathit{x}\right|(x,y)\in {\mathit{\gamma }}_{\mathrm{d}},z_b(x,y,t)\le z\le h(x,y,t)\mathit{\}}.\end{array}\end{array}$$

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, $\mathit{h}=(h_x,h_y,-\mathrm{1}{)}^{\mathrm{T}}$ describes the spatial gradient of the ice surface (in two vertical dimensions $\mathit{h}=(h_x,-\mathrm{1}{)}^{\mathrm{T}}$). The strain rate D and the viscosity η are given by

$$\begin{array}{}\text{(2)}& \begin{array}{r}\mathbf{D}\left(\mathit{u}\right)=\frac{\mathrm{1}}{\mathrm{2}}(\mathrm{\nabla }\mathit{u}+\mathrm{\nabla }{\mathit{u}}^{\mathrm{T}}),\mathit{\eta }\left(\mathit{u}\right)=\frac{\mathrm{1}}{\mathrm{2}}{A}^{-\frac{\mathrm{1}}{n}}\left(\mathrm{tr}{\mathbf{D}}^{\mathrm{2}}\right(\mathit{u}){)}^{\frac{\mathrm{1}-n}{\mathrm{2}n}},\end{array}\end{array}$$

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

$$\begin{array}{}\text{(3)}& \mathit{\sigma }(\mathit{u},p)=\mathrm{2}\mathit{\eta }\mathbf{D}\left(\mathit{u}\right)-\mathbf{I}p,\end{array}$$

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 $\mathbf{T}\mathit{\sigma }\mathit{n}=-\mathit{\beta }\mathbf{T}\mathit{u}$, where a projection (Petra et al., 2012) on the tangential plane of Γ_b is denoted by $\mathbf{T}=\mathbf{I}-\mathit{n}\otimes \mathit{n}$, where the Kronecker outer product between two vectors a and c or two matrices A and C is defined by

$$(\mathit{a}\otimes \mathit{c}{)}_{ij}=a_i{c}_j,(\mathbf{A}\otimes \mathbf{C}{)}_{ijkl}=A_{ij}{C}_{kl}.$$

The friction coefficient has a general form $\mathit{\beta }(\mathit{u},\mathit{x},t)=C(\mathit{x},t)f\left(\mathit{u}\right)$, where C(x,t) is independent of the velocity u and f(u) represents some linear or nonlinear function of u. For instance, $f\left(\mathit{u}\right)=\Vert \mathit{u}{\Vert }^{m-\mathrm{1}}$ with the norm $\Vert \mathit{u}\Vert =(\mathit{u}\cdot \mathit{u}{)}^{\frac{\mathrm{1}}{\mathrm{2}}}$ introduces a Weertman-type friction law (Weertman, 1957) on ω with a Weertman friction coefficient $C(\mathit{x},t)\ge \mathrm{0}$ and an exponent parameter m>0. Common choices of m are $\frac{\mathrm{1}}{\mathrm{3}}$ 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

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& h_t+\mathit{h}\cdot \mathit{u}=a_{\mathrm{s}},\mathrm{on}{\mathrm{\Gamma }}_{\mathrm{s}},t\ge \mathrm{0},\\ & h(\mathit{x},\mathrm{0})=h_{\mathrm{0}}\left(\mathit{x}\right),\mathit{x}\in \mathit{\omega },h(\mathit{x},t)=h_{\mathit{\gamma }}(\mathit{x},t),\mathit{x}\in {\mathit{\gamma }}_{\mathrm{u}},\\ & -\mathrm{\nabla }\cdot \mathit{\sigma }(\mathit{u},p)=-\mathrm{\nabla }\cdot \left(\mathrm{2}\mathit{\eta }\right(\mathit{u}\left)\mathbf{D}\right(\mathit{u}\left)\right)+\mathrm{\nabla }p=\mathit{\rho }\mathit{g},\\ & \mathrm{\nabla }\cdot \mathit{u}=\mathrm{0},\mathrm{in}\mathrm{\Omega }\left(t\right),\\ & \mathit{\sigma }\mathit{n}=\mathbf{0},\mathrm{on}{\mathrm{\Gamma }}_{\mathrm{s}},\\ & \mathbf{T}\mathit{\sigma }\mathit{n}=-Cf\left(\mathbf{T}\mathit{u}\right)\mathbf{T}\mathit{u},\mathit{n}\cdot \mathit{u}=\mathrm{0},\mathrm{on}{\mathrm{\Gamma }}_{\mathrm{b}},\end{array}\end{array}$$

where $h_{\mathrm{0}}\left(\mathit{x}\right)>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:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& z_{bt}+u{z}_{bx}+v{z}_{by}-w=a_{\mathrm{b}},\mathrm{on}{\mathrm{\Gamma }}_{\mathrm{w}},t\ge \mathrm{0},\\ & \mathit{\sigma }\mathit{n}=-p_{\mathrm{w}}\mathit{n},\mathrm{on}{\mathrm{\Gamma }}_{\mathrm{w}},\end{array}\end{array}$$

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_{\mathrm{w}}=-{\mathit{\rho }}_{\mathrm{w}}g{z}_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(\mathit{x},t)=z_b(\mathit{x},t)-b\left(\mathit{x}\right).$$

On Γ_w, we have $d(\mathit{x},t)>\mathrm{0}$ and $\mathit{n}\cdot \mathit{\sigma }\mathit{n}+p_{\mathrm{w}}=\mathrm{0}$ and on Γ_b, we have $d(\mathit{x},t)=\mathrm{0}$ and $\mathit{n}\cdot \mathit{\sigma }\mathit{n}+p_{\mathrm{w}}<\mathrm{0}$. The complementarity relation at the ice base Γ_b∪Γ_w can be summarized as

$$d(\mathit{x},t)\ge \mathrm{0},\mathit{n}\cdot \mathit{\sigma }\mathit{n}+p_{\mathrm{w}}\le \mathrm{0},d(\mathit{x},t)(\mathit{n}\cdot \mathit{\sigma }\mathit{n}+p_{\mathrm{w}})=\mathrm{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

$$\begin{array}{}\text{(6)}& \mathit{u}{|}_{{\mathrm{\Gamma }}_{\mathrm{u}}}={\mathit{u}}_{\mathrm{u}},\mathit{u}{|}_{{\mathrm{\Gamma }}_{\mathrm{d}}}={\mathit{u}}_{\mathrm{d}},\end{array}$$

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 $\mathit{u}{|}_{{\mathrm{\Gamma }}_{\mathrm{u}}}=\mathbf{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.

2.2 Shallow-shelf approximation

The three-dimensional FS problem (4) in Ω can be reduced to a two-dimensional, horizontal problem with $\mathit{x}=(x,y)\in \mathit{\omega }$ by adopting the SSA, in which only $\mathit{u}=(u_{\mathrm{1}},u_{\mathrm{2}}{)}^{\mathrm{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

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& \mathit{\eta }\left(\mathit{u}\right)=\\ & {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{A}^{-\frac{\mathrm{1}}{n}}{\left(u_{\mathrm{1}x}^{\mathrm{2}}+u_{\mathrm{2}y}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}(u_{\mathrm{1}y}+u_{\mathrm{2}x}{)}^{\mathrm{2}}+u_{\mathrm{1}x}{u}_{\mathrm{2}y}\right)}^{\frac{\mathrm{1}-n}{\mathrm{2}n}}\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{A}^{-\frac{\mathrm{1}}{n}}{\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathbf{B}:\mathbf{D}\right)}^{\frac{\mathrm{1}-n}{\mathrm{2}n}},\end{array}\end{array}$$

where $\mathbf{B}\left(\mathit{u}\right)=\mathbf{D}\left(\mathit{u}\right)+\mathrm{\nabla }\cdot \mathit{u}\mathbf{I}$, with $\mathrm{\nabla }\cdot \mathit{u}=\mathrm{tr}\mathbf{D}\left(\mathit{u}\right)$. 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 $\mathbf{A}:\mathbf{C}={\sum }_{ij}{\mathbf{A}}_{ij}{\mathbf{C}}_{ij}$. The vertically integrated stress tensor ς(u) (cf. 3 for the FS model) is given by

$$\begin{array}{}\text{(8)}& \mathit{\varsigma }\left(\mathit{u}\right)=\mathrm{2}H\mathit{\eta }\mathbf{B}\left(\mathit{u}\right),\end{array}$$

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 $\mathit{\beta }(\mathit{u},\mathit{x},t)=C(\mathit{x},t)f\left(\mathit{u}\right)=C(\mathit{x},t)\Vert \mathit{u}{\Vert }^{m-\mathrm{1}}$ with a friction coefficient $C(\mathit{x},t)\ge \mathrm{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

$$\begin{array}{}\text{(9)}& \begin{array}{rl}& H_t+\mathrm{\nabla }\cdot \left(\mathit{u}H\right)=a,t\ge \mathrm{0},\mathit{x}\in \mathit{\omega },\\ & h(\mathit{x},\mathrm{0})=h_{\mathrm{0}}\left(\mathit{x}\right),\mathit{x}\in \mathit{\omega },h(\mathit{x},t)=h_{\mathit{\gamma }}(\mathit{x},t),\mathit{x}\in {\mathit{\gamma }}_{\mathrm{u}},\\ & \mathrm{\nabla }\cdot \mathit{\varsigma }-Cf\left(\mathit{u}\right)\mathit{u}=\mathit{\rho }gH\mathrm{\nabla }h,\mathit{x}\in \mathit{\omega },\\ & \mathit{n}\cdot \mathit{u}(\mathit{x},t)=u_{\mathrm{u}}(\mathit{x},t),\mathit{x}\in {\mathit{\gamma }}_{\mathrm{u}},\mathit{n}\cdot \mathit{u}(\mathit{x},t)=u_{\mathrm{d}}(\mathit{x},t),\\ & \mathit{x}\in {\mathit{\gamma }}_{\mathrm{d}},\\ & \mathit{t}\cdot \mathit{\varsigma }\mathit{n}=-C_{\mathit{\gamma }}{f}_{\mathit{\gamma }}(\mathit{t}\cdot \mathit{u})\mathit{t}\cdot \mathit{u},\mathit{x}\in {\mathit{\gamma }}_g,\\ & \mathit{t}\cdot \mathit{\varsigma }\mathit{n}=\mathrm{0},\mathit{x}\in {\mathit{\gamma }}_{\mathrm{w}}.\end{array}\end{array}$$

Above, t is the tangential vector on $\mathit{\gamma }={\mathit{\gamma }}_{\mathrm{u}}\cup {\mathit{\gamma }}_{\mathrm{d}}$ such that $\mathit{n}\cdot \mathit{t}=\mathrm{0}$ and $a=a_{\mathrm{s}}-a_{\mathrm{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 $\mathit{\gamma }={\mathit{\gamma }}_g\cup {\mathit{\gamma }}_{\mathrm{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

$$H_{\mathrm{f}}=-{\mathit{\rho }}_{\mathrm{w}}{z}_b/\mathit{\rho }.$$

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

$$\begin{array}{}\text{(10)}& \mathit{\varsigma }\left(\mathit{u}\right)\cdot \mathit{n}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\rho }g{H}^{\mathrm{2}}\left(\mathrm{1}-{\displaystyle \frac{\mathit{\rho }}{{\mathit{\rho }}_{\mathrm{w}}}}\right)\mathit{n},\mathit{x}\in {\mathit{\gamma }}_{\mathrm{d}}.\end{array}$$

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 derived 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 ${\mathrm{\Gamma }}_{\mathrm{b}}=[\mathrm{0},x_{\mathrm{GL}}],{\mathrm{\Gamma }}_{\mathrm{w}}=(x_{\mathrm{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₁. The forward equations thus become

$$\begin{array}{}\text{(11)}& \begin{array}{rl}& H_t+(uH{)}_x=a,\mathrm{0}\le t\le T,\mathrm{0}\le x\le L,\\ & h(x,\mathrm{0})=h_{\mathrm{0}}\left(x\right),h(\mathrm{0},t)=h_L\left(t\right),\\ & (H\mathit{\eta }{u}_x{)}_x-Cf(u)u-\mathit{\rho }gH{h}_x=\mathrm{0},\\ & u(\mathrm{0},t)=u_{\mathrm{u}}\left(t\right),u(L,t)=u_{\mathrm{d}}\left(t\right),\end{array}\end{array}$$

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{\left[{\displaystyle \frac{\mathit{\rho }gH(L,t)}{\mathrm{4}}}(\mathrm{1}-{\displaystyle \frac{\mathit{\rho }}{{\mathit{\rho }}_{\mathrm{w}}}})\right]}^n.$$

Assuming that u>0 and u_x>0, the viscosity becomes $\mathit{\eta }=\mathrm{2}{A}^{-\frac{\mathrm{1}}{n}}{u}_x^{\frac{\mathrm{1}-n}{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 longitudinal stress can be ignored in the steady-state solution (see Schoof (2007)). With a sliding law in the form $f\left(u\right)=u^{m-\mathrm{1}}$ and the thickness given at x_GL, (11) thus reduces to

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& (uH{)}_x=a,\mathrm{0}\le x\le x_{\mathrm{GL}},\\ & H\left(x_{\mathrm{GL}}\right)=H_{\mathrm{GL}},\\ & -C{u}^m-\mathit{\rho }gH{h}_x=\mathrm{0},\\ & u\left(\mathrm{0}\right)=\mathrm{0}.\end{array}\end{array}$$

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

$$\begin{array}{}\text{(13)}& \begin{array}{rl}& H\left(x\right)={\left(H_{\mathrm{GL}}^{m+\mathrm{2}}+{\displaystyle \frac{m+\mathrm{2}}{m+\mathrm{1}}}{\displaystyle \frac{C{a}^m}{\mathit{\rho }g}}(x_{\mathrm{GL}}^{m+\mathrm{1}}-x^{m+\mathrm{1}})\right)}^{\frac{\mathrm{1}}{m+\mathrm{2}}},\\ & \mathrm{0}\le x\le x_{\mathrm{GL}},\\ & H\left(x\right)=H_{\mathrm{GL}},x_{\mathrm{GL}}<x<L,\\ & u\left(x\right)={\displaystyle \frac{ax}{H}},\mathrm{0}\le x\le x_{\mathrm{GL}},u\left(x\right)={\displaystyle \frac{ax}{H_{\mathrm{GL}}}},\\ & x_{\mathrm{GL}}<x<L.\end{array}\end{array}$$

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=\mathrm{1.6}\times {\mathrm{10}}^{\mathrm{6}}$ m on a single slope bed defined by $b\left(x\right)=\mathrm{720}-\frac{\mathrm{778.5}x}{\mathrm{7.5}\times {\mathrm{10}}^{\mathrm{5}}}$ and lifts from it at the grounding line position $x_{\mathrm{GL}}=\mathrm{1.035}\times {\mathrm{10}}^{\mathrm{6}}$ m. As x approaches x_GL, H decreases to become H_GL at x_GL in Fig. 2b.

Table 1 The model parameters.

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

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

3 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 ℱ:

$$\begin{array}{}\text{(14)}& \mathcal{F}=\underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathrm{\Gamma }}_{\mathrm{s}}}{\int }F(\mathit{u},h)\text{d}\mathit{x}\text{d}t.\end{array}$$

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 ℱ, the forward solution $(\mathit{u},p,h)$ to (4) or (u,h) to (9), and the adjoint solution $(\mathit{v},q,\mathit{\psi })$ or (v,ψ) to the adjoint FS and adjoint SSA equations (both derived in the following and in Appendix A), we introduce a Lagrangian $\mathcal{L}(\mathit{u},p,h;\mathit{v},q,\mathit{\psi };b,C)$. The Lagrangian for the FS equations is

$$\begin{array}{}\text{(15)}& \begin{array}{rl}& \mathcal{L}(\mathit{u},p,h;\mathit{v},q,\mathit{\psi };C)=\\ & \underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathrm{\Gamma }}_{\mathrm{s}}}{\int }F(\mathit{u},h)+\mathit{\psi }(h_t+\mathit{h}\cdot \mathit{u}-a)\text{d}\mathit{x}\text{d}t\\ & +\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathit{\omega }}{\int }\underset{b}{\overset{h}{\int }}-\mathit{v}\cdot (\mathrm{\nabla }\cdot \mathit{\sigma }(\mathit{u},p\left)\right)-\mathit{\rho }\mathit{g}\cdot \mathit{v}\\ & -q\mathrm{\nabla }\cdot \mathit{u}\text{d}\mathit{x}\text{d}t,\end{array}\end{array}$$

with the Lagrange multipliers $\mathit{v},q,$ and ψ corresponding to the forward equations for $\mathit{u},p,$ and h. The effect of perturbations δC and δb in C and b on ℱ is given by the perturbation δℒ, viz.

$$\begin{array}{}\text{(16)}& \begin{array}{r}\mathit{\delta }\mathcal{F}=\mathit{\delta }\mathcal{L}=\mathcal{L}(\mathit{u}+\mathit{\delta }\mathit{u},p+\mathit{\delta }p,h+\mathit{\delta }h;\mathit{v}+\mathit{\delta }\mathit{v},q+\mathit{\delta }q,\\ \mathit{\psi }+\mathit{\delta }\mathit{\psi };b+\mathit{\delta }b,C+\mathit{\delta }C)-\mathcal{L}(\mathit{u},p,h;\mathit{v},q,\mathit{\psi };b,C).\end{array}\end{array}$$

Examples of F(u,h) in (14) are $F=\Vert \mathit{u}-{\mathit{u}}_{\mathrm{obs}}{\Vert }^{\mathrm{2}}$, and $F=|h-h_{\mathrm{obs}}{|}^{\mathrm{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). Another example is $F(\mathit{u},h)=\frac{\mathrm{1}}{T}{u}_{\mathrm{1}}(\mathit{x},t)\mathit{\delta }(\mathit{x}-{\mathit{x}}_{*})$ with the Dirac delta δ at x_* to measure the time-averaged horizontal velocity u₁ at x_* on the ice surface Γ_s with

$$\begin{array}{}\text{(17)}& \mathcal{F}=\underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathrm{\Gamma }}_{\mathrm{s}}}{\int }F(\mathit{u},h)\text{d}\mathit{x}\text{d}t={\displaystyle \frac{\mathrm{1}}{T}}\underset{\mathrm{0}}{\overset{T}{\int }}{u}_{\mathrm{1}}({\mathit{x}}_{*},t)\text{d}t,\end{array}$$

where T is the duration of the observation at Γ_s. If the horizontal velocity is observed at $(x_{*},t_{*})$ then $F(\mathit{u},h)=u_{\mathrm{1}}(\mathit{x},t)\mathit{\delta }(\mathit{x}-{\mathit{x}}_{*})\mathit{\delta }(t-t_{*})$ and

$$\begin{array}{}\text{(18)}& \mathcal{F}=\underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathrm{\Gamma }}_{\mathrm{s}}}{\int }F(\mathit{u},h)\text{d}\mathit{x}\text{d}t=u_{\mathrm{1}}({\mathit{x}}_{*},t_{*}).\end{array}$$

The sensitivity in ℱ 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₁ 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, ℱ is the objective function in a minimization procedure and $F=\Vert \mathit{u}-{\mathit{u}}_{\mathrm{obs}}{\Vert }^{\mathrm{2}}$. The change δℱ in ℱ is of interest when C and b are changed during the solution procedure. In order to minimize ℱ, δC and δb are chosen such that δℱ<0 and ℱ 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 δℱ→0 and ℱ 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 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 ℒ for the FS equations is given in (15) and (A15) in Appendix A. To determine $(\mathit{v},q,\mathit{\psi })$ in (15), the following adjoint problem is solved:

$$\begin{array}{}\text{(19)}& \begin{array}{rl}& {\mathit{\psi }}_t+\mathrm{\nabla }\cdot \left(\mathit{u}\mathit{\psi }\right)-\mathit{h}\cdot {\mathit{u}}_z\mathit{\psi }=F_h+F_{\mathit{u}}\cdot {\mathit{u}}_z,\mathrm{on}{\mathrm{\Gamma }}_{\mathrm{s}},\\ & \mathrm{0}\le t\le T,\\ & \mathit{\psi }(\mathit{x},T)=\mathrm{0},\mathit{\psi }(\mathit{x},t)=\mathrm{0},\mathrm{on}{\mathrm{\Gamma }}_{\mathrm{d}},\\ & -\mathrm{\nabla }\cdot \stackrel{\mathrm{̃}}{\mathit{\sigma }}(\mathit{v},q)=-\mathrm{\nabla }\cdot \left(\mathrm{2}\stackrel{\mathrm{̃}}{\mathit{\eta }}\right(\mathit{u})\star \mathbf{D}(\mathit{v}\left)\right)+\mathrm{\nabla }q=\mathbf{0},\\ & \mathrm{\nabla }\cdot \mathit{v}=\mathrm{0},\mathrm{in}\mathrm{\Omega }\left(t\right),\\ & \stackrel{\mathrm{̃}}{\mathit{\sigma }}(\mathit{v},q)\mathit{n}=-(F_{\mathit{u}}+\mathit{\psi }\mathit{h}),\mathrm{on}{\mathrm{\Gamma }}_{\mathrm{s}},\\ & \mathbf{T}\stackrel{\mathrm{̃}}{\mathit{\sigma }}(\mathit{v},q)\mathit{n}=-Cf\left(\mathbf{T}\mathit{u}\right)\left(\mathbf{I}+{\mathbf{F}}_b\left(\mathbf{T}\mathit{u}\right)\right)\mathbf{T}\mathit{v},\mathrm{on}{\mathrm{\Gamma }}_{\mathrm{b}},\\ & \mathit{n}\cdot \mathit{v}=\mathrm{0},\mathrm{on}{\mathrm{\Gamma }}_{\mathrm{b}},\end{array}\end{array}$$

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

$$F_{\mathit{u}}={\left({\displaystyle \frac{\partial F}{\partial u_{\mathrm{1}}}},{\displaystyle \frac{\partial F}{\partial u_{\mathrm{2}}}},{\displaystyle \frac{\partial F}{\partial u_{\mathrm{3}}}}\right)}^{\mathrm{T}},F_h={\displaystyle \frac{\partial F}{\partial 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 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 δℱ=0 in (22); see below.

The adjoint viscosity and adjoint stress are

$$\begin{array}{}\text{(20)}& \begin{array}{rl}\stackrel{\mathrm{̃}}{\mathit{\eta }}\left(\mathit{u}\right)& =\mathit{\eta }\left(\mathit{u}\right)\left(\mathcal{I}+\frac{\mathrm{1}-n}{n\mathbf{D}\left(\mathit{u}\right):\mathbf{D}\left(\mathit{u}\right)}\mathbf{D}\left(\mathit{u}\right)\otimes \mathbf{D}\left(\mathit{u}\right)\right),\\ \stackrel{\mathrm{̃}}{\mathit{\sigma }}(\mathit{v},q)& =\mathrm{2}\stackrel{\mathrm{̃}}{\mathit{\eta }}\left(\mathit{u}\right)\star \mathbf{D}\left(\mathit{v}\right)-q\mathbf{I};\end{array}\end{array}$$

cf. also Petra et al. (2012). For the rank four-tensor ℐ, ℐ_ijkl=1 only when $i=j=k=l$; otherwise ℐ_ijkl=0. The ⋆ operation in (20) between a rank four-tensor 𝒜 and a rank two-tensor (viz. a matrix) C is defined by $(\mathcal{A}\star \mathbf{C}{)}_{ij}={\sum }_{kl}{\mathcal{A}}_{ijkl}{C}_{kl}$. In general, F_b(Tu) in (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, $f\left(\mathbf{T}\mathit{u}\right)=\Vert \mathbf{T}\mathit{u}{\Vert }^{m-\mathrm{1}}$,

$$\begin{array}{}\text{(21)}& {\mathbf{F}}_b\left(\mathbf{T}\mathit{u}\right)={\displaystyle \frac{m-\mathrm{1}}{\mathbf{T}\mathit{u}\cdot \mathbf{T}\mathit{u}}}\left(\mathbf{T}\mathit{u}\right)\otimes \left(\mathbf{T}\mathit{u}\right).\end{array}$$

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

$$\begin{array}{}\text{(22)}& \mathit{\delta }\mathcal{F}=\mathit{\delta }\mathcal{L}=\underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathrm{\Gamma }}_{\mathrm{b}}}{\int }f\left(\mathbf{T}\mathit{u}\right)\mathbf{T}\mathit{u}\cdot \mathbf{T}\mathit{v}\mathit{\delta }C\text{d}\mathit{x}\text{d}t.\end{array}$$

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

3.1.1 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 $u_{\mathrm{1}\ast }=u_{\mathrm{1}}({\mathit{x}}_{*},t_{*})$ is observed at $({\mathit{x}}_{*},t_{*})$ at the ice surface as in (18) and that $t_{*}<T$, then

$$\begin{array}{}\text{(23)}& u_{\mathrm{1}*}=\mathcal{F}=\underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathrm{\Gamma }}_{\mathrm{s}}}{\int }F\left(\mathit{u}\right)\text{d}\mathit{x}\text{d}t,\end{array}$$

with $F\left(\mathit{u}\right)=u_{\mathrm{1}}\mathit{\delta }(\mathit{x}-{\mathit{x}}_{*})\mathit{\delta }(t-t_{*}),F_{u_{\mathrm{1}}}=\mathit{\delta }(\mathit{x}-{\mathit{x}}_{*})\mathit{\delta }(t-t_{*}),F_{u_{\mathrm{2}}}=F_{u_{\mathrm{3}}}=\mathrm{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 $({\mathit{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 $\mathit{\psi }(\mathit{x},T)=\mathrm{0}$ as the corresponding final condition. Obviously, the solution for $t_{*}<t\le T$ is $\mathit{\psi }(\mathit{x},t)=\mathrm{0}$ and $\mathit{v}(\mathit{x},t)=\mathbf{0}$. Letting eⁱ denote the unit vector with 1 in the ith component, the boundary condition in (19) becomes $\stackrel{\mathrm{̃}}{\mathit{\sigma }}(\mathit{v},q)\mathit{n}=-{\mathit{e}}^{\mathrm{1}}\mathit{\delta }(\mathit{x}-{\mathit{x}}_{*})\mathit{\delta }(t-t_{*})-\mathit{\psi }\mathit{h}$ at $t=t_{*}$. For $t<t_{*}$, $\stackrel{\mathrm{̃}}{\mathit{\sigma }}(\mathit{v},q)\mathit{n}=-\mathit{\psi }\mathit{h}$. Since ψ is small for $t<t_{*}$ (see Sect. 3.1.4), the dominant part of the solution is $\mathit{v}(\mathit{x},t)={\mathit{v}}_{\mathrm{0}}\left(\mathit{x}\right)\mathit{\delta }(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 Γ_s is, according to (22), given by

$$\begin{array}{}\text{(24)}& \begin{array}{rl}\mathit{\delta }{u}_{\mathrm{1}*}& =\mathit{\delta }\mathcal{L}=\underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathrm{\Gamma }}_{\mathrm{b}}}{\int }f\left(\mathbf{T}\mathit{u}\right)\mathbf{T}\mathit{u}\cdot \mathbf{T}\mathit{v}\mathit{\delta }C\text{d}\mathit{x}\text{d}t\\ & \approx \underset{{\mathrm{\Gamma }}_{\mathrm{b}}}{\int }f\left(\mathbf{T}\mathit{u}\right)\mathbf{T}\mathit{u}\cdot \mathbf{T}{\mathit{v}}_{\mathrm{0}}\mathit{\delta }C(\mathit{x},t_{*})\text{d}\mathit{x}.\end{array}\end{array}$$

This implies that the perturbation δu_1∗ mainly depends on δC at time t_* and that contributions from previous $\mathit{\delta }C(\mathit{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

$$\begin{array}{rl}& F\left(h\right)=h(x,t)\mathit{\delta }(\mathit{x}-{\mathit{x}}_{*})\mathit{\delta }(t-t_{*}),F_h=\mathit{\delta }(\mathit{x}-{\mathit{x}}_{*})\mathit{\delta }(t-t_{*}),\\ & F_{\mathit{u}}=\mathbf{0}.\end{array}$$

The solution of the adjoint equation (19) with $\stackrel{\mathrm{̃}}{\mathit{\sigma }}(\mathit{v},q)\mathit{n}=-\mathit{\psi }\mathit{h}$ at Γ_s for v(x,t) is non-zero since $\mathit{\psi }(\mathit{x},t)\ne \mathrm{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 $\mathit{\delta }C(\mathit{x},t)=\mathit{\delta }{C}_{\mathrm{0}}\left(\mathit{x}\right)\cos (\mathrm{2}\mathit{\pi }t/\mathit{\tau })$, viz. a time-dependent perturbation added to a stationary time average C(x), with $\mathrm{0}<\mathit{\delta }{C}_{\mathrm{0}}\le 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\mathit{\tau },n=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{\dots }$ (winter) and $(n+\mathrm{1}/\mathrm{2})\mathit{\tau }$, $n=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{\dots }$ (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

$$\begin{array}{}\text{(25)}& \begin{array}{ll}\mathit{\delta }{h}_{*}& =\mathit{\delta }\mathcal{L}={\int }_{\mathrm{0}}^T\underset{{\mathrm{\Gamma }}_{\mathrm{b}}}{\int }f\left(\mathbf{T}\mathit{u}\right)\mathbf{T}\mathit{u}\cdot \mathbf{T}\mathit{v}\mathit{\delta }C(\mathit{x},t)\text{d}\mathit{x}\text{d}t\\ & \approx \mathcal{J}{\int }_{\mathrm{0}}^{t_{*}}\cos (\mathrm{2}\mathit{\pi }t/\mathit{\tau })\text{d}t=\mathcal{J}\frac{\mathit{\tau }}{\mathrm{2}\mathit{\pi }}\sin (\mathrm{2}\mathit{\pi }{t}_{*}/\mathit{\tau }),\end{array}\end{array}$$

where the integral 𝒥 is

$$\begin{array}{}\text{(26)}& \mathcal{J}=\underset{{\mathrm{\Gamma }}_{\mathrm{b}}}{\int }f\left(\mathbf{T}\mathit{u}\right)\mathbf{T}\mathit{u}\cdot \mathbf{T}\mathit{v}\mathit{\delta }{C}_{\mathrm{0}}\text{d}\mathit{x}.\end{array}$$

Obviously, from the properties of the cosine function, the friction perturbation δC is large at $t_{*}=\mathrm{0},\mathit{\tau }/\mathrm{2},\mathit{\tau },\mathrm{\dots },$ and vanishes at $t_{*}=\mathit{\tau }/\mathrm{4},\mathrm{3}\mathit{\tau }/\mathrm{4},\mathrm{\dots }$. Yet, (25) shows that $\mathit{\delta }{h}_{*}=\mathrm{0}$ at $t_{*}=\mathrm{0},\mathit{\tau },\mathrm{\dots }$ (during maximal friction in the winter) and at $t_{*}=\mathit{\tau }/\mathrm{2},\mathrm{3}\mathit{\tau }/\mathrm{2},\mathrm{\dots }$ (during minimal friction in the summer), while $\mathit{\delta }{h}_{*}\ne \mathrm{0}$ when δC=0 at $t_{*}=\mathit{\tau }/\mathrm{4},\mathrm{3}\mathit{\tau }/\mathrm{4},\mathrm{\dots }$ in the spring and the fall. The response in h by changing C is delayed in phase by $\mathit{\pi }/\mathrm{2}$ or in time by $\mathit{\tau }/\mathrm{4}=\mathrm{0.25}$ years. This is in contrast to the observation of u₁ 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 $\mathit{\delta }h(\mathit{x},t_{*})=\mathrm{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 $\mathit{\delta }C(\mathit{x},t)=\mathit{\delta }{C}_{\mathrm{0}}\left(\mathit{x}\right)s\left(t\right)$, where s(t)=1 in the time interval [t₀,t₁] and 0 otherwise. Then with 𝒥 in (26), δh_* in (25) is

$$\mathit{\delta }{h}_{*}\approx \mathcal{J}\underset{\mathrm{0}}{\overset{t_{*}}{\int }}s\left(t\right)\text{d}t=\left\{\begin{array}{ll}\mathrm{0},& t_{*}\le t_{\mathrm{0}},\\ \mathcal{J}(t_{*}-t_{\mathrm{0}}),& t_{\mathrm{0}}<t_{*}<t_{\mathrm{1}},\\ \mathcal{J}(t_{\mathrm{1}}-t_{\mathrm{0}}),& t_{*}>t_{\mathrm{1}}.\end{array}\right.$$

The effect of the basal perturbation successively increases in the elevation when $t_{*}>t_{\mathrm{0}}$ and stays at a higher level for $t_{*}>t_{\mathrm{1}}$ when δC has vanished.

3.1.2 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, $\mathit{\delta }C(x,t)=\mathrm{0.01}C\cos \left(\mathrm{2}\mathit{\pi }t\right)$ in $x\in [\mathrm{0.9},\mathrm{1.0}]\times {\mathrm{10}}^{\mathrm{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₁ and h are observed at $x\in [\mathrm{0.85},\mathrm{1.02}]\times {\mathrm{10}}^{\mathrm{6}}$ m. The steady-state solution of the forward equation with the GL located at $x_{\mathrm{GL}}=\mathrm{1.035}\times {\mathrm{10}}^{\mathrm{6}}$ m is perturbed by δu₁ and δh when C is perturbed by δC as expressed in formulas $\mathit{\delta }{u}_{\mathrm{1}}=u_{\mathrm{1}}(C+\mathit{\delta }C)-u_{\mathrm{1}}\left(C\right)$ and $\mathit{\delta }h=h(C+\mathit{\delta }C)-h\left(C\right)$. 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₁ and δh in the grounded part of the ice sheet, specifically at $x_{*}=\mathrm{0.85},\mathrm{0.9},\mathrm{0.95},\mathrm{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 $[\mathrm{0.9},\mathrm{1.0}]\times {\mathrm{10}}^{\mathrm{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₀ in (24) is negative, and an increase in the friction, δC>0, leads to a decrease in the velocity, and δC<0 increases the velocity in all panels of Fig. 3. The velocities δu₁ and the surface elevations δh are separated by a phase shift in time, $\mathrm{\Delta }\mathit{\varphi }=\mathit{\pi }/\mathrm{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_{*}=\mathrm{0.9}\times {\mathrm{10}}^{\mathrm{6}}$ to 1.0×10⁶, explaining why the shift changes sign in the red dashed lines shown in the two lower panels of Fig. 3.

Figure 3 Observations at $x_{*}=\mathrm{0.85},\mathrm{0.9},\mathrm{0.95},\mathrm{1.0},\mathrm{1.02}\times {\mathrm{10}}^{\mathrm{6}}$ m with FS in time $t\in [\mathrm{0},\mathrm{10}]$ of δu₁ (solid blue) and δh (dashed red) with perturbation δC(t)=0.01Ccos (2πt) for $x\in [\mathrm{0.9},\mathrm{1.0}]\times {\mathrm{10}}^{\mathrm{6}}$ m. Notice the different scales on the y axes.

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 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 $({\mathit{v}}^i,q^i,{\mathit{\psi }}^i),i=\mathrm{1},\mathrm{\dots },d,$ be the steady-state solution to (19) when u_i is observed at $\overline{\mathit{x}}$ and $F_{\mathit{u}}={\mathit{e}}^i\mathit{\delta }(\mathit{x}-\stackrel{\mathrm{‾}}{\mathit{x}})$. These are solutions to the sensitivity problem. We will show that the adjoint solution and the variation δℱ of the inverse problem can be expressed in $({\mathit{v}}^i,q^i,{\mathit{\psi }}^i)$. The perturbation δC is chosen such that δℱ<0 in each step in the iterative solution of the inverse problem. Then the objective function ℱ decreases stepwise toward the minimum.

It is shown in Appendix B that

$$\left(\underset{\mathit{\omega }}{\int }\sum _{i=\mathrm{1}}^d{w}_i\left(\overline{\mathit{x}}\right){\mathit{v}}^i\text{d}\stackrel{\mathrm{‾}}{\mathit{x}},\underset{\mathit{\omega }}{\int }\sum _{i=\mathrm{1}}^d{w}_i\left(\stackrel{\mathrm{‾}}{\mathit{x}}\right)q^i\text{d}\stackrel{\mathrm{‾}}{\mathit{x}},\underset{\mathit{\omega }}{\int }\sum _{i=\mathrm{1}}^d{w}_i\left(\stackrel{\mathrm{‾}}{\mathit{x}}\right){\mathit{\psi }}^i\text{d}\stackrel{\mathrm{‾}}{\mathit{x}}\right)$$

is a solution of (19) with arbitrary weights $w_i\left(\stackrel{\mathrm{‾}}{\mathit{x}}\right),i=\mathrm{1},\mathrm{\dots },d,$ when

$$\begin{array}{}\text{(27)}& F_{\mathit{u}}=\underset{\mathit{\omega }}{\int }\sum _{i=\mathrm{1}}^d{w}_i\left(\stackrel{\mathrm{‾}}{\mathit{x}}\right){\mathit{e}}^i\mathit{\delta }(\mathit{x}-\stackrel{\mathrm{‾}}{\mathit{x}})\text{d}\stackrel{\mathrm{‾}}{\mathit{x}}=\sum _{i=\mathrm{1}}^d{w}_i\left(\stackrel{\mathrm{‾}}{\mathit{x}}\right){\mathit{e}}^i.\end{array}$$

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

$$\begin{array}{}\text{(28)}& \mathit{\delta }\mathcal{F}=\underset{{\mathrm{\Gamma }}_{\mathrm{b}}}{\int }f\left(\mathbf{T}\mathit{u}\right)\mathbf{T}\mathit{u}\cdot \mathbf{T}\left(\underset{\mathit{\omega }}{\int }\sum _{i=\mathrm{1}}^d{w}_i\left(\stackrel{\mathrm{‾}}{\mathit{x}}\right){\mathit{v}}^i\text{d}\stackrel{\mathrm{‾}}{\mathit{x}}\right)\mathit{\delta }C\text{d}\mathit{x}.\end{array}$$

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

$$\begin{array}{}\text{(29)}& \mathcal{F}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\underset{\mathit{\omega }}{\int }\Vert \mathit{u}\left(\mathit{x}\right)-{\mathit{u}}_{\mathrm{obs}}\left(\mathit{x}\right){\Vert }^{\mathrm{2}}\text{d}\mathit{x},F_{\mathit{u}}=\mathit{u}\left(\mathit{x}\right)-{\mathit{u}}_{\mathrm{obs}}\left(\mathit{x}\right).\end{array}$$

The weights in (27) for the inverse problem are $w_i\left(\mathit{x}\right)=u_i\left(\mathit{x}\right)-u_{\mathrm{obs},i}\left(\mathit{x}\right)$. Let $\stackrel{\mathrm{̃}}{\mathit{v}}$ denote a weighted sum of the solutions of the sensitivity problem vⁱ over the whole domain ω:

$$\begin{array}{}\text{(30)}& \stackrel{\mathrm{̃}}{\mathit{v}}\left(\mathit{x}\right)=\underset{\mathit{\omega }}{\int }\sum _{i=\mathrm{1}}^d\left(u_i\right(\stackrel{\mathrm{‾}}{\mathit{x}})-u_{\mathrm{obs},i}(\stackrel{\mathrm{‾}}{\mathit{x}}\left)\right){\mathit{v}}^i\text{d}\stackrel{\mathrm{‾}}{\mathit{x}}.\end{array}$$

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

$$\begin{array}{}\text{(31)}& \mathit{\delta }\mathcal{F}=\underset{{\mathrm{\Gamma }}_{\mathrm{b}}}{\int }f\left(\mathbf{T}\mathit{u}\right)\mathbf{T}\mathit{u}\cdot \mathbf{T}\stackrel{\mathrm{̃}}{\mathit{v}}\left(\mathit{x}\right)\mathit{\delta }C\text{d}\mathit{x}.\end{array}$$

The same construction of the solution is possible when h_obs is given. Then d=1, $F\left(h\right)=\frac{\mathrm{1}}{\mathrm{2}}(h-h_{\mathrm{obs}}{)}^{\mathrm{2}}$, and $F_h=w=h-h_{\mathrm{obs}}$.

We have investigated the relation between the sensitivity problem and the inverse problem. By solving d sensitivity problems with $F_{\mathit{u}}={\mathit{e}}^i\mathit{\delta }(\mathit{x}-\stackrel{\mathrm{‾}}{\mathit{x}}),i=\mathrm{1},\mathrm{\dots },d,$ to obtain their adjoint solutions $({\mathit{v}}^i,q^i,{\mathit{\psi }}^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 ℱ of a change in C as in (31). In practice, if we are interested in solving the inverse problem and determine δℱ in (28) in order to iteratively compute the optimal solution with a gradient method, then we solve (19) directly with $F_{\mathit{u}}=\mathit{u}-{\mathit{u}}_{\mathrm{obs}}$ or $F_h=h-h_{\mathrm{obs}}$ to obtain $\stackrel{\mathrm{̃}}{\mathit{v}}$ without computing d vectors vⁱ. Taking $\mathit{\delta }C=-f\left(\mathbf{T}\mathit{u}\right)\mathbf{T}\mathit{u}\cdot \mathbf{T}\stackrel{\mathrm{̃}}{\mathit{v}}$ in (31) guarantees that δℱ<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 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

$$\begin{array}{}\text{(32)}& (u_{\mathrm{1}}\mathit{\psi }{)}_x=F_h+(\mathit{h}\mathit{\psi }+F_{\mathit{u}})\cdot {\mathit{u}}_z,z=h,\mathrm{0}\le x\le L.\end{array}$$

The velocity from the forward equation is $\mathit{u}(x,z)=(u_{\mathrm{1}},u_{\mathrm{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₁ is observed and let g(x)=1 if h is observed. Then the adjoint solution is

$$\begin{array}{}\text{(33)}& \mathit{\psi }\left(x\right)=\left\{\begin{array}{l}{\displaystyle }-{\displaystyle \frac{g\left(x_{*}\right)}{u_{\mathrm{1}}\left(x\right)}}\exp \left(-{\int }_x^{x_{*}}{\displaystyle \frac{\mathit{h}\cdot {\mathit{u}}_z\left(y\right)}{u_{\mathrm{1}}\left(y\right)}}\text{d}y\right),\\ \mathrm{0}\le x\le x_{*},\\ \mathrm{0},x_{*}<x\le L.\end{array}\right.\end{array}$$

So, this solution has a jump $-g\left(x_{*}\right)/u_{\mathrm{1}}\left(x_{*}\right)$ at x_*.

With a small $\mathit{h}\cdot {\mathit{u}}_z\left(y\right)\approx \mathrm{0}$ in (33), an approximate solution is $\mathit{\psi }\left(x\right)\approx -g\left(x_{*}\right)/u_{\mathrm{1}}\left(x\right)$. If u₁ is observed and $g\left(\mathit{x}\right)=u_{\mathrm{1}z}\approx \mathrm{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_{\mathrm{1}z}(x,h)=\mathrm{0}$ (Greve and Blatter, 2009; Hutter, 1983). When these approximations are accurate then u_1z will be small. Consequently, when u₁ 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 observed and $\mathit{\psi }\approx \mathrm{1}/u_{\mathrm{1}}\left(x_{*}\right)$ in (33).

3.2 Adjoint equations based on SSA

Starting from (9), a Lagrangian ℒ 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 ℒ 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 δℒ; for details, see Appendix A.

The adjoint SSA equations read

$$\begin{array}{}\text{(34)}& \begin{array}{rl}& {\mathit{\psi }}_t+\mathit{u}\cdot \mathrm{\nabla }\mathit{\psi }+\mathrm{2}\mathit{\eta }\mathbf{B}\left(\mathit{u}\right):\mathbf{D}\left(\mathit{v}\right)-\mathit{\rho }gH\mathrm{\nabla }\cdot \mathit{v}\\ & +\mathit{\rho }g\mathit{v}\cdot \mathrm{\nabla }b=F_h,\mathrm{in}\mathit{\omega },\mathrm{0}\le t\le T,\\ & \mathit{\psi }(\mathit{x},T)=\mathrm{0},\mathrm{in}\mathit{\omega },\mathit{\psi }(\mathit{x},t)=\mathrm{0},\mathrm{on}{\mathit{\gamma }}_{\mathrm{w}},\\ & \mathrm{\nabla }\cdot \stackrel{\mathrm{̃}}{\mathit{\varsigma }}\left(\mathit{v}\right)-Cf\left(\mathit{u}\right)(\mathbf{I}+{\mathbf{F}}_{\mathit{\omega }}(\mathit{u}\left)\right)\mathit{v}-H\mathrm{\nabla }\mathit{\psi }\\ & =-F_{\mathit{u}},\mathrm{in}\mathit{\omega },\\ & \mathit{t}\cdot \stackrel{\mathrm{̃}}{\mathit{\varsigma }}\left(\mathit{v}\right)\mathit{n}=-C_{\mathit{\gamma }}{f}_{\mathit{\gamma }}(\mathit{t}\cdot \mathit{u})(\mathrm{1}+F_{\mathit{\gamma }}(\mathit{t}\cdot \mathit{u}\left)\right)\mathit{t}\cdot \mathit{v},\\ & \mathrm{on}{\mathit{\gamma }}_g,\\ & \mathit{t}\cdot \stackrel{\mathrm{̃}}{\mathit{\varsigma }}\left(\mathit{v}\right)\mathit{n}=\mathrm{0},\mathrm{on}{\mathit{\gamma }}_{\mathrm{w}},\\ & \mathit{n}\cdot \mathit{v}=\mathrm{0},\mathrm{on}\mathit{\gamma },\end{array}\end{array}$$

where the adjoint viscosity $\stackrel{\mathrm{̃}}{\mathit{\eta }}$ and adjoint stress $\stackrel{\mathrm{̃}}{\mathit{\varsigma }}$ are (cf. 20 for the case of FS)

$$\begin{array}{}\text{(35)}& \begin{array}{rl}& \stackrel{\mathrm{̃}}{\mathit{\eta }}\left(\mathit{u}\right)=\mathit{\eta }\left(\mathit{u}\right)\left(\mathcal{I}+{\displaystyle \frac{\mathrm{1}-n}{n\mathbf{B}\left(\mathit{u}\right):\mathbf{D}\left(\mathit{u}\right)}}\mathbf{B}\left(\mathit{u}\right)\otimes \mathbf{D}\left(\mathit{u}\right)\right),\\ & \stackrel{\mathrm{̃}}{\mathit{\varsigma }}\left(\mathit{v}\right)=\mathrm{2}H\stackrel{\mathrm{̃}}{\mathit{\eta }}\left(\mathit{u}\right)\star \mathbf{B}\left(\mathit{v}\right).\end{array}\end{array}$$

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\left(\mathit{u}\right)=\Vert \mathit{u}{\Vert }^{m-\mathrm{1}}$ and $f_{\mathit{\gamma }}(\mathit{t}\cdot \mathit{u})=|\mathit{t}\cdot \mathit{u}{|}^{m-\mathrm{1}}$ (cf. also Appendix A1),

$${\mathbf{F}}_{\mathit{\omega }}\left(\mathit{u}\right)={\displaystyle \frac{m-\mathrm{1}}{\mathit{u}\cdot \mathit{u}}}\mathit{u}\otimes \mathit{u},F_{\mathit{\gamma }}=m-\mathrm{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 ℒ is changed by (note that the weight in front of δC in Eq. 36 is actually the same as in Eq. 22)

$$\begin{array}{}\text{(36)}& \begin{array}{rl}\mathit{\delta }\mathcal{L}& =\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathit{\omega }}{\int }\left(\mathrm{2}\mathit{\eta }\mathbf{B}\right(\mathit{u}):\mathbf{D}(\mathit{v})+\mathit{\rho }g\mathit{v}\cdot \mathrm{\nabla }h+\mathrm{\nabla }\mathit{\psi }\cdot \mathit{u})\mathit{\delta }b\\ & -f\left(\mathit{u}\right)\mathit{u}\cdot \mathit{v}\mathit{\delta }C\text{d}\mathit{x}\text{d}t-\underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathit{\gamma }}_g}{\int }{f}_{\mathit{\gamma }}(\mathit{t}\cdot \mathit{u})\mathit{t}\cdot \mathit{u}\mathit{t}\\ & \cdot \mathit{v}\mathit{\delta }{C}_{\mathit{\gamma }}\text{d}s\text{d}t.\end{array}\end{array}$$

The same perturbations in $\mathit{\delta }C,\mathit{\delta }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 δℒ 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 δℒ 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 adjoint equation derived in MacAyeal (1993) is identical to the stress equation in (34), if H is constant, F_ω=0 (e.g., m=1), and $\stackrel{\mathrm{̃}}{\mathit{\eta }}\left(\mathit{u}\right)=\mathit{\eta }\left(\mathit{u}\right)$.

3.2.1 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₁ be independent of y and setting v₂=0 and C_γ=0. To simplify the notation, we also let v=v₁. The adjoint equations for v and ψ follow either from simplifying (34) or from the Lagrangian and (11) and read as follows:

$$\begin{array}{}\text{(37)}& \begin{array}{rl}& {\mathit{\psi }}_t+u{\mathit{\psi }}_x+(\mathit{\eta }{u}_x-\mathit{\rho }gH)v_x+\mathit{\rho }g{b}_{x}v=F_h,\\ & \mathrm{0}\le t\le T,\mathrm{0}\le x\le L,\\ & \mathit{\psi }(x,T)=\mathrm{0},\mathit{\psi }(L,t)=\mathrm{0},\\ & ({\displaystyle \frac{\mathrm{1}}{n}}\mathit{\eta }H{v}_x{)}_x-Cmf(u)v-H{\mathit{\psi }}_x=-F_{\mathrm{u}},\\ & v(\mathrm{0},t)=\mathrm{0},v(L,t)=\mathrm{0}.\end{array}\end{array}$$

Note that the viscosity above is multiplied by a factor $\mathrm{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):

$$\begin{array}{}\text{(38)}& \begin{array}{rl}\mathit{\delta }\mathcal{L}& =\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathrm{0}}{\overset{L}{\int }}({\mathit{\psi }}_{x}u+v_x\mathit{\eta }{u}_x+v\mathit{\rho }g{h}_x)\mathit{\delta }b\\ & -vf\left(u\right)u\mathit{\delta }C\text{d}x\text{d}t.\end{array}\end{array}$$

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

$$\begin{array}{}\text{(39)}& \begin{array}{rl}& w_b(x,t)={\mathit{\psi }}_{x}u+v_x\mathit{\eta }{u}_x+v\mathit{\rho }g{h}_x,\\ & w_C(x,t)=-vf\left(u\right)u.\end{array}\end{array}$$

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\in [\mathrm{0},x_{\mathrm{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

$$\begin{array}{}\text{(40)}& \begin{array}{rl}& u{\mathit{\psi }}_x-\mathit{\rho }gH{v}_x=F_h,\mathrm{0}\le x\le x_{\mathrm{GL}},\\ & {\mathit{\psi }}_x\left(\mathrm{0}\right)=\mathrm{0},\mathit{\psi }\left(x_{\mathrm{GL}}\right)=\mathrm{0},\\ & -Cm{u}^{m-\mathrm{1}}v-H{\mathit{\psi }}_x=-F_{\mathrm{u}},\\ & v\left(\mathrm{0}\right)=\mathrm{0}.\end{array}\end{array}$$

3.2.2 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_*,

$$\mathcal{F}=\underset{\mathrm{0}}{\overset{L}{\int }}u\left(x\right)\mathit{\delta }(x-x_{*})\text{d}x=u_{*},F_{\mathrm{u}}=\mathit{\delta }(x-x_{*}),F_h=\mathrm{0},$$

and the adjoint solutions are

$$\begin{array}{}\text{(41)}& \begin{array}{rl}\mathit{\psi }\left(x\right)& ={\displaystyle \frac{C{a}^m{x}_{*}}{\mathit{\rho }g{H}_{*}^{m+\mathrm{3}}}}\left(x_{\mathrm{GL}}^m-x^m\right),x_{*}<x\le x_{\mathrm{GL}},\\ \mathit{\psi }\left(x\right)& =-{\displaystyle \frac{\mathrm{1}}{H_{*}}}+{\displaystyle \frac{C{a}^m{x}_{*}}{\mathit{\rho }g{H}_{*}^{m+\mathrm{3}}}}\left(x_{\mathrm{GL}}^m-x_{*}^m\right),\mathrm{0}\le x<x_{*},\\ v\left(x\right)& ={\displaystyle \frac{a{x}_{*}}{\mathit{\rho }g{H}_{*}^{m+\mathrm{3}}}}{H}^m,x_{*}<x\le x_{\mathrm{GL}},\\ v\left(x\right)& =\mathrm{0},\mathrm{0}\le x<x_{*},\end{array}\end{array}$$

where ψ(x) and v(x) have discontinuities at the observation point x_*. The perturbation of the Lagrangian (38) is with the Heaviside step function ℋ(x) and the Dirac delta δ(x) (cf. Appendix F):

$$\begin{array}{}\text{(42)}& \begin{array}{rl}\mathit{\delta }{u}_{*}& =\mathit{\delta }\mathcal{L}=\underset{\mathrm{0}}{\overset{x_{\mathrm{GL}}}{\int }}{w}_b\mathit{\delta }b+w_C\mathit{\delta }C\text{d}x\\ & =\underset{\mathrm{0}}{\overset{x_{\mathrm{GL}}}{\int }}({\mathit{\psi }}_{x}u+v_x\mathit{\eta }{u}_x+v\mathit{\rho }g{h}_x)\mathit{\delta }b-v{u}^m\mathit{\delta }C\text{d}x\\ & =\underset{x_{*}^{-}}{\overset{x_{\mathrm{GL}}}{\int }}{\displaystyle \frac{a{x}_{*}{H}^m}{H_{*}^{m+\mathrm{3}}}}\left[(m+\mathrm{1})H_x\mathcal{H}(x-x_{*})\right.\\ & \left.+H\mathit{\delta }(x-x_{*})\right]\mathit{\delta }b-{\displaystyle \frac{a{x}_{*}(ax{)}^m}{\mathit{\rho }g{H}_{*}^{m+\mathrm{3}}}}\mathit{\delta }C\text{d}x\\ & ={\displaystyle \frac{u_{*}\mathit{\delta }{b}_{*}}{H_{*}}}-{\displaystyle \frac{u_{*}}{\mathit{\rho }g{H}_{*}^{m+\mathrm{2}}}}\underset{x_{*}}{\overset{x_{\mathrm{GL}}}{\int }}C(ax{)}^m\\ & \left((m+\mathrm{1}){\displaystyle \frac{\mathit{\delta }b}{H}}+{\displaystyle \frac{\mathit{\delta }C}{C}}\right)\text{d}x,\end{array}\end{array}$$

or, after scaling with u_*,

$$\begin{array}{}\text{(43)}& \begin{array}{rl}{\displaystyle \frac{\mathit{\delta }{u}_{*}}{u_{*}}}& ={\displaystyle \frac{\mathit{\delta }{b}_{*}}{H_{*}}}-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }g{H}_{*}^{m+\mathrm{2}}}}\underset{x_{*}}{\overset{x_{\mathrm{GL}}}{\int }}C(ax{)}^m\\ & \left((m+\mathrm{1}){\displaystyle \frac{\mathit{\delta }b}{H}}+{\displaystyle \frac{\mathit{\delta }C}{C}}\right)\text{d}x.\end{array}\end{array}$$

The relation in (43) between the relative perturbations $\mathit{\delta }b/H,\mathit{\delta }C/C$, and $\mathit{\delta }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-\mathit{\delta }b,b+\mathit{\delta }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 ${\mathit{x}}_{\ast ,i},i=\mathrm{1},\mathrm{2},\mathrm{\dots }$ due to perturbations δC_j at discrete points ${\mathit{x}}_j,j=\mathrm{1},\mathrm{2},\mathrm{\dots }$ 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

$$\mathit{\delta }{u}_{\mathrm{1}i}=\sum _j{W}_{Cij}\mathit{\delta }{C}_j.$$

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₁ 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., 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.

Figure 4 The analytical solutions of ψ in (40) of the observations of (a) u and (b) h at different locations $x_{*}=\mathrm{0.25}\times {\mathrm{10}}^{\mathrm{6}},\mathrm{0.5}\times {\mathrm{10}}^{\mathrm{6}},\mathrm{0.7}\times {\mathrm{10}}^{\mathrm{6}}$, and 0.9×10⁶ m.

Figure 5 The analytical solutions of v in (40) of the observations of (a) u and (b) h at different locations $x_{*}=\mathrm{0.25}\times {\mathrm{10}}^{\mathrm{6}},\mathrm{0.5}\times {\mathrm{10}}^{\mathrm{6}},\mathrm{0.7}\times {\mathrm{10}}^{\mathrm{6}}$, and 0.9×10⁶ m.

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_{\mathrm{GL}})$ in Fig. 6a and (42). The same conclusion holds true for perturbations in b because w_b<0 in $(x_{*},x_{\mathrm{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.

Figure 6 The analytical solution of the weights $w_C=-v{u}^m$ on δC in (38) for (a) u and (b) h observed at $x_{*}=\mathrm{0.25}\times {\mathrm{10}}^{\mathrm{6}},\mathrm{0.5}\times {\mathrm{10}}^{\mathrm{6}},\mathrm{0.7}\times {\mathrm{10}}^{\mathrm{6}}$, and 0.9×10⁶ m.

Figure 7 The analytical solution of weights $w_b={\mathit{\psi }}_{x}u+v_x\mathit{\eta }{u}_x+v\mathit{\rho }g{h}_x$ on δb in (38) for (a) u and (b) h observed at $x_{*}=\mathrm{0.25}\times {\mathrm{10}}^{\mathrm{6}},\mathrm{0.5}\times {\mathrm{10}}^{\mathrm{6}},\mathrm{0.7}\times {\mathrm{10}}^{\mathrm{6}}$, and 0.9×10⁶ m.

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_{*}\to x_{\mathrm{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_{\mathrm{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, $\mathit{\delta }{u}_{*}\approx u_{*}\mathit{\delta }{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 in $[x_{\mathrm{0}},x_{\mathrm{1}}]\in [x_{*},x_{\mathrm{GL}}]$ and evaluate the integral in (42) to obtain

  $$\begin{array}{}\text{(44)}& \begin{array}{rl}\mathit{\delta }{u}_{*}& =-{\displaystyle \frac{u_{*}}{\mathit{\rho }g{H}_{*}^{m+\mathrm{2}}}}\underset{x_{\mathrm{0}}}{\overset{x_{\mathrm{1}}}{\int }}(ax{)}^m\mathit{\delta }{C}_{\mathrm{0}}\text{d}x\\ & =-{\displaystyle \frac{a^m{u}_{*}}{(m+\mathrm{1})\mathit{\rho }g{H}_{*}^{m+\mathrm{2}}}}(x_{\mathrm{1}}^{m+\mathrm{1}}-x_{\mathrm{0}}^{m+\mathrm{1}})\mathit{\delta }{C}_{\mathrm{0}}.\end{array}\end{array}$$

  The same δu_* is observed with a constant perturbation in $[x_{\mathrm{2}},x_{\mathrm{3}}]\in [x_{*},x_{\mathrm{GL}}]$ with the amplitude $\mathit{\delta }{C}_{\mathrm{0}}(x_{\mathrm{1}}^{m+\mathrm{1}}-x_{\mathrm{0}}^{m+\mathrm{1}})/(x_{\mathrm{3}}^{m+\mathrm{1}}-x_{\mathrm{2}}^{m+\mathrm{1}})$.

• 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 $\mathit{\delta }C=\mathit{ϵ}\cos (kx/x_{\mathrm{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 $\mathit{\lambda }=\mathrm{2}\mathit{\pi }{x}_{\mathrm{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

  $$\begin{array}{}\text{(45)}& \begin{array}{rl}\mathit{\delta }{u}_{*}& =-\underset{x_{*}}{\overset{x_{\mathrm{GL}}}{\int }}\mathit{ϵ}{\displaystyle \frac{a^{\mathrm{2}}{x}_{*}}{\mathit{\rho }g{H}_{*}^{\mathrm{4}}}}x\cos \left({\displaystyle \frac{kx}{x_{\mathrm{GL}}}}\right)\text{d}x\\ & =-\mathit{ϵ}{\displaystyle \frac{a^{\mathrm{2}}{x}_{*}}{\mathit{\rho }g{H}_{*}^{\mathrm{4}}}}{\displaystyle \frac{x_{\mathrm{GL}}^{\mathrm{2}}}{k}}\left(\sin \left(k\right)-{\displaystyle \frac{x_{*}}{x_{\mathrm{GL}}}}\sin \left({\displaystyle \frac{k{x}_{*}}{x_{\mathrm{GL}}}}\right)\right.\\ & +\left.{\displaystyle \frac{\mathrm{1}}{k}}\left(\cos \left(k\right)-\cos \left({\displaystyle \frac{k{x}_{*}}{x_{\mathrm{GL}}}}\right)\right)\right).\end{array}\end{array}$$

  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 $\mathrm{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 depends 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 $\mathit{\delta }b=\mathit{ϵ}\cos (kx/x_{\mathrm{GL}})$, then any perturbation at x_* is propagated to the surface by $\frac{u_{*}\mathit{\delta }{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.

3.2.3 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=\mathit{\delta }(x-x_{*})$, the expressions for ψ and v satisfying (40) are

$$\begin{array}{}\text{(46)}& \begin{array}{rl}\mathit{\psi }\left(x\right)& ={\displaystyle }-{\displaystyle \frac{C{a}^{m-\mathrm{1}}}{\mathit{\rho }g{H}_{*}^{m+\mathrm{1}}}}\left(x_{\mathrm{GL}}^m-x^m\right),x_{*}<x\le x_{\mathrm{GL}},\\ \mathit{\psi }\left(x\right)& ={\displaystyle }-{\displaystyle \frac{C{a}^{m-\mathrm{1}}}{\mathit{\rho }g{H}_{*}^{m+\mathrm{1}}}}\left(x_{\mathrm{GL}}^m-x_{*}^m\right),\mathrm{0}\le x<x_{*},\\ v\left(x\right)& ={\displaystyle }-{\displaystyle \frac{H^m}{\mathit{\rho }g{H}_{*}^{m+\mathrm{1}}}},x_{*}<x\le x_{\mathrm{GL}},\\ v\left(x\right)& =\mathrm{0},\mathrm{0}\le x<x_{*}.\end{array}\end{array}$$

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 $(\frac{\mathrm{1}}{n}\mathit{\eta }H{v}_x{)}_x$ is neglected in the simplified Eq. (40) but is of importance at x_*. A correction $\widehat{\mathit{\psi }}$ of ψ at x_* in (46) is therefore introduced to satisfy ${\left(\frac{\mathrm{1}}{n}\mathit{\eta }H{v}_x\right)}_x-H{\widehat{\mathit{\psi }}}_x=\mathrm{0}$. With $v_x\left(x_{*}\right)=-\mathit{\delta }(x-x_{*})/\left(\mathit{\rho }g{H}_{*}\right)$, the correction is $\widehat{\mathit{\psi }}\left(x\right)=-\mathit{\delta }(x-x_{*}){\mathit{\eta }}_{*}/\left(n\mathit{\rho }g{H}_{*}\right)$. The solution ψ is updated at each x_* in Fig. 4b, with $\widehat{\mathit{\psi }}$ 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 $\widehat{\mathit{\psi }}$:

$$\begin{array}{}\text{(47)}& \begin{array}{rl}{\displaystyle \frac{\mathit{\delta }{h}_{*}}{H_{*}}}& =\underset{x_{*}^{-}}{\overset{x_{\mathrm{GL}}}{\int }}-{\displaystyle \frac{u{\mathit{\eta }}_{*}}{n\mathit{\rho }g{H}_{*}^{\mathrm{2}}}}{\mathit{\delta }}_x(x-x_{*})\mathit{\delta }b\text{d}x\\ & +\underset{x_{*}}{\overset{x_{\mathrm{GL}}}{\int }}{\displaystyle \frac{C(ax{)}^m}{\mathit{\rho }g{H}_{*}^{m+\mathrm{2}}}}\left((m+\mathrm{1}){\displaystyle \frac{\mathit{\delta }b}{H}}+{\displaystyle \frac{\mathit{\delta }C}{C}}\right)\text{d}x\\ & ={\displaystyle \frac{a{\mathit{\eta }}_{*}}{n\mathit{\rho }g{H}_{*}^{\mathrm{2}}}}{\left(x{\displaystyle \frac{\mathit{\delta }b}{H}}\right)}_x\left(x_{*}\right)+{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }g{H}_{*}^{m+\mathrm{2}}}}\underset{x_{*}}{\overset{x_{\mathrm{GL}}}{\int }}C(ax{)}^m\\ & \left((m+\mathrm{1}){\displaystyle \frac{\mathit{\delta }b}{H}}+{\displaystyle \frac{\mathit{\delta }C}{C}}\right)\text{d}x,\end{array}\end{array}$$

where $a(x\mathit{\delta }b/H{)}_x(x_{*})=(u\mathit{\delta }b{)}_x\left(x_{*}\right)$ represents the x-derivative of uδb evaluated at x_*. When δb=0 then δu_* in (43) and $\mathit{\delta }{h}_{*}=\mathit{\delta }{H}_{*}$ in (47) satisfy $\mathit{\delta }{u}_{*}{H}_{*}=-\mathit{\delta }{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):

$$\begin{array}{}\text{(48)}& \begin{array}{rl}\mathit{\delta }{h}_{*}& =\underset{\mathrm{0}}{\overset{x_{\mathrm{GL}}}{\int }}({\mathit{\psi }}_{x}u+v_x\mathit{\eta }{u}_x+v\mathit{\rho }g{h}_x)\mathit{\delta }b-v{u}^m\mathit{\delta }C\text{d}x\\ & =\underset{\mathrm{0}}{\overset{x_{\mathrm{GL}}}{\int }}{w}_b\mathit{\delta }b+w_C\mathit{\delta }C\text{d}x.\end{array}\end{array}$$

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 $\mathit{\delta }b/H$, and the first term in (47) depends on the derivative of $ax\mathit{\delta }b/H=u\mathit{\delta }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) and the left and right columns of Figs. 6 and 7, the conclusions (i), (ii), (iv), (v), and (vi) in Sect. 3.2.2 from (42) and (43) for δu_* are valid also for δh_* in (47).

3.2.4 The two-dimensional time-dependent adjoint solution

Finally, the time-dependent adjoint equation (37) is investigated. Equation (37) is solved numerically for the same MISMIP 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):

$$\begin{array}{}\text{(49)}& C(x,t)=C_{\mathrm{0}}(\mathrm{1}+\mathit{\kappa }\cos (\mathrm{2}\mathit{\pi }t\left)\right),\mathrm{0}<\mathit{\kappa }<\mathrm{1}.\end{array}$$

Apparently, C has its highest value at $t=n,n=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{\dots }$, i.e., the winter, and its lowest value at $t=n+\mathrm{1}/\mathrm{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_{*}=\mathrm{9}\times {\mathrm{10}}^{\mathrm{5}}$ m for 0.1 year in the four seasons starting from the summer of the 10th year, e.g., in the summer $(t_{*}=\mathrm{9.5})$, the fall $(t_{*}=\mathrm{9.75})$, the winter $(t_{*}=\mathrm{10})$, and the spring $(t_{*}=\mathrm{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³ 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_{*}=\mathrm{9}\times {\mathrm{10}}^{\mathrm{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 $[\mathrm{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⁵ 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 $[\mathrm{0},t_{*})$, the smaller the effect of it is at t_*. After 5 years a perturbation can hardly be detected at the surface.

Figure 8 The adjoint weights for the observations at $x_{*}=\mathrm{9}\times {\mathrm{10}}^{\mathrm{5}}$ m of the four seasons. (a) w_C for the observation of u. (b) w_b for the observation of u. (c) w_C for the observation of h. (d) w_b for the observation of h.

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

Figure 9 The adjoint weights at x_* in the four seasons of the 10th year with seasonally varying friction coefficient. The black dashed line is a reference solution without seasonal variations which is observed at $t_{*}=\mathrm{9.75}$. (a) w_C for the observation of u. (b) w_b for the observation of u. (c) w_C for the observation of h. (d) w_b for the observation of h.

However, when h is observed, the effects of δC_* and δb_* are not visible directly because $w_{C\ast }\approx \mathrm{0}$ and $w_{b*}\approx \mathrm{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_{*}=\mathrm{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\ast }\exp (-(T-t)/\mathit{\tau })$ 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

$$\mathit{\delta }C(x,t)=\mathit{\delta }{C}_{\mathrm{0}}\cos \left(\mathrm{2}\mathit{\pi }ft\right)\mathit{\delta }(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

$$\begin{array}{}\text{(50)}& \begin{array}{rl}\mathit{\delta }{h}_{*}& =\underset{\mathrm{0}}{\overset{t_{*}}{\int }}{w}_{C\ast }\exp (-(t_{*}-t)/\mathit{\tau })\mathit{\delta }{C}_{\mathrm{0}}\cos \left(\mathrm{2}\mathit{\pi }ft\right)\text{d}t\\ & =\left({\displaystyle \frac{\cos \left(\mathrm{2}\mathit{\pi }f{t}_{*}\right)+\mathrm{2}\mathit{\pi }\mathit{\tau }f\sin \left(\mathrm{2}\mathit{\pi }f{t}_{*}\right)-e^{-t_{*}/\mathit{\tau }}}{\mathrm{4}{\mathit{\pi }}^{\mathrm{2}}\mathit{\tau }{f}^{\mathrm{2}}+{\mathit{\tau }}^{-\mathrm{1}}}}\right)w_{C*}\mathit{\delta }{C}_{\mathrm{0}};\end{array}\end{array}$$

cf. (45). If the frequency is high, f≫1, then $\mathit{\delta }{h}_{*}\propto \mathrm{1}/f$ and high-frequency perturbations are damped efficiently. At certain times of observation t_* when $\sin \left(\mathrm{2}\mathit{\pi }f{t}_{*}\right)=\mathrm{0}$, the damping is even stronger with $\mathit{\delta }{h}_{*}\propto \mathrm{1}/f^{\mathrm{2}}$. If the frequency is low, f≪1, then $\mathit{\delta }{h}_{*}\propto \mathit{\tau }$ and the change in h_* is insensitive to the frequency. The same conclusions hold true for δb where decennial perturbations seem more realistic.

4 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., $\mathit{\delta }{u}_{\mathrm{1}}=u_{\mathrm{1}}(C+\mathit{\delta }C)-u_{\mathrm{1}}\left(C\right)$ 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\ne 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.

Appendix A: Derivation of the adjoint equations

A1 Adjoint viscosity and friction in SSA

The adjoint viscosity $\stackrel{\mathrm{̃}}{\mathit{\eta }}\left(\mathit{u}\right)$ in SSA in (20) is derived as follows. The SSA viscosity for u and u+δu is

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}& \mathit{\eta }(\mathit{u}+\mathit{\delta }\mathit{u})\approx \mathit{\eta }\left(\mathit{u}\right)(\mathrm{1}+{\displaystyle \frac{\mathrm{1}-n}{\mathrm{2}n\widehat{\mathit{\eta }}}}((\mathrm{2}{u}_{\mathrm{1}x}+u_{\mathrm{2}y})\mathit{\delta }{u}_{\mathrm{1}x}\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(u_{\mathrm{1}y}+u_{\mathrm{2}x})\mathit{\delta }{u}_{\mathrm{2}x}+(\mathrm{2}{u}_{\mathrm{2}y}+u_{\mathrm{1}x})\mathit{\delta }{u}_{\mathrm{2}y}\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(u_{\mathrm{1}y}+u_{\mathrm{2}x})\mathit{\delta }{u}_{\mathrm{1}y}\left)\right).\end{array}\end{array}$$

Determine ℬ(u) such that

$$\mathit{\varrho }(\mathit{u},\mathit{\delta }\mathit{u})\mathbf{B}\left(\mathit{u}\right)=\mathcal{B}\left(\mathit{u}\right)\star \mathbf{B}\left(\mathit{\delta }\mathit{u}\right).$$

First note that

$$\begin{array}{rl}& \mathbf{B}\left(\mathit{u}\right):\mathbf{D}\left(\mathit{\delta }\mathit{u}\right)=\left(\mathbf{D}\right(\mathit{u})+\mathrm{\nabla }\cdot \mathit{u}\mathbf{I}):\mathbf{D}\left(\mathit{\delta }\mathit{u}\right)=\mathbf{D}\left(\mathit{u}\right):\mathbf{D}\left(\mathit{\delta }\mathit{u}\right)\\ & +(\mathrm{\nabla }\cdot \mathit{u})(\mathrm{\nabla }\cdot \mathit{\delta }\mathit{u})=\mathbf{D}\left(\mathit{u}\right):\left(\mathbf{B}\right(\mathit{\delta }\mathit{u})-\mathrm{\nabla }\cdot \mathit{\delta }\mathit{u}\mathbf{I})\\ & +(\mathrm{\nabla }\cdot \mathit{u})(\mathrm{\nabla }\cdot \mathit{\delta }\mathit{u})=\mathbf{D}\left(\mathit{u}\right):\mathbf{B}\left(\mathit{\delta }\mathit{u}\right).\end{array}$$

Then use the ⋆ operator to define ℬ:

$$\begin{array}{rl}& {\displaystyle \frac{\mathrm{1}-n}{\mathrm{2}n\widehat{\mathit{\eta }}}}\sum _{kl}{B}_{kl}\left(\mathit{u}\right)D_{kl}\left(\mathit{\delta }\mathit{u}\right)B_{ij}\left(\mathit{u}\right)\\ & ={\displaystyle \frac{\mathrm{1}-n}{\mathrm{2}n\widehat{\mathit{\eta }}}}\sum _{kl}{D}_{kl}\left(\mathit{u}\right)B_{kl}\left(\mathit{\delta }\mathit{u}\right)B_{ij}\left(\mathit{u}\right)\\ & =\sum _{kl}{\mathcal{B}}_{ijkl}\left(\mathit{u}\right)D_{kl}\left(\mathit{\delta }\mathit{u}\right)=(\mathcal{B}\star D{)}_{ij}.\end{array}$$

Thus, let

$$\begin{array}{rl}& {\mathcal{B}}_{ijkl}={\displaystyle \frac{\mathrm{1}-n}{\mathrm{2}n\widehat{\mathit{\eta }}}}{B}_{ij}\left(\mathit{u}\right)D_{kl}\left(\mathit{u}\right),\\ & {\stackrel{\mathrm{̃}}{\mathit{\eta }}}_{ijkl}\left(\mathit{u}\right)=\mathit{\eta }\left(\mathit{u}\right)({\mathcal{I}}_{ijkl}+{\mathcal{B}}_{ijkl}(\mathit{u}\left)\right),\end{array}$$

or in tensor form

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}& \mathcal{B}={\displaystyle \frac{\mathrm{1}-n}{n\mathbf{B}\left(\mathit{u}\right):\mathbf{D}\left(\mathit{u}\right)}}\mathbf{B}\left(\mathit{u}\right)\otimes \mathbf{D}\left(\mathit{u}\right),\\ & \stackrel{\mathrm{̃}}{\mathit{\eta }}\left(\mathit{u}\right)=\mathit{\eta }\left(\mathit{u}\right)\left(\mathcal{I}+\mathcal{B}\right).\end{array}\end{array}$$

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 $\mathit{\xi }=\mathit{u},\mathit{\zeta }=\mathit{v},c=C,$ and F=F_ω and at γ_g with $\mathit{\xi }=\mathit{t}\cdot \mathit{u},\mathit{\zeta }=\mathit{t}\cdot \mathit{v},c=C_{\mathit{\gamma }},f=f_{\mathit{\gamma }},$ and $\mathbf{F}=F_{\mathit{\gamma }},$ we arrive at the adjoint friction term cf(ξ)(I+F(ξ))ζ, where

$$\begin{array}{}\text{(A3)}& \mathbf{F}\left(\mathit{\xi }\right)={\displaystyle \frac{m-\mathrm{1}}{\mathit{\xi }\cdot \mathit{\xi }}}\mathit{\xi }\otimes \mathit{\xi }.\end{array}$$

A2 Adjoint equations in SSA

The Lagrangian for the SSA equations is with the adjoint variables $\mathit{\psi },\mathit{v},q$

$$\begin{array}{}\text{(A4)}& \begin{array}{rl}& \mathcal{L}(\mathit{u},h;\mathit{v},\mathit{\psi };b,C_{\mathit{\gamma }},C)=\\ & \underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathit{\omega }}{\int }F(\mathit{u},h)+\mathit{\psi }(h_t+\mathrm{\nabla }\cdot (\mathit{u}H)-a)\text{d}\mathit{x}\text{d}t\\ & +\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathit{\omega }}{\int }\mathit{v}\cdot \mathrm{\nabla }\cdot \left(\mathrm{2}H\mathit{\eta }\mathbf{B}\right(\mathit{u}\left)\right)-Cf\left(\mathit{u}\right)\mathit{v}\cdot \mathit{u}-\mathit{\rho }gH\mathit{v}\cdot \mathrm{\nabla }h\text{d}\mathit{x}\text{d}t\\ & =\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathit{\omega }}{\int }F(\mathit{u},h)+\mathit{\psi }(h_t+\mathrm{\nabla }\cdot (\mathit{u}H)-a)\text{d}\mathit{x}\text{d}t\\ & +\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathit{\omega }}{\int }-\mathrm{2}H\mathit{\eta }\left(\mathit{u}\right)\left(\mathbf{D}\right(\mathit{v}):\mathbf{D}(\mathit{u})+\mathrm{\nabla }\cdot \mathit{u}\mathrm{\nabla }\cdot \mathit{v})\\ & -Cf\left(\mathit{u}\right)\mathit{v}\cdot \mathit{u}-\mathit{\rho }gH\mathit{v}\cdot \mathrm{\nabla }h\text{d}\mathit{x}\text{d}t-\underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathit{\gamma }}_g}{\int }{C}_{\mathit{\gamma }}{f}_{\mathit{\gamma }}(\mathit{t}\cdot \mathit{u})\mathit{t}\\ & \cdot \mathit{u}\mathit{t}\cdot \mathit{v}\text{d}s\text{d}t\end{array}\end{array}$$

after partial integration and using the boundary conditions. The perturbed SSA Lagrangian is split into the unperturbed Lagrangian and three integrals:

$$\begin{array}{}\text{(A5)}& \begin{array}{rl}& \mathcal{L}(\mathit{u}+\mathit{\delta }\mathit{u},h+\mathit{\delta }h;\mathit{v}+\mathit{\delta }\mathit{v},\mathit{\psi }+\mathit{\delta }\mathit{\psi };b+\mathit{\delta }b,C_{\mathit{\gamma }}+\mathit{\delta }{C}_{\mathit{\gamma }},\\ & C+\mathit{\delta }C)=\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathit{\omega }}{\int }F(\mathit{u}+\mathit{\delta }\mathit{u},h+\mathit{\delta }h)\text{d}\mathit{x}\text{d}t+\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathit{\omega }}{\int }(\mathit{\psi }+\mathit{\delta }\mathit{\psi })\\ & (h_t+\mathit{\delta }{h}_t+\mathrm{\nabla }\cdot \left((\mathit{u}+\mathit{\delta }\mathit{u})(H+\mathit{\delta }H)\right)-a)\text{d}\mathit{x}\text{d}t\\ & +\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathit{\omega }}{\int }-\mathrm{2}(H+\mathit{\delta }H)\mathit{\eta }(\mathit{u}+\mathit{\delta }\mathit{u})\mathbf{D}(\mathit{v}+\mathit{\delta }\mathit{v}):\mathbf{B}(\mathit{u}+\mathit{\delta }\mathit{u})\\ & -(C+\mathit{\delta }C)f(\mathit{u}+\mathit{\delta }\mathit{u})(\mathit{u}+\mathit{\delta }\mathit{u})\cdot (\mathit{v}+\mathit{\delta }\mathit{v})-\mathit{\rho }g(H+\mathit{\delta }H)\\ & \mathrm{\nabla }(h+\mathit{\delta }h)\cdot (\mathit{v}+\mathit{\delta }\mathit{v})\text{d}\mathit{x}\text{d}t-\underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathit{\gamma }}_g}{\int }(C_{\mathit{\gamma }}+\mathit{\delta }{C}_{\mathit{\gamma }})f_{\mathit{\gamma }}\\ & (\mathit{t}\cdot (\mathit{u}+\mathit{\delta }\mathit{u}\left)\right)\mathit{t}\cdot (\mathit{u}+\mathit{\delta }\mathit{u})\mathit{t}\cdot (\mathit{v}+\mathit{\delta }\mathit{v})\text{d}s\text{d}t\\ & =\mathcal{L}(\mathit{u},h;\mathit{v},\mathit{\psi };b,C_{\mathit{\gamma }},C)+I_{\mathrm{1}}+I_{\mathrm{2}}+I_{\mathrm{3}}.\end{array}\end{array}$$

The perturbation in ℒ is

$$\begin{array}{}\text{(A6)}& \mathit{\delta }\mathcal{L}=I_{\mathrm{1}}+I_{\mathrm{2}}+I_{\mathrm{3}}.\end{array}$$

Terms of order 2 or more in δℒ are neglected. Then the first term in δℒ satisfies

$$\begin{array}{}\text{(A7)}& \begin{array}{rl}{I}_{\mathrm{1}}& =\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathit{\omega }}{\int }F(\mathit{u}+\mathit{\delta }\mathit{u},h+\mathit{\delta }h)-F(\mathit{u},h)\text{d}\mathit{x}\text{d}t\\ & =\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathit{\omega }}{\int }{F}_{\mathit{u}}\mathit{\delta }\mathit{u}+F_h\mathit{\delta }h\text{d}\mathit{x}\text{d}t.\end{array}\end{array}$$

Using partial integration, Gauss' formula, and the initial and boundary conditions on u and H and $\mathit{\psi }(\mathit{x},T)=\mathrm{0},\mathit{x}\in \mathit{\omega },$ and $\mathit{\psi }(\mathit{x},t)=\mathrm{0},\mathit{x}\in {\mathit{\gamma }}_{\mathrm{w}},$ in the second integral we have

$$\begin{array}{}\text{(A8)}& \begin{array}{rl}{I}_{\mathrm{2}}& =\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathit{\omega }}{\int }\mathit{\delta }\mathit{\psi }(h_t+\mathrm{\nabla }\cdot (\mathit{u}H)-a)\\ & +\mathit{\psi }(\mathit{\delta }{h}_t+\mathrm{\nabla }\cdot (\mathit{\delta }\mathit{u}H)+\mathrm{\nabla }\cdot (\mathit{u}\mathit{\delta }H\left)\right)\text{d}\mathit{x}\text{d}t\\ & =\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathit{\omega }}{\int }\mathit{\delta }\mathit{\psi }(h_t+\mathrm{\nabla }\cdot (\mathit{u}H)-a)\text{d}\mathit{x}\text{d}t+\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathit{\omega }}{\int }\\ & -{\mathit{\psi }}_t\mathit{\delta }h-H\mathrm{\nabla }\mathit{\psi }\cdot \mathit{\delta }\mathit{u}-\mathrm{\nabla }\mathit{\psi }\cdot \mathit{u}\mathit{\delta }H\text{d}\mathit{x}\text{d}t.\end{array}\end{array}$$

The first integral after the second equality vanishes since h is a weak solution and I₂ is

$$\begin{array}{}\text{(A9)}& \begin{array}{rl}{I}_{\mathrm{2}}& =\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathit{\omega }}{\int }-({\mathit{\psi }}_t+\mathit{u}\cdot \mathrm{\nabla }\mathit{\psi })\mathit{\delta }h-H\mathrm{\nabla }\mathit{\psi }\cdot \mathit{\delta }\mathit{u}\\ & +\mathit{u}\cdot \mathrm{\nabla }\mathit{\psi }\mathit{\delta }b\text{d}\mathit{x}\text{d}t.\end{array}\end{array}$$

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

$$\begin{array}{}\text{(A10)}& \begin{array}{rl}{I}_{\mathrm{3}}=& I_{\mathrm{31}}+I_{\mathrm{32}},\\ I_{\mathrm{31}}=& \underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathit{\omega }}{\int }-\mathrm{2}(H+\mathit{\delta }H)\mathit{\eta }(\mathit{u}+\mathit{\delta }\mathit{u})\mathbf{D}(\mathit{v}+\mathit{\delta }\mathit{v})\\ & :\mathbf{B}(\mathit{u}+\mathit{\delta }\mathit{u})-(C+\mathit{\delta }C)f(\mathit{u}+\mathit{\delta }\mathit{u})(\mathit{u}+\mathit{\delta }\mathit{u})\\ & \cdot (\mathit{v}+\mathit{\delta }\mathit{v})\\ & -\mathit{\rho }g(H+\mathit{\delta }H)\mathrm{\nabla }(h+\mathit{\delta }h)\cdot (\mathit{v}+\mathit{\delta }\mathit{v}))\text{d}\mathit{x}\text{d}t\\ & -\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathit{\gamma }}{\int }(C_{\mathit{\gamma }}+\mathit{\delta }{C}_{\mathit{\gamma }})f_{\mathit{\gamma }}(\mathit{t}\cdot (\mathit{u}+\mathit{\delta }\mathit{u}\left)\right)\mathit{t}\\ & \cdot (\mathit{u}+\mathit{\delta }\mathit{u})\mathit{t}\cdot (\mathit{v}+\mathit{\delta }\mathit{v})\text{d}s\text{d}t\\ =& I_{\mathrm{311}}+I_{\mathrm{312}}-I_{\mathrm{313}},\end{array}\end{array}$$

where

$$\begin{array}{}\text{(A11)}& \begin{array}{rl}{I}_{\mathrm{311}}=& \underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathit{\omega }}{\int }-\mathrm{2}H\mathbf{D}\left(\mathit{v}\right):\left(\mathit{\eta }\right(\mathit{u}+\mathit{\delta }\mathit{u}\left)\mathbf{B}\right(\mathit{u}+\mathit{\delta }\mathit{u}\left)\right)+\mathrm{2}H\mathbf{D}\left(\mathit{v}\right)\\ & :\left(\mathit{\eta }\right(\mathit{u}\left)\mathbf{B}\right(\mathit{u}\left)\right)\text{d}\mathit{x}\text{d}t\\ =& \underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathit{\omega }}{\int }-\mathrm{2}H\mathbf{D}\left(\mathit{v}\right):\left(\stackrel{\mathrm{̃}}{\mathit{\eta }}\right(\mathit{u})\star \mathbf{B}(\mathit{\delta }\mathit{u}\left)\right)\text{d}\mathit{x}\text{d}t\\ I_{\mathrm{312}}=& \underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathit{\omega }}_g}{\int }-\mathit{\delta }Cf\left(\mathit{u}\right)\mathit{u}\cdot \mathit{v}\text{d}\mathit{x}\text{d}t+\underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathit{\omega }}_g}{\int }-C\left(f\right(\mathit{u}+\mathit{\delta }\mathit{u})\mathit{v}\\ & \cdot (\mathit{u}+\mathit{\delta }\mathit{u})-f\left(\mathit{u}\right)\mathit{v}\cdot \mathit{u})\text{d}\mathit{x}\text{d}t\\ =& \underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathit{\omega }}_g}{\int }-\mathit{\delta }Cf\left(\mathit{u}\right)\mathit{u}\cdot \mathit{v}+Cf\left(\mathit{u}\right)(\mathbf{I}+{\mathbf{F}}_{\mathit{\omega }}(\mathit{u}\left)\right)\mathit{\delta }\mathit{u}\cdot \mathit{v}\text{d}\mathit{x}\text{d}t\\ I_{\mathrm{313}}=& \underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathit{\gamma }}_g}{\int }(C_{\mathit{\gamma }}+\mathit{\delta }{C}_{\mathit{\gamma }})\left(f_{\mathit{\gamma }}\right(\mathit{t}\cdot (\mathit{u}+\mathit{\delta }\mathit{u}))\mathit{t}\cdot \mathit{v}\mathit{t}\cdot (\mathit{u}+\mathit{\delta }\mathit{u})\\ & -f_{\mathit{\gamma }}(\mathit{t}\cdot \mathit{u})\mathit{t}\cdot \mathit{v}\mathit{t}\cdot \mathit{u})\text{d}s\text{d}t\\ =& \underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathit{\gamma }}_g}{\int }(C_{\mathit{\gamma }}+\mathit{\delta }{C}_{\mathit{\gamma }})\left(f_{\mathit{\gamma }}\right(\mathit{t}\cdot \mathit{u})\mathit{t}\cdot \mathit{u}\mathit{t}\cdot \mathit{v}+C_{\mathit{\gamma }}{f}_{\mathit{\gamma }}(\mathit{t}\cdot \mathit{u})\\ & (\mathbf{I}+{\mathbf{F}}_{\mathit{\gamma }}(\mathit{t}\cdot \mathit{u}\left)\right)\mathit{t}\cdot \mathit{\delta }\mathit{u}\mathit{t}\cdot \mathit{v}\text{d}s\text{d}t\\ I_{\mathrm{32}}=& \underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathit{\omega }}{\int }-\mathit{\rho }gH\mathrm{\nabla }h\cdot \mathit{v}-\mathrm{2}\mathit{\eta }\mathbf{D}\left(\mathit{v}\right):\mathbf{B}\left(\mathit{u}\right)\mathit{\delta }H-\mathit{\rho }g\mathrm{\nabla }h\\ & \cdot \mathit{v}\mathit{\delta }H-\mathit{\rho }gH\mathit{v}\cdot \mathrm{\nabla }\mathit{\delta }h\text{d}\mathit{x}\text{d}t\\ =& \underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathit{\omega }}{\int }-\mathit{\rho }gH\mathrm{\nabla }h\cdot \mathit{v}-\left(\mathrm{2}\mathit{\eta }\mathbf{D}\right(\mathit{v}):\mathbf{B}(\mathit{u})+\mathit{\rho }g\mathrm{\nabla }h\cdot \mathit{v})\\ & \mathit{\delta }H+\mathit{\rho }g\mathrm{\nabla }\cdot \left(H\mathit{v}\right)\mathit{\delta }h\text{d}\mathit{x}\text{d}t.\end{array}\end{array}$$

Collecting all the terms in (A7), (A9), and (A10), the first variation of ℒ is

$$\begin{array}{}\text{(A12)}& \begin{array}{rl}\mathit{\delta }\mathcal{L}=& I_{\mathrm{1}}+I_{\mathrm{2}}+I_{\mathrm{3}}\\ =& \underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathit{\omega }}{\int }{F}_{\mathit{u}}\mathit{\delta }\mathit{u}-\mathrm{2}HD\left(\mathit{v}\right):\left(\stackrel{\mathrm{̃}}{\mathit{\eta }}\right(\mathit{u})\star \mathbf{B}(\mathit{\delta }\mathit{u}\left)\right)-H\mathrm{\nabla }\mathit{\psi }\\ & \cdot \mathit{\delta }\mathit{u}\text{d}\mathit{x}\text{d}t-\underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathit{\omega }}_g}{\int }Cf\left(\mathit{u}\right)(\mathbf{I}+{\mathbf{F}}_{\mathit{\omega }}(\mathit{u}\left)\right)\mathit{v}\cdot \mathit{\delta }\mathit{u}\text{d}\mathit{x}\text{d}t\\ & -\underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathit{\gamma }}_g}{\int }{C}_{\mathit{\gamma }}{f}_{\mathit{\gamma }}(\mathit{t}\cdot \mathit{u})(\mathbf{I}+{\mathbf{F}}_{\mathit{\gamma }}(\mathit{t}\cdot \mathit{u}\left)\right)\mathit{t}\cdot \mathit{v}\mathit{t}\cdot \mathit{\delta }\mathit{u}\text{d}s\text{d}t\\ & -\underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathit{\gamma }}_g}{\int }\mathit{\delta }{C}_{\mathit{\gamma }}{f}_{\mathit{\gamma }}(\mathit{t}\cdot \mathit{u})\mathit{t}\cdot \mathit{u}\mathit{t}\cdot \mathit{v}\text{d}s\text{d}t\\ & +\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathit{\omega }}{\int }(F_h-({\mathit{\psi }}_t+\mathit{u}\cdot \mathrm{\nabla }\mathit{\psi }+\mathrm{2}\mathit{\eta }\mathbf{D}\left(\mathit{v}\right):\mathbf{B}\left(\mathit{u}\right)\\ & -\mathit{\rho }g\mathrm{\nabla }b\cdot \mathit{v}+\mathit{\rho }gH\mathrm{\nabla }\cdot \mathit{v}\left)\right)\mathit{\delta }h\text{d}\mathit{x}\text{d}t\\ & +\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathit{\omega }}{\int }-\mathit{\delta }Cf\left(\mathit{u}\right)\mathit{v}\cdot \mathit{u}+\left(\mathrm{2}\mathit{\eta }\mathbf{D}\right(\mathit{v}):\mathbf{B}(\mathit{u})\\ & +\mathit{\rho }g\mathrm{\nabla }h\cdot \mathit{v}+\mathit{u}\cdot \mathrm{\nabla }\mathit{\psi })\mathit{\delta }b\text{d}\mathit{x}\text{d}t.\end{array}\end{array}$$

The forward solution $({\mathit{u}}^{*},p^{*},h^{*})$ and adjoint solution $({\mathit{v}}^{*},q^{*},{\mathit{\psi }}^{*})$ satisfying (9) and (34) are inserted into (A4) resulting in

$$\begin{array}{}\text{(A13)}& \begin{array}{rl}& \mathcal{L}({\mathit{u}}^{*},p^{*};{\mathit{v}}^{*},q^{*};h^{*},{\mathit{\psi }}^{*};b,C_{\mathit{\gamma }},C)\\ & =\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathit{\omega }}{\int }F({\mathit{u}}^{*},h^{*})\text{d}\mathit{x}\text{d}t.\end{array}\end{array}$$

Then (A12) yields the variation in ℒ in (A13) with respect to perturbations $\mathit{\delta }b,\mathit{\delta }{C}_{\mathit{\gamma }},$ and δC in $b,C_{\mathit{\gamma }},$ and C:

$$\begin{array}{}\text{(A14)}& \begin{array}{rl}\mathit{\delta }\mathcal{L}=& \underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathit{\omega }}{\int }\left(\mathrm{2}\mathit{\eta }\mathbf{D}\right({\mathit{v}}^{*}):\mathbf{B}({\mathit{u}}^{*})+\mathit{\rho }g\mathrm{\nabla }{h}^{*}\cdot {\mathit{v}}^{*}\\ & +{\mathit{u}}^{*}\cdot \mathrm{\nabla }{\mathit{\psi }}^{*})\mathit{\delta }b\text{d}\mathit{x}\text{d}t-\underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathit{\gamma }}_g}{\int }\mathit{\delta }{C}_{\mathit{\gamma }}{f}_{\mathit{\gamma }}(\mathit{t}\cdot {\mathit{u}}^{*})\mathit{t}\\ & \cdot {\mathit{u}}^{*}\mathit{t}\cdot {\mathit{v}}^{*}\text{d}s\text{d}t-\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathit{\omega }}{\int }\mathit{\delta }Cf\left({\mathit{u}}^{*}\right){\mathit{v}}^{*}\cdot {\mathit{u}}^{*}\text{d}\mathit{x}\text{d}t.\end{array}\end{array}$$

A3 Adjoint equations in FS

The FS Lagrangian is as in (15):

$$\begin{array}{}\text{(A15)}& \begin{array}{rl}& \mathcal{L}(\mathit{u},p,h;\mathit{v},q,\mathit{\psi };C)=\underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathrm{\Gamma }}_{\mathrm{s}}}{\int }F(\mathit{u},h)+\mathit{\psi }(h_t+\mathit{h}\\ & \cdot \mathit{u}-a_{\mathrm{s}})\text{d}\mathit{x}\text{d}t+\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathit{\omega }}{\int }\underset{b}{\overset{h}{\int }}-\mathit{v}\cdot (\mathrm{\nabla }\cdot \mathit{\sigma }(\mathit{u},p))-q\mathrm{\nabla }\\ & \cdot \mathit{u}-\mathit{\rho }\mathit{g}\cdot \mathit{v}\text{d}\mathit{x}\text{d}t=\underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathrm{\Gamma }}_{\mathrm{s}}}{\int }F(\mathit{u},h)+\mathit{\psi }(h_t+\mathit{h}\\ & \cdot \mathit{u}-a)\text{d}\mathit{x}\text{d}t+\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathit{\omega }}{\int }\underset{b}{\overset{h}{\int }}\mathrm{2}\mathit{\eta }(\mathit{u}\left)\mathbf{D}\right(\mathit{v}):\mathbf{D}(\mathit{u})-p\mathrm{\nabla }\\ & \cdot \mathit{v}-q\mathrm{\nabla }\cdot \mathit{u}-\mathit{\rho }\mathit{g}\cdot \mathit{v}\text{d}\mathit{x}\text{d}t+\underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathrm{\Gamma }}_{\mathrm{b}}}{\int }Cf\left(\mathbf{T}\mathit{u}\right)\mathbf{T}\mathit{u}\\ & \cdot \mathbf{T}\mathit{v}\text{d}\mathit{x}\text{d}t.\end{array}\end{array}$$

In the same manner as in (A5), the perturbed FS Lagrangian is

$$\begin{array}{}\text{(A16)}& \begin{array}{rl}& \mathcal{L}(\mathit{u}+\mathit{\delta }\mathit{u},p+\mathit{\delta }p;\mathit{v}+\mathit{\delta }\mathit{v},q+\mathit{\delta }q;h+\mathit{\delta }h,\mathit{\psi }+\mathit{\delta }\mathit{\psi };\\ & C+\mathit{\delta }C)=\mathcal{L}(\mathit{u},p,h;\mathit{v},q,\mathit{\psi };C)+I_{\mathrm{1}}+I_{\mathrm{2}}+I_{\mathrm{3}}.\end{array}\end{array}$$

Terms of order 2 or more in $\mathit{\delta }\mathit{u},\mathit{\delta }\mathit{v},$ and δh are neglected. The first integral I₁ in (A16) is

$$\begin{array}{}\text{(A17)}& \begin{array}{rl}{I}_{\mathrm{1}}& =\underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathrm{\Gamma }}_{\mathrm{s}}}{\int }F\left(\mathit{u}\right(\mathit{x},h+\mathit{\delta }h,t)+\mathit{\delta }\mathit{u},h+\mathit{\delta }h)\\ & -F\left(\mathit{u}\right(\mathit{x},h,t),h)\text{d}\mathit{x}\text{d}t=\underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathrm{\Gamma }}_{\mathrm{s}}}{\int }{F}_{\mathit{u}}(\mathit{\delta }\mathit{u}+{\mathit{u}}_z\mathit{\delta }h)\\ & +F_h\mathit{\delta }h\text{d}\mathit{x}\text{d}t.\end{array}\end{array}$$

Partial integration, the conditions $\mathit{\psi }(\mathit{x},T)=\mathrm{0}$ and $\mathit{\psi }(\mathit{x},t)=\mathrm{0}$ at Γ_s, and the fact that h is a weak solution simplify the second integral:

$$\begin{array}{}\text{(A18)}& \begin{array}{rl}{I}_{\mathrm{2}}& =\underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathrm{\Gamma }}_{\mathrm{s}}}{\int }\mathit{\delta }\mathit{\psi }(h_t+\mathit{h}\cdot \mathit{u}-a_{\mathrm{s}})+\mathit{\psi }(\mathit{\delta }{h}_t+\mathit{u}\cdot \mathit{\delta }\mathit{h}\\ & +{\mathit{u}}_z\cdot \mathit{h}\mathit{\delta }h+\mathit{h}\cdot \mathit{\delta }\mathit{u})\text{d}\mathit{x}\text{d}t=\underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathrm{\Gamma }}_{\mathrm{s}}}{\int }\mathit{\delta }\mathit{\psi }(h_t+\mathit{h}\cdot \mathit{u}\\ & -a_{\mathrm{s}})\text{d}\mathit{x}\text{d}t+\underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathrm{\Gamma }}_{\mathrm{s}}}{\int }(-{\mathit{\psi }}_t-\mathrm{\nabla }\cdot \left(\mathit{u}\mathit{\psi }\right)\\ & +\mathit{h}\cdot {\mathit{u}}_z\mathit{\psi })\mathit{\delta }h+\mathit{h}\cdot \mathit{\delta }\mathit{u}\mathit{\psi }\text{d}\mathit{x}\text{d}t.\end{array}\end{array}$$

Define $\mathrm{\Xi },\mathit{\xi },$ and Υ to be

$$\begin{array}{}\text{(A19)}& \begin{array}{rl}\mathrm{\Theta }(\mathit{u},p;\mathit{v},q;C)=& \mathrm{2}\mathit{\eta }\left(\mathit{u}\right)\mathbf{D}\left(\mathit{v}\right):\mathbf{D}\left(\mathit{u}\right)-p\mathrm{\nabla }\cdot \mathit{v}\\ & -q\mathrm{\nabla }\cdot \mathit{u}-\mathit{\rho }\mathit{g}\cdot \mathit{v},\\ \mathit{\theta }(\mathit{u};\mathit{v};C)=& Cf\left(\mathbf{T}\mathit{u}\right)\mathbf{T}\mathit{u}\cdot \mathbf{T}\mathit{v},\\ \mathrm{{\rm Y}}(\mathit{u},p;\mathit{v},q)=& -\mathit{v}\cdot (\mathrm{\nabla }\cdot \mathit{\sigma }(\mathit{u},p\left)\right)-q\mathrm{\nabla }\\ & \cdot \mathit{u}-\mathit{\rho }\mathit{g}\cdot \mathit{v}.\end{array}\end{array}$$

Then a weak solution, (u,p), for any (v,q) satisfying the boundary conditions, fulfills

$$\begin{array}{}\text{(A20)}& \begin{array}{rl}& \underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathit{\omega }}{\int }\underset{b}{\overset{h}{\int }}\mathrm{\Theta }(\mathit{u},p;\mathit{v},q;C)\text{d}\mathit{x}\text{d}t-\underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathrm{\Gamma }}_{\mathrm{b}}}{\int }\mathit{\theta }(\mathit{u};\mathit{v};C)\\ & \text{d}\mathit{x}\text{d}t=\mathrm{0}.\end{array}\end{array}$$

The third integral in (A16) is

$$\begin{array}{}\text{(A21)}& \begin{array}{rl}{I}_{\mathrm{3}}=& I_{\mathrm{31}}+I_{\mathrm{32}},\\ I_{\mathrm{31}}=& \underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathit{\omega }}{\int }\underset{b}{\overset{h}{\int }}\mathrm{\Theta }(\mathit{u}+\mathit{\delta }\mathit{u},p+\mathit{\delta }p;\mathit{v}+\mathit{\delta }\mathit{v},q+\mathit{\delta }q;\\ & C+\mathit{\delta }C)\text{d}\mathit{x}\text{d}t-\underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathrm{\Gamma }}_{\mathrm{b}}}{\int }\mathit{\theta }(\mathit{u}+\mathit{\delta }\mathit{u};\mathit{v}+\mathit{\delta }\mathit{v};\\ & C+\mathit{\delta }C)\text{d}\mathit{x}\text{d}t,\\ I_{\mathrm{32}}=& \underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathit{\omega }}{\int }\underset{h}{\overset{h+\mathit{\delta }h}{\int }}\mathrm{{\rm Y}}(\mathit{u},p;\mathit{v},q)\text{d}\mathit{x}\text{d}t.\end{array}\end{array}$$

The integral I₃₁ is expanded as in (A10) and (A11) or Petra et al. (2012) using the weak solution, Gauss' formula, and the definitions of the adjoint viscosity and adjoint friction coefficient in Appendix A1. When $b<z<h$ we have $\mathrm{{\rm Y}}(\mathit{u},p;\mathit{v},q)=\mathrm{0}$. If Υ is extended smoothly in the positive z direction from z=h, then with $z\in [h,h+\mathit{\delta }h]$ for some constant c>0 we have $\left|\mathrm{{\rm Y}}\right|\le c\mathit{\delta }h$. Therefore,

$$\begin{array}{rl}& \left|\underset{h}{\overset{h+\mathit{\delta }h(x,t)}{\int }}\mathrm{{\rm Y}}\right(\mathit{u},p;\mathit{v},q\left)dz\right|\le \underset{h}{\overset{h+\mathit{\delta }h(x,t)}{\int }}sup\left|\mathrm{{\rm Y}}\right|dz\\ & \le c\left|\mathit{\delta }h\right(x,t{)}^{\mathrm{2}}|,\end{array}$$

and the bound on I₃₂ in (A21) is

$$\begin{array}{}\text{(A22)}& \left|I_{\mathrm{32}}\right|\le ct\left|\mathit{\omega }\right|max\left|\mathit{\delta }h\right(x,t){|}^{\mathrm{2}},\end{array}$$

where $\left|\mathit{\omega }\right|$ is the area of ω. This term is a second variation in δh which is neglected and I₃=I₃₁.

The first variation of ℒ is then

$$\begin{array}{}\text{(A23)}& \begin{array}{rl}\mathit{\delta }\mathcal{L}=& I_{\mathrm{1}}+I_{\mathrm{2}}+I_{\mathrm{3}}\\ =& \underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathrm{\Gamma }}_{\mathrm{s}}}{\int }(F_{\mathit{u}}+\mathit{\psi }\mathit{h})\cdot \mathit{\delta }\mathit{u}\text{d}\mathit{x}\text{d}t+\underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathrm{\Gamma }}_{\mathrm{s}}}{\int }(F_h+F_{\mathit{u}}{\mathit{u}}_z\\ & -({\mathit{\psi }}_t+\mathrm{\nabla }\cdot (\mathit{u}\mathit{\psi })-\mathit{h}\cdot \mathit{u}\mathit{\psi }))\mathit{\delta }h\text{d}\mathit{x}\text{d}t\\ & +\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathit{\omega }}{\int }\underset{b}{\overset{h}{\int }}\mathrm{2}\mathbf{D}\left(\mathit{v}\right):\left(\stackrel{\mathrm{̃}}{\mathit{\eta }}\right(\mathit{u})\star \mathbf{D}(\mathit{\delta }\mathit{u}\left)\right)-\mathit{\delta }p\mathrm{\nabla }\cdot \mathit{v}\\ & -q\mathrm{\nabla }\cdot \mathit{\delta }\mathit{u}\text{d}\mathit{x}\text{d}t+\underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathrm{\Gamma }}_{\mathrm{b}}}{\int }Cf\left(\mathbf{T}\mathit{u}\right)(\mathbf{I}+{\mathbf{F}}_b(\mathit{u}\left)\right)\mathbf{T}\mathit{v}\\ & \cdot \mathbf{T}\mathit{\delta }\mathit{u}\text{d}\mathit{x}\text{d}t+\underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathrm{\Gamma }}_{\mathrm{b}}}{\int }\mathit{\delta }Cf\left(\mathbf{T}\mathit{u}\right)\mathbf{T}\mathit{u}\cdot \mathbf{T}\mathit{v}\text{d}\mathit{x}\text{d}t.\end{array}\end{array}$$

With the forward solution $({\mathit{u}}^{*},p^{*},h^{*})$ and the adjoint solution $({\mathit{v}}^{*},q^{*},{\mathit{\psi }}^{*})$ satisfying (4) and (19), the first variation with respect to perturbations δC in C is (cf. A14)

$$\begin{array}{}\text{(A24)}& \mathit{\delta }\mathcal{L}=\underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathrm{\Gamma }}_{\mathrm{b}}}{\int }f\left(\mathbf{T}{\mathit{u}}^{*}\right)\mathbf{T}{\mathit{u}}^{*}\cdot \mathbf{T}{\mathit{v}}^{*}\mathit{\delta }C\text{d}\mathit{x}\text{d}t.\end{array}$$

Appendix B: The adjoint solution in the inverse and sensitivity problems

Assume that $({\mathit{v}}^i,q^i,{\mathit{\psi }}^i),i=\mathrm{1},\mathrm{\dots },d,$ solves adjoint FS equation (19) in the steady state with observation of u_i with d=2 or 3

$$\begin{array}{}\text{(B1)}& \begin{array}{rl}& \mathcal{F}=u_i\left(\mathit{x}\right)=\underset{\mathit{\omega }}{\int }{u}_i\mathit{\delta }(\mathit{x}-\stackrel{\mathrm{‾}}{\mathit{x}})\text{d}\stackrel{\mathrm{‾}}{\mathit{x}},F_{\mathit{u}}={\mathit{e}}^i\mathit{\delta }(\mathit{x}-\stackrel{\mathrm{‾}}{\mathit{x}}),\\ & i=\mathrm{1},\mathrm{\dots },d,\end{array}\end{array}$$

or observation of h with d=1

$$\begin{array}{}\text{(B2)}& \mathcal{F}=h\left(\mathit{x}\right)=\underset{\mathit{\omega }}{\int }h\mathit{\delta }(\mathit{x}-\stackrel{\mathrm{‾}}{\mathit{x}})\text{d}\stackrel{\mathrm{‾}}{\mathit{x}},F_h=\mathit{\delta }(\mathit{x}-\stackrel{\mathrm{‾}}{\mathit{x}}).\end{array}$$

Introduce the weight functions $w_i\left(\mathit{x}\right),i=\mathrm{1},\mathrm{\dots },d$. It follows from (19) that $\left(w_i\right(\stackrel{\mathrm{‾}}{\mathit{x}}\left){\mathit{v}}^i\right(\mathit{x}),w_i(\stackrel{\mathrm{‾}}{\mathit{x}}\left)q^i\right(\mathit{x}),w_i(\stackrel{\mathrm{‾}}{\mathit{x}}\left){\mathit{\psi }}^i\right(\mathit{x}\left)\right)$ is a solution with $F_{\mathit{u}}=w_i\left(\stackrel{\mathrm{‾}}{\mathit{x}}\right){\mathit{e}}^i\mathit{\delta }(\mathit{x}-\stackrel{\mathrm{‾}}{\mathit{x}})$ or $F_h=w\left(\stackrel{\mathrm{‾}}{\mathit{x}}\right)\mathit{\delta }(\mathit{x}-\stackrel{\mathrm{‾}}{\mathit{x}})$. Therefore, also

$$\begin{array}{}\text{(B3)}& \left(\underset{\mathit{\omega }}{\int }{w}_i\left(\stackrel{\mathrm{‾}}{\mathit{x}}\right){\mathit{v}}^i\text{d}\stackrel{\mathrm{‾}}{\mathit{x}},\underset{\mathit{\omega }}{\int }{w}_i\left(\stackrel{\mathrm{‾}}{\mathit{x}}\right)q^i\text{d}\stackrel{\mathrm{‾}}{\mathit{x}},\underset{\mathit{\omega }}{\int }{w}_i\left(\stackrel{\mathrm{‾}}{\mathit{x}}\right){\mathit{\psi }}^i\text{d}\stackrel{\mathrm{‾}}{\mathit{x}}\right)\end{array}$$

is a solution with $F_{\mathit{u}}={\int }_{\mathit{\omega }}{w}_i\left(\stackrel{\mathrm{‾}}{\mathit{x}}\right){\mathit{e}}^i\mathit{\delta }(\mathit{x}-\stackrel{\mathrm{‾}}{\mathit{x}})\text{d}\stackrel{\mathrm{‾}}{\mathit{x}}=w_i\left(\mathit{x}\right){\mathit{e}}^i$ or $F_h={\int }_{\mathit{\omega }}w\left(\stackrel{\mathrm{‾}}{\mathit{x}}\right)\mathit{\delta }(\mathit{x}-\stackrel{\mathrm{‾}}{\mathit{x}})\text{d}\stackrel{\mathrm{‾}}{\mathit{x}}=w\left(\mathit{x}\right)$. A sum over $i,i=\mathrm{1},\mathrm{\dots },d,$ of each integral in (B3) is also a solution.

Consider a target functional ℱ for the steady-state solution with a weight vector $\mathit{w}\left(\stackrel{\mathrm{‾}}{\mathit{x}}\right)$ with components $w_i\left(\stackrel{\mathrm{‾}}{\mathit{x}}\right)$ multiplying δuⁱ in the first variation of ℱ. Using (22), δℱ is

$$\begin{array}{}\text{(B4)}& \begin{array}{rl}\mathit{\delta }\mathcal{F}& =\underset{\mathit{\omega }}{\int }\mathit{w}\left(\stackrel{\mathrm{‾}}{\mathit{x}}\right)\cdot \mathit{\delta }\mathit{u}\text{d}\stackrel{\mathrm{‾}}{\mathit{x}}\\ & =\underset{\mathit{\omega }}{\int }\sum _{i=\mathrm{1}}^d{w}_i\left(\stackrel{\mathrm{‾}}{\mathit{x}}\right)\mathit{\delta }{u}^i\text{d}\stackrel{\mathrm{‾}}{\mathit{x}}\\ & =\underset{\mathit{\omega }}{\int }\sum _{i=\mathrm{1}}^d{w}_i\left(\stackrel{\mathrm{‾}}{\mathit{x}}\right)\underset{{\mathrm{\Gamma }}_{\mathrm{b}}}{\int }f\left(\mathbf{T}\mathit{u}\right)\mathbf{T}\mathit{u}\cdot \mathbf{T}{\mathit{v}}^i\mathit{\delta }C\text{d}\mathit{x}\text{d}\stackrel{\mathrm{‾}}{\mathit{x}}\\ & =\underset{{\mathrm{\Gamma }}_{\mathrm{b}}}{\int }f\left(\mathbf{T}\mathit{u}\right)\mathbf{T}\mathit{u}\cdot \mathbf{T}\left(\underset{\mathit{\omega }}{\int }\sum _{i=\mathrm{1}}^d{w}_i\left(\stackrel{\mathrm{‾}}{\mathit{x}}\right){\mathit{v}}^i\text{d}\stackrel{\mathrm{‾}}{\mathit{x}}\right)\mathit{\delta }C\text{d}\mathit{x}.\end{array}\end{array}$$

Appendix C: Steady-state solution of the adjoint height equation in the FS model

In a two-dimensional vertical ice with $\mathit{u}(x,z)=(u_{\mathrm{1}},u_{\mathrm{3}}{)}^T$, the stationary equation for ψ in (19) is

$$\begin{array}{}\text{(C1)}& (u_{\mathrm{1}}\mathit{\psi }{)}_x=F_h+(\mathit{h}\mathit{\psi }+F_{\mathit{u}})\cdot {\mathit{u}}_z,z=h,\mathrm{0}\le x\le L.\end{array}$$

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₁ is observed at Γ_s then $F(\mathit{u},h)=u_{\mathrm{1}}\left(x\right)\mathit{\chi }\left(x\right)$, $F_{\mathit{u}}=\left(\mathit{\chi }\right(x),\mathrm{0}{)}^T$, and F_h=0. The weight χ on u₁ may be a Dirac delta, a Gaussian, or a constant in a limited interval. On the other hand, if $F(\mathit{u},h)=h\left(x\right)\mathit{\chi }\left(x\right)$ then F_h=χ(x) and F_u=0.

Let g(x)=u_1z(x) when F_u≠0 and let g(x)=1 when F_h≠0. Then by (32),

$$\begin{array}{}\text{(C2)}& (u_{\mathrm{1}}\mathit{\psi }{)}_x-\mathit{h}\cdot {\mathit{u}}_z\mathit{\psi }=g(x\left)\mathit{\chi }\right(x).\end{array}$$

The solution to (C2) is

$$\begin{array}{}\text{(C3)}& \begin{array}{rl}& \mathit{\psi }\left(x\right)=-{\displaystyle \frac{\mathrm{1}}{u_{\mathrm{1}}\left(x\right)}}\underset{x}{\overset{x_{*}}{\int }}\exp \left(-\underset{x}{\overset{\mathit{\xi }}{\int }}{\displaystyle \frac{\mathit{h}\cdot {\mathit{u}}_z\left(y\right)}{u_{\mathrm{1}}\left(y\right)}}\text{d}y\right)g\left(\mathit{\xi }\right)\mathit{\chi }\left(\mathit{\xi }\right)\text{d}\mathit{\xi },\\ & \mathrm{0}\le x<x_{*},\mathit{\psi }\left(x\right)=\mathrm{0},x_{*}<x\le L.\end{array}\end{array}$$

In particular, if $\mathit{\chi }\left(x\right)=\mathit{\delta }(x-x_{*})$ then $\mathcal{F}=u_{\mathrm{1}}\left(x_{*}\right)$ or $\mathcal{F}=h\left(x_{*}\right)$, and the multiplier is

$$\begin{array}{}\text{(C4)}& \begin{array}{rl}& \mathit{\psi }\left(x\right)=-{\displaystyle \frac{g\left(x_{*}\right)}{u_{\mathrm{1}}\left(x\right)}}\exp \left(-\underset{x}{\overset{x_{*}}{\int }}{\displaystyle \frac{\mathit{h}\cdot {\mathit{u}}_z\left(y\right)}{u_{\mathrm{1}}\left(y\right)}}\text{d}y\right),\\ & \mathrm{0}\le x<x_{*},\end{array}\end{array}$$

which has a jump $-g\left(x_{*}\right)/u_{\mathrm{1}}\left(x_{*}\right)$ at x_*.

Appendix D: Simplified SSA equations

The forward and adjoint SSA equations in (12) and (40) are solved analytically. The conclusion from the thickness equation in (12) is that

$$\begin{array}{}\text{(D1)}& u\left(x\right)H\left(x\right)=u\left(\mathrm{0}\right)H\left(\mathrm{0}\right)+ax=ax,\end{array}$$

since u(0)=0. Solve the second equation in (12) for u on the bedrock with x≤x_GL and insert into (D1) using the assumptions for x>0 that b_x≪H_x and h_x≈H_x to have

$$\begin{array}{}\text{(D2)}& {\displaystyle \frac{\mathit{\rho }g}{C}}{H}^{m+\mathrm{1}}{H}_x={\displaystyle \frac{\mathit{\rho }g}{C(m+\mathrm{2})}}(H^{m+\mathrm{2}}{)}_x=-(ax{)}^m.\end{array}$$

The equation for H^m+2 for x≤x_GL is integrated from x to x_GL such that

$$\begin{array}{}\text{(D3)}& \begin{array}{rl}H\left(x\right)& ={\left(H_{\mathrm{GL}}^{m+\mathrm{2}}+\frac{m+\mathrm{2}}{m+\mathrm{1}}\frac{C{a}^m}{\mathit{\rho }g}(x_{\mathrm{GL}}^{m+\mathrm{1}}-x^{m+\mathrm{1}})\right)}^{\frac{\mathrm{1}}{m+\mathrm{2}}},\\ u\left(x\right)& ={\displaystyle \frac{ax}{H}},H_x={\displaystyle }-{\displaystyle \frac{C{a}^m}{\mathit{\rho }g}}{\displaystyle \frac{x^m}{H^{m+\mathrm{1}}}}.\end{array}\end{array}$$

For the floating ice at x>x_GL, ρgHh_x=0, implying that h_x=0 and H_x=0. Hence, H(x)=H_GL. The velocity increases linearly beyond the grounding line

$$\begin{array}{}\text{(D4)}& u\left(x\right)=ax/H\left(x\right)=ax/H_{\mathrm{GL}},x>x_{\mathrm{GL}}.\end{array}$$

By including the viscosity term in (11) and assuming that H(x) is linear in x, a more accurate formula is obtained for u(x) on the floating ice in (6.77) of Greve and Blatter (2009).

Appendix E: Jumps in ψ and v in SSA

Multiply the first equation in (40) by H and the second equation by u to eliminate ψ_x. We get

$$\begin{array}{}\text{(E1)}& -Cm{u}^{m}v-\mathit{\rho }g{H}^{\mathrm{2}}{v}_x=H{F}_h-u{F}_{\mathrm{u}}.\end{array}$$

Use the expression for u and H_x in (D3). Then

$$\begin{array}{}\text{(E2)}& \mathit{\rho }gH(m{H}_{x}v-H{v}_x)=H{F}_h-u{F}_{\mathrm{u}},\end{array}$$

or equivalently

$$\begin{array}{}\text{(E3)}& {\left({\displaystyle \frac{v}{H^m}}\right)}_x=-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }g{H}^{m+\mathrm{2}}}}(H{F}_h-u{F}_{\mathrm{u}}).\end{array}$$

The solutions ψ(x) and v(x) of the adjoint SSA equation (37) have jumps at the observation point x_*. For x close to x_* in a short interval $[x_{*}^{-},x_{*}^{+}]$ with $x_{*}^{-}<x_{*}<x_{*}^{+}$, integrate (E3) to receive

$$\begin{array}{}\text{(E4)}& \underset{x_{*}^{-}}{\overset{x_{*}^{+}}{\int }}{\left({\displaystyle \frac{v}{H^m}}\right)}_x\text{d}x=-\underset{x_{*}^{-}}{\overset{x_{*}^{+}}{\int }}{\displaystyle \frac{H{F}_h-u{F}_{\mathrm{u}}}{\mathit{\rho }g{H}^{m+\mathrm{2}}}}\text{d}x.\end{array}$$

Since H is continuous and u and v are bounded, when $x_{*}^{-}\to x_{*}^{+}$, then

$$\begin{array}{}\text{(E5)}& \begin{array}{rl}& v\left(x_{*}^{+}\right)-v\left(x_{*}^{-}\right)=\\ & -{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }g{H}_{*}^{\mathrm{2}}}}\left(H_{*}\underset{x_{*}^{-}}{\overset{x_{*}^{+}}{\int }}{F}_h\text{d}x-u_{*}\underset{x_{*}^{-}}{\overset{x_{*}^{+}}{\int }}{F}_{\mathrm{u}}\text{d}x\right).\end{array}\end{array}$$

A similar relation for ψ can be derived:

$$\begin{array}{}\text{(E6)}& \mathit{\psi }\left(x_{*}^{+}\right)-\mathit{\psi }\left(x_{*}^{-}\right)={\displaystyle \frac{\mathrm{1}}{H_{*}}}\underset{x_{*}^{-}}{\overset{x_{*}^{+}}{\int }}{F}_{\mathrm{u}}\text{d}x.\end{array}$$

With F_u=0 and F_h=0 for $x<x_{*}$ and $v\left(\mathrm{0}\right)={\mathit{\psi }}_x\left(\mathrm{0}\right)=\mathrm{0}$, we find that

$$\begin{array}{}\text{(E7)}& v\left(x\right)={\mathit{\psi }}_x\left(x\right)=\mathrm{0},\mathit{\psi }\left(x\right)=\mathit{\psi }\left(x_{*}^{-}\right),\mathrm{0}\le x<x_{*}.\end{array}$$

If $F(u,h)=u\mathit{\delta }(x-x_{*})$, then by (E5) and (E6)

$$\begin{array}{}\text{(E8)}& v\left(x_{*}^{+}\right)={\displaystyle \frac{u_{*}}{\mathit{\rho }g{H}_{*}^{\mathrm{2}}}},\mathit{\psi }\left(x_{*}^{+}\right)-\mathit{\psi }\left(x_{*}^{-}\right)={\displaystyle \frac{\mathrm{1}}{H_{*}}},\end{array}$$

and if $F(u,h)=h\mathit{\delta }(x-x_{*})$, then

$$\begin{array}{}\text{(E9)}& v\left(x_{*}^{+}\right)=-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }g{H}_{*}}},\mathit{\psi }\left(x_{*}^{+}\right)-\mathit{\psi }\left(x_{*}^{-}\right)=\mathrm{0}.\end{array}$$

Appendix F: Analytical solutions in SSA

By Appendix E, v(x)=0 for $\mathrm{0}\le x<x_{\ast }$. Use equations in (40) with H_x in (D3) for $x_{*}<x\le x_{\mathrm{GL}}$ to have

$${\displaystyle \frac{v_x}{v}}=-{\displaystyle \frac{axCm{u}^{m-\mathrm{1}}}{\mathit{\rho }g{H}^{\mathrm{3}}}}=-{\displaystyle \frac{Cm{u}^m}{\mathit{\rho }g{H}^{\mathrm{2}}}}={\displaystyle \frac{m{H}_x}{H}}.$$

Let $\mathcal{H}(x-x_{*})={\int }_{-\mathrm{\infty }}^{x-x_{*}}\mathit{\delta }\left(s\right)ds$ be the Heaviside step function at x_*. Then

$$\begin{array}{}\text{(F1)}& v\left(x\right)=C_{v}H\left(x{)}^m\mathcal{H}\right(x-x_{*}),\mathrm{0}\le x\le x_{\mathrm{GL}}.\end{array}$$

To satisfy the jump condition in (E8) and (E9), the constant C_v is

$$\begin{array}{}\text{(F2)}& C_v=\left\{\begin{array}{rl}{\displaystyle \frac{a{x}_{*}}{\mathit{\rho }g{H}_{*}^{m+\mathrm{3}}}},& F(u,h)=u\mathit{\delta }(x-x_{*}),\\ -\frac{\mathrm{1}}{\mathit{\rho }g{H}_{*}^{m+\mathrm{1}}},& F(u,h)=h\mathit{\delta }(x-x_{*}).\end{array}\right.\end{array}$$

Combine (F1) with the relation ${\mathit{\psi }}_x=(F_{\mathrm{u}}-Cm{u}^{m-\mathrm{1}}v)/H$ and integrate from x to x_GL to obtain

$$\begin{array}{}\text{(F3)}& \mathit{\psi }\left(x\right)=C_v{a}^{m-\mathrm{1}}C\left(x_{\mathrm{GL}}^m-x^m\right),x_{*}<x\le x_{\mathrm{GL}}.\end{array}$$

With the jump condition in (E8) and (E9), ψ(x) at $\mathrm{0}\le x<x_{*}$ is

$$\begin{array}{}\text{(F4)}& \mathit{\psi }\left(x\right)=\left\{\begin{array}{l}-\frac{\mathrm{1}}{H_{*}}+\frac{C{a}^m{x}_{*}}{\mathit{\rho }g{H}_{*}^{m+\mathrm{3}}}\left(x_{\mathrm{GL}}^m-x_{*}^m\right),\\ F(u,h)=u\mathit{\delta }(x-x_{*}),\\ -\frac{C{a}^{m-\mathrm{1}}}{\mathit{\rho }g{H}_{*}^{m+\mathrm{1}}}(x_{\mathrm{GL}}^m-x_{*}^m),\\ F(u,h)=h\mathit{\delta }(x-x_{*}).\end{array}\right.\end{array}$$

The weight for δC in the functional δℒ in (38) is non-zero for $x_{*}<x\le x_{\mathrm{GL}}$:

$$\begin{array}{}\text{(F5)}& -v{u}^m=-C_v(ax{)}^m.\end{array}$$

Use (F1) and (40) in (38) to determine the weight for δb in δℒ:

$$\begin{array}{}\text{(F6)}& \begin{array}{rl}& {\mathit{\psi }}_{x}u+v_x\mathit{\eta }{u}_x+v\mathit{\rho }g{h}_x=\mathit{\rho }g(Hv{)}_x+F_h\\ & =C_v\mathit{\rho }g{H}^m\left[(m+\mathrm{1})H_x\mathcal{H}(x-x_{*})+H\mathit{\delta }(x-x_{*})\right]\\ & +F_h.\end{array}\end{array}$$