Introduction
Alaska glaciers, including the Juneau Icefield of Southeast Alaska, have lost mass rapidly in the late 20th century and early 21st century, as reported using both remote sensing data (e.g. Larsen and others, Reference Larsen2015; Berthier and others, Reference Berthier, Larsen, Durkin, Willis and Pritchard2018) and in situ, mass balance measurements (e.g. Pelto and others, Reference Pelto, Kavanaugh and McNeil2013; McNeil and others, Reference McNeil2020). This mass loss has contributed significantly to global sea level rise (e.g. Gardner and others, Reference Gardner2013) and the ongoing retreat of mountain glaciers affects the hydrology of impacted regions (e.g. Frans and others, Reference Frans, Istanbulluoglu, Lettenmaier, Fountain and Riedel2018; Chesnokova and others, Reference Chesnokova, Baraër, Laperriére-Robillard and Huh2020). Strengthening our understanding of dynamic glacier processes in a warming climate broadens our understanding of present glacier conditions and improves our ability to model future changes. Seismic reflection imaging of the glacier bed – including mapping bed topography and bed conditions –has been useful in elucidating changes and dynamics of other Alaska valley glaciers in the Juneau Icefield (Zechmann and others, Reference Zechmann2018) and elsewhere (e.g. Nolan and Echelmeyer, Reference Nolan and Echelmeyer1999; Babcock and Bradford, Reference Babcock and Bradford2014). In this study, we use active-source seismic reflection and refraction data from Lemon Creek Glacier in the Juneau Icefield to develop a new model of the glacier thickness and bed topography for a portion of the central section of the glacier. Our results will be useful in ongoing studies of subglacial hydrology and subglacial materials.
Lemon Creek Glacier is a small (~9.7 km2 (McNeil and others, Reference McNeil2020)) temperate valley glacier located in a maritime climate (Criscitiello and others, Reference Criscitiello, Kelly and Tremblay2010) at the southern end of the Juneau Icefield (Fig. 1). The glacier's small size belies its glaciological importance; as one of the target glaciers of the long-running Juneau Icefield Research Program, one of the ‘reference glaciers’ of the World Glacier Monitoring Service, and the U.S. Geological Survey ‘benchmark glacier’, Lemon Creek Glacier is the site of the second longest record (1953 to present) of in situ mass balance measurements in North America (Pelto and others, Reference Pelto, Kavanaugh and McNeil2013). Consequently, the glacier's history of retreat and thinning has been thoroughly studied and is well constrained (e.g. Marcus and others, Reference Marcus, Chambers, Miller and Lang1995; Sapiano and others, Reference Sapiano, Harrison and Echelmeyer1998; Criscitiello and others, Reference Criscitiello, Kelly and Tremblay2010; Pelto and others, Reference Pelto, Kavanaugh and McNeil2013; O'Neel and others, Reference O'Neel2019; McNeil and others, Reference McNeil2020). This long-running data set makes Lemon Creek an excellent laboratory in which to explore the links between climate, glacier dynamics, and glacier mass-loss in subarctic mountain glaciers.
Main figure: Landsat image acquired 7 July 2018 (P-058, R-019) showing the location of Lemon Creek Glacier (circled) at the southern end of the Juneau Icefield in Southeast Alaska. Inset: Location of the Juneau Icefield within Western North America. (Landsat 8 image courtesy of the U.S. Geological Survey.)
Prior geophysical study
Despite its long history of glaciological study, there have been a paucity of published geophysical studies focused on quantifying the subglacial topography and ice thickness of Lemon Creek Glacier, particularly in recent decades. The most recent ice-volume (and thus thickness) estimates come from large-scale global or regional studies (Farinotti and others, Reference Farinotti2019). Such studies may be effective in aggregate, but may underestimate ice thickness in the case of small mountain glaciers (Pelto and others, Reference Pelto, Maussion, Menounos, Radić and Zeuner2020).
Prior focused geophysical study has resulted in three published models of the subglacial topography for the whole of Lemon Creek Glacier (Fig. 2). The first model of subglacial topography was presented in Thiel and others (Reference Thiel, LaChapelle and Behrendt1957), and was derived from four gravity transects of the glacier conducted in the summer of 1956 (lines BB′, CC′, DD′ and EE′ in Fig. 3). An update of the Thiel and others (Reference Thiel, LaChapelle and Behrendt1957) model appeared in Miller (Reference Miller1972); in addition to incorporating the gravity transects described by Thiel and others (Reference Thiel, LaChapelle and Behrendt1957), this model included three seismic transects in the upper portion of the glacier (AA′, BB′′ and CC′ in Fig. 3). This model, which we will refer to as the ‘1972 model’, was unchanged from the model of Thiel and others (Reference Thiel, LaChapelle and Behrendt1957) in areas not subject to new data from seismic surveys. The third model – which we will refer to as the ‘1975 model’ – was presented in Miller (Reference Miller1975) and also utilized the gravity transects of Thiel and others (Reference Thiel, LaChapelle and Behrendt1957), additional gravity lines collected in 1971 and 1973 (LA and LB in Fig. 3), and a single seismic line collected in 1968 (SA in Fig. 3).
Two previous models of Lemon Creek's ice thickness as adapted from Miller (Reference Miller1972) (a) and Miller (Reference Miller1975) (b) plotted on the approximate 1972 and 2017 extents of Lemon Creek Glacier. Solid lines are contours of ice thickness in 50 m intervals. In the 1972 model (a), the green dashed line represents the axis of maximum depth and in the 1975 model (b), the dashed lines represent 25 m depth contours. In each panel, the light-gray dotted lines present the 50 m contours of the other model in order to ease direct comparison. The ‘axis of maximum depth’ in the 1972 model is inherited from Thiel and others (Reference Thiel, LaChapelle and Behrendt1957), and represents a line connecting the deepest flow-perpendicular points along the glacier. Lemon Creek Glacier flows in an approximately S to N direction, with the terminus in the NW (top right) corner of the figure. The model of Miller (Reference Miller1972) was digitized from Figure 77 of that study, and the model of Miller (Reference Miller1975) was digitized from Figure 66 of that study. The figure is projected in WGS84/UTM Zone 8 North (EPSG:32608).
Geophysical survey transects as digitized from Miller (Reference Miller1972) (Figure 75 of that work) and Miller (Reference Miller1975) (Figure 65 of that work) compared to the seismic line of this study. BB′, CC′, DD′, and EE′ are gravity transects first described in Thiel and others (Reference Thiel, LaChapelle and Behrendt1957). AA′, BB′′ and CC′are seismic transects described by Miller (Reference Miller1972) from work by Poulter and others (Reference Poulter, Prather and Shaw1967) (AA′ and CC′) and Prather and others (Reference Prather, Schoen, Classen and Miller1968) (BB′′). Lines LA and LB are described by Miller (Reference Miller1975) as additional gravity surveys by Shaw and others (Reference Shaw, Hinze and Asher1972). Miller (Reference Miller1975) places the seismic line of Prather and others (Reference Prather, Schoen, Classen and Miller1968) at the line labeled SA, which is inconsistent with the location given by Miller (Reference Miller1972) (BB′′). We were unable to obtain a copy of the Prather and others (Reference Prather, Schoen, Classen and Miller1968) conference abstract and are unable to reconcile the discrepancy between the locations for that line given in the 1972 and 1975 models. Thiel and others (Reference Thiel, LaChapelle and Behrendt1957) used a slightly different nomenclature for the gravity transects (AA′–DD′) than the nomenclature of Miller (Reference Miller1972) that we use here (BB′–EE′ for the transects of Thiel and others (Reference Thiel, LaChapelle and Behrendt1957)); however the locations of those surveys are consistent amongst Thiel and others (Reference Thiel, LaChapelle and Behrendt1957), Miller (Reference Miller1972) and Miller (Reference Miller1975). The figure is projected in WGS84/UTM Zone 8 North (EPSG:32608).
Qualitatively, the 1972 (Fig. 2a) model suggested basal topography that was quasi-symmetric across the central flow line, with the thickest ice near the southern headwall in the cross-flow center of the glacier. The map view, depth-contour presentation of this model suggested sidewalls that are 50–75% as steep as the immediately adjacent, non-glacierized valley slopes and not ‘U-shaped’ as is typical of glacial erosion. Although Miller (Reference Miller1972) presented cross sections that are more typical of glacial valleys in another figure, these cross sections seem to be inconsistent with the map-view presentation of the bed shown elsewhere in Miller (Reference Miller1972) and Miller and Pelto (Reference Miller and Pelto1999). This model appeared to be the preferred model for more recent publications (e.g. Miller and Pelto, Reference Miller and Pelto1999; Pelto and others, Reference Pelto, Kavanaugh and McNeil2013).
The 1975 model (Fig. 2b) also incorporated the seismic imaging of Prather and others (Reference Prather, Schoen, Classen and Miller1968), as well as additional gravity surveys carried out in 1971 and 1973. While the 1972 and 1975 models agreed somewhat in the range on ice thicknesses expected for Lemon Creek Glacier during the 1950s to 1970s, there are major differences. Most notably, the 1975 model suggested a secondary cirque and overdeepening near the southwestern headwall of the glacier in a region where the 1972 model predicted relatively shallow ice. When compared to the 1972 model, the 1975 model also suggested somewhat steeper sidewalls through the central portion of the glacial trough (with slopes comparable to slopes seen in the non-glacierized portions of the valley), representing a more typical ‘U-shaped’ trough. The 1975 model showed a central trough that shallows as it extended upglacier, while thicker ice in this trough extended to somewhat lower elevations than in the 1972 model.
Several questions arise from the differences in these models and from the details contained in the source publication. There are conflicting reports on the nature of past seismic studies of Lemon Creek Glacier: both Miller (Reference Miller1972) and Miller (Reference Miller1975) agreed that a single seismic line was collected in 1968, with a citation to Prather and others (Reference Prather, Schoen, Classen and Miller1968). However, in the accompanying figures (Figure 75 in Miller (Reference Miller1972) and Figure 65 of Miller (Reference Miller1975)) they provided different locations for this line (Fig. 3: BB′′ from the 1972 model, and SA from the 1975 model), using the same citation of Prather and others (Reference Prather, Schoen, Classen and Miller1968). Furthermore, the two additional seismic lines described by Miller (Reference Miller1972) were not included by Miller (Reference Miller1975), nor did Miller (Reference Miller1975) address this data or why it was excluded from the 1975 model. These contradictions cannot be reconciled from the texts of the earlier reports, and we were unable to obtain copies of the works cited in those seismic studies (Poulter and others, Reference Poulter, Prather and Shaw1967; Prather and others, Reference Prather, Schoen, Classen and Miller1968; Shaw and others, Reference Shaw, Hinze and Asher1972). Additionally, the profile presented in Thiel and others (Reference Thiel, LaChapelle and Behrendt1957) for the gravity transect of CC′ (Fig. 2 of that study, labeled BB′ in their nomenclature) is identical to the profile presented by Miller (Reference Miller1972) (Figure 76 of that study, line CC′), but reported by Miller (Reference Miller1972) to be a seismic profile.
All of the published studies agreed on the locations of land-based gravity studies – originally described by Thiel and others (Reference Thiel, LaChapelle and Behrendt1957) – that have been used by both studies to constrain the thickness and subglacial topography of the lower two-thirds of the glacier. But the two most recent models to utilize those surveys are many decades old produced incompatible models of subglacial topography, even in areas from which both models seemed to be exclusively using data from those gravity transects.
Aside from any discussion of the findings of those prior models, Lemon Creek Glacier has changed significantly since their publication. Although the glacier was already recognized to be in a generally negative mass balance regime at the time of the earlier models (Heusser and Marcus, Reference Heusser and Marcus1964), the annualized rate of mass loss at Lemon Creek Glacier has accelerated since 1980 (McNeil and others, Reference McNeil2020). This ongoing mass loss has been accompanied by a glacier-averaged 37 m decrease in the surface elevation, and a loss of ~20% of its surface area since 1979 (McNeil and others, Reference McNeil2020). This reduced surface area has caused Lemon Creek's tributary branches to separate from the glacier (McNeil and others, Reference McNeil2020). Thus, the 1972 and 1975 models are fundamentally outdated – as would be any model from that time period – and do not represent the glacier in its current configuration. This is clearly illustrated in Figure 2, where the existing models provide depth estimates of up to 100 m in areas that are no longer ice covered as of 2017.
The significant glaciological changes that have affected Lemon Creek over the last 45 years complicate efforts to translate published models of ice thickness into models of subglacial topography. Because of these ongoing changes and the ambiguities and disagreements in published models of ice thickness and subglacial topography, our goal in this paper is to observe the Lemon Creek's bed topography using high-resolution seismic imaging along a ~1 km long flow-parallel line near the center of the glacier. These new data and imaging results will provide clarity regarding the subglacial topography and current ice thickness of the glacier. Our work will support future modeling and geophysical work for which reliable estimates of ice volumes, hydrological gradients and internal stresses are key.
Deployment and data acquisition
In this paper, we present updated bed topography for a section of Lemon Creek Glacier derived from analysis of seismic data collected during June and July of 2017, while the glacier surface was covered with seasonal snow. In addition, we will discuss in some detail the procedures and practices we used in the deployment of the nodes and in our active-source shooting. We hope that a record of the challenges we faced in this deployment as well as the value of some of our attempts to mitigate them will be useful for the planning of future nodal deployments and seismic gun surveys on glaciers.
We deployed 51 Magseis Fairfield Z-Land, Generation 2, 5 Hz three-component nodes (Ringler and others, Reference Ringler, Anthony, Karplus, Holland and Wilson2018), programed to record the vertical component only. The nodes were set to record continuously for the duration of the deployment at 1000 Hz using a pre-set gain of 12 dB. We fitted all of the nodes with aluminum plates about 50 cm in diameter to improve coupling and to mitigate the effects of melt out and subsequent deleveling. We surveyed the station locations (Table S1) at the time of deployment. We surveyed most stations – including all stations discussed in detail in this study–with a Topcon GMS2 Global Positioning System (GPS) unit providing better than sub-meter horizontal accuracy. For a small number of stations surveyed in the later phase of the deployment we used a less accurate GPS (Garmin eTrex20). We report positions and elevations in the WGS84 horizontal datum and ellipsoid, and projected into the UTM Zone 8 North coordinate system.
We used a Betsy seismic gun, provided by the Seismic Source Facility at the University of Texas at El Paso, as the seismic source for this study. The Betsy gun generates a seismic source by firing 8-gauge blanks. We weighted the Betsy gun with two 5 gal (~19 L) jugs of water to improve coupling between with the glacier surface. We initiated shots via hammer strike and timed them using a Verifi-i GPS Synchronizer. We surveyed shot locations (Table S2) at the time of shooting with the Topcon GPS and the same projection as described above.
Our field work was focused on the upper reaches of the glacier and consisted of two phases. In the first phase of the study, we deployed nodes and shot along a 1-km long active-source seismic line incorporating 51 nodal seismometers at 20 m spacing and 52 shot locations at 20 m spacing (Fig. 4, Tables S1, S2). Each of the 52 shot locations featured between one and five shots for a total of 98 shots. Shooting locations between and at the ends of the seismic line resulted in offsets ranging from 10 to 1010 m. Each receiver was buried in snow at a depth of ~0.3 m. The shots were located at the midpoints between the receivers, with one shot off each end of the line. The line was oriented at 352°C (roughly flow-parallel) and near the cross-sectional center of the glacier.
Satellite image of Lemon Creek Glacier showing the locations of seismic nodes for Phases 1 and 2 of this study. Active-source shots for the seismic line discussed in detail in this study (Phase 1) occupied the same footprint as Phase 1 nodes (~5 m). Background image is a high-resolution (0.44 m) orthorectified satellite image provided by the USGS as part of the benchmark-monitoring program at Lemon Creek Glacier (McNeil and others, Reference McNeil2019), dated 28 August 2016. The figure is projected in WGS84/UTM Zone 8 North (EPSG:32608).
The section of the glacier containing our seismic line is an area where the previously published models of the glacier's bed topography differ notably (Fig. 2) and predominantly over a region of the glacier outside of the area directly addressed by previous geophysical transects (Fig. 3). Our seismic line also roughly straddles Lemon Creek's typical equilibrium line altitude (ELA) (Sapiano and others, Reference Sapiano, Harrison and Echelmeyer1998; Mernild and others, Reference Mernild2013; Pelto and others, Reference Pelto, Kavanaugh and McNeil2013) and lies just above an area of increased glacier slope associated with an icefall (Pelto and others, Reference Pelto, Kavanaugh and McNeil2013).
For the second phase of the study, we redeployed 41 of the nodal seismometers across the upper half of the glacier. We redeployed 34 of these nodes as part of 7 cross-flow seismic lines with station spacings of 300–330 m, and along-flow line spacings of ~400 m at elevations ranging from 1050 to 1200 m (Fig. 4). We releveled and reburied ten stations from the phase-one seismic line in the same positions that we deployed them in during the first phase of the study. These remaining ten stations had station spacings of 100 m, plus one at 150 m. After we redeployed the nodes, we left them to record passive seismicity for 17 days. After 17 days, we made an additional 79 shots at 25 locations, including 47 shots at 14 locations along the cross-flow lines, and 32 shots at the 11 locations still occupied by nodes along the phase-one seismic line. Many of the nodes had experienced significant melt out during the interval between the two phases of shooting and had to be releveled.
Data analysis and results
We focus this paper on analysis of the data from phase one of our deployment and data collection: the densely spaced 1-km seismic line. We completed initial analyses using the data from the later phases of the deployment, including depth-to-bed picks at near offsets for the distributed shots and sparse, long-offset CMP imaging for the shots along the seismic line in July. However, our near offset picks were unreliable due to source noise on near-zero offset shots (discussed further later), and the later long offset shots did not add significant details to our common midpoint (CMP) image. Future work may glean more information from these data.
In our studies of the data collected along our dense, 1-km active-source seismic line we focus on two main analyses. We first used the arrivals from our shots in order to perform a joint refraction/reflection tomographic inversion to derive a P-wave velocity model using the Tomo2d software package (Korenaga and others, Reference Korenaga2000). We then used that tomographic model in processing our data into a CMP reflection image that provides a record of ice thickness along our seismic line.
Velocity modeling
Because the depth of Lemon Creek Glacier is not well constrained, and because there is a direct trade off between reflector depth and seismic velocity (both in general and in the parameterization in Tomo2d), we used an initial ‘brute stack’ CMP reflection image with minimal processing to obtain an initial estimate of the thickness of the glacier ice (Fig. 5a).
Comparison of three CMP stack images of Lemon Creek Glacier: (a) Automatic gain control and depth migration has been applied, but no further processing has been done. (a) Application of a Stolt Migration using the initial velocity model as shown in Table 1. (c) Final CMP image after application of our final velocity model (as shown in Table 1). In these profiles, 0 m identifies the location of the farthest north/downglacier Betsy gun shot. The vertical-axis elevations are relative to the WGS84 ellipsoid. These images have a vertical exaggeration of 1.5×.
We generated travel times by manually picked arrivals of 585 reflected and 1737 refracted P-waves to perform tomographic velocity modeling across the portion of the glacier covered by our seismic line. We used the Tomo2d software package for our modeling (Korenaga and others, Reference Korenaga2000). Tomo2d allows for joint reflection and refraction tomography over a user-defined mesh and reflective layer, and generates models through perturbation and minimization of travel time misfits. Although Tomo2d allows for joint inversion for both reflector depth and velocity, we defined our inversion settings such that perturbations to the initial basal reflector were smaller than perturbations to velocity. We made this choice because we believed the basal reflector from our brute-stack imaging was more accurate than our initial assumed seismic velocity, and we found in sensitivity tests that our inversion parameters performed better in resolving velocity than reflector depth.
We performed our velocity modeling over a model space that was discretized horizontally at 100 m intervals and vertically at 25 m intervals. Additionally, we parameterized the model to include an additional mesh-node at 5 m depth in order to better account for very low velocities in the seasonal snow layer near the glacier surface. Our large horizontal grid spacing was designed in order to maximize the utility of our somewhat sparse and predominantly vertical-traveling reflected-arrival coverage (Fig. 6), while still allowing us somewhat fine vertical resolution. Our initial model is shown in Table I and was informed by the previously published velocity estimates for Lemon Creek Glacier from Miller (Reference Miller1972).
Example shot gathers and modeled ray paths for three shooting locations along our dense seismic line. In keeping with previous figures, 0 m is defined as the NNE (downglacier) end of the seismic line. Shots are located at 1010 m (left panels), 420 m (center panels) and 0 m (right panels) along the seismic line. Shot gather traces have been normalized by the RMS amplitude of each, have not been filtered, and have not been corrected for instrument response. Raypaths are shown for only the first-arriving example of each phase, and are shown without vertical exaggeration. Ice (blue) and bedrock/till (dark gray) are based on the results of this study. Dashed lines within the ice indicate ice depths of 100–300 m
Comparison of initial and final model for P-wave seismic velocity modeling
Initial velocity model is based on previously published velocity findings of Miller (Reference Miller1972).
Final 1-D velocity model and standard deviation are calculated via weighted averaging of a 2-D velocity model, where weights are assigned based on ray coverage of each grid-node.
Although Tomo2D allows for joint inversion of reflected and refracted P-wave arrivals, we first modeled the shallow subsurface of the glacier using refracted arrivals and then modeled the reflected arrivals separately. We found that this separate modeling produced more stable results, in particular in the shallowest regions of the glacier where we observed both the highest vertical velocity gradients and the highest potential for lateral variations as well. For both models, we saw significant root mean square (RMS) misfit reduction: from 3.1 to 1.2 ms for the refracted arrivals, and from 6.1 to 3.3 ms for the reflected arrivals. For our final model, we combined the results of these two models (Table 1).
We compiled the results of our tomographic modeling into the one-dimensional (1-D) model shown in Table 1. The 1-D model we present is a result of a weighted average of a 2-D model, where overall ray coverage determined the weight assigned to each node, resulting in the final model and standard deviation of velocities in each layer of the 2-D model is shown in Table 1. The variability within the 2-D model is distributed throughout the model, and the model lacks any obvious horizontal features that would seem to imply changes in the structure of the glacier. We also smoothed vertical velocity variations below the maximum penetration depth of our refracted arrivals (~75 m), as ray coverage below those depths consists only of near-vertically traveling reflected phases. Our velocity model is not significantly different than earlier published estimates for the seismic velocity of Lemon Creek Glacier (Miller, Reference Miller1972).
We utilized the final velocity model (Table 1) in our further reflection imaging (CMP stacking) of the glacier. Our processing and interpretations of this imaging ultimately resulted in our final estimates of ice thickness and bed topography.
Reflection imaging
We based our reflection imaging on recordings from the 51 receivers (Table S1; stations 1007–1057) and 52 shot locations (Table S2; shotpoints 7–58) at elevations ranging from ~1080 to ~1150 m (WGS-84 elevations) that formed phase one of our summer 2017 deployment. We performed initial quality control to identify shots with either timing errors or contamination by other seismic signals. We stacked 83 (of 98) shots that passed quality control checks by shotpoint for further processing. We stacked and processed our data at the full recorded sample rate of 1000 Hz. Examples of several representative stacked shot gathers and ray coverage are shown in Figure 6.
After stacking, we performed further processing steps using the Landmark SeisSpace software suite. We defined our CMPs every 10 m along the length of the line for 101 CMP stacking locations (Table S3). We first applied a Butterworth bandpass filter from 5 to 150 Hz to our stacked shot gathers and then applied an automatic gain control with an operating window of 400 ms. For our initial unmigrated and migrated stacks, we applied normal moveout corrections (NMO) using a constant velocity of 3700 ms−1 and then performed a time to depth conversion using the same constant velocity.
Our initial, unmigrated brute-stack CMP image reveals multiple intersecting reflectors generating the first reflected arrivals (Fig. 5a). We interpret the first reflector(s) as reflections off the ice-bed interface, this reflector contains a number of ‘bowtie’ structures that are likely the result of a depression in the topography of the reflector (Sheriff and Geldart, Reference Sheriff and Geldart1995). We believe the later arriving reflectors to be reflections of the sidewalls of the glacial trough. Our image lacks any interpretable structure above the basal reflector.
We performed further processing on the unmigrated image in order to remove the ‘bowtie’ structures by applying a basic Stolt migration (Stolt, Reference Stolt1978), which we applied using the SeisSpace processing software to generate the CMP in Figure 5b. We used this image to obtain an initial estimate of ice thickness to use as a constraint for our tomographic modeling.
After deriving a velocity model based on refractions and reflections as described in the previous section we used the updated velocity model for NMO correction and time-to-depth conversion for the reflection image. In the resulting CMP image, the signal-to-noise ratio of the stacked traces is improved. After migrating this latter image with a Stolt migration (Fig. 5c), we were able to clearly image the basal topography of the glacier. Our picks of the basal reflector in the final CMP image, with an estimated error of 3 ms in our depth-to-bed picks. This translates to a potential error in our depth picks of ±5 m.
Our imaging reveals the presence of a significant overdeepening in the glacier's subglacial topography that is defined by two step-wise decreases in the elevation of the bed (Fig. 7, Table S3). From a riegel at ~150 m along-track distance, the elevation of the bed decreases by ~40 m, reaching an initial low-point at ~300 m along track. Following this initial low, the elevation of the bed remains mostly flat until ~500 m along track. At this point, the elevation of the bed begins to decrease again, and decreases by additional ~75 m to reach its nadir at ~775 m along track. The bed subsequently increases in elevation slightly to the end of the seismic line at 1010 m along track. Within this overdeepening, the glacier thickens from ~210 m at the riegel to ~250 m after the initial decrease in the basal elevation. The ice remains at this thickness for ~200 m along track until the second decrease in basal elevation. Upglacier of this second decrease the glacier thickens further, reaching a maximum thickness of 359 m at the nadir of the basal topography. When compared to our brute-stack CMP (Fig. 5b), our velocity model enhances the visibility of the shallower, flatter portion of the overdeepening and extends the down-glacier edge of the deepest portion of the overdeepening to the north (in the down-glacier direction).
Top panel (a): Depth to basal reflector from this study overlain on the glacier's 2017 extent. Glacier flows from South to North. Bottom panel (b): Profile view of the Lemon Creek's basal topography as found in this study (solid) and earlier studies (dashed) plotted without vertical exaggeration. Solid circles (1972 model) and solid squares (1975 model) indicate the locations where our seismic line intersects bed-depth contours for the previous models. Open circles (1972 model) and squares (1975 model) with question marks indicate the locations where our seismic line intersects the likely maximum depth of each model. This ‘axis of maximum depth’ is explicitly defined in Miller (Reference Miller1972). For the Miller (Reference Miller1975) model, it is assumed to be the midpoint between the deepest indicated contours. The depth at these points is assumed to be just less than the value that would have required an additional depth contour in the source model. Basal topography for 2017 is based on the ice thickness we derive from seismic-reflection imaging and a glacier surface derived from our surveyed station and shot elevations. Basal topography for the earlier models is based on digitized and georeferenced isopach depth maps of those models as well as a glacier surface derived from a 1974 DEM publicly available as part of the USGS Benchmark Glaciers program (McNeil and others, Reference McNeil2019). The figure is projected in WGS84/UTM Zone 8 North (EPSG:32608), and elevations are displayed relative to the WGS84 ellipsoid.
Discussion
Data collection
In assessing our deployment and data collection procedures, we report positive outcomes for the use of the Betsy gun. We were satisfied with the effectiveness of the Betsy gun as a seismic source for use on a midsummer, wet-snow-covered mountain glacier. Refracted and first-reflected P-waves were visible across the majority of our 1-km seismic line with sufficient signal-to-noise ratio that we were able to pick them reliably in geophone recordings. Although we did not experience any of the pore-fluid effects reported by an earlier study using a Betsy gun on a glacier (Baker and others, Reference Baker2003), at near offsets we did find the traces to be somewhat noisy, which complicated our efforts by obscuring near-zero offset reflections and made our attempts to pick vertical reflections at distributed shotpoints unreliable. The noise at these near-offset shots was centered at ~175 Hz. Although relatively narrow band, the signal created a ‘ringing’ in the traces that obscured reflected arrivals, even after careful filtering.
Additionally, our later distributed shooting suggests that in this environment, Betsy gun sources are not recorded at nodes with offsets of >1 km. However, the portability, cost and potential ease of permitting for the Betsy gun make it an attractive option as a source for small-scale active-source studies of glaciers.
The Fairfield nodes performed well in the cool, extremely wet environment of the glacier when buried beneath the snow. In addition to the water-saturated environment of the seasonal snow, the nodes experienced heavy rainfall and/or 100% humidity on many days of the deployment. While buried, the nodes’ internal temperature sensors recorded diurnal variations of 1°C to 2°C; these variations increased to as much as 12°C after the nodes began to melt out.
All of the nodes melted out of their shallow snow burial after 10–14 days of deployment. Some of these nodes became off level after melt out and tilted by up to 45°C; although we did not take systematic measurements of all of the tilted sensors, typical tilts were much less. Prior to the second phase of shooting, we releveled all of the nodes that remained along the 1 km seismic line. We observed that the snow plates appeared to cause extra tilt beyond what may have occurred without the snow plates, perhaps due to heat conduction facilitated by the aluminum plates. Consequently, we do not recommend using thermally conducting snow plates on node deployments where snow or ice melt out are potential concerns.
P-wave velocity model
Our P-wave velocity model is in large part consistent with the previously published velocities of Miller (Reference Miller1972) (Table 1). Our final model is fairly robust, as varying the initial models by ±5% results in functionally similar final models. Our initial model uses a slightly lower velocity than that of Miller (Reference Miller1972) in order to ensure that our modeling is sensitive to the observed travel times and that we are truly resolving velocity variations at the depths displayed in Table 1.
In our 2-D model we see the greatest lateral variability at a depth of 5 m, as represented in the standard deviation for each layer in Table 1. This lateral variation may be a result of real variations in the thickness and/or the character of the summer snow, the shallow ice structure, or the increased sensitivity of this layer to short-offset arrivals for which small absolute-errors in travel time estimates are fractionally larger. However, this variability is distributed in such a way that does not suggest any clear physical or structural changes and neither our survey nor our velocity modeling was optimized to image this layer of the glacier.
Seismic imaging and updated depth estimates
The depth and shape of the basal reflector we map differs significantly from the previously published basal topography of both the Miller (Reference Miller1972) and Miller (Reference Miller1975) models are shown in Figure 2. Most notably, both models underestimate the maximum thickness of the ice by as much as 150 m (40% of the overall thickness; Fig. 7). Both the 1972 and 1975 models show a maximum depth of slightly more than 200 m. The 1972 model has contours at 50 m intervals, for a maximum ice thickness of <250 m, while the 1975 model has contours at 25 m intervals, for a maximum ice thickness of <225 m. This contrasts with the maximum thickness we find of 359 m. This difference in ice thickness is sufficiently large that it is exceedingly unlikely to be the result of any differences in migration, velocity models, or processing of the seismic data, and it is particularly surprising given of thinning the glacier has experienced in recent decades (Sapiano and others, Reference Sapiano, Harrison and Echelmeyer1998; Criscitiello and others, Reference Criscitiello, Kelly and Tremblay2010). Indeed, a comparison of our surveyed elevations (Tables S1, S2) to the elevations from a 1974 digital elevation model of the glacier (McNeil and others, Reference McNeil2019) provide a mean elevation reduction of ~35 m along the length of our seismic line. This reduction in elevation is an order of magnitude higher than the expected uplift from ongoing glacio-isostatic adjustment in the region (Larsen and others, Reference Larsen, Motyka, Freymueller, Echelmeyer and Ivins2005), thus we expect to find ~35 m of thinning at the location of our shot profile since the time of the earlier geophysical surveys, not the significantly thicker ice that we observe.
Except for the divergent estimates of maximum depth, the 1972 model is more qualitatively consistent with the depths we observe across the area covered by our seismic line. At the downglacier end of our seismic line, the predicted depth of the 1972 model is consistent with our observations, and the predicted deepest point along our seismic line is somewhat similar, probably falling within the errors inherent in digitizing and georeferencing that model from a journal figure. However, it is surprising that this model seems to have missed the significant overdeepening we see at the upglacier end of our line, as this is an area that Miller (Reference Miller1972) reported to have covered with a seismic survey (Fig. 3, CC′) in addition to the earlier gravity survey of Thiel and others (Reference Thiel, LaChapelle and Behrendt1957) along that same line.
The 1975 model diverges more starkly from what we observe across our seismic line. In addition to underestimating the thickness of the ice by ~150 m, the 1975 model proposes that the thickest section of the glacier extends downglacier of the area where we see the thickest ice, thus missing the overdeepening we observe in shape as well as depth. The 1975 model also proposes shallowing of the bed across the area where we see a significant deepening of the bed.
Our imaging of Lemon Creek Glacier also differs from ice thickness estimates recently published for Lemon Creek as part of a broader, global effort (Farinotti and others, Reference Farinotti2019) by at least 100 m over the deepest portion of the glacier. The thickness estimates of Farinotti and others (Reference Farinotti2019) use multiple flow-modeling procedures to estimate thickness based on glacier surface topography and climate. That our thickness estimate differs from those results suggests that Lemon Creek Glacier contains complexities that are not well modeled through such broad-scale modeling efforts. We therefore infer that more site-specific data is important for reliable modeling of mountain glaciers.
Our study captures the basal topography of a small-section of Lemon Creek Glacier, but suggests a thicker glacier with more varied bed topography than previous models have captured. This finding is important and has implications for geophysical and climatological studies of the glacier. Most simply, thicker ice means a larger ice volume and affects estimates of future and evolving glacier dynamics. This is particularly important on a glacier with a well recorded history of negative mass balance as thicker ice implies a slower response to climate forcings (Harrison and others, Reference Harrison, Elsberg, Echelmeyer and Krimmel2001). If this is indeed the case, the previously observed record of negative mass balance at the glacier suggests that the glacier may be further out of equilibrium than previously suggested, and that future losses may be greater than previously predicted.
From the perspective of glacier dynamics; thicker ice and the presence of retrograde bed topography affect hydrologic gradients (e.g. Hubbard and Nienow, Reference Hubbard and Nienow1997), the likelihood of channelized versus distributed drainage configurations (Flowers, Reference Flowers2008; Vore and others, Reference Vore, Bartholomaus, Winberry, Walter and Amundson2019), and subglacial sediment distribution (Creyts and others, Reference Creyts, Clarke and Church2013). Thus studies of the subglacial hydrology of the glacier require reliable models of bed topography and ice thickness, and the differences between the bed topography we observe and the bed topography in earlier models are significant in such hydrological studies. Furthermore, basal sediment distributions are significant in a wide range of basal dynamic processes (Cuffey and Paterson, Reference Cuffey and Paterson2010).
Finally, independent measurements of glacier surface velocity during 2017 at a site along the seismic profile revealed relatively steady 4.3 cmd−1 of motion over 1.5 months midsummer (Bartholomaus and others, Reference Bartholomaus2019). Ice thickness at this location, according to our measurements, is ~300 m, which, according to Glen's flow law, would lead us to expect a deformation velocity of 4.9 cmd−1 using a parabolic shape factor to account for traction along the glacier side walls (Cuffey and Paterson, Reference Cuffey and Paterson2010). Thus, our current bed mapping implies that, during the majority of the summer, no basal motion occurs at this site along the seismic profile. However, if one assumed that the earlier bed models were correct, ice thickness of only ~ 190 m at this site would be expected, implying deformation velocities of 1.2 cmd−1. Acceptance of the earlier ice thickness maps would then lead us to assume a dominant contribution of background, hydrologically invariant basal motion to the total, measured surface motion. The new bed mapping reported here clearly has a critically important impact on interpretations of the dynamics of Lemon Creek Glacier.
Conclusions
We present a new velocity model and seismic reflection image of the bed of Lemon Creek Glacier in the Juneau Icefield, Alaska. Our field acquisition techniques combining a seismic gun and nodal seismometers make for a cost-effective and flexible method of performing active source surveys on steep, wet mountain glaciers.
Existing models of the subglacial topography for Lemon Creek Glacier have significant disagreements and represent the glacier prior to ~50 years of thinning, retreat and area reduction associated with persistently negative mass balance. Thus, the existing models of the glacier represent the glacier in a significantly different glaciological context than the current situation. We compare our new model and images in detail with previous geophysical results. The model of Miller (Reference Miller1972), which has been used by recent publications on Lemon Creek (e.g. Pelto and others, Reference Pelto, Kavanaugh and McNeil2013), is more compatible with our imaging in qualitative shape, but underpredicts the thickness of the glacier across our line by as much as 150 m (~40% of the overall thickness of the glacier), as well as overestimating the elevation of the subglacial topography. It also lacks overdeepenings and areas of retrograde bed slope, all of which have significant glacier dynamic implications.
Lemon Creek Glacier is the site of a long-running glaciological data collection project, and accurate estimates of thickness are required for analysis of future and historical glaciological data. Thus, our new results will inform and improve the analyses and models of the glaciers long-term behavior by improving our knowledge of the glacier's basal topography and characteristics. Our data shows a much (up to 150 m) thicker glacier than previous models. This represents a considerable difference, and significantly affects the dynamics and ice volume estimates of the glacier, as well as understanding the past and future evolution of the glacier in light of its ongoing thinning and retreat. Our new profile allows us to infer that, during the majority of the summer, basal motion does not contribute appreciably to the velocity measured at the glacier surface.
Our imaging differs sufficiently from existing models such that we believe further radar or seismic imaging should be employed in order to obtain an updated model of the bed topography and ice thickness of Lemon Creek Glacier. This glacier is the site of long-running, high-quality glaciological observations, and is well suited as a case study for many studies of glacier geophysics; however, many of these methods require more expansive, reliable and accurate estimates of ice thickness than currently exist. Consequently, we suggest Lemon Creek as a priority for expanded geophysical imaging to better constrain ice thickness and subglacial topography across the whole of the glacier.
Supplementary Material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2021.32.
Acknowledgements
We thank the Juneau Icefield Research Program (JIRP) for their support with facilities at Camp 17 and logistical planning for the field work. We thank Emily Graves, Joachim Schalk and Celeste Labedz for their participation in the field work. We thank the Incorporated Research Institutions for Seismology Seismic Source Facility at the University of Texas at El Paso for their assistance with the seismic sources. We thank two reviewers for their helpful commentary, and the Journal's editorial volunteers and staff for their invaluable role in supporting this publication. This work was partially funded by Marianne Karplus's startup funds at the University of Texas at El Paso.
Open access
