An Assessment of Geophysical Survey Techniques for Characterising the Subsurface Around Glacier Margins, and Recommendations for Future Applications

Geophysical surveys provide an efficient and non-invasive means of studying subsurface conditions in numerous sedimentary settings. In this study, we explore the application of three geophysical methods to a proglacial environment, namely ground penetrating radar (GPR), seismic refraction and multi-channel analysis of surface waves (MASW). We apply these geophysical methods to three glacial landforms with contrasting morphologies and sedimentary characteristics, and we use the various responses to assess the applicability and limitations of each method for these proglacial targets. Our analysis shows that GPR and seismic (refraction and MASW) techniques can provide spatially extensive information on the internal architecture and composition of moraines, but careful survey designs are required to optimise data quality in these geologically complex environments. Based on our findings, we define a number of recommendations and a potential workflow to guide future geophysical investigations in analogous settings. We recommend the initial use of GPR in future studies of proglacial environments to inform (a) seismic survey design and (b) the selection of seismic interpretation techniques. We show the benefits of using multiple GPR antenna frequencies (e.g., 25 and 100 MHz) to provide decimetre scale imaging in the near surface (e.g., < 15 m) while also enabling signal penetration to targets at up to ∼40 m depth (e.g., bedrock). This strategy helps to circumvent changes in radar signal penetration resulting from variations in substrate conductivity or abundant scatterers. Our study also demonstrates the importance of combining multiple geophysical methods together with ground-truthing through sedimentological observations to reduce ambiguity in interpretations. Implementing our recommendations will improve geophysical survey practice in the field of glacial geology and allow geophysical methods to play an increasing role in the interpretation of glacial landforms and sediments.

Introduction

Near-surface geophysics offers considerable potential for non-invasive imaging of the subsurface within Earth Science, but the use of geophysical methods in terrestrial (onshore) proglacial environments remains relatively limited. Where these methods have been applied in glacial geological studies, they have mostly used ground penetrating radar (GPR) (e.g., Spagnolo et al., 2014; Tonkin et al., 2017; Lovell et al., 2019; Fitzsimons and Howarth, 2020; Woodard et al., 2020) and few have applied and compared or combined multiple geophysical surveying approaches (e.g., Parkes et al., 2013; Dobiński et al., 2017; Kunz and Kneisel, 2020). However, an integrated, multi-method geophysical surveying approach can provide a range of complementary but distinct data (e.g., on sedimentary architecture, sediment thickness, depth to bedrock) for characterising proglacial environments.

Detailed knowledge of both the subsurface and surficial geomorphology in glacial environments is required to provide crucial empirical data on past glacier behaviour over different timescales. Critically, such evidence is used to address outstanding questions surrounding the current and future response of glaciers to rapid climate change, by providing context for current observations (e.g., Hannesdóttir et al., 2015; Chandler et al., 2016a; Chandler et al., 2016b; Åkesson et al., 2017; Weber et al., 2019). Frequently, glacial-geological interpretations of past glacier dynamics rely on either a “geomorphology only” approach, where the interpretations are based purely on surface morphology (e.g., Beedle et al., 2009), or sedimentological data from a small number of glacial landforms (e.g., Krüger 1995; Lukas 2007). Subsurface glacial-geological investigations are necessarily selective due to the often-limited availability of exposures through glacial landforms and sedimentary sequences, the impracticality of manually creating exposures through often large glacial landforms, and restrictions on excavations (e.g., in environmentally protected areas). Geophysics could alleviate these issues by imaging extensive glacial landform and sediment suites, thereby providing considerably more spatially-extensive data (e.g., Busby and Merritt, 1999; Jakobsen and Overgaard, 2002; Hauck and Kneisel, 2008; Lukas and Sass, 2011). This can establish a broader geological context for highly localised sedimentological datasets, aiding studies of past glacier behaviour (Schrott and Sass, 2008; Van Dam, 2012).

Although geophysical applications in glacial geology have increased in number during the past 2 decades (e.g., Bennett et al., 2004; Hauck and Kneisel, 2008; Stokes et al., 2008; Spagnolo et al., 2014), there is currently a lack of literature explicitly focusing on the use of multiple geophysical techniques at glacial landforms within active forelands (e.g., Lecomte et al., 2008; McClymont et al., 2011). The limited application of near-surface geophysics to contemporary proglacial environments potentially reflects the complexity and diversity of sediments and landforms in these dynamic environments (e.g., Evans, 2003; Benn and Evans, 2010; Hambrey and Glasser, 2012), which makes geophysical investigations challenging (due to the non-uniqueness of geophysical interpretations). The performance and applicability of different geophysical surveying methods and associated parameters and processing steps are, however, difficult to assess due to the limited number of previous studies and the lack of a systematic assessment. This knowledge gap will be addressed in this paper.

We present the application and outcomes of multiple geophysical surveying techniques in an active proglacial environment and develop a framework for best practice in future applications. To do this, we test and assess the performance of three geophysical surveying methods: GPR, seismic refraction, and multichannel analysis of surface waves (MASW). This testing is undertaken at three contrasting sites in the dynamic proglacial environment of Midtdalsbreen (Norway), which contains a range of glacial landform and sediment types (e.g., Andersen and Sollid 1971; Sollid and Bjørkenes 1978; Reinardy et al., 2013; Weber et al., 2019). In the following section, we provide the motivation behind the site selection, survey designs, and parameter choices, along with detailed descriptions of the procedures undertaken in the investigation at Midtdalsbreen. We then present the results of the geophysical surveys at each of the individual survey sites, before synthesising the collective geophysical datasets. These results are used to inform recommendations for future geophysical investigations in proglacial, and other analogous environments.

Methodology

Field Site Selection

To assess the performance of GPR, seismic refraction, and MASW for characterising proglacial landforms and environments, we aimed to apply the methods in an active proglacial environment that (a) contains a range of glacial landforms and sediment types and (b) is easily accessible (e.g., by snowmobile) with the various geophysical equipment. Based on these criteria, we selected the foreland of Midtdalsbreen, an outlet glacier of the Hardangerjøkulen icefield in Norway (60^o34′N, 7^o27′E; Figure 1). Midtdalsbreen has a dynamic ice-marginal to proglacial environment that contains a variety of glacial sediment-landform associations, including recessional (“annual”) moraines, controlled moraine, de-iced hummocky moraine, flutes, variable till cover (thin/discontinuous to thick/continuous), glaciofluvial sediments, and localised areas of supraglacial debris and boulder spreads (e.g., Andersen and Sollid, 1971; Sollid and Bjørkenes, 1978; Reinardy et al., 2013, 2019; Weber et al., 2019). The glacier snout overlies a subsurface variably composed of wet sediment (∼10 m thick), fractured bedrock, and permafrost (Killingbeck et al., 2020), with bedrock comprised of granitic gneiss, phyllites, and schists (Murray and Dowdeswell, 1992; Willis et al., 2011; Åkesson et al., 2017). The foreland of Midtdalsbreen is easily accessible from the Finse Alpine Research Centre (either via snowmobile or on foot in snow-free conditions), which was ideal for our test study.

FIGURE 1

FIGURE 1. (A) Study area location within Norway, red marker, (B) Satellite image of Hardangerjøkulen with the location of C outlined in red, (C) Satellite image of the study site showing the locations of the geophysical survey lines and sediment excavations (E1–E7) and the approximate position of the 2018 glacier margin (dashed line). In this paper, GPR data are used to exemplify the characteristic outcomes of the techniques at the contrasting sites so only a subset of them is presented. LIA = Little Ice Age. (D) Photograph of section E1, facing West, (E) Photograph of test pit E3, facing South, (F) Photograph of section E4, facing East.

For the geophysical surveys, our goal was to test the performance of the three methods for surveying proglacial landforms by applying them to sites with (a) contrasting internal sedimentary and structural properties (e.g., coarse, cobble-rich material vs. better sorted, fine-grained sediments) and (b) contrasting surface morphologies (e.g., simple linear ridge vs. more complex hummocky topography). To achieve this, we selected three sites in the Midtdalsbreen foreland (Figure 1C) after undertaking initial geomorphological and sedimentological observations:

1. Hummocky Site: An area of hummocks that consist of glaciotectonised sandy sediments with brittle and ductile deformation structures (Figure 1D). The site is located to the north west of the glacier, on the ice proximal side of the prominent Little Ice Age moraine (Figure 1C). The hummocks are ∼2–8 m wide and up to 4 m in height, with occasional surface boulders up to 2 m long. Bedrock is exposed between the hummocks.

2. Lateral Moraine Site: A broad lateral moraine ridge with sand, gravel and widely spaced boulders at its surface. The lateral moraine is situated close to the present-day western ice margin (Figure 1C), and it is oriented parallel to the dominant ice flow direction (from NNE to SSW). The moraine ridge is ∼10 m wide and has a relatively flat crest ∼2 m above the surrounding foreland.

3. Terminal Moraine Site: A pronounced terminal moraine ridge with an open blockwork of medium sized (10–80 cm) cobbles and boulders over diamicton at its surface (Figure 1F). This prominent moraine was formed in 2001 and can be traced across the central foreland (Figure 1C; Reinardy et al., 2013). The moraine is up to ∼5 m high and ∼10 m wide, and it has an asymmetric form with a steeper (∼30°) ice-distal slope. A meltwater stream dissects the moraine ∼100 m from its western end, providing exposure E4 (Figure 1C).

The geophysical and sedimentological data presented here were collected across four field campaigns in 2018 and 2019. Site identification took place in September 2018, when the foreland was snow-free, and was guided by geomorphological mapping of the Midtdalsbreen foreland (Sollid and Bjørkenes, 1978; Reinardy et al., 2013). Geophysical data were acquired when the foreland was snow-covered, allowing transport of equipment by snowmobile. All surveys were undertaken in spring 2019, except for the seismic survey at the Terminal Moraine Site, which took place in April 2018. Sites were revisited in June 2019 to ground truth the geophysical interpretation.

Ground Penetrating Radar

GPR is widely used for shallow environmental investigations, offering high resolution imaging with relatively portable apparatus (Jol and Bristow 2003). It uses pulses of high frequency (commonly 10–2,500 MHz) electromagnetic (EM) waves to image subsurface boundaries between materials with contrasting EM properties whereby the EM waves are recorded when they return to the surface. The reflectivity of a boundary is defined by contrasts in the propagation velocity of the GPR wavelet either side of that boundary, itself influenced by contrasts in the bulk dielectric permittivity. In the proglacial environment, changes in dielectric permittivity and velocity are caused by variations in sediment grain size, water saturation, ice content, or the presence of bedrock (Supplementary Table S1). Depth of sampling is influenced by electrical conductivity, with high conductivities limiting sample depth because GPR energy is absorbed. This can be overcome to some extent by using antennas of lower centre-frequency, since high-frequency wavelet components are preferentially absorbed.

Acquisition

We undertook GPR common-offset surveys using a MALÅ ProEx system and unshielded Rough Terrain Antennas (RTA) with nominal centre-frequencies of 100 and 25 MHz at all three sites (Figure 1C). The bistatic RTA system has transmitting and receiving antennas orientated in-line, separated by 2.2 and 6.2 m for the 100 and 25 MHz systems, respectively. For a target at 10 m depth, antennas with ∼100 MHz centre-frequencies may ordinarily be used. However, a lower-frequency antenna was also used in this study, due to the heterogeneity and electrical conductivity of substrates in proglacial environments, which can lead to scattering of high frequency signals. To simultaneously record the geographical location and elevation of the GPR data, we deployed a GPS receiver (Hemisphere A101) with ±2.5 m horizontal accuracy, mounted on a backpack, 4.9 m (100 MHz) and 8.6 m (25 MHz) ahead of the RTA midpoints.

