Acessibilidade / Reportar erro

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

Abstract

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 density log from a close-by borehole. We performed that correction by aligning a synthetic 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.

Key words
ice core; GPR; velocity model; dielectric property; well logging

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. 1973TINGA WR, VOSS W & BLOSSEY D. 1973. Generalized approach to multiphase dielectric mixture theory. J Appl Phys 44(9): 3897-3902., Kovacs et al. 1995KOVACS A, GOW AJ & MOREY RM. 1995. The in-situ dielectric constant of polar firn revisited. Cold Reg Sci Technol 23(3): 245-256., Sugiyama et al. 2010SUGIYAMA S, ENOMOTO H, FUJITA S, FUKUI K, NAKAZAWA F & HOLMLUND P. 2010. Dielectric permittivity of snow measured along the route traversed in the Japanese–Swedish Antarctic Expedition 2007/08. Ann Glaciol 51(55): 9-15.).

In this paper, we use GPR data collected at (84S, 79.5W) 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 1971JIRACEK G & BENTLEY CR. 1971. Velocity of electromagnetic waves in Antarctic ice. Antarctic snow and ice studies II 16: 199-208., Jezek et al. 1978JEZEK KC, CLOUGH JW, BENTLEY CR & SHABTAIE S. 1978. Dielectric permittivity of glacier ice measured in situ by radar wide-angle reflection. J Glaciol 21(85): 315-329., Kovacs et al. 1995KOVACS A, GOW AJ & MOREY RM. 1995. The in-situ dielectric constant of polar firn revisited. Cold Reg Sci Technol 23(3): 245-256.). 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 2001YILMAZ O. 2001. Seismic Data Analysis: Processing, Inversion, and Interpretation of Seismic Data. 8801 South Yale Suite 500 Tulsa OK 74137 USA: SEG.).

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. 1998NIELSEN NPV, CARSTENSEN JM & SMEDSGAARD J. 1998. Aligning of single and multiple wavelength chromatographic profiles for chemometric data analysis using correlation optimised warping. J Chromatogr A 805(1): 17-35., Pravdova et al. 2002PRAVDOVA V, WALCZAK B & MASSART D. 2002. A comparison of two algorithms for warping of analytical signals. Anal Chim Acta 456(1): 77-92., Tomasi et al. 2004TOMASI G, VAN DEN BERG F & ANDERSSON C. 2004. Correlation optimized warping and dynamic time warping as preprocessing methods for chromatographic data. J Chemom 18(5): 231-241.). Notwithstanding DTW has begun associated with speech recognition (Gilbert et al. 2010GILBERT JM, RYBCHENKO SI, HOFE R, ELL SR, FAGAN MJ, MOORE RK & GREEN P. 2010. Isolated word recognition of silent speech using magnetic implants and sensors. Med Eng Phys 32(10): 1189-1197., Rabiner et al. 1978RABINER L, ROSENBERG A & LEVINSON S. 1978. Considerations in dynamic time warping algorithms for discrete word recognition. IEEE Trans Signal Process 26(6): 575-582., Sakoe & Chiba 1978SAKOE H & CHIBA S. 1978. Dynamic programming algorithm optimization for spoken word recognition. IEEE Trans Signal Process 26(1): 43-49.) it has proved to be useful in several other applications. DTW has been used in seismic applications to automate well-to-well log correlation (White & Simm 2003WHITE R & SIMM R. 2003. Tutorial: Good practice in well ties. First Break 21(10)., Munoz & Hale 2012MUNOZ A & HALE D. 2012. Automatically tying well logs to seismic data. CWP-725 URL http://cwp.mines.edu/ Documents/cwpreports/CWP-764.pdf.
http://cwp.mines.edu/ Documents/cwprepor...
).

The fact of density and permittivity being functions of the air and ice volumetric contents 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 (84S, 79.5W 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 2010CUFFEY K & PATERSON W. 2010. Physics of Glaciers, 4th Edn.), which can be expressed as a deviation from the density of pure glacier ice ρice (Herron & Langway Jr 1980HERRON MM & LANGWAY JR CC. 1980. Firn densification: an empirical model. J Glaciol 25: 373-385.),

dρ(z)dz=C0(ρ(z)ρice).(1)
A solution to this equation is
ρ̂(z)=C1+C2exp(C0z),(2)
where the Ci,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 1980HERRON MM & LANGWAY JR CC. 1980. Firn densification: an empirical model. J Glaciol 25: 373-385.) to three (Ligtenberg et al. 2011LIGTENBERG S, HEILSEN M & VAN DE BROEKE M. 2011. An improved semi-empirical model for the densification of Antarctic firn. Cryosphere 5(4): 809-819.) distinct stages of firn densification: ρ0.55,0.550.83,0.83gcm3. In the last stage pores close off and the firn is transformed into ice with ρice=0.917gcm3 (Cuffey & Paterson 2010CUFFEY K & PATERSON W. 2010. Physics of Glaciers, 4th Edn.). 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 C0=0.23±0.17,C1=0,C2=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.

