Reanalysis of Polythermal Glacier Thermal Structure Using Radar Diffraction Focusing

Ground‐penetrating radar (GPR) is widely used on polythermal glaciers to image bed topography and detect internal scatter due to water inclusions in temperate ice. The glaciological importance of this is twofold: bed topography is a primary component for modeling the long‐term evolution of glaciers and ice sheets, and the presence of temperate ice and associated englacial water significantly reduces overall ice viscosity. Englacial water has a direct influence on radar velocity, which can result in incorrect observations of bed topography due to errors in depth conversion. Assessment of radar velocities often requires multi‐offset surveys, yet these are logistically challenging and time consuming to acquire, hence techniques to extract velocity from common‐offset data are required. We calculate englacial radar velocity from common offset GPR data collected on Von Postbreen, a polythermal glacier in Svalbard. We first separate and enhance the diffracted wavefield by systematically assessing data coherence. We then use the focusing metric of negative entropy to deduce a migration velocity field and produce a velocity model which varies spatially across the glacier. We show that this velocity field successfully differentiates between areas of cold and temperate ice and can detect lateral variations in radar velocity close to the glacier bed. This velocity field results in consistently lower ice depths relative to those derived from a commonly assumed constant velocity, with an average difference of 4.9 ± 2.5% of local ice depth. This indicates that diffraction focusing and velocity estimation are crucial in retrieving correct bed topography in the presence of temperate ice.

Here, we present a method for the automated estimation of a 2D englacial radar velocity field to provide an improved understanding of the temperate-ice distribution within a polythermal glacier using a series of 2D, common offset GPR data sets. The underpinning method of inverting from diffractions in common-offset data has been developed in seismic applications Preine et al., 2020;Schwarz & Gajewski, 2017), yet has previously been applied only manually to radioglaciological data (Bradford, 2005). We apply the workflow described here to field data collected on a polythermal glacier in Svalbard, and detect areas of low radar-wave velocity coincident with high signal scattering within the data, in accordance with previous studies on similar polythermal glaciers . We then use the resultant radar-wave velocity field to estimate the improvements in glacier bed topography, compared to the assumption of a constant representative ice velocity (Navarro et al., 2014).

Physical Setting
Von Postbreen (78°25.07 N and 17°43.27 E) is a polythermal glacier (Sevestre et al., 2015) located in the Tempelfjorden region of Svalbard. It is approximately 15 km in length and has an area of 168 km 2 (König et al., 2014). The upper tributaries include Fimbulisen to the south, Potpeschniggbreen to the North-East, and Phillippbreen to the north. The lower ablation area consists of two tongues, the southernmost of which is the widest at 2 km width on average and forms the focus of this study. The glacier is classed as surge-type and is currently in a long quiescent phase, typical of many Svalbard glaciers (J. A. Dowdeswell et al., 1991), having last surged in 1870 (De Geer, 1910). Crevasse-squeezed ridges are present in the modern-day forefield and are noted as a sign of historic surging (Farnsworth et al., 2016). The glacier was historically confluent with Tunabreen with a continuous calving front in Tempelfjorden, but has steadily retreated throughout the twentieth century and is now entirely land-terminating (Flink et al., 2015;Sweeting & Groom, 1956). While Von Postbreen has been long quiescent, neighboring Tunabreen has surged regularly and is well-studied (How et al., 2019). The confluent minor glacier Bogebreen surged in 1980 (J. Dowdeswell et al., 1984).
Geophysical surveys on Von Postbreen were first undertaken in 1980 using an airborne radio echo sounding system with a central frequency of 60 MHz (Drewry et al., 1980), sounding both the bed and internal scatter within the ice as part of extensive surveys over the archipelago (J. Dowdeswell et al., 1984). Bamber (1987) applied a reflection coefficient-based analysis of the uppermost internal scattering horizon on similar glaciers across Svalbard from the same surveys, and concluded that it represents a layer of temperate ice with a water content of approximately 3%. In 2012 and 2015, Sevestre et al. (2015) undertook an extensive GPR campaign at 100 MHz, successfully mapping the distribution of temperate ice within the glacier, but encountering strong scattering throughout the deepest ice, precluding the ability to pick the bed confidently throughout.

Data Acquisition
Data were collected on the 13 and 15 March 2018, using a 25 MHz PulseEkko Pro GPR system mounted on plastic pulks towed behind a snowmobile. We acquired eight transverse profiles across the ablation zone, each approximately 2 km long and spaced by 200 m, and a 13 km longitudinal profile extending to the upper accumulation zone (Figure 1). Antennas were mounted in a parallel endfire configuration in line with the direction of travel; antenna dipole centers were separated by an offset of 5 m, which in this study we refer to as a common-offset acquisition. A 4 ns sample interval (sampling frequency 250 MHz) was used which gives a good reconstruction of the received waveform, and a trace length of 1,125 samples was used. The snowmobile was driven at 10-15 km/ hr, giving an average trace spacing in the raw data of 0.5 m. Traces were located with a Leica GS10 differential GPS rover feeding real time kinematic corrections directly to the GPR system, the base station for which was located on Przybyllokfjellet to the south of the glacier at an elevation of 720 m above sea level. The base station location was corrected using the Canadian Spatial Reference System Precise Point Positioning service following a 30 min static recording time on a previous day.