At the Hummocky Site, we conducted surveys on a 100 m by 70 m grid, with a line spacing of 10 m. Additional profiles were taken over prominent features, including along the seismic survey line and over the location of the logged sedimentary section (E1; Figure 1C). At the Lateral Moraine and Terminal Moraine sites, we surveyed a combination of long-profiles (along the moraine crestlines) and cross-profiles (perpendicular to the moraine crestlines), with cross-profiles surveyed at 15 and 30 m line spacings for the 100 and 25 MHz systems, respectively (Figure 1C). Additionally, we collected profiles parallel to the crestline at ∼4 m from the ice-proximal slope and ∼6 m from the ice-distal slope. All surveys were conducted by towing the RTAs on foot, at a pace of ∼1 m/s, as exposed surface boulders and snow-free moraine crests were impassable with a snowmobile.

Processing

We processed the GPR data using ReflexW (version 8.5.6; Sandmeier, 2016) and Matlab (version R2018a; MATLAB 2018) software. Initial processing followed a common sequence (Supplementary Figure S1):

1. dewow filter using 30 and 80 ns windows for the 100 and 25 MHz data, respectively;

2. static corrections, to shift first-breaks to the expected arrival time of the air-wave (7.3 ns for the 100 MHz antennas, 20.7 ns for the 25 MHz);

3. automatic gain control (AGC), using a window length of 50 ns for both frequencies; and

4. bandpass filtering, using trapezoidal (Ormsby) filter, with corner frequencies of 5–20–240–480 MHz for 100 MHz data, and 5–10–40–80 MHz for 25 MHz data).

To account for variable trace spacings resulting from variations in tow speed, trace spacings were then regularised (to 0.25 and 0.50 m, respectively, for the 100 and 25 MHz data) using a linear interpolation in Matlab and the recorded GPS coordinates. The process also corrected the horizontal offset between the RTA midpoint and the GPS antenna, which is located ahead of the RTA midpoint (see Acquisition). Thereafter, spatial filtering was undertaken to suppress horizontal noise in the data (e.g., direct waves and antenna ringing) using various “subtracting average” filters. All filter types and gain functions were selected based on qualitative examination of the outcomes after extensive parameter testing.

Migration was undertaken using a Kirchhoff algorithm, requiring an estimate of the subsurface radar velocity. Common midpoint surveying was not possible as the RTA antennas are inseparable; thus, velocities were defined using a migration velocity analysis routine (Yilmaz, 2001; Neal, 2004), with the optimal migration velocity chosen based on the degree of image focusing. The Kirchhoff time migration was conducted in Matlab, using an algorithm that honours topographic variation (Allroggen et al., 2014) and prevents incorrect focusing of energy as a result of changes to surface slope (Supplementary Figure S1F). The final velocities were also used for depth estimation. These depth estimates are likely to be accurate to ±20% or less, due to uncertainties in the velocity estimates, which could be high for the heterogeneous sedimentary sequences typically found in proglacial environments (cf. Carrick Utsi, 2017).

Following Neal (2004), we interpreted the processed radar data to identify prominent reflections and classify radargram characteristics in terms of reflector shape, dip, continuity, and configuration. These interpretations were checked for consistency, albeit with the change in resolution, between the 100 and 25 MHz datasets. The vertical resolution of a GPR antenna under ideal conditions is considered to be a quarter of the received wavelength (Sheriff and Geldart, 1983) and wavelength is equal to the propagation velocity divided by the signal frequency. Here, using a propagation velocity of 0.1 m/ns, we calculated the vertical resolution for the 100 and 25 MHz wavelets to be approximately 0.25 and 1 m, respectively (Davis and Annan, 1989).

Seismic Refraction

Seismic refraction methods determine both the thickness and seismic velocity of target layers by transmitting seismic energy (usually compressional (P-) waves) into the ground and subsequently recording the return signal with a spread of geophones (Steeples 2005). The seismic refraction method uses the arrival time of the first onset (“first-break”) of seismic energy at each geophone (Reynolds, 2011). By plotting arrival time against source-receiver offset, the simplest assessments of seismic velocity use the reciprocal of the slope of first-break times as a velocity estimate. This is justified where horizontal, planar interfaces between discrete layers are present and at each, the velocity of P-waves (v_p) increases. In glacial environments, strong velocity contrasts (Supplementary Table S2) include hard bedrock underlying sediment (e.g., Sass 2006; McClymont et al., 2011) or the transition from unfrozen to frozen ground (e.g., Hauck et al., 2004), making them suitable for detection with seismic refraction methods.

Acquisition

Our seismic surveys were planned to ensure that the data would be appropriate for refraction and MASW processing: We acquired seismic data using a Geometrics GEODE system, with 48 vertical-component geophones and a 6 kg sledgehammer-and-plate source. Previous studies have shown that a sledgehammer source provides enough energy for surveying targets at 10–30 m depth (Schrott and Sass, 2008) and this option produces a good low frequency (<30 Hz) signal, which is required for MASW (Rossi et al., 2018). To ensure that both refractions and surface waves were captured, each record was 1 s long, sampled at 0.125 ms and logged using geophones with low natural frequencies (4.5 and 10 Hz; Ivanov et al., 2008).

As a general rule, the geophone array should be at least five times as long as the depth to the deepest target refractor (e.g., bedrock) such that refractions from that target can be observed as first arrivals in the seismic dataset (Redpath, 1973). Sediment thicknesses across the foreland were expected to be up to 10 m, based on previous work by Reinardy et al. (2013), Reinardy et al. (2019) and Killingbeck et al. (2019). Thus, the seismic surveys were designed with geophone array lengths of ∼50 m or more. During 2018, surveys at the Terminal Moraine Site used geophones (denoted G 1–G 48) with a natural frequency of 10 Hz, planted at 1 m intervals. Seismic shots were taken at each geophone location (denoted S_n, where n is the distance from G 1) and 1 m away from each end of the line, forming an array length of 49 m (Supplementary Figure S2). In some cases, geophone planting was difficult and inconsistent due to exposed, snow-free moraine crests: geophones G 17–G 22 were partially exposed, and up to 0.3 m laterally offset from the survey line, hence greater uncertainty was assigned to their first-break picks during analysis.

Prior to the 2019 seismic data acquisition, we modelled the subsurface for both sites using expected depths to interfaces and substrate velocities. Based on knowledge of sediment types in the foreland (Reinardy et al., 2013; Reinardy et al., 2019), initial results from the 2018 seismic data, and typical seismic velocities (Supplementary Table S2), we predicted substrate v_p values of 1,000–1,800 m/s (partially saturated sediments), 2,000–3,500 m/s (frozen sediments) and 4,400–6,400 m/s (bedrock) at each of the sites. We then used these modelled values to design the 2019 seismic surveys. Our survey strategy allowed for depth sampling to at least 15 m, while also ensuring close geophone spacing to minimise the risk of missing arrivals from thin layers (Lankston, 1990) and to conform to topographic restrictions.

The 2019 surveys at the Hummocky and Lateral Moraine sites (Figure 1C) used 48 vertical component geophones with 4.5 Hz natural frequency (lower than in 2018, to benefit MASW analysis), installed at 1.5 m intervals and a spread-length of 70.5 m. Line locations were selected to ensure that all geophones would be planted in packed snow (having examined signal-to-noise ratios in the 2018 survey). We took offset shots 1.5 and 30 m from each end of the spread to extend the depth of investigation and provide reverse cover of the deepest target refractor (Reynolds, 2011). To sample direct waves for the calculation of velocity in the uppermost unit (v₁), we installed short-offset spreads of 24 geophones at 0.5 m spacing at each end of the survey line. We recorded at least ten ‘clean’ hammer strikes (i.e., that made a firm contact with the impact plate, with minimal plate movement) at each shot location, for later stacking and improvement of signal-to-noise ratio (Processing). The repeat shots also compacted fresh snow beneath the impact plate, improving coupling and the consistency of the source waveform.

Processing

Seismic refraction processing and interpretation was undertaken using ReflexW and Matlab. We first removed trigger errors, resulting from disparities between the recording time initiation and the hammer blow by considering the zero-offset time of the air wave, found by linear regression on arrivals at other offsets, as the actual shot instant (Figure 2B), and adjusting arrival times to reflect this. Individual shots were then stacked (between 5 and 15 records at each location) to boost signal-to-noise ratio, with noisy records and those with inconsistent source waveforms discarded. We undertook manual first-break picking in ReflexW (e.g., Figures 2A,B), with picks interpreted as two-layer (Terminal Moraine Site) and three-layer (Lateral Moraine Site) cases. The slope of straight-line segments fitted to picks of first arrival times against shot-receiver distance defined the reciprocal of the P-wave velocity, v_p, of each layer (Supplementary Figure S3). By evaluating slope using the approach defined in York et al. (2004), velocity estimates incorporated the uncertainty in both time picking and geophone location.

FIGURE 2

FIGURE 2. (A) First 70 ms of stacked seismic shot gather acquired at the Lateral Moraine Site, with the shot location S −1.5. First-break picks are marked with blue dashes. (B) First 230 ms of reverse seismic shot gather at the Lateral Moraine Site with shot location S 72. Sections are displayed in positive standard polarity. First-break picks are marked with orange dashes. Grey areas highlight surface wave arrivals and geophone offsets used in MASW analysis are defined with “MASW S” and “MASW N” for the southern and northern ends of the line, respectively.

The thickness of the uppermost layer (H₁) or depth to the first refractor (Z₁) at the Terminal and Lateral Moraine Sites were calculated using

$\begin{array}{c}{H}_1 = \frac{t_{i1}{v}_1{v}_2}{2\sqrt{\left(v_2^2 - v_1^2\right)}} \end{array}(1)$

where t_i1 is the intercept time of the second linear segment, v₁ is v_p in layer 1 and v₂ is v_p in layer 2 (Supplementary Figure S3). At the Lateral Moraine Site, the thickness of the second layer (H₂) was calculated using

$\begin{array}{c}{H}_2=\left[t_{i2}- \frac{2H_1\cos \left(i_{c\mathrm{1,3}}\right)}{v_1}\right]\frac{v_2}{2\cos \left(i_{c\mathrm{2,3}}\right)}\end{array}(2)$

where i_c1,3 = sin(v₁/v₃); i_c2,3 = sin(v₂/v₃); and t_i2 is the intercept time of the third linear segment (Palmer 1986). H₁ and H₂ were then summed to give the overall depth to the second refractor (Z₂).

The variation in the first-break pick slope at the Hummocky Site suggested significant subsurface heterogeneity, incompatible with this analysis, so seismic refraction interpretation was not performed on those data. Even for the other sites, systematic variations in the arrival time of the deepest refraction imply undulating boundaries, suggesting that our preliminary analysis could be enhanced via (e.g.,) the plus-minus or TimeDepth methods (Hagedoorn 1959; Hawkins 1961), which use delay-times to calculate the depth to the refractor beneath each geophone. Implementation of such methods is beyond the scope of this paper, but examples of delay-time application in glacial and formerly glaciated environments can be found in Hausmann et al. (2007) and Ross et al. (2019).

Multichannel Analysis of Surface Waves

Although most seismic surveys focus on the propagation of body waves to interpret subsurface properties, MASW uses the vertical component of surface waves, known as Rayleigh waves or ‘groundroll’. Rayleigh waves travel along the near surface in an elliptical motion where the depth of penetration is determined by its wavelength (λ) (Biot, 1962). For any given frequency (f) the wavelength is

$\begin{array}{c}\lambda =\frac{v_{phase}\left(f\right)}{f}\end{array}(3)$

where v_phase is the phase velocity of a given frequency component (Stokoe et al., 1994). The relationship between f and v_phase is called dispersion. This dispersive behaviour is used in MASW inversions to produce 1-D profiles of the depth variation of shear (S) wave velocity (v_s), beneath the geophone spread (Supplementary Figure S4; Aki, 1957; Park et al., 1999), albeit with a further assumption of lateral homogeneity beneath the spread, which we address in Processing.

Acquisition

