A synthetic study of acoustic full waveform inversion to improve seismic modelling of firn

Abstract

The density structure of firn has implications for hydrological and climate modelling and for ice shelf stability. The firn structure can be evaluated from depth models of seismic velocity, widely obtained with Herglotz-Wiechert inversion (HWI), an approach that considers the slowness of refracted seismic arrivals. However, HWI is appropriate only for steady-state firn profiles and the inversion accuracy can be compromised where firn contains ice layers. In these cases, Full Waveform Inversion (FWI) can be more successful than HWI. FWI extends HWI capabilities by considering the full seismic waveform and incorporates reflected arrivals, thus offering a more accurate estimate of a velocity profile. We show the FWI characterisation of the velocity model has an error of only 1.7% for regions (vs. 4.2% with HWI) with an ice slab (20 m thick, 40 m deep) in an otherwise steady-state firn profile.

1. Introduction and motivation

The structure and properties of firn are important to understand for numerous glaciological applications (e.g. Kinar and Pomeroy, 2007; Diez and others, 2014; Schlegel and others, 2019). 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, 2011), and firn densification has been implicated in catastrophic ice-shelf collapse (Glasser and others, 2008). 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, 2017; Zechmann and others, 2018).

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, 1962; MacFerrin and others, 2019, Miller and others, 2022) 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, 2016; MacFerrin and others, 2019; Culberg and others, 2021). 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, 2019).

Borehole sampling can determine the thickness of firn layers (e.g., Hubbard and others, 2016) 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, 1979; King and Jarvis, 1991, 2007; Booth and others, 2013; Hollmann and others, 2021) 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, 2009) 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, 1907; Wiechert, 1910; Slichter, 1932; Nowack, 1990). HWI obtains a velocity-depth model by considering the slowness (the reciprocal of velocity) of refracted seismic arrivals (e.g., Thiel and Ostenso, 1961, Rege and Godio, 2011; Diez and others, 2013). 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.1 FWI motivation

The limitations of HWI can potentially be overcome with FWI. Where HWI only considers the slowness of first-arrivals, FWI considers both slowness, amplitude and phase information for refracted arrivals (beyond just the first arrival) and reflected wavelets. In so doing, it iteratively improves the match between the recorded and model seismic data predicted by an evolving distribution of subsurface seismic properties. The incorporation of reflected arrivals into FWI also means that it is, in principle, sensitive to more complex velocity structures than can be detected by HWI. For the case of firn structures, this means that FWI should be able to reconstruct the extent of ice slabs by considering reflections from their bounding interfaces.

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, 1980) 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 (1972) 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, 2009; Zhang and Yao, 2013). 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, 1967). 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, 2013), provided by Center for Wave Phenomena at Colorado School of Mines, to implement this approach (Aragao and Sava, 2020). 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.

3. Results

Synthetic seismic data are forward modelled (Pearce, 2022) from the ‘True’ velocity profile. These data are used to initiate the HWI and FWI, and their output velocity models are compared to the ‘True’ model (Fig. 1). Given its limiting assumptions, the model derived from HWI (blue) cannot resolve the base of the ice slab from the seismic data. Instead, it reaches its maximum velocity (~3750 m s⁻¹) at 40 m depth. This corresponds to the top of the ice slab, but this velocity is incorrectly propagated to the base of the model (Fig. 1a). By contrast, the FWI velocity model (red) detects both contacts of the ice slab: the resolution of its upper contact is lower than HWI, but the increase and deeper decrease of velocity correctly indicates the ice slab's presence. By 80 m depth, the FWI velocity model is comparable to that of the ‘True’ model, again indicating that the ice slab extends from depths 40–80 m. The relative performance of HWI and FWI is benchmarked by defining a percentage error from the ‘True’ velocity model (Fig. 1b), with FWI providing the more accurate representation throughout. Our current models assume a maximum wavelet frequency of 60 Hz; extending the bandwidth at the high frequency end would facilitate improved resolution but at the potential cost of inversion instabilities such as cycle skipping (Hu and others, 2018).

Fig. 1. (a) Velocity model outputs from HWI (Red) and FWI (Blue) compared to the True model (Black). The HWI generates a velocity profile that stops increasing in velocity at the depth of the ice slab top (40 m). This velocity is then extrapolated to the base of the model. FWI produces a velocity increase at 30 m depth and a decrease at 70 m depth, indicating the ice slab's location. (b) The absolute percentage error between the inversion models for HWI (Red) and (FWI) relative to the True model. The average difference for the total depth shows FWI reduces the error by 2.5%. A maximum difference of 20% for the HWI velocity model is observed in the top 40 m.

Figure 2 shows examples of the seismic traces generated by HWI (red) and FWI (blue), compared to the True model (black). Refractions in near-offset traces (Fig. 2a) are well-represented by both methods. However, at 500 m offset (Fig. 2b), the FWI trace is much closer to that of the reference model, while the over-estimate of velocity in the HWI case predicts a first-arrival that precedes that in the reference trace. Furthermore, FWI is able to model reflected arrivals, and Figure 2c shows an enlarged section of the near-offset traces, showing the insensitivity of HWI to a reflected arrival.

Fig. 2. Seismic data forward modelled from the True (Black), HWI (Red) and FWI (Blue). (a) All data with LMO applied with a velocity of 3800 m s⁻¹. (b) Comparison of trace from an offset of 5 m, the first arrival is modelled well by both the HWI and FWI, the green panel indicates the zoomed in section of figure (d). (c) Comparison of trace from an offset of 500 m. FWI closely reproduces the true data, while HWI poorly represents the true data. (d) Zoom of the reflection at 5 m offset produced by the ice slab visible at near offsets.

4. Future research priorities

The presence of ice slabs influences drainage and meltwater runoff across glaciers and ice masses (MacFerrin and others, 2019). 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, 2014). 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, 2004). Hubbard and others (2016) 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, 2019) 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, 2019) 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, 2018). 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 (2010), King and Jarvis (1991)).

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 (2014), Borisov and others (2020), Irnaka and others (2022)). 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 (2022), Kan and others (2023)).