Home Contact Us Site Map  
 
       
    next up previous contents
Next: 6.7 Packages Related to Up: 6.6 Sea Ice Packages Previous: 6.6.3 SHELFICE Package   Contents

Subsections


6.6.4 STREAMICE Package

Authors: Daniel Goldberg

6.6.4.1 Introduction

Package ``streamice'' provides a dynamic land ice model for MITgcm.

6.6.4.2 Equations Solved

As of now, the model tracks only 3 variables: $ x$ -velocity ($ u$ ), $ y$ -velocity ($ v$ ), and thickness ($ h$ ). There is also a variable that tracks coverage of fractional cells, discussed...

By default the model solves the Shallow Shelf approximation (SSA) for velocity. The SSA is appropriate for floating ice (ice shelf) or ice flowing over a low-friction bed (e.g. MacAyeal, 1989). The SSA consists of the $ x$ -momentum balance:

$\displaystyle \partial_x(h\nu(4\dot{\varepsilon}_{xx}+2\dot{\varepsilon}_{yy})) + \partial_y(2h\nu\dot{\varepsilon}_{xy}) - \tau_{bx} = \rho g h s_x$ (6.88)

the $ y$ -momentum balance:

$\displaystyle \partial_x(2h\nu\dot{\varepsilon}_{xy}) + \partial_y(h\nu(4\dot{\varepsilon}_{yy}+2\dot{\varepsilon}_{xx})) - \tau_{by} = \rho g h s_y.$ (6.89)

From the velocity field, thickness evolves according to the continuity equation:

$\displaystyle h_t + \nabla\cdot(h\vec{u}) = -\dot{b},$ (6.90)

Where $ \dot{b}$ is a basal mass balance (e.g. melting due to contact with the ocean), positive where there is melting. Surface mass balance is not considered, since pkg/streamice is intended for localized areas (tens to hundreds of kilometers) over which the integrated effect of surface mass balance is generally small. Where ice is grounded, surface elevation is given by

$\displaystyle s = R + h,$ (6.91)

where $ R(x,y)$ is the bathymetry, and the basal elevation $ b$ is equal to $ R$ . If ice is floating, then the assumption of hydrostasy and constant density gives

$\displaystyle s = (1-\frac{\rho}{\rho_w} h,$ (6.92)

where $ \rho_w$ is a representative ocean density, and $ b=-(\rho/\rho_w)h$ . Again by hydrostasy, floation is assumed wherever

$\displaystyle h \leq -\frac{\rho_w}{\rho}R$ (6.93)

is satisfied. Floatation criteria is stored in float_frac_streamice, equal to 1 where ice is at floatation.

The strain rates $ \varepsilon_{ij}$ are generalized to the case of orthogonal curvilinear coordinates, to include the "metric" terms that arise when casting the equations of motion on a sphere or projection on to a sphere (see pkg/SEAICE, 6.6.2.4.8 of the MITgcm documentation). Thus

$\displaystyle \dot{\varepsilon}_{xx} =$ $\displaystyle u_x + k_1 v,$    
$\displaystyle \dot{\varepsilon}_{yy} =$ $\displaystyle v_y + k_1 u,$    
$\displaystyle \dot{\varepsilon}_{xy} =$ $\displaystyle \frac{1}{2}(u_y+v_x) + k_1 u + k_2 v.$    

$ \nu$ has the form arising from Glen's law

$\displaystyle \nu = \frac{1}{2}A^{-\frac{1}{n}}\left(\dot{\varepsilon}_{xx}^2+\... ...y}+\dot{\varepsilon}_{xy}^2+\dot{ \varepsilon}_{min}^2\right)^{\frac{1-n}{2n}},$ (6.94)

though the form is slightly different if a hybrid formulation is used. Whether $ \tau_b$ is nonzero depends on whether the floatation condition is satisfied. Currently this is determined simply on an instantaneous cell-by-cell basis (unless subgrid interpolation is used), as is the surface elevation $ s$ , but possibly this should be rethought if the effects of tides are to be considered. $ \vec{\tau}_b$ has the form

