Glacial seismology

The vast majority of ice on Earth today exists in the polar regions (figures 2–6). The volume and geographic distribution of ice has varied tremendously over Earth history on time scales ranging from years to millions of years. This variation principally arises from nonlinear feedback processes occurring between the atmospheric, oceanic, and cryospheric elements of the climate system. Earth's current geographical, geological, and climatological state is amenable to dramatic cryospheric changes, and during the past few million years of Earth history the planet has undergone repeated and extensive glacial- interglacial oscillations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The glacial environment of Earth has exerted profound influence on the evolution of life, and study of glacial conditions and effects across Earth history from recent epochs to deep time is a dynamic field within Earth science, e.g. Macdonald et al (2010). To set the stage for the present glacial cryosphere, we characterize the glaciological history to the Quaternary period (composed of the Pleistocene and Holocene epochs) spanning the most recent approximately 2.5 My of Earth history. During the Pleistocene (2.588–0.0117 My before present), the volume and extent of ice on Earth underwent dramatic and repeated glacial/interglacial cycles, each typically lasting between approximately 40 to 100 kYr. These variations are extensively and widely documented by ice core temperature proxy measurements (e.g. oxygen isotope ratios), glacial landforms, marine sediments, and other paleoclimatological proxies. Aspects of glacial/interglacial cycles have been convincingly linked to insolation variations, and associated feedbacks, arising from periodic astronomically forced orbital and axial inclination changes, including eccentricity, obliquity and precession, referred to as Milankovich Cycles. However, the detailed processes controlling glaciation and glacial cycles, which represent a complex and highly nonlinear interaction between astronomical, atmospheric, oceanic, solid Earth, biological, and other forcing mechanisms, remain an active area of investigation.

The maximal extent of glacial ice during the most recent global glacial advance, referred to as last glacial maximum (LGM; e.g. Clark et al (2009)), occurred approximately 26.5–19.5 Kyr before present. At LGM, extensive ice sheets, with thicknesses of up to several kilometers, covered much of Fennoscandia and North America and, to a lesser extent, southernmost South America. Significantly greater glacial volumes than today also existed in Antarctica and across high mountain regions. Earth presently resides in the interglacial Holocene epoch (0.0117 My present). The International Commission on Stratigraphy (the entity responsible for official geologic time scale designations) is presently considering whether recent and ongoing anthropogenic changes to the planet, including widespread cryospheric melting due to increasing atmospheric and ocean temperatures primarily attributed to increasing atmospheric CO₂ and other anthropogenic changes warrant the designation of an Anthropocene epoch or age for the latest Holocene. Despite its lack of formal geological designation, the term 'Anthropocene' is currently in wide colloquial, and sometimes scholarly, use to designate the present geological time.

Seismically relevant glacial features range from millimeters (e.g. crystal scale) to thousands of kilometers (ice sheets). Large (volume greater that 1 km³) glaciers are currently present in North America, Asia, South America, Europe, and Antarctica, as well as higher latitude oceanic island landmasses (notably New Zealand and numerous sub-Antarctic and Arctic islands) (figures 2–6). The ice sheets of Antarctica (≈$13.5 \times 10^6~{\rm{km}^2}$) and Greenland (≈$1.74 \times 10^6 ~{\rm{km}^2}$) currently constitute approximately 97% of Earth's glacial area and more than 99% of its volume. Equivalent average global sea level rise for the present land-based Greenland and Antarctic ice masses are approximately 7 and 57 m, respectively (Gregory et al 2004, Fretwell 2013). Non-ice sheet glaciers are concentrated in mountainous regions of the higher-latitude southern and (especially) northern hemispheres. Although non-ice sheet glaciers constitute only a small fraction of the global glacial cryosphere, they have large ecological, climate, and hydrological roles in the natural and human resource landscape (Fountain et al 2012).

The high elevation topographic regions of the ice sheets are slowly flowing domes that transport ice via sheet flow from ice divides, e.g. Fretwell (2013). The largest glacial systems terminate at the oceans. These marine-terminating systems ultimately deliver ice seaward via outlet glaciers (topographically bounded glaciers) and ice stream (ice-flanked streams of faster moving ice) systems (Bentley 1987). In northern Canada, Greenland and Antarctica, some large marine-extending glacial systems terminate in low-topography floating ice shelves. Many other large glaciers, however, terminate more abruptly at the ocean as grounded or floating calving fronts, which may be abutted by an extensive iceberg mélange of calved ice. In both cases, the ocean can exert strong influences on terminus behavior and stability.

Glacial ice under Earth pressure and temperature conditions (0.1–35 MPa (0–4 km ice thickness) and $-80$ to $0~^\circ$C.) resides exclusively in the ice-$\mathrm {I_h}$ (hexagonal) polymorph (figure 7). The melting under pressure (regelation) of ice near $0^\circ$C has important consequences for the hydrologic and frictional characteristics of glacial beds, e.g. Weertman (1957) and Clarke (1987) (figure 7). Pure $\mathrm {I_h}$ ice at $0~^\circ$C and one atmosphere of pressure has a density of 916.79 km ${\rm m}^{-3}$, which is 91.69% that of liquid water under identical conditions (Lide 2004). The density of terrestrial in situ ice, however, varies widely based on the presence of atmospheric bubbles, rock fragments, and other impurities. The melting-point temperature of ice as a function of pressure can be approximately expressed as a Clausius–Clapeyron relationship

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle T_m = T_0 - \gamma (P - P_0) \nonumber \end{align} \tag{ 1 }$$

