1. Introduction and motivation

The structure and properties of firn are important to understand for numerous glaciological applications (e.g. Kinar and Pomeroy, Reference Kinar and Pomeroy2007; Diez and others, Reference Diez2014; Schlegel and others, Reference Schlegel2019). For ice shelves, the loss of firn air content is linked to changes in surface elevation that could be misdiagnosed as evidence of basal melting (Holland and others, Reference Holland2011), and firn densification has been implicated in catastrophic ice-shelf collapse (Glasser and others, Reference Glasser and Scambos2008). Effective geophysical characterisation of the subglacial environment of ice masses requires firn effects to be removed, e.g., to correct observations of basal seismic and/or radar reflectivity for energy loss through firn layers (e.g., Kulessa and others, Reference Kulessa2017; Zechmann and others, Reference Zechmann2018).

Ice slabs form in firn where surface melt occurs, including the coastal regions of Antarctica and Greenland. With a lateral extent of tens-to-hundreds of metres, ice slabs make the shallow firn column impermeable (Benson, Reference Benson1962; MacFerrin and others, Reference MacFerrin2019, Miller and others, Reference Miller2022) and can increase its local density from between 400–800 kg m^–3 (typical of firn) to that of pure glacier ice, 830–917 kg m⁻³ (e.g. Hubbard and others, Reference Hubbard2016; MacFerrin and others, Reference MacFerrin2019; Culberg and others, Reference Culberg, Schroeder and Chu2021). Without explicitly accounting for the effects of ice slabs, meltwater runoff can be underestimated: in Greenland, for example, regional climate models suggest runoff to be underestimated by almost 60% if ice slab formation is excluded (MacFerrin and others, Reference MacFerrin2019).

Borehole sampling can determine the thickness of firn layers (e.g., Hubbard and others, Reference Hubbard2016) but geophysical techniques facilitate measurement away from borehole control and can constrain other physical properties of the firn column. Several seismic techniques have been developed for firn characterisation (Kirchner and Bentley, Reference Kirchner and Bentley1979; King and Jarvis, Reference King and Jarvis1991, Reference King and Jarvis2007; Booth and others, Reference Booth2013; Hollmann and others, Reference Hollmann, Treverrow, Peters, Reading and Kulessa2021) but these make assumptions (e.g., that firn density always increases with depth) that limit their applicability. In this study, we report on the application of acoustic Full Waveform Inversion (FWI) techniques (Virieux and Operto, Reference Virieux and Operto2009) as a means of improving seismic surveying for constraining the physical properties of firn. We use synthetic seismic data to highlight promising results for the FWI approach and establish the next steps that would broaden the applicability of the method for firn studies.

2. Seismic methods to model FIRN and ice

Firn structures are detectable with seismic refraction techniques since, as firn is compacted and densified, its seismic velocity increases with depth. Seismic velocity models can therefore be linked, via empirical models, to firn density. A commonly used approach for characterising the seismic velocity through firn from controlled- source seismic Data are Herglotz-Wiechert inversion (HWI) (Herglotz, Reference Herglotz1907; Wiechert, Reference Wiechert1910; Slichter, Reference Slichter1932; Nowack, Reference Nowack1990). HWI obtains a velocity-depth model by considering the slowness (the reciprocal of velocity) of refracted seismic arrivals (e.g., Thiel and Ostenso, Reference Thiel and Ostenso1961, Rege and Godio, Reference Rege and Godio2011; Diez and others, Reference Diez, Eisen, Hofstede, Bohleber and Polom2013). HWI is practical for many firn applications since its key assumption – that seismic velocity increases gradually with depth – is fulfilled by many steady-state firn profiles. However, this assumption is violated for more complex firn profiles, including those which include so-called ‘ice slabs’ (locations in which water has infiltrated and refrozen). In such cases, the ice slab represents a high velocity anomaly, with a velocity reduction likely at its lower interface on returning to unmodified firn. Consequently, when HWI is applied to firn profiles containing infiltration ice, the boundaries of ice slabs are improperly represented, leading to errors in the final velocity model and any derivative density estimate.

2.2 Progress on FWI development

We have been exploring the use of FWI to recover the extent of an ice slab located within a firn column. Synthetic firn velocity profiles were generated using the Herron-Langway (HL) (Herron and Langway, Reference Herron and Langway1980) accumulation model with the critical, pure ice and surface snow densities set as 550, 917 and 400 kg m^–3 respectively, a 10 m depth temperature of −30°C and an accumulation rate of 0.2 m w.e. a⁻¹. These parameters predict that firn density will increase from 400 to 830 kg m^–3; and continues to increase to a maximum ice density, here considered as 917 kg m^–3 at which the firn-ice transition would occur at a predicted depth of 200 m. This density profile is used to model seismic velocity with the Kohnen (Reference Kohnen1972) empirical expression. To represent the ice slab, velocities in the depth interval 40–60 m are increased to 3800 m s⁻¹. Throughout this study we consider the synthetic velocity model with ice as the ‘True’ model; models implied by HWI and FWI are benchmarked against this reference model. The FWI algorithm models the acoustic wave equation with the finite-difference (FD) method, due to its simplicity and efficiency compared to other techniques available to solve partial differential equations (Virieux and Operto, Reference Virieux and Operto2009; Zhang and Yao, Reference Zhang and Yao2013). To simplify the source characterisation in the FD approach the source wavelet is modelled using a Ricker wavelet with a peak frequency of 60 Hz. To ensure modelling stability, data are recorded with a time sampling of 0.001 s, for 1 s of propagation (Courant and others, Reference Courant, Friedrichs and Lewy1967). In our FWI, we use a gradient approach to define the fit between recorded and modelled data. The objective function (OF) used is a least-squares (LS) formulation that minimises the sum of the square of the difference between the observed (d) and predicted (p’) datasets:

