Skip to main content Accessibility help


  • Access
  • Open access



      • Send article to Kindle

        To send this article to your Kindle, first ensure is added to your Approved Personal Document E-mail List under your Personal Document Settings on the Manage Your Content and Devices page of your Amazon account. Then enter the ‘name’ part of your Kindle email address below. Find out more about sending to your Kindle. Find out more about sending to your Kindle.

        Note you can select to send to either the or variations. ‘’ emails are free but can only be sent to your device when it is connected to wi-fi. ‘’ emails can be delivered even when you are not connected to wi-fi, but note that service fees apply.

        Find out more about the Kindle Personal Document Service.

        Active seismic studies in valley glacier settings: strategies and limitations
        Available formats

        Send article to Dropbox

        To send this article to your Dropbox account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Dropbox.

        Active seismic studies in valley glacier settings: strategies and limitations
        Available formats

        Send article to Google Drive

        To send this article to your Google Drive account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Google Drive.

        Active seismic studies in valley glacier settings: strategies and limitations
        Available formats
Export citation


Subglacial tills play an important role in glacier dynamics but are difficult to characterize in situ. Amplitude Variation with Angle (AVA) analysis of seismic reflection data can distinguish between stiff tills and deformable tills. However, AVA analysis in mountain glacier environments can be problematic: reflections can be obscured by Rayleigh wave energy scattered from crevasses, and complex basal topography can impede the location of reflection points in 2-D acquisitions. We use a forward model to produce challenging synthetic seismic records in order to test the efficacy of AVA in crevassed and geometrically complex environments. We find that we can distinguish subglacial till types in moderately crevassed environments, where ‘moderate’ depends on crevasse spacing and orientation. The forward model serves as a planning tool, as it can predict AVA success or failure based on characteristics of the study glacier. Applying lessons from the forward model, we perform AVA on a seismic dataset collected from Taku Glacier in Southeast Alaska in March 2016. Taku Glacier is a valley glacier thought to overlay thick sediment deposits. A near-offset polarity reversal confirms that the tills are deformable.


Glaciers and ice sheets are changing rapidly due to global warming. Accurate projections of future sea level rise require improved models of glacier evolution. The response of glaciers to climate change depends strongly on the glacier/bed interface. For many glaciers, the strength and texture of subglacial sediments (tills) determine how quickly a glacier can slide and what seasonal changes can be expected in glacier behavior (Iverson and others, 1994, 1998; Walder and Fowler, 1994; Tulaczyk, 1999; Truffer and others, 2001; Hewitt, 2011). Tills are also important for glaciers that reach the ocean (tidewater glaciers) because changes in till distribution can allow tidewater glaciers to advance independently of climate, or to retreat very quickly with small amounts of warming (Alley, 1991; Post and Motyka, 1995; Motyka and Begét, 1996; Motyka and others, 2006; Brinkerhoff and others, 2017).

High water content (>40% porosity) tills are soft and deformable (Iverson and others, 1998; Iverson, 2010). They are referred to as ‘dilatant’,1 while lower porosity ($\lesssim\! 30 {\%}$) tills are stiffer, nondeforming and referred to as ‘dewatered’. Dilatant tills can contribute to glacier motion by deforming (Tulaczyk, 1999) or by facilitating glacier sliding (Iverson and others, 1995). This explains why, in ice-sheet settings, dilatant tills underlie areas of the fast-sliding ice, while dewatered tills exist under slow-moving ice (Alley and others, 1987; Anandakrishnan and others, 1998; Vaughan and others, 2003; Peters and others, 2006, 2008; Peters, 2009; Christianson and others, 2014; Luthra and others, 2016). Walter and others (2014) and Hart and others (2011) also suggest that tills can change seasonally from dilatant to dewatered states and vice-versa, due to changes in subglacial water pressure.

Observations of these subglacial sediments are important but difficult to obtain. One powerful method for determining subglacial conditions is Amplitude Variation with Angle (AVA) analysis of seismic reflection data (Sheriff and Geldart, 1995; Aki and Richards, 2002). This method has been used successfully on the Antarctic and Greenland ice sheets to distinguish dilatant tills from dewatered tills (Anandakrishnan, 2003; Peters and others, 2006, 2007; Peters, 2009; Dow and others, 2013; Christianson and others, 2014; Luthra and others, 2016). In controlled-source AVA surveys, seismic energy is recorded at a number of receivers (geophones). The records include reflections from the glacier bed as well as various types of noise.

The seismic reflectivity of the glacier bed is described by the Zoeppritz Equations (Zoeppritz, 1919; Aki and Richards, 2002), which are functions of contrasts in the density and seismic wavespeed across a reflective interface. By evaluating reflected amplitudes and plotting them against incidence angle, a series of diagnostic AVA curves (Peters and others, 2008; Booth and others, 2016) can be defined (Fig. 1). The Zoeppritz Equations assume a specular interface and that both layers are homogeneous, isotropic and have thicknesses greater than one wavelength. Though glacier ice can be anisotropic and inhomogeneous, we assume that this results in negligible errors in AVA analysis.

Fig. 1. Reflectivity curves and curve ranges for interfaces between glacier ice and various materials. Seismic parameters used to produce the ranges for till are listed in Table 1.

Table 1. Values used to produce the curves in Fig. 1 (Morgan, 1969; Hamilton, 1976; Clarke and others, 2008; Christianson and others, 2014).

AVA analysis can yield sediment compressional and shear wave velocities, as well as sediment density. These seismic parameters can be used to infer other characteristics of the glacier bed. For example, the compressional wave velocity (α) depends on sediment strength and water content. Shear wave velocity (β) also depends on sediment strength and is more sensitive to water content than α. Sediment density (ρ) is also a good indicator of water content. Table 1 shows typical seismic velocity and density ranges for various till types found beneath glacier ice.

AVA analysis requires that the primary reflection from the bed have a sufficiently high signal-to-noise ratio (SNR) to allow estimation of the polarity of the wavelet and to measure its peak or root mean square (RMS) amplitude. In addition, AVA analysis is more robust if the bed reflection multiple can be identified and measured (Peters and others, 2006, 2007, 2008; Peters, 2009; King and others, 2008; Booth and others, 2012; Christianson and others, 2014; Luthra and others, 2016). Uncertainties in bed dip and strike at the reflecting point map to uncertainties in incidence angles. Finally, high-amplitude Rayleigh waves (surface waves or colloquially ‘ground-roll’) and their scattering from surface crevasses can obscure parts of both the bed reflection and its multiple.

To avoid uncertainties in strike and dip resulting from glacier bed geometries, some studies have taken advantage of local flat spots. King and others (2008) performed reflectivity surveys at Midtre Lovenbreen, a small valley glacier in Svalbard. This survey used primary and multiple reflections from a surface-parallel planar part of the bed to determine a value for normal incidence reflectivity that was used to compare the strengths of bed arrivals elsewhere on the glacier, and revealed transitions between frozen talus and bedrock. Babcock and Bradford (2014) performed a seismic reflection survey on Bench Glacier, another valley glacier, with a steep, undulating bed. Like King and others (2008), they focused on obtaining data from a small flat area and stacked multiple wavelets to increase SNR for a full waveform inversion analysis of a thin basal ice layer.

In the absence of detectable multiple reflections, AVA analysis is possible with higher uncertainty. Dow and others (2013) and Anandakrishnan (2003) have applied AVA principles to noisy datasets by adjusting the usual AVA workflow. Dow and others (2013) worked with seismic data from Russell Glacier in West Greenland, which lacked a bed reflection multiple due to crevasse dispersion noise. They were still able to obtain a range of acceptable AVA curves by inverting for the source amplitude. Anandakrishnan (2003) analyzed subglacial sediments using a dataset with unreliable amplitudes (due to the style of seismic acquisition, which used towed snow streamers instead of geophones) from the upper part of the Kamb Ice Stream, using the incidence angle of an observed reflectivity reversal to constrain the AVA curve.

