Hydrology Solution - SHAKTI

Description

SHAKTI (Subglacial Hydrology and Kinetic, Transient Interactions) is a transient subglacial hydrology model that has flexible geometry and treats the entire domain with one set of governing equations, allowing for any drainage configuration to arise, which can include efficient (channelized) and inefficient (distributed) features. [Sommers2018]

Equations

The SHAKTI model is based upon governing equations that describe conservation of water mass, evolution of the system geometry, basal water flux (approximate momentum equation), and internal melt generation (approximate energy equation).

• Continuity equation (water mass balance):

  

  where is subglacial gap height, is the volume of water stored englacially per unit area of bed, is basal water flux, is melt rate, and is the input rate of water from the englacial to subglacial system.

• Basal gap dynamics (subglacial geometry):

  

  where is the subglacial gap height, is melt rate, is the ice flow law parameter, is the flow law exponent, is the overburden pressure of ice, is water pressure, is a dimensionless parameter governing opening by sliding, and is sliding velocity. According to this equation, the subglacial gap height evolves with time by: opening by both melt and sliding over bumps on the bed, and closing due to ice creep.

• Basal water flux (approximate momentum equation):

  

  where is subglacial gap height, is gravitational acceleration, is kinematic viscosity of water, is a dimensionless parameter controlling the nonlinear transition from laminar to turbulent flow (for turbulent flow, the flux becomes proportional to the square root of the head gradient), is the Reynolds number, and is hydraulic head:

  

Equation (3) is a key piece of our model formulation, in that it allows for a spatially and temporally variable hydraulic transmissivity in the system, and facilitates easeful transition between laminar and turbulent flow regimes. Most existing models prescribe a hydraulic conductivity and assume the flow to be turbulent everywhere. Our momentum equation acts to “unify” different drainage modes in a single model.

• Internal melt generation (energy balance at the bed):

  

  where is latent heat of fusion of water, is geothermal flux, is basal sliding velocity, is the stress exerted by the bed onto the ice, is basal water flux (discharge), is hydraulic head, is the pressure melt coefficient (Clapeyron constant), is the heat capacity of water, is density of water, and is water pressure. In words, melt is produced through a combination of geothermal flux, frictional heat due to sliding, and heat generated through internal dissipation (where mechanical energy is converted to thermal energy), minus the heat consumed due to changes in water pressure. We note that this form of the energy equation assumes that all heat produced is converted locally to melt and is not advected downstream. We assume that the ice and liquid water are consistently at the pressure melting point temperature. While these assumptions may not be strictly valid under certain real conditions that may have interesting implications, we leave that discussion for future work.

Following Werder et al. (2013), the englacial storage volume is a function of water pressure:

where is the englacial void ratio (zero for no englacial storage).

Equations (1), (2), (3), and (5) are combined to form a nonlinear partial differential equation (PDE) in terms of hydraulic head, :

With no englacial storage (), Eq. (7) takes the form of an elliptic PDE.

Defining the hydraulic “transmissivity”:

Equation (7) can be written more compactly as:

or simply as:

where the forcing () is a function of the relevant dependent variables. Within each time step, this nonlinear PDE is solved using a Picard iteration to establish the head () distribution.

Boundary conditions

Boundary conditions can be applied as either prescribed head (Dirichlet) conditions or as flux (Neumann) conditions. We typically apply a Dirichlet boundary condition of atmospheric pressure at the edge of the ice sheet, and Neumann boundary conditions (no flux or prescribed flux, which can be constant or time-varying) on the other boundaries of the subglacial drainage domain.

Model parameters

The parameters relevant to the SHAKTI (hydrologyshakti) solution can be displayed by running:

>> md.hydrology

• md.hydrology.head subglacial hydrology water head (m)
• md.hydrology.gap_height height of gap separating ice to bed (m)
• md.hydrology.bump_spacing characteristic bedrock bump spacing (m)
• md.hydrology.bump_height characteristic bedrock bump height (m)
• md.hydrology.englacial_input liquid water input from englacial to subglacial system (m/yr)
• md.hydrology.moulin_input liquid water input from moulins (at the vertices) to subglacial system (m³/s)
• md.hydrology.reynolds Reynolds’ number
• md.hydrology.neumannflux water flux applied along the model boundary (m²/s)
• md.hydrology.spchead water head constraints (NaN means no constraint) (m)
• md.hydrology.relaxation under-relaxation coefficient for nonlinear iteration
• md.hydrology.storage englacial storage coefficient (void ratio)

Running a simulation

To run a transient stand-alone subglacial hydrology simulation, use the following commands:

md.transient = deactivateall(md.transient);
md.transient.ishydrology = 1;

%Change hydrology class to SHAKTI
md.hydrology = hydrologyshakti();

%Set model paramters here

md = solve(md, 'Transient');

References