Methodology
The approach can be described in five main steps. We first pre-process the acquired radar data to remove significant noise sources. We then apply a diffraction enhancement algorithm to remove planar signals and noise in the data, and generate coherent diffractions. Then, following an approach previously developed and applied by Fomel et al. (2007) to assess the velocity from a post-stack seismic data set, we apply diffraction focusing to the radar profiles to generate time-domain englacial radar-wave velocity fields. We modify this approach by applying a novel focusing measure of the local signal entropy as a function of focusing velocity. We then undertake parameter testing on shaping regularization of the radar-wave-velocity field to assess the tradeoff between overfitting, with associated non-physical local velocities, and oversmoothing, whereby local anomalies are lost. Finally, we use the estimated velocity field to (a) derive a model of water content distribution across the glacier, and (b) migrate and depth-convert the radar data to derive a correct bed topography model.

Pre-Processing
Data pre-processing was required before application of the velocity analysis algorithm. Traces were re-sorted to 1 m intervals using associated GPS locations. A static time shift was applied such that t = 0 ns at the first break of the direct arrival, followed by a 2 MHz high-pass (dewow) filter to remove low frequency or static, DC noise.
Due to the small offset between transmitter and receiver, imposed by logistical constraints, significant ringing noise was experienced through approximately the first 300 ns of each trace. To remove this noise, we used a modification of the singular value decomposition (SVD) filter, developed by Freire and Ulrych (1988) for seismic data and applied to helicopter-borne radar data over Glacier de la Plaine Morte, Switzerland, by Grab et al. (2018). This approach decomposes the radargram D of m traces and n data points to where U is a m × m orthogonal matrix consisting of the eigenvectors u i of the covariance matrix of DD T , (the left singular vectors), S is a diagonal matrix whose diagonal elements are the singular values σ i of D in descending order, and V is a n × n matrix of the eigenvectors v i of D T D (the right singular vectors). Equation 1 shows that D can be reconstructed by a σ i -weighted summation of the orthogonal eigenimages . The ringing noise constitutes a highly correlated signal between traces which can be described well by a low number of high-amplitude eigenvalues σ i . Filtering is applied over S, scaling high-amplitude values of σ i and retaining low values of σ i which constitute the uncorrelated internal reflections. To generate a modified vector ̂ , and the filtered radargram is reconstructed by This approach is successful in removing near-surface ringing, yet in the case of a flat or low-dip bed reflection, it can result in an additional attenuation of the bed reflection amplitude. We therefore modify this approach by applying a smooth amplitude taper to the data between t = 0 ns and a user-defined two-way travel time (twtt), such that we generate two matrices D = D ns + D bed , where D ns captures the near-surface and upper scattering layer and D bed represents the deeper radargram including the bed, with data partitioned in twtt by a tapered transition. SVD filtering (Equations 1 and 2) is applied to D ns , and the output radargram is reconstructed by = + . Parameter testing is required on the scale of the pre-filter data partition (separating D ns and D bed ) which we generally define as a 100-sample plateau with a 200-sample cosine taper, and the filter of eigenvalues, which we scale from zero using a cosine taper over the first 20 values.

Retrieving the Diffracted Wavefield
Estimation of a full-ice-depth velocity field requires retrieval and focusing of diffractions throughout the ice column. The bed reflection is high-amplitude and specular in nature relative to diffractions within the ice. As such, there is little change with migration velocity perturbation, hence little variation in a time-velocity-domain common image gather. Suppression of the bed reflection is therefore essential to see the diffracted signals clearly, and thus enable retrieval of a full-depth average velocity and interval velocity of the deepest regions of ice.
We use a coherent summation and subtraction scheme (Schwarz, 2019) to extract diffraction-only data from the pre-processed radar lines. This is a three-step process: first, the slope of arriving wavefronts is estimated by a directional stacking optimization problem; second, planar wavefronts are subtracted from the raw data; and third, the diffractions are further enhanced by applying local time and amplitude corrections.
To remove planar wavefronts, the coherence of arriving wavefronts is estimated as a directional stacking optimization problem. For every data point D(x 0 , t 0 ), an estimation of the inclination of the most coherent arrival can be represented as an optimization problem using the semblance norm (Neidell & Taner, 1971), in which the horizontal slowness p acts as the optimization parameter where the index i refers to the ith trace within a spatial aperture, and δt is a time window over which summation is performed.
To ensure reliable convergence even in the case of strong interference and noise contamination, the problem expressed in Equation 3 is solved within a local aperture spanning n traces using a global differential evolution optimizer that is constrained to low wavenumbers (Storn & Price, 1997). Directional stacking on the original data is then applied using the resultant slope field, leaving strong, near-planar reflections which are subtracted from the raw data to leave high slope angle arrivals (Schwarz & Gajewski, 2017). Used alone, this step provides an approximately diffraction-only data set, yet much of the high-frequency scatter is retained in temperate ice, resulting in imperfect and noisy hyperbolas. In addition, the derived low-wavenumber wavefield may encode sufficiently accurate waveforms, but is likely to suffer from amplitude deviations and temporal offsets, which can compromise the subsequent step of noise attenuation.
To enhance the diffracted wavefield further and suppress high-frequency scatter, we follow Schwarz (2019) in applying a least squares optimization approach to subtract the previously obtained low-wavenumber estimate from the input data. The optimization for the local amplitude scaling α 0 and a local time shift τ 0 minimizes the misfit where D low is the coherent reflection model estimated in the previous step. To remove planar wavefronts and estimate the local slope, we used an aperture of 20 m, a coherence window of 20 ns and a maximum angle of incidence of 3°. Diffraction separation and coherent stacking were then performed with a wider aperture of 100 m, with a temporal adaptation window of 100 ns (δt in Equation 4). The result can be considered to be a diffraction-only data set with contaminating noise and planar reflections suppressed. We apply this approach to our data, and use this diffraction-only data set going forward for velocity analysis.