If a reflectivity reversal cannot be observed, then a wide range of glacier bed types is possible. Richards (1988) observed a temporal change from positive to negative bed reflection polarity in near-offset seismic data during the surge of Variegated Glacier, which he attributed to changes in basal sediments. This method cannot always distinguish between saturated and dewatered sediments because both sediment types can have negative reflectivities at near offsets (Fig. 1). Richards (1988) also kept the seismic line as far from crevasses as possible to minimize groundroll noise.

An additional consideration in AVA is the error associated with uncertainties in seismic wave attenuation in ice. Reflectivity calculations require a knowledge of seismic quality factor Q, which is the inverse of internal friction, a material property proportional to the fraction of energy a wave loses per cycle as it travels through a material. Low Q values correspond to higher attenuation. Q decreases as temperature and level of material fracture increase, and attenuation calculations also become more sensitive to changes in Q as Q decreases. AVA surveys on temperate, crevassed glaciers will be more vulnerable than cold ice sheets to errors in Q. Holland and Anandakrishnan (2009) suggest a modification to the AVA method that minimizes this error; however, in many cases it will not be applicable to valley glaciers as it requires a bed reflection multiple. Holland and Anandakrishnan (2009) also describe how to obtain source amplitude from the direct wave. Valley glaciers may not be good candidates for this technique because wave surface interactions could make this method unreliable if seismic lines are short.

To date, AVA has been used successfully in regions with thick ice (which separates the surface wave from the reflection); in regions with few crevasses (which significantly reduces backscatter noise); in regions where the bed and surface are relatively flat (reducing uncertainties in incidence angle estimates); and in Polar glaciers (where the cold ice has low and relatively constant attenuation with depth). In this paper, we explore the utility of the method when those conditions may not be met.

In our study, we produce synthetic seismic gathers from thin, geometrically complex, crevassed glaciers and test the usual AVA workflow, as well as the methods of Dow and others (2013) (source amplitude inversion) and Anandakrishnan (2003) (polarity crossing angle), on the synthetic data to see how well each of these methods can recover the input seismic parameters or, more broadly, the input till type. We will additionally apply these methods to a real seismic dataset we collected from Taku Glacier, a valley glacier in Southeast Alaska.


2.1. Taku glacier ice-sediment dynamics

Taku Glacier is a tidewater glacier in Southeast Alaska, currently in its advancing phase and protected from subaqueous melt and calving by a sediment shoal at its terminus (Fig. 2). The glacier experiences tidewater advance and retreat cycles that can be asynchronous with climate fluctuations and are controlled by bedrock shape and sediment dynamics and modulated by climate (Post, 1975; Brinkerhoff and others, 2017). Taku Glacier has undergone many cycles of advance and retreat during the past 3000 years (Motyka and Begét, 1996). The most recent retreat began in 1750 CE (Post and Motyka, 1995) and ended in the early to mid-1800s CE (Nolan and others, 1995).

Fig. 2. Taku Glacier terminus, showing the approximate location of the 2016 seismic survey. Historical terminus locations are shown in red (Motyka and others, 2006). Imagery is from 2010 (Google Earth).

In its advancing phase Taku Glacier is stabilized by a proglacial moraine, which shields the ice from the tides and warm water of Taku Inlet (Kuriger and others, 2006). During this phase, the glacier advances by excavating subglacial sediments and expelling them in debris flows to form this moraine (Motyka and others, 2006). As the glacier excavates sediments, its bed becomes overdeepened. Once sediment loss and/or climate change triggers a retreat, the glacier will rapidly lose mass as it calves into deeper and deeper water. In its retracted phase, Taku Glacier leaves behind a fjord (Taku Inlet) that eventually becomes filled in with outwash sediments from the glacier and fluvial sediments from the Taku River (Nolan and others, 1995).