$\displaystyle \vec{\tau}_b = C (\vert\vec{u}\vert^2+u_{min}^2)^{\frac{m-1}{2}}\vec{u}.$ (6.95)

Again, the form is slightly different if a hybrid formulation is to be used. The scalar term multiplying $ \vec{u}$ is referred to as $ \beta$ below.

The momentum equations are solved together with appropriate boundary conditions, discussed below. In the case of a calving front boundary condition (CFBC), the boundary condition has the following form:


Here $ \vec{n}$ is the normal to the boundary, and $ R(x,y)$ is the bathymetry.

6.6.4.2.1 Hybrid SIA-SSA stress balance

The SSA does not take vertical shear stress or strain rates (e.g., $ \sigma_{xz}$ , $ \partial u/\partial z$ ) into account. Although there are other terms in the stress tensor, studies have found that in all but a few cases, vertical shear and longitudinal stresses (represented by the SSA) are sufficient to represent glaciological flow. streamice can allow for representation of vertical shear, although the approximation is made that longitudinal stresses are depth-independent. The stress balance is referred to as "hybrid" because it is a joining of the SSA and the Shallow Ice Approximation (SIA), which only accounts only for vertical shear. Such hybrid formulations have been shown to be valid over a larger range of conditions than SSA (Schoof and Hindmarsh 2010, Goldberg 2011).

In the hybrid formulation, $ \overline{u}$ and $ \overline{v}$ , the depth-averaged $ x-$ and $ y-$ velocities, replace $ u$ and $ v$ in (6.88), (6.89), and (6.90), and gradients such as $ u_x$ are replaced by $ (\overline{u})_x$ . Viscosity becomes

$\displaystyle \nu = \frac{1}{2}A^{-\frac{1}{n}}\left(\dot{\varepsilon}_{xx}^2+\... ...1 }{4}u_z^2+\frac{1}{4}v_z^2+\dot{\varepsilon}_{min}^2\right)^{\frac{1-n}{2n}}.$ (6.98)

In the formulation for $ \tau_b$ , $ u_b$ , the horizontal velocity at $ u_b$ is used instead. The details are given in Goldberg (2011).

6.6.4.3 Numerical Scheme

\resizebox{5.0in}{3.0in}{\includegraphics{s_phys_pkgs/figs/stencil.eps}}

\resizebox{5.0in}{3.0in}{\includegraphics{s_phys_pkgs/figs/mask_cover.eps}}

6.6.4.3.1 Stress/momentum equations

The stress balance is solved for velocities using a very straightforward, structured-grid finite element method. Generally a finite element method is used to deal with irregularly shaped domains, or to deal with complicated boundary conditions; in this case it is the latter, as explained below. (NOTE: this is not meant as a finite element tutorial, and various mathematical objects are defined nonrigorously. It is simply meant as an overview of the numerical scheme used.) In this method, the numerical velocities $ u$ and $ v$ (or $ \overline{u}$ , $ \overline{v}$ if a hybrid formulation) have a bilinear shape within each cell. For instance, referring to Fig. 6.6.4.3, at a point ($ x,y$ ) within cell ($ i,j$ ), and assuming a rectangular mesh, the $ x-$ velocity $ u$ would be given by

$\displaystyle u(x,y) =$ $\displaystyle (1-\zeta_1)\times(1-\zeta_2)\times u_{i,j}\ +$ (6.99)
  $\displaystyle (\zeta_1)\times(1-\zeta_2)\times u_{i+1,j}\ +$ (6.100)
  $\displaystyle (1-\zeta_1)\times(\zeta_2)\times u_{i,j+1}\ +$ (6.101)
  $\displaystyle (\zeta_1)\times(\zeta_2)\times u_{i+1,j+1},$ (6.102)

where

$\displaystyle \zeta_1 = x-x_i, \ \ \zeta_2 = y-y_j$ (6.103)