Time Migration Velocity Analysis
We follow the approach of Fomel (2003) to estimate time-migration velocity field using local focusing velocities for each hyperbola in the data set. This step, and subsequent processing, are undertaken in Madagascar, an open-source geophysical processing package (Madagascar, 2021). An initial Stolt F-K migration is applied using a velocity of 0.1 m/ns, which is below the minimum expected value for radio-wave transit through glacier ice. We then apply velocity continuation using a velocity interval of 0.005 m/ns up to 0.2 m/ns to produce a 3D data volume D(t, x, v). To constrain the focusing function, we chose a velocity discretization of 0.005 m/ns as a balance between velocity precision and computational demand; this follows the discretization used in the manual MVA approach of Bradford (2005).
We measure the level of local focusing through use of the statistical measure of negative entropy. First introduced by Shannon (1948) as a measure of information order in communications technology, the concept of negative entropy as a measure of signal focusing was extended and subsequently applied to diffraction imaging by De Vries and Berkhout (1984) by calculating the measure of focusing for a constant-velocity migrated data set to estimate a constant migration velocity for the whole data set. The concept can be extended to a 2D v(x, z) problem through the use of moving windows. This can be achieved by first applying a Hilbert transform to generate amplitude data a ij of amplitude a at pixel i, j, and calculating negative entropy focusing (S) by for each location throughout the data, where is the data amplitude normalized by the mean over an N1 × N2 moving window, and N = N 1 N 2 , the number of cells in the moving window. Each window is tapered with a Gaussian kernel such that only the central pixels of the window are equal to the input data. This approach using moving windows is highly computationally intensive, and we therefore calculate S using an alternate approach, by scaling all-positive amplitude data using a 2D automatic gain control (AGC) filter on which we calculate local negative entropy by where g(x, t) is the AGC gain function (Yilmaz, 2001) and a is the Hilbert-transform derived amplitude as previously. This results in one calculation to estimate g(t) within a moving window, followed by a point-wise calculation of S (Equation 7), and is significantly more computationally efficient. The result is a statistical measure of the degree of signal focusing for each location within the data; this metric is calculated for each constant migration panel. Figure 2 compares the performance of negative entropy to the continuous and windowed negative entropy function, semblance, and kurtosis measures of focusing. The diffraction is modeled at a depth of 100 m within a homogeneous medium of v = 0.165 m/ns, representing an anomaly within a cold ice medium. The diffraction is focused using a velocity range of 0.13-0.19 m/ns. We find that the negative entropy measure gives an improved metric of focusing. This is principally because negative entropy is less susceptible to the presence of multiple wavelet cycles typical of GPR data, as discussed by Booth et al. (2010), in addition to being less sensitive to regions of noise and crossing diffraction tails in the data. This is beneficial in the application of automated picking algorithms as in the subsequent step, where the focusing metric is used as an objective function to maximize 10.1029/2021JF006382 7 of 19 (Toldi, 1989), as it avoids the potential of small-scale local maxima which may lead to an incorrect velocity pick in the case of semblance or kurtosis (Figures 2c and 2d).
In practical terms, this seeks to find a surface in the trace-twtt axis through the 3D focusing volume generated in the previous step, which maximizes the sum of focusing across the surface. This generates a surface which is constrained by peaks in the focusing metric and enables a controlled interpolation between. The picking algorithm is undertaken in a two-step process. The first step maximizes the function in Equation 9 for each trace to give a coarse, strongly laterally variant velocity field. The maximum vertical gradient of the picked velocity profile is dictated by defining the picking gate, that is, the maximum variation in velocity between time samples, and a starting velocity is provided giving the surface velocity. In the absence of near-surface diffractions we estimate this velocity the velocity of radar through cold ice with a small air fraction, the calculation for which is discussed in Section 3.5. A second step undertakes smoothing through shaping regularization in trace and time coordinates , to limit the gradient. We then apply a time backshift Δt = 5.6 × 10 −9 s to the data, to account for the time difference between the wavelet maximum and first break . The derived velocity field is considered to be an average between the surface and the scattering point, which we henceforth refer to as V RMS . A smooth Dix inversion (Dix, 1955;Fomel, 2007), is then applied to V RMS to derive the local velocity, which we refer to as interval velocity or V local .

