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.

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.

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.

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.

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.

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 λ = v p h a s e ( f ) f (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.

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.

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

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

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.

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

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.

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.

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.

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.

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

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

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.

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

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.