MASW in this study is based on the same datasets collected for seismic refraction (see Acquisition). The depth of investigation (Z) provided by the surface waves is dependent on their wavelength ($\lambda$). Depending on subsurface conditions, Z can be given by $\lambda /2$ (e.g., Park et al., 1999) or $\lambda /3$ (Gazetas, 1982). The wavelengths surveyed are controlled by survey geometry:

- the geophone spacing (d_x) controls the minimum horizontal wavelength (λ_min) recorded,

$\begin{array}{c}{\lambda }_{min} = 2d_x\end{array}(4)$

and hence, the minimum resolvable depth of investigation (Z_min), and.

- the geophone spread length (D) is directly related to the maximum horizontal wavelength (λ_max) recorded,

$\begin{array}{c}{\lambda }_{max}=D\end{array}(5)$

hence D controls the maximum resolvable depth of investigation (Z_max) (Park et al., 1999).

To account for subsurface heterogeneities, we estimated Z_min and Z_max using $\lambda /3$ (Gazetas, 1982), suggesting that spread lengths between 10 and 70 m would provide sensitivity to a 1–20 m depth range.

Processing

We performed MASW on the processed seismic shots using the open-source software Geopsy (version 2.9.1.; Wathelet, 2011). Dispersion curves were obtained by transforming the seismic data from the time-offset domain to the frequency-phase velocity domain using a linear frequency-wavenumber (f-k) transform (March and Bailey, 1983). Differing dispersion patterns from the near and far offset geophones, at the Hummocky and Lateral Moraine sites, indicated a change in subsurface structure along the lines. To avoid violating the assumption of lateral homogeneity, we processed near offset and far offset traces separately to create independent dispersion curves for each end of the seismic lines (offsets in Table 1). The range of geophones used at each end was guided by the characteristics of the seismic shot gather and features in the co-located GPR data. Dispersion images were calculated using a progressively smaller geophone offset range, with the optimal range representing a trade-off between noise in the image and smearing of laterally-variable structure.

TABLE 1

TABLE 1. MASW dispersion curve inputs and resolutions for the three sites.

The dispersion curves were defined by picking the maximum response for each frequency (Supplementary Figure S4). In this study, dispersion curve picking was restricted to the fundamental mode only: the higher modes are not dominant and cannot be reliably picked. The frequency ranges of the dispersion curves, clearly identified above the background noise level, and their associated Z_min and Z_max, are given in Table 1 for each of the selected geophone ranges.

We inverted the picked dispersion curves in the Geopsy Dinver module to produce 1-D v_s profiles. The Dinver module uses the inversion process developed by Wathelet et al. (2004), in which Sambridge’s (1999) Neighbourhood Algorithm is implemented. The initial input model includes constraints applied to the Earth properties v_p, v_s, density (ρ), and unit thicknesses, for the expected layers to Z_max (Xia et al., 1999; Killingbeck et al., 2018). The inversion process then seeks to minimise the difference between the modelled and observed dispersion curves (Park et al., 1999), leading to the definition of a v_s depth model and its uncertainty. We constrained model ranges for v_p using the results from seismic refraction and contact depths using the GPR data. For each picked dispersion curve, we ran the inversion process with at least 10 different initial models, ranging from 3 to 10 subsurface layers, to assess the non-uniqueness of the inversion. Each run of the inversion process tested at least 90,000 models. All the calculated v_s profiles and dispersion curves were displayed with a colour representative of their misfit value to highlight the models with the closest fit to the observed data.

Sedimentological Observations

To ground truth the geophysical data, we undertook sedimentological investigations in selected locations. These sedimentological observations supplement more extensive information on the sedimentology of the Midtdalsbreen foreland, which has been published elsewhere (Reinardy et al., 2013; Reinardy et al., 2019). Further information on sediment properties (e.g., grain size distribution and sediment strength) could be determined by collecting sediment samples for lab testing; however, this was beyond the scope of this study.

We manually excavated exposures (E1–E7; Figures 1C–F) at each site to identify dominant sediment types and measure depths to sedimentary boundaries. Where possible, we created sediment sections perpendicular to moraine ridges to enable identification and description of sedimentary structures. Exposure E1 (Figure 1D) was oriented north-west-south-east (325⁰–145⁰), perpendicular to the long axis of a hummock, enabling us to examine its internal sediment structure and texture. At the Lateral Moraine Site, where natural cross-sectional exposures were unavailable, we dug two pits along the moraine crest, reaching ∼0.8 m in depth (Figure 1E). The Terminal Moraine Site is dissected by a meltwater stream, which exposes 2 m of sediment in the moraine (Exposure E4, Figure 1F) ∼1 m below the moraine crest. We supplemented observations of E4 with excavations of three pits along the crestline of the terminal moraine, reaching depths of ∼0.5 m.

Results

Lateral Moraine Site

In test pits E2 and E3 (Figures 1C,E), the upper 0.8 m consists of highly saturated, sandy silt with occasional gravel-sized clasts (1–10 cm long). High water saturation led to continual collapse of the pit walls, meaning it was not possible to dig below 0.8 m. A sand layer, possibly partially frozen or highly consolidated, was recognised at the base of both pits. E2 at the bottom (North) end of the latero-terminal moraine, had a higher concentration of gravel near the surface compared to E3 and was less saturated.

GPR velocities at the Lateral Moraine Site range from 0.07 to 0.13 m/ns, with the highest velocities observed at the northern end of Profile 1. Depth conversions are performed with the representative velocity of 0.10 m/ns, although their accuracy is likely no better than ±20%. In Profiles 1 and 2, continuous reflectors are apparent between 3 and 5 m depth in the 100 MHz radargrams (Figures 3A,B,E,F) and 8–14 m depth in the 25 MHz radargrams (Figures 3C,D,G,H). Both reflectors have a similar geometry, inclining steeply towards the surface just north of a relict meltwater channel (Figures 3B–D), which we observed during summer fieldwork. This channel drained snowmelt, leading to very wet ground conditions at the time of the survey, which may explain the attenuated signal below. There are clear changes in attenuation, apparent for both antenna frequencies at coincident locations. On the cross-profiles (e.g., Profile 2a, b; Figures 3E−H), attenuation is higher on the ice-proximal side of the moraine, with a distinct change directly below the moraine crest. On the profiles running transverse to the moraine crestline (e.g., Profile 1a, b; Figures 3A−D), high attenuation is observed at the northern end of the moraine, where the moraine flattens out at its lowest elevation. In the centre of Profiles 1a and b (Figures 3A–C), low attenuation allows for penetration to later than 300 ns (∼15 m; 100 MHz) and 800 ns (∼40 m; 25 MHz).

FIGURE 3

FIGURE 3. Processed GPR data from the Lateral Moraine Site. (A) Profile 1a: 100 MHz longitudinal profile at the crest of the moraine and (B) interpretation. (C) Profile 1b: 25 MHz longitudinal profile parallel to moraine crest and (D) interpretation. (E) Profile 2a: 100 MHz cross-sectional profile over the moraine and (F) interpretation. (G) Profile 2b: 25 MHz cross sectional profile over the moraine and (H) interpretation. (I) Satellite image with diagram of approximate GPR profile locations at the Lateral Moraine Site. N, S, P, and D indicate North, South, ice proximal, and ice distal, respectively. In B, D, F & H, solid lines indicate the snow-ground interface. Dashed lines mark continuous reflections, dotted lines delineate discontinuous reflections, and dot-dashed lines define changes in radargram character. G1 and G48 indicate the approximate locations of seismic geophones 1 and 48. X and yellow lines indicate the approximate crossover location of Profiles 1 and 2 on each of the radargrams.

The slope-intercept method applied to the seismic refraction data gives P-wave velocities of v₁ = 1,860 ± 60 m/s, v₂ = 2,420 ± 20 m/s, and v₃ = 4,460 ± 20 m/s, expressed with 95% confidence assuming first-break timing errors range from ±0.125 ms at near offsets to ±2 ms at far offsets, and geophone placement errors of ±0.2 m. Refractor depths are 1.5 ± 0.7 m for the interface between layers 1–2 (Z₁) and 8–13 m between layers 2–3 (Z₂). Z₁ is likely an average depth value for a refractor that shallows to the south, and the upper limit for Z₂ has been defined assuming bedrock below with a maximum velocity of 6,400 m/s.

Dispersion curve inversions from either end of the spread produced similar v_s profiles (Figures 4C–F). Both suggest a 4-layer structure with a velocity inversion at 2 ± 1 m. Shear velocity then increases at 4 ± 1 m depth, followed by a further increase from 600 to 1,400 m/s to 2,200–3,000 m/s at 13 ± 5 m depth (South; Figure 4C) and from 500 to 900 m/s to 1,000–2,200 m/s at 9 ± 3 m depth (North; Figure 4F).

FIGURE 4

FIGURE 4. MASW dispersion curves and inversion results from the Lateral Moraine Site: (A–C) use geophones G 2 to G 16 and the southern shotpoint, at −1.5 m, and (D–F) use geophones G 45 to G 25 and the northern shotpoint, at 72 m. Picked dispersion curves for each site are shown in A & D. All modelled dispersion curves tested in the inversion are shown in B & E, with corresponding misfits plotted as the colour scale, high misfit models progressively overlain by lower misfit models, and picked values displayed in black. C & F show all v_s profiles associated with the modelled dispersion curves and corresponding misfits, with minimum misfit v_s profile (black dashed line). The grey regions in C & F show approximate resolution limits defined by Z_min and Z_max (Table 1).

Combined, our data suggest three horizons at this site. The change in substrate firmness observed at ∼0.8 m in the test pits coincides approximately, within resolution limits, with the shallowest boundary detected with the two seismic methods. However, this is not observed in the GPR data due to direct wave arrivals obscuring any structure in the uppermost metre. GPR and MASW profiles suggest a boundary is present at 3–5 m, but this did not produce observable refracted wave arrivals. All three geophysical methods indicate that a horizon is present at ∼12 m depth and, although depth resolution differs across the datasets, seismic refraction places it between 8 and 13 m.

Terminal Moraine Site

Exposures E5–E7 (Terminal Moraine Site) reveal that the upper 40–50 cm of the subsurface consists of a clast-rich diamicton with a silty sand matrix. At a depth of 56 cm (E5), 42 cm (E6), and 38 cm (E7), the matrix becomes coarser, with increased concentrations of gravel. The very firm character of the lithofacies restricted further manual excavation. Sediments exposed in the stream cut (E4, Figure 1F) were similar to those identified in the test pits (E5–E7), but less firm. A gradual coarsening from pebbles to cobbles, and a decrease in firmness, was identified at ∼2.5 m below the moraine crest.

Based on migration velocity analysis, GPR velocities at the Terminal Moraine Site appeared to range from 0.07 m/ns to 0.17 m/ns. A velocity of 0.08 m/ns best-focused the majority of the image, hence it was used in depth conversions. Strong reflections are observed at ∼1–2 m depth on Profile 3a (Figures 5A,B) and at ∼10–14 m depth on Profile 3b (Figures 5C,D). A discontinuous reflector in the 100 MHz profile (3a), coincides (within resolution limits) with a distinct change in radargram character on the 25 MHz profile (3b) at about ∼3–5 m depth. Both frequencies have a low amplitude response directly below the moraine crest, ∼15 m wide and down to ∼6 m depth. Continuous, sub-horizontal reflectors are truncated on either side of this zone and the signal is attenuated below (Figures 5E–H). Penetration depths are smaller but more consistent laterally than at the Hummocky and Lateral Moraine sites. The 100 MHz profiles generally have a maximum penetration to 200 ns (∼8 m) and the 25 MHz profiles to 500 ns (∼20 m).

FIGURE 5