and e.g. $ u_{i,j+1}$ is the nodal value of $ u$ at the top right corner of the cell. In finite element terms, the functions multiplying a nodal value of $ u$ are the shape functions $ \phi $ corresponding to that node, e.g. $ \phi_{ij}=\zeta_1\zeta_2$ in cell ($ i,j$ ). Note that $ h_{ij}$ is defined at the center of each cell. In the code velocities are U_STREAMICE and V_STREAMICE, and thickness is H_STREAMICE.

If we take a vector-valued function $ \Phi=(\phi,\psi)$ which is in the ``solution space'' of ($ u,v$ ) (i.e. $ \phi $ and $ \psi$ are bilinear in each cell as defined above, and are zero where velocities are imposed), take the dot product with by (6.88,6.89), and integrate by parts over the domain $ D$ , we get the weak form of the momentum equation:

$\displaystyle -$ $\displaystyle \int_D \left(h\nu\phi_x (4\dot{\varepsilon}_{xx}+2\dot{\varepsilo... ...xy} + h\nu\psi_y (4\dot{\varepsilon}_{yy}+2\dot{\varepsilon}_{xx}) \right) dA =$    
  $\displaystyle \int_D \Phi\vec{\tau}_d - \int_{\Gamma} \frac{1}{2}g \left(\rho h^2 - \rho_w R^2\right) \Phi\cdot\vec{n} dS,$ (6.104)

where $ \Gamma$ is the part of the domain boundary where a calving front stress condition is imposed, and $ \vec{\tau}_d$ is the driving stress $ \rho g h \nabla s$ . Note that the boundary integral does not involve viscosity, or any velocity-dependent terms. This is the advantage of using a finite-element formulation: viscosity need not be represented at boundaries, and so no complicated one-sided derivative expressions need to be used.

Let $ (U,V)$ be the vector of all nodal values $ u_{i,j},v_{i,j}$ that determine $ (u,v)$ . If we assume that $ \nu$ and $ \beta$ are independent of $ (U,V)$ , then for any $ (\phi,\psi)$ , the left hand side of (6.104) is a linear functional of $ (U,V)$ , i.e. a linear operator that returns a scalar. We can choose $ \phi $ and $ \psi$ so that they are nonzero only at a single node; and we can evaluate (6.104) for each such $ \phi_{ij}$ and $ \psi_{ij}$ , giving a linearly independent set of equations for $ U$ and $ V$ . The left hand side of (6.104) for each $ \phi_{ij}$ and $ \psi_{ij}$ is evaluated in the subroutine STREAMICE_CG_ACTION() in the source file streamice_cg_functions.f. The right hand side is evaluated in STREAMICE_DRIVING_STRESS(). The latter is not a straightforward calculation, since the thickness $ h$ and the surface $ s$ are not expressed by piecewise bilinear functions, and is discussed below.

STREAMICE_CG_ACTION() represents the action of a matrix on the vector represented by $ (U,V)$ and STREAMICE_DRIVING_STRESS() gives the right hand side of a linear system of equations. This information is enough to solve the linear system, and this is done with the method of conjugate gradients, implemented (with simple Jacobi preconditioning) in STREAMICE_CG_SOLVE().

The full nonlinear system is solved in STREAMICE_VEL_SOLVE(). This subroutine evaluates driving stress and then calls a loop in which a sequence of linear systems is solved. After each linear solve, the viscosity $ \nu$ (the array visc_streamice) and basal traction $ \beta$ (tau_beta_eff_streamice) is updated from the new iterate of $ u$ and $ v$ . The iteration continues until a desired tolerance is reached or the maximum number of iterations is reached.

Rather than calling STREAMICE_CG_ACTION() each iteration of the conjugate gradient, the sparse matrix can be constructed explicitly. This is done if STREAMICE_CONSTRUCT_MATRIX is #defined in STREAMICE_OPTIONS.h. This speeds up the solution because STREAMICE_CG_ACTION() contains a number of nested loops related to quadrature points and interactions between nodes of a cell.

The driving stress $ \tau_d=\rho g \nabla s$ is not as straightforward to deal with in the weak formulation, as $ h$ and $ s$ are considered constant in a cell. Furthermore, in the ice shelf the driving stress can be written as

$\displaystyle \vec{\tau}_d = \frac{1}{2}\rho g (1-\frac{\rho}{\rho_w})\nabla(h^2).$ (6.105)

Among other things, this (along with equations 6.96 and 6.97) implies that for an unconfined shelf (one for which the only boundaries are a calving front and a grounding line), the stress state along the grounding line depends only on the location at depth of the grounding line. The numerical representation of velocities must respect this. Thus the driving stress contribution in the $ x-$ momentum equation at node $ (i,j)$ is given as follows.

$\displaystyle \tau_{d,x}(i,j) = \begin{cases}\frac{1}{4}\rho g (\gamma_{ij} h_{... ...ta x_{f,ij}} g (\rho h_{ij}^2 - \rho_w b_{ij}^2) & \mbox{case III}. \end{cases}$ (6.106)

where

$\displaystyle \gamma_{ij} = \sqrt{\frac{rA_{ij}}{\Delta x_{f,ij}}},$ (6.107)

$ rA_{ij}$ is the area of cell ($ i,j$ ) and $ \Delta x_{f,ij}$ is the width of the cell at its $ y-$ midpoint. Cases I, II, and III refer to the situations when (I) both cells ($ i,j$ ) and ($ i-1,j$ ) are in the domain, (II) only cell ($ i-1,j$ ) is in the domain, and (III) only cell ($ i,j$ ) is in the domain. (The fact that $ \tau_{d,x}(i,j)$ is being calculated means that the implied boundary is a calving stress boundary, as explained below.) The above expression gives the contribution from the cells above node ($ i,j$ ) (in the $ y-$ direction). A similar expression gives the contribution from cells ($ i-1,j-1$ ) and ($ i,j-1$ ), and corresponding expressions approximate driving stress in the $ y-$ momentum equation. In the ice shelf, the case I expressions give approximations for (6.105) which are discretely integrable; in grounded ice the expression need not be integrable, so any errors associated with varying grid metrics can be subsumed into the discretization error.

6.6.4.3.2 Orthogonal curvilinear coordinates

If the grid is not rectangular, but given in orthogonal curvilinear coordinates, the evaluation of (6.104) is not exact, since the gradient of the basis functions are approximated, as we do not have complete knowledge of the coordinate system, only the grid metrics. Still, we consider nodal basis functions that are equal to 1 at a single node, and zero elsewhere in a cell.

A cell ($ i,j$ ) has separation $ \Delta x_{i,j}$ at the bottom and $ \Delta x_{i,j+1}$ at the top (referred to as dxG(i,j) and dxG(i,j+1) in the code). So, for instance, for a cell given by $ [\xi_1,\xi_2]\times[\eta_1,\eta_2]$ , the $ \xi$ -derivative of the basis function centered at ( $ \xi_1,\eta_1$ ) is approximated by

$\displaystyle \left(\frac{1-q}{\Delta x_{i,j}}\right) + \left(\frac{q}{\Delta x_{i,j+1}}\right)$ (6.108)

where $ q$ is the $ y$ -coordinate of the quadrature point being used. (The code currently uses 2-point gauss quadrature to evaluate the integrals in (6.104); with a rectangular grid and cellwise-constant $ \beta$ and $ \nu$ , this is exact.)

Basis function derivatives and quadrature weights are stored in Dphi and grid_jacq_streamice. Both are initialized in STREAMICE_INIT_PHI, called only once.

6.6.4.3.3 Boundary conditions and masks

The computational domain (which may be smaller than the array/grid as defined by SIZE.h and GRID.h) is determined by a number of mask arrays within the STREAMICE package. They are

  • $ hmask$ (STREAMICE_hmask): equal to 1 (ice-covered), 0 (open ocean), 2 (partly-covered), or -1 (out of domain)
  • $ umask$ (STREAMICE_umask): equal to 1 (an ``active'' velocity node), 3 (a Dirichlet node), or 0 (zero velocity)
  • $ vmask$ (STREAMICE_vmask): similar to umask
  • $ ufacemaskbdry$ (STREAMICE_ufacemaskbdry): equal to -1 (interior face), 0 (no-slip), 1 (no-stress), 2 (calving stress front), 3 (Dirichlet), or 4 (flux input boundary); when 4, then u_flux_bdry_SI must be initialized, through binary or parameter file
  • $ vfacemaskbdry$ (STREAMICE_vfacemaskbdry): similar to ufacemaskbdry
$ hmask$ is defined at cell centers, like $ h$ . $ umask$ and $ vmask$ are defined at cell nodes, like velocities. $ ufacemaskbdry$ and $ vfacemaskbdry$ are defined at cell faces, like velocities in a $ C$ -grid - but unless STREAMICE_GEOM_FILE_SETUP is #defined in STREAMICE_OPTIONS.h, the values are only relevant at the boundaries of the grid.

Essentially $ hmask$ controls the domain; it is specified at initialization and does not change unless calving front advance is allowed. $ ufacemaskbdry$ and $ vfacemaskbdry$ are defined at initialization as well. Masks are either initialized through binary files or through the text file data.streamice (see the streamice tutorial)

The values of $ umask$ and $ vmask$ determine which nodal values of $ u$ and $ v$ are involved in the solve for velocities. Before each call to STREAMICE_VEL_SOLVE(), $ umask$ and $ vmask$ are re-initialized. Values are only relevant if they border a cell where $ hmask=1$ . Furthermore, if a cell face is a ``no-slip'' face, both $ umask$ and $ vmask$ are set to zero at the face's nodal endpoints. If the face is a Dirichlet boundary, both nodal endpoints are set to 3. If $ ufacemaskbdry=1$ on the $ x-$ boundary of a cell, i.e. it is a no-stress boundary, then $ vmask$ is set to 1 at both endpoints while $ umask$ is set to zero (i.e. no normal flow). If the face is a flux boundary, velocities are set to zero (see ???). For a calving stress front, $ umask$ and $ vmask$ are 1: the nodes represent active degrees of freedom.

With $ umask$ and $ vmask$ appropriately initialized, STREAMICE_VEL_SOLVE can proceed rather generally. Contributions to (6.104) are only evaluated if $ hmask=1$ in a given cell, and a given nodal basis function is only considered if $ umask=1$ or $ vmask=1$ at that node.

6.6.4.3.4 Thickness evolution

(6.90) is solved in the subroutine STREAMICE_ADVECT_THICKNESS, similarly to the advection routines in MITgcm. Mass fluxes are evaluated, first in the $ x$ -direction. This is done in the generic subroutine STREAMICE_ADV_FLUX_FL_X. Flux velocity in the $ x-$ direction at face ($ i,j$ ) are generated by averaging $ u_{i,j}$ and $ u_{i,j+1}$ . Assuming the flux velocity is positive, if $ hmask_{i-2,j},\ mask{i-1,j}$ and $ hmask_{i,j}$ are equal to 1, then flux thickness, i.e. the interpolation of $ h$ at face ($ i,j$ ), is through a minmod limiter as in the generic_advdiff package. If these values are not available, then a zero-order upwind flux is used. The exception is when STREAMICE_ufacemask(i,j) is equal to 4; then u_flux_bdry_SI(i,j) is used for the flux. Fluxes are then differenced to update $ h$ in cells that are active ($ hmask=1$ ); a similar procedure follows for the $ y-$ direction.

6.6.4.3.5 Ice front advance

Flux out of the domain across Dirichlet boundaries is ignored, and flux at no-flow or no-stress boundaries is zero. However, fluxes that cross calving boundaries are stored, and if STREAMICE_move_front=.true., then STREAMICE_ADV_FRONT is called after all cells with $ hmask=1$ have their thickness updated. In this subroutine, cells with $ hmask=0$ or $ hmask=2$ with nonzero fluxes entering their boundaries are processed.

The algorithm is based on Albrecht (2011). In this scheme, if cell ($ i,j$ ) fits the criteria, a reference thickness $ h_{ref}(i,j)$ is found, defined as an average over the thickness of all neighboring cells with $ hmask=1$ that are contributing mass to ($ i,j$ ). The total mass input over the time step to ($ i,j$ ) is calculated as $ V_{in}(ij)$ . This is added to the volume of ice already in the cell, which is nonzero only if $ hmask_{ij}=2$ at the beginning of the time step, and is equal to

$\displaystyle V_{cell}(i,j) = area_{ij} h_{ij}.$ (6.109)

(If $ hmask_{ij}$ is not equal to 1, then $ h_{ij}$ has not yet been updated in the current time step.) $ area_{ij}$ is stored in the array area_shelf_streamice. The total volume, $ V_{tot}(i,j)$ , is then equal to $ V_{cell}(i,j)+V_{in}(i,j)$ . Next an effective thickness

$\displaystyle h_{eff}(i,j) = \frac{V_{ij}}{rA_{ij}}$ (6.110)

is found. If $ h_{eff}(ij)<h_{ref}(i,j)$ (but nonzero) then
  • $ hmask_{ij}$ is set to 2,
  • $ h_{ij}$ is set to $ h_{ref}(ij)$ ,
  • $ area_{ij}$ is set to $ \frac{V_{tot}(i,j)}{h_{ref}(i,j)}$ .
If, on the other hand, $ h_{eff}(ij)\geq h_{ref}(i,j)$ ,
  • $ hmask_{ij}$ is set to 1,
  • $ h_{ij}$ is set to $ h_{ref}(ij)$ ,
  • $ area_{ij}$ is set to $ rA_{ij}$ ,
  • Excess flux $ q_{ex}(i,j) = V_{tot}(i,j) - h_{eff}(i,j)rA_{ij}$ is found.
If $ q_{ex}(i,j)$ is zero, nothing more needs to happen. If it is positive, then
  • adjacent cells that are not completely covered ($ hmask=0$ or 1) are searched for.
  • if there are no suitable cells, $ q_{ex}(i,j)$ is put back into cell ($ i,j$ ) ($ h_{ij}$ is set to $ h_{ref}(i,j)+\frac{q_{ex}(i,j)}{rA_{ij}}$ ).
  • if there are suitable cells, $ q_{ex}$ is divided uniformly among the corresponding faces of the receiving cells, and the algorithm begins again.
The last step is a recursive one; in theory the recursion could continue indefinitely, but in practice at most one recursive call is done. The decision to spread the excess flux uniformly is not the most physically motivated choice, but with appropriate time steps the excess flux value is expected to be small, and is more a matter of mass/volume conservation. This algorithm, while complicated, is necessary for realistic front advance. If newly-covered cells were given full coverage instead of partial coverage and volume were conserved, the front would quickly diffuse.

If calve_to_mask=.true., this sets a limit to how far the front can advance, even if advance is allowed. The front will not advance into cells where the array calve_mask is not equal to 1. This mask must be set through a binary input file.

No calving parameterisation is implemented in STREAMICE. However, front advancement is a precursor for such a development to be added.

6.6.4.3.6 Surface/basal mass balance

After updating thickness in fully- and partially-covered cells, surface and basal mass balances are applied to nonempty cells. Currently surface mass balance is a uniform value, set through the parameter streamice_adot_uniform. Basal melt rates are stored in the array BDOT_streamice, but are multiplied through by float_frac_streamice before being applied.

next up previous contents
Next: 6.7 Packages Related to Up: 6.6 Sea Ice Packages Previous: 6.6.3 SHELFICE Package   Contents
mitgcm-support@mitgcm.org
Copyright © 2006 Massachusetts Institute of Technology Last update 2018-01-23