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 OpenADgenerated ice model adjoint and to provide synchronous iceocean coupling through the SHELFICE package. It solves a set of dynamic equations appropriate for floating iceshelf flow as well as icestream and slower icesheet 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. Compiletime 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.
CPP Flag Name 
Default 
Description 

#define 
explicit construction of matrix for Picard iteration for velocity 

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

#undef 
use package array for rLow rather than model 

#undef 
use files rather than parameters in STREAMICE_PARM03 to configure boundaries 

#undef 
enable interface to PETSc for velocity solver matrix solve 

#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 runtime by setting useSTREAMICE to .TRUE.
in file data.pkg
.
8.6.5.2.3. Runtime parmeters: general flags and parameters¶
Runtime 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.23.
Parameter 
Default 
Description 

910 
the (uniform) density of land ice (kg/m^{3}) 

1024 
the (uniform) density of ocean (kg/m^{3}) 

3 
Glen’s Flow Law exponent (nondim.) 

1e12 
minimum strain rate in Glen’s Law (\(\varepsilon_0\), yr^{1}) 

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

0 
exponent in nonlinear sliding law (nondim.) 

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

TRUE 
lower CG tolerance when nonlinear residual decreases by fixed factor 

2000 
maximum iterations in linear solve 

0 
as above when coupled with pkg/shelfice 

1e6 
tolerance of nonlinear residual for velocity (relative to initial) 

100 
maximum Picard iterations in solve for velocity 

0 
as above when coupled with pkg/shelfice 

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

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

FALSE 
terminate velocity iteration based on relative change per iteration 

TRUE 
terminate velocity iteration based on residual 

