Governing Equations

Preliminaries

WAVI.jl uses a Cartesian co-ordinate system $\mathbf{x}=(x,y,z)$, with $z$ positive upwards; the corresponding velocity components are $\mathbf{u}=(u,v,w)$. We use bar notation to denote depth averaged quantities, for example:

$$\begin{array}{}\text{(1)}& \overline{f}=\frac{1}{h}{\int }_{z=b(x,y)}^{z=s(x,y,t)}f\mathrm{d}z\end{array}$$

is the depth average of the quantity $f$. Here, $t$ is the denotes time, $b(x,y)$ is the (known) bed elevation (measured positive upwards), and $s(x,y,t)$ is the surface elevation.

We assume that the ice is in hydrostratic equilibrium, so that regions are with $h<-({\rho }_i/{\rho }_w)b$ are floating, and regions with $h\ge -({\rho }_i/{\rho }_w)b$ are grounded, where ${\rho }_i$ and ${\rho }_w$ the ice and ocean density, respectively. Where the ice is grounded, we have $s=h+b$, while where the ice is floating, the hydrostratic assumption enforces $s=(1-{\rho }_i/{\rho }_w)h$.

WAVI.jl solves equations describing conservation of momentum and conservation of mass for $\overline{\mathbf{u}}(x,y,t)=(\overline{u}(x,y,t),\overline{v}(x,y,t))$, the depth averaged velocity components in the $(x,y)$ directions, respectively, and the ice thickness $h(x,y,t)$

Conservation of Momentum

Conservation of momentum requires that the $\overline{u}$ and $\overline{v}$ satisfy ([1]):

$$\begin{array}{}\text{(2)}& \frac{\partial }{\partial x}(4\overline{\eta }h\frac{\partial \overline{u}}{\partial x}+2\overline{\eta }h\frac{\partial \overline{v}}{\partial y}))+\frac{\partial }{\partial y}(\overline{\eta }h\frac{\partial \overline{v}}{\partial x}+\overline{\eta }h\frac{\partial \overline{u}}{\partial y})-{\tau }_{b,x}& ={\rho }_{i}gh\frac{\partial s}{\partial x},\text{(3)}& \frac{\partial }{\partial y}(4\overline{\eta }h\frac{\partial \overline{v}}{\partial y}+2\overline{\eta }h\frac{\partial \overline{u}}{\partial x}))+\frac{\partial }{\partial x}(\overline{\eta }h\frac{\partial \overline{u}}{\partial y}+\overline{\eta }h\frac{\partial \overline{v}}{\partial x})-{\tau }_{b,y}& ={\rho }_{i}gh\frac{\partial s}{\partial y},\end{array}$$

where ${\rho }_i$ is the ice density, $g$ is the gravitational acceleration, ${\tau }_b=({\tau }_{b,x},{\tau }_{b,y})$ is the basal drag in the $(x,y)$ directions, and $\eta$ is the ice viscosity, defined implicity in terms of the velocity components (the strain components are themselves functions of $\eta$, see below):

$$\begin{array}{}\text{(4)}& \eta =\frac{B}{2}{[{\left(\frac{\partial \overline{u}}{\partial x}\right)}^2+{\left(\frac{\partial \overline{v}}{\partial y}\right)}^2+\frac{\partial \overline{u}}{\partial x}\frac{\partial \overline{v}}{\partial y}+\frac{1}{4}{(\frac{\partial \overline{u}}{\partial y}+\frac{\partial \overline{u}}{\partial x})}^2+\frac{1}{4}{\left(\frac{\partial \overline{u}}{\partial z}\right)}^2+{\left(\frac{\partial \overline{v}}{\partial z}\right)}^2+{ϵ}^2]}^{\frac{1-n}{2n}}.\end{array}$$

Here $n$ is the exponent in a nonlinear Glen flow law, $ϵ$ is a regularization parameter that prevents the viscosity becoming unbounded at small strain rates (for small strain rates, $\eta$ is constant, corresponding to a linear rheology), and $B(x,y,z)$ is a temperature-dependent coefficient that determines the stiffness of the ice.

Boundary Conditions

The momentum equations are solved alongside boundary conditions at the lateral boundary of the ice sheet,

$$\begin{array}{}\text{(5)}& -\frac{1}{2}{\rho }_w{h}_w^2{\hat{n}}_x=2\overline{\eta }h(2\frac{\partial \overline{u}}{\partial x}+\frac{\partial \overline{v}}{\partial y}){\hat{n}}_y-\frac{1}{2}{\rho }_{i}g{h}^2{\hat{n}}_x+\overline{\eta }h(\frac{\partial \overline{u}}{\partial y}+\frac{\partial \overline{v}}{\partial x}){\hat{n}}_y,\text{(6)}& -\frac{1}{2}{\rho }_w{h}_w^2{\hat{n}}_y=2\overline{\eta }h(2\frac{\partial \overline{v}}{\partial y}+\frac{\partial \overline{u}}{\partial x}){\hat{n}}_y-\frac{1}{2}{\rho }_{i}g{h}^2{\hat{n}}_y+\overline{\eta }h(\frac{\partial \overline{u}}{\partial y}+\frac{\partial \overline{v}}{\partial x}){\hat{n}}_x,\end{array}$$

which impose continuity of depth-integrated momentum there. In $\text{(5)}$–$\text{(6)}$, $h_w=max(h-s+\zeta ,0)$ is the thickness of ice below the water level, where $\zeta$ is the sea level with respect to $z=0$, and $\hat{\mathbf{n}}=({\hat{n}}_x,{\hat{n}}_y)$ is the normal to the lateral boundary.

In addition, a Robin boundary condition at the bed linearly relates the basal stress ${\tau }_b=({\tau }_{b,x},{\tau }_{b,y})$ to the basal velocity ${\mathbf{u}}_b$ via a multiplicative drag coefficient $\beta$:

$$\begin{array}{}\text{(7)}& {\tau }_b=\beta {\mathbf{u}}_b.\end{array}$$

(A no-stress condition at the surface is also implicit in the derivation of~$\text{(2)}$–$\text{(2)}$.)

Conservation of Mass

For a given depth-averaged velocity $\overline{\mathbf{u}}$, accumulation rate $a(x,y,t)$ (positive for ice gain), and basal melt rate $m(x,y,t)$ (positive for ice loss), conservation of ice mass requires that the ice thickness $h$ satisfies

$$\begin{array}{}\text{(8)}& \frac{\partial h}{\partial t}=a-m-\mathrm{\nabla }.\left(h\mathbf{u}\right).\end{array}$$