8.6.5. STREAMICE Package

Author: Daniel Goldberg

8.6.5.1. Introduction

Package STREAMICE provides a dynamic land ice model for MITgcm. It was created primarily to develop a TAF- and OpenAD-generated ice model adjoint and to provide synchronous ice-ocean coupling through the SHELFICE package. It solves a set of dynamic equations appropriate for floating ice-shelf flow as well as ice-stream and slower ice-sheet flow. It has been tested at the scale of one or several ice streams, but has not been tested at the continental scale.

8.6.5.2. STREAMICE configuration

8.6.5.2.1. Compile-time options

pkg/streamice can be included on at compile time in the packages.conf file by adding a line streamice (see Section 8.1.1).

Parts of the pkg/streamice code can be enabled or disabled at compile time via CPP flags. These options are set in STREAMICE_OPTIONS.h.

Table 8.23 CPP flags used by pkg/streamice.

CPP Flag Name

Default

Description

STREAMICE_CONSTRUCT_MATRIX

#define

explicit construction of matrix for Picard iteration for velocity

STREAMICE_HYBRID_STRESS

#undef

use L1L2 formulation for stress balance (default shallow shelf approx.)

USE_ALT_RLOW

#undef

use package array for rLow rather than model

STREAMICE_GEOM_FILE_SETUP

#undef

use files rather than parameters in STREAMICE_PARM03 to configure boundaries

ALLOW_PETSC

#undef

enable interface to PETSc for velocity solver matrix solve

STREAMICE_COULOMB_SLIDING

#undef

enable basal sliding of the form (8.48)

8.6.5.2.2. Enabling the package

Once it has been compiled, pkg/streamice is switched on/off at run-time by setting useSTREAMICE to .TRUE. in file data.pkg.

8.6.5.2.3. Runtime parmeters: general flags and parameters

Run-time parameters are set in file data.streamice (read in streamice_readparms.F). General pkg/streamice parameters are set under STREAMICE_PARM01 as described in Table 8.24.

Table 8.24 Run-time parameters and default values (defined under STREAMICE_PARM01 namelist)

Parameter

Default

Description

streamice_density

910

the (uniform) density of land ice (kg/m3)

streamice_density_ocean_avg

1024

the (uniform) density of ocean (kg/m3)

n_glen

3

Glen’s Flow Law exponent (non-dim.)

eps_glen_min

1e-12

minimum strain rate in Glen’s Law (\(\varepsilon_0\), yr-1)

eps_u_min

1e-6

minimum speed in nonlinear sliding law (\(u_0\), m/yr)

n_basal_friction

0

exponent in nonlinear sliding law (non-dim.)

streamice_cg_tol

1e-6

tolerance of conjugate gradient of linear solve of Picard iteration for velocity

streamice_lower_cg_tol

TRUE

lower CG tolerance when nonlinear residual decreases by fixed factor

streamice_max_cg_iter

2000

maximum iterations in linear solve

streamice_maxcgiter_cpl

0

as above when coupled with pkg/shelfice

streamice_nonlin_tol

1e-6

tolerance of nonlinear residual for velocity (relative to initial)

streamice_max_nl_iter

100

maximum Picard iterations in solve for velocity

streamice_maxnliter_cpl

0

as above when coupled with pkg/shelfice

streamice_nonlin_tol_fp

1e-6

tolerance of relative change for velocity iteration (relative to magnitude)

streamice_err_norm

0

type of norm evaluated for error (\(p\) in \(p\)-norm; 0 is \(\infty\))

streamice_chkfixedptconvergence

FALSE

terminate velocity iteration based on relative change per iteration

streamice_chkresidconvergence

TRUE

terminate velocity iteration based on residual

streamicethickInit

FILE

method by which to initialize thickness (FILE or PARAM)

streamicethickFile

' '

thickness initialization file, in meters (rather than parameters in STREAMICE_PARM03)

streamice_move_front

FALSE

allow ice shelf front to advance

streamice_calve_to_mask

FALSE

if streamice_move_front TRUE do not allow to advance beyond streamice_calve_mask

streamicecalveMaskFile

' '

file to initialize streamice_calve_mask

streamice_diagnostic_only

FALSE

do not update ice thickness (velocity solve only)

streamice_CFL_factor

0.5

CFL factor which determine maximum time step for thickness sub-cycling