Migration and Depth Conversion
With the assumption of a constant velocity field, Stolt (Hambrey et al., 2005;Murray et al., 1997;Navarro et al., 2014;Saintenoy et al., 2013) or Kirchhoff migrations (Arcone et al., 1995;Schannwell et al., 2014;Sevestre et al., 2015) have been commonly applied to similar glacier GPR surveys on polythermal glaciers. In exploration seismic imaging, depth migration or reverse time migration (RTM) is undertaken when there are lateral variations in interval velocity due to the effect of energy refraction with depth (Yilmaz, 2001). Some studies have used RTM with airborne GPR data over temperate ice to reconstruct bed topography, where two constant velocity layers are required due to the elevated platform (Grab et al., 2018;Langhammer et al., 2019), yet this approach is not widely implemented due to the wavelength and model scales required. Assuming a migration aperture of 200 m, while we expect a degree of lateral variation in the velocity field, the magnitude and scale of velocity variations is minimal on the scale of the migration algorithm aperture as a result of the regularization used in the velocity picking step. For this reason, we apply Kirchhoff time migration followed by a 1D, trace-by-trace depth conversion, both using the derived V RMS model.

Water Content Inversion
Inversion for water content requires the assumption of a geophysical mixing model. Previous studies have used the Looyenga mixing model (Macheret et al., 1993;Murray et al., 2007), which makes the assumption of isotropic and spherical inclusions in a two-phase mixture, or the Complex Refractive Index Model (CRIM; Bradford et al., 2009;Brown et al., 2017) which allows inclusion of an arbitrary number of phase contributions by estimating a time-averaged slowness weighted using the relative volumetric contributions of each phase (Greaves et al., 1996). More advanced mixing models include the formulation of Giordano (2005), which considers the case where water is held within disc-shaped inclusions with a preferred orientation resulting in an anisotropy in velocity. This was used by Bradford et al. (2013) to quantify anisotropy on Bench Glacier, Alaska. Analysis of the bed pick at survey line intersections in our data collected on Von Postbreen shows most intersections are within a 50 ns (≈4 m at 0.166 mns −1 ) time difference, showing little evidence of significant anisotropy, hence we use the CRIM equation.
We treat the glacier as a three-phase material whereby pore spaces are made up of water or air pockets, and as such the water content by volume can be estimated by where ϕ w and ϕ a are the volumetric water and air fraction, respectively, v is the estimated local velocity, and v i , v w , and v a are the velocity of propagation through ice, water, and air.
The fraction of air in glacier ice, held within pockets and bubbles, decreases as a function of depth due to closure as a result of overburden pressure (Cuffey & Paterson, 2010), and must be estimated for Equation 10 to remove a possible depth-dependent bias in water content. Using the assumption of zero deviatoric stress, Bradford et al. (2009) balance the gas pressure within bubbles with the hydrostatic pressure due to an ice overburden to show that the volumetric air fraction ϕ a as a function of depth z, calculated at a depth discretized by Δz, is where g is acceleration due to gravity (9.81 ms −2 ), ρ the density of ice (918 kg/m 3 ) and K = ϕ 0 R/M, the volumetric contribution of air in ice at the surface scaled by R, the ideal gas constant and M, the molar volume of air at atmospheric pressure. β′ = 9.8 × 10 −8 K Pa −1 is a constant representing the rate of change of melting point with pressure, and P 0 is the atmospheric pressure at the surface. ϕ a(0) is the air fraction at the surface, which is assumed to be 0.1.