FIGURE 5. Processed GPR data from the Terminal Moraine Site. (A) Profile 3a: Eastern end (located on E) of the 100 MHz profile along the crest of the prominent terminal moraine and (B) interpretation. (C) Profile 3b: 25 MHz longitudinal profile, parallel to moraine crest, and (D) interpretation. (E) Profile 4a: 100 MHz transverse profile and (F) interpretation. (G) Profile 4b: 25 MHz transverse profile and (H) interpretation. (I) Satellite image with approximate GPR profile locations defined. W, E, P, and D indicate West, East, ice proximal, and ice distal, respectively. In B, D, F, H, solid lines indicate the snow-ground interface. Dashed lines mark continuous reflections, dotted lines delineate discontinuous reflections, and dot-dashed lines define changes in radargram character. G1 and G48 indicate the approximate locations of seismic geophones 1 and 48. X and yellow lines indicate the approximate crossover location of Profiles 3 and 4 on the radargrams.

A two-layer interpretation of the seismic refraction data from the Terminal Moraine Site (Figure 6) suggests an interface at between 0.8 and 4.0 m (Z₁), with a depth of ∼1.6 m most likely. In the upper layer, v_p is 2,460 ± 70 m/s and in the lower layer it increases to 2,920 ± 30 m/s. Velocities are stated with 95% confidence limits, with the same errors assumed as for the Lateral Moraine Site, excepting the poorly planted geophones G17−G22 being allocated ±0.5 m placement error. The large range of the possible depths for Z₁ is due to the assumptions of a planar, horizontal boundary for Eq. 1.

FIGURE 6

FIGURE 6. First 150 ms TWT of seismic gathers acquired at the Terminal Moraine Site, displayed in positive standard polarity (A) Forward shot from shot location S 1 and (B) reverse shot from shot location S 46. First-break picks marked by blue dashes and reverberant surface wave arrivals (potentially guided waves) highlighted in grey.

Dispersion curve analysis was undertaken separately for the shots at the ends of the spread to check for lateral heterogeneity (Figure 7). On the shot gathers (Figures 6A,B), high amplitude, monochromatic, reverberant wave trains are dominant beyond ∼8 m offsets (highlighted in grey). This character can be caused by waveguide effects, e.g., a subsurface layer with a large acoustic impedance contrast at its upper and lower interfaces. The continuous reflection of seismic energy within this layer gives rise to the observed resonant character of the waves, known as guided waves (Roth et al., 1998). An initial interpretation for the waveguide suggests a seismically slow layer (550–600 m/s) between ∼4 and 15 m depth, bounded by high impedance contrasts.

FIGURE 7

FIGURE 7. MASW dispersion curves and inversion results from the Terminal Moraine Site: (A–C) use geophones G 6 to G 48 and shotpoint S 1 (shot direction: west to east; W−E) and (D–F) use geophones G 41 to G 1 and shotpoint S 46 (shot direction: east to west; E−W). Picked dispersion curves for each site are shown in A & D. All modelled dispersion curves tested in the inversion are shown in B & E, with corresponding misfits plotted as the colour scale, high misfit models progressively overlain by lower misfit models, and picked values displayed in black. C & F show all v_s profiles associated with the modelled dispersion curves and corresponding misfits, with minimum misfit v_s profile (black dashed line). The grey regions in C & F show approximate resolution limits defined by Z_min and Z_max (Table 1).

The maximum depth of the pits (E5−E7) at the Terminal Moraine Site is below the minimum resolution of the geophysical data, so comparison in horizon depths is not possible between these datasets. Depths to horizons in the GPR and MASW data correspond well within resolution limits, with three horizons identified at about 1–2, 3–5, and 10–15 m depth. The depth to the refractor identified through seismic refraction analysis has high associated uncertainty but may align with the shallowest aforementioned horizon. The increase in grain size, identified at 2.5 m depth at section E4 does not coincide with any of these boundaries, however, E4 is not coincident with the seismic survey line or the displayed GPR profiles.

Hummocky Site

Exposure E1 revealed laterally alternating, well-sorted, sand-rich, and silt-rich lithofacies adjacent to a gravel unit at the ice-distal side of the hummock, with sharp, vertical contacts between the lithofacies units (Figure 1D). The internal structure of the hummock is complex with structures indicative of ductile and brittle deformation (e.g., folds and faults). Thus, these sediments are interpreted as glaciotectonites, i.e., sediments that have been deformed by subglacial shearing but retain structural characteristics of the parent material (Evans et al., 2006; Evans, 2018).

The radar velocities, identified through migration velocity analysis, range from 0.07 to 0.12 m/ns. For guideline depth estimates, a velocity of 0.09 m/ns was assumed for Profile 5 and 0.10 m/ns for Profile 6. Most of the 100 MHz profiles (e.g., Profiles 5a and 6) are characterised by discontinuous, in some cases chaotic, reflections (Figures 8A–E). Where Profile 6 (Figures 8E,F) crosses the logged moraine (E1), the 100 MHz data suggests a structurally complex subsurface. However, the resolution is insufficient to represent the numerous fine-scale deformation structures within the glaciotectonites observed in section E1. At the location of the seismic spread (Profiles 5a and b; Figures 8A−D), a dipping reflection in the 100 MHz data reaches the surface approximately 30 m south of G 1. Immediately south of this, there is a zone of low amplitude responses (identifiable in Profiles 5a and b), followed by a higher amplitude response and a break in the reflectors at the location of a meltwater channel, identified during summer fieldwork. The signal penetration varies across the area for both antenna frequencies. Maximum signal penetration changes abruptly on some profiles from 200 ns (∼9 m) to more than 500 ns (∼23 m) along 100 MHz profiles and 400 ns (∼18 m) to 800 ns (∼36 m) on 25 MHz profiles. It appears that attenuation is particularly high where the hummocky moraines are present, suggesting a high clay content, the presence of impure water or both.

FIGURE 8

FIGURE 8. Processed GPR data from the Hummocky Site. (A) Profile 5a: 100 MHz profile along the seismic survey line and (B) its interpretation. (C) Profile 5b: 25 MHz profile along the seismic survey line and (D) its interpretation. (E) Profile 6: 100 MHz profile over section E1 (Figure 1C). (F) Enlarged view of the moraine from Profile 6. (G) Satellite image with diagram of GPR profile orientations at the Hummocky Site. In B, D & F, solid lines show snow-ground interface, dashed and dotted lines delineate continuous and discontinuous reflections, respectively. Dot-dashed boxes outline columns of high amplitude reflections. G1 and G48 indicate the approximate locations of geophones 1 and 48. X and yellow lines indicate approximate crossover location between Profiles 5 and 6.

The refracted, reflected, and surface wave arrivals on the shot gathers from the Hummocky Site do not have simple arrival time and waveform characteristics (Figures 9A,B). Rayleigh wave propagation is disrupted at the green lines in Figures 9A,B, and first break travel-times do not increase with offset (purple lines on Figures 9A,B) and instead imply a refractor with significant topography under the middle of the spread (note, there is no surface topographic expression of this at the site). The seismic response and GPR Profiles 5a and b (Figures 8A–C) suggest that the structure is laterally heterogeneous and that refracting interfaces are non-planar. Since assumptions in the slope-intercept method are likely violated at this site, the interpretation is not performed and a more sophisticated algorithm is currently being explored. MASW (Figures 10A–D) was also impacted by the subsurface heterogeneities at the site, since it could only be performed for a restricted range of offsets. This limited the bandwidth of the picked dispersion curve to 20–60 Hz (Figures 10A–D), and hence restricted the vertical resolution limits of the v_s profiles to 2–10 m depth (Table 1). The best-fit calculated v_s models from both ends of the spread indicate v_s is around 400–500 m/s from 2 to 9 m depth (Figures 10C–F).

FIGURE 9

FIGURE 9. First 250 ms TWT of seismic gathers acquired at the Hummocky Site, displayed in positive standard polarity. Strong airwave components and potential back-scattered surface waves are highlighted and geophone offsets used in MASW analysis defined. Vertical dashed green lines mark the offsets at which a change in arrival characteristics occurs. (A) Forward shot taken at shotpoint S −1.5. (B) Reverse shot taken at shotpoint S 72.

FIGURE 10

FIGURE 10. MASW dispersion curves and inversion results from the Hummocky Site: (A–C) use geophones G 16 to G 23 and the ice proximal shotpoint, at −1.5 m and (D–F) use geophones G 47 to G 35 and the ice distal shotpoint, at 72 m. Picked dispersion curves for each site are shown in A & D. All modelled dispersion curves tested in the inversion are shown in B & E, with corresponding misfits plotted as the colour scale, high misfit models progressively overlain by lower misfit models, and picked values displayed in black. C & F show all v_s profiles associated with the modelled dispersion curves and corresponding misfits, with minimum misfit v_s profile (black dashed line). The grey regions in C & F show approximate resolution limits defined by Z_min and Z_max (Table 1).

Unlike at the Lateral and Terminal Moraine sites, the complex subsurface architecture at the Hummocky Site limited combined analysis of horizon depths. The GPR profiles and seismic sections indicate subsurface heterogeneity and structural complexity on a larger spatial scale (metre to tens of metres) than can be observed at section E1 (centimetre to tens of centimetres). Under the seismic survey line, probable refractor topography is corroborated by concave upwards GPR reflectors around the meltwater channel, 40–60 m from the southern end of Profile 5a (Figure 8B).

Discussion

Synthesis of Acquired Geophysical Data

Our results demonstrate that geophysical techniques can deliver spatial information on the structure and depth of substrates in proglacial environments which traditional methods alone cannot access. However, some glacial sediment-landform associations are more amenable to geophysical survey than others: survey designs must be adapted for particular site conditions and target structures. Here, we compare the results from the Hummocky and Lateral Moraine sites to show the strong influence of structural complexity on the suitability of the applied techniques, and to explore the effects of surface sediment type we use results from the Lateral and Terminal Moraine sites. We then offer recommendations for future geophysical survey campaigns in proglacial, and analogous (e.g., periglacial) settings.

Influence of Subsurface Structural Complexity

Sediment excavation at the Hummocky Site showed that the hummocks contain glaciotectonised sediments, which are structurally complex due to high stress and strain histories (Evans et al., 2006; Evans, 2018). GPR profiles across the area (Figures 8A–F) support this complexity and show that subsurface conditions are also highly variable on a larger scale (Figure 11I). The GPR surveys over the Lateral Moraine Site revealed a simpler subsurface structure with continuous, subparallel reflectors (Figures 3A−D) detailing a stratified substrate with moderately undulating boundaries (Figure 11A). The contrasting subsurface conditions at these two sites led to considerably different seismic responses (Figures 2A,B, 8A,B). The v_p values and, to a lesser extent, unit thicknesses at the Lateral Moraine Site could be approximated using the slope-intercept interpretation method owing to its simple structure (Hoffmann and Schrott, 2003). By contrast, the heterogeneity at the Hummocky Site violates the plane-layer assumption of the slope-intercept method and requires more advanced interpretation techniques, such as Wavepath Eikonal Traveltime (WET) (Woodward, 1992), the General Reciprocal Method (GRM) (Palmer, 1981), or Full Waveform Inversion (FWI) (Tarantola, 1986), to approximate the substrate v_p values. It must be considered that associated interpretations would still be ambiguous (Whiteley and Eccleston, 2006) and improved precision in results should not necessarily be regarded as increased certainty (Van Overmeeren, 2001).

FIGURE 11

FIGURE 11. Summary diagrams of the subsurface conditions beneath the seismic spread locations at the Lateral Moraine Site (A–D), the Terminal Moraine Site (E–H), and the Hummocky Site (I–K). (A), (E) & (I) Key features identified from GPR radargrams. Dotted lines depict dominant reflections from 100 MHz profiles (Figures 4B, 6B, 9B), dashed lines from 25 MHz profiles (Figures 4D, 6D, 9D), dot-dashed lines indicate a strong change in radargram character, and the solid lines delineate the surface topography. Grey boxes show approximate areas over which MASW was performed with grey line segments corresponding to shading in respective MASW plots (C, D; G, H & J, K). (B) & (F) Seismic refraction results at the Lateral Moraine Site and Terminal Moraine Site, respectively, positioned at centre of seismic line location. (C), (D), (G), (H), (J) & (K) The v_s models with the lowest misfit (solid line) and error bounds displayed as solutions with a 10% fit to the data (dashed lines) from MASW analysis over the areas defined in A, E, and I. Grey shading shows identified layers; darker represents higher velocity. Small triangles show ends of seismic spread locations.

