Introduction
The englacial hydrologic system controls water transport from the surface to the base of glaciers, regulating the sensitivity of glaciers and ice sheets to environmental influences. The geometry of water transport pathways is poorly understood, but it determines the efficacy of water transport to the glacier base (Fountain and others, Reference Fountain, Schlichting, Jansson and Jacobel2005) and thus strongly affects subglacial drainage and ice flow dynamics (Schoof, Reference Schoof2010; Bartholomew and others, Reference Bartholomew2012; Banwell and others, Reference Banwell, Hewitt, Willis and Arnold2016). In this study, we investigate how the detection of fluid resonance in subsurface cracks (quasi-planar cavities at interfaces or arising from fractures, with small opening compared to the other two dimensions) and conduits (long channels with roughly equant cross-sectional dimensions) may be used to interrogate hidden water transport pathways in glaciers and ice sheets.
The englacial system has been studied through various direct and indirect methods. Direct speleological observations (e.g. Vatne, Reference Vatne2001), in situ monitoring of boreholes with cameras (e.g. Fountain and others, Reference Fountain, Schlichting, Jansson and Jacobel2005), as well as indirect imaging from various types of ice-penetrating radar (e.g. Catania and others, Reference Catania, Neumann and Price2008; Badgeley and others, Reference Badgeley2017), suggest that the englacial system can be composed of cracks, conduits or possibly both. However, these methods are limited in their scope and generally fail to inform water transport dynamics or internal geometry such as is predicted from theory (e.g. Röthlisberger, Reference Röthlisberger1972; Shreve, Reference Shreve1972; Spring and Hutter, Reference Spring and Hutter1982; Gulley and others, Reference Gulley, Benn, Müuller and Luckman2009). Therefore, we look to other methods that may provide more insight.
Unsteady fluid motion has been used to probe the subglacial hydrologic system for decades. For example, slug tests and borehole studies have been conducted to infer basal structures (Stone and others, Reference Stone, Clarke and Ellis1997; Kulessa and others, Reference Kulessa, Hubbard, Williamson and Brown2005), subglacial storage and subglacial drainage efficiency (Iken and others, Reference Iken, Fabri and Funk1996; Kulessa and others, Reference Kulessa, Hubbard, Williamson and Brown2005; Kavanaugh and others, Reference Kavanaugh, Moore, Dow and Sanders2010; Rada and Schoof, Reference Rada and Schoof2018). More recently, seismic and geodetic studies have used transients in internal and basal water pressures on hour to day timescales, to infer basal sliding and ice flow dynamics (Iken, Reference Iken1972; Engelhardt and Kamb, Reference Engelhardt and Kamb1997; Kavanaugh and Clarke, Reference Kavanaugh and Clarke2001; Kavanaugh and others, Reference Kavanaugh, Moore, Dow and Sanders2010; Andrews and others, Reference Andrews2018).
Oscillatory fluid motions, generally at much higher frequencies, have also been leveraged to infer subsurface transport structures. Fluid oscillations in a confining structure lead to resonance at distinct frequencies or ‘eigenmodes’ that can be related to excitation mechanisms and the geometry of a resonator. This phenomenon has seen application in oil and gas exploration and volcanology (Aki and others, Reference Aki, Fehler and Das1977; Holzhausen and Egan, Reference Holzhausen and Egan1986; Chouet, Reference Chouet1988; Chen and others, Reference Chen, Quan and Harris1996; Mondal, Reference Mondal2010; Liang and others, Reference Liang, Karlstrom and Dunham2020), as well as some studies in the glacial environment (Lipovsky and Dunham, Reference Lipovsky and Dunham2015; Roeoesli and others, Reference Roeoesli, Walter, Ampuero and Kissling2016; Podolskiy, Reference Podolskiy2020). Here, we aim to advance this particular field of research, by investigating the limits to which fluid resonance can be used to define a variety of englacial geometries.
Because the problem of fluid resonance in elastic-walled transport structures has been studied in other contexts, we present a synthesis of relevant results and an overview of resonant modes using the englacial system to constrain parameters. Through this investigation we find one mode – a coupled oscillation between a conduit and branching crack – of particular interest, and dedicate much of the subsequent text to understanding this mode. In particular, we explore the limits to which the coupled mode provides information about the geometry of cracks branching from a central conduit connected to the surface, considering the influence of the flow regime within the conduit, crack dip and the excitation function that initiates wave motion. Finally, we apply our results to published high-frequency pressure time series data, and discuss the limits of subsurface crack detectability using fluid resonance in both alpine and ice-sheet environments.
Theory
We model the englacial system as a vertical tube with a circular cross section that may be connected to one or more dipping rectangular cracks, as seen in Fig. 1. Generalizations to conduits with variable radii, non-circular cross-sectional geometry or along-axis curvature are fairly straightforward extensions, but are not explored here. We assume that the fluid in the englacial system is pure water (no bubbles or sediment). Additionally, we neglect any energy dissipation from elastic waves emitted from the crack and tube walls and assume that the surrounding ice is homogeneous (e.g. Tsai and Rice, Reference Tsai and Rice2010). This last assumption is not strictly valid for cracks at the glacier bed between ice and till or bedrock, but we do not consider bi-material interfaces or inhomogeneous material properties here.
Fig. 1. Idealized model of the englacial system showing multiple cracks branching from a central conduit. Cracks are defined by two length scales L x and L y and opening w 0. Conduit sections are defined by a length L and a radius R. Cracks may also be dipping at an angle θ with respect to the conduit axis. The free surface is denoted by a red triangle. h(t) is the height of the water surface in reference to an unperturbed fluid surface at atmospheric pressure.
Governing equations: fluid-filled conduit
In the conduit, we solve the Navier–Stokes equations for flow of a viscous, compressible fluid in a pipe. We linearize around a static background state upon which small amplitude perturbations are superimposed at long wavelengths compared to crack openings and the conduit radius. Thus we can consider the fluid pressure to be uniform in the radial direction and fluid velocity to be axisymmetric and unidirectional. Finally, we assume a linearized equation of state for water in the conduit, giving rise to a constant fluid sound speed. The description and numerical implementation of equations in this section essentially follows Liang and others (Reference Liang, Karlstrom and Dunham2020), which is the modeling framework on which our study is built.
We use cross-sectionally averaged governing equations for fluid motion in a conduit with a constant radius R and length L, which could apply to distinct conduit sections as depicted in Fig. 1. Cross-sectional averaging is defined by the operator 〈〉, acting on a generic function Φ of vertical position z, radial coordinate r and time t:
We will generally omit explicit listing of function arguments except when necessary in what follows. Using this width-averaged description we solve the following governing equations in the conduit sections:
where u z = 〈v z〉 is the cross-sectionally averaged vertical velocity, ρ is the fluid density with ρ0 its constant and uniform background value, p is the fluid pressure, τ is the viscous stress tensor and K f is the bulk modulus of water (numerical values for parameters are listed in Table 1). For our numerical experiments we assume fluid flow is fully developed, and the vertical component of the cross-sectionally averaged shear stress term becomes
with μ as the dynamic viscosity of water. The limits for which this assumption is valid will be evaluated and discussed later for specific resonant modes.
Table 1. Values for glacial parameters
To supplement the governing Eqns (2–4), we impose no slip boundary conditions at the conduit walls to ensure no fluid is transferred out or into the system. At the top of the conduit, the water surface is allowed to move freely, resulting in the linearized boundary condition:
where g is the gravity, P source is an externally imposed source function and h(t) is the displacement of the fluid surface around hydrostatic equilibrium. h(t) can be calculated using the cross-sectionally averaged fluid velocity at the surface:
Boundary conditions at the conduit base, assuming a basal connection to a crack as depicted in Fig. 1, can be stated as
where p t is the pressure in the crack evaluated at the conduit junction, q is the volumetric flux from the conduit into the crack and A c = πR 2 is the conduit cross-sectional area. p t can be calculated by considering the mass balance at the base of the conduit
where V c is the volume of the crack and ρ is fluid density. For long period oscillatory perturbations (much longer than L x/c 0 with L x a crack length and c 0 fluid sound speed), fluid compressibility effects are negligible compared to crack volume changes and Eqn (9) can be written as
In Eqn (10), the storativity C t is the characteristic response of the crack volume to a uniform opening pressure, neglecting inertia and viscous pressure drops within the crack. In general, it depends on both fluid and crack bulk moduli (Segall, Reference Segall2010). However, the bulk modulus of water, K f, is generally much larger than K c and we neglect its effects here. Crack storativity may then be calculated as
The storativity of a symmetric crack can be represented in terms of the effective elastic modulus $G^\ast = G/( 1-\nu )$
of ice, with shear modulus G and Poisson's ratio ν, as well as crack length L x and a constant κ
The scaling factor κ depends on the aspect ratio of the crack, a purely geometric quantity. It is found by calculating dV c/dp numerically using the displacement discontinuity method described above, combining Eqns (11) and (12) and solving for κ. dV c/dp is found by applying a uniform pressure to the crack and calculating change in volume. For glacial parameters (Table 1), κ = 0.4814 ± 0.0006 with respect to the variable crack length. Note that Eqn (12) is strictly valid only for cracks in a homogeneous elastic half space. To find pressure at the crack mouth p t we integrate Eqn (10) in time. For low-frequency perturbations the fluid is essentially incompressible, and p t can then be written using Eqn (7) as
For high-frequency perturbations, there is no flux into the crack and the conduit bottom becomes a constant velocity boundary condition, u z = 0. When cracks intersect the conduit between the surface and its base, conservation of fluid momentum and mass at the crack mouth gives rise to additional frequency dependent junction conditions.
Governing equations: fluid-filled crack
For the crack, we consider a 3-D tabular cavity defined by the coordinate system x, y and ξ, where ξ represents the direction of opening for the crack as seen in Fig. 1. We use width-averaged equations for fluid flow in a 3-D crack (Lipovsky and Dunham, Reference Lipovsky and Dunham2015; Liang and others, Reference Liang, Karlstrom and Dunham2020):
where w is the perturbed crack opening in reference to the crack center-line (perpendicular to the unperturbed opening dimension w 0). Width-averaged velocities u x and u y in the crack are expressed as
with i = x, y. Similar to the conduit, we assume fully developed flow in the crack for our numerical experiments, in which case the shear stress terms become
To close our system of crack equations, we apply no slip boundary conditions on the crack faces and zero normal velocity boundary conditions at the edges of the crack.
The non-local response of crack walls to pressure perturbations is modeled with quasi-static elasticity, using the displacement discontinuity method (Crouch and others, Reference Crouch, Starfield and Rizzo1983; Okada, Reference Okada1985). This involves discretizing the crack into rectangular elements subject to opening-mode displacements, which are stitched together to ensure displacement and stress continuity (Liang and others, Reference Liang, Karlstrom and Dunham2020). A symmetric and positive definite stiffness matrix K s then relates crack opening w(x, y, t) and pressure p(x, y, t) on the mesh to evaluate Eqn (16), as
We also use K s to calculate the crack effective bulk modulus, which is defined as
where V c,0 = L xL yw 0 is a reference volume and dV c/dp is approximated numerically as discussed above.
Coupling between conduit and crack
To simulate fluid flow between conduit and crack elements, we couple the conduit and crack system of equations through pressure continuity and a discontinuous jump in the velocity due to the flux q into and out of the crack across each coupling point z c:
where A f = 2πR w 0 is the cross-sectional area of the crack intersecting the conduit assuming crack dip is 0 °. For dipping cracks the effective area increases with increasing dip as explained below.
Next, we incorporate the volume exchange between components and crack elasticity into the governing equations. Assuming wavelengths of the perturbation are long compared to the crack opening, and combining Eqns (16), (17) and (20) we get
where ${1}/{\bar{K}} = {1}/{K_{\rm c}} + {1}/{K_{\rm f}}$
and q is the fluid flux into the crack from the conduit. The fluid flux is modeled as a point source excitation with a delta function δ(x − x c)δ(y − y c) at the coupling point (x c, y c) of the conduit in the crack coordinate system. This delta function is scaled by the calculated energy transfer from the conduit into the crack seen in the right-hand-side of Eqn (24).
If a crack is dipping at some angle θ with respect to the conduit, the energy exchange between the conduit and crack will change as the cross section becomes elliptical instead of circular. We account for this area change using an effective circular radius (Tang and Cheng, Reference Tang and Cheng1993). If we consider the semi-minor axis of the ellipse equal to the conduit radius, the semi-major axis can be written as
using the ellipticity $e_{\rm E} = {\sqrt {R_1^2 - R^2}}/{R_1}$
, we can write the effective radius as
where E(e E, π/2) is the complete elliptic integral of the second kind: $E( k ,\; \, \theta ) = \int _{0}^{\theta } \sqrt {1-k^2 \sin ^2( \theta ') } {\rm d}\theta '$
.
Excitation of wave motion
Natural excitation of oscillatory flows in the englacial environment could arise from impulsive sources (Gräff and others, Reference Gräff, Walter and Lipovsky2019) or continuous sources (Roeoesli and others, Reference Roeoesli, Walter, Ampuero and Kissling2016). We largely focus on the impulsive case for this study and use a Gaussian pressure pulse applied at the water surface to perturb our system
with amplitude P 0 and half width ΔT. Due to the linearity of our problem, we assume unit forcing P 0 = 1 Pa in this study. A wavelength scale for the Gaussian function is calculated as λex = c T × ΔT where ΔT is the full width at half maximum of the Gaussian in the time domain (Eqn (27)) and c T is the tube wave speed (defined in the next section). We note that Eqn (27) is used for its relative simplicity, but other impulsive excitation functions might be used to narrow the input frequency range, such as a chirp signal.
Numerical implementation
We solve the coupled system of equations using a 6th-order summation-by-parts finite difference scheme in space and a 4th-order Runge–Kutta method in time (O'Reilly and others, Reference O'Reilly, Dunham and Nordström2017). Boundary conditions in the conduit are weakly enforced using simultaneous approximation terms (Del Rey Fernández and others, Reference Del Rey Fernández, Hicken and Zingg2014; Erickson and others, Reference Erickson, O'Reilly and Nordström2019) and boundary conditions in the crack are strongly enforced. Pressure and velocity in the conduit are solved on regular collocated grids, while the crack requires velocity grids to be staggered with respect to the pressure grids to prevent a coordinate singularity at the crack center (Prochnow and others, Reference Prochnow, O'Reilly, Dunham and Petersson2017). To maintain stability and high-order accuracy, the viscosity term in Eqns (3), (14) and (15) are solved implicitly in time (Liang and others, Reference Liang, Karlstrom and Dunham2020).
Results
Interpreting resonant oscillations
The characteristics of eigenmodes depend on the geometry of the resonating structure, the fluid properties and the source–time function of excitation (e.g. Biot, Reference Biot1952; Krauklis, Reference Krauklis1962). We focus on the characteristic frequencies of eigenmodes here as a discriminator, noting that eigenmode damping rate and resulting time-dependent deformation patterns in the ice, if detectable, provide additional information.
In tube-like geometries, interface waves or ‘tube waves’ can involve either local or non-local elastic coupling to the fluid (Burridge and others, Reference Burridge, Kostek and Kurkjian1993; Sheriff and Geldart, Reference Sheriff and Geldart1995). In this study, we focus on long wavelength excitation compared to the conduit radius, in which conduit walls deform quasi-statically and a non-dispersive compressional wave results. The tube wave travels along the fluid–solid interface at a slightly slower speed than the acoustic wave speed (Biot, Reference Biot1952). For a straight circular borehole, this speed is
In this study, we neglect the effects of borehole along-axis curvature and cross-sectional ellipticity. However, natural englacial conduits can be sinuous along their axis and elliptical in cross section. These effects modify the long wavelength tube wave speed by changing the compliance of the borehole (Norris, Reference Norris1990; Burridge and others, Reference Burridge, Kostek and Kurkjian1993).
Tube waves in the conduit will resonate at ‘organ pipe’ frequencies. These resonant frequencies depend on the tube wave speed, the boundary conditions – whether the ends of the pipe are considered open or closed – and the conduit length (Lighthill, Reference Lighthill1978)
Here, n = 1, 2, 3… and m = 1, 3, 5… are integer harmonic numbers. For example, a vertical englacial conduit connected to a large cavity could resonate like an organ pipe open at both ends since both the conduit top and bottom are free pressure boundaries. Conversely, a conduit terminating at the base of a glacier is equivalent to an organ pipe with one end open (free surface at the top of the conduit) and one end closed (solid boundary at the bottom of the conduit).
Crack waves, also known as Krauklis waves, are a type of interface wave that manifest in fluid-filled cracks where one dimension is small compared to the other two (Krauklis, Reference Krauklis1962; Ferrazzini and Aki, Reference Ferrazzini and Aki1987). Crack waves involve coupled fluid and elastic solid motion and are excited when perturbation wavelengths exceed the scale length for elastic coupling (e.g. Dunham and Ogden, Reference Dunham and Ogden2012). For high-frequency perturbations, the crack will respond rigidly to perturbations and crack waves will not be excited. At low frequencies, crack walls become more compliant and crack waves will manifest resulting in lower phase velocities with decreasing frequency. Crack wave dispersion manifests as unequal spacing of resonant harmonics in a crack. Lipovsky and Dunham (Reference Lipovsky and Dunham2015) estimated the frequency spacing between crack wave modes as f n/f 1 = n 3/2.
We calculate crack wave frequencies using the normalized crack impedance, which is a transfer function between pressure and fluid velocity at the crack mouth
where $\hat p$
and $\hat u_z$
are the Fourier transformed pressure and velocity. The Fourier transform and its inverse are defined as
Note that we define the normalized crack impedance as the inverse of the crack transfer function in Liang and others (Reference Liang, O'Reilly, Dunham and Moos2017). Frequencies where |F(ω)| is maximized represent the most efficient exchange between the crack and conduit for crack waves.
Coupled wave motion
In addition to crack and tube wave modes, coupled resonant modes exist between the conduit and crack. One mode, identified as the ‘conduit reservoir mode’ in Liang and others (Reference Liang, Karlstrom and Dunham2020), will be referred to as the ‘coupled mode’ here. For the coupled mode, fluid may be considered incompressible in the conduit, allowing the entire fluid column to oscillate under the influence of gravity in and out of the elastic crack. Liang and others (Reference Liang, Karlstrom and Dunham2020) derived this resonant frequency accounting for fluid density stratification in the conduit. We assume constant density here, however, we will generalize the modeling of this mode by analytically solving for non-fully developed flow regimes in the conduit.
Starting from our linearized momentum balance in the conduit, we integrate Eqn (3) in the incompressible limit from z = 0 to L:
For fully developed flow, the shear stress term $\langle \nabla \cdot \tau \rangle _z$
is Eqn (5) and does not depend on conduit length or frequency. However, in the boundary layer limit this shear stress term is more complicated. To determine the shear stress term, we use the solution of Womersley (Reference Womersley1955) for flow in a tube subject to a driving pressure variation ∂p/∂z that is harmonic in time with frequency ω. The resulting velocity is
where $\alpha ( \omega ) = R \sqrt {{\omega \rho _0}/{\mu }}$
and J i are Bessel functions of the first kind, order i. Using the identity $\int _{0}^{\psi } \psi ' J_0( \psi ') {\rm d}\psi ' = \psi J_1( \psi )$
, the cross-sectionally averaged velocity is
where for notational convenience ψ(ω) = α(ω) i 3/2. Shear stress may then be computed using the Bessel function recurrence relation (2m/ψ) J m(ψ) = J m+1(ψ) + J m−1(ψ) as
Assuming quasi-static opening of the crack, if we substitute Eqns (6) and (13) into Eqn (33) and rewrite in terms of the cross-sectionally averaged fluid displacement (Eqn (7)), we obtain the equation of a damped harmonic oscillator
with natural frequency
and damping coefficient
For an under-damped harmonic oscillator the resulting resonant frequency can be calculated as
and the quality factor can be calculated as the energy stored divided by the energy lost in one oscillation
If we consider the fully developed limit ψ ≪ 1, so J m(ψ) → (ψ/2) m(1/(Γ(m + 1))) with Γ(x) the Gamma function, Eqn (39) loses its frequency dependence:
If we neglect viscosity, Eqn (37) becomes a simple harmonic oscillator and the coupled mode frequency reduces to the natural frequency (Eqn (38)).
This coupled mode is derived for a crack intersecting at the base of a conduit, however a similar oscillation can arise in a system with a crack intersecting anywhere along the conduit. In that case, we must account for partitioning of flow between the crack and the conduit based on the relative cross-sectional areas of the crack opening and the conduit. It is standard to express this energy exchange in terms of reflection and transmission coefficients (Lighthill, Reference Lighthill1978). For our system we use reflection and transmission coefficients similar to those derived in Liang and others (Reference Liang, O'Reilly, Dunham and Moos2017):
where χ = (c 0/A f)/(c T/A c). These reflection and transmission coefficients are frequency dependent, implying selective interactions of the crack with oscillatory flow in the conduit. In addition, Eqn (26) implies that A f increases with crack dip, such that |R| → 1 and flux into the crack is maximized as dip angle θ → 90°.
Numerical experiments
Here, we present and analyze model results based on the theory and numerical approach outlined in the previous sections to explore which resonant modes might manifest in typical englacial systems, and how appropriate interpretation of frequency spectra could lead to inference of hidden transport geometries. We first present our analysis of crack resonance in isolation to connect with previously published study, then we study wave motion in three coupled englacial geometries: a conduit connected to a basal crack, a conduit connected to a crack located somewhere between the surface and the base and finally a multi-crack system.
We excite wave motion in our numerical experiments by initiating a pressure pulse at the water surface, then analyze fluid pressure and particle velocity time series at locations in the conduit and in the crack. Frequency domain analysis permits an interpretation of the resulting power spectrum to infer the resonant frequencies in the system. All model parameters for our englacial system can be found in Table 1 and were informed by studies of conduits and cracks in the englacial system (Anandakrishnan and Alley, Reference Anandakrishnan and Alley1997; Vatne, Reference Vatne2001; Petrenko and Whitworth, Reference Petrenko and Whitworth2002; Fountain and others, Reference Fountain, Schlichting, Jansson and Jacobel2005; Catania and others, Reference Catania, Neumann and Price2008; West and others, Reference West, Larsen, Truffer, O'Neel and LeBlanc2010).
Crack resonance
We first focus on fluid resonance in cracks alone. Lipovsky and Dunham (Reference Lipovsky and Dunham2015), assuming a symmetric 2-D crack, showed that the period and quality factor (energy loss over an oscillation cycle) for the fundamental resonant mode of a crack may be used to infer crack length and opening. We reproduced their calculation in Fig. 2A (red contours), demonstrating that an extension to two length dimensions (a 3-D tabular crack) is comparable in the symmetric case (Fig. 2A). We note that the over-damped region is determined from the Lipovsky and Dunham (Reference Lipovsky and Dunham2015) analysis of quality factor. We do not consider boundary layer effects in the crack and thus generally underpredict quality factor for the crack in our model. Across this range of crack dimensions, fundamental frequencies derived from our 3-D simulations are up to two times smaller than in the 2-D simulations presented in Lipovsky and Dunham (Reference Lipovsky and Dunham2015).
Fig. 2. Fundamental crack resonant frequencies for symmetric and asymmetric cracks. (A) Fundamental mode of a symmetric tabular crack, comparing the 2-D model (one length dimension and opening) of Lipovsky and Dunham (Reference Lipovsky and Dunham2015) (red contours) to the 3-D model studied here (blue contours). The over-damped region is determined from Lipovsky and Dunham (Reference Lipovsky and Dunham2015); no resonance occurs in this part of parameter space. (B) Variation of the fundamental mode for an asymmetric 3-D crack for two different crack length dimensions L x and L y. Where red contours represent a crack thickness of 0.005 m and blue contours represent results for a crack thickness of 0.05 m. Undulations in the contours reflect finite numerical resolution of crack storativity, leading to $\sim \!5 \percnt$
uncertainty in crack dimensions.
In the englacial environment, cracks are likely to be asymmetric in their length dimensions. Figure 2B shows that crack asymmetry can result in non-unique fundamental frequencies. For example, if we consider a crack length of L x = 10 m, and a crack opening between 5 and 50 mm, the fundamental frequency could vary between 0.4 and 30 Hz depending on L y. This non-uniqueness implies that higher order eigenmodes are required to constrain both length scales of an asymmetric crack.
Basal crack
We next consider a coupled geometry with a conduit connected to a crack at its base. This models a moulin or borehole connected to a basal water cavity. For the calculations and model results shown below we consider a 100 m conduit with a 10 cm radius, similar to a borehole, and a symmetric basal crack of length L x = 5 m with an opening of 1 cm. We excite the system with a scale wavelength (λex = c T/ΔT) of 10 m or ΔT = 0.0075 s (Eqn (27)), following the crack resonance condition from Lipovsky and Dunham (Reference Lipovsky and Dunham2015).
Figure 3 shows a spatial time series of the water particle velocity (u z) in response to a pressure perturbation at the water surface, where velocity is scaled by the characteristic fluid pressure magnitude, v scaled = (ρ0 c/P 0) v z. This figure highlights how the surface pressure perturbation affects the overall fluid motion in the conduit through time. Multiple wave modes are evident but somewhat complicated to disentangle in the time domain. However, in the frequency domain these modes and their harmonics are easily distinguished. In Fig. 4, we show pressure time series and spectra associated with the simulation in Fig. 3. Figures 4A and 4B show how pressure changes through time in the conduit and the crack. Figure 4D shows the fast Fourier transform (FFT) of a pressure time series taken at 50 m down the conduit shown in Fig. 4C. In this time series, we see organ pipe modes and harmonics beginning at 5.6 Hz, three crack wave modes at 41, 68 and 103 Hz, and the coupled mode at 0.75 Hz corresponding to the fundamental eigenmode of the system. Crack modes correspond to the peaks in the normalized crack impedance defined in Eqn (31), seen in Fig. 4E.
Fig. 3. Fluid particle velocity (u z) throughout the conduit in time, for a basal crack system with conduit length L = 100 m, radius R = 10 cm, a symmetric crack of length L x = 5 m and opening of w 0 = 1 cm. The center of the Gaussian excitation pulse is at 0 s model time with a scaled wavelength of 10 m. Following excitation we follow the propagating tube wave down the conduit until it reaches the basal crack at 100 m. When the tube wave reaches the crack, most of the energy is reflected and subsequently resonates at the fundamental organ pipe frequency. Energy that enters the crack is transmitted back into the conduit in the form of crack waves (first seen from ~0.1–0.3 s). Although wave motion becomes more complex through time, we note a full conduit velocity polarity change at ~1–1.2 s. This bulk motion in the fluid shows the long period coupled mode in the basal crack system.
Fig. 4. Spatial time series and frequency spectra for resulting wave motion in an englacial geometry with a conduit connected to a basal crack (assuming fully developed flow). (A) Similar to Fig. 3 but showing pressure time series throughout the conduit–crack network. (B) Spatial time series in the crack. (C) Pressure time series in the conduit taken at 50 m depth noted as the dashed line in panel (A). (D) Frequency spectrum of the time series shown in (C). In this frequency spectra, we note the coupled frequency at 0.75 Hz, crack wave modes at frequencies 41, 68 and 103 Hz, and many harmonics of the fundamental organ pipe frequency at 5.6 Hz. (E) Normalized crack impedance |F| for the basal crack, where |F| maxima correspond to observed crack wave frequencies in the conduit.
The power spectrum for the basal crack case shows that the coupled mode has the highest amplitude of all modes shown in the conduit. High amplitude modes are most likely to be measurable in the field, as has been observed in the volcanic setting (Chouet and Dawson, Reference Chouet and Dawson2013). We therefore now focus on this coupled mode as a possible tool for englacial geometric inference and consider both fully developed and boundary layer effects. We first note that the inviscid frequency Eqn (38) has two limits associated with distinct restoring forces for the oscillation. If A c/ρ0 C t g ≫ 1, the frequency of the coupled mode does not depend on gravity and equals $f_{{\rm el}} = ( {1}/{2\pi }) \sqrt {( {A_c/\rho _0 C_t}) /{L}}$
. For A c/ρ0 C t g ≪ 1, Eqn (38) reduces to a gravity-dependent limit, $f_{\rm g} = \sqrt {g/L}$
where crack elasticity no longer matters. Figure 5 shows how the coupled mode frequency and the damping rate scale for the full range of englacial parameters described in Table 1. The x-axis of Fig. 5A illustrates data collapse as the wavelength associated with the elastic limit frequency, λel, is $\lesssim \!1$
. For λel/L ≫ 1, the data collapse onto the gravity-dependent frequency f g plotted as the dashed lines in Fig. 5A. The gravity-dominated limit does not depend on crack parameters and thus it is not possible to determine crack geometry when in this limit.
Fig. 5. Coupled mode frequency and damping rate regimes. (A) Dominant restoring force as a function of frequency f for the coupled mode in a basal crack system. In the elastic limit, the coupled mode frequency is $f_{{\rm el}} = ( {1}/{2\pi }) \sqrt {{( A_c/\rho _0 C_t) }/{L}}$
, while in the gravity limit frequency follows $f_{g} = ( {1}/{2\pi }) \sqrt {{g}/{L}}$
. The x-axis shows the wavelength λel = c T/f el normalized by conduit length L. The dashed lines represent the gravity limit frequency for three different conduit lengths of 10, 100 and 1000 m. This highlights where the data collapse onto the gravity limit. (B) Quality factor of the coupled mode is controlled by the ratio of mode period and momentum diffusion time across the conduit T visc = R 2ρ0/μ. The y-axis shows cross-sectionally averaged shear stress in the conduit (Eqn (39)) normalized by the fully developed limit γdev = 4μ/ρ0 R 2 = 4/T visc. This illustrates the importance of a boundary layer treatment of viscous damping for assessing quality factors in the glacial environment.
Quality factor for the coupled mode (Eqn (41)) depends on resonant frequency ω and frequency-dependent damping rate γ(ω) but not on flow inside the crack. Crack geometry enters the coupled mode only through the crack storativity, which does not constrain opening w 0. Dissipation is thus governed by flow in the conduit, where the low viscosity of water implies that flow may not be fully developed at the periods of interest. Figure 5B shows that the coupled mode for glacial parameters mostly falls in the boundary layer limit. Fully developed flow is only a valid assumption for small conduit radii and long conduit lengths.
We can determine the crack length from the coupled mode using both the quality factor and frequency. We first note that for very large cracks our model assumption of a uniform opening pressure in the crack breaks down. Using Eqns (40) and (41), we take 2-D slices through the 3-D parameter space (L, L x, R) to plot crack length versus period and conduit length, and contour both the frequency and quality factor as shown in Fig. 6. Alternatively, if one knew the conduit radius or conduit length and the frequency or quality factor, the crack length could be determined using a graphical approach similar to Fig. 6, with parameters tailored to the specific case.
Fig. 6. Coupled mode frequency and quality factor for the range of glacial parameters in Table 1. (A and B) Frequency contours for R = 0.1 and 0.5 m, and L = 100 and 1000 m. (C and D) Quality factor contours for R = 0.1 and 0.5 m, and L = 100 and 1000 m. In all panels, transitions to the gravity-dominated limit can be seen where the frequency and quality factor lose crack length dependence. Panel (B) in particular shows no frequency variation for crack lengths above the 0.06 Hz (red) and 0.02 Hz (blue) contour. Conduit lengths and radii were chosen as proxies for an alpine and ice-sheet environment.
Middle crack
Next we examine a crack located between the base of the conduit and the surface. Using the same parameters as the basal cavity case, we place a horizontal crack 30 m from the base of the conduit. Figures 7A and 7B show these results for the fully developed limit. Similar to the basal crack geometry, we see organ pipe modes (in both conduit sections, due to reflection and transmission of energy past the crack mouth) as well as crack wave modes and a coupled mode. The organ pipe mode associated with the top conduit section is ~8 Hz and is associated with a pipe open at both ends (Eqn (29)). The organ pipe frequency associated with the bottom conduit section is ~9 Hz, and is equivalent to a pipe open at one end and closed at the other (Eqn (30)). The difference in these modes is due to the no slip boundary condition at the conduit base versus the pressure boundary condition at the crack and conduit interface. We also see a coupled mode associated with flow from upper section of the conduit in and out of the crack as the fundamental mode. We do not see a coupled mode associated with the bottom section of the conduit because it is closed to flow in the incompressible limit. Finally, we see the same crack modes as those in the basal case since the crack dimensions have not changed. The crack impedance used to identify these modes is the same as that in Fig. 4D.
Fig. 7. Spatial time series and frequency spectra for two middle crack models with the same geometric parameters and forcing as the coupled basal crack simulation in Fig. 4, except that the crack is located at 70 m down the conduit and varies in dip between the two simulations. (A) Spatial time series for a crack dipping at 0°. (C) Spatial time series for a crack dipping at 70°. (B) Frequency spectra for a crack dipping at 0°. (D) Frequency spectra for a crack dipping at 70°. The thick black line in (A) and (C) represents the location of the crack and the gray dashed lines represent the location of the time series in the upper and lower conduits that are Fourier transformed to create the spectra in (B) and (D).
The middle crack geometry provides an ideal test bed for the role of crack dip. Englacial cracks are often observed to be dipping at fairly steep angles, ~ 70° (Fountain and others, Reference Fountain, Schlichting, Jansson and Jacobel2005). Figure 8 shows how dip angle affects the reflection and transmission of waves, using the same crack geometry as all previous experiments. The normalized crack impedance for this crack is shown in Fig. 4E and is used to calculate the reflection coefficient from Eqn (43), using effective radius for the crack area from Eqn (26).
Figures 8B and 8C show that the reflection coefficient increases with increasing dip angle, becoming largest for angles $\gtrsim \!75^{\circ }$
. The reflection coefficient is frequency dependent and controlled by the crack impedance. Fluid flux into the crack is determined as jump in flux across the crack, which is affected by the reflection and transmission coefficients as well as crack geometry.
These effects can be seen in Fig. 7 which compares a horizontal and 70° dipping crack at 70 m depth. Increased flux in the dipping case is reflected off of the crack tip and re-emitted to the conduit as stronger upward and downward propagating waves compared to the horizontal crack. The crack has a small length of 5 m compared to the incoming perturbation wavelength of ~ 10 m, resulting in two-way transit within the crack before the incoming wave has past. The increase in flux into and out of the dipping crack results in increased amplitude of the lower conduit organ pipe modes, coupled mode and crack modes (Figs 7B and D).
Multiple branching cracks
Finally, we present results from a multi-crack system to illustrate how the basic elements studied individually interact in a more complex englacial network. We consider three cracks located at depths of 40, 90 m and at the base of a 150 m long conduit. The cracks at 40 and 90 m have a length of 5 m whereas the crack at 150 m has a length of 10 m. The top crack at 40 m is dipping at a 70° angle and has an opening of 1 mm compared to 10 mm for the other two. All other model parameters are the same as previous experiments.
We use this example to illustrate the role of the pressure source–time function on eigenmode excitation. Figure 9 shows the results of three different multi-crack experiments, differing only in the form of P source(t) in Eqn (6). The plots in the left column are for a system forced impulsively with ΔT = 0.05 s or 130 m wavelength, and the plots in the center column are forced with ΔT = 0.0075 s (10 m wavelength). The plots in the right column were forced continually with an excitation function comprised of 2000 sine waves with random initial phase and a small amount of numerical noise, normalized to unit maximum. The resulting quasi-white noise source–time function is a model for continuous excitation, for example from falling water in a moulin.
Fig. 9. Excitation of resonant modes from different source–time functions in a three crack system where the cracks are at depths of 40, 90 and 150 m. All other model parameters are the same as previous experiments and assuming fully developed flow. Left column: A long wavelength Gaussian excitation, with a scaled wavelength of 130 m. (A) Forcing function for the long wavelength model run. (B) Spatial time series for the long wavelength model run. (C) Frequency spectra for the long wavelength model run showing the FFT of time series taken at 20, 65 and 120 m. Middle column: A short wavelength excitation, with a scaled wavelength of 10 m. (D) Forcing function for the short wavelength model run. (E) Spatial time series for the short wavelength model run. (F) Frequency spectra for the short wavelength model run showing the FFT of time series taken at 20, 65 and 120 m. Right column: A simulated white noise continuous excitation, comprised of 2000 sine waves with random initial phase and a small amount of numerical noise. (G) Forcing function for the continuous forcing model run. (H) Spatial time series for the continuous forcing model run. (I) Frequency spectra for the continuous forcing model run showing the FFT of time series taken at 20, 65 and 120 m.
The conduit pressure frequency spectrum for each multi-crack system exhibits the same set of eigenmodes, excited to various degrees according to the forcing frequency spectra. The short wavelength excitation results in many resonant modes due to a broad range of excitation frequencies, while the long wavelength excitation preferentially excited low-frequency modes. The continuous forcing, while generating complex motions in the time domain, is still composed of a superposition of the same clearly identifiable eigenmodes.
In general, the amplitude of eigenmode excitation can be predicted directly from the frequency spectra of the source–time function given an appropriate transfer function (Karlstrom and Dunham, Reference Karlstrom and Dunham2016). However we do not formulate such a transfer function here, or an explicit dependence on transport network complexity (e.g. number of cracks) as in Fig. 9. It is easy to imagine englacial geometries that are effectively unresolvable through fluid resonance measurements, such as a case when branching cracks are spaced more closely than the wavelength of tube waves, but we leave this problem to future study.
Discussion
Increasing resolution and availability of measurement technologies have opened up new pathways toward discovery in most sub-disciplines of glaciology. Studies of water movement through glaciers and ice sheets are no exception having seen a burst of activity driven by the proliferation of seismologic and geodetic measurement in recent years (e.g. Kavanaugh and others, Reference Kavanaugh, Moore, Dow and Sanders2010; Roeoesli and others, Reference Roeoesli, Walter, Ampuero and Kissling2016; Andrews and others, Reference Andrews2018; Rada and Schoof, Reference Rada and Schoof2018). The interpretation of such measurements relies on models. We have developed a quantitative framework to understand a class of unsteady fluid motions within glaciers and these resonant oscillations are controlled by the geometry and material properties of the hidden englacial transport network. Such signals could be detected by a range of geophysical measurements, including water pressure and seismicity time series. In this section we comment on the resolving power of the resonant modes of coupled conduit–crack networks, focusing on a coupled mode and conditions appropriate for alpine glaciers and ice sheets that are motivated by previously published data.
Detecting branching cracks from a conduit
We have investigated whether the geometry of branching cracks is encoded in low-frequency eigenmodes of a coupled conduit–crack network. This is motivated by the reasonable expectation that lower frequency signals will be easiest to observe due to less attenuation and scattering for seismic measurements. However, we also recognize that low-frequency eigenmodes suffer some non-uniqueness when it comes to parameter resolving power. Therefore, we must assess whether sufficient geometric information is contained in low-frequency eigenmodes to be interesting.
Low-frequency eigenmodes and in particular the coupled mode constrain branching crack length in some circumstances, as illustrated in Fig. 6. Figure 5 shows that the coupled mode exhibits two regimes, a gravity-dominated regime and an elasticity-dominated regime. Crack geometry can be determined when elasticity matters, as seen in Fig. 6 where frequency and quality factor have power law dependencies over a broad range of crack lengths less than L x ~ 100 m. The maximum detectable crack length increases for increasing borehole radius and decreasing conduit length, however it does not much exceed 100 m for reasonable glacial parameters. In addition, for quality factors Q < 0.5, the coupled mode is over-damped and will not generally encode a diagnostic period and quality factor. We find that the coupled mode is over-damped for a small portion of the glacial parameter space where conduit radii are small and conduit lengths are on the ice-sheet scale. Crack dip also plays an important role, by controlling the degree to which oscillatory flow in the conduit interacts with branching cracks. As shown in Fig. 8, the reflection coefficient increases with increasing dip angle. As the magnitude of the reflection coefficient tends to 1, more flow enters the crack to excite resonance. Figure 8B shows dip angle impacts higher frequency oscillations more than lower frequencies, but generally cracks with high dip angles will be more connected hydraulically with the conduit during fluid oscillation.
Applications to the field
These limitations notwithstanding, our study suggests that excitation of resonant modes in coupled conduit–crack networks may provide a tool to infer englacial geometries in a range of realistic conditions. Appropriate instrumentation of glacial conduits, whether natural or human-made, is often more challenging than the equivalent problem in industry due to the dynamic conditions experienced on glaciers. However, the ability to directly instrument the fluid, such as with high-frequency pressure sensors or cameras, makes the detection of fluid resonance possibly simpler than for the equivalent volcanic problem. Field applications of the model developed here likely require a case-by-case examination. We illustrate how this might be accomplished by examining two contrasting glacial settings, where high-frequency pressure time series inside a vertical borehole already exist.
As a first example, we consider the results of Gräff and others (Reference Gräff, Walter and Lipovsky2019), who conducted experiments in a drilled borehole on Rhonegletscher Glacier in the Swiss Alps. Although their passive observations are not identical to our numerical experiments forced by a surface pressure pulse, they would exhibit the same eigenmodes of the system that we model all else constant. At the drilling site, Rhonegletscher Glacier was 117 m thick and was penetrated by hot water drilling with a hole of radius 10 cm. Gräff and others (Reference Gräff, Walter and Lipovsky2019) noted both impulsive and a more sustained excitation of crack waves in their borehole, resulting in high-frequency pressure time series qualitatively similar to our numerical experiments. They recorded several distinct spectral peaks that were inferred to represent resonant crack wave modes of a basal water-filled crack.
We propose that the range of spectral peaks observed by Gräff and others (Reference Gräff, Walter and Lipovsky2019) can be interpreted according to the framework introduced here for a conduit and basal crack. According to this model, the lowest frequency mode observed by Gräff and others (Reference Gräff, Walter and Lipovsky2019) will be the coupled mode, with higher frequency peaks being associated with crack and organ pipe modes. The relative amplitude of these modes will vary with the excitation function, which was not well constrained but might include both stick slip at the glacier bed and surface forcing. We therefore focus primarily on interpreting the eigenmodes themselves.
Using Eqn (40) and ice properties shown in Table 1, we calculate the crack length to be ~4.4 m for a coupled mode to match the ~ 1 Hz lowest frequency in (Gräff and others, Reference Gräff, Walter and Lipovsky2019, Fig. 4). Our modeled quality factor of ~ 60 does not match with that reported of ~ 1.5, but does if we artificially increase viscous dissipation in the conduit by ~1000. This discrepancy suggests that our model does not account for all sources of dissipation present in the natural system. The next highest mode at ~ 2.3 Hz is consistent with an open-closed organ pipe mode for a 107 m water column in the conduit (the water surface was ~ 10 below the ice surface), assuming sufficient uncertainty in the estimated tube wave speed due variable ice properties to account for slight discrepancies between calculated organ pipe frequencies and observations. The ~ 8 Hz mode is consistent with a harmonic of this mode. We interpret the ~ 4 Hz mode to be the fundamental mode of the crack, since it cannot be explained by an organ pipe mode. However, we are unable to identify subsequent crack wave modes. Similar to Gräff and others (Reference Gräff, Walter and Lipovsky2019) we use Fig. 7 from Lipovsky and Dunham (Reference Lipovsky and Dunham2015) to determine a crack length using our interpreted fundamental frequency of 4 Hz. We determine a crack length of 5.8 m, which is remarkably close to the length inferred independently from the fundamental mode. Note that we cannot use results shown in Fig. 2 without modeling non-fully developed flow in the crack to accurately determine quality factor contours. Crack opening is predicted to be ~ 1 mm by the crack wave mode, but is unresolved by the coupled mode (crack storativity is sensitive to crack length only as in Eqn (12)).
Our interpretation differs from Gräff and others (Reference Gräff, Walter and Lipovsky2019), who associate the lowest observed frequency in the power spectrum with the fundamental crack eigenmode. This predicts a 19.8 m basal crack length. These authors also argue that the 2.3 ± 0.22 Hz peak does not reflect tube wave resonance, because it was not measured by a nearby surface seismometer. We do not attempt to assess this here, but argue that the relative self-consistency offered by our model for the full range of spectral peaks observed by Gräff and others (Reference Gräff, Walter and Lipovsky2019), despite its simplicity, is interesting and potentially compelling. Our results are consistent with a borehole connected to a thin, 4–5 m long, subglacial crack. In addition, if we consider the sustained events to be similar to the continuous forcing example in Fig. 9, differences between spectral content of an impulsive versus continuous surface forcing are in qualitative agreement with variations in spectral peak amplitude. In particular, reduction in amplitude of an open-closed organ pipe mode for the longer forcing wavelengths present in a continuous source is consistent with a coincident increase in amplitude of the coupled mode.
This result provides a measure of confidence that our model can reasonably be applied in the alpine environment. Limited availability of high-frequency data from ice sheets precludes a direct assessment of our methods in thick ice. But we consider representative experiments on the Greenland ice sheet as a template for whether englacial features should be detectable in that environment. In 2014 and 2015, drilled boreholes in SE Greenland near Isunnguata Sermia outlet glacier penetrated ~ 700 m thick ice (Meierbachtol and others, Reference Meierbachtol, Harper and Humphrey2018), and were instrumented with 4 Hz pressure sensors at the glacier bed. For the majority of their sampling period the conduit was filled with ~600 m of water. Assuming a conduit radius similar to the boreholes of 15 cm (Meierbachtol and others, Reference Meierbachtol, Harper and Humphrey2013), if a coupled mode is excited, we could constrain a crack length up to ~70 m before the gravity-dominated regime limits inference of crack dimensions.
Although no resonant oscillations in the pressure time series were reported by Meierbachtol and others (Reference Meierbachtol, Harper and Humphrey2018), transient pressure pulses with a duration of ~ 10–20 s were reported. The signals qualitatively resemble over-damped waves and might be interpreted as such. Based on the geometry of the conduit, the coupled mode might not be over-damped if excited. Therefore, we consider the possibility of this over-damped pulse originating from the basal water layer, given that the pressure transducer was located at base of the borehole. Figure 7 from Lipovsky and Dunham (Reference Lipovsky and Dunham2015) predicts that a 100 m long crack with an opening of ~ 1 mm would result in a 0.1 Hz over-damped oscillation. Table 1 from Meierbachtol and others (Reference Meierbachtol, Harper and Humphrey2018) shows most pressure pulse recovery times are around 10 s, which matches the predicted periods. We recognize that, in practice, over-damped signals alone are not likely to result in strong constraints on englacial or subglacial geometry. But their existence could be used to motivate future experiments.
These examples suggest that low-frequency eigenmodes represent a promising tool for inferring englacial geometry in ice-sheet settings as well as alpine glaciers. The predictive power of the coupled mode is maximized for measurements in short conduit lengths (L < 100 m) and conduit radii larger than 50 cm (in order to detect cracks up to ~100 m in length). In a controlled active source setting, one can control conduit radius (through drilling) which could help in designing experiments. Either continuous or impulsive forcing will excite low-frequency modes, and modulated frequency sweeping of input forcing could be leveraged to seek signatures of matched resonance (Liang and others, Reference Liang, O'Reilly, Dunham and Moos2017) or be otherwise tuned to detect cracks of different sizes.
It is also possible that excitation of resonant signals could be used to probe time-evolving englacial or subglacial geometry. For example, seasonal evolution of subglacial drainage pathways (Gimbert and others, Reference Gimbert, Tsai, Amundson, Bartholomaus and Walter2016) or lake drainage (e.g. Bigelow and others, Reference Bigelow2020) occur on timescales long enough that one would expect to see systematic changes in resonant periods and quality factors as subsurface geometry evolves. We leave further investigation of these possibilities to future study.
Concluding remarks
Englacial geometry is an elusive missing piece of the glacial hydrologic system. In this study, we have investigated the limits for which fluid resonance could be used to gain insight into englacial geometries, through the analysis of wave propagation in generic coupled conduit–crack networks. We presented an overview of expected modes in the englacial system and focused on resonant modes that may be sensitive to englacial geometry. We showed that although there are many modes excited in complex englacial water transport networks, a resonant mode originating from the coupled oscillation of fluid into and out of a branching crack is the fundamental mode of the system and provides information on crack and conduit geometry. Future extensions of the modeling framework presented here could include more systematic characterization of coupled conduit and crack eigenmodes, examining boundary layer flow in the crack, adaptations to curved and cross-sectionally varying conduits, consideration of anisotropic elastic media, and modeling of seismic radiation in the ice associated with fluid motions. However, even in its present state, we believe that applications of the framework presented here to seismic and borehole water pressure studies could provide new insights on englacial transport and the englacial to subglacial connection.
Acknowledgements
The authors thank the editor and reviewers Brad Lipovsky and Dominik Gräff for their constructive comments on this paper. We thank C. Liang for writing the underlying code associated with time domain simulations presented here, and thank Eric M. Dunham, Alan Rempel, Brittany A. Erickson and Josh Crozier for discussions. L.K. acknowledges support from NSF EAR-1624431 and EAR-2036980.
Open access