Uncertainty Analysis
Uncertainty in water content can be estimated with propagation of errors using partial derivatives (Bradford et al., 2009). We consider contributions from uncertainty in air fraction σ a , uncertainty in Dix-inverted interval velocity σ v , and uncertainty in cold ice velocity . σ a and σ v are estimated through partial differentiation of the CRIM equation (Equation 10); and the total error, assuming the error in v and error in ϕ a are uncorrelated, is The uncertainties in the radar velocity in water v w and in air v a are negligible, as these materials have well-defined real dielectric constants of 81 and 1 (Daniels, 2004), respectively, hence their contribution to uncertainty can be ignored for this analysis. The velocity of cold ice v i is estimated using the derived velocity above the scattering region and below 50 m depth (i.e., in cold ice and where ϕ a < 0.02) and is estimated from the distribution of the derived velocity. Consistent with the uncertainty estimations of ϕ a from Bradford et al. (2009), we hold the uncertainty in air fraction as = 0.5 which allows a significant variation in air fraction estimation given the simplified nature of ϕ a estimation.
We now wish to estimate the uncertainty in local velocity, first considering the effects of the width of peaks in the negative entropy peaks in the velocity panels (as in Figure 2b), and second considering the influence of propagating these uncertainties through the Dix velocity inversion step. The uncertainty in V RMS is estimated by assessing the width at half maximum of the negative entropy peaks derived from Equation 7. The negative entropy response in Figure 2b shows a width at half maximum of ±0.075 m/ns, which we use as an uncertainty going forward. Note that, as this is an uncertainty in V RMS , this can be used directly as an estimate of the uncertainty in bed topography correction. To estimate the uncertainty in interval velocity we must consider the implications of propagating this error in V RMS through a smooth Dix inversion algorithm. To do so, we use a Monte Carlo approach to velocity analysis using an example trace, similar to Booth et al. (2011), but using a continuous function of V RMS . We generate an ensemble of 100,000 traces with added random Gaussian noise with a standard deviation σ = 0.075 mns −1 , matching that of the V RMS range estimated above, smoothed over a scale factor of 10 pixels. We then undertake Dix inversion on all 100,000 noisy V RMS trajectories independently and calculate a 1D depth-dependent estimate of ensemble standard deviation. Due to the computational cost of this approach in 2D, we repeat this in 1D for a single trace and use the standard deviation in local velocity as representative for the wider data set; an alternative approach to uncertainty analysis with a lower computational cost would be desirable for operational use, however this approach is sufficient for this study.

Manual Validation of Velocity Field
We validate the results of velocity estimation by manually estimating the V RMS for survey line 5, taken across the glacier. We use the approach of Booth and Pringle (2016) for calculating the semblance response of individual hyperbolas and picking the semblance response of the first peak of the wavelet. Uncertainties on velocity are estimated using the width of the semblance response at 50% of the semblance peak. We pick semblance peaks for clear hyperbolas at two depths within the ice: along the uppermost surface of the scattering layer, and at or immediately above the ice-bed interface. For picks close to bed, we pick hyperbolas with no high-amplitude features immediately above, which may indicate multiples from sub-or englacial channels (Stuart, 2003) and hence focus at a lower V RMS than primary diffractions. We use the coherence-enhanced diffraction-only data set, as data with only pre-processing applied do not contain sufficient diffraction hyperbolas to enable this approach.

Diffraction Focusing
We apply the preprocessing, SVD filtering, and coherent wavefield separation and stacking to each acquired survey line and present the results from a radargram taken across-glacier in Figure 3. Figures 3a and 3b show the significant near-surface ringing in the pre-processed data, which was effectively removed by SVD filtering (Figure 3c) without lateral smearing or degradation of near-surface hyperbolas ("A" in Figures 3b and 3d). The SVD approach additionally handled changes in the noise characteristics along the survey line well ("B" in Figures 3b and 3d), and was successful in preserving the bed reflection throughout.
Coherent wavefield separation and stacking enhanced diffractions at the bed while removing planar bed returns, enabling estimation of V RMS for the full ice depth throughout much of the survey area. Figures 3e and 3f show the resultant diffracted wavefield. Coherent diffractions can be detected throughout the ice column within the scattering region (highlighted as "A" and "C" in Figures 3b, 3d, and 3f), and planar reflections have been removed. Some residual noise in the near-surface is additionally further suppressed in the diffraction-only data set. Figure 4 shows three example negative entropy panels from a segment of line 05 across the glacier, showing the correspondence of negative entropy peaks to prominent diffractions throughout the data. As there are no diffractions to estimate a starting velocity in the upper sections of ice, we use a starting pick velocity of 0.173 mns −1 , estimated through using Equations 11 and 10 (estimating v with ϕ w = 0) to represent the velocity in the uppermost ice across the glacier. Note that due to smoothing across traces, the picked V RMS profiles do not pass through all peaks in the panel but nevertheless provide a representative V RMS across the regularization aperture.