Taku Glacier is currently advancing over fjord sediment deposits (Post and Motyka, 1995; McNeil, 2016) and offers an opportunity to study sediments under a tidewater glacier terminus. Taku Glacier experiences low strain rates compared with other Alaskan tidewater glaciers (Pfeffer and others, 2000; O'Neel and others, 2003; Truffer and others, 2009), resulting in less crevassing which allows us to perform seismic reflection surveys on its surface. We performed such a survey in March 2016 in order to obtain observations of Taku subglacial sediment properties using AVA analysis. The seismic line was located in the ablation area, ~1 km from the terminus and oriented perpendicular to glacier flow (Fig. 2).

Taku Glacier most likely overlies thick sediments in the area of our seismic line. Radar data from Motyka and others (2006) show that the glacier bed in the area of our 2016 survey lies at an elevation of ~−90 m at its deepest, 20 m below the fjord bottom mapped in 1890. Bathymetric maps (Post and Motyka, 1995) show that the deposition rate in the area of our seismic line was ~ 0.3 m a−1 from 1890 to 1937. This would extrapolate to a 1750 fjord bottom elevation of ~−110 m.

If we assume this, our seismic line is still >20 m above the 1750 fjord floor, which itself probably was not a bedrock surface. Marine seismic surveys in similar fjords show that the bedrock surface can be deeper than 300 m below sea level (Post and Motyka, 1995). If the same were true of Taku Inlet, the subglacial sediment layer would be at least 220 m thick at the location of our seismic survey.

2.2. Taku glacier sediment samples

Small sediment samples were recovered with a gravity corer in August 2015 from two boreholes at the site of the 2016 seismic survey. They consisted of sandy clays with water contents of 15–24% and 16–26%. Sample densities could not be obtained, so the porosities of these samples are unknown. Upper porosity limits of 34–45% and 35–47% are calculated based on a solid fraction density of 2600 kg m−3, a reasonable density for the local bedrock, which consists of tonalite (Gehrels and Berg, 1992).

Our calculated porosity ranges indicate a dilatant till, though uncertainties in water content do not rule out a dewatered till. The sample porosities could fall into the dewatered range if we unintentionally added water during the extraction process and/or overestimated solid fraction density. The sediment could also have experienced dewatering due to the gravity corer method (Talalay, 2013) to the effect that our range underestimates the sample porosity.

2.3. Taku glacier seismic experiment

To perform the seismic survey, we deployed 120 40 Hz geophones spaced at 5 m and buried in a vertical orientation 1 m below the surface. The geophones were installed at the southern end of the seismic line and sources were initiated along the entire 930 m long line. Sources were spaced at 10 m intervals and consisted of 120 g charges of Kinepak installed in boreholes which were backfilled with snow. Boreholes were 5 m deep (the top 1 m was snow with ice below that).

From this survey, we obtained 94 shot gathers. Seismic data processing for viewing the data and picking amplitudes was minimal and consisted of a gain applied to correct for energy loss during spherical spreading. Further data processing was performed to produce a seismic image of the glacier bed. We applied a normal moveout correction using an ice velocity of 3640 m s−1, then stacked common midpoint gathers to produce the seismic image. We derived our ice velocity from the direct wave instead of the glacier bed reflection, as the bed reflection had a lower SNR and exhibited geometrical complexity. The unmigrated image shows that the ice depth increases to the north and varies from ~160 m to 200 m. Based on these depths, the maximum source-receiver angle would be ~70°.

Figure 3 shows the shape of the glacier bed obtained from the stacked unmigrated seismic image, as well as the modeled ray tracing from one shot. This figure illustrates the effects of a rough bed on the propagation of seismic energy. Because the glacier bed is not planar, a sampling of the glacier bed is not uniform and incidence angles do not always increase with offset. Note also that the glacier bed shape in the figure is only an estimation. Our stacked common midpoint gathers may not share the same depth points and arrivals could also originate from outside the vertical plane with the seismic line. As such, it is not straightforward to define the reflection point of either a seismic reflection or its multiple.

Fig. 3. The 2016 survey geometry with ray tracing for one shot. Only rays for every 2nd receiver are shown, for clarity.

The Taku Glacier data suffer from unwanted signal interference. Ground roll contamination strongly affects bed reflections from traces up to 200 m offset from the energy source (Fig. 4).

Fig. 4. An example of a divergence-compensated, but otherwise raw seismic record from the 2016 Taku survey. Amplitudes are multiplied by travel time to correct for spherical spreading. The left lower panel shows a clear bed reflection. The right lower panel shows a signal that could be a bed reflection multiple, though wavelets are oddly-shaped due to groundroll interference.


In order to better understand the Taku Glacier seismic data, we first test AVA analysis methods on synthetic datasets. This allows us to gauge the ability of reflectivity methods to estimate subglacial velocity and density in the presence of noise. We are also able to investigate the specific effects of different noise sources and processing techniques.

The forward model uses a specified glacier bed, surface and crevasse geometry, a set of seismic parameters for the bed, an acquisition geometry similar to our field experiment and a realistic source wavelet. We assume that the glacier ice has a uniform seismic velocity, that there is no firn and that the subglacial material is uniform and thicker than 1/4 of the seismic wavelength (~10 m at the center frequency of 100 Hz), so that thin-layer effects do not distort the reflection wavelets (Widess, 1973). In reality, dilatant till layers are often thinner (Iverson and others, 1995; Porter and Murray, 2001; Evans and others, 2006; Reinardy and others, 2011; Booth and others, 2012), but a consideration of thin layer effects is beyond the scope of this paper.

We choose some typical values for dilatant till seismic parameters: ρ = 1800 kg m−3, based on a porosity of 50% and a solid fraction density of 2600 kg m−3, α = 1700 m s−1 and β = 200 m s−1, which are all within the range of observed values in Table 1.

We generate synthetic data by assuming that the seismic ray paths can be described by Snell's Law and that the seismic energy can be described by spherical spreading losses, losses due to internal friction, and losses due to reflection(s). We further assume that reflections are from specular interfaces (either crevasses or the bed). We calculate the travel path between a source and a receiver along two possible ray paths: a reflection from the bed and a surface path that interacts with (possibly many) crevasses. Reflection losses are calculated with a Zoeppritz equation Matlab script (Krebes and Margrave, 1991; Aki and Richards, 2002). We account for attenuation from internal friction by convolving the source wavelet with a frequency-dependent constant seismic quality factor (Q) impulse response (Kjartansson, 1979; Margrave, 1999).

We use a 1 m digital surface model (DSM) of a deglaciated valley as a model input. We choose the Green Lakes Valley in the Colorado Front Range, which Anderson (2014) used to illustrate a typical lumpy glacier valley. Two lakes cause flat spots to appear in the DSM; we transform these into depressions. For the ice surface we use a parabolic sheet inclined by 2° relative to the DSM, with a glacier outline defined by the intersection of the ice surface with the DSM. The DSM (760 m × 1990 m) was resampled to 10 m resolution to conserve memory and computing power. We add a chevron pattern of crevasses at the glacier edges and flow-perpendicular crevasses at the glacier midpoint. Crevasses are spaced ~15 m apart. We find from satellite images (Google Earth) that valley glaciers tend to have crevasses spaced at ~10 m to 20 m intervals.

We base our choices on forward model parameters (wavelet shapes, relative amplitudes, frequencies and anelastic attenuation values) on the 2016 Taku seismic data (Fig. 4). The following describes how we arrived at these parameter values.

3.1. The source wavelet

We design two source wavelets, one parameterized to resemble the Taku Glacier compressional wave arrivals (Fig. 5) and another to have the appearance of Rayleigh waves from the Taku Glacier dataset. Input parameters include frequency, wave damping and the number of nodes. We use Berlage wavelets (Aldridge, 1990), though other wavelet types could be substituted.

Fig. 5. The Berlage source wavelet. (a) The plain Berlage wavelet. (b) The wavelet with windowed Gaussian-random noise added; the red dashed line shows the window shape. (c) The same wavelet affected by a seismic quality factor impulse response to simulate anelastic attenuation from travel through 200 m of ice. (d) A direct arrival wavelet from the Taku Glacier dataset, recorded 200 m from the shot.

To simulate real data we add windowed white noise to the source wavelet. The noise window has zero amplitude at the start of the first arrival and ramps parabolically up to a maximum amplitude over the first wavelet half-cycle to remain constant for the next two periods. After that, its amplitude halves every two periods. We then use a highpass filter (above 50 Hz) to improve the similarity between the spectra of our real and synthetic wavelets. The resulting waveform is similar to that observed in Taku data (Fig. 5).

3.2. Seismic quality factor

Seismic quality factor Q is the inverse of internal friction, a material property proportional to the fraction of energy a wave loses per cycle as it travels through a material. We require a value for Q to calculate seismic wave attenuation in ice. Our choice for modeling Q is based roughly on observations of the Taku Glacier dataset.

We calculate the seismic quality factor values from stacked common offset gathers of the direct wave. Q is calculated independently for every offset, using the spectral ratio method, in which the frequency spectra of near and far offset direct arrivals are compared (Gusmeroli and others, 2010).

We find that the seismic quality factor varies with shot offset, first increasing and then leveling off. This reflects the fact that direct waves traveling a longer distance sample deeper glacier ice with higher Q-values. Seismic quality factor decreases with degree of material fracture unless the material is fully saturated with water. Thus, we can expect a glacier surface to have low seismic quality factor values, with Q increasing with depth (Gusmeroli and others, 2010; Babcock and Bradford, 2014) as voids in the ice become smaller and more water-saturated.

In order to find the seismic quality factor of the deeper ice using the direct wave data, we create a forward model to calculate the average seismic quality factor, Q a, in the Fresnel volume. We assume that the ice thickness is divided into a lower layer where seismic quality factor is constant (Q i) and an upper layer of thickness c = 30 m (equal to the maximum vertical extent of crevassing) where Q z (depth-dependent seismic quality factor) increases linearly from a surface value (Q s) to Q i (1). This is equivalent to a power law decrease in internal friction.

(1)$$Q_{\rm z} = \left\{ {\matrix{ {Q_{\rm s} + \left( {\displaystyle{{Q_{\rm i} - Q_{\rm s}} \over c}} \right)z} & {{\rm }0 \le z < c} \cr {Q_{\rm i}} & {{\rm }z \ge c} \cr } } \right.$$

Q is found from the average internal friction of the Fresnel zone, which is estimated by numerically integrating $Q_{ {\rm z}}^{-1}$ over the elliptical Fresnel zone of the surface wave and dividing by total Fresnel zone volume. A model with Q s = 30 and Q i = 230 provides a reasonable fit to the observed Q-values.

The Q-value of the groundroll (Q r) is also calculated following the spectral ratio method. We use a constant seismic quality factor when we calculate Rayleigh wave attenuation (Q r = 12).

3.3. Reflection raytracing

We construct bed reflections using a two-layer raytracing algorithm with a 3D layer interface. Five thousand rays are emanated from each shot location (from angle ranges −90° to +90°) and traced to the bed of the glacier and then to their emergent position on the glacier surface. The density of ray coverage is usually sufficient to provide each receiver position with at least one emergent ray, to within 10 m tolerance.

If the search returns multiple rays for a given geophone, it bins the rays by bed incident location and averages the ray arrival times so that only one arrival per bin is recorded. This model assumes that a lack of rays reaching a geophone is due to an insufficient density of modeled rays, so if no rays are returned to a geophone within the search radius, the nearest surface-incident ray is chosen. To conserve computing speed, we chose not to increase the number of modeled rays.

We then produce the reflection wavelet scaled by bed reflectivity, geometric spreading and anelastic attenuation. Since bed reflected rays interact with all layers of ice, we calculate a bulk average Q-value using (1).

A bed reflection multiple is modeled by continuing to trace the surface-incident ray to the glacier bed again and back to the surface. We record the longer travel time and transform the multiple wavelets based on Q and spherical spreading accordingly. We also scale the wavelet by the product of the reflectivities of its two reflections with the bed and its reflection with the surface. We approximate the reflectivity of the ice/air interface as −1; due to acquisition geometries, recorded multiples will have ice-air reflection incidence angles ${\vskip 10pt\lt \atop \vskip -12pt\approx}\! 40^{\circ }$ (see Fig. 1).

3.4. Reflections from surface features

The model includes backscattered signals from crevasses when it calculates direct wave and Rayleigh wave arrivals. Backscattered direct waves and surface waves are clearly visible in the Taku seismic data (see Fig. 4).

We choose a general ‘reflectivity’ c r (proportion of backscattered vs transmitted energy, converted to amplitude) of crevasses in order to produce a noise pattern, loosely based on that observed in the Taku Glacier data, that we deem realistic and desire to test. We choose a c r value of 0.3, and hold it constant for every crevasse reflection in the model regardless of incident angle or crevasse size. Note that in reality, c r depends on the ray incident angle with the crevasse and the size of the crevasse relative to the wave Fresnel zone (Benjumea and Teixidó, 2001); we ignore these considerations because our c r is just a crude approximation of highly variable values.

We determine Rayleigh and direct wave arrival times using a 2-D raytracing model, treating crevasses as planar reflectors. This is an approximation that allows us to produce different levels of surface wave noise and different arrival patterns of surface wave noise.

Ray amplitudes decrease as they are transmitted past or reflected off of crevasses. The rays that propagate past the crevasse have amplitudes equal to $\sqrt {a_{ {\rm i}}^2 - (c_{ {\rm r}} a_{ {\rm i}})^2}$, where a i is the amplitude of the incident ray. The amplitude of the reflected ray is c ra i.

Rayleigh waves and direct compressional waves also reflect off of glacier sidewalls. We use a 2-D adaptation of the 3-D raytracing algorithm to produce sidewall reflections. Rayleigh wave reflectivity is set to c r for simplicity, while direct wave reflectivity is determined via the Zoeppritz equations and the assumption that the sidewall material is identical to the basal material.

Arrivals from sidewall and surface reflections are sorted in the same way as the bed arrivals. We add scaled Berlage wavelets according to modeled arrivals times and correct the wavelets for spherical spreading and attenuation due to Q a of the wave Fresnel zone.

3.5. Seismic record assembly

Bed reflections, sidewall reflections, primary waves with crevasses and Rayleigh waves with crevasses are calculated separately. All component simulations are sampled at the same temporal sampling interval, therefore no further interpolation is required when assembling them into the full synthetic record. They are simply added together. Very low amplitude nearfield white noise is added to this record to add further realism to the model. Finally, we simulate variability in shot-geophone coupling by multiplying each trace by a factor chosen at random between 0.6 and 1.0. The coupling range is arbitrary and loosely based on observations from the Taku Glacier data.


We perform AVA analysis on the synthetic seismograms using the procedures discussed below. We then discuss our ability to retrieve the till properties prescribed in the forward model.

4.1. Incidence angle and depth

A normal moveout correction is applied to the model results, using an ice velocity calculated from the first breaks of the stacked direct wave. Common-midpoint gathers are stacked to produce a seismic image. We assume that the first breaks of the stacked section represent the bed cross section directly under the seismic line. Incidence angles and locations for each reflected wave are derived from a forward model of raypaths using this bed shape.

4.2. Shot-geophone coupling

We determine the RMS amplitudes of the direct arrival wavelets. Amplitudes are corrected for spherical spreading by multiplying by x −1, where x is the direct wave raypath length. Direct waves are also multiplied by e ax to correct for anelastic attenuation. a is found from

(2)$$a = \displaystyle{{\pi f} \over {\alpha _{{\rm ice}}Q_{\rm a}}},$$

where α ice is the compressional wave velocity of the ice and f is the dominant frequency of the wavelet. We determine Q s and Q i from the synthetic data using the same process we employed with the 2016 Taku data and find Q a for each direct wave using (1) and integrating $Q_{ {\rm z}}^{-1}$ over the Fresnel zone volume.

Attenuation-corrected direct wave amplitudes are normalized to their average and these normalization factors are assumed to correct for shot-geophone coupling variability. We multiply the amplitudes of each raw trace by its corresponding entry in the normalization vector.

4.3. The reflectivity curve

We pick the amplitudes of the bed reflection wavelets and near-offset bed reflection multiple wavelets. With these amplitudes and raylengths and with source amplitude A 0, we can calculate bed reflectivity R using:

(3)$$R = \displaystyle{{A_1} \over {A_0{\mkern 1mu} \gamma }}{\mkern 1mu} e^{ax},$$

where A 1 is the bed reflection amplitude, x is the raypath length and a is calculated from (2) using the average Q-value of the entire ice thickness and the center frequency of the bed reflection. The factor γ is a geometrical correction term,

(4)$$\gamma = \displaystyle{1 \over x}\cos \left( \theta \right)$$

where θ is the angle at which the seismic wave reaches the receiver.

We calculate A 0 using

(5)$$A_0 = \displaystyle{{A_1^2 } \over {A_2}}{\mkern 1mu} \displaystyle{{\gamma _2} \over {\gamma _1^2 }}{\mkern 1mu} {\rm e}^{2a_1x_1 - a_2x_2}$$

from Peters (2009), where A 2 is the amplitude of the bed reflection multiple. Equation (5) requires that A 1 and A 2 have similar incidence angles (we require incidence angles to be within 5° of each other). We calculate γ, a and x separately for the reflection and multiple.

This method has the potential to produce an error due to geophone coupling variability and the fact that A 1 and A 2 sample different parts of the seismic interface, unless A 1 and A 2 are zero-offset and co-trace. However, due to shot proximity and groundroll noise, we are unable to find useable zero offset, co-trace A 1 and A 2 arrival pairs in real data and all synthetic runs except the Flat run.

The reflectivity of every wavelet is calculated using (3). These reflectivities yield the AVA curve when plotted against the incidence angle. We invert for ρ, α and β using a grid search to find the best Zoeppritz curve fit in the least-squares sense (Booth and others, 2012). Grid search spacing is Δα = 20 m s−1, Δβ = 20 m s−1 and Δρ = 20 kg m−3.

The grid search is restricted to parameter combinations that are physically plausible. α, β and ρ combinations must lay within the range of a dilatant till, a dewatered till, or a consolidated till (see Table 1 for acceptable ranges).

We also test the use of a frequency bandpass filter as well as an FK filter. The bandpass filter has a lower cutoff of 60 Hz, a plateau between 120 Hz and 300 Hz, and a higher cutoff at 600 Hz. Such a filter has worked well to reduce groundroll noise in the Taku 2016 dataset, although not to the point where seismic reflection energy can be recovered at all offsets.

The FK filter is designed to remove signals with velocities less than the compressional wave velocity in ice and also includes a frequency bandpass filter with a lower cutoff of 100 Hz, a plateau between 140 Hz and 260 Hz, and an upper cutoff of 300 Hz. This is more successful than a simple bandpass filter at revealing a greater angle range of the bed reflection signal, but has the disadvantage of affecting reflection wavelet amplitudes to a greater extent.

4.4. Inverting for the source amplitude

We attempt performing AVA without the bed reflection multiple by following the methods of Dow and others (2013). For every combination of α, β and ρ, we compare the modeled Zoeppritz curve with simulated reflectivity curves calculated from the bed reflection amplitudes (binned by incidence angle) and a range of possible A 2 values. The tested range of A 2 values are equally spaced from zero to half of a reference A 1 value. This reference amplitude is equal to the normal-incidence reflection amplitude, or if that is not available, the maximum reflection amplitude. We use the range of A 2 values to calculate corresponding A 0 values using the reference reflection amplitude and (5). Next, we calculate simulated reflectivity curves from each A 0 value.

Simulated AVA curves are rejected if normal incidence reflectivity exceeds 0.6 (the maximum for any type of ice/bed interface), or if the absolute value of reflectivity for any angle exceeds 1. To allow for some data error, we add a buffer of 0.1 to both of these values. We calculate simulated curve misfits (by summing squared residuals divided by the number of datapoints, then taking the square root) and assign the smallest misfit to the grid cell for the tested α, β and ρ combination.

4.5. Crossing angle analysis

Anandakrishnan (2003) estimated seismic parameters using reflectivity crossing angle and we test his methods here. A grid search finds α, β and ρ based on the angle at which the phase reversal occurs. In source amplitude inversion, incorrect calculation of attenuation alters the curve shape and changes the results. Crossing angle analysis avoids this problem and furthermore allows us to skip Q calculation and coupling correction.

In the crossing angle inversion process, we define the misfit as the gap between the observed and calculated crossing angles. We also define the maximum acceptable misfit as either the width of an angle bin or, if angle bins adjacent to the crossing are empty, half the width of the gap plus the width of an angle bin.

4.6. Acceptable misfit

In order to characterize an acceptable range of till parameter combinations (α, β, and ρ) in AVA analysis (crossing angle analysis excluded), we calculate an envelope of acceptable Zoeppritz curve fits. We do not perform a rigorous data error analysis here, as the nature of coupling corrections and reflectivity calculations in AVA results in errors that are systematic and non-Gaussian. Instead, we want to define an error range that approximates how far up or down we can shift the best-fit AVA curve before it no longer passes between datapoints; the ‘highest’ and ‘lowest’ it can range defines the envelope.

To fall within the envelope, Zoeppritz curves must satisfy a maximum acceptable error value E max. E max is determined from the best fit Zoeppritz curve and the maximum data residual as follows.

The best-fit curve error (our minimum error E min) is equal to the sum of squares of the differences between the observed reflectivities R d and the best-fit curve R m:

(6)$$E_{{\rm min}}^2 = \sum_{i=1}^{n}\left(R_{{\rm d}}(i) - R_{{\rm m}}(i)\right)^2.$$

Here, n is the number of datapoints.

We now find the maximum residual, h, in the dataset,

(7)$$h = \max{\left(\vert R_{{\rm d}}(i) - R_{{\rm m}}(i) \vert \right)},$$

then we shift the best-fit curve up by h (approximating the top of our envelope) and re-calculate the error (6). We then define a maximum acceptable error as

(8)$$E_{{\rm max}}^2 = \sum_{i=1}^{n}\left(R_{{\rm d}}(i) - R_{{\rm m}}(i)+h\right)^2.$$

Multiplying the terms in (8), we obtain

(9)$$E_{{\rm max}}^2 = \sum_{i=1}^{n} (R_{{\rm d}}(i) - R_{{\rm m}}(i))^2 + \sum_{i=1}^{n}(R_{{\rm d}}(i) - R_{{\rm m}}(i))h + nh^2.$$

Assuming that the middle term is negligible because

(10)$$\sum_{i=1}^{n}R_{{\rm d}} \approx \sum_{i=1}^{n}R_{{\rm m}},$$

Equation 9 reduces to

(11)$$E_{{\rm max}}^2 = \sum_{i=1}^{n} (R_{{\rm d}}(i) - R_{{\rm m}}(i))^2 + nh^2 = E_{{\rm min}}^2 + nh^2.$$

We convert this to a maximum misfit value by dividing $E_{ {\rm max}}^2$ by the number of datapoints and taking the square root:

(12)$$\sigma _{{\rm max}} = \sqrt {\displaystyle{{E_{{\rm max}}^2 } \over n}} ,$$

where σ max is the maximum misfit.

We now examine the results from our grid search over α, β and ρ and designate all combinations with misfits smaller than σ max as acceptable.

4.7. Acoustic impedance

The α, β and ρ grid search results may also be expressed as an acoustic impedance (Z) vs β parameter space, where Z = αρ. For easier viewing of our results, we translate the α, β and ρ combination misfits into a Z vs β misfit plot. A problem arises at this step as numerous combinations of α and ρ produce the same Z, yet result in different AVA curves (and curve misfits), even as β is held constant. We address this problem by using the lowest misfit combination of α and ρ for each β, Z pair.

Tested Z values from our grid search are not uniformly spaced, so we perform a resampling of the Z vs β misfit function by binning misfits by Z and taking the mean misfit value, then applying a 2-D Gaussian filter to the misfit plot. Note that this causes the Z vs β plots (Fig. 6) to show best-fit β values and acceptable ranges that differ slightly from the α, β and ρ grid search results (Fig. 7 and 8). We use only α, β and ρ grid search results when reporting best-fit values and acceptable ranges for β.

Fig. 6. Parameter ranges returned by AVA analysis of model runs and Taku Glacier data, showing best fit values (dots) and acceptable ranges (whiskers). Subscripts following chart labels refer to the following AVA methods: no subscript, source amplitude inversion; A 0, source amplitude calculation; x, crossing angle analysis; fk, source amplitude inversion with an FK filter applied; xfk, crossing angle analysis with an FK filter; bandpass, source amplitude inversion with a bandpass filter.

Fig. 7. AVA results from synthetic seismic records. Tables (left) show best-fit parameter combinations and acceptable ranges, which result in the red dashed curves and gray-dashed curve envelopes in the reflectivity vs incidence angle plots (center); data reflectivities or crossing angles appear as blue dots, and the green curve is the model input. To the right are acoustic impedance (Z) vs β misfit plots; boxes labeled A, B and C encompass the dilatant till, dewatered till and consolidated till ranges, respectively. White lines mark the range of acceptable Z vs β combinations and white crosses mark the best-fit Z, β pair.

Fig. 8. AVA curve fits to reflectivities calculated from the Taku Glacier seismic record.


We perform three model runs which are distinguished from one another by input geometry (Table 2). AVA analysis is performed on these synthetic datasets (Flat, GL-long and GL-trans), as well as the real data from Taku Glacier.

Table 2. Model runs.

5.1. Deep, flat glacier run (‘Flat’)

In the Flat model, the bed reflection and the multiple are easily identified and it is possible to complete a full AVA analysis. This returns seismic parameters that are close to the input parameters and lie within dilatant till ranges (Fig. 7 Subplot A). Traditional AVA produces a misfit function with a well-defined minimum, resulting in a narrow solution range. A 0 inversion also results in a well-defined, but somewhat broader minimum which correctly establishes the till as dilatant (Fig. 7 Subplot B).

When we invert for the crossing angle only, we also obtain acceptable parameter values that fall completely within the dilatant till range. Misfit plots do not center on a local minimum, but rather a trough. The relationship between acoustic impedance Z and β is well-constrained, but there is little variation in misfit value along that line.

5.2. Longitudinal Green Lakes Valley run (‘GL-long’)

This model run uses the Green Lakes Valley geometry and a seismic line that is parallel to the glacier axis and perpendicular to crevasses (Fig. 9). The survey samples a relatively flat part of the glacier bed, keeping the raypaths simple. In the output seismic gathers (see Fig. 10 for an example), the multiple is obscured and so is part of the bed reflection. A crossing angle becomes visible when we apply an FK filter to the data; crossing angle analysis reveals the tills as dilatant (Fig. 7 Subplot D). Without an FK filter, no crossing angle is visible and source amplitude inversion yields parameters that extend into the dewatered till range (Fig. 6 and Fig. 7 Subplot C.

Fig. 9. GL-long (magenta) and GL-trans (red) survey setup. Geophones are marked as dots. Examples of raypaths (red and magenta lines) are shown emanating from shots (yellow asterisks) and reflecting off of the Green Lakes Valley surface and returning to geophones (red and magenta dots). Cyan lines represent modeled crevasses.

Fig. 10. An example of divergence-compensated, but otherwise raw seismic records from the GL-long and the GL-trans synthetic data. Amplitudes are multiplied by travel time to correct for spherical spreading. Insets show close-ups of the bed reflection and the expected location of the normal-incidence bed reflection multiple. In both seismic gathers, the ground roll signal completely obscures the multiple.

5.3. Transverse Green Lakes Valley run (‘GL-trans’)

A geometry with a transverse seismic line (parallel to crevasses, Fig. 9) produces a seismogram that looks very different from its longitudinal counterpart (Fig. 10). In GL-trans, the bed reflection is affected by overlaying hyperbolic signals where the direct wave has reflected off of crevasses. These signals resemble the bed reflection.

No multiple is visible, so we are limited to crossing angle analysis or A 0 inversion to perform AVA on this dataset. Applying an FK filter does not improve bed reflection SNR enough to see a polarity reversal, so A 0 inversion is our only course of action. Unfortunately, after A 0 inversion, possible parameter combinations exist in the dewatered till range (Fig. 6 and Fig. 7 Subplot E). Thus AVA fails to uniquely identify the till as dewatered or dilatant.

5.4. Taku Glacier 2016 survey data analysis

A signal that resembles a bed reflection multiple appears in the Taku Glacier dataset but applying AVA using this ‘multiple’ results in reflectivities greater than one. Thus we believe that this signal is a primary reflection from a large parallel crevasse or from a different part of the glacier bed and a multiple reflection is not visible.

Instead of traditional AVA, we perform source amplitude inversion and crossing angle analysis. The source amplitude inversion method is not precise enough to reveal till types in this case (Fig. 8 Subplot A), but crossing angle analysis produces acceptable parameter combinations that fall completely within the dilatant till range (Fig. 6 and Fig. 8 Subplot B).


6.1. Influence of noise

The largest impact on data quality was increased surface wave noise from crevasse backscattering. The amount of noise varied between GL-long and GL-trans due to differences in crevasse orientation producing different patterns of noise arrival times. GL-long produces a seismic record that allows AVA analysis, whereas GL-trans data fails AVA analysis. In GL-long, destructive interference resulting from modeled crevasse spacing caused decreased noise (Fig. 10). Crevasses in GL-trans are parallel to the seismic line and arrivals from them constructively interfere. Furthermore, crevasses parallel to the seismic line produce arrival times in hyperbolic patterns which resemble and could occur within a wavelength of the bed reflection, making picking difficult.

6.2. Ice thickness and AVA

Our modeling approach has allowed us to analyze geometry requirements for a successful AVA survey, assuming the seismic source and ice qualities of our Taku Glacier survey. With this model, we can predict the quality of the seismic survey before entering the field.

We observe in our model results that, at some critical ice thickness, the direct wave and/or its reflections from nearby crevasses dominates the bed reflection signal even before Rayleigh waves arrive. This critical ice thickness is dependent on compressional wavelength, the amount and duration of noise that trails the wavelet, the character of crevasse dispersion noise and the reflectivity of the glacier bed. In the best-case scenario (no crevasses, no noise trailing the direct wavelets) and with our model parameters, ice must be >46 m thick to avoid interference with the bed reflection at 60° offset, or >26 m thick for bed reflections up to 40° offset. These thicknesses will vary from survey to survey according to the shape of the direct wavelet, which depends on the frequency spectrum and impulse duration of the seismic energy source.

At this minimum thickness, modeled random noise trailing the direct wave still overlaps with the reflection signal, obscuring the shape of the bed reflection wavelet at angles with very low reflectivities. To avoid bed reflection contamination at 60° offset with the highest-amplitude part of our modeled noise tail, the ice thickness would need to be >115 m. This number varies between seismic surveys depending on the length and amplitude of noise associated with the direct arrival.

Bed reflections appear clearly in GL-trans despite the fact that ice thickness is only 70 m, because crevasse signals fortunately interfere destructively. Without this destructive interference, the ice would have to be much thicker for the bed reflection to arrive after crevasse noise is sufficiently attenuated. Such crevasse signals seen in the Taku Glacier data last up to 0.026 sec after the direct wave arrival, which results in a minimum glacier thickness up to ~180 m for a 60° offset. Within a 40° reflection angle, this minimum thickness decreases to 100 m. In order to distinguish a reflection at incidence angles past the groundroll signal (>25°) the glacier would need to be 76 m thick. This assumes the worst case crevasse spacing for the direct compressional signal wavelength (even and continuous at increments of 22-28 m). With fewer crevasses, thinner ice is acceptable.

In a scenario with many reflectors parallel to the seismic line, surpassing the minimum thickness does not guarantee a usable bed reflection because unwanted reflection hyperbolas could overlay the glacier bed signal at any depth.

6.3. Quality factor inversions

Underestimating or overestimating Q had little impact on AVA uncertainty in our experimental runs. Its effect on reflectivity calculations was superceded by other sources of error. However, Q errors could become more important in other survey scenarios, such as when Q is exceptionally low or if AVA can only be accomplished via source amplitude inversion, thus requiring an accurate reflectivity trend. Q has the potential to introduce false trends if a coupling correction is applied to the data. The effect of overestimating attenuation, for example, is that the correction vector amplifies bed returns as receiver offset increases.

To examine the importance of Q uncertainties, we re-ran AVA analysis using the model inputs for Q i and Q s (Fig. 11). The bars labeled ‘Normal’ correspond to the source amplitude inversion procedure using unfiltered data, incidence angles calculated by raytracing, and seismic quality factors calculated from the seismic line. ‘Known Q’ bars correspond to AVA results using the Q values from the forward model (Q i = 230 and Q s = 30). Bar charts show metrics for accuracy and precision. Accuracy is quantified by the percentage of acceptable parameter combinations that lie outside the range of a dilatant till (first row) and the misfit between the best-fit curve and the input curve over tested incidence angles (third row). The second row quantifies precision and shows the size of the acceptable parameter combination range as a percentage of all tested parameter combinations. Smaller values for all three metrics correspond to higher accuracy or precision.

Fig. 11. Metrics for accuracy and precision for synthetic model runs. Lower values indicate higher success. Columns show different model runs. Top row: percent of parameter combinations that lie outside of the dilatant till range. Middle row: percent of allowable parameter combinations out of all tested parameter combinations. Bottom row: the misfit between the best fit and the input AVA curves. Results are from source amplitude inversion unless indicated otherwise. Bars in the same plot represent different improvements or modifications to the AVA analysis, either by FK or bandpass filtering the data, substituting in the original Q or ϕ values, or performing crossing angle analysis (abbreviated as ‘X’) with or without a bandpass or FK filter.

Using correct Q values causes very little change in AVA results for all models using source amplitude inversion. Results were not very sensitive to Q-value calculations, even though calculated Q-values for GL-long and GL-trans differed significantly from the prescribed values (Table 3).

Table 3. Calculated seismic quality factors from model runs.

6.4. Reconstructions of incidence angles: the effects of incorrect ϕ

The GL-trans run is the most geometrically complex and produces the least-reliable raytracing results. Depth points and incidence angles calculated using the brute stack are within a median distance of 2° and 16 m (26% of the glacier thickness) of the real forward model depth point and incidence angle values, respectively. The 75th percentile incidence angle difference is smaller than our 5° angle bins, so it seems unlikely that incidence angle miscalculations could affect our results. Fig. 11 shows that GL-trans AVA accuracy (measured as misfit between input and best-fit curves) and precision (measured as mean R range or percent of parameter combinations outside of the dilatant till range) does not improve with incidence angle accuracy.

The GL-trans profile overlays ice that is only 70 m deep and yet our errors in calculated depth points range up to a few tens of meters. This would be problematic if we desired to bin bed reflections by depth point to look for spatial variability in till qualities. For GL-trans, we would need to choose a bin size of ~30% of the ice thickness or greater to reflect this uncertainty.

6.5. Effects of survey ϕ range

In this experiment, we found that crossing angle analysis is the most reliable tool for confirming a subglacial dilatant till. Crossing angle analysis is limited by the incidence angle range of observable bed reflection polarities, however. Dilatant tills have polarity reversals at low incidence angles ($\phi {\vskip 10pt\lt \atop \vskip -12pt\approx}30^{\circ }$; ϕ = 20° for our input parameters), where the bed reflection can be obscured by groundroll waves. However, a wide range of glacier bed materials, including dilatant tills, also exhibit polarity reversals at greater ϕ values (70° for our modeled dilatant till). Extending a seismic survey line past this incidence angle could increase the ability to constrain glacier bed type. A single observation of positive polarity past ϕ > 70° can rule out dewatered tills and consolidated tills as the subglacial layer, for example (Fig. 1), though it cannot distinguish between dilatant tills and bedrock. Placing a single receiver at a very far offset (${\vskip 10pt\gt \atop \vskip -12pt\approx}$6 times the ice thickness from the farthest seismic source for ϕ > 70°) could be an inexpensive and useful addition to an AVA survey with a short line.

6.6. Effects of bandpass and FK filtering

FK filtering was the only processing method that allowed us to recover input till type in the GL-long model run, as its suppression of surface wave noise allowed us to locate a polarity reversal in the bed reflection for crossing angle analysis. This polarity reversal clearly established the till as dilatant. We did not experience the same success when applying an FK filter to GL-trans, as GL-trans had more backscattered noise from compressional waves, which could not be filtered out without also removing the bed reflection.

Bandpass filtering, though less effective at removing groundroll noise than FK filtering, carries some advantages over FK filtering. An FK filter requires a close geophone spacing; the 5 m spacing of the Taku survey and the 10 m spacing of the Greenlakes surveys were sufficiently small, but the 60 m Flat survey spacing ruled out an FK filter as an option. Bandpass filtering, though unnecessary for the Flat data and leading to no advantage in accuracy or precision, could still be applied. Bandpass filtering was also applied to the GL-trans and GL-long model runs, though did not improve source amplitude inversion results and did not reveal near-offset returns. In GL-long, an FK filter did improve accuracy in the source amplitude inversion method and allowed for a crossing angle analysis, as mentioned previously. No processing methods were able to improve results in the case of GL-trans–unless noise was removed completely from the synthetic run (see the ‘Noise-free’ bar).

We found that filtering was not required for a successful crossing angle analysis in the case of Taku Glacier; though groundroll covered up many returns at near offsets, the shape of the glacier bed (having a dip in the middle of the line) ensured that a sufficient number of incidence angles were sampled to perform a crossing angle analysis–though one of the incidence angle bins closest to the reversal only included two datapoints. Bandpass and especially FK filtering were better at populating angle bins with bed polarity observations.

6.7. Shot-geophone coupling corrections

We found that our AVA results were not very sensitive to the application of shot-geophone coupling corrections. Such corrections would be more useful in a case where shot strength and geophone coupling were more variable. It is possible that errors in attenuation calculations could cause the coupling correction to introduce false amplitude trends, so the researcher should check the coupling correction vector for unrealistic offset-dependent trends.


This study affirms the notion that deeper and less-crevassed glaciers with simple bed geometries are the best candidates for AVA studies (Fig. 12). At the other end of the spectrum, very thin (<100 m) heavily-crevassed glaciers cannot yield data suitable for AVA analysis.

Fig. 12. Conceptualization of AVA survey quality based on ice thickness and degree of crevasse noise. Each blue dot marks a reported reflectivity survey. Red dots are from modeled surveys. Degree of crevasse noise is from remarks made by the author, or we estimate it from photographs or satellite imagery of the studied glacier.

Glaciers in the middle of the spectrum–without workable bed reflection multiples, but with bed reflections that are not overwhelmed by crevasse noise–can be successful AVA candidates. A near-offset crossing angle confirms a dilatant till, as long as crossing angle can be recovered after FK or bandpass filtering. Source amplitude inversion is a reasonable method if source amplitude is uncertain and it is not possible to identify glacier bed type using crossing angle analysis. However, for source amplitude inversion to work, bed reflection SNR must be high enough to reveal an amplitude trend. In the case of a high-quality bed reflection but an uncertain A 0, source amplitude inversion may be as accurate as AVA. It avoids the possibility of unquantified errors in the source amplitude calculation skewing results.

This forward model serves as a planning tool for seismic surveys on thin, crevassed, geometrically-complex glaciers. Crevasse locations can be obtained from satellite imagery of the glacier surface. The glacier bed can be constrained using existing radar data or, barring that, an estimated glacier bed that combines a typical glacier width/depth ratio with the maximum amount of basal topography that can be reasonably expected. If AVA is successful from the simulated data despite the high amplitude of basal topography, then it should be successful with the study glacier. The model could be run with the hypothesized till parameters and perhaps repeated with the local bedrock to ensure the two can be distinguishable from each other. Model results can also indicate the best use of resources in the seismic survey. If the bed reflection is obscured at near offsets due to crevasse noise, shots should be located farther away from the geophone line. However, this will make it impossible to use a direct calculation of the source amplitude.

Our synthetic survey results show that geometry uncertainties do not significantly impact AVA results. Because of this, we are confident that the Taku Glacier AVA results are not misleading, even with only an estimate of the glacier bed shape. Our brute stack did not show a basal topography that was as rough as the GL-trans topography, so our ϕ calculation errors could not have exceeded the GL-trans errors. It is probable that the Taku Glacier incidence angle calculation errors were smaller than the size of our angle bins, so the topography beneath the Taku Glacier seismic line could not have caused a till type misidentification.

Crossing angle analysis shows that Taku Glacier sediments are within the range we have defined as dilatant and as a consequence, are deformable. A deformable till under Taku Glacier has consequences for its terminus evolution. Till deformation allows faster evacuation of sediments from beneath Taku Glacier, as deforming till creeps towards subglacial channels, where it can be eroded fluvially–an important mechanism for sediment evacuation (Motyka and others, 2006) and tidewater glacier dynamics (Brinkerhoff and others, 2017).


JZ performed most of the writing of this paper, most of the model development, led the seismic survey portion of the fieldwork performed on Taku Glacier in 2016, and performed all analysis of the 2016 seismic dataset with the advice and supervision of MT. MT, AG, JA and CL conceived of the study and participated in the field work. AB proposed and produced a preliminary version of the forward model and advised on AVA processing strategies. AG helped with seismic data processing. All authors contributed to the writing of the manuscript.


The Taku 2016 survey was made possible by input and assistance from Thomas Hart (blaster), Aurora Roth, and Andy Aschwanden. IRIS PASSCAL provided seismic recording equipment for the Taku 2016 seismic survey. Thanks to Peter Burkett, Sridhar Anandakrishnan and Kiya Riverman for instrument support, and to Bernard Coakley, Bernard Hallet, Esther Babcock and some anonymous reviewers for their input and advice. Also thanks to the Consortium for Research in Elastic Wave Exploration Seismology (CREWES) for providing Matlab scripts that facilitated the seismic forward modeling. This project is funded by National Science Foundation grant number 130 4899 and supported in part by a University of Alaska Fairbanks Center for Global Change Student Research Grant with funds from the Cooperative Institute for Alaska Research.

1 Dilatant till is defined herein as ‘simply [a till] that has expanded by enlarging pore spaces with water’ (Clark, 2010).


Aki, K and Richards, PG (2002) Quantitative Seismology, 2nd Edn. San Francisco: University Science Books
Aldridge, GF (1990) Short note: the Berlage wavelet. Geophysics, 55(11), 15081511
Alley, RB (1991) Sedimentary processes may cause fluctuations of tidewater glaciers. Ann. Glaciol., 15, 119124
Alley, RB, Blankenship, DD and Bentley, CR (1987) Till beneath Ice Stream B: 3. Till deformation: evidence and implications. J. Geophys. Res., 92(B9), 89218929
Anandakrishnan, S (2003) Dilatant till layer near the onset of streaming flow of Ice Stream C, West Antarctica, determined by AVO (amplitude vs offset) analysis. Ann. Glaciol., 36, 283286
Anandakrishnan, S, Blankenship, DD, Alley, RB and Stoffa, PL (1998) Influence of subglacial geology on the position of a West Antarctic ice steam from seismic observations. Nature, 394(September), 6265
Anderson, RS (2014) Evolution of lumpy glacial landscapes. Geology, (August), 42(8), 679682
Babcock, E and Bradford, J (2014) Quantifying the basal conditions of a mountain glacier using a targeted full-waveform inversion: Bench Glacier, Alaska, USA. J. Glaciol., 60(224), 12211231
Benjumea, B and Teixidó, T (2001) Seismic reflection constraints on the glacial dynamics of Johnsons Glacier, Antarctica. J. Appl. Geophys., 46(1), 3144
Booth, AD and 6 others (2012) Thin-layer effects in glaciological seismic amplitude-versus-angle (AVA) analysis: Implications for characterising a subglacial till unit, Russell Glacier, West Greenland. Cryosphere, 6, 909922
Booth, AD, Emir, E and Diez, A (2016) Approximations to seismic AVA responses: validity and potential in glaciological applications. Geophysics, 81(1), WA1WA11
Brinkerhoff, DJ, Truffer, M and Aschwanden, A (2017) Sediment transport drives tidewater glacier periodicity. Nature Comms., 8(90), 18 DOI 10.1038/s41467-017-00095-5
Christianson, K and 7 others (2014) Dilatant till facilitates ice-stream flow in northeast Greenland. Earth Plan. Sc. Lett., 401, 5769
Clark, CD (2010) Emergent drumlins and their clones: from till dilatancy to flow instabilities. J. Glaciol., 51(200), 10111025
Clarke, BG, Hughes, DB and Hashemi, S (2008) Physical characteristics of subglacial tills. Geotechnique, 58(1), 6776
Dow, CF and 5 others (2013) Seismic evidence of mechanically weak sediments underlying Russell Glacier, West Greenland. Ann. Glaciol., 54(64), 135141
Evans, DJA, Phillips, ER, Hiemstra, JF and Auton, CA (2006) Subglacial till: Formation, sedimentary characteristics and classification. Earth-Sc. Rev., 78(1–2), 115176
Gehrels, GE and Berg, HC (1992) Geologic map of Southeastern Alaska, U.S Geological Survey map I-1867
Gusmeroli, A and 5 others (2010) Instruments and Methods: Seismic wave attenuation in the uppermost glacier ice of Storglaciären, Sweden. J. Glaciol., 56(196), 249256
Hamilton, EL (1976) Shear-wave velocity versus depth in marine sediments: a review. Geophysics, 41(5), 985996
Hart, JK, Rose, KC and Martinez, K (2011) Subglacial till behaviour derived from in situ wireless multi-sensor subglacial probes: Rheology, hydro-mechanical interactions and till formation. Quat. Sc. Rev., 30(1–2), 234247
Hewitt, IJ (2011) Modelling distributed and channelized subglacial drainage: the spacing of channels. J. Glaciol., 57(202), 302314
Holland, CW and Anandakrishnan, S (2009) Subglacial seismic reflection strategies when source amplitude and medium attenuation are poorly known. J. Glaciol., 55(193), 931937
Iverson, NR (2010) Shear resistance and continuity of subglacial till: Hydrology rules. J. Glaciol., 56(200), 11041114
Iverson, NR, Jansson, P and Hooke, RL (1994) In-situ measurement of the strength of deforming subglacial till. J. Glaciol., 40(136), 497503
Iverson, NR, Hanson, B, Hooke, R and Jansson, P (1995) Flow mechanism of glaciers on soft beds. Science, 267(5194), 8081
Iverson, NR, Hooyer, TS and Baker, RW (1998) Ring-shear studies of till deformation: Coulomb-plastic behaviour and distributed strain in glacier beds. J. Glaciol., 44, 634642
King, EC, Smith, AM, Murray, T and Stuart, GW (2008) Glacier-bed characteristics of midtre Lovénbreen, Svalbard, from high-resolution seismic and radar surveying. J. Glaciol., 54(184), 145156
Kjartansson, E (1979) Constant Q-wave propagation and attenuation. J. Geophys. Res., 84, 47374748
Krebes, E and Margrave, GF (1991) Zoeppritz [Computer software], CREWES Toolbox Version 1272, syntraces, Available at
Kuriger, EM, Truffer, M, Motyka, RJ and Bucki, AK (2006) Episodic reactivation of large-scale push moraines in front of the advancing Taku Glacier, Alaska. J. Geophys. Res.: Earth Surface, 111, 113 F01009, DOI 10.1029/2005JF000385
Luthra, T, Anandakrishnan, S, Winberry, JP, Alley, RB and Holschuh, N (2016) Basal characteristics of the main sticky spot on the ice plain of Whillans Ice Stream, Antarctica. Earth Plan. Sc. Lett., 440, 1219
Margrave, GF (1999) Einar [Computer software], CREWES Toolbox Version 1272, Qtools,Available at
McNeil, C (2016) Surface mass balanche reanalysis of Taku and Lemon Creek Glaciers, Alaska: 1946–2015. (PhD thesis)
Morgan, NA (1969) Physical properties of marine sediments as related to seismic velocities. Geophysics, 34(4), 529545
Motyka, RJ and Begét, JE (1996) Taku Glacier, southeast Alaska, U.S.A.: Late Holocene history of a tidewater glacier. Arctic Alpine Res., 28, 4251
Motyka, RJ, Truffer, M, Kuriger, EM and Bucki, AK (2006) Rapid erosion of soft sediments by tidewater glacier advance: Taku Glacier, Alaska, USA. Geophys. Res. Lett., 33(24), 15
Nolan, M, Motkya, RJ, Echelmeyer, K and Trabant, DC (1995) Ice-thickness measurements of Taku Glacier, Alaska, USA, and their relevance to its recent behavior. J. Glaciol., 54(184), 145156
O'Neel, SO, Echelmeyer, KA and Motyka, RJ (2003) Short-term variations in calving of a tidewater glacier : LeConte Glacier, Alaska, U.S.A. J. Glaciol., 41(139), 541553
Peters, LE (2009) A seismic investigation of basal conditions in glaciated regions. (PhD thesis, Pennsylvania State University)
Peters, LE and 6 others (2006) Subglacial sediments as a control on the onset and location of two Siple Coast ice streams, West Antarctica. J. Geophys. Res.: Solid Earth, 111, 114
Peters, LE, Anandakrishnan, S, Alley, RB and Smith, AM (2007) Extensive storage of basal meltwater in the onset region of a major West Antarctic ice stream. Geology, 35, 251254
Peters, LE and 5 others (2008) Seismic detection of a subglacial lake near the South Pole, Antarctica. Geophys. Res. Lett., 35(L23501), 18 DOI 10.1029/2008GL035704
Pfeffer, WT, Cohn, J, Meier, M and Krimmel, RM (2000) Alaskan Glacier beats a dramatic retreat. Eos Trans. AGU, 81(48), 577
Porter, PR and Murray, T (2001) Mechanical and hydraulic properties of till beneath Bakaninbreen, Svalbard. J. Glaciol., 47(157), 167175
Post, A (1975) Preliminary hydrology and historic terminal changes of Columbia Glacier, Alaska, U.S. Geolological Survey Hydrologic Investigations Atlas, (map HA-559)
Post, A and Motyka, RJ (1995) Taku and LeConte Glaciers, Alaska: Calving-Speed of Late-Holocene Asynchronous advances and Retreats. Phys. Geogr., 16, 5982
Reinardy, BTI and 5 others (2011) Streaming flow of an Antarctic Peninsula palaeo-ice stream, both by basal sliding and deformation of substrate. J. Glaciol., 57(204), 596608
Richards, MA (1988) Seismic evidence for a weak basal layer during the 1982 surge of Variegated Glacier, Alaska, USA. J. Glaciol., 34(116), 111120
Sheriff, RE and Geldart, LP (1995) Exploration Seismology, 2nd Edn. New York: Cambridge University Press
Talalay, PG (2013) Subglacial till and bedrock drilling. Cold Reg. Sc. Tech., 86, 142166
Truffer, M, Echelmeyer, KA and Harrison, WD (2001) Implications of till deformation on glacier dynamics. J. Glaciol., 47(156), 123134
Truffer, M, Motyka, RJ, Hekkers, M, Howat, IM and King, MA (2009) Terminus dynamics at an advancing glacier: Taku Glacier, Alaska. J. Glaciol., 55(194), 10521060
Tulaczyk, S (1999) Ice sliding over weak, fine-grained tills: dependence of ice-till interactions on till granulometry. Geol. Soc. Am. Spec. Pap., 337, 159177
Vaughan, DG, Smith, AM, Nath, PC and Le Meur, E (2003) Acoustic impedence and basal shear stress beneath four Antarctic ice streams. Ann. Glaciol., 36, 225232
Walder, JS and Fowler, AC (1994) Channelized subglacial drainage over a deformable bed. J. Glaciol., 40(134), 315
Walter, F, Chaput, J and Lüthi, MP (2014) Thick sediments beneath Greenland's ablation zone and their potential role in future ice sheet dynamics. Geology, 42(6), 487490
Widess, MB (1973) How thin is a thin bed? Geophysics, 38(6), 11761180
Zoeppritz, K (1919) VII b. Über Reflexion und Durchgang seismischer Wellen durch Unstetigkeitsflächen. Nachrichten von der Gesellschaft der Wissenschaften zu Göttingen, Mathematisch-Physikalische Klasse, 1919, 6684