Basal Water Modeling

From Interactive System for Ice sheet Simulation
Jump to: navigation, search

Basal water distribution and its relation to basal traction is a key focus of modern ice sheet model development. As there are a paucity of data relating to basal processes, generic formulations are implemented in CISM, allowing a number of scenarios to be investigated quickly, and setting the stage for future efforts.

Basal processes involve a complex coupling of

  • ice and the latent heat associated with its transition to water,
  • rheological properties of underlying till or bedrock,
  • transfer of potential energy to the margins, and
  • friction.

This poorly understood dynamic is responsible for the streaming, or fast flow of ice. Streaming behavior is responsible for ice flowing an order of magnitude or more faster than surrounding ice, and 90% of ice flow into the oceans.

This work will focus on the interplay between the thermodynamics of basal melt and the transport of basal water. To quickly establish a basis to work from we will realize the following scheme, which will quickly enable the numerical experiments specified in assessment to be preformed.


Approach to Modeling Basal Water

CISM already has a working and well validated means of computing the temperature distribution within ice sheets. The values are computed with a finite difference discretization of a vertically rescaled diffusion-advection equation. The results of the temperature calculation are coupled to the ice dynamics through; (1) the effective viscosity in the ice, and (2) in the cases where basal melting is occurring, there is a user-definable interface to sliding.

Once basal melt rates are computed from the temperature distribution, the melt water is distributed according to water movement in the direction of decreasing hydraulic potential \phi. The hydraulic potential gradient (or water pressure gradient) is written as,

\nabla \phi = \rho_i g \nabla z_s + (\rho_w - \rho_i) g \nabla z_b - \nabla N,

where \rho_w is the water density, and z_b is the lower surface elevation of the ice mass.

A continuity equation for basal water flow is written

\frac{\partial w}{\partial t} = -\nabla \cdot \left (\mathbf{v}_w \, w
\right ) + \dot{m}_b,

where w is the thickness of the water layer, \mathbf{v}_w the vertically averaged velocity of water in the layer and \dot{m}_b the basal melting rate.

In steady-state, the basal melting rate must balance the water flux divergence, or

\dot{m}_b = \nabla \cdot \left (\mathbf{v}_w \, w \right ).

In this simplest case, basal hydrology is represented in terms of the sub-glacial water flux, which is a function of the basal melt rate.

The steady-state basal water flux \mathbf{\psi}_w = (\mathbf{v}_w \, w) is obtained by integrating the basal melt rate \dot{m}_b over the entire ice sheet, beginning at the maximum potential, and proceeding in the direction of decreasing hydraulic potential.

In practice, this can be done with the algorithm for balance flux of ice sheets described in detail by Budd 1996. A improvement upon that scheme is found in Quinn 1991, where fluxes are distributed across neighbors in a more computationally efficient and flux conservative manner.

Basal sliding is often taken as being proportional to the basal drag (or driving stress) and inversely proportional to the effective pressure, which is defined as the ice overburden pressure minus the basal water pressure. Weertman 1964 relates basal sliding velocity to the driving stress and a roughness factor. The presence of a thin water film underneath an ice sheet will reduce friction and thus bed roughness, hence facilitating basal sliding. In that sense it is possible to involve basal water layer thickness as inversely proportional to roughness Johnson:2002. However, water layer thickness can only be determined if the sub-glacial water velocity is known. We can assume water velocity constant over the entire basin, so that the sub-glacial water flux \mathbf{\psi}_w can be inversely related to bed roughness.

Hence, a reasonable method of estimating ice velocity based upon basal conditions is formulated. Our intention will be to iteratively run the model with various sliding relations that are driven by basal water until the surface velocity measurements are matched, at which point an empirical relation will have been reached that can be used to explore a host of scenarios (See Section assessment).

Algorithmic Approach

The Static Algorithm

Compute the effective pressure

Measurements from boreholes show that basal water is typically close to, but slightly less than the ice overburden pressure (Kamb 01). This indicates that some fraction of the overburden is being supported by the bed. If N represents pressure of the ice being supported by the bed, then

p_w = \rho_i g H - N.

A crude parametrization of this is provided by the relation, provided by (Alley 96)

N = \frac{c}{w},

where w is the basal water thickness, and c the proportionality that is a function of bed roughness.

To complete the parametrization, let's assume that that p_w=p_i=\rho_i g H if w> 0.2 m, and that the effective pressure is close to over burden if the w< .005 m.

Compute the pressure potential

The water will move down the potential gradient, \nabla \phi. This potential is given by

\phi =\rho_w g z_b + \rho_i g H - N.

After some rearrangement the potential gradient is

\nabla \phi = \rho_i g \nabla z_s + (\rho_w - \rho_i) g \nabla z_b - \nabla N,

Route water down the potential gradient

Quinn 1999 provides the algorithm used. The method is

  1. Sort the data from high to low pressure potential.
  2. The flux into the highest potential cell is equal to the melt rate. The flux out is equal to the flux in, and assigned to lower potential neighbors, proportionate to the downhill slope.
  3. The procedure is continued, in the order of the descending sort, the only difference is that now each cell has a flux in equal to the melt rate plus the flux in from up potential.

Compute thickness from flux

The fluxes through each cell provide a basis for determining the water thickness. The flux is

Q = V w,

hence if the velocity is known, then the water depth can be computed.

The velocity is assumed to be Gauckler-Manning,

V = \frac{k}{n}  R_h ^\frac{2}{3} \cdot S^\frac{1}{2},


V is the cross-sectional average velocity (ft/s, m/s)
k is a conversion constant equal to 1.486 for U.S. customary units or 1.0 for SI units
n is the Gauckler-Manning coefficient (independent of units)
R_h is the hydraulic radius (ft, m)
S is the slope of the water surface or the linear hydraulic head loss (ft/ft, m/m) (S=h_f/L)

A Dynamic Algorithm

Consider the water balance in a grid cell and denote the storage  S . Conservation dictates that

\frac{dS}{dt} = I - Q

with I inflow and Q outflow. A complete accounting for the inflow term is

I = f_s + f_n + f_g

where f_s is melt in the grid cell, f_n flux in from upstream, and f_g the flow into the cell from below, due to changes in the sediment distribution.

The velocities in the system should be described by the Manning equation

\mathbf{V} = \frac{1}{n}R^{2/3} s^{1/2} = \frac{1}{n} \left( \frac{Wh}{W + 2h} \right)^{2/3} s^{1/2}