Velocity Picking and Validation
We perform a sensitivity analysis on the picked and Dix-inverted velocity field to select regularization scales by testing the results of a range of distance and twtt smoothing windows, after Nicolson et al. (2014). Figure 5 shows the result of testing a range of regularization windows for a single cross-line data set. Small regularization radii result in overfitting of Dix-inverted radar velocity (to the upper left of Figure 5). Cases with velocities below 0.1 and above 0.18 m/ns are rejected as these correspond to physically unreasonable velocities and are likely to be as a result of overfitting, and rapid transitions in average englacial water content are unlikely over short length scales. At large regularization scales, the velocity of the upper ice column reduces indicating an oversmoothing of low velocities into the higher-velocity upper layer, and a loss of distinctive cold/temperate ice distribution. In general, we find that a large smoothing radius is required to minimize rapid local transitions in velocity, which resulted in extreme local velocities. We used a window size of 100 m laterally and 50 samples in time (0.2 μs) across all lines for consistency, corresponding to panel (i) in Figure 5.
We manually pick the location of the apex of 13 diffractions across the top of the scattering layer and 11 diffractions at or immediately (<100 ns twtt) above the bed throughout the data. Figure 6 shows the extracted V RMS from the bed pick compared to the manually derived velocities. The uncertainty in Figure 6 is the uncertainty in V RMS as derived in Section 3.6. The manually derived picks do show some variations on short length scales which are not fully resolved by the automated approach (such as at 200 m in Figure 6b) but in general, the diffraction focusing-derived velocity field closely matches the manually derived picks to within the estimated uncertainty. Figure 6a presents an elevation-corrected radargram, produced assuming a constant englacial velocity field, along the full 13 km centerline profile through Von Postbreen. The profile excludes only the uppermost 900 m and lowermost 1.1 km of the glacier, where we could not gather data. On top of the radargram, we have superimposed a manually picked upper surface of the scattering zone, a surface that elsewhere has been taken as the interface between cold and warm ice (Bamber, 1987;Irvine-Fynn et al., 2011;Schannwell et al., 2014) and hence is hereafter  (Schwarz, 2019). Planar reflections are removed and the diffractions are enhanced ("A" and "C") and ringing is successfully removed through much of the profile ("B").

The Thermal Structure of Von Postbreen
termed the cold/temperate-ice transition zone (CTZ). This shows that the glacier has a typical Svalbard-or Scandinavian-type thermal structure, whereby temperate ice makes up much of the ice column in the accumulation area and is advected downglacier to the ablation zone, where a layer of cold ice is present close to the surface across the ablation zone (Blatter & Hutter, 1991;Irvine-Fynn et al., 2011).
Data collected in the upper accumulation area depict scattering with an upper surface between 50 and 100 m below the surface across the upper 2 km of Von Postbreen. We are unable to determine if this extends to the surface as ringing in the near surface has resulted in data clipping, hence SVD filtering in pre-processing would be unable to retrieve this near-surface scatter. Such scattering is expected in the near surface however due to heterogeneities within the snow and firn layers in the accumulation zone, and was present in surveys undertaken by Sevestre et al. (2015). The lower ablation region consists of a layer of cold ice of approximately 80 m thickness overlying temperate ice in the lower ablation zone (up to 2 km in Figure 7). Our centerline does not capture the glacier front, which is approximately 1.1 km from the survey start point. Figure 7a shows the CTZ to be coincident with the bed at the survey start point, and it is likely therefore that the glacier is cold-based for the lowermost 1 km.
There is significant topography in the CTZ along the glacier centerline, with high-amplitude scattering reaching close to the surface at 2 and 5 km in Figure 7a. These are likely to correspond to crevasses or moulins whereby surface meltwater is able to access the englacial hydrological system during the summer melt, potentially enhancing local englacial melt. No clear surface features were apparent when undertaking the survey. Further topography is seen at approximately 10 km, where a small high in CTZ topography is seen up-glacier of a significant bed topography high.
The diffraction focusing-derived radar-wave interval velocity field (V local ) along the same profile (Figure 7b) resembles structurally the contrasts in scattering depicted by Figure 7a. The upper ice column is characterized by higher velocities, generally 0.16 m/ns, while the lower ice column consists of lower velocities, generally ranging between 0.14 and 0.15 m/ns. The distribution of velocities above and below the manually derived CTZ is shown in Figure 8, reflecting the ability of the approach to differentiate the radar-wave velocities in the respective thermal regimes. Several regions of very low velocity (down to 0.13 m/ns) are found in the lower glacier (at approximately 3.8 and 5.4 km). The mean velocity below the scattering interface is 0.150 ± 0.004 m/ns, with the mean above being 0.159 ± 0.006 m/ns, although with a significant negative skew with a peak at 0.165 m/ns. This peak corresponds well with typical values for cold ice of 0.165-0.168 m/ns (Bradford et al., 2009;Navarro et al., 2014). Figure 7c shows the englacial water content distribution derived from Figure 7b. We estimate the uncertainty in englacial water content, following Equation 15, to be ±1.6%. Some regions of cold ice within the upper ice column have non-zero water content estimated, which may be a reflection of this uncertainty.

Glacier Bed Topography Reconstruction
We compare the results of undertaking bed topography estimation using an assumed englacial velocity of 0.166 mns −1 , as used by Navarro et al. (2014) to our newly derived radar velocity. Figures 9a shows the glacier-wide retrieved bed depth using picked velocities. Figure 9b shows the difference between the bed topography of Figure 9a and bed topography retrieved using the constant velocity of 0.166 mns −1 . The histogram in Figure 10 shows the distribution of differences between bed topography reconstructions across the glacier. We observe that ice depths derived using the picked velocity are consistently lower, with an average difference of 4.9 ± 2.5%, equivalent to 12.0 ± 8.4 m between constant-and picked-velocity bed topography.

