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 \(\phi\), water sheet thickness \(h\), and water channel cross-sectional area \(S\).

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):

    \(\frac{e_v}{\rho_w g}\frac{\partial\phi}{\partial t} + \frac{\partial h}{\partial t} - \nabla\cdot\boldsymbol{q} - m_b = 0,\) where: \(e_v\) is the englacial void ratio, \(\rho_w\) is water density (kg~m-3), \(g\) is gravitational acceleration (kg~m-3), \(\phi\) is the hydraulic potential (Pa), and \(h\) is the sheet thickness (m). The water discharge \(\boldsymbol{q}\) (m2~s-1) is given by:

    \(\boldsymbol{q}=-k_s\;h^{\alpha_s}\left|\nabla\phi\right|^{\beta_s-2}\nabla\phi,\) where \(k_s\) is the sheet conductivity (m~s~kg-1), and \(\alpha_s=\)5/4 and \(\beta_s=\)3/2 are two constant exponents. Finally, the melt source term \(m_b\) (m~s-1) is given by:

    \(m_b=\frac{G+|\boldsymbol{\tau}_b\cdot\boldsymbol{u}_b|}{\rho_{i}L},\) where \(G\) is the geothermal heat flux (W~m-2), \(|\boldsymbol{\tau}_b\cdot\boldsymbol{u}_b|\) is the frictional heating (W~m-2), for basal stress \(\boldsymbol{\tau}_b\) (Pa), \(\rho_i\) is ice density (kg~m-3), and \(L\) is latent heating (J~kg-1).

  • Sheet thickness:

    \(\frac{\partial h}{\partial t} = w_s - v_s.\) Here, \(w_s\) is the cavity opening rate due to sliding over bed topography (m~s-1), given by:

    \(w_s\left(h\right) =\begin{array}{ll}\displaystyle \frac{\left|\boldsymbol{u}_b\right|}{l_r}\left(h_r-h\right), & \text{if } h<h_r\\\\0, & \text{otherwise,}\end{array}\) where \(h_r\) and \(l_r\) are both constants (m), and \(\boldsymbol{u}_b\) is the basal sliding velocity vector (provided by the ice flow model). The cavity closing rate due to ice creep \(v_s\) (m~s-1), is given by:

    \(v_s\left(h,\phi\right) = \frac{2A}{n^n}h\left|N\right|^{n-1}N,\) where \(A\) is the basal flow parameter (Pa-3~s-1), related to the ice hardness by \(B=A\)-1/3, \(n\) is the Glen flow relation exponent, and \(N= \phi_0-\phi\) is the effective pressure. The overburden hydraulic potential is given by \(\phi_0 = \phi_m+p\), with the ice pressure \(p = \rho_i g H\) and elevation potential \(\phi_m = \rho_w g b\), all of which are given in units of Pa.