For MASW, problems of heterogeneity were mitigated by undertaking separate v_s inversions at each end of the seismic line (Figures 10C,D, 11J,K). This approach worked well at the Lateral Moraine Site where well-defined fundamental modes could be picked up to 80 and 100 Hz (Figures 4A–D). However, at the Hummocky Site, the restricted geophone offsets and extensive spatial variability of the subsurface conditions led to reduced precision and restricted bandwidths (to 60 Hz) of the dispersion curve picks (Park et al., 2001; Figures 10A–D). MASW inversions can also benefit from constrained layer boundaries, to reduce ambiguity in inversions. The simplicity of structure at the Lateral Moraine Site meant that such constraints were applicable from both the seismic refraction and GPR datasets, resulting in well-constrained v_s models. However, constraints were unavailable for the more complex Hummocky Site, hence the corresponding MASW inversions are less precise. Nevertheless, the non-uniqueness of MASW results, compared to GPR and other seismic methods, must always be considered to avoid over-interpretation of the resultant v_s profiles.

While the seismic methods applied here were poorly suited to the complex glaciotectonites found at the Hummocky Site, both frequencies in the GPR surveys provided useful images of subsurface structures (e.g., Figures 8A,C,E). A joint interpretation of these GPR data allowed shallow (<10 m depth) structures to be defined at sub-metre scale spatial resolution (e.g., Figure 8F), and broader-scale structure to be imaged to ∼35 m depth. This combined interpretation highlights the value of multi-frequency GPR surveying (e.g., Booth et al., 2009).

The terrain at the Hummocky Site does not provide a clear continuous ridge to guide profile orientations, as is the case at the Lateral Moraine Site. Instead, the adoption of a survey grid (Figure 8G) aided spatial positioning of the profiles during processing and meant that the complex structures could be viewed from orthogonal orientations, reducing the risk of misinterpretation of out-of-plane reflections (Neal, 2004). The interpretation of complex 3-D structures could be aided by undertaking 3-D GPR surveys, but this requires a dense grid spacing with spatial sampling at a maximum of a quarter wavelength (e.g., 0.25 m trace spacing for 100 MHz and 1 m for 25 MHz surveys), corresponding positional accuracy, and more computationally demanding 3-D migration (Grasmueck et al., 2005). Therefore, 3-D surveying is more labour and cost intensive than traditional 2-D surveying (Lehmann and Green, 1999). However, in structurally complex situations, 3-D GPR surveys may prove to be a more valuable use of field time than attempting advanced seismic investigations or potentially misinterpreting sparse 2-D GPR profiles.

Influence of Surface Composition

Structurally, the Lateral Moraine and Terminal Moraine sites provide comparable conditions for geophysical techniques, but their contrasting surface materials present different challenges to subsurface interpretation. The open blockwork of coarse gravel, cobbles, and occasional boulders at the surface of the Terminal Moraine Site led to signal scattering and abundant diffraction hyperbolae in the unmigrated GPR data (Supplementary Figures S1A−E). The simplified velocity analysis used here cannot represent the complex velocity structure of the moraine, so the scattered signal is poorly migrated in the uppermost 25 ns (100 MHz; e.g., Figures 5A–F) and 100 ns (25 MHz; e.g., Figures 5C–H) of the radargrams. This signal scattering may be the cause of the restricted GPR investigation depth at this site, reaching a maximum of ∼25 m, compared to the ∼40 m reached at the Lateral Moraine Site. At the Lateral Moraine Site, the changes in signal penetration (Figure 3C) are more likely due to spatial variations in liquid water content, with particularly wet conditions observed in the field at the northern end of Profiles 1a and b (Figures 3A–D).

Visual inspection of the seismic data at the Lateral Moraine Site (Figures 2A,B) and the Terminal Moraine Site (Figures 6A,B) reveals differences in the main arrival types. A large portion of the seismic energy at the Terminal Moraine Site appears as monochromatic, reverberant wave trains (potentially guided waves; Figures 6A,B), which exhibit little change across the gather. In contrast, at the Lateral Moraine Site, the surface wave component arrives as Rayleigh waves which change character (i.e., amplitude and frequency) from trace to trace (Figure 2B). The MASW process has been developed to utilise Rayleigh wave energy and it is still not clear how well it performs with guided waves, so the results from the Terminal Moraine Site should be treated with caution. Despite the differences in arrival types, the inversion processes define similar subsurface structures for forward and reverse shots at both sites (Figures 11C,D,G,H), with depths to unit boundaries corresponding, within resolution limits, to prominent reflectors identified on coincident GPR profiles (Figures 11A–E, respectively). The first arrivals (i.e., direct and critically refracted waves) at the Lateral Moraine Site are identifiable across the full seismic records, out to 72 m offset (Figures 2A,B). At the Terminal Moraine Site, first arrivals are poorly defined beyond ∼25 m offset (Figures 6A,B). This is likely to have been influenced by inadequately planted geophones in the frozen, blocky, exposed surface and hence, a high level of noise as the geophones had minimal protection from the wind. Therefore, refraction results have higher associated uncertainty at the Terminal Moraine Site than the Lateral Moraine Site.Where co-located measurements of v_p and v_s are obtained from seismic refraction analysis and MASW, the subsurface can be interpreted with reduced ambiguity. Due to the different influences on P− and S wave propagation, v_s is less sensitive to fluid fill than v_p. Therefore, co-located v_s and v_p profiles can reveal subsurface interfaces that would be missed using a standalone seismic interpretation, for example, the water table will not be evident on a v_s profile but would exhibit as a contrast in v_p with depth. At the Lateral Moraine Site the v_p/v_s ratio increases from 2 to 6 at ∼2 m depth, which is typical for unconsolidated saturated sediments overlying frozen or unsaturated sediments (e.g., Zimmerman and King 1986; Schön 2015). However, due to overlap in the seismic velocities of bedrock and frozen substrate, additional data are required to make a definitive interpretation of the half spaces identified at the Lateral and Terminal Moraine sites (Killingbeck et al., 2020).

The identification of frozen ground has frequently been achieved using techniques that investigate changes in substrate resistivity (Hauck and Kneisel, 2008). This includes electrical techniques, for example, direct current (DC) resistivity soundings (e.g., Yoshikawa et al., 2006) and electrical resistivity tomography (ERT), and electromagnetic approaches, such as time-domain (TDEM) and frequency-domain (FDEM) electromagnetic soundings (e.g., Maurer and Hauck, 2007; Killingbeck et al., 2020). The contribution of such techniques to foreland investigations is beyond the scope of this study, but it is clear that they can be a valuable addition to geophysical campaigns in proglacial and analogous environments (e.g., Hauck et al., 2011).

Although there are limitations on the spatial resolution and applicability of the geophysical techniques presented here, particularly where complex sedimentary sequences are present, the results from the Lateral and Terminal Moraine sites show how combined GPR and seismic analysis can characterise the internal properties of proglacial landforms when surveys are appropriately designed. Results from the Hummocky Site highlight the challenges presented by the proglacial environment. The complexity of the site and centimetre scale structures observed at section E1 shows that detailed logging of sedimentary sections remains essential for genetic interpretations, and as such we do not provide interpretations for our three sites here. By contrast, at the Lateral and Terminal Moraine sites, the ground conditions hampered excavation efforts hence the additional insight from the geophysical surveys was particularly valuable. As a rule, however, we would always recommend interpreting geophysical data alongside sedimentological observations (where possible) to reduce ambiguity and provide detailed information on the subsurface.

Recommendations for Geophysical Investigations in Analogous Settings

Although sedimentological studies can identify individual lithofacies and sedimentary structures, the information they provide is spatially limited. Sediment sections are valuable at the fine scale but are highly localised, often requiring extrapolation from just one or two sections, whereas geophysical methods are useful at a broader scale and to greater depths (i.e., entire landform-suites). Through comparison of the outcomes at the three proglacial sites studied here, and with guidance from published results from other sediment studies, we provide a framework for future geophysical investigations in analogous settings (Figure 12) and outline the applicability of each technique to common targets (Table 2). We recommend a survey approach that:

- visits the survey sites in snow-free conditions to identify profiles with minimal topography (e.g., along moraine ridges) and enough space for a seismic spread length ∼5 times the expected target depth (determined, where possible, from published literature),

- uses 2-D GPR results to guide the seismic survey design (i.e., survey location, spread length, geophone spacing, shotpoint positions),

- performs refraction analysis and MASW on seismic data (where appropriate) to provide complimentary v_p and v_s profiles from single seismic surveys,

- co-locates GPR and seismic surveys to constrain seismic inversions and enable combined interpretation,

- acquires GPR data with multiple antenna frequencies, including one ≤50 MHz system to counter high signal attenuation in glacial deposits,

- acquires GPR surveys along orthogonal profiles, where possible guided by landform orientation; if highly complex architecture is present, focus field effort on a 3-D GPR survey,

- undertakes geophysical surveys when snow cover is present (if possible), for easier data acquisition,

- includes direct observations (e.g., sediment logging) to enable past process interpretations (where possible).

FIGURE 12

FIGURE 12. Suggested workflow for undertaking combined GPR, seismic refraction, and MASW surveys in a proglacial environment.

TABLE 2

TABLE 2. A summary of which techniques are suggested for providing information on typical proglacial investigation targets.

The following section explains the motivation for these suggestions.

The contrasting data quality acquired from the three sites at Midtdalsbreen highlights the importance of a comprehensive desk study prior to any near-surface geophysical campaign in proglacial environments to aid survey planning. In particular, a key consideration is the viability and applicability of certain geophysical survey methods, depending on surface morphology, surface material, and predicted subsurface composition at the target sites. To this end, remotely sensed data (e.g., aerial photographs, satellite imagery) should be used to guide initial target selection and to make geomorphological interpretations of these often-remote sites. . Where possible, in situ reconnaissance should be undertaken in snow-free conditions prior to geophysical surveying in order to (a) assess site conditions (e.g., surface sediment types), (b) identify obstacles to geophysical techniques (e.g., large boulders or meltwater streams), and (c) guide survey designs as well as equipment and parameter selection (e.g., GPR antenna frequencies, GPR sampling intervals, seismic array lengths, geophone spacing). Acquiring, processing, and interpreting GPR data prior to seismic surveying can help target seismic acquisitions and lead to more successful investigations. Radargrams can be used to locate a laterally homogeneous and representative subsurface and inform the seismic survey design with respect to expected target depths and subsurface architecture (e.g., with resolution and offset considerations in mind). For example, when highly heterogeneous substrates are present, such as at the Hummocky Site, simple slope-intercept seismic refraction interpretation is not appropriate; more sophisticated interpretation techniques are required and surveys must be planned accordingly. For example, WET demands a dense shot spacing to ensure adequate ray-path coverage (Whiteley and Eccleston, 2006). It should be considered that some sites are inherently unsuitable for the seismic refraction method (e.g., where no substantial velocity contrasts are present or v_p decreases with depth) and applying more sophisticated interpretation tools in these situations will not yield more reliable results (Whiteley and Eccleston, 2006). Where a complex subsurface (e.g., highly deformed glacial sedimentary sequences) is identified through initial sedimentological observations and GPR surveys, a focused 3-D GPR survey may be more appropriate.