The Velocity Structure of Von Postbreen
In Figure 7, we show the retrieved velocity structure with the interpretation of the upper scattering layer, which is representative of the interface between cold and temperate ice in a polythermal glacier. Figure 8 shows the distribution of local velocities below the interpreted CTZ. These both show significant variations in local velocity Figure 5. Testing picking regularization parameters using an example cross-line data set. Panels (a-c) show smoothing with the smallest (5-sample) window in the time-direction, with each subsequent row showing smoothing radii of 25, 50, and 100 samples. Panels (a, d, g, j) show smoothing with the smallest (25 m) window in the trace orientation (x-direction), with subsequent columns showing smoothing over 50 and 100 m smoothing. A low distance smoothing results in physically unreasonable velocity models with very rapid lateral changes in velocity. A low time smoothing parameter increases the chance of velocity spikes outside the expected range of 1.4-1.7 × 10 8 m/s after Dix velocity inversion. Data are clipped before the bed reflection to avoid fitting to subglacial diffractions. below the CTZ along the length of the centerline profile, with a range of 0.131-0.164 mns −1 . This is representative of lateral variations in water content through the glacier. Of particular note are areas of very low velocity at or down-glacier of regions of enhanced scattering (approximately 2 and 5 km in Figure 7a), where the upper scattering layer extends into the upper ice column. Given the low prevalence of crevassing across much of the glacier surface, this is likely to be a small moulin whereby water is accessing the glacier interior over the summer melt season and either (a) being held within the ice as part of a disconnected englacial drainage system or (b) enhancing local melt through release of latent heat. Our data were collected in spring, before significant seasonal melt began, hence this enhanced englacial scatter is likely to have been persistent through winter and represent seasonal storage of water. Similar undulations in the upper surface of the scattering zone have been observed in the presence of surface crevassing and moulins at Hansbreen, Svalbard (Moore et al., 1999), and seasonal changes in the vertical distribution of water content have been observed, whereby a decreased scattering intensity at the CTZ may be attributed to enhanced drainage through the summer melt season (Jania et al., 2005). A similar conceptual model of overwinter englacial water storage is proposed by Irvine-Fynn et al. (2006) using repeat GPR surveys over the polythermal Stagnation glacier, Canada, suggesting that the englacial piezometric surface is elevated in the early melt season due to an inefficient englacial drainage system. This may be consistent with observations of Hodson et al. (2005), who hypothesize that over-winter subglacial storage of water is necessary to complete the annual runoff and mass balance budget of Midre Lovenbreen, a smaller polythermal glacier with a low relative proportion of temperate ice by depth (Sevestre et al., 2015).
The presence of enhanced scattering through much of the region immediately below the CTZ may indicate a low hydraulic connectivity with the deeper subglacial hydrological network. At several locations in Figure 7, the inferred water content close to the bed is lower than in the overlying mid-ice column (for example, at approximately 2, 6, and 8 km in Figures 7b and 7c). Such a multi-layer structure has previously been observed using temperature profile measurements (Jania et al., 1996) and from GPR CMP velocity analysis (Murray et al., 2000), both inferring a 4-layer structure with temperate ice of water content of 3 − 5% overlying a low water content region close to the bed. Our results partly follow this structure, but not continuously throughout the glacier, breaking down in some regions of low scattering intensity (e.g., at 7 km in Figure 7b). The wider variations in water content imply that the glacier thermal structure may not be well-represented by sparse CMP-derived estimates of 1D velocity profiles.