streamice_adjDump

0.0

frequency (s) of writing of adjoint fields to file (TAF only)

streamicebasalTracConfig

UNIFORM

method by which to initialize basal traction (FILE or UNIFORM)

streamicebasalTracFile

' '

basal trac initialization file (see Units of input files for units)

C_basal_fric_const

31.71

uniform basal traction value (see Units of input files for units)

streamiceGlenConstConfig

UNIFORM

method by which to initialize Glen’s constant (FILE or UNIFORM)

streamiceGlenConstFile

' '

Glen’s constant initialization file (see Units of input files for units)

B_glen_isothermal

9.461e-18

uniform Glen’s constant value (see Units of input files for units)

streamiceBdotFile

' '

file to initialize time-indep melt rate (m/yr)

streamiceBdotTimeDepFile

' '

file to initialize time-varying melt rate (m/yr), based on streamice_forcing_period

streamiceTopogFile

' '

topography initialization file (m); requires #define USE_ALT_RLOW

streamiceHmaskFile

' '

streamice_hmask initialization file; requires #define STREAMICE_GEOM_FILE_SETUP

streamiceuFaceBdryFile

' '

streamice`STREAMICE_ufacemask_bdry` initialization file; requires #define STREAMICE_GEOM_FILE_SETUP

streamicevFaceBdryFile

' '