Channel model equations

  • Channel discharge (along mesh edges):

    \(\frac{\partial Q}{\partial s} + \frac{\Xi-\Pi}{L}\left(\frac{1}{\rho_i} - \frac{1}{\rho_w}\right) - v_c - m_c = 0,\) where \(Q\) is the channel discharge (m3~s-1), which evolves with respect to the horizontal coordinate along the channel \(s\), \(\Xi\) is the channel potential energy dissipated per unit length and time (W~m-1), \(\Pi\) is the channel pressure melting/refreezing (W~m-1), \(v_c\) is the channel closing rate (m2~s-1) and \(m_c\) is the source term (m2~s-1). The discharge \(Q\) is defined as:

    \(Q= -\underbrace{k_cS^{\alpha_c}\left|\frac{\partial\phi}{\partial s}\right|^{\beta_c-2}}_{K_c}\frac{\partial\phi}{\partial s},\) where \(k_c\) is the channel conductivity, and \(\alpha_c=\)3 and \(\beta_c=\)2 are two exponents. The term \(v_c\) is the closing of the channels by ice creep, and is given by:

    \(v_c\left(S,\phi\right) = \frac{2A}{n^n}S\left|N\right|^{n-1}N,\\\) where \(S\) is the channel cross-sectional area (m2). Finally, \(m_c\), the channel source term balancing the flow of water out of the adjacent sheet, is:

    \(m_c = \boldsymbol{q}\cdot\boldsymbol{n}|_{\partial{\Omega_{i1}} } +\boldsymbol{q}\cdot\boldsymbol{n}|_{\partial{\Omega_{i2}} }.\) where \(\boldsymbol{n}\) is the normal to the channel edge.

  • Channel dissipation of potential energy:

    \(\Xi(S,\phi)=\Biggl|Q\frac{\partial\phi}{\partial s}\Biggr| +\left|l_cq_c\frac{\partial\phi}{\partial s}\right|,\) where \(l_c\) is the channel width (m), and \(q_c\) is the discharge in the sheet (flowing in the direction of the channel; m2~s-1), given by:

    \(q_c\left(h,\phi\right)= -\underbrace{k_s h^{\alpha_s} \left|\frac{\partial\phi}{\partial s}\right|^{\beta_s-2}}_{K_s}\frac{\partial\phi}{\partial s},\) with \(k_s\), \(\alpha_s\), and \(\beta_s\) as given above.

  • Channel melt/refreeze rate:

    \(\Pi(S,\phi)=-c_tc_w\rho_w(Q+f l_c q_c)\frac{\partial\phi-\partial\phi_m}{\partial s},\) Here, \(c_t\) is the Clapeyron slope (K~Pa-1), \(c_w\) is the specific heat capacity of water (J~kg-1~K-1), and \(f\) is a switch parameter that accounts for the fact that the channel area cannot be negative, turning off the sheet flow for refreezing as \(S\rightarrow0\), i.e.:

    \[f =\left\{\begin{array}{ll}1, & \text{if }S>0 \text{ or } q_c\partial(\phi-\phi_m)\partial s>0\\0, & \text{otherwise}\end{array}\right.\]

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

    \[\frac{\partial S}{\partial t} = \frac{\Xi - \Pi}{\rho_i L} - v_c.\]

Boundary conditions

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

\(\phi=\phi_D \quad\text{on} \quad\partial\Omega_D,\) where \(\phi_D\) is a specific potential, and the Neumann boundary condition is:

\(\frac{\partial\phi}{\partial n}=\Phi_N \quad\text{on} \quad\partial\Omega_N,\) corresponding to the specific discharge

\(q_N=-k_s h^{\alpha_s}|\nabla\phi|^{\beta_s-2}\Phi_N.\) 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 \(h\) and \(S\) 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 (\(c_t\)) [K Pa-1]
  • md.hydrology.sheet_conductivity: sheet conductivity (\(k\)) [m7/4 kg-1/2]
  • md.hydrology.cavity_spacing: cavity spacing (\(l_r\)) [m]
  • md.hydrology.bump_height: typical bump height (\(h_r\)) [m]
  • md.hydrology.ischannels: Do we allow for channels? 1: yes, 0: no
  • md.hydrology.channel_conductivity: channel conductivity (\(k_c\)) [m3/2 kg-1/2]
  • md.hydrology.spcphi: Hydraulic potential Dirichlet constraints [Pa]
  • md.hydrology.neumannflux: water flux applied along the model boundary (m2 s-1)
  • md.hydrology.moulin_input: moulin input (\(Q_s\)) [m3 s-1]
  • md.hydrology.englacial_void_ratio: englacial void ratio (\(e_v\))
  • 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

  • Ian J. Hewitt. Modelling distributed and channelized subglacial drainage: the spacing of channels. J. Glaciol., 57(202):302-314, 2011.

  • J. F. Nye. Water flow in glaciers: jokulhlaups, tunnels and veins. J. Glaciol., 17(76):181-207, 1976.

  • H. Rothlisberger. Water pressure in intra-and subglacial channels. J. Glaciol., 11(62):177-203, 1972.

  • Mauro A. Werder, Ian J. Hewitt, Christian G. Schoof, and Gwenn E. Flowers. Modeling channelized and distributed subglacial drainage in two dimensions. J. Geophys. Res., 118:1-19, 2013.