The Nature of Scattering Bodies
Having elected to use the CRIM methodology of calculating englacial water content (Bradford et al., 2009;Brown et al., 2017), we assumed a 3-phase isotropic internal structure of englacial water content. The precise nature of englacial water is likely to be more complex than this. Some previous studies have assumed that water held within ice does so within intra-ice crystal veins, yet time-domain reflectometry of temperate glacier ice cores has shown that such veins only have a low capacity to store water, with an inferred upper bound of storage of <0.03% at 2° (West et al., 2007). Much water has instead been hypothesized to be stored predominantly in hydraulically disconnected macroscale pores, or voids, such as those imaged in the upper temperate ice of Storglaciaren using borehole video (Fountain et al., 2005). There, the majority of water storage and flow was found to be through fractures in the lower regions of ice, generally in the form of subglacial crevasses, at a similar orientation to the surface strain-related features and at sub-vertical (>70°) angles. This is consistent with Bradford et al. (2013), who showed through azimuthal GPR CMP surveys that velocity anisotropy is possible close to the bed of the glacier as a result of such subglacial features. Velocity anisotropy can be described by the retrieval of two dielectric constants parallel and perpendicular to the fracture orientation, known as ϵ ∥ and ϵ ⊥ , respectively. This approach suggests an anisotropic mixing model should be used, for example, as developed by Giordano (2005), which enables estimation of a ϵ ⊥ and ϵ ∥ for a distribution or arbitrarily oriented bodies of defined aspect ratio.
We can test for anisotropy by considering the differences in twtt at survey line crossovers. Our survey grid gives a total of 16 crossover locations, with 8 at near-perpendicular and 8 at approximately 50° intersection angles. Figure 11 shows the results of crossover analysis on depth-converted data for all 16 locations using both constant and picked velocity fields. While all of the intersections are derived from regions with very small variations in bed topography, two outliers are present (highlighted in the constant-velocity plot of Figure 11), which are the result of slight cross-line variations in topography, resulting in a twtt bias. Our results show very little anisotropy at radar line crossover points. Average mistie between crossovers for data with constant velocity migration applied is 16 ± 32 ns, with all mis-ties below 61 ns, showing no systematic bias, suggesting a low anisotropy between orientations. This implies that water may be held  in three possible configurations: (a) within macroscopic, disconnected pores, as predicted by West et al. (2007), with isotropic behavior, (b) within a series of subglacial fractures with a near-random distribution relative to the survey line such that the bulk dielectric constant has no azimuthal variation, or c), within fractures intersecting the survey lines by ≈45° such that ϵ ∥ ≈ ϵ ⊥ . Of these, (a) would support use of a spherical inclusion mixing model such as Looyenga or the CRIM model, and (b) and (c) would support treating the inclusions as fractures or disks as in Giordano (2005) with different distributions of orientations. Of these, (b) is unlikely as crevasses will generally be uniform in orientation (perpendicular to the direction of bulk ice deformation after Cuffey and Paterson (2010)), and (c) is unlikely given that anisotropy is not observed for profiles intersecting at both ∼90° and ∼50°. Brown et al. (2017) interpret a similar result with TMz (parallel antennas) and TEz (perpendicular antennas) oriented common offset surveys in Greenland, and as such use isotropic mixing models. Air voids or pockets are unlikely to be a source of scattering within the deeper ice as they would be likely to close rapidly due to overburden pressure (Cuffey & Paterson, 2010), and scattering is unlikely to be debris-induced due to the lack of surface debris.

Radar-Wave Velocity and Bed Topography
Our derived velocity structure results in a consistently reduced bed depth compared to the assumption of a cold ice glacier with a velocity of 0.166 mns −1 (Figure 9). This is expected in the sense that other authors have shown bed depths can be overestimated when assuming constant englacial velocities in polythermal ice (e.g., Murray et al., 2000). For many such glaciers, a practical workaround has been to assume that high radio-wave velocities in the firn counteract the lower radio-wave velocities through the temperate ice, such that a constant velocity can be applied in most instances (Jania et al., 2005;Navarro et al., 2014). While absolute glacier volume estimates are within a suitable error range using this approach (Lapazaran et al., 2016), the relative error in bed depth along the glacier profile remains unconstrained, resulting in possible underestimates where firn/cold ice dominates and overestimates in regions of predominantly temperate ice. The effect of such a systematic error is unclear when mass-conservation approaches to glacier volume estimation are used, such as in retrieving archipelago-wide bed topography estimates in Svalbard by Fürst et al. (2018). Our data from Von Postbreen show a consistently lower velocity due to temperate ice through the full accumulation zone ice depth, similar to previous observations at this site (Sevestre et al., 2015), in addition to Kongsvegen (Björnsson et al., 1996) and Hansbreen (Jania et al., 2005). As such, we cannot assume that the higher velocities expected within snow and firn will be sufficient to counter low velocities of temperate ice.

Conclusions
We have applied a novel approach to estimate the radar wave velocity within polythermal glaciers to derive a continuous profile of englacial water content. We have achieved this by firstly extracting the diffracted wavefield using a coherent diffraction extraction technique, followed by automated extraction of diffraction focusing velocities. This study represents the first such automated approach to this problem, building on previously developed manual approaches. The quantification of englacial water content has for the most part relied on the use of common mid-point or common source-point surveys, which are logistically challenging, hence only a sparse representation of spatially variant glacier-wide water content distribution can be achieved. Common-offset data are far more readily available, and can be acquired over a sufficient spatial scale to well-represent the glacier-wide distribution of water content.  We find that the velocities derived using this approach typically result in ice depths which are lower than those obtained using a commonly assumed constant velocity, with an average difference of 4.9 ± 2.5% of local ice depth, or 12.0 ± 8.4 m in absolute terms, for data collected over Von Postbreen. We have shown that the spatial distribution of low radar-wave velocity regions, and the corresponding volumetric contributions of englacial water content, are more complex than can be described by interpretation of the CTZ alone. For example, regions of high englacial water content can be found in the vicinity of near-surface scatter, which may represent enhanced local melt or seasonal storage. We have demonstrated that understanding lateral variations of englacial radar velocity is essential, and that diffraction focusing and velocity estimation are crucial in retrieving correct bed topography in the presence of temperate ice. The workflow developed for this paper provides a mechanism for capturing these englacial radar velocity variations that, with careful tuning of scale factors in the velocity picking and smoothing stages, can be applied to GPR data acquired from diverse radar systems and survey settings.