streamice`STREAMICE_vfacemask_bdry`` initialization file; requires #define STREAMICE_GEOM_FILE_SETUP

streamiceuMassFluxFile

' '

mass flux at \(u\)-faces init. file (m2/yr); requires #define STREAMICE_GEOM_FILE_SETUP

streamicevMassFluxFile

' '

mass flux at \(v\)-faces init. file (m2/yr); requires #define STREAMICE_GEOM_FILE_SETUP

streamiceuFluxTimeDepFile

' '

time-depend. mass flux at \(u\)-faces file (m2/yr); requires #define STREAMICE_GEOM_FILE_SETUP

streamicevFluxTimeDepFile

' '

time-depend. mass flux at \(v\)-faces file (m2/yr); requires #define STREAMICE_GEOM_FILE_SETUP

streamiceuNormalStressFile

' '

calving front normal stress parm along \(u\)-faces (non-dim.; see Boundary Stresses)

streamicevNormalStressFile

' '

calving front normal stress parm along \(v\)-faces (non-dim.; see Boundary Stresses)

streamiceuShearStressFile

' '

calving front normal stress parm along \(u\)-faces (non-dim.; see Boundary Stresses)

streamicevShearStressFile

' '

calving front normal stress parm along \(v\)-faces (non-dim.; see Boundary Stresses)

streamiceuNormalTimeDepFile

' '

time-dependent version of streamiceuNormalStressFile

streamicevNormalTimeDepFile

' '

time-dependent version of streamicevNormalStressFile

streamiceuShearTimeDepFile

' '

time-dependent version of streamiceuShearStressFile

streamicevShearTimeDepFile

' '

time-dependent version of streamicevShearStressFile

streamice_adot_uniform

0

time/space uniform surface accumulation rate (m/yr)

streamice_forcing_period

0

file input frequency for streamice time-dependent forcing fields (s)

streamice_smooth_gl_width

0

thickness range parameter in basal traction smoothing across grounding line (m)

streamice_allow_reg_coulomb

FALSE

use regularized Coulomb sliding (8.48). Requires STREAMICE_COULOMB_SLIDING CPP option.

8.6.5.2.4. Configuring domain through files

The STREAMICE_GEOM_FILE_SETUP CPP option allows versatility in defining the domain. With this option, the array streamice_hmask must be initialized through a file (streamiceHmaskFile) as must streamice_ufacemask_bdry and streamice_vfacemask_bdry (through streamiceuFaceBdryFile and streamicevFaceBdryFile) as well as u_flux_bdry_SI and v_flux_bdry_SI, volume flux at the boundaries, where appropriate (through streamiceuMassFluxFile and streamicevMassFluxFile). Thickness must be initialized through a file as well (streamicethickFile); streamice_hmask is set to zero where ice thickness is zero, and boundaries between in-domain and out-of-domain cells (according to streamice_hmask) are no-slip by default.

When using this option, it is important that for all internal boundaries, streamice_ufacemask_bdry and streamice_vfacemask_bdry are -1 (this will not be the case if streamiceuFaceBdryFile and streamicevFaceBdryFile are undefined).

In fact, if streamice_hmask is configured correctly, streamice_ufacemask_bdry and streamice_vfacemask_bdry can be set uniformly to -1, UNLESS there are no-stress or flux-condition boundaries in the domain. Where streamice_ufacemask_bdry and streamice_vfacemask_bdry are set to -1, they will be overridden at (a) boundaries where streamice_hmask changes from 1 to -1 (which become no-slip boundaries), and (b) boundaries where streamice_hmask changes from 1 to 0 (which become calving front boundaries).

An example of domain configuration through files can be found in verification/halfpipe_streamice. By default, verification/halfpipe_streamice is compiled with STREAMICE_GEOM_FILE_SETUP undefined, but the user can modify this option. The file verification/halfpipe_streamice/input/data.streamice_geomSetup represents an alternative version of verification/halfpipe_streamice/input/data.streamice in which the appropriate binary files are specified.

8.6.5.2.5. Configuring domain through parameters

For a very specific type of domain the boundary conditions and initial thickness can be set via parameters in data.streamice. Such a domain will be rectangular. In order to use this option, the STREAMICE_GEOM_FILE_SETUP CPP flag should be undefined.

There are different boundary condition types (denoted within the parameter names) that can be set:

  • noflow: \(x\)- and \(y\)-velocity will be zero along this boundary.

  • nostress: velocity normal to boundary will be zero; there will be no tangential stress along the boundary.

  • fluxbdry: a mass volume flux is specified along this boundary, which becomes a boundary condition for the thickness advection equation (see Equations Solved). Velocities will be zero. The corresponding parameters flux_bdry_val_NORTH, flux_bdry_val_SOUTH, flux_bdry_val_EAST and flux_bdry_val_WEST then set the values.

  • CFBC: calving front boundary condition, a Neumann condition based on ice thickness and bed depth, is imposed at this boundary (see Equations Solved).

Note the above only apply if there is dynamic ice in the cells at the boundary in question. The boundary conditions are then set by specifying the above conditions over ranges of each (north/south/east/west) boundary. The division of each boundary should be exhaustive and the ranges should not overlap. Parameters to initialize boundary conditions (defined under STREAMICE_PARM03 namelist) are listed in Table 8.25.

Table 8.25 Parameters to initialize boundary conditions (defined under STREAMICE_PARM03 namelist)

Parameter

Default

Description

min_x_noflow_NORTH

0

western limit of no-flow region on northern boundary (m)

max_x_noflow_NORTH

0

eastern limit of no-flow region on northern boundary (m)

min_x_noflow_SOUTH

0

western limit of no-flow region on southern boundary (m)

max_x_noflow_SOUTH

0

eastern limit of no-flow region on southern boundary (m)

min_y_noflow_EAST

0

southern limit of no-flow region on eastern boundary (m)

max_y_noflow_EAST

0

northern limit of no-flow region on eastern boundary (m)

min_y_noflow_WEST

0

southern limit of no-flow region on western boundary (m)

max_y_noflow_WEST

0

northern limit of no-flow region on eastern boundary (m)

min_x_nostress_NORTH

0

western limit of no-stress region on northern boundary (m)

max_x_nostress_NORTH

0

eastern limit of no-stress region on northern boundary (m)

min_x_nostress_SOUTH

0

western limit of no-stress region on southern boundary (m)

max_x_nostress_SOUTH

0

eastern limit of no-stress region on southern boundary (m)

min_y_nostress_EAST

0

southern limit of no-stress region on eastern boundary (m)

max_y_nostress_EAST

0

northern limit of no-stress region on eastern boundary (m)

min_y_nostress_WEST

0

southern limit of no-stress region on western boundary (m)

max_y_nostress_WEST

0

northern limit of no-stress region on eastern boundary (m)

min_x_fluxbdry_NORTH

0

western limit of flux-boundary region on northern boundary (m)

max_x_fluxbdry_NORTH

0

eastern limit of flux-boundary region on northern boundary (m)

min_x_fluxbdry_SOUTH

0

western limit of flux-boundary region on southern boundary (m)

max_x_fluxbdry_SOUTH

0

eastern limit of flux-boundary region on southern boundary (m)

min_y_fluxbdry_EAST

0

southern limit of flux-boundary region on eastern boundary (m)

max_y_fluxbdry_EAST

0

northern limit of flux-boundary region on eastern boundary (m)

min_y_fluxbdry_WEST

0

southern limit of flux-boundary region on western boundary (m)

max_y_fluxbdry_WEST

0

northern limit of flux-boundary region on eastern boundary (m)

min_x_CFBC_NORTH

0

western limit of calving front condition region on northern boundary (m)

max_x_CFBC_NORTH

0

eastern limit of calving front condition region on northern boundary (m)

min_x_CFBC_SOUTH

0

western limit of calving front condition region on southern boundary (m)

max_x_CFBC_SOUTH

0

eastern limit of calving front condition region on southern boundary (m)

min_y_CFBC_EAST

0

southern limit of calving front condition region on eastern boundary (m)

max_y_CFBC_EAST

0

northern limit of calving front condition region on eastern boundary (m)

min_y_CFBC_WEST

0

southern limit of calving front condition region on western boundary (m)

max_y_CFBC_WEST

0

northern limit of calving front condition region on eastern boundary (m)

flux_bdry_val_SOUTH

0

volume flux per width entering at flux-boundary on southern boundary (m2/a)

flux_bdry_val_NORTH

0

volume flux per width entering at flux-boundary on southern boundary (m2/a)

flux_bdry_val_EAST

0

volume flux per width entering at flux-boundary on southern boundary (m2/a)

flux_bdry_val_WEST

0

volume flux per width entering at flux-boundary on southern boundary (m2/a)

8.6.5.3. Description

8.6.5.3.1. Equations Solved

The model solves for 3 dynamic variables: \(x\)-velocity (\(u\)), \(y\)-velocity (\(v\)), and thickness (\(h\)). There is also a variable that tracks coverage of fractional cells, discussed in Ice front advance.

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) [Mac89]). The SSA consists of the \(x\)-momentum balance:

(8.43)\[\partial_x(h\nu(4\dot{\varepsilon}_{xx}+2\dot{\varepsilon}_{yy})) + \partial_y(2h\nu\dot{\varepsilon}_{xy}) - \tau_{bx} = \rho g h \frac{\partial s}{\partial x}\]

the \(y\)-momentum balance:

(8.44)\[\partial_x(2h\nu\dot{\varepsilon}_{xy}) + \partial_y(h\nu(4\dot{\varepsilon}_{yy}+2\dot{\varepsilon}_{xx})) - \tau_{by} = \rho g h \frac{\partial s}{\partial y}\]

where \(\rho\) is ice density, \(g\) is gravitational acceleration, and \(s\) is surface elevation. \(\nu\), \(\tau_{bi}\) and \(\dot{\varepsilon}_{ij}\) are ice viscosity, basal drag, and the strain rate tensor, respectively, all explained below.

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

(8.45)\[h_t + \nabla\cdot(h\vec{u}) = \dot{a}-\dot{b}\]

Where \(\dot{b}\) is a basal mass balance (e.g., melting due to contact with the ocean), positive where there is melting. This is a field that can be specified through a file. At the moment surface mass balance \(\dot{a}\) can only be set as uniform. Where ice is grounded, surface elevation is given by

\[s = R + h\]

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

\[s = (1-\frac{\rho}{\rho_w}) h,\]

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

\[h \leq -\frac{\rho_w}{\rho}R\]

is satisfied. Floatation criteria is stored in float_frac_streamice, equal to 1 where ice is grounded, and equal to 0 where ice is floating.

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 Finite-volume discretization of the stress tensor divergence). Thus

\[\begin{split}\begin{aligned} \dot{\varepsilon}_{xx} = & u_x + k_1 v, \notag \\ \dot{\varepsilon}_{yy} = & v_y + k_1 u, \notag \\ \dot{\varepsilon}_{xy} = & \frac{1}{2}(u_y+v_x) + k_1 u + k_2 v. \notag \end{aligned}\end{split}\]

\(\nu\) has the form arising from Glen’s law

(8.46)\[\nu = \frac{1}{2}A^{-\frac{1}{n}}\left(\dot{\varepsilon}_{xx}^2+\dot{\varepsilon}_{yy} ^2+\dot{\varepsilon}_{xx}\dot{\varepsilon}_{yy}+\dot{\varepsilon}_{xy}^2+\dot{ \varepsilon}_{min}^2\right)^{\frac{1-n}{2n}}\]

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

(8.47)\[\vec{\tau}_b = C (|\vec{u}|^2+u_{min}^2)^{\frac{m-1}{2}}\vec{u}.\]

Again, the form is slightly different if a hybrid formulation is to be used, and the velocity refers to sliding velocity (\(u_b\)).

An alternative to the above “power law” sliding parameterization can be used by defining the STREAMICE_COULOMB_SLIDING CPP option and setting streamice_allow_reg_coulomb to .TRUE.:

(8.48)\[\vec{\tau}_b = C\frac{|u|^{m}N}{2\left[C^{1/m}|u|+(0.5N)^{1/m}\right]^{m}}u^{-1}\vec{u}\]

where \(u\) is shorthand for the regularized norm in (8.47) (or for \(u_b\) if a hybrid formulation is used). \(m\) is the same exponent as in (8.47). \(N\) is effective pressure:

(8.49)\[N = \rho g (h - h_f),\]

with \(h_f\) the floatation thickness

\[h_f = max\left(0,-\frac{\rho_w}{\rho}R\right),\]

where \(R\) is bed elevation. This formulation was used in the MISMIP+ intercomparison tests [ADCD+16]. (8.49) assumes complete hydraulic connectivity to the ocean throughout the domain, which is likely only true within a few tens of kilometers of the grounding line. With this sliding relation, Coulomb sliding is predominant near the grounding line, with the yield strength proportional to height above floatation. Further inland sliding transitions to the power law relation in (8.47).

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:

(8.50)\[(h\nu(4\dot{\varepsilon}_{xx}+2\dot{\varepsilon}_{yy}))n_x + (2h\nu\dot{\varepsilon}_{xy})n_y = \frac{1}{2}g \left(\rho h^2 - \rho_w b^2\right)n_x\]
(8.51)\[(2h\nu\dot{\varepsilon}_{xy})n_x + (h\nu(4\dot{\varepsilon}_{yy}+2\dot{\varepsilon}_{xx}))n_y = \frac{1}{2}g \left(\rho h^2 - \rho_w b^2\right)n_y.\]

Here \(\vec{n}\) is the normal to the boundary, and \(b\) is ice base.

8.6.5.3.2. 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. pkg/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 accounts only for vertical shear. Such hybrid formulations have been shown to be valid over a larger range of conditions than SSA (Goldberg 2011) [Gol11].

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

\[\nu = \frac{1}{2}A^{-\frac{1}{n}}\left(\dot{\varepsilon}_{xx}^2+\dot{\varepsilon}_{yy} ^2+\dot{\varepsilon}_{xx}\dot{\varepsilon}_{yy}+\dot{\varepsilon}_{xy}^2+\frac{1 }{4}u_z^2+\frac{1}{4}v_z^2+\dot{\varepsilon}_{min}^2\right)^{\frac{1-n}{2n}}\]

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

8.6.5.3.3. Ice front advance

By default all mass flux across calving boundaries is considered lost. However, it is possible to account for this flux and potential advance of the ice shelf front. If streamice_move_front is TRUE, then a partial-area formulation is used.

The algorithm is based on Albrecht et al. (2011) [AMH+11]. In this scheme, for empty or partial cells adjacent to a calving front, a reference thickness \(h_{ref}\) is found, defined as an average over the thickness of all neighboring cells that flow into the cell. The total volume input over a time step is added to the volume of ice already in the cell, whose partial area coverage is then updated based on the volume and reference thickness. If the area coverage reaches 100% in a time step, then the additional volume is cascaded into adjacent empty or partial cells.

If streamice_calve_to_mask is 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 streamice_calve_mask is not equal to 1. This mask must be set through a binary input file to allow the front to advance past its initial position.

No calving parameterization is implemented in pkg/streamice. However, front advancement is a precursor for such a development to be added.

8.6.5.3.4. Units of input files

The inputs for basal traction (streamicebasalTracFile, C_basal_fric_const) and ice stiffness (streamiceGlenConstFile, B_glen_isothermal) require specific units. For ice stiffness (A in (8.46)), \(B=A^{-1/n}\) is specified; or, more accurately, its square root \(A^{-1/(2n)}\) is specified (this is to ensure positivity of B by squaring the input). The units of streamiceGlenConstFile and B_glen_isothermal are \(\mathrm{Pa}^{1/2}\ \mathrm{yr}^{1/(2n)}\) where \(n\) is n_glen.

streamicebasalTracFile and C_basal_fric_const initialize the basal traction (C in (8.47)). Again \(C^{1/2}\) is directly specified rather than C to ensure positivity. The units are \(\mathrm{Pa}^{1/2} (\mathrm{m }\ \mathrm{yr}^{-1})^{n_b}\) where \(n_b\) is n_basal_friction.

8.6.5.4. Numerical Details

STREAMICE stencil

Figure 8.14 Grid locations of thickness (h), velocity (u,v), area, and various masks.

STREAMICE masks

Figure 8.15 Hypothetical configuration, detailing the meaning of thickness and velocity masks and their role in controlling boundary conditions.

The momentum balance is solved via iteration on viscosity (Goldberg 2011 [Gol11]). At each iteration, a linear elliptic differential equation is solved via a finite-element method using bilinear basis functions. The velocity solution “lives” on cell corners, while thickness “lives” at cell centers (Figure 8.14). The cell-centered thickness is then evolved using a second-order slope-limited finite-volume scheme, with the velocity field from the previous solve. To represent the flow of floating ice, basal stress terms are multiplied by an array float_frac_streamice, a cell-centered array which determines where ice meets the floation condition.

The computational domain of pkg/streamice (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 pkg/streamice. 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_ufacemask_bdry): equal to -1 (interior face), 0 (no-slip), 1 (no-stress), 2 (calving stress front), or 4 (flux input boundary); when 4, then u_flux_bdry_SI must be initialized, through binary or parameter file

  • \(vfacemaskbdry\) (streamice_vfacemask_bdry): 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 one sets #define STREAMICE_GEOM_FILE_SETUP in STREAMICE_OPTIONS.h, the values are only relevant at the boundaries of the grid.

The values of \(umask\) and \(vmask\) determine which nodal values of \(u\) and \(v\) are involved in the solve for velocities. These masks are not configured directly by the user, but are re-initialized based on streamice_hmask, streamice_ufacemask_bdry and streamice_vfacemask_bdry at each time step. Figure 8.15 demonstrates how these values are set in various cells.

With \(umask\) and \(vmask\) appropriately initialized, subroutine streamice_vel_solve.F can proceed rather generally. Contributions 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.

8.6.5.5. Additional Features

8.6.5.5.1. PETSc

There is an option to use PETSc for the matrix solve component of the velocity solve, and this has been observed to give a 3- or 4-fold improvement in performance over the inbuilt conjugate gradient solver in a number of cases. To use this option, the CPP option ALLOW_PETSC must be defined, and MITgcm must be compiled with the -mpi flag (see Section 3.5.4). However, often a system-specific installation of PETSc is required. If you wish to use PETSc with pkg/streamice, please contact the author.

8.6.5.5.2. Boundary Stresses

The calving front boundary conditions (8.50) and (8.51) are intended for ice fronts bordering open ocean. However, there may be reasons to apply different Neumann conditions at these locations, e.g., one might want to represent force associated with ice melange, or to represent parts of the ice shelf that are not resolved, as in Goldberg et al. (2015) [GHJS15]. The user can then modify these boundary conditions in the form

\[(h\nu(4\dot{\varepsilon}_{xx}+2\dot{\varepsilon}_{yy}))n_x + (2h\nu\dot{\varepsilon}_{xy})n_y = \frac{1}{2}g \left(\rho h^2 - \rho_w b^2\right)n_x + \sigma n_x + \tau n_y\]
\[(2h\nu\dot{\varepsilon}_{xy})n_x + (h\nu(4\dot{\varepsilon}_{yy}+2\dot{\varepsilon}_{xx}))n_y = \frac{1}{2}g \left(\rho h^2 - \rho_w b^2\right)n_y + \sigma n_y + \tau n_x\]

In these equations, \(\sigma\) and \(\tau\) represent normal and shear stresses at the boundaries of cells. They are not specified directly, but through coefficients \(\gamma_{\sigma}\) and \(\gamma_{\tau}\):

\[\sigma = \frac{1}{2}g \left(\rho h^2 - \rho_w b^2\right)\gamma_{\sigma}\]
\[\tau = \frac{1}{2}g \left(\rho h^2 - \rho_w b^2\right)\gamma_{\tau}\]

\(\gamma_{\sigma}\) is specified through streamiceuNormalStressFile, streamicevNormalStressFile, streamiceuNormalTimeDepFile, streamicevNormalTimeDepFile and \(\gamma_{\tau}\) is specified through streamiceuShearStressFile, streamicevShearStressFile, streamiceuShearTimeDepFile, and streamicevShearTimeDepFile. Within the file names, the u and v determine whether the values are specified along horizontal (\(u\)-) faces and vertical (\(v\)-) faces. The values will only have an effect if they are specified along calving front boundaries (see Configuring domain through files).

8.6.5.6. Adjoint

The STREAMICE package is adjointable using both TAF (Goldberg et al. 2013 [GH13]) and OpenAD (Goldberg et al. 2016 [GNHU16]). In OpenAD, the fixed-point method of [Chr94] is implemented, greatly reducing the memory requirements and also improving performance when PETSc is used.

Verification experiments with both OpenAD and TAF are located in the verification/halfpipe_streamice (see below).

8.6.5.7. Key Subroutines

Top-level routine: streamice_timestep.F (called from model/src/do_oceanic_phys.F)

   CALLING SEQUENCE
...
 streamice_timestep (called from DO_OCEANIC_PHYS)
 |
 |-- #ifdef ALLOW_STREAMICE_TIMEDEP_FORCING
 |    STREAMICE_FIELDS_LOAD
 |   #endif
 |
 |--#if (defined (ALLOW_STREAMICE_OAD_FP))
 |    STREAMICE_VEL_SOLVE_OPENAD
 |  #else
 |    STREAMICE_VEL_SOLVE
 |    |
 |    |-- STREAMICE_DRIVING_STRESS
 |    |
 |    | [ITERATE ON FOLLOWING]
 |    |
 |    |-- STREAMICE_CG_WRAPPER
 |    |   |
 |    |   |-- STREAMICE_CG_SOLVE
 |    |       #ifdef ALLOW_PETSC
 |    |        STREAMICE_CG_SOLVE_PETSC
 |    |       #endif
 |    |
 |    |-- #ifdef STREAMICE_HYBRID_STRESS
 |         STREAMICE_VISC_BETA_HYBRID
 |        #else
 |         STREAMICE_VISC_BETA
 |        #endif
 |
 |-- STREAMICE_ADVECT_THICKNESS
 |   |
 |   |-- STREAMICE_ADV_FRONT
 |
 |-- STREAMICE_UPD_FFRAC_UNCOUPLED
 |

8.6.5.8. STREAMICE diagnostics

Diagnostics output is available via the diagnostics package (Packages II - Diagnostics and I/O). Available output fields are summarized in the following table:

----------------------------------------------------------------------------
<-Name->|Levs|  mate |<- code ->|<--  Units   -->|<- Tile (max=80c)
----------------------------------------------------------------------------
SI_Uvel |  1 |       |UZ      L1|m/a             |Ice stream x-velocity
SI_Vvel |  1 |       |VZ      L1|m/a             |Ice stream y-velocity
SI_Thick|  1 |       |SM      L1|m               |Ice stream thickness
SI_area |  1 |       |SM      L1|m^2             |Ice stream cell area coverage
SI_float|  1 |       |SM      L1|none            |Ice stream grounding indicator
SI_hmask|  1 |       |SM      L1|none            |Ice stream thickness mask
SI_usurf|  1 |       |SM      L1|none            |Ice stream surface x-vel
SI_vsurf|  1 |       |SM      L1|none            |Ice stream surface y-vel
SI_ubase|  1 |       |SM      L1|none            |Ice stream basal x-vel
SI_vbase|  1 |       |SM      L1|none            |Ice stream basal y-vel
SI_taubx|  1 |       |SM      L1|none            |Ice stream basal x-stress
SI_tauby|  1 |       |SM      L1|none            |Ice stream basal y-stress
SI_selev|  1 |       |SM      L1|none            |Ice stream surface elev

8.6.5.9. Experiments and tutorials that use streamice

The verification/halfpipe_streamice experiment uses pkg/streamice.