where T₀ and P₀ are the temperature and pressure of pure water at the triple point ($373.16~^\circ$K and $611.73$ Pa). The constant γ can range appreciably (e.g. $7.42 \times 10^{-5 \circ}$K ${\rm kPa}^{-1}$ for pure water versus $9.87 \times 10^{-5 \circ}$K ${\rm kPa}^{-1}$ for air-saturated water (Harrison 1975).

Zoom In Zoom Out Reset image size

The density of freshly fallen snow in cold regions is typically 50–70 kg ${\rm m}^{3}$. In the absence of surface or sub-surface melt, snow compacts over tens to hundreds of meters, first into compact snow (firn) and eventually into ice with increasing burial depth and pressure, eventually approaching the density of glacial ice with depth (approximately 820–920 kg ${\rm m}^{-3}$). Glacial ice is a nearly pure to highly impure polycrystalline solid composed of $\mathrm {I_h}$ crystals in varying orientations and sizes (typically from sub-mm to many cm).

At time scales exceeding seismic (elastic) periods (i.e. longer than many minutes), glacial ice viscously deforms under deviatoric stress as a strain rate-softening material. The large-strain deformation of glacial ice at these time scales may be grossly characterized by the Glen-Nye law, e.g. Glen (1955),

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \dot{\epsilon}_{ij}=\frac{1}{2 \eta (T, ~ S)} d_{ij} \label{Glen_eq} \nonumber \end{align} \tag{ 2 }$$

where η is (temperature- and stress-dependent) power law viscosity

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \eta (T, ~ S) = \frac{1}{2 C (T) S^{n-1}} \nonumber \end{align} \tag{ 3 }$$

with a rate constant $C(T)$ (dependent on both temperature and ice microstructure and impurities), and $n \approx 3$, C exhibits strong temperature dependence, for example, varying from approximately $7 \times 10^{-15}$ to $4 \times 10^{-18} ~{\rm s}^{-1}$ between 0 and $-50~^\circ$C (Cuffey and Paterson 2010).

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle d_{ij} = \sigma_{ij} - \sigma_{kk}/3 \delta_{ij} \nonumber \end{align} \tag{ 4 }$$

is the deviatoric stress tensor (the stress tensor, σ minus the mean isotropic stress), and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle S = (1/2) (d_{ij} d_{ij}) \nonumber \end{align} \tag{ 5 }$$

is the squared maximum shear stress (i.e. the second tensor invariant of d).

While (2) provides a useful formulation for fundamental transport, glacial flow is also influenced by topography, multiple intra- and inter-crystalline creep mechanisms, crystal size distribution, basal or englacial melting, water pathways and inclusions, recrystallization, and temperature profile. Resulting complex strain histories produce spatiotemporally variable, and anisotropic glacial rheology, flow, and fracture structures, e.g. Budd and Jacka (1989). In many cases internal glacial structures can be imaged in substantial detail by land-sited, airborne, or spaceborne ice penetrating radars operating in the MF, HF and VHF frequency bands, e.g. Bogorodsky et al (1985) and Kofman et al (2010).

The broad thermal structure of glaciers may be usefully categorized into warm, cold, and polythermal cases, e.g. Benn and Evans (2014). Warm glaciers occur where summer melt rates are sufficiently high. As a result, the entire ice column is close to the freezing point due to a combination of snow insulation in the winter, and melt, infiltration, refreezing, and latent heat release during summer. Cold glaciers are substantially below the freezing point from surface to base, with very limited or no surface or basal melting, and occur in very cold and dry environments, such as the Dry Valleys of Antarctica. Polythermal glaciers include elements of both warm and cold glacial ice, with their heterogeneous thermal structure reflecting surface and basal conductive and advective heat sources and transport processes (Irvine-Fynn et al 2011), as well as the influence of strain or frictional heating. The Greenland and Antarctic ice caps are broadly polythermal, and basal thermal condition there may be substantially affected by pressure melting, frictional sliding, and geothermal heat flux, while the surface may be uniformly cold, or experience seasonal melt (especially in Greenland).

2.1.1. Seismic waves in ice.

Seismic (elastic) waves in ice, and in solid media generally, are propagating, reversible (or nearly so) elastic strains at frequencies from sub-millihertz to greater than kilohertz. Except in near-source circumstances, seismic strains are usually very small ($<\sim$$10^{-6}$). In this case, seismic waves are well characterized by linear elasticity theory described by a Hooke's law constitutive relationship of the form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sigma_{ij} = c_{ijkl} \epsilon_{kl} \label{hooke} \nonumber \end{align} \tag{ 6 }$$

where σ is the (symmetric) stress tensor and is the corresponding (symmetric, and thus irrotational) strain tensor, both expressed in a common, predefined coordinate system. The stress tensor element $\sigma_{ij}$ indicates the x_j force component per unit area acting on the x_i face of an infinitessimal cube oriented along $(i, ~ j, ~ k)$ Cartesian coordinate directions. The force vector per area (traction) on a general cube face defined by a unit normal vector ${\bf {\hat{n}}}$ is thus given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tau_i = \sigma_{ij}{n_j}. \nonumber \end{align} \tag{ 7 }$$

Note that seismological convention is for extensional strain (and thus extensional stress) defined as positive. Note that this is the opposite to the convention adopted in some other (e.g. engineering) fields.

The seismic strain tensor is defined in terms of displacement, ${\bf u}$, Lagrangian spatial derivatives

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \epsilon_{ij} = (\partial u_i / \partial x_j + \partial u_j / \partial x_i)/2 . \label{strainten} \nonumber \end{align} \tag{ 8 }$$

Following general seismological practice extensional strain is defined as positive. ${\bf c}$ in (6) is the elasticity (or stiffness) tensor, with elements given by the appropriate elastic constants. Once symmetry and thermodynamic conservation considerations are taken into account, the most general physical (triclinic) $\mathbf c$ contains 21 independent constants (e.g. Aki and Richards (2002)). For isotropic media, (6) can be expressed using just two elastic constants as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sigma_{ij} = \lambda \Theta \delta_{ij} + 2 \mu \epsilon_{ij} \label{stressstrainiso} \nonumber \end{align} \tag{ 9 }$$

where λ and μ are the Lamé (elastic) parameters (where μ is the rigidity and λ approaches the bulk modulus as μ approaches zero), and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Theta = \epsilon_{ii} = \nabla \cdot {\bf u} \nonumber \end{align} \tag{ 10 }$$

is the trace of the strain tensor (the dilitation).

The hexagonal symmetry of ${\rm I_h}$ of non-randomly oriented ice crystals and/or thin (relative to seismic wavelength) layering due to deposition, flow, and/or compaction, or aligned inclusions of liquid water (Bradford et al 2013) can introduce appreciable bulk elastic rheological anisotropy within glacial ice of up to several percent, e.g. Blankenship and Bentley (1987), Diez et al (2014) and Diez and Eisen (2015) (figure 9). If the elastic properties of a bulk medium are azimuthally isotropic (i.e. anisotropic with a vertical rotational symmetry axis), then the associated elastic anisotropy is referred to as vertical transversely isotropic (VTI). The constitutive relationship for linear elasticity in a VTI solid (6) requires five distinct elastic constants. If the ice contains appreciable crystalline heterogeneity or impurities, such as englacial rock, air, or water, seismic wave propagation for wavelengths substantially longer than size scale of the bulk elastic behavior will be governed by effective (e.g. Voigt-, Reuss-, or Hill-averaged; e.g. Man and Huang (2011)) elastic moduli and density.

The existence of seismic waves in general linear elastic media can be deduced from the equation of motion (Navier's equation) for continuous media expressed in terms of stress, displacement, ${\bf u}$, and density, ρ

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial \sigma_{ij}}{\partial x_j} = \rho \frac{\partial^2 u_i}{\partial t^2} . \label{eqmo} \nonumber \end{align} \tag{ 11 }$$

Specific derivation of seismic wave equations arising from (11) requires the adoption of an appropriate constitutive relationship (6). For isotropic media, using (9), wave equations for elastic wave propagation of P (longitudinal) and S (transverse) (body) waves can be readily derived from (11) (appendix). The nondispersive propagation velocities of P and S waves in this situation are

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle v_{\rm P} = \sqrt{\frac{\lambda + 2 \mu}{\rho}} \label{wavespeed1} \nonumber \end{align} \tag{ 12 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle v_{\rm S} = \sqrt{\frac{\mu}{\rho}} . \label{wavespeed2} \nonumber \end{align} \tag{ 13 }$$

where λ and μ are the Lamé elastic moduli of (9), e.g. Stein and Wysession (1991).

The propagation velocities of P and S waves, in ice are temperature dependent, largely reflecting increasing elastic moduli at lower temperatures. Approximate linear regressions for seismic body wave speeds in isotropic compact glacial ice, compiled from laboratory and in situ measurements of seismic velocities, are (Kohnen 1974)

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle v_{\rm P} = -(2.30 \pm 0.17)T + 3795 ~ {\rm m~s^{-1}} \label{v_P} \nonumber \end{align} \tag{ 14 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle v_{\rm S} = -(1.2 \pm 0.58)T + 1915 ~ {\rm m~s^{-1}} \label{v_S} \nonumber \end{align} \tag{ 15 }$$

where T is temperature in °C (figure 8).

Zoom In Zoom Out Reset image size

These velocities are in the general range of shallow crustal rocks. This ensures appreciable elastic wave coupling between the cryospheric and other parts of the solid Earth. However, in addition to the significantly lesser density, the Poisson's ratio for glacial ice

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \nu = \frac{1}{2} \frac{(v_{\rm P}/v_{\rm S})^2 -2}{(v_{\rm p}/v_{\rm S})^2-1} \approx 0.33 \nonumber \end{align} \tag{ 16 }$$

corresponding to $v_{\rm p}/v_{\rm S} \approx 2$, is substantially higher than for typical crustal rocks, for which, typically, $\nu \approx 0.25$ and $v_{\rm p}/v_{\rm S} \approx 1.73$. Similarly, the bidirectional transmission of elastic wave energy across ice-water interfaces is also significant.

The snow-firn-glacial ice transition thickness can exceed 100 m in cold low-accumulation regions (van den Broeke 2008). This transition zone is manifested as a progressive increase in density and elastic stiffness that produces an increasing seismic velocity with depth where v_P and v_S increase from a few hundred m/s near the surface to glacial velocities (14) and (15), e.g. Albert (1998), Albert and Bentley (1998) and Diez et al (2016). This strong seismic velocity gradient creates a near-surface waveguide for shallowly dipping higher frequency (i.e. $>$1 Hz) seismic energy (Anthony et al 2015) (figure 9). For glacial systems that are strongly affected by surface melt and ablation, the transition zone will be considerably more complex, or may be absent (in the case of glacial ice exposed by ablation).

Zoom In Zoom Out Reset image size

Seismic sources associated with cryospheric processes exhibit characteristic source durations from thousands of seconds to tens of milliseconds, resulting in seismic radiation across a similar range of periods, e.g. summarized in figure 14 of Podolskiy and Walter (2016). Seismic signals are most commonly detected and studied using seismographic or strainmeter records of surface or subsurface (borehole) ground motion from seismometers installed at fixed (relative to the ice or land surface) geographic locations. Technologies and scientific analysis methods in general solid-Earth and glacial seismology are highly similar and cross-influential, with modifications for glaciological studies usually reflecting specific environmental and/and logistical considerations. As a result, the seismological and glaciological scientific communities and their geophysical instrumentation and facilities have developed increasing levels of collaboration in recent years (Aster and Simons 2015).

Seismic motion is most commonly detected using an electromagnetic inertial seismometer, which may be passive or incorporate active force feedback, damped mass-spring seismometer that produces a vertical- or three-component time series (seismogram) of ground motion. The (typically) analog seismometer output is digitized, time-stamped (usually to a GPS-synchronized clock), and recorded locally and/or telemetered. Modern portable and observatory-based broadband seismographs have useful ground motion responses that can range from many hundreds of Hz to periods of many hundreds of seconds, with a pass band amplitude response that is typically nearly flat to ground velocity within linear bounds. The instrumental noise floor is commonly designed as much as possible to lie below the natural noise floor for seismically 'quiet' sites on Earth (Peterson 1993). Seismic background noise is highly site-specific, and is dictated by a wide range of wind, ocean and other environmental processes, e.g. Dìaz (2016) that produce both quasi-continuous and transient seismic wavefield excitation.

The past two decades have produced a remarkable expansion in seismological data from the cryosphere, collected using both portable (typically seasonal to 1- to 2-year) and long term instrument deployments, such as those of the global seismographic network (GSN) (Butler et al 2004). Important community instrument pools and networks that have seen wide use in the glacial seismology (and broader seismological) research community notably include Incorporated Research Institutions for Seismology (IRIS) Consortium PASSCAL, Polar Instrumentation and UNAVCO geodetic facilities (Aster et al 2005, Nyblade et al 2012) (NSF), global seismographic network (Park et al 2005), GEOSCOPE (Roult et al 2010), GFZ Potsdam GIPP pool (Haberland and Ritter 2016), and SEIS-UK (NERC) facilities. For ground motions at periods beyond the passband of typical seismographic systems (e.g. many hundreds of seconds to thousands of seconds and beyond), larger (e.g. greater than mm-range) surface displacement fields, as well as secular glacial motions may be usefully recorded by geodetic (particularly GPS/GNSS) instruments. Integration of geodetic and seismic observations has been widely utilized to study the full bandwidth of deformation and displacement in seismogenic glaciological processes (Nyblade et al 2012).

The detection and analysis of seismic signals depends on the level of seismic background noise arising from a range of environmental processes. At periods between hundreds of seconds and several seconds, Earth's background seismic noise is dominated the oceanic microseism. The microseism consists of both body and, especially, Rayleigh surface waves, and is globally detectable at all suitably quiet seismic sites (Peterson 1993). Microseism excitation results from energy transfer from ocean gravity waves into seismic waves, principally at or near the coasts, via two distinct mechanisms. The 'primary' microseism is excited by direct ocean swell-coastal interactions, while the more powerful, frequency-doubled 'secondary' microseism is generated by standing-wave components of the ocean gravity wave field), e.g. Longuet-Higgins (1950), Douglas (1977), Ardhuin et al (2011) and Ardhuin et al (2015) that create pressure fluctuations at the ocean floor. Persistent global oceanic excitation of the elastic infrasound field in the atmosphere also occurs, e.g. Ishihara et al (2015). The global microseism noise spectrum has a power spectral density mode dictated by secondary microseism energy between approximately 5 and 10 s period, and a secondary peak or shoulder, representing primary microseism energy, at twice this period. At still shorter periods, seismic noise is commonly dominated by more local wind, water wave, and/or anthropogenic processes (Larose et al 2015). Figure 10 displays a representative probability density function of noise power spectral density for an Antarctic icecap sited station near the Amundsen-Scott Station (U.S.) near South Pole.

Zoom In Zoom Out Reset image size

Seismic signals from impulsive glaciological processes that exhibit clear P and S phases may be detected and located using standard seismological methods that have been extensively developed by the earthquake seismology community. In this case, determining source location and origin time is a well-posed nonlinear problem that may easily be solved using iterative methods provided the seismic velocity structure can be well estimated and the network of recording stations is sufficiently close to the (3-dimensional) source hypocenter. Given a sufficiently large and well-distributed seismographic network, it is feasible and common to jointly estimate both velocity structure and source locations as an inverse problem, e.g. Zhang and Thurber (2003) and Aster et al (2013) If seismographs geographically encircle the source region, but none are sufficiently close to lying directly above the source, the geographic location (epicenter) may be well determined, but source depth many be poorly constrained. If signal-to-noise ratios are high, event detection is commonly and robustly performed by envelope comparison of short-term and long-term energy averages, e.g. Withers et al (1998). However, where signal-to-noise ratios are low, or where signals are emergent or otherwise do not show clearly interpretable seismic phases due to an emergent or complex source process, or due to complex wave propagation between source and receiver, more advanced methods may be required to detect and locate a seismic source. In this case, incisive amplitude and timing information to assist with source location and other analysis may be extracted from spectrograms, e.g. Martin et al (2010a) and Bartholomaus et al (2012). Coherent arrays of seismographs (in conjunction with individual stations or as networks of arrays), when deployed, facilitate powerful beamforming methodologies that can constrain source location (or azimuth), as well as provide local measurements of apparent wave speed, wave type, and local seismic structure (Diez et al 2016), even for complex sources, e.g. Richardson et al (2012), Köhler et al (2016) and Bromirski et al (2017). For glacial seismic sources that produce a surface manifestation, camera, video, infrasound, and other ancillary instrumentation has obvious utility, e.g. Richardson et al (2010) and Bartholomaus et al (2012).

First proposed by VanWormer and Berg (1970) as a mechanism to explain seismic energy originating from glaciers, it is now clear that seismogenic stick-slip is common within the glacial environment. Basal glacial stick-slip processes have been observed on both large ice sheets and on smaller mountain glaciers across the globe. These diverse environments host asperity patches that span several orders of magnitude in scale, ranging from O(1 m) to O(100 km) (figure 12).

Where recorded with sufficient azimuthal coverage to resolve the pattern of seismic radiation (figure 11), the slip mechanism as well as the size of the slip patch can be estimated using earthquake seismological methods. Analysis of polarity or moment tensor inversion has confirmed that energy radiating for these asperities can be modeled as arising from low-angle thrust faulting (double couple mechanism), consistent with topside motion in the down-glacier direction (Anandakrishnan and Bentley 1993, Walter et al 2009). The seismic moment (A.15) can be estimated from seismograms using the long-period asymptote of the seismic spectrum, while the fall off in spectral amplitude at higher frequency can be used to constrain the slip patch size, allowing the slip, d, to be calculated from (A.15). The spectrum is often fit using

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Omega(\,f)=\frac{\Omega_0}{1+(\,f/f_c)^2} \label{brune} \nonumber \end{align} \tag{ 17 }$$

where is f_c is the spectral corner frequency and $\Omega_0$ is the long-period spectral-amplitude at frequencies less than f_c. An estimate of the radius, R, (and thus, the slip area, A) of an assumed circular asperity can be estimated using

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle R = 0.32 v_{\rm P} / f_c \nonumber \end{align} \tag{ 18 }$$

where v_P is the P wave speed for material bounding the fault (Brune 1970).

Zoom In Zoom Out Reset image size

Smaller events (magnitude 0–1), such as those observed beneath both alpine glaciers (Walter et al 2009, Allstadt and Malone 2014, Helmstetter et al 2015b) as well as ice streams (Blankenship et al 1987, Smith 2006, Winberry et al 2013), are most revealingly studied using local (typically up to a few km distant) seismographs due to their small energies. In such cases, dozens to hundreds of distinct bed asperities may be often mapped over relatively small regions (e.g. few square kilometers). Figure 12(a), shows a recent example from the Rutford Ice Stream. Slightly larger slip patches O(100–1000) m, may produce magnitude up to  ≈2 events that are readily detectable at regional distances (up to a few hundred km) (Bannister et al 2002, Danesi et al 2007, Zoet et al 2012). Figure 12(b) shows a notable example from the David glacier in Antarctica.

Zoom In Zoom Out Reset image size

The Whillans Ice Stream in Antarctica is home to largest glacial stick-slip events yet noted (Bindschadler et al 2003) (figure 12(c)). During these events, the entire  ≈100 km downstream portion of the ice stream moves in approximately twice daily 'slow earthquakes'; 30 min slip events individual seismic moment magnitudes of $M_w \approx 7$. These events were initially detected geodetically via GPS measurements of displacement collected on the glacier surface, but were later found to radiate long-period seismic energy (40–80 s) that can be observed at across Antarctica at distances of up to 1000 km (Wiens et al 2008, Winberry et al 2011). While a single event, each rupture has subsequently been shown to represent the linked breakage of several distinct large asperities ($>$5 km radius). This behavior is distinct relative to smaller events, which appear to be isolated asperities that generate several packages of energy observed in far-field seismograms (Pratt et al 2014). Additionally, the ability to jointly observe the Whillans Ice Stream slip events via GPS and seismic instrumentation deployed directly above the slip surface has revealed more detailed evolution of rupture velocity (Winberry et al 2011, Walter et al 2011, 2015b).

While subglacial stick-slip has been documented across a broad range of spatial scales, a relatively simple model of the process is supported by the observations. In particular, the generation of nearly identical waveforms associated with repeating events indicates repeated slip over isolated subglacial asperities (figure 13) (Smith 2006, Zoet et al 2012, Winberry et al 2013, Allstadt and Malone 2014). We will motivate discussion of subglacial stick-slip mechanics with a relatively simple model that captures its salient features. More detailed and complete summaries of the frictional stick-slip instabilities may be found elsewhere in the seismological literature (Scholz 1998, Marone 1998). We then discuss how the bi-material nature of the subglacial interface (i.e. ice-bedrock or ice-sediment) influences subglacial stick-slip behavior.

Zoom In Zoom Out Reset image size

As is frequently done, we will utilize a simple spring slider-block analog to discuss the stick-slip process. The spring is extended at constant rate, resulting in a linearly increasing elastic stress $\sigma_e$. If the block is initially at rest, motion will be inhibited until $\sigma_e$ exceeds the frictional strength of the block-surface contact, with frictional strength of this interface defined as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tau=\sigma^{\prime}_n c \nonumber \end{align} \tag{ 19 }$$

where $\sigma^{\prime}_n$ is effective normal stress, the applied normal stress $\sigma_n$ minus pore pressure p, and c is the coefficient of friction, and reflects the sensitivity of basal traction to subglacial hydrologic conditions. This relatively simple relationship between τ and $\sigma^{\prime}_n$ can describe basal traction beneath both soft-bedded glaciers (Iverson 2010) and debris laden basal ice sliding over hard beds (Zoet et al 2013).

However, c is not generally constant, but is instead a function of sliding speed, with the simplest case being a static value, c_s, while the block is at rest and dynamic value, c_d, once the block is moving. If friction increases at the onset of motion (velocity strengthening behavior), the block will slide smoothly. However, if friction drops at the onset of motion (velocity weakening behavior), an instability summarized in figure 14 may arise.

Zoom In Zoom Out Reset image size

First consider the case where stress levels at Point A in figure 14 initiate slip of a previously motionless block. At the onset of slip, $\sigma_e$ will begin to decrease as the spring uncompresses due to motion of the block. However, simultaneously, c will decrease as the system transitions to its dynamic value over some finite distance $\mathcal{L}$. A relatively short time after slip initiation $\sigma_e > \tau$. This force imbalance creates an acceleration of the block. The block will continue to accelerate until Point B, after which $\sigma_e < \tau$, and the block will begin to decelerate. The block will eventually come to rest at Point C, at which point c will return to its static value. Extension in spring will afterwards accumulate until $\sigma_e$ initiates another slip event. This process will produce periodic slip events with a recurrence interval of

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \delta t=2\frac{\sigma^{\prime}_n (c_s - c_d)}{\dot{\sigma_e}} \label{rieq} \nonumber \end{align} \tag{ 20 }$$

This relatively simple model has been implemented (with vary degrees of sophistication) to understand glacial stick-slip at all scales (Bindschadler et al 2003, Winberry et al 2009b, Sergienko et al 2009, Lipovsky and Dunham 2016). It is worth noting that in the glaciological literature, the term 'stick-slip' has additionally been used to describe transitions between periods of minimal and relatively fast subglacial motion that are triggered by diurnal subglacial water-pressure fluctuations (Fischer and Clarke 1997, Boulton et al 2001, Knight 2002, Damsgaard et al 2016). This hydrologic mechanism should not be confused with the rapidly-acting frictional mechanism discussed here that is the driver for seismogenic glacial stick-slip systems.

In subglacial systems, variations in pore pressure are known to drive significant variations in $\sigma^{\prime}_n$ over time scales ranging from minutes to years. Importantly for stick-slip in glacial systems, the instability of a slip patch with fixed frictional and geometrical properties will be suppressed at effective normal stresses $\sigma^{\prime}_n$ below a critical threshold given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sigma^{\prime}_{\rm crit} \leqslant \frac{K\mathcal{L}}{c_s - c_d} \label{normcrit} \nonumber \end{align} \tag{ 21 }$$

where K is the spring stiffness (Scholz 1998).

Consider the stress evolution when $\sigma^{\prime}_n=\sigma^{\prime}_{\rm crit}$ (figure 14). In this scenario, an initially motionless block begins to move in response to stress conditions represented at Point D. As with the previous case, τ will begin to drop, however, due to the value of $\sigma^{\prime}_n$, the rate of stress reduction occurs at the same rate as the reduction of $\sigma_e$, thus the force balance requires no rapid acceleration between point D and E. When point E is reached, τ remains constant and the block will slide at the rate of spring extension v. However, large perturbations to the background stressing-rate can cause a jump from stable to stick-slip behavior.

In natural systems, the slip behavior is controlled by the elasticity of the materials bounding the slip patch and frictional properties of the interface. We now discuss how the glacial setting influences each of these parameters. From the discussion above, it is clear that the stiffness of system K is fundamental to determining sliding stability. For a slip patch, K dictates the change in stress with slip and is defined as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle K=\mu/L \nonumber \end{align} \tag{ 22 }$$

where μ is the shear modulus of the material bounding the slip interface and L is the length of the slip patch. However, when asperities occur between two strongly contrasting materials μ is not well defined, such as when ice ($\mu_{\rm ice} \approx 3 \times 10^{9}$ Pa) moves over either bedrock ($\mu_{\rm rock} \approx 5 \times 10^{10}$ Pa) or sediment, ($\mu_{\rm sed} \approx 1 \times 10^{6} - 10^{9}$ Pa) (Iverson et al 1999, Lipovsky and Dunham 2016). When the material contrast is significant, elastic strain will preferentially be accommodated in the more compliant material (figure 15). Thus, for the case of a hard bed asperity, since ice is more compliant than typical bedrock by about an order of magnitude, strain will preferentially be accommodated within the ice at a bedrock contact. During a slip event under such circumstances, the strain in the ice be released due to elastic rebound during the slip. However, for an asperity atop weak sediments, the situation is less clear. The elasticity of sediment is dependent on the particular properties of the material, and its effective normal stress. However, at levels of $\sigma^{\prime}_n$ that promote sliding, $\mu_{\rm sed}$ will often be significantly lower than $\mu_{\rm ice}$ and the sediment will accommodate the majority of the inter-event strain (Iverson et al 1999, Lipovsky and Dunham 2016). For both soft and hard bed cases however, net motion only occurs in the ice.

Zoom In Zoom Out Reset image size

It is worth noting that due to the length of the central Whillans Ice Stream asperity ($>$10 km) relative to the ice thickness (≈1 km), this simple model of strain accumulation does not hold. On Whillans Ice Stream, it appears that most inter-event strain is accommodated upstream and downstream of the main sticky-spot (Winberry et al 2014).

This material dependence on strain partitioning has a significant implications for seismic estimation of seismic slip and patch size. To convert the long-period component of the seismogram spectrum into the seismic moment an estimate of the shear modulus is required. Similarly, rupture velocity v_r, which controls the duration of the slip event, will be a function of μ, influencing the estimation of slip patch size. Thus, without constraints on material properties of the asperity, authors usually bound estimates for seismic moment, patch size, slip amount with 'hard' and 'soft' bed estimates (Anandakrishnan and Bentley 1993).

Velocity-weakening friction is also a fundamental element for the above stick-slip mechanism. For hard bedded glaciers, insights emerging from limited laboratory studies are showing that the existence of velocity weakening behavior is controlled by ice temperature, basal debris content, and sliding rate (Zoet et al 2013, McCarthy et al 2017). At relatively high temperatures, near the pressure melting point, the production of meltwater reduces steady-state friction by lubricating contacts and promotes velocity strengthening behavior for both clean and debris-laden ice sliding over bedrock. Below the pressure melting point, velocity weakening behavior is favored by increased basal debris contents, while for relatively debris free basal ice, lower temperatures ($<$$-10$ $^{\circ}\rm C$) appear to promote slip instabilities. For glaciers overriding a soft bed, basal motion may be accommodated by the glacier sliding over the sediment or dominated by sediment deformation near the ice bottom ($<$10 m), thus, in many ways similar to gouge in tectonic fault systems. Velocity weakening friction is a common attribute of many granular materials, and significant efforts have been undertaken to constrain the rheological and frictional behavior of subglacial sediments. While these studies reveal that some subglacial sediments display velocity weakening behavior similar to fault gouges, most laboratory results on subglacial sediments tend to display a slight velocity strengthening (or neutral) behavior (Iverson et al 1998, Tulaczyk et al 2000, Rathbun et al 2008). These observations have prompted consideration of velocity-weakening mechanisms associated with pore-pressure reductions in the lee of basal debris as it is dragged through the underlying sediment (Thomason and Iverson 2008) and/or co-seismic pumping and pore pressure change during slip events.

In subglacial environments, resistance to motion is controlled by rheological, stress, and hydrologic conditions. Thus, the stick-slip seismicity offers a window into the temporal and spatial variability of subglacial conditions. Below, we highlight a few of the ways in which stick-slip behavior can be used to illuminate subglacial physical conditions.

3.3.1. Frictional properties of stick-slip.

3.3.2. Hydrologic modulation of seismicity.

3.3.3. Tidal and stressing rate modulation of seismicity.

It is clear from equation (20) that stressing rate variations $\dot\sigma_e$ will also modulate seismicity patterns. Tidal modulation of the force budget near the grounding lines of marine-terminating glaciers and ice sheets is the best studied driver of variable $\dot\sigma_e$. Tides modulate the glacial force budget through a combination of flexural and hydrologic perturbations that propagate upstream. The speed and magnitude with which a grounding line perturbation propagates upstream is dictated to large degree by subglacial conditions and thus, several studies have demonstrated how along flow observations of tidal influences on seismicity and velocity records can be used to infer subglacial conditions, potentially the appropriate sliding laws (Anandakrishnan and Alley 1997, Gudmundsson 2011, Walker et al 2013). However, in some regions where the surface velocity is tidally modulated, asperities appear to be unperturbed, indicating that the sticky-spots are strong relative to the size of the tidal perturbation (a few kPa) (Aalgeirsdóttir et al 2008, Smith et al 2015).

In the Ross Sea region of Antarctica, a semi-diurnal tide results in peak grounding line flow speeds on the falling tide (Anandakrishnan et al 2003) (figure 16). Analysis of $\delta t$ shows a clear modulation by the tide for both the large Whillans Ice Stream events as well as the smaller David Glacier events. Winberry et al (2009b) and Zoet et al (2012). As expected, higher stressing rates promote a decrease in $\delta t$ by promoting failure of the asperity earlier in the stick-slip cycle for the David Glacier events. On the Whillans Ice Stream, increasing flow rates near the grounding line on falling tide lead to shorter $\delta t$ (Bindschadler et al 2003, Winberry et al 2014). It is also observed in both settings that event size scales with inter-event time. This implies that c_s is not solely a function on $\sigma^{\prime}_n$, but is a time-dependent feature. This well known behavior is called 'healing' in the rock mechanics literature and is consistent with more complex friction laws (Marone 1998). Several mechanism likely control healing rates in the subglacial environment, but potential processes include increasing contact area between either side of the slip patch with time due to creep of the ice (Zoet et al 2013, McCarthy et al 2017) or pore-pressure diffusion following slip events (Iverson 2010).

Zoom In Zoom Out Reset image size

The Greenland ice sheet, surrounded by coastal mountains, is primarily drained by a small number of very large and fast moving marine-terminating outlet glaciers, some with terminus thicknesses of up to 1 km or more. Enderlin et al (2014) recently estimated that calving from only five Greenland glaciers accounted for approximately 50% (just 15 accounted for 77% and a single glacier, Jakobshavn Isbrae, drains approximately 5 to 8%; Bindshadler (1984)) of icecap mass delivered to the ocean by glacial flow ($739 \pm 29$ Gt) between 2000 and 2014.

The past decade has witnessed dramatic changes in the Greenland ice sheet due to atmospheric and oceanic forcings attributed to climate change. The proportion of ice delivered to the ocean by surface melt and sub- and en-glacial flow has been estimated by Enderlin et al (2014) to now be comparable or greater to the glacial (calving) flux, and to be increasing, (32% to 58% between 2010 and 2014) primarily in response to a climate-change driven increase in surface temperature. Another recent change to the Greenland cryosphere has been the thinning, acceleration and terminus retreat, punctuated by ice mélange and other seasonally variations (Amundson et al 2010) and of large outlet glaciers (for example, over 10 km of recent retreat at Jakobshavn Isbrae; e.g. Joughin et al (2008)).

Greenland and Antarctic outlet glaciers can generate seismic events with source moment magnitudes of up to M_w 5.1 and source durations of 60 s or longer (by comparison, the typical source duration of an M_w 5.0 tectonic earthquake is around 2 s). These events are globally observable as seismic signals with unusually strong Rayleigh waves for their magnitude between 10 and 150 s period (Ekstrom et al 2003, Ekstrom 2006, Ekstrom et al 2006, Tsai and Ekstrom 2007, Tsai et al 2008, Nettles and Ekstrom 2010, Chen et al 2011, Veitch and Nettles 2012). Early source modeling using seismographic data from global distances showed that teleseismic surface waves from these events were inconsistent with tectonic earthquake-like (double couple) source mechanisms, but could be modeled using sub-horizontally directed single-forces applied to Earth's surface in the mass-distance product range (for Greenland) of 0.1 to $2.0 \times 10^{14}$ kg m (Tsai and Ekstrom 2007). Further global monitoring established their approximate locations at the termini of large tidewater glaciers in both Greenland and Antarctica (but preferentially in Greenland) (e.g. Nettles and Ekstrom (2010)). The spatiotemporal evolution of large Greenland glacial earthquakes show seasonal and secular changes, in particular a northward advancement along the west coast of Greenland over the past two decades, that correlate with atmosphere and ocean warming, suggesting large-scale glaciological response to changing climate (Ekstrom et al 2006, Veitch and Nettles 2012).

Investigations into the source processes of Greenland glacial earthquakes, utilizing seismology, GPS positioning, satellite, and local video observations show that these events arise from glacier, water column, and ocean bed forces applied during the detachment and capsizing of very large ($\approx \!\!1 \times 10^9~{\rm m}^3$) gravitationally unstable newly calved icebergs (Amundson et al 2008, Walter et al 2012, Sergeant et al 2016, Murray et al 2015; (figures 18 and 19). Because larger events are seismically globally observable, their location and magnitudes can be used as valuable proxy measurements to monitor this type of calving. The displacement of large water volumes during these events can additionally generate significant local seiches that apply additional forces to the solid Earth and generate ocean waves that are readily detectable up to 50 km away (Amundson et al 2008). Recent seismic modeling at periods between 10 and 100 s by Sergeant et al (2016) of a particularly well studied megaiceberg calving event (which incorporated two large icebergs) at Jakobshaven showed that at least four important seismogenic mechanisms occurred during the event; collapse initiation, capsizing iceberg-to-terminus contact (distinct for the two associated icebergs), and ice mélange motion (interpreted as sub-horizontal forces generated for 20 min after calving due to collisional or shear forces within the mélange (Amundson et al 2010) and/or along the terminus margin). It appears that the complexity and variability of large glacial calving events results in highly variable seismic excitation, and recent studies incorporating video and other corroborative observations indicate that waveform inversions for moment rate functions of these events may have a limited ability to produce reliable estimates of ice loss (Walter et al 2012, Sergeant et al 2016).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Closer inspection from more locally sited seismographs has further revealed that large fast-moving Greenland glaciers can also produce long-duration (10–40 min) 'glacial rumbling' events, such as observed at Jakobshavn to occur approximately once very two days during the summer melt season by Rial et al (2009). These particular events, recorded by a seismographic array at a range near 50 km, were found to commonly culminated in a large-amplitude subevents showing distinct P and S phases that were found to be consistent with fault plane shear (i.e. a double-couple source mechanism) that was geographically associated with a bend in the feeding ice stream. The events had seismic moment magnitudes between approximately 4 and 4.3. However, the sense of slip was not able to be robustly constrained due to a lack of azimuthal source coverage by the available seismic array. Rial et al (2009) interpreted these rumbling events as progressive upstream (5–10 km from terminus) ice–ice or ice-bed stick-slip failure, typically induced by calving events at the terminus.

In addition to the major outlets of Antarctica and Greenland, large tidewater glaciers are found in other glaciated regions, including Alaska, Arctic Canada, Patagonia, and Svalbard. The stability of many of these glacial systems appears to be dominated by surface mass balance (Larsen et al 2015), but exceptions exist, including examples of surging mediated by hydrothermodynamic feedback process linked to summer melt noted in Svalbard (Schellenberger et al 2017). Calving dynamics are essential to understanding the full range of tidewater systems on all scales. Particularly well-documented in southeast Alaska, e.g. O'Neel et al (2007), it has long been recognized that calving seismicity may be a useful proxy for calving rate (Qamar 1988, Bartholomaus et al 2012, 2015b, Köhler et al 2015, Köhler et al 2016). Longer term records have been used to understand seasonal fluctuations, surge-type behavior, and dynamic changes during retreat (Walter et al 2010). The seismic and hydroacoustic source mechanisms associated iceberg calving are complex, and include shear and normal forces during detachment and gravitationally driven slip, impact forces of falling ice with the water surface, drag and buoyancy forces associated with ice deceleration, and air/water cavity collapse.

Significant recent progress in understanding seismogenic calving processes has been achieved through simultaneous video/time lapse imagery and seismic observations. Moderate sized calving events involving calved ice volumes of tens to hundreds of thousands of $\textrm m^3$ (i.e. much smaller than the megacalving events of Greenland; figures 18 and 19), and source durations of several to 1000 s), generate local to regionally detectable seismic energy between 0.5 to tens of Hz, that may be strongly peaked in the short period (e.g. 1–3 Hz; O'Neel et al (2007) and Richardson et al (2010)) band (figure 20). A common phenomenon observed to be associated with the generation of this narrowband energy is the generation of air cavity collapse-driven subvertical Worthington jets occurring several seconds after the iceberg impacts the water column (Bartholomaus et al 2012).

Zoom In Zoom Out Reset image size

A long-term goal of glacial seismology is to quantitatively estimate calving rates of marine terminating glaciers at temporal resolution that is not possible with remote-sensing observations. However, at the present state of knowledge, complexity and diversity of calving events (e.g. events with a high degree of subareal fragmentation versus those that are largely submarine) makes the inference of calved ice volumes from purely seismological observables challenging. Calibration for specific glaciers and glacial settings may substantially improve accuracy of these estimates. For example, recent work has shown that relatively simple approaches (generalized linear models) can be used to develop scaling relationships between calving observations (either visual (Bartholomaus et al 2015b) or satellite-derived (Köhler et al 2016)) and seismic attributes, in particular the duration of seismic signal associated with a discrete calving event. Once calibrated, seismic records can then be used to study calving variability over time-windows (daily, seasonal, yearly) during which other observations are not available.

Highly harmonic spectral characteristics of some iceberg-associated seismic signals were initially noted from seismic data recorded at Mount Erebus Volcano, Ross Island (MacAyeal et al 2008), in seismic signals recorded on Pacific islands (Talandier et al 2002, 2006) created by coupled hydroacoustic (T) phases (Talandier and Okal 1998, Hanson and Bowman 2006), and in seismograms recorded near the Neumayer Antarctic research station on the Princess Martha Coast of Queen Maud Land (Müller et al 2005). Müller et al (2005) noted the highly harmonic seismic signals (figure 22) associated with the local passage of drifting tabular icebergs, and hypothesized a mechanism involving forced fluid flow and excited resonances. The true mechanism of iceberg harmonic tremor was determined by MacAyeal et al (2008), who exploited a multi-year assemblage of wind and tidally driven, repeatedly colliding large tabular icebergs calved from the Ross Ice Shelf near Ross Island and McMurdo (U.S.) and Scott (New Zealand) Stations (figure 21). This situation facilitated the deployment of iceberg-sited seismographs and GPS positioning instruments to establish both relative iceberg motions and to record associated near-source seismic signals. Iceberg harmonic tremor signals were revealed to arise from repetitive stick-slip events occurring at ice–ice or ice-ground contacts (MacAyeal et al 2015) (figure 22). The harmonic features of these signals arise because impulsive subevents, sufficiently regularly spaced in time, produce harmonic spectra with a fundamental frequency that is equal to the subevent frequency (i.e. a periodic sequence of delta functions has a Fourier spectrum that is a sequence of delta functions in the frequency domain (Bracewell 2000)). The repetitive stick-slip nature of this source process was only revealed when iceberg-sited seismographs were deployed very close (several km or less) to the source zone, because multipathing and attenuation obscure the individual subevents in more distantly recorded seismograms (figure 22).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Complex seismic and hydroacoustic signals have also been found to be generated during sustained iceberg-seafloor collisions. Martin et al (2010b) analyzed iceberg-sited and remote seismographic recordings of the grounding and breakup of B15A in October 2015 as it collided with an undersea shoal near Cape Adare, Antarctica (figures 23 and 24). This relatively well-studied grounding and breakup event generated powerful harmonic and chaotic signals that were observed as far away as South Pole station (≈$17.9^\circ$ distant; Martin et al (2010b)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Ice shelves and ice tongues, the floating extensions of ice sheets and glaciers, exert fundamental control on the grounded ice through buttressing. A reduction in buttressing due to collapse or weakening of ice shelves results in the acceleration of inland grounded ice, e.g. Scambos et al (2004), thus accelerating sea level rise. The structural integrity of ice shelves is regulated by the distribution and size of crevasses and rifts (a full depth fracture on an ice shelf). Recent catastrophic ice shelf collapses have been observed to be abetted by the pooling of supra-glacial melt water that promotes hyrofracturing through preexisting fractures and subsequent gravitationally driven disintegration (MacAyeal et al 2003, Scambos et al 2009, Pollard et al 2015), and ice shelves are presently thinning globally in response to changes in ocean temperatures (Paolo et al 2015, Furst et al 2016, Marsh et al 2016, Jeong et al 2016).

Because of importance to glacial stability and sea-level rise, significant efforts are being undertaken to improve understanding of ice shelf structure, dynamics, and stability, including focused seismic studies that illuminate controls on their strain and fracture (rifting and calving) evolution. Ocean swells from distant storms create flexural elastic waves and attendant strains in ice shelves (Bromirski et al 2015) and large icebergs (MacAyeal et al 2006, 2015) that may in some cases promote fracturing (Sergienko 2010). Tsunamis from distant earthquakes have also been observed to spur the calving of large tabular icebergs from ice shelves (Brunt et al 2011). While external environmental factors appear secondary in influencing ice shelf rift propagation (Bassis et al 2007, 2008, Bassis et al 2005), near grounding lines tidal stresses have been shown to exert significant control on fracture rates (Osten-Woldenburg 1990, Barruol et al 2013, Lombardi et al 2016, Hulbe et al 2016) and microseismicity (Podolskiy et al 2016). This is significant, since many ice shelf crevasses originate at grounding lines and are then advected seaward with flow of the ice shelf, where they are increasingly influenced by the ocean.

Recent seismographic deployments atop ice shelves and tabular icebergs (Nicolas et al 1973, Okal and MacAyeal 2006, Bassis et al 2008, MacAyeal et al 2009, Bromirski and Stephen 2012, Zhan et al 2014, Bromirski et al 2015, Diez et al 2016, Bromirski et al 2017), and, in a few cases, on Arctic sea ice (Läderach and Schlindwein 2011) have observed elastic and gravity wave phenomena over periods from thousands of seconds to frequencies of many Hz. Mechanical and other perturbations to ice shelves are of significant importance to assessing the stability of these features. It has recently been discovered that ice shelf displacements also couple into acoustic-gravity atmospheric waves in the 3–10 h period range (Godin and Zabotin 2016).

Harmonic strain from ocean gravity waves may be important to ice shelf stability, e.g. Sergienko (2010), Bromirski and Stephen (2012) Oceanographic mechanical forcings include wind-driven waves and ocean swell, infragravity waves (ocean gravity waves with periods longer than 100 s), tides, ocean gravity wave signals from regional calving events MacAyeal et al (2009), and tsunamis (Brunt et al 2011). Ocean waves (Paolo et al 2015) and tides (Bromirski et al 2015) couple with the sub-ice shelf ocean water column and overlying floating shelf to produce a propagating elastic (compressional, shear, and flexural; e.g. Press et al (1951) and Sergienko (2010)) strain field within the floating glacial ice. Recent large-scale deployment of broadband seismographs across the Ross Ice Shelf has shown that coupled gravity wave/flexural wave strains can penetrate the ice shelf/water cavity system to the ice shelf grounding line (Bromirski et al 2015) (figure 25).

Zoom In Zoom Out Reset image size

Expressing the strain tensor elements (8) in terms of the displacement vector, $\mathbf u$ and expanding in terms of displacement terms lead to the general equation of motion for isotropic media, expressed in vector formulation as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \partial^2 \mathbf u / \partial t^2 &= (\lambda + 2 \mu) \nabla (\nabla \cdot {\bf u}) + \mu \nabla^2 {\bf u} \nonumber \\ &= (\lambda + 2 \mu) \nabla (\nabla \cdot {\bf u}) - \mu \nabla \times \nabla \times {\bf u} . \nonumber \end{align} \tag{ A.1 }$$

Decomposing the displacement field ${\bf u}(x, ~t)$ into the sum of the gradient of a (zero curl) scalar potential and the curl of a (zero divergence) vector potential via a Helmholtz decomposition gives

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho \frac{\partial^2}{\partial t^2} (\nabla \Phi + \nabla \times {\Psi}) =&\, (\lambda + 2 \mu) \nabla(\nabla \cdot (\nabla \Phi + \nabla \times {\Psi})) \nonumber \\ &- \mu \nabla \times (\nabla \times (\nabla \Phi + \nabla \times {\Psi})) \nonumber \end{align} \tag{ A.2 }$$

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle =(\lambda + 2 \mu) \nabla (\nabla \cdot \nabla \Phi) - \mu \nabla \times (\nabla \times \nabla \times {\Psi}) . \nonumber \end{align*}$$

We can simplify this expression by noting that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle -\nabla \times (\nabla \times \nabla \times {\Psi})& = - \nabla (\nabla \cdot (\nabla \times {\Psi})) \nonumber \\ &\quad+ \nabla^2 (\nabla \times {\Psi}) = \nabla^2 (\nabla \times {\Psi}) \nonumber \end{align} \tag{ A.3 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \nabla \cdot (\nabla \Phi) = \nabla^2 \Phi \nonumber \end{align} \tag{ A.4 }$$

so that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho \frac{\partial^2}{\partial t^2} (\nabla\Phi + \nabla \times {\Psi}) = (\lambda + 2 \mu) \nabla (\nabla^2 \Phi) + \mu \nabla^2 (\nabla \times {\Psi}) . \label{em3} \nonumber \end{align} \tag{ A.5 }$$

Grouping the potential terms in (A.5) gives

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \nabla \left((\lambda + 2 \mu) \nabla^2 \Phi - \rho \frac{\partial^2 \Phi}{\partial t^2} \right) = \nabla \times \left(\mu \nabla^2 {\Psi} - \rho \frac{\partial^2 {\Psi}} {\partial t^2} \right) . \label{sepeq} \nonumber \end{align} \tag{ A.6 }$$

The left hand side of (A.6) is the gradient of a function of Φ, while the right hand side is the curl of a function of Ψ. Setting the two bracketed terms to a constant separates the elastodynamic equation of motion into two wave equations. Choosing a constant of zero (i.e. no source terms) gives the wave equations

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \nabla^2 \Phi ({\bf x}, ~ t) = \frac{\rho}{\lambda + 2 \mu} \frac{\partial^2 \Phi ({\bf x}, ~ t)} {\partial t^2} \nonumber \end{align} \tag{ A.7 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \nabla^2 {\Psi} ({\bf x}, ~ t) = \frac{\rho}{\mu} \frac{\partial^2 {\Psi} ({\bf x}, ~ t)} {\partial t^2} . \nonumber \end{align} \tag{ A.8 }$$

that characterize the two fundamental (P and S) seismic body waves in an isotropic, homogeneous elastic medium.

Because wave fields are most often observed in the far field, plane wave modeling is frequently employed for studying seismic phenomena. A harmonic Helmholz potential plane-wave scalar displacement potential is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Phi ({\bf x}, ~ t) = \Phi_0 ({\bf x}, ~ t) = A {\rm e}^{\imath (\omega t \pm {\bf k} \cdot {\bf x})} . \nonumber \end{align} \tag{ A.9 }$$

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle =A {\rm e}^{\imath (\omega t \pm (k_x, ~ k_y, ~ k_z) \cdot (x, ~ y, ~ z))} . \nonumber \end{align*}$$

where ${\bf k}$ is the vector (propagation direction) wave number. The P wave displacement field due to the displacement potential (A.9) is thus

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tr}{{\rm tr}} \displaystyle {\bf u} = \triangledown \Phi = \pm {\rm i} A (k_x \Phi \hat{x} + k_y \Phi \hat{y} + k_z \Phi \hat{z}) \equiv B {\rm e}^{\imath (\omega t \pm {\bf k} \cdot {\bf x})} \hat{k} . \nonumber \end{align} \tag{ A.10 }$$

where $B = \pm i A \vert {\bf k} \vert$. This is a harmonic displacement disturbance with all displacement in the propagation direction $\mp \hat{k}$ and an amplitude proportional to $A \vert {\bf k} \vert$. The P wave displacement field thus consists of once-per-wavelength compressions and rarefactions. The dilatation of the P wave as a function of space and time is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \nabla \cdot {\bf u} = \mp A \vert {\bf k} \vert ^2 {\rm e}^{\imath (\omega t \mp {\bf k} \cdot {\bf x})}, \nonumber \end{align} \tag{ A.11 }$$

which does not vary perpendicular to the $\hat{k}$ direction. Volumetric strain propagating in the P wave is thus caused by shortening in the $\hat{k}$ direction.

The corresponding displacement field from a harmonic Helmholtz potential plane-wave vector displacement potential Ψ is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\bf u} = \nabla \times {\Psi} = \nabla \times (A_x, ~ A_y, ~ A_z) {\rm e}^{\imath (\omega t \pm {\bf k} \cdot {\bf x}) } \nonumber \end{align} \tag{ A.12 }$$

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle = \left(\frac{\partial \Psi_z}{\partial y} - \frac{\partial \Psi_y}{\partial z} \right) \hat{x} + \left(\frac{\partial \Psi_x}{\partial z} - \frac{\partial \Psi_z}{\partial x} \right) \hat{y} + \left(\frac{\partial \Psi_y}{\partial x} - \frac{\partial \Psi_x}{\partial y} \right) \hat{z} \nonumber \end{align*}$$

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle =&\, \pm {\rm i} \left((A_z k_y - A_y k_z) \hat{x} + (A_x k_z - A_z k_x) \hat{y}\right. \nonumber \\ &\left.+ (A_y k_x - A_x k_y) \hat{z} \right) {\rm e}^{\imath (\omega t \pm {\bf k} \cdot {\bf x})} \nonumber \end{align*}$$

The displacement in the S wave is perpendicular to the propagation direction, $\hat{k}$, and the dilatation, consistent with infinitessimal shear strain, is zero everywhere and at all times

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \nabla \cdot {\bf u} = \nabla \cdot (\nabla \times {\Psi}) = 0 \nonumber \end{align} \tag{ A.13 }$$

so that the S wave propagates solely shear strain.

A horizontally layered medium of uniform elastic layers is often usefully employed in seismology, Reflection and transmission (scattering) coefficients for like (P–P and S–S) and converted (P-S and S-P) plane seismic waves for this case are readily derivable by matching displacement and traction boundary conditions at media interfaces for ansatz plane harmonic wave solutions (obeying Snell's Law, as dictated by associated wave speeds defined by (12) and (13)). In the case of plane isotropic layers and plane waves, the seismic wavefield can decomposed into two independent wave systems. The first of these wave system is a system of S waves with particle motions that are perpendicular to the ray plane (referred to as S_H). The second wave system is the coupled system of P and S waves with particle motions lying within the ray plane (referred to as $P-S_V$). Seismic transmission and reflection coefficients can be expressed algebraically in closed forms in terms of Snell's law angles and material density and elastic parameters, and a complete set of reflection and transmission coefficient relationships for isotropic media can be found in Aki and Richards (2002). However, in elastically anisotropic media $P-S_V$ and S_H are no longer decoupled, and are thus not generally separable in this manner. Transmission and reflection coefficients for VTI media may be found in Graebner (1992), and linearized expressions for more general anisotropic media can be found in Behura and Tsvankin (2009).

An important system of coupled P and S waves is the phenomenon of seismic surface and other boundary waves associated with a free surface or other inhomogeneous interface. Such waves can be theoretically derived as coupled propagating and evanescent components of the wavefield. The interface between the solid Earth and the atmosphere can be practically considered as a free surface under most circumstances, and for horizontally layered media, surface waves can be derived using the stress-free boundary condition and corresponding reflection coefficients for a free surface, e.g. Stein and Wysession (1991). Coupled P and S surface waves polarized in the $P-S_V$ (ray) plane are generically referred to as Rayleigh waves, and surface waves comprised of S waves with particle motions perpendicular to the ray plane are generically referred to as Love waves. As with systems of non-evanescent waves, the Rayleigh (P-SV) and Love wave (S_H) components of the seismic wavefield are fully decoupled only in the case of horizontally layered isotropic media. Surface waves in layered media will exhibit dispersion, which is commonly utilized for structural inversion on all scales. Surface waves are also preferentially excited by shallower sources, and may serve as a discriminant for depth.

Seismic waves are excited when transient forces are applied at the boundaries of or internally within an elastic medium. Important sources in the cryospheric environment include examples of both types. Internal localized (so that they may be approximated by a point) sources within an elastic medium are typically characterized by a second-order seismic moment tensor, ${\bf M}$, where they corresponding scalar moment is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M_0 = \sqrt{\frac{\Sigma_{ij} M_{ij}^2}{2}}. \nonumber \end{align} \tag{ A.14 }$$

The seismic moment tensor elements $M_{ij}$ represent the relative magnitudes and signs of opposing force couples oriented in the x_i directions operating in plane specified by the normal vector x_j. Diagonal moment tensor elements thus represent opposing normal forces that are applied normal to principal planes, and off-diagonal elements represent shear couples. Because internal source process (e.g. collapse, explosion, or fault slip) cannot impart a net change in angular momentum, the moment tensor is symmetric. M₀ has units of force times distance (i.e. N-m in SI units).

A practically important moment tensor type is that of a double-couple, which describes equivalent force couples corresponding to a slipping mode II (in-plane shear) crack and is thus appropriate for slip on a fault or at a glacier bed. The scalar moment of a double couple source can be usefully calculated as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M_0 = \mu A d \label{momeq} \nonumber \end{align} \tag{ A.15 }$$

where A is the slip patch area, d is the average fault slip, and μ is the rigidity of the surrounding elastic medium. Isotropic volumetric expansion (e.g. explosion) and force couple components of a point seismic source can be represented by a diagonal moment tensor. Additionally, single (non-couple) forces, not represented in the moment tensor, can also be important elements of a seismic source characterizations, such as in modeling reaction forces related to the advection of gas or fluid, or gravitational or otherwise externally driven boundary force, such as may occur during glacial calving processes.

Seismic source sizes are commonly reported using magnitude scales, which are calculated from an appropriately scaled logarithmic formula that incorporates geometric spreading and otherwise appropriately corrected seismic amplitude or other robust observables, e.g. Shearer (2006). Magnitude estimates are typically constrained by a network of seismic stations that measure the seismic radiation over an appropriate angular extent. There are many magnitude scales in historical usage, each derived from different observables (e.g. amplitudes of body waves or surface waves within particular period bands). The most physically fundamental magnitude scale is one based on the scalar moment, referred to moment magnitude, which is now in wide general use. For a scalar moment M₀ in SI units, the corresponding moment magnitude is given by Hanks and Kanamori (1979)

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M_w = \frac{2}{3} \left(\log_{10} M_0 - 9.05 \right) . \label{mommag} \nonumber \end{align} \tag{ A.16 }$$

(A.16) provides reasonable consistency with historic magnitude scales over a wide range of sources sizes, but does not exhibit saturation or other technical problems at very large magnitudes, as can occur with other magnitude scales, e.g. Shearer (2006). For example, fault slip of 0.01 m on a fault plane measuring 100 by 100 m in ice with a medium rigidity characterized as $2 \times {10}^{10}$ Pa corresponds to a scalar seismic moment of $2 \times {10}^{12}$ N m and a corresponding moment magnitude of $M_w \approx 2.17$.

For a general double-couple seismic source representation with $\hat{n} = (n_x, ~ n_y, ~ n _z)$ specifying the normal to the fault plane and $\hat{d} = (d_x, ~ d_y, ~ d_z)$ specifying the unit fault slip vector, the moment tensor is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{dcmom} {\bf M} &= M_0 (n_i d_j+n_j d_i) \nonumber \\ &= M_0 \left(\begin{array}{@{}ccc@{}} 2 n_x d_x & n_x d_y + n_y d_x & n_x d_z + n_z d_x \nonumber \\ n_yd_x + n_x d_y & 2 n_y d_y & n_y d_z + n_z d_y \nonumber \\ n_z d_x + n_x d_z& n_z d_y + n_y d_z & 2 n_z d_z \nonumber \\ \end{array} \right) . \nonumber \end{align} \tag{ A.17 }$$

Seismic sources induce both near-field (non-propagating, static) and far-field (propagating as elastic waves and transitory) strain within an elastic medium. Near-field strain terms may become important when seismographic or other measurements of strain are recorded very close to a source region and/or for large seismic events, and may best be observed geodetically due to their static and sometimes large-amplitude components (e.g. using GPS or other positional instrumentation). In the more commonly relevant situation for glacial seismology, far-field, case, P and S wave displacement amplitudes will decay due to geometric spreading as the inverse of distance and be proportional to applied force amplitudes for single forces and proportional to the time derivative of applied couple amplitudes (or to the moment rate function in the case of extended time sources) for force couples (Aki and Richards 2002).

The orthogonal eigenvectors of (A.17) are referred to as the $\mathbf t$ (tension), $\mathbf p$ (pressure), and $\mathbf b$ (null) axes, and have eigenvalues M₀, $-M_0$, and 0, respectively. The $\mathbf t$ and $\mathbf p$ axes are the centers of the compressional and dilitational quadrants of the P wave radiation pattern from the double couple source (the maximum and minimum P wave amplitude radiation directions). Conversely, the $\mathbf b$ axis is the axis of zero radiated P- and S- wave energy, which occurs at the intersection of the fault plane and its (indistinguishable from P wave polarity information) counterpart, the auxiliary plane. The fault and auxiliary planes of a double-couple solution lie at $45^\circ$ to the $\mathbf p$ and $\mathbf t$ axes, and the intersection of the fault and auxiliary plane denotes the $\mathbf b$ axis (figure A1). Analytic expressions for seismic radiation terms for isotropic media can be found in Aki and Richards (2002).

Zoom In Zoom Out Reset image size