Figure 1
Panel (a) shows the density ρ(z) as dots and the model expressed by equation (2), as a solid curve. Panel (b) shows the decimated density series as a solid line.

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 𝒢. 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. 2013MARTINS SS, TRAVASSOS JM & SIMOES JC. 2013. Interpolating aliased end-on GPR data. In: Explor Geophys, p. 1868-1872.).

Now estimate a velocity model for the 𝒢 by automatic velocity picking on its semblance spectrum (Fomel 2009FOMEL S. 2009. Velocity analysis using AB semblance. Geophys Prospect 57(3): 311-321., Yilmaz 2001YILMAZ O. 2001. Seismic Data Analysis: Processing, Inversion, and Interpretation of Seismic Data. 8801 South Yale Suite 500 Tulsa OK 74137 USA: SEG.). Use that model to obtain a stacking velocity model to produce the stacked trace for the 𝒢 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 2001YILMAZ O. 2001. Seismic Data Analysis: Processing, Inversion, and Interpretation of Seismic Data. 8801 South Yale Suite 500 Tulsa OK 74137 USA: SEG.),

v(t)vn,n1=(vr,n2tnvr,n12tn1tntn1)1/2(3)
where vr,n and vr,n1 are the RMS velocities and tn and tn1 are the zero-offset traveltimes for the n and n1 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 vr,n estimates, (ii) are inversely proportional to the n-layer thickness, the numerator in equation (3) as well as of the ratio (vr,n1vn), and to errors in tn (Hajnal & Sereda 1981HAJNAL Z & SEREDA I. 1981. Maximum uncertainty of interval velocity estimates. Geophysics 46(11): 1543-1547.).

Figure 2
Panel (b) shows an automatic velocity picking on the semblance spectrum from gather 𝒢 shown in Panel (a). Panel (c) shows the stacked trace for gather 𝒢.

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.2mns.

Figure 3
The 1-D interval velocity model v(t) estimated from gather 𝒢.

Methods

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 1965EVANS S. 1965. Dielectric properties of ice and snow–a review. J Glaciol 5(42): 773-792., Sihvola 2000SIHVOLA A. 2000. Mixing rules with complex dielectric coefficients. Subsurface sensing technologies and applications 1(4): 393-415.). εe can be disregarded when the loss tangent, tanδ =ε eε e1, δ 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 2000SIHVOLA A. 2000. Mixing rules with complex dielectric coefficients. Subsurface sensing technologies and applications 1(4): 393-415.). 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

v=cεe,(4)
where c is the speed of light in vacuum, the εe being independent of frequency on the GPR plateau 13GHz. 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. 1973TINGA WR, VOSS W & BLOSSEY D. 1973. Generalized approach to multiphase dielectric mixture theory. J Appl Phys 44(9): 3897-3902., Sihvola & Kong 1988SIHVOLA AH & KONG JA. 1988. Effective permittivity of dielectric mixtures. IEEE Trans Geosci Remote Sens 26(4): 420-429., Sihvola 2000SIHVOLA A. 2000. Mixing rules with complex dielectric coefficients. Subsurface sensing technologies and applications 1(4): 393-415.) reads

εe=3fεh(εiεh)(εi+2εh)f(εiεh)+εh,(5)
where the subscript (h) refers to the host material and (i) to the material in the spherical inclusions and f=ViVh is the fractional volume of the inclusions in relation to the host volume. Note that the fractional volume is related to the total porosity
ϕ=f1+f.(6)
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 50.5dB100m, in the temperature range [040]C (Evans 1965EVANS S. 1965. Dielectric properties of ice and snow–a review. J Glaciol 5(42): 773-792.); 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 𝒢. Assuming a total porosity range ϕ0[0.1,0.2] (Fourteau et al. 2020FOURTEAU K, ARNAUD L, FAïN X, MARTINERIE P, ETHERIDGE DM, LIPENKOV V & BARNOLA JM. 2020. Historical porosity data in polar firn. Earth Syst Sci Data 12(2): 1171-1177.), 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: 1 an ice matrix with air inclusions and 2 an air matrix with ice inclusions, which represent two limiting approaches based on the bonding of ice particles (Sugiyama et al. 2010SUGIYAMA S, ENOMOTO H, FUJITA S, FUKUI K, NAKAZAWA F & HOLMLUND P. 2010. Dielectric permittivity of snow measured along the route traversed in the Japanese–Swedish Antarctic Expedition 2007/08. Ann Glaciol 51(55): 9-15.). 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 1 and 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 2003WHITE R & SIMM R. 2003. Tutorial: Good practice in well ties. First Break 21(10).).

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 1980HERRON MM & LANGWAY JR CC. 1980. Firn densification: an empirical model. J Glaciol 25: 373-385.). 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 14 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,

