Hydrology Solution - GlaDS

Description

The two-dimensional Glacier Drainage System model (GlaDS, [Werder2013]) couples a distributed water sheet model — a continuum description of a linked cavity drainage system [Hewitt2011] — with a channelized water flow model — modeled as R channels [Rothlisberger1972,Nye1976]. The coupled system collectively describes the evolution of hydraulic potential , water sheet thickness , and water channel cross-sectional area .

Sheet model equations

• Mass conservation: The mass conservation equation describes water storage changes over longer timescales (dictated by cavity opening due to sliding) as well as shorter timescales (e.g. due to surface melt water input):

  

  where: is the englacial void ratio, is water density (kg~m^-3), is gravitational acceleration (kg~m^-3), is the hydraulic potential (Pa), and is the sheet thickness (m). The water discharge (m²~s^-1) is given by:

  

  where is the sheet conductivity (m~s~kg^-1), and 5/4 and 3/2 are two constant exponents. Finally, the melt source term (m~s^-1) is given by:

  

  where is the geothermal heat flux (W~m^-2), is the frictional heating (W~m^-2), for basal stress (Pa), is ice density (kg~m^-3), and is latent heating (J~kg^-1).

• Sheet thickness:

  

  Here, is the cavity opening rate due to sliding over bed topography (m~s^-1), given by:

  

  where and are both constants (m), and is the basal sliding velocity vector (provided by the ice flow model). The cavity closing rate due to ice creep (m~s^-1), is given by:

  

  where is the basal flow parameter (Pa^-3~s^-1), related to the ice hardness by ^-1/3, is the Glen flow relation exponent, and is the effective pressure. The overburden hydraulic potential is given by , with the ice pressure and elevation potential , all of which are given in units of Pa.

Channel model equations

• Channel discharge (along mesh edges):

  

  where is the channel discharge (m³~s^-1), which evolves with respect to the horizontal coordinate along the channel , is the channel potential energy dissipated per unit length and time (W~m^-1), is the channel pressure melting/refreezing (W~m^-1), is the channel closing rate (m²~s^-1) and is the source term (m²~s^-1). The discharge is defined as:

  

  where is the channel conductivity, and 3 and 2 are two exponents. The term is the closing of the channels by ice creep, and is given by:

  

  where is the channel cross-sectional area (m²). Finally, , the channel source term balancing the flow of water out of the adjacent sheet, is:

  

  where is the normal to the channel edge.

• Channel dissipation of potential energy:

  

  where is the channel width (m), and is the discharge in the sheet (flowing in the direction of the channel; m²~s^-1), given by:

  

  with , , and as given above.

• Channel melt/refreeze rate:

  

  Here, is the Clapeyron slope (K~Pa^-1), is the specific heat capacity of water (J~kg^-1~K^-1), and is a switch parameter that accounts for the fact that the channel area cannot be negative, turning off the sheet flow for refreezing as , i.e.:

  

• Cross-sectional channel area (defined along mesh edges):

  

Boundary conditions

Boundary conditions for the evolution of hydraulic potential are applied on the domain boundary , as either a prescribed pressure or water flux. The Dirichlet boundary condition is:

where is a specific potential, and the Neumann boundary condition is:

corresponding to the specific discharge

Channels are defined only on element edges and are not allowed to cross the domain boundary, so we do not require flux conditions. Since the evolution equations for and do not contain their spatial derivatives, we do not require any boundary conditions for their evolution equations.

Model parameters

The parameters relevant to the GlaDS (hydrologyglads) solution can be displayed by running:

>> md.hydrology

• md.hydrology.pressure_melt_coefficient: Pressure melt coefficient () [K Pa^-1]
• md.hydrology.sheet_conductivity: sheet conductivity () [m^7/4 kg^-1/2]
• md.hydrology.cavity_spacing: cavity spacing () [m]
• md.hydrology.bump_height: typical bump height () [m]
• md.hydrology.ischannels: Do we allow for channels? 1: yes, 0: no
• md.hydrology.channel_conductivity: channel conductivity () [m^3/2 kg^-1/2]
• md.hydrology.spcphi: Hydraulic potential Dirichlet constraints [Pa]
• md.hydrology.neumannflux: water flux applied along the model boundary (m² s^-1)
• md.hydrology.moulin_input: moulin input () [m³ s^-1]
• md.hydrology.englacial_void_ratio: englacial void ratio ()
• md.hydrology.requested_outputs: additional outputs requested?
• md.hydrology.melt_flag: User specified basal melt? 0: no (default), 1: use md.basalforcings.groundedice_melting_rate

Running a simulation

To run a transient standalone subglacial hydrology simulation, use the following commands:

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

%Change hydrology class to GlaDS;
md.hydrology = hydrologyglads();

%Set model parameters here;
md = solve(md, 'Transient');

References