The joint application of seismic refraction and MASW can greatly aid subsurface interpretation when initial GPR surveys reveal a layered subsurface with continuous reflectors and only minor lateral heterogeneities, such as along the lateral moraine (Figures 3A−D). Field efforts can be maximised if the seismic survey design is compatible with both refraction and MASW analysis (e.g., by deploying geophones with a low natural frequency). Performing refraction and MASW analysis on the data can deliver depths to horizons and unit elastic properties (v_p, v_s, and their derivatives, e.g., v_p/v_s ratio) that compliment structural and electrical information delivered by coincident GPR profiles. If field sampling for density measurements can be undertaken, further physical properties such as shear modulus (rigidity) can be defined (Pegah and Liu, 2016; Clarke, 2018). Knowledge of such properties can be useful for engineering purposes or, for example, to aid our understanding of landform stress histories or preservation potential (Clarke, 2017). If the ground surface is rocky or frozen, performing geophysical surveys when the ground is snow covered can improve geophone coupling for seismic surveys and provide a smooth surface for GPR surveying. For seismic surveys, wind noise can be reduced by planting geophones in snow pits, where snow depth is sufficient (e.g., Rossi et al., 2018).

In this study, we have shown the benefits of incorporating multiple antenna frequencies in geophysical campaigns. The GPR results from all three sites demonstrate the value of low frequency (≤50 MHz) GPR surveys in proglacial environments to provide adequate depth return where liquid water, abundant scatterers (e.g., large cobbles), and fine sediments are commonly present. However, the low frequency systems lack the fine scale resolution achieved with higher frequency antennas: the 100 MHz data revealed decimetre scale lithofacies that could not be resolved in the 25 MHz datasets. Thus, combining low frequency (≤50 MHz) and higher frequency (≥100 MHz) is beneficial in proglacial environments. Nonetheless, the combination of GPR frequencies used in this study was not able to image the level of detail observed in the sediment exposures. Greater resolution could be achieved with higher frequency (e.g., 200 MHz) GPR systems with the caveats of poor penetration depth and abundant signal scattering (e.g., Lukas and Sass, 2011; Pellicer and Gibson, 2011). Such higher frequency GPR antennas are thus best suited to sites where sand and gravel are dominant. Regardless of antenna frequency, subsurface interpretation also benefits from accurate migration and, by extension, accurate positional and velocity information to provide more representative structures and better vertical accuracy (Blindow et al., 2007; Booth et al., 2010). To enable process-form reconstructions in glacial geological studies, accurate migration is important, as the internal structural architecture of landforms (and thus reflector geometry) provides critical information on depositional and deformation histories of landforms (Evans and Benn, 2021). However, migration is not routinely used in glacial geological studies (e.g. Midgley et al., 2013). Accurately representing the internal structural architecture is particularly significant in the interpretation of moraines (Benn and Lukas, 2006; Benediktsson et al., 2010; Benediktsson et al., 2015), as these landforms (a) may contain sedimentary structures that do not conform to the surface slopes (e.g., overturned folds, thrusts) and (b) have morphologies that are the most challenging in migration (c.f. Allroggen et al., 2014). Thus, topographic migration with representative signal velocities and accurate elevation data is necessary for these landforms.

In this study, velocity estimates are determined using iterative migration. Where the subsurface is approximately laterally homogeneous (e.g., along the crest of the Terminal Moraine Site), these estimates can be improved through the implementation of separable antennas for common midpoint surveys (where possible). Where the subsurface is heterogeneous (e.g., at the Hummocky Site), multiple common midpoint surveys could be applied to determine typical velocities for the area. However, these will not fully represent the complex velocity structure, so the uncertainty in depth estimates will remain high (Carrick Utsi, 2017). Even with improved velocity accuracy and a higher frequency GPR system, the level of detail and unambiguous material properties provided by sedimentological observations cannot be achieved. Therefore, glacial sedimentological investigations are required, both for ground truthing geophysical data and for interpreting depositional processes at the individual landform scale.

Conclusion

A combined geophysical approach to the investigation of glacial sediments can provide valuable information on subsurface material properties and architecture if appropriate survey practices are applied. Sediment-landform associations in proglacial environments have traditionally been studied using standard glacial geological methods, but a combination of problems (e.g., limited moraine exposures) often precludes extensive subsurface investigations in proglacial settings. Through applying GPR, seismic refraction, and MASW across the foreland of Midtdalsbreen (Norway) we have highlighted some of the strengths and limitations of these techniques in investigations of glacial landforms and outlined a suggested framework for their application in foreland characterisation, and analogous studies. The success of the individual methods applied here varies depending on subsurface composition and structural complexity. Combined analysis of the datasets is required to overcome the inherent non-uniqueness in the results and allow sedimentological interpretations to be drawn.

GPR provides laterally extensive information on subsurface architecture in the investigated proglacial settings. Combining low and mid-range antenna frequencies overcame the resolution-depth trade off, enabling delineation of near surface unit boundaries and identification of interfaces at up to 40 m depth. Seismic surveys can be challenging in the proglacial setting, and the geological complexity of targets can undermine assumptions for simple seismic interpretations. We suggest undertaking GPR surveys and sediment studies prior to seismic investigations, so that survey practices and interpretation methods can be optimised for the specific site conditions. Where seismic methods are viable, applying both refraction and surface wave analysis to the seismic data optimises field efforts and can produce v_p and v_s values for the substrates, to further reduce the ambiguity of the possible interpretations. The recommendations outlined here can help future analogous studies to optimise geophysical investigations, leading to greater insights into past glacier behaviour and ongoing changes in the surrounding landscape.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: Figshare: https://doi.org/10.6084/m9.figshare.c.5828600.v1.

Author Contributions

Fieldwork preparation, data collection, and analysis were performed by HW, BR, AB, and SK with additional data analysis performed by RC. The first draft of the manuscript was written by HW and all authors commented on subsequent versions of the manuscript. All authors read and approved the final manuscript.

Funding

The research leading to these results has received funding from the European Union’s Horizon 2020 project INTERACT, under grant agreement No 730938; The Bolin Centre for Climate Research; Albert and Maria Bergström's Foundation; and the Gerard de Geer Foundation. Funding from the aforementioned bodies covered the costs of four field campaigns across the years 2018 and 2019, including lodging, transport, hiring of personnel and fuel costs. Benjamin Chandler was funded by a grant from The Leverhulme Trust (MULTIPLEX; SAS-2019-075). Open access publication fees are covered by the Stockholm University Library.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

Niklas Allroggen provided updates to a topographic migration algorithm for our dataset; Richard Gyllencreutz supplied the GPS equipment; Erika Leslie and marit Tangen gave invaluable logistical help at Finse Alpine Research Station; Kjell Magne Tangen provided support and local knowledge during field campaigns; Sven Lukas gave sedimentological input in the field; and Emma Pearce helped with seismic data collection. The manuscript benefited from input from Dominik Gräff and, Lukas Arenson and the editor Alun Hubbard.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2022.734682/full#supplementary-material

References

Åkesson, H., Nisancioglu, K. H., Giesen, R. H., and Morlighem, M. (2017). Simulating the Evolution of Hardangerjøkulen Ice Cap in Southern Norway since the Mid-holocene and its Sensitivity to Climate Change. The Cryosphere. 11, 281–302. doi:10.5194/tc-11-281-2017

CrossRef Full Text | Google Scholar

Aki, K. (1957). Space and Time Spectra of Stationary Stochastic Waves, with Special Reference to Microtremors. Bull. Earthq. Res. Inst. 35, 415–456.

Google Scholar

Allroggen, N., Tronicke, J., Delock, M., and Böniger, U. (2014). Topographic Migration of 2D and 3D Ground‐penetrating Radar Data Considering Variable Velocities. Near Surf. Geophys. 13, 253–259. doi:10.3997/1873-0604.2014037

CrossRef Full Text | Google Scholar

Andersen, J. L., and Sollid, J. L. (1971). Glacial Chronology and Glacial Geomorphology in the Marginal Zones of the Glaciers, Midtdalsbreen and Nigardsbreen, South Norway. Norsk Geografisk Tidsskrift - Norwegian J. Geogr. 25, 1–38. doi:10.1080/00291957108551908

CrossRef Full Text | Google Scholar

Baker, P. L. (1991). Response of Ground-Penetrating Radar to Bounding Surfaces and Lithofacies Variations in Sand Barrier Sequences. Exploration Geophys. 22, 19–22. doi:10.1071/EG991019

CrossRef Full Text | Google Scholar

Beedle, M. J., Menounos, B., Luckman, B. H., and Wheate, R. (2009). Annual Push Moraines as Climate Proxy. Geophys. Res. Lett. 36, L20501. doi:10.1029/2009GL039533

CrossRef Full Text | Google Scholar

Benediktsson, Í. Ö., Schomacker, M. D., Geiger, A. J., Ingólfsson, Ó., and Guðmundsdóttir, E. R. (2015). Architecture and Structural Evolution of an Early Little Ice Age Terminal Moraine at the Surge-type Glacier Múlajökull, Iceland. J. Geophys. Res. Earth Surf. 120, 1895–1910. doi:10.1002/2015JF003514

CrossRef Full Text | Google Scholar

Benediktsson, Í. Ö., Schomacker, H., and Ingólfsson, Ó. (2010). The 1890 Surge End Moraine at Eyjabakkajökull, Iceland: a Re-assessment of a Classic Glaciotectonic Locality. Quat. Sci. Rev. 29, 484–506. doi:10.1016/j.quascirev.2009.10.004

CrossRef Full Text | Google Scholar

Benn, D., and Evans, D. (2010). Glaciers and Glaciation. 2nd ed. London: Routledge.

Google Scholar

Benn, D. I., and Lukas, S. (2006). Younger Dryas Glacial Landsystems in North West Scotland: an Assessment of Modern Analogues and Palaeoclimatic Implications. Quat. Sci. Rev. 25, 2390–2408. doi:10.1016/j.quascirev.2006.02.015

CrossRef Full Text | Google Scholar

Bennett, M. R., Huddart, D., Waller, R. I., Cassidy, N., Tomio, A., Zukowskyj, P., et al. (2004). Sedimentary and Tectonic Architecture of a Large Push Moraine: a Case Study from Hagafellsjökull-Eystri, Iceland. Sediment. Geology. 172, 269–292. doi:10.1016/j.sedgeo.2004.10.002

CrossRef Full Text | Google Scholar

Biot, M. A. (1962). Mechanics of Deformation and Acoustic Propagation in Porous Media. J. Appl. Phys. 33, 1482–1498. doi:10.1063/1.1728759

CrossRef Full Text | Google Scholar

Blindow, N., Eisenburger, D., Illich, B., Petzold, H., and Richter, T. (2007). “Ground Penetrating Radar,” in “Ground Penetrating Radar,” in Environmental Geology: Handbook of Field Methods and Case Studies. Editors K. Knödel, G. Lange, and H.-J. Voigt (Berlin, Heidelberg: Springer Berlin Heidelberg), 283–335. doi:10.1007/978-3-540-74671-3_10:

CrossRef Full Text | Google Scholar

Booth, A. D., Clark, R. A., Hamilton, K., and Murray, T. (2010). Multi-offset Ground Penetrating Radar Methods to Image Buried Foundations of a Medieval Town wall, Great Yarmouth, UK. Archaeol. Prospect. 17, a–n. doi:10.1002/arp.377

CrossRef Full Text | Google Scholar

Booth, A. D., Endres, A. L., and Murray, T. (2009). Spectral Bandwidth Enhancement of GPR Profiling Data Using Multiple-Frequency Compositing. J. Appl. Geophys. 67, 88–97. doi:10.1016/j.jappgeo.2008.09.015

CrossRef Full Text | Google Scholar

Busby, J. P., and Merritt, J. W. (1999). Quaternary Deformation Mapping with Ground Penetrating Radar. J. Appl. Geophys. 41, 75–91. doi:10.1016/S0926-9851(98)00050-0

CrossRef Full Text | Google Scholar

Carrick Utsi, E. (2017). Survey Strategies, In. Ground Penetrating Radar, 73–82. doi:10.1016/b978-0-08-102216-0.00007-2

CrossRef Full Text | Google Scholar

Chandler, B. M. P., Evans, D. J. A., and Roberts, D. H. (2016a). Characteristics of Recessional Moraines at a Temperate Glacier in SE Iceland: Insights into Patterns, Rates and Drivers of Glacier Retreat. Quat. Sci. Rev. 135, 171–205. doi:10.1016/j.quascirev.2016.01.025