ρ iρ i^ =1σ (ρ i)(ρ iρ ¯ ),(7)
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 ρî 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) is invertible, allowing to go back to the original values whenever needed.

The warping of series ρ̂ε, the sample, onto ρ̂w, the reference, is achieved by constructing a warp path between them W=(p1,p2,,pK), where each path element pk 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 NK2N1 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, p1=(1,1) and pK=(N,N), all elements considered. Second, the ordering of the two aligned sequences has to be preserved, pk(i,i)pk+1(î,î)iî(i+1)and iî(i+1). Third the elements of W are monotonically spaced in the independent variable, thus preventing big jumps, pk(i,i)pk+1(î,î)(iî)0and(iî)0.

The process of warping seeks a minimum-distance,

Dp=12Nmin{k=1Kd(pk,pk+1)},(8)
where d(pk,pk+1) is the distance between two contiguous elements. Dp attains its minimum when the sample is correctly warped onto the reference signal (Sakoe & Chiba 1978SAKOE H & CHIBA S. 1978. Dynamic programming algorithm optimization for spoken word recognition. IEEE Trans Signal Process 26(1): 43-49.). 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. 1998NIELSEN NPV, CARSTENSEN JM & SMEDSGAARD J. 1998. Aligning of single and multiple wavelength chromatographic profiles for chemometric data analysis using correlation optimised warping. J Chromatogr A 805(1): 17-35., Pravdova et al. 2002PRAVDOVA V, WALCZAK B & MASSART D. 2002. A comparison of two algorithms for warping of analytical signals. Anal Chim Acta 456(1): 77-92., Tomasi et al. 2004TOMASI G, VAN DEN BERG F & ANDERSSON C. 2004. Correlation optimized warping and dynamic time warping as preprocessing methods for chromatographic data. J Chemom 18(5): 231-241.). An integer slack parameter limits the range of possible segments l, initially set to unity, 𝔰=1. We apply the COW algorithm on all N/l segments through dynamic programming (Nielsen et al. 1998NIELSEN NPV, CARSTENSEN JM & SMEDSGAARD J. 1998. Aligning of single and multiple wavelength chromatographic profiles for chemometric data analysis using correlation optimised warping. J Chromatogr A 805(1): 17-35., Pravdova et al. 2002PRAVDOVA V, WALCZAK B & MASSART D. 2002. A comparison of two algorithms for warping of analytical signals. Anal Chim Acta 456(1): 77-92., Tomasi et al. 2004TOMASI G, VAN DEN BERG F & ANDERSSON C. 2004. Correlation optimized warping and dynamic time warping as preprocessing methods for chromatographic data. J Chemom 18(5): 231-241.). 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.

Figure 4
Standartized series ρ̃v, in red, and ρ̃w, in black. Panel (a) shows the original series and the relative shifts which should be applied to ρ̃v in order to align it to ρ̃w, the vector lengths proportional to the amount of shift. Panel (b) shows the aligned series.

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.550.83,0.83gcm3. 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 𝒢. That longer profile led to an independent 2-D velocity model, used in producing a pre-stack time-migrated profile (Yilmaz 2001YILMAZ O. 2001. Seismic Data Analysis: Processing, Inversion, and Interpretation of Seismic Data. 8801 South Yale Suite 500 Tulsa OK 74137 USA: SEG.), which revealed a velocity variation v2D[0.16,0.24]mns (Martins & Travassos 2015MARTINS SS & TRAVASSOS JM. 2015. Interpolating wide-aperture ground-penetrating radar beyond aliasing. Geophysics 80(2): H13-H22.), 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.

Figure 5
Permittivity estimates εe(z) at their corrected depths.

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

ACKNOWLEDGMENTS

This work was fully supported by the Programa Antártico Brasileiro (PROANTAR), through the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) and Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES). The two main sources come from the INCT da Criosfera (CAPES Proj. 88887.136384/2017-00) and PROANTAR/CNPq (Proj. 442755/2018-0).