(1)$$f = \displaystyle{1 \over 2}\vert {\boldsymbol {\,p}^{\prime}}-{\boldsymbol d}\vert _2^2 $$

We use the Madagascar framework (Fomel and others, Reference Fomel, Sava, Vlad, Liu and Bashkardin2013), provided by Center for Wave Phenomena at Colorado School of Mines, to implement this approach (Aragao and Sava, Reference Aragao and Sava2020). To minimise the impact of cycle skipping, we start iterations for the low frequency wavefield, filtered between 3 and 10 Hz, until the objective function updates suggest model convergence. Thereafter, the frequency content is progressively increased in 10 Hz bands, to a maximum of 60 Hz.

4. Future research priorities

The presence of ice slabs influences drainage and meltwater runoff across glaciers and ice masses (MacFerrin and others, Reference MacFerrin2019). In the case of ice shelves, this meltwater increase and decrease in permeability can lead to a reduction in ice shelf stability (Munneke and others, Reference Munneke, Ligtenberg, Van Den Broeke and Vaughan2014). This process was observed on Larsen B Ice Shelf, where firn compaction, meltwater ponding and hydrofracturing were strongly implicated in the shelf's rapid disintegration in 2002 (Scambos and others, Reference Scambos, Bohlander, Shuman and Skvarca2004). Hubbard and others (Reference Hubbard2016) detected a 40 m thick ice slab in the upstream reaches of Larsen C Ice Shelf (LCIS) in borehole Optical Televiewer (OPTV) data, interpreted as being the accumulation of episodic refreezing of meltwater ponds. Here, the firn zone is 10°C warmer and 170 kg m^–3 denser than undisturbed firn in the surrounding area. Regional geophysical surveys suggested that the ice slab is at least sixteen kilometres across and several kilometres long; while GPR surveys could constrain thickness, and seismic velocity models (Kulessa and others, Reference Kulessa2019) were consistent with pure ice, neither method could establish both the full depth extent and velocity anomaly of the slab. FWI methods show promise for this application and offer a means of monitoring the development of similar processes at other sites. This potential therefore motivates the future application of FWI to similarly stratigraphically-complex areas of firn, such as interior areas of the Greenland Ice Sheet (MacFerrin and others, Reference MacFerrin2019) and LCIS.

FWI was applied to legacy field datasets from Antarctica's Pine Island Glacier, but none are currently compliant with a stable inversion. Hence, the promising results from synthetic trials reported herein have not yet been validated for field data. The high frequency components of explosive sources used within glaciology lead to cycle skipping when the starting model is not close enough to the true model. Field data have additional difficulties vs synthetic data for FWI, due to the assumption of an acoustic wavefield. The effect of neglecting the elastic component in field data leads to a loss in resolution and accuracy in the recovered velocity model (Agudo and others, Reference Agudo, Da Silva, Warner and Morgan2018). Attenuation was not accounted for, even though firn is known to be highly attenuative and has an effect on the wavelet with offset (e.g. Eisen and others (Reference Eisen2010), King and Jarvis (Reference King and Jarvis1991)).

To obtain a compatible field dataset, particular care must be taken to ensure that the source wavelet is sufficiently consistent and rich in low-frequency content, to facilitate stable FWI. Many examples of terrestrial field data examples exist that satisfy these criteria, allowing for successful FWI and thus improved near-surface velocity models (e.g. Adamczyk and others (Reference Adamczyk, Malinowski and Malehmir2014), Borisov and others (Reference Borisov, Gao, Williamson and Tromp2020), Irnaka and others (Reference Irnaka, Brossier, Métivier, Bohlen and Pan2022)). The extension of this practice into glaciology motivates the acquisition of FWI-compliant data over an ice slab target, and LCIS offers a promising opportunity for this given its apparent complexity and the wealth of data already available for validation. FWI over an area such as this would provide greater constraints on the structure of the firn column, allowing for a more comprehensive reconstruction of seismic properties than currently feasible with existing seismic techniques. If our acoustic FWI algorithm is deemed to be effective, it could in principle undergo further development to recover the anisotropic and/or elastic properties of the firn and deeper ice column (Li and Alkhalifah (Reference Li and Alkhalifah2022), Kan and others (Reference Kan, Chevrot and Monteiller2023)).