CrossRef Full Text | Google Scholar

Chandler, B. M. P., Evans, D. J. A., and Roberts, D. H. (2016b). Recent Retreat at a Temperate Icelandic Glacier in the Context of the Last ∼80 Years of Climate Change in the North Atlantic Region. arktos. 2, 1–13. doi:10.1007/s41063-016-0024-1

CrossRef Full Text | Google Scholar

Clarke, B. G. (2017). Engineering of Glacial Deposits. London, United Kingdom: CRC Press 10.1201/9781315149356.

Google Scholar

Clarke, B. G. (2018). The Engineering Properties of Glacial Tills. Geotechnical Res. 5, 262–277. doi:10.1680/jgere.18.00020

CrossRef Full Text | Google Scholar

Davis, J. L., and Annan, A. P. (1989). Ground-penetrating Radar for High-Resolution Mapping of Soil and Rock Stratigraphy1. Geophys. Prospect. 37, 531–551. doi:10.1111/j.1365-2478.1989.tb02221.x

CrossRef Full Text | Google Scholar

Dobiński, W., Grabiec, M., and Glazer, M. (2017). Cold-temperate Transition Surface and Permafrost Base (CTS-PB) as an Environmental axis in Glacier-Permafrost Relationship, Based on Research Carried Out on the Storglaciären and its Forefield, Northern Sweden. Quat. Res. 88, 551–569. doi:10.1017/qua.2017.65

CrossRef Full Text | Google Scholar

Evans, D. J. A. (2003). Glacial Landsystems. London: Arnold.

Google Scholar

Evans, D. J. A. (2018). Till: A Glacial Process Sedimentology. The Cryosphere Scenice Series. Chichester: Wiley.

Google Scholar

Evans, D. J. A., and Benn, D. I. (2021). “A Practical Guide to the Study of Glacial Sediments,”. Editors D. J. A. Evans, and D. I. Benn. 2nd ed (London: Quaternary Research Association).

Google Scholar

Evans, D. J. A., Phillips, E. R., Hiemstra, J. F., and Auton, C. A. (2006). Subglacial till: Formation, Sedimentary Characteristics and Classification. Earth-Science Rev. 78, 115–176. doi:10.1016/j.earscirev.2006.04.001

CrossRef Full Text | Google Scholar

Fitzsimons, S., and Howarth, J. (2020). Development of Push Moraines in Deeply Frozen Sediment Adjacent to a Cold‐based Glacier in the McMurdo Dry Valleys, Antarctica. Earth Surf. Process. Landforms. 45, 622–637. doi:10.1002/esp.4759

CrossRef Full Text | Google Scholar

Gazetas, G. (1982). Vibrational Characteristics of Soil Deposits with Variable Wave Velocity. Int. J. Numer. Anal. Methods Geomech. 6, 1–20. doi:10.1002/nag.1610060103

CrossRef Full Text | Google Scholar

Grasmueck, M., Weger, R., and Horstmeyer, H. (2005). Full-resolution 3D GPR Imaging. Geophysics. 70, K12–K19. doi:10.1190/1.1852780

CrossRef Full Text | Google Scholar

Hagedoorn, J. G. (1959). The Plus-Minus Method of Interpreting Seismic Refraction Sections*. Geophys. Prospect. 7, 158–182. doi:10.1111/j.1365-2478.1959.tb01460.x

CrossRef Full Text | Google Scholar

Hambrey, M. J., and Glasser, N. F. (2012). Discriminating Glacier thermal and Dynamic Regimes in the Sedimentary Record. Sediment. Geology. 251-252, 1–33. doi:10.1016/j.sedgeo.2012.01.008

CrossRef Full Text | Google Scholar

Hannesdóttir, H., Ađalgeirsdóttir, G., Jóhannesson, T., Guđmundsson, S., Crochet, P., Ágústsson, H., et al. (2015). Downscaled Precipitation Applied in Modelling of Mass Balance and the Evolution of Southeast Vatnajökull, Iceland. J. Glaciol. 61, 799–813. doi:10.3189/2015JoG15J024

CrossRef Full Text | Google Scholar

Hauck, C., Böttcher, M., and Maurer, H. (2011). A New Model for Estimating Subsurface Ice Content Based on Combined Electrical and Seismic Data Sets. The Cryosphere. 5, 453–468. doi:10.5194/TC-5-453-2011

CrossRef Full Text | Google Scholar

Hauck, C., Isaksen, K., Vonder Mühll, D., and Sollid, J. L. (2004). Geophysical Surveys Designed to Delineate the Altitudinal Limit of Mountain Permafrost: An Example from Jotunheimen, Norway. Permafrost Periglac. Process. 15, 191–205. doi:10.1002/ppp.493

CrossRef Full Text | Google Scholar

C. Hauck, and C. Kneisel. (Editors) (2011). A New Model for Estimating Subsurface Ice Content Based on Combined Electrical and Seismic Data Sets. Cambridge: Cambridge University Press.

Google Scholar

Hausmann, H., Krainer, K., Brückl, E., and Mostler, W. (2007). Internal Structure and Ice Content of Reichenkar Rock Glacier (Stubai Alps, Austria) Assessed by Geophysical Investigations. Permafrost Periglac. Process. 18, 351–367. doi:10.1002/ppp.601

CrossRef Full Text | Google Scholar

Hawkins, L. V. (1961). The Reciprocal Method of Routine Shallow Seismic Refraction Investigations. Geophysics. 26, 806–819. doi:10.1190/1.1438961

CrossRef Full Text | Google Scholar

Hecht, S. (2003). Differentiation of Loose Sediments with Seismic Refraction Methods Potentials and Limitations Derived from Case Studies. Geophysical Applications in Geomorphology. Zeitschift fur Geomorphologie. Editors L. Schrott, A. Hördt, and R. Dikau 132, 89–102.

Google Scholar

Hoffmann, T., and Schrott, L. (2003). "Determing Sediment Thickness of Talus Slopes and valley Fill Deposits Using Seismic Refraction - A Comparision of 2D Interpretation Tools. "Geophysical Applications in Geomorphology. Zeitschift fur Geomorphologie. Editors L. Schrott, A. Hördt, and R. Dikau 132, 71–87.

Google Scholar

Ivanov, J., Miller, R. D., and Tsoflias, G. (2008). Some Practical Aspects of MASW Analysis and Processing. Symp. Appl. Geophys. Eng. Environ. Probl. Proc., 1186–1198. doi:10.4133/1.2963228

CrossRef Full Text | Google Scholar

Jakobsen, P. R., and Overgaard, T. (2002). Georadar Facies and Glaciotectonic Structures in Ice Marginal Deposits, Northwest Zealand, Denmark. Quat. Sci. Rev. 21, 917–927. doi:10.1016/S0277-3791(01)00045-2

CrossRef Full Text | Google Scholar

Jol, H. M., and Bristow, C. S. (2003). GPR in Sediments: Advice on Data Collection, Basic Processing and Interpretation, a Good Practice Guide. Geol. Soc. Lond. Spec. Publications. 211, 9. doi:10.1144/gsl.sp.2001.211.01.02

CrossRef Full Text | Google Scholar

Killingbeck, S. F., Booth, A. D., Livermore, P. W., Bates, C. R., and West, L. J. (2020). Characterisation of Subglacial Water Using a Constrained Transdimensional Bayesian Transient Electromagnetic Inversion. Solid Earth. 11, 75–94. doi:10.5194/se-11-75-2020

CrossRef Full Text | Google Scholar

Killingbeck, S. F., Booth, A. D., Livermore, P. W., West, L. J., Reinardy, B. T. I., and Nesje, A. (2019). Subglacial Sediment Distribution from Constrained Seismic Inversion, Using MuLTI Software: Examples from Midtdalsbreen, Norway. Ann. Glaciol. 60, 206–219. doi:10.1017/aog.2019.13

CrossRef Full Text | Google Scholar

Killingbeck, S. F., Livermore, P. W., Booth, A. D., and West, L. J. (2018). Multimodal Layered Transdimensional Inversion of Seismic Dispersion Curves With Depth Constraints. Geochem. Geophys. Geosyst. 19, 4957–4971. doi:10.1029/2018GC008000

CrossRef Full Text | Google Scholar

Krüger, J. (1995). Origin, Chronology and Climatological Significance of Annual-Moraine Ridges at Myrdalsjökull, Iceland. The Holocene. 5, 420–427. doi:10.1177/095968369500500404

CrossRef Full Text | Google Scholar

Kunz, J., and Kneisel, C. (2020). Glacier-Permafrost Interaction at a Thrust Moraine Complex in the Glacier Forefield Muragl, Swiss Alps. Geosciences. 10, 205. doi:10.3390/geosciences10060205

CrossRef Full Text | Google Scholar

Lankston, R. W. (1990). 3. High-Resolution Refraction Seismic Data Acquisition and Interpretation.Geotechnical and Environmental Geophysics. 1, 45, 74. doi:10.1190/1.9781560802785.ch3

CrossRef Full Text | Google Scholar

Lecomte, I., Thollet, I., Juliussen, H., and Hamran, S.-E. (2008). Using Geophysics on a Terminal Moraine Damming a Glacial lake: The Flatbre Debris Flow Case, Western Norway. Adv. Geosci. 14, 301–307. doi:10.5194/adgeo-14-301-2008

CrossRef Full Text | Google Scholar

Lehmann, F., and Green, A. G. (1999). Semiautomated Georadar Data Acquisition in Three Dimensions. Geophysics. 64, 719–731. doi:10.1190/1.1444581

CrossRef Full Text | Google Scholar

Lovell, H., Livingstone, S. J., Boston, C. M., Booth, A. D., Storrar, R. D., and Barr, I. D. (2019). Complex Kame belt Morphology, Stratigraphy and Architecture. Earth Surf. Process. Landforms. 44, 2685–2702. doi:10.1002/esp.4696

CrossRef Full Text | Google Scholar

Lukas, S. (2007). Early-Holocene Glacier Fluctuations in Krundalen, South central Norway: Palaeoglacier Dynamics and Palaeoclimate. The Holocene. 17, 585–598. doi:10.1177/0959683607078983

CrossRef Full Text | Google Scholar

Lukas, S., and Sass, O. (2011). The Formation of Alpine Lateral Moraines Inferred from Sedimentology and Radar Reflection Patterns: a Case Study from Gornergletscher, Switzerland. Geol. Soc. Lond. Spec. Publications. 354, 77–92. doi:10.1144/SP354.5

CrossRef Full Text | Google Scholar

March, D. W., and Bailey, A. D. (1983). A Review of the Two-Dimensional Transform and its Use in Seismic Processing, First Break. 1, 9–21. doi:10.3997/1365-2397.1983001

CrossRef Full Text | Google Scholar

MATLAB (2018). Version 9.4.0.813654 (R2018a). Natick, Massachusetts: The MathsWorks Inc.

Google Scholar

Maurer, H., and Hauck, C. (2007). Geophysical Imaging of alpine Rock Glaciers. J. Glaciol. 53, 110–120. doi:10.3189/172756507781833893

CrossRef Full Text | Google Scholar

McClymont, A. F., Roy, J. W., Hayashi, M., Bentley, L. R., Maurer, H., and Langston, G. (2011). Investigating Groundwater Flow Paths within Proglacial Moraine Using Multiple Geophysical Methods. J. Hydrol. 399, 57–69. doi:10.1016/j.jhydrol.2010.12.036

CrossRef Full Text | Google Scholar

Midgley, N. G., Cook, S. J., Graham, D. J., and Tonkin, T. N. (2013). Origin, Evolution and Dynamic Context of a Neoglacial Lateral-Frontal Moraine at Austre Lovénbreen, Svalbard. Geomorphology. 198, 96–106. doi:10.1016/j.geomorph.2013.05.017

CrossRef Full Text | Google Scholar