REFERENCES

  • CUFFEY K & PATERSON W. 2010. Physics of Glaciers, 4th Edn.
  • EVANS S. 1965. Dielectric properties of ice and snow–a review. J Glaciol 5(42): 773-792.
  • FOMEL S. 2009. Velocity analysis using AB semblance. Geophys Prospect 57(3): 311-321.
  • FOURTEAU K, ARNAUD L, FAïN X, MARTINERIE P, ETHERIDGE DM, LIPENKOV V & BARNOLA JM. 2020. Historical porosity data in polar firn. Earth Syst Sci Data 12(2): 1171-1177.
  • GILBERT JM, RYBCHENKO SI, HOFE R, ELL SR, FAGAN MJ, MOORE RK & GREEN P. 2010. Isolated word recognition of silent speech using magnetic implants and sensors. Med Eng Phys 32(10): 1189-1197.
  • HAJNAL Z & SEREDA I. 1981. Maximum uncertainty of interval velocity estimates. Geophysics 46(11): 1543-1547.
  • HERRON MM & LANGWAY JR CC. 1980. Firn densification: an empirical model. J Glaciol 25: 373-385.
  • JEZEK KC, CLOUGH JW, BENTLEY CR & SHABTAIE S. 1978. Dielectric permittivity of glacier ice measured in situ by radar wide-angle reflection. J Glaciol 21(85): 315-329.
  • JIRACEK G & BENTLEY CR. 1971. Velocity of electromagnetic waves in Antarctic ice. Antarctic snow and ice studies II 16: 199-208.
  • KOVACS A, GOW AJ & MOREY RM. 1995. The in-situ dielectric constant of polar firn revisited. Cold Reg Sci Technol 23(3): 245-256.
  • LIGTENBERG S, HEILSEN M & VAN DE BROEKE M. 2011. An improved semi-empirical model for the densification of Antarctic firn. Cryosphere 5(4): 809-819.
  • MARTINS SS & TRAVASSOS JM. 2015. Interpolating wide-aperture ground-penetrating radar beyond aliasing. Geophysics 80(2): H13-H22.
  • MARTINS SS, TRAVASSOS JM & SIMOES JC. 2013. Interpolating aliased end-on GPR data. In: Explor Geophys, p. 1868-1872.
  • MUNOZ A & HALE D. 2012. Automatically tying well logs to seismic data. CWP-725 URL http://cwp.mines.edu/ Documents/cwpreports/CWP-764.pdf
    » http://cwp.mines.edu/ Documents/cwpreports/CWP-764.pdf
  • NIELSEN NPV, CARSTENSEN JM & SMEDSGAARD J. 1998. Aligning of single and multiple wavelength chromatographic profiles for chemometric data analysis using correlation optimised warping. J Chromatogr A 805(1): 17-35.
  • PRAVDOVA V, WALCZAK B & MASSART D. 2002. A comparison of two algorithms for warping of analytical signals. Anal Chim Acta 456(1): 77-92.
  • RABINER L, ROSENBERG A & LEVINSON S. 1978. Considerations in dynamic time warping algorithms for discrete word recognition. IEEE Trans Signal Process 26(6): 575-582.
  • SAKOE H & CHIBA S. 1978. Dynamic programming algorithm optimization for spoken word recognition. IEEE Trans Signal Process 26(1): 43-49.
  • SIHVOLA A. 2000. Mixing rules with complex dielectric coefficients. Subsurface sensing technologies and applications 1(4): 393-415.
  • SIHVOLA AH & KONG JA. 1988. Effective permittivity of dielectric mixtures. IEEE Trans Geosci Remote Sens 26(4): 420-429.
  • SUGIYAMA S, ENOMOTO H, FUJITA S, FUKUI K, NAKAZAWA F & HOLMLUND P. 2010. Dielectric permittivity of snow measured along the route traversed in the Japanese–Swedish Antarctic Expedition 2007/08. Ann Glaciol 51(55): 9-15.
  • TINGA WR, VOSS W & BLOSSEY D. 1973. Generalized approach to multiphase dielectric mixture theory. J Appl Phys 44(9): 3897-3902.
  • TOMASI G, VAN DEN BERG F & ANDERSSON C. 2004. Correlation optimized warping and dynamic time warping as preprocessing methods for chromatographic data. J Chemom 18(5): 231-241.
  • WHITE R & SIMM R. 2003. Tutorial: Good practice in well ties. First Break 21(10).
  • YILMAZ O. 2001. Seismic Data Analysis: Processing, Inversion, and Interpretation of Seismic Data. 8801 South Yale Suite 500 Tulsa OK 74137 USA: SEG.

Publication Dates

  • Publication in this collection
    30 May 2022
  • Date of issue
    2022

History

  • Received
    4 June 2021
  • Accepted
    8 Dec 2021
Academia Brasileira de Ciências Rua Anfilófio de Carvalho, 29, 3º andar, 20030-060 Rio de Janeiro RJ Brasil, Tel: +55 21 3907-8100 - Rio de Janeiro - RJ - Brazil
E-mail: aabc@abc.org.br