Notes about the flow-line SSA¶

Using the same notation as in the rest of the manual, the flow-line case the shallow shelf approximation reads

(64)¶\[\left(4 \nu H u' \right)' = \rho_i\, g\, H\, h'.\]

Here \(u\) is the ice velocity, \(\nu = \nu(u')\) is the ice viscosity, \(H\) is the ice thickness, and \(h\) is the ice surface elevation.

Let \(\alpha = (1 - \rhoi / \rhow)\) and note that \(h = \alpha\, H\) whenever ice is floating, so

(65)¶\[\left(4 \bar \nu H u' \right)' = \frac12 \alpha\, \rho_i\, g\, (H^2)'.\]

We can easily integrate this, getting

(66)¶\[4 \bar \nu H u' = \frac12 \alpha\, \rho_i\, g\, H^2.\]

Note

This is a non-linear first-order ODE for the ice velocity \(u\).

An observation¶

Consider a boundary-value problem for the velocity \(u\) of a floating ice shelf on an interval \([0, L]\):

(67)¶\[ \begin{align}\begin{aligned}\left(4 \nu H u' \right)' &= \rho_i\, g\, H\, h',\\u(0) &= u_{0},\\\left. \left( 4\nu H u' \right) \right|_{x=L} &= \tau_{\text{stat}}(L),\\\tau_{\text{stat}}(x) &= \frac12 \rhoi\, g\, \alpha\, H(x)^2.\end{aligned}\end{align} \]

Now consider a similar BVP on \([0, a]\) for some positive \(a < L\):

(68)¶\[ \begin{align}\begin{aligned}\left(4 \nu H v' \right)' &= \rho_i\, g\, H\, h',\\v(0) &= u_{0},\\\left. \left( 4\nu(v') H v' \right) \right|_{x=a} &= \tau_{\text{stat}}(a).\end{aligned}\end{align} \]

Because (67) is a first-order ODE, the solutions \(u\) of (67) and \(v\) of (68) coincide on \([0, a]\), i.e. \(u(x) = v(x)\) for all \(x \in [0, a]\).

Note

This implies that the velocity \(u(x)\) of a floating flow-line ice shelf modeled by (67) is not sensitive to ice geometry perturbations at locations downstream from \(x\).

Discrete analog of this property¶

Let \(N > 2\) be an integer and define the \(N\)-point grid \(x_{j} = 0 + j\dx\), \(\dx = L/(N-1)\).

Discretizing (67) on this grid results in a non-linear system with \(N-1\) unknowns (call this system \(S_{L}\)).

Let \(a = L - \dx\) and write down a system of equations by discretizing (68) on the grid consisting of the first \(N-1\) of \(x_{j}\). This system (call it \(S_{a}\)) has \(N-2\) unknowns.

The property we would like our discretization to have is this:

If \(u_{1},\dots,u_{N-1}\) is the solution of \(S_{L}\) and \(v_{1},\dots,v_{N-2}\) solves \(S_{a}\), then \(u_{k} = v_{k}\) for all \(k = 1,\dots,N-2\).

Note that the first \(N-3\) equations in \(S_{L}\) and \(S_{a}\) are the same.

Discretization¶

Let \([X]_{j}\) be an approximation of \(X\) at a grid location \(j\).

The standard approach within the domain is to use centered finite differences and linear interpolation to approximate staggered-grid values, i.e. averaging values at immediate regular grid point neighbors:

(69)¶\[ \begin{align}\begin{aligned}{}[f]_{i + \frac12} &= \frac12 (f_i + f_{i+1}),\\{}[f']_{i} &= \frac1{\dx}([f]_{i+\frac12} + [f]_{i-\frac12}).\end{aligned}\end{align} \]

Interior¶

(70)¶\[\frac1{\dx} \left( [4\nu(u')H u']_{i+\frac12} - [4\nu(u')H u']_{i-\frac12} \right) = [\rhoi\, g\, \alpha\, H\, H']_{i}.\]

Ice front¶

Let \(n\) be the last grid point with non-zero ice thickness, i.e. assume that \(H_{k} = 0\) for \(k > n\).

The implementation of the stress boundary condition at the ice front amounts to adding one more equation (see (63)):

(71)¶\[{}[4\nu(u') H u']_{n+\frac12} = [\tau_{\text{stat}}]_{n+\frac12}.\]

We can then combine (71) with (70) (with \(i\) replaced by \(n\)) to get the discretization at the ice front:

(72)¶\[ \begin{align}\begin{aligned}\frac1{\dx}\left( [\tau_{\text{stat}}]_{n+\frac12} - [4\nu(u')H u']_{n-\frac12} \right) &= \rhoi\, g\, \alpha\, [H\, H']_{n},\\- [4\nu(u')H u']_{n-\frac12} &= \dx\rhoi\, g\, \alpha\, [H\, H']_{n} - [\tau_{\text{stat}}]_{n+\frac12}.\end{aligned}\end{align} \]

Choosing FD approximations¶

Assuming that the ice front is at \(n\)

If we assume that the ice front is at \(n\), the last equation in the system looks like (72):

(73)¶\[- [4\nu(u')H u']_{n-\frac12} = \dx\rhoi\, g\, \alpha\, [H\, H']_{n} - [\tau_{\text{stat}}]_{n+\frac12}.\]

Assuming that the ice front is at \(n+1\)

If we assume that the ice front is at \(n+1\), the \(n\)-th equation looks like a “generic” interior equation (70) and we have one more equation ((73)) shifted by \(1\):

(74)¶\[ \begin{align}\begin{aligned}\frac1{\dx} \left( [4\nu(u')H u']_{n+\frac12} - [4\nu(u')H u']_{n-\frac12} \right) &= [\rhoi\, g\, \alpha\, H\, H']_{n}.\\- [4\nu(u')H u']_{(n+1)-\frac12} &= \dx\rhoi\, g\, \alpha\, [H\, H']_{n+1} - [\tau_{\text{stat}}]_{n+1+\frac12}.\end{aligned}\end{align} \]

Note that the second equation in (74) is the same as (73), but with the index shifted by \(1\). Both correspond to locations at the ice front.

The goal

We want to choose FD approximations \([\tau_{\text{stat}}]_{*}\) and \([H\, H']_{*}\) in a way that would make it possible to obtain (73) by transforming equations (74).

We propose using constant extrapolation to approximate \(H_{n+\frac12}\)

(75)¶\[\begin{split}{}[H]_{i+\frac12} &= \begin{cases} H_{i}, &i \text{ is at the ice front}\\ \frac12(H_{i} + H_{i+1}), &\text{otherwise}. \end{cases}\end{split}\]

This gives us the following approximation of derivatives:

\[\begin{split}{}[H']_{i} &= \begin{cases} \frac1{\dx}(H_{i} - [H]_{i-\frac12}), &i \text{ is at the ice front}\\ \frac1{\dx}([H]_{i+\frac12} - [H]_{i-\frac12}), &\text{otherwise}. \end{cases}\end{split}\]

After substituting (69) this becomes

(76)¶\[\begin{split}{}[H']_{i} &= \begin{cases} \frac1{2\dx}(H_{i} - H_{i-1}), &i \text{ is at the ice front}\\ \frac1{2\dx}(H_{i+1} - H_{i-1}), &\text{otherwise}. \end{cases}\end{split}\]

Note

The ice front case in (76) is the one half of the standard one-sided finite-difference approximation of \(H'\).

Checking if (75) is the right choice

Consider the first equation in (74) and note that it corresponds to the case in which \(n\) is not at the ice front.

\[\frac1{\dx} \left( [4\nu(u')H u']_{n+\frac12} - [4\nu(u')H u']_{n-\frac12} \right) = [\rhoi\, g\, \alpha\, H\, H']_{n},\]

Multiplying by \(\dx\) and moving one of the terms to the right hand side, we get

(77)¶\[ \begin{align}\begin{aligned}- [4\nu(u')H u']_{n-\frac12} &= \dx[\rhoi\, g\, \alpha\, H\, H']_{n} - [4\nu(u')H u']_{n+\frac12}.\\- [4\nu(u')H u']_{n-\frac12} &= \dx\, \rhoi\, g\, \alpha\, H_{n}\frac{H_{n+1} - H_{n-1}}{2\dx} - [4\nu(u')H u']_{n+\frac12}.\\&= \frac12\, \rhoi\, g\, \alpha\, H_{n}\, (H_{n+1} - H_{n-1}) - [4\nu(u')H u']_{n+\frac12}.\end{aligned}\end{align} \]

Now consider the second equation in (74). Note that here \(n+1\) is at the ice front.

(78)¶\[ \begin{align}\begin{aligned}- [4\nu(u')H u']_{(n+1)-\frac12} &= \dx\rhoi\, g\, \alpha\, [H\, H']_{n+1} - [\tau_{\text{stat}}]_{n+1+\frac12}.\\&= \dx\rhoi\, g\, \alpha\, H_{n+1}[H']_{n+1} - \frac12 \rhoi\, g\, \alpha\, [H]_{n+1+\frac12}^{2}.\\&= \dx\rhoi\, g\, \alpha\, H_{n+1}\frac{H_{n+1} - H_{n}}{2\dx} - \frac12 \rhoi\, g\, \alpha\, H_{n+1}^{2}\\&= \frac12\, \rhoi\, g\, \alpha( H_{n+1}(H_{n+1} - H_{n}) - H_{n+1}^{2} )\\&= \frac12\, \rhoi\, g\, \alpha( H_{n+1}^{2} - H_{n+1}H_{n} - H_{n+1}^{2} )\\&= -\frac12\, \rhoi\, g\, \alpha\, H_{n+1}H_{n}.\end{aligned}\end{align} \]

Put together, (77) and (78) read as follows:

(79)¶\[ \begin{align}\begin{aligned}- [4\nu(u')H u']_{n-\frac12} &= \frac12\, \rhoi\, g\, \alpha\, H_{n}\, (H_{n+1} - H_{n-1}) - [4\nu(u')H u']_{n+\frac12},\\- [4\nu(u')H u']_{n+\frac12} &= -\frac12\, \rhoi\, g\, \alpha\, H_{n+1}H_{n}.\end{aligned}\end{align} \]

Substituting the second equation into the first produces

(80)¶\[ \begin{align}\begin{aligned}- [4\nu(u')H u']_{n-\frac12} &= \frac12\, \rhoi\, g\, \alpha\, H_{n}\, (H_{n+1} - H_{n-1}) - \frac12\, \rhoi\, g\, \alpha\, H_{n+1}H_{n}\\&= \frac12\, \rhoi\, g\, \alpha\, (H_{n}\, (H_{n+1} - H_{n-1}) - H_{n+1}H_{n})\\&= -\frac12\, \rhoi\, g\, \alpha\, H_{n}\, H_{n-1}\end{aligned}\end{align} \]

Compare (80) to (78). Note that they are the same, except for the index shift. In other words, (80) is the same as (73), as desired.

This confirms that finite difference approximations (75) and (76) result in a discretization with the property we seek:

modeled ice velocity at a given location \(x\) along a flow-line ice shelf is not sensitive to geometry perturbations downstream from \(x\).