A firn dielectric log depth-tied to an ice core on the West Antarctica Ice Sheet

: We have estimated a 1-D permittivity model from a 100m long variable offset GPR in the West Antarctic ice sheet. That model inherits the inaccuracies in depth from the velocity model, which should be corrected before attempting to correlate it with the densitylogfroma close-byborehole.Weperformedthatcorrectionbyaligningasynthetic ice density derived from a Maxwell Garnett two-phase mixture model to the ice core density measurements through dynamic time warping. The shifts to bring the permittivity estimates to their proper depths suggest a direct correlation of radar-derived data to borehole depths may suffer from noise to an unknown degree. The present methodology is within reach of a standard GPR survey, having at least one variable offset gather.


INTRODUCTION
Estimating the dielectric properties of ice is essential for glaciological research employing both ground and satellite radar on glaciers, ice sheets, and snow cover in polar regions. Those properties have been measured both in the field and in laboratories, with several studies indicating that permittivity is controlled by density, and this dependence follows the dielectric mixture theory (Tinga et al. 1973, Kovacs et al. 1995, Sugiyama et al. 2010. In this paper, we use GPR data collected at (84 • S, 79.5 • W) in the West Antarctic ice sheet to estimate the permittivity of the firn layer. The use of variable offset GPR for estimating ice permittivity begun a long time ago, (e.g. Jiracek & Bentley 1971, Jezek et al. 1978, Kovacs et al. 1995. The permittivity is accessed through the estimation of a velocity model, which bears significant uncertainties, accuracy in-depth being a particular challenge to achieve (Yilmaz 2001).
We corrected the depth inaccuracies of our permittivity estimates by tying them to the ice cores density measurements using Dynamic Time Warping (DTW) to align the former with the latter. The DTW aligns the permittivity estimates (sample) onto the density, a reference by piecewise linear stretching and compression of the warping segments with variable lengths (Nielsen et al. 1998, Pravdova et al. 2002, Tomasi et al. 2004. Notwithstanding DTW has begun associated with speech recognition (Rabiner et al. 1978, Sakoe & Chiba 1978, Gilbert et al. 2010 it has proved to be useful in several other applications. DTW has been used in seismic applications to automate well-to-well log correlation ( of firn is in the core of the methodology herein. For that matter the DTW was a tool to seek similarities between the two series with minute human intervention. The method is at the reach of a standard GPR survey.

THE DATA SETS The Density data
The density data was measured in ice cores from a 100.5m borehole dug at (84 • S, 79.5 • W at 1285m asl) during the austral summer 2011/12 in west Antarctic ice sheet. The density record comes from cores limited to the first 97.1m, with a fixed diameter of 0.85m and lengths up to .8m, weighted in situ with an nominal accuracy of 0.1g.
The densification of snow is a diagenetic process where pore space is reduced, and grains are transformed into a solid ice mass; the grains reaming at constant density, implying the porosity of polar firn decreases with depth. Accordingly the estimated density displays a standard, monotonically increasing data series with depth, (e.g. Cuffey & Paterson 2010), which can be expressed as a deviation from the density of pure glacier ice ρ ice (Herron & Langway Jr 1980), A solution to this equation iŝ where the C i , i = 0, 1, 2, are constants. The empirical model (1) ranges from the surface to the zone of pore close-off, assuming two (Herron & Langway Jr 1980) to three (Ligtenberg et al. 2011) distinct stages of firn densification: ρ 0.55, 0.55 -0.83, 0.83 g /cm 3 . In the last stage pores close off and the firn is transformed into ice with ρ ice = 0.917 g /cm 3 (Cuffey & Paterson 2010). Our density estimates display scatter at near surface < 10m and at depths > 50m so we rather used one single stage for the entire density log, with C 0 = 0.23 ± 0.17, C 1 = 0, C 2 = 0.32 ± 0.09.
Lower the total number of density estimates by decimating them by two, followed by applying an anti-alias low-pass filter to smooth the density oscillation. Figure 1 shows both the original ρ (z) and its decimated counterpart. From now on, we are going to use only the decimated density series.

The GPR data
We collected the data with a PULSE EKKO Pro GPR with two 100MHz unshielded antennae mounted in broadside perpendicular configuration, each of which on one sled. Positioning was done with differential post-processed GPS positioning with an accuracy of ≈ 0.05m. We concentrate here on a 100m common-source gather 107m long with 1111 traces, yielding around step size of 0.1m, we are going to refer to as G. Its time window is 2000ns, sampled with 0.8ns, which is enough to span over the entire snow/firn stratigraphy probed with the ice cores. Further details on the acquisition parameters can be found elsewhere (Martins et al. 2013).
Now estimate a velocity model for the G by automatic velocity picking on its semblance spectrum (Yilmaz 2001, Fomel 2009). Use that model to obtain a stacking velocity model to produce the stacked trace for the G shown in Figure 2. As we deal with horizontal layers, approximate the Root-Mean-Square (RMS) velocity by the obtained stacking velocity. Estimate the interval velocity for the subsurface layers through Dix formula (equation 3) to construct an interval 1-D velocity model (Yilmaz 2001), v (t) ≡ v n,n-1 = v 2 r,n t n -v 2 r,n-1 t n-1 t n -t n-1  where v r,n and v r,n-1 are the RMS velocities and t n and t n-1 are the zero-offset traveltimes for the n and n -1 reflectors, which define the layer n. It is known that interval estimates in the velocity model of equation (3) have inaccuracies due to a series of reasons, ranging from using stacking for the RMS velocities to eventual effects from a given gather spread and subsurface complexities. In particular the inaccuracies in v (t) are a function of: (i) error in the v r,n estimates, (ii) are inversely proportional to the n-layer thickness, the numerator in equation (3) as well as of the ratio ( v r,n-1/vn), and to errors in t n (Hajnal & Sereda 1981).
Velocities in model v (t) are exact for a horizontally layered medium, where a multi-layer conversion between time-and depth-domain, v (t) → v (z), is theoretically straightforward. Nevertheless, v (t) can over or underestimate the actual velocity; the uncertainty in the calculated interval velocity increases with depth and is inversely proportional to layer thickness, even when the input velocity and time variables are errors remain constant. Figure 3 shows too low or too high interval velocity estimates, evident at earlier and later times as well as in the range [600, 1000] ns. We choose to be rather conservative in discarding any non-physical velocity estimates so as not to interfere too much in our data as well as to display their variable reliability with depth; we have retained all estimates falling within the interval 0.05 < v (t) < 0.2 m /ns.

The firn effective medium
The effective permittivity is a complex parameter, ε e = ε e + ıε e , its imaginary part ε e being related to the losses due to electrical conduction, dielectric relaxation, dielectric resonance and other processes such as scattering (Evans 1965, Sihvola 2000. ε e can be disregarded when the loss tangent, tan δ = ε e/ε e 1, δ being related to the phase angle between the total current density and the incident electric field. For dry snow and firn ε e ε is dependent mainly on the volume fraction of ice and air, usually expressed as a mixing formula of its constituent permittivities and their fractional volumes, e.g., (Sihvola 2000). Therefore we regard a dry firn layer here as an idealized mixture model of ice and air with an effective permittivity ε e , which reflects its bulk properties, see equation (5) below. The relationship between velocity and the dielectric permittivity is given by where c is the speed of light in vacuum, the ε e being independent of frequency on the GPR plateau 1 -3GHz. For the sake of simplicity we drop the subscript of ε e . A basic dielectric mixing formula estimates an effective permittivity of the mixture, using its constituent permittivities and their fractional volumes. In this case, the mixture's electromagnetic response is as if the discrete medium were homogeneous, a frequency-dependent simplification valid only if inclusions are small compared to the wavelength. Following the dielectric mixture theory, we adopt a simple two-phase mixture model consisting of a background host medium with embedded spherical inclusions, both dielectrically isotropic. The Maxwell Garnett formula for the effective permittivity (Tinga et al. 1973, Sihvola & Kong 1988, Sihvola 2000 reads where the subscript ( h ) refers to the host material and ( i ) to the material in the spherical inclusions and f = V i /V h is the fractional volume of the inclusions in relation to the host volume. Note that the fractional volume is related to the total porosity The underlying assumption in equation (5) is that the mixture's effective permittivity is a function of the materials and the inclusions geometries. For GPR frequencies the relative permittivity of ice we adopt ε ice = 3.17, with an attenuation of ≈ 5 -0.5 dB /100m, in the temperature range [0 -40] • C (Evans 1965); for air ε air = 1. From the interval velocities v (t), equation (3) we estimate a permittivity ε (t) through relation (4), which which can be fed to the left hand side of equation (5). In this way equation (5) becomes a function of the fractional volume f only, which can be estimated for entire time, or depth, range spanned by gather G. Assuming a total porosity range φ 0 ∈ [0.1, 0.2] (Fourteau et al. 2020), we derive a first guess for f and obtain an estimate for it using the Powell hybrid method. We then use the estimates of f for two distinct models: M 1 an ice matrix with air inclusions and M 2 an air matrix with ice inclusions, which represent two limiting approaches based on the bonding of ice particles (Sugiyama et al. 2010). It is reasonable to assume the real fractional volume with depth, or time, should be contained between the fractional volumes of those two limiting models.
The estimates for f and by consequence for φ, allow for mapping the velocity v onto a velocity-derived "density" ρ ε = φρ air + (1φ) ρ ice . Such density ρ ε (t; z) may be then compared with the measured ρ (z) from the ice cores, an entirely independent dataset. As models M 1 and M 2 represent limiting ideal cases, there is no reason to favour one over another, so we adopt herein their average value for f . The above mentioned noise in the velocity model v (t), refer to Figure 3, propagates through the estimation of the models' volumetric fraction, equation (5), resulting in some inaccuracies in f . Again we were conservative in discarding only non-physical estimates, the ones who fell outside the interval 0 < f (z) < 1, totaling 13% of the estimates.

Tying the velocity-derived density to well density
In this section we tie the velocity-derived ρ v (t) to the well density ρ w (z) in order to have a time-depth mapping transforming ρ v (t) → ρ v (z). Once this is achieved we can correct v (z), equation (3), to their corrected depths and, by direct consequence of equation (4) to ε (z). We resort to a DTW algorithm to align the two series ρ v (t) and ρ w (z) through a sequence of shifts on the independent variable, which are a globally optimal solution to a non-linear optimization problem with linear inequality constraints, a problem with many local minima. The procedure we adopt here is akin to using a sonic tool to tie a well-log synthetic seismogram to a seismic data set twice. Well ties are used to produce a time-depth function enabling the correlation of properties measured in well-logs with the reflectors in seismic images (White & Simm 2003).
Before applying the DTW, it is necessary to go through a few processing steps to turn the data amenable to the algorithm. Start by computing first guess depths by using the interval velocities, equation (3), and the time at the layer boundaries. With those first guess depths map the velocity model in time into a model in depth, v (t) → v (z). Now transform velocity to permittivity and finally to density, ρ v (z). The decimated density record ρ w (z) with its somewhat ironed out oscillations should pair better with the significantly smoother ρ v (z).
The simple model of equation (5), a purely scattering model, is far from being able to deal with the complexities of compaction and densification of snow with depth. To make ρ v (z) comparable with ρ w (z) extract the firn densification from the density series subtracting the density estimates from the empirical model of equation (1). This quantityρ w should express the density anomaly to the normal firn densification based on Sorge's law which states the change in air space of the firn is linearly related to the weight of the snow above (Herron & Langway Jr 1980). The two data series, ρ v (z) and ρ w (z), are independent and so are their resolution and noise contents. The resolution of ρ w (z) is ≈ .8m, varying little with depth, whereas the resolution of ρ v (z) varies with time and noise contents, being related to 1 /4 of the wavelength or ≈ 0.4m. The noise contents of both series is difficult to estimate, but it is reasonable to assume it should be comparatively higher for ρ v (z) than for ρ w (z). This is a positive factor as we use the latter as the reference for pairing the two series with the DTW.
Begin by making both data series, ρ v (z) and ρ w (z), comparable though a standardizing procedure, where ρ i ≡ ρ v (z) , ρ w (z), ρ is an average, ς (ρ i ) is a standard deviation; the indexes i = 1, . . . , N become now the independent variable. The two standardized series ρ i have the same number of data points as the original series, are zero-mean and have a unit standard deviation. The standardization process also minimizes eventual y-axis discrepancies between the two series, minimizing the possibility of misalignment by the DTW algorithm. The mapping (7)   allowing to go back to the original values whenever needed.
The warping of seriesρ ε , the sample, ontô ρ w , the reference, is achieved by constructing a warp path between them W = (p 1 , p 2 , . . . , p K ), where each path element p k is linked to the two series indexes i, i , for the N elements inρ ε and ρ w , respectively. The warp path length W is bound to N ≤ K ≤ 2N -1 and subject to several criteria to assure continuity and monotonicity. First the warp path start and end at the first and last elements of the two sequences, p 1 = (1, 1) and p K = (N, N), all elements considered. Second, the ordering of the two aligned sequences has to be preserved, p k i, i → p k+1 i, i ⇒ i ≤ i ≤ (i + 1) and i ≤ i ≤ i + 1 . Third the elements of W are monotonically spaced in the independent variable, thus preventing big jumps, The process of warping seeks a minimum-distance, where d p k , p k+1 is the distance between two contiguous elements. D p attains its minimum when the sample is correctly warped onto the reference signal (Sakoe & Chiba 1978). We perform the DTW with an algorithm that uses a correlation optimized warping, or COW, which aligns the sample onto the reference by piecewise linear stretching and compression of warping segments of variable lengths l (Nielsen et al. 1998, Pravdova et al. 2002, Tomasi et al. 2004). An integer slack parameter limits the range of possible segments l, initially set to unity, s = 1. We apply the COW algorithm on all N/l segments through dynamic programming (Nielsen et al. 1998, Pravdova et al. 2002, Tomasi et al. 2004. A complete description of the DTW and COW algorithms is well beyond the scope of this work; the reader is kindly referred to the literature cited herein.
Apply the DTW to the standardized versions of both velocity-derived and well density,ρ v andρ w , the reference. The algorithm applies shifts in the independent variable, stretching and squeezingρ v to matchρ w until it finds an optimal solution. The result is aρ v with decreased depth inaccuracies due to the ice core data, allowing for a more accurate relation of permittivity to density. Figure 4 shows theρ v before and after being tied toρ w , showing the optimal alignment of the two series, as well depicting the shifts as vectors. The algorithm did an excellent job in aligning the two series, more conspicuously in the interval 40 < z < 85m. It is essential to note the quality of the alignment notwithstanding ρ v being significantly smoother thanρ w , a fact that harmed the shallower data, z < 40m. The algorithm performed poorer toward the end of the borehole, z > 85m, whereρ v andρ w display a vertical shift, probably resulting from inaccuracies on the former. These problems do not diminish the results' quality; it would be too optimistic to expect good results throughout. At least within the depth range 40 < z < 85m the errors in depth estimates ofρ v should now be more uniformly distributed. The magnitudes of the shift vectors shown in Figure 4 suggest the original errors in depth estimates were not only highly variable but also significantly larger than the depth resolution of the individual ice cores, ∼ = .8m.

RESULTS
The first guess depths derived from the interval velocity model v (t) are now corrected using the DTW shifts from the warped results, thus producing permittivity estimates at corrected depths, ε e (z). From now on, the permittivity is depth-corrected or somewhat tied to the well depths. Figure 5 shows ε e (z) for almost the entire depth interval of the well, but for the shallowest depths. Nevertheless the ε e (z) estimates span over the three distinct stages of firn densification, corresponding to densities ρ 0.55, 0.55 -0.83, 0.83 g /cm 3 . The permittivity estimates are relatively high throughout, but we are reasonably confident on them as we do have an independent velocity model obtained along with a 500m profile on the same site, where it is first 100m are coincident with gathering G. That longer profile led to an independent 2-D velocity model, used in producing a pre-stack time-migrated profile (Yilmaz 2001), which revealed a velocity variation v 2D ∈ [0.16, 0.24] m /ns (Martins & Travassos 2015), which corresponds to a permittivity variation of ε ∈ [3.5, 1.6]. This is a marginally shorter interval than the one reported herein for the 1-D model. Therefore most of our estimates ε v (z) should be correct but for a few peaks and valleys resulting from noise in the data.

CONCLUSIONS
The present paper dealt with two independent data sets measured on the firn layer of the ice cover at a site located in the West Antarctic ice sheet. The two datasets were acquired during the same fieldwork, comprising density measurements of ice cores and firn permittivity estimates from a 1-D velocity model from a multi-offset GPR profile obtained near the borehole.
While density estimates inaccuracies are more or less confined within a relatively narrow range due to errors in estimating cores volumes and weights, with a depth resolution of 0.85m, the radar velocity v (t), or for that matter the directly derived permittivity ε (t), are notoriously prone to inaccuracies arising from the simplifying assumptions in modeling an unknown subsurface complexity. Those inaccuracies become more intense and variable with the transformation of time to depth, ε (t) → ε (z), making a direct correlation with density ρ (z) difficult.
Nevertheless, density and permittivity are functions of the air and ice volumetric contents of firn, suggesting they should bear some direct or inverse similarity in their respective morphology. Moreover, density is a direct measure, whereas the permittivity is an estimate from a remote measurement at the surface, with a markedly distinct sensibility to localized changes in the volumetric contents of ice and air. All this considered renders a correlation between the two series a problematic task. A manual attempt of doing that correlation leaves too much room for disposition to be reliable.
We tackled correlating the two series based on their morphology using a DTW algorithm, which warps the permittivity by having the density as a reference. This is done entirely by mathematical optimization, with no human intervention but for a few trial-and-error choices on a pair of parameters. A necessary by-product of warping is a direct mapping of permittivity to a set of corrected depths.
The DTW algorithm did a good job in aligning the two series, notwithstanding their widely different noise contents and depth resolution, which is more conspicuously for intermediate depths, 40 < z < 85m. The results shown in this work showed the necessary shifts to permittivity depths are not only highly variable but also significantly larger than the depth resolution of density series, ∼ = .8m. The worst results were toward the end of the borehole, z > 85m.
We bring this section to a close, saying the anything but apparent alignment obtained herein, Figure 4, is the main contribution of this paper. It is definitively improbable one may have reached comparable results based on personal judgment only. Moreover, the needed shifts to bring the permittivity estimates to their proper depths suggest other results obtained by directly correlating radar-derived data to borehole depths may suffer from noise to an unknown degree. As a last remark, the methodology presented here is at the reach of a standard GPR survey, typically performed with a few variable offset and more numerous fixed-offset gathers.