FILE 
method by which to initialize thickness ( 

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

FALSE 
allow ice shelf front to advance 

FALSE 
if streamice_move_front TRUE do not allow to advance beyond streamice_calve_mask 

' ' 
file to initialize streamice_calve_mask 

FALSE 
do not update ice thickness (velocity solve only) 

0.5 
CFL factor which determine maximum time step for thickness subcycling 

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

UNIFORM 
method by which to initialize basal traction ( 

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

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

UNIFORM 
method by which to initialize Glen’s constant ( 

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

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

' ' 
file to initialize timeindep melt rate (m/yr) 

' ' 
file to initialize timevarying melt rate (m/yr), based on streamice_forcing_period 

' ' 
topography initialization file (m); requires #define USE_ALT_RLOW 

' ' 
streamice_hmask initialization file; requires #define STREAMICE_GEOM_FILE_SETUP 

' ' 
streamice`STREAMICE_ufacemask_bdry` initialization file; requires #define STREAMICE_GEOM_FILE_SETUP 

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

' ' 
mass flux at \(u\)faces init. file (m^{2}/yr); requires #define STREAMICE_GEOM_FILE_SETUP 

' ' 
mass flux at \(v\)faces init. file (m^{2}/yr); requires #define STREAMICE_GEOM_FILE_SETUP 

' ' 
timedepend. mass flux at \(u\)faces file (m^{2}/yr); requires #define STREAMICE_GEOM_FILE_SETUP 

' ' 
timedepend. mass flux at \(v\)faces file (m^{2}/yr); requires #define STREAMICE_GEOM_FILE_SETUP 

' ' 
calving front normal stress parm along \(u\)faces (nondim.; see Boundary Stresses) 

' ' 
calving front normal stress parm along \(v\)faces (nondim.; see Boundary Stresses) 

' ' 
calving front normal stress parm along \(u\)faces (nondim.; see Boundary Stresses) 

' ' 
calving front normal stress parm along \(v\)faces (nondim.; see Boundary Stresses) 

' ' 
timedependent version of streamiceuNormalStressFile 

' ' 
timedependent version of streamicevNormalStressFile 

' ' 
timedependent version of streamiceuShearStressFile 

' ' 
timedependent version of streamicevShearStressFile 

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

0 
file input frequency for streamice timedependent forcing fields (s) 

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

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 indomain and outofdomain cells (according to streamice_hmask) are noslip 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 nostress or fluxcondition 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 noslip 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.24.
Parameter 
Default 
Description 

0 
western limit of noflow region on northern boundary (m) 

0 
eastern limit of noflow region on northern boundary (m) 

0 
western limit of noflow region on southern boundary (m) 

0 
eastern limit of noflow region on southern boundary (m) 

0 
southern limit of noflow region on eastern boundary (m) 

0 
northern limit of noflow region on eastern boundary (m) 

0 
southern limit of noflow region on western boundary (m) 

0 
northern limit of noflow region on eastern boundary (m) 

0 
western limit of nostress region on northern boundary (m) 

0 
eastern limit of nostress region on northern boundary (m) 

0 
western limit of nostress region on southern boundary (m) 

0 
eastern limit of nostress region on southern boundary (m) 

0 
southern limit of nostress region on eastern boundary (m) 

0 
northern limit of nostress region on eastern boundary (m) 

0 
southern limit of nostress region on western boundary (m) 

0 
northern limit of nostress region on eastern boundary (m) 

0 
western limit of fluxboundary region on northern boundary (m) 

0 
eastern limit of fluxboundary region on northern boundary (m) 

0 
western limit of fluxboundary region on southern boundary (m) 

0 
eastern limit of fluxboundary region on southern boundary (m) 

0 
southern limit of fluxboundary region on eastern boundary (m) 

0 
northern limit of fluxboundary region on eastern boundary (m) 

0 
southern limit of fluxboundary region on western boundary (m) 

0 
northern limit of fluxboundary region on eastern boundary (m) 

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

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

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

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

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

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

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

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

0 
volume flux per width entering at fluxboundary on southern boundary (m^{2}/a) 

0 
volume flux per width entering at fluxboundary on southern boundary (m^{2}/a) 

0 
volume flux per width entering at fluxboundary on southern boundary (m^{2}/a) 

0 
volume flux per width entering at fluxboundary on southern boundary (m^{2}/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 lowfriction bed (e.g., Macayeal (1989) [Mac89]). The SSA consists of the \(x\)momentum balance:
the \(y\)momentum balance:
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:
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
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
where \(\rho_w\) is a representative ocean density, and \(b=(\rho/\rho_w)h\). Again by hydrostasy, floation is assumed wherever
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 Finitevolume discretization of the stress tensor divergence). Thus
\(\nu\) has the form arising from Glen’s law
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 cellbycell 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
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.
:
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:
with \(h_f\) the floatation thickness
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:
Here \(\vec{n}\) is the normal to the boundary, and \(b\) is ice base.
8.6.5.3.2. Hybrid SIASSA 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 depthindependent. 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 depthaveraged \(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
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 partialarea 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¶
The momentum balance is solved via iteration on viscosity (Goldberg 2011 [Gol11]). At each iteration, a linear elliptic differential equation is solved via a finiteelement method using bilinear basis functions. The velocity solution “lives” on cell corners, while thickness “lives” at cell centers (Figure 8.14). The cellcentered thickness is then evolved using a secondorder slopelimited finitevolume 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 cellcentered 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 (icecovered), 0 (open ocean), 2 (partlycovered), 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 (noslip), 1 (nostress), 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 reinitialized 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 4fold 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 systemspecific 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
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}\):
\(\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 fixedpoint 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¶
Toplevel 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 L1m/a Ice stream xvelocity
SI_Vvel  1  VZ L1m/a Ice stream yvelocity
SI_Thick 1  SM L1m Ice stream thickness
SI_area  1  SM L1m^2 Ice stream cell area coverage
SI_float 1  SM L1none Ice stream grounding indicator
SI_hmask 1  SM L1none Ice stream thickness mask
SI_usurf 1  SM L1none Ice stream surface xvel
SI_vsurf 1  SM L1none Ice stream surface yvel
SI_ubase 1  SM L1none Ice stream basal xvel
SI_vbase 1  SM L1none Ice stream basal yvel
SI_taubx 1  SM L1none Ice stream basal xstress
SI_tauby 1  SM L1none Ice stream basal ystress
SI_selev 1  SM L1none Ice stream surface elev
8.6.5.9. Experiments and tutorials that use streamice¶
The verification/halfpipe_streamice experiment uses pkg/streamice.