Murray, T., and Dowdeswell, J. A. (1992). Water Throughflow and the Physical Effects of Deformation on Sedimentary Glacier Beds. J. Geophys. Res. 97, 8993–9002. doi:10.1029/92JB00409

CrossRef Full Text | Google Scholar

Neal, A. (2004). Ground-penetrating Radar and its Use in Sedimentology: Principles, Problems and Progress. Earth-Science Rev. 66, 261–330. doi:10.1016/j.earscirev.2004.01.004

CrossRef Full Text | Google Scholar

O'Connell, R. J., and Budiansky, B. (1974). Seismic Velocities in Dry and Saturated Cracked Solids. J. Geophys. Res. 79, 5412–5426. doi:10.1029/JB079i035p05412

CrossRef Full Text | Google Scholar

Palmer, D. (1981). An Introduction to the Generalized Reciprocal Method of Seismic Refraction Interpretation. Geophysics. 46, 1508–1518. doi:10.1190/1.1441157

CrossRef Full Text | Google Scholar

Palmer, D. (1986). Refraction Seismics: The Lateral Resolution of Structure and Seismic Velocity. London, United Kingdom: Geophysical Press. Available at: https://books.google.se/books?id=_Ui4AAAAIAAJ.

Google Scholar

Park, C. B., Miller, R. D., and Xia, J. (2001). “Offset And Resolution of Dispersion Curve in Multichannel Analysis of Surface Waves (Masw),” in 14th EEGS Symposium on the Application of Geophysics to Engineering and Environmental Problems (European Association of Geoscientists & Engineers), 1186–1198. doi:10.3997/2214-4609-pdb.192.SSM_4

CrossRef Full Text | Google Scholar

Park, C. B., Miller, R. D., and Xia, J. (1999). Multichannel Analysis of Surface Waves. Geophysics. 64, 800–808. doi:10.1190/1.1444590

CrossRef Full Text | Google Scholar

Parkes, A. A., Stimpson, I. G., and Waller, R. I. (2013). Geophysical Surveying to Determine the Structure and Origin of the Woore Moraine, Shropshire, UK. Ann. Glaciol. 54, 97–104. doi:10.3189/2013AoG64A112

CrossRef Full Text | Google Scholar

Pegah, E., and Liu, H. (2016). Application of Near-Surface Seismic Refraction Tomography and Multichannel Analysis of Surface Waves for Geotechnical Site Characterizations: A Case Study. Eng. Geology. 208, 100–113. doi:10.1016/j.enggeo.2016.04.021

CrossRef Full Text | Google Scholar

Pellicer, X. M., and Gibson, P. (2011). Electrical Resistivity and Ground Penetrating Radar for the Characterisation of the Internal Architecture of Quaternary Sediments in the Midlands of Ireland. J. Appl. Geophys. 75, 638–647. doi:10.1016/j.jappgeo.2011.09.019

CrossRef Full Text | Google Scholar

Redpath, B. B. (1973). Seismic Refraction Exploration for Engineering Site Investigations. United States. doi:10.2172/4409605

CrossRef Full Text | Google Scholar

Reinardy, B. T. I., Booth, A. D., Hughes, A. L. C., Boston, C. M., Åkesson, H., Bakke, J., et al. (2019). Pervasive Cold Ice within a Temperate Glacier - Implications for Glacier thermal Regimes, Sediment Transport and Foreland Geomorphology. The Cryosphere. 13, 827–843. doi:10.5194/tc-13-827-2019

CrossRef Full Text | Google Scholar

Reinardy, B. T. I., Leighton, I., and Marx, P. J. (2013). Glacier thermal Regime Linked to Processes of Annual Moraine Formation at Midtdalsbreen, Southern Norway. Boreas. 42, a–n. doi:10.1111/bor.12008

CrossRef Full Text | Google Scholar

Reynolds, J. M. (2011). An Introduction to Applied and Environmental Geophysics. 2nd ed. Chichester, United Kingdom: Wiley-Blackwell.

Google Scholar

Ross, N., Brabham, P., and Harris, C. (2019). The Glacial Origins of Relict 'pingos', Wales, UK. Ann. Glaciol. 60, 138–150. doi:10.1017/aog.2019.40

CrossRef Full Text | Google Scholar

Rossi, G., Accaino, F., Boaga, J., Petronio, L., Romeo, R., and Wheeler, W. (2018). Seismic Survey on an Open Pingo System in Adventdalen Valley, Spitsbergen, Svalbard. Near Surf. Geophys. 16, 1–15. doi:10.3997/1873-0604.2017037

CrossRef Full Text | Google Scholar

Roth, M., Holliger, K., and Green, A. G. (1998). Guided Waves in Near-Surface Seismic Surveys. Geophys. Res. Lett. 25, 1071–1074. doi:10.1029/98GL00549

CrossRef Full Text | Google Scholar

Sambridge, M. (1999). Geophysical Inversion with a Neighbourhood Algorithm-I. Searching a Parameter Space. Geophys. J. Int. 138, 479–494. doi:10.1046/j.1365-246X.1999.00876.x

CrossRef Full Text | Google Scholar

Sandmeier, K. J. (2016). GPR and Seismic Data Processing Software - Sandmeier. Available at https://www.sandmeier-geo.de/.

Google Scholar

Sass, O. (2006). Determination of the Internal Structure of alpine Talus Deposits Using Different Geophysical Methods (Lechtaler Alps, Austria). Geomorphology. 80, 45–58. doi:10.1016/j.geomorph.2005.09.006

CrossRef Full Text | Google Scholar

Schön, J. H. (2015). “Elastic Properties,” in Physical Properties Of Rocks: Developments in Petroleum Science. Editor J. H. Schön (Amsterdam, Netherlands: Elsevier), 167–268. doi:10.1016/B978-0-08-100404-3.00006-8

CrossRef Full Text | Google Scholar

Schrott, L., and Sass, O. (2008). Application of Field Geophysics in Geomorphology: Advances and Limitations Exemplified by Case Studies. Geomorphology. 93, 55–73. doi:10.1016/j.geomorph.2006.12.024

CrossRef Full Text | Google Scholar

Sheriff, R. E., and Geldart, L. P. (1983). Exploration Seismology. Volume 1: History, Theory, and Data Acquisition. United States: Cambridge University Press.

Google Scholar

Sollid, J. L., and Bjørkenes, A. (1978). Genesis of Moraines at Midtdalsbreen, South Norway (1:5000). Norsk Geografiske Oppmåling: Norsk Geogrfisk Tidskrift.

Google Scholar

Spagnolo, M., King, E. C., Ashmore, D. W., Rea, B. R., Ely, J. C., and Clark, C. D. (2014). Looking through Drumlins: Testing the Application of Ground-Penetrating Rada. J. Glaciol. 60, 1126–1134. doi:10.3189/2014JoG14J110

CrossRef Full Text | Google Scholar

Steeples, D. W. (2005). Shallow Seismic Methods. Springer Netherlands: Hydrogeophysics, 215–251. doi:10.1007/1-4020-3102-5_8

CrossRef Full Text | Google Scholar

Stokes, C. R., Lian, O. B., Tulaczyk, S., and Clark, C. D. (2008). Superimposition of Ribbed Moraines on a Palaeo-Ice-Stream Bed: Implications for Ice Stream Dynamics and Shutdown. Earth Surf. Process. Landforms. 33, 593–609. doi:10.1002/esp.1671

CrossRef Full Text | Google Scholar

Stokoe, K. H., Wright, S. G., Bay, J. A., and Roesset, J. M. (1994). “Characterization of Geotechical Sites by SASW Method,”. Editor R. D. Woods (ISSMFE Technical Committee #10).

Google Scholar

Tarantola, A. (1986). A Strategy for Nonlinear Elastic Inversion of Seismic Reflection Data. Geophysics. 51, 1893–1903. doi:10.1190/1.1442046

CrossRef Full Text | Google Scholar

Tonkin, T. N., Midgley, N. G., Graham, D. J., and Labadz, J. C. (2017). Internal Structure and Significance of Ice-Marginal Moraine in the Kebnekaise Mountains, Northern Sweden. Boreas. 46, 199–211. doi:10.1111/bor.12220

CrossRef Full Text | Google Scholar

Van Dam, R. L. (2012). Landform Characterization Using Geophysics-Recent Advances, Applications, and Emerging Tools. Geomorphology. 137, 57–73. doi:10.1016/j.geomorph.2010.09.005

CrossRef Full Text | Google Scholar

Van Overmeeren, R. A. (2001). Hagedoorn’s Plus-Minus Method: the beauty of Simplicity. Geophys. Prospect. 49, 687–696. doi:10.1111/j.1365-2478.1964.tb01888.x

CrossRef Full Text | Google Scholar

Wathelet, M. (2011). Geopsy home: Software Applications for Ambient Vibration Techniques. Available at http://www.geopsy.org/.

Google Scholar

Wathelet, M., Jongmans, D., and Ohrnberger, M. (2004). Surface-wave Inversion Using a Direct Search Algorithm and its Application to Ambient Vibration Measurements. Near Surf. Geophys. 2, 211–221. doi:10.3997/1873-0604.2004018

CrossRef Full Text | Google Scholar

Weber, P., Boston, C. M., Lovell, H., and Andreassen, L. M. (2019). Evolution of the Norwegian Plateau Icefield Hardangerjøkulen since the 'Little Ice Age'. The Holocene. 29, 1885–1905. doi:10.1177/0959683619865601

CrossRef Full Text | Google Scholar

Whiteley, R. J., and Eccleston, P. J. (2006). Comparison of Shallow Seismic Refraction Interpretation Methods for Regolith Mapping. Exploration Geophys. 37, 340–347. doi:10.1071/EG06340

CrossRef Full Text | Google Scholar

Willis, I. C., Fitzsimmons, C. D., Melvold, K., Andreassen, L. M., and Giesen, R. H. (2012). Structure, Morphology and Water Flux of a Subglacial Drainage System, Midtdalsbreen, Norway. Hydrol. Process. 26, 3810–3829. doi:10.1002/hyp.8431

CrossRef Full Text | Google Scholar

Woodard, J. B., Zoet, L. K., Benediktsson, Í. Ö., and Iverson, A. (2020). Insights into Drumlin Development from Ground-Penetrating Radar at Múlajökull, Iceland, a Surge-type Glacier. J. Glaciol. 66, 822–830. doi:10.1017/jog.2020.50

CrossRef Full Text | Google Scholar

Woodward, M. J. (1992). Wave‐equation Tomography. Geophysics. 57, 15–26. doi:10.1190/1.1443179

CrossRef Full Text | Google Scholar

Xia, J., Miller, R. D., and Park, C. B. (1999). Estimation of Near‐surface Shear‐wave Velocity by Inversion of Rayleigh Waves. Geophysics. 64, 691–700. doi:10.1190/1.1444578

CrossRef Full Text | Google Scholar

Yilmaz, Ö. (2001). Seismic Data Processing and Seismic Data Analysis. 2nd ed. Tulsa: Society of Exploration Geophysicists.

Google Scholar

York, D., Evensen, N. M., Martı́nez, M. L., and De Basabe Delgado, J. (2004). Unified Equations for the Slope, Intercept, and Standard Errors of the Best Straight Line. Am. J. Phys. 72, 367–375. doi:10.1119/1.1632486

CrossRef Full Text | Google Scholar

Yoshikawa, K., Leuschen, C., Ikeda, A., Harada, K., Gogineni, P., Hoekstra, P., et al. (2006). Comparison of Geophysical Investigations for Detection of Massive Ground Ice (Pingo Ice). J. Geophys. Res. 111, E06S19. doi:10.1029/2005JE002573

CrossRef Full Text | Google Scholar

Zimmerman, R. W., and King, M. S. (1986). The Effect of the Extent of Freezing on Seismic Velocities in Unconsolidated Permafrost. Geophysics. 51, 1285–1290. doi:10.1190/1.1442181

CrossRef Full Text | Google Scholar