A framework for attenuation method selection evaluated with ice-penetrating radar data at South Pole Lake

Abstract

All radar power interpretations require a correction for attenuative losses. Moreover, radar attenuation is a proxy for ice-column properties, such as temperature and chemistry. Prior studies use either paired thermodynamic and conductivity models or the radar data themselves to calculate attenuation, but there is no standard method to do so; and, before now, there has been no robust methodological comparison. Here, we develop a framework meant to guide the implementation of empirical attenuation methods based on survey design and regional glaciological conditions. We divide the methods into the three main groups: (1) those that infer attenuation from a single reflector across many traces; (2) those that infer attenuation from multiple reflectors within one trace; and (3) those that infer attenuation by contrasting the measured power from primary and secondary reflections. To assess our framework, we introduce a new ground-based radar survey from South Pole Lake, comparing selected empirical methods to the expected attenuation from a temperature- and chemistry-dependent Arrhenius model. Based on the small surveyed area, lack of a sufficient calibration surface and low reflector relief, the attenuation methods that use multiple reflectors are most suitable at South Pole Lake.

1. Introduction

The measured power of a reflected radar wave preserves information from both reflector properties and path effects. Thus, for ice-penetrating radar, returning waveforms provide insight into the nature of subsurface interfaces as well as the ice column. Attenuation is generally the strongest of the path effects acting in polar ice that is well below the pressure-melting point. Radio-frequency attenuation is mostly due to the dielectric relaxations that occur due to lattice defects (Fletcher, 1970). Therefore, the attenuation rate, or the energy absorption as a function of travel path length, depends on the defect concentration and ice temperature (Corr and others, 1993; MacGregor and others, 2007; Stillman and others, 2013). Empirically-derived attenuation measurements in ice generally have high uncertainty, inhibiting comprehensive, ice-sheet scale, radar power interpretation (Jordan and others, 2016).

In glaciological applications, attenuation estimates are most commonly used as a correction for the power returned from the ice–bed interface (Matsuoka, 2011). Only after this as well as other corrections discussed below, can the measured bed power be interpreted in terms of bed properties; for example, roughness (MacGregor and others, 2013; Christianson and others, 2016) or wetness (Peters and others, 2005; Christianson and others, 2016), which both affect physical glacier processes, such as basal slip (Weertman, 1957; Lliboutry, 1968; Schoof, 2002). More recently, attenuation measurements have also been used to infer properties within the ice column, such as ice temperature (MacGregor and others, 2015) or water in the firn (Chu and others, 2018b).

Unfortunately, disentangling attenuation from reflector characteristics is challenging. Even carefully designed experiments, such as common-midpoint surveys, which eliminate the effects of spatial variability in ice temperature, impurity concentration and basal reflectivity by isolating differences in reflected power as a function of the travel path in the ice, often suffer from geometric (angle-dependent) effects. Common-midpoint surveys can therefore result in erroneously high attenuation rates (Holschuh and others, 2016).

Thermomechanical ice-flow models provide an alternative approach to determining attenuation through model-derived ice temperature (Matsuoka and others, 2010, 2012; Chu and others, 2018a) or data-model inversion (Jordan and others, 2016). However, model results will inherently be limited by the chosen model physics, numerical implementation and input data. Often, attenuation models only consider the ice temperature, ignoring the effects of reflector geometry and impurities (c.f. MacGregor and others, 2007, 2009, 2012).

Many empirical attenuation methods are used throughout the glaciological literature. Here, we explore the variety of methods, address their limiting assumptions and develop a framework for method selection that uses survey design characteristics and ice-dynamic setting to determine expectations for method performance. We test the framework by applying the suitable methods to a newly collected common-offset survey from South Pole Lake (SPL), showing that the variety of methods can yield varied results even when applied to the same input data. Finally, we discuss the benefits and limitations of each method. We aim for a comprehensive assessment that will inform future attenuation calculations, associated uncertainty, and best practices for radar survey design in glaciological studies.

2. Attenuation estimation methods

2.1 Empirical attenuation derivation

Every empirical attenuation method described here is based on a regression of diminishing radar power with depth. Measured radar power can be broken into its constitutive components with the radar power equation, generally presented in the glaciology literature using a simplified form which assumes ray theory can describe the radar wave (Bogorodsky and others, 1985),

(1)$$P_{\rm r} = P_{\rm t}\displaystyle{{A^2Rq} \over {4\pi {\lpar {2z} \rpar }^2B}}{\rm e}^{{-}2\alpha z}\comma \;$$

where P _r is the returning power measured at the receiver, P _t is the transmitted power, A is the antenna gain, R is the reflection coefficient (which depends on the magnitude of the permittivity contrast across the interface, the incidence angle of the wave, and the interface specularity (Schroeder and others, 2015)), q is refractive focusing (Bogorodsky and others, 1985, equation 3.8), z is the reflector depth, 4π(2z)² is a spherical spreading term (assuming the target is large relative to the first Fresnel zone), B is birefringence loss, and α is the dielectric absorption coefficient (attenuation coefficient). It is worth noting that the plane-wave approximation, necessary to implement Eqn (1), is only accurate for depths well below the antenna separation distance (Arcone, 1995). In the rare cases where near-surface power interpretations are desired, alternative spreading and refractive focusing corrections could be made for a more specific radiation pattern and firn density profile.

To simplify Eqn (1), we correct the power for spreading and refractive focusing, P _c = P _r4π(2z)²/q, group all the system properties into one term, S = P _tA ², and rewrite as

(2)$$P_{\rm c} = SR{\rm e}^{{-}2\alpha z}$$

ignoring any effects of birefringence for which corrections require polarimetric data (Fujita and others, 2006; Matsuoka and others, 2009). Birefringence losses only need to be considered when there is substantial phase lag between propagating waves of different polarizations, a phenomenon that requires both a well-developed crystal orientation fabric in ice and a relatively high-frequency radar system. To linearize the equation, power is commonly converted to decibels [dB],

(3)$$\lsqb {P_{\rm c}} \rsqb = \lsqb S \rsqb + \lsqb R \rsqb -2\alpha z\lpar {10{\rm lo}{\rm g}_{10}\lpar {\rm e} \rpar } \rpar \comma \;\;$$

where [X] = 10 log₁₀(X). Then, the one-way attenuation rate, in dB km⁻¹, is defined N = α(10 log₁₀(e)) (Winebrenner and others, 2003) so,

(4)$$\lsqb {P_{\rm c}} \rsqb = \lsqb S \rsqb + \lsqb R \rsqb -2Nz.$$

To solve for N, the depth derivative is carried out by linear regression of Eqn (4), with the regression slope being equal to the depth-averaged attenuation rate over the chosen sampling window. This procedure ignores any source or reflectivity changes with depth. Typically, the regression is done by either ordinary least squares (Jacobel and others, 2009) or weighted least squares using the uncertainty of the reflected power as weights (MacGregor and others, 2015); however, problems which have two measured variables (i.e. depth and power both subject to uncertainty) can also be solved by regression with errors in variables and reporting the uncertainty to some confidence interval, as we do here (Appendix).

In principle, empirical attenuation calculations all follow this same standard derivation, where some set of observations for the measured power and depth are isolated in Eqn (4) and regressed to calculate the attenuation rate (e.g. Fig. 1). In practice, the precise methods for arriving at a calculated value change in order to sufficiently isolate the depth–power relationship by satisfying the critical assumptions that (∂[S]/∂z) = (∂[S]/∂x) = 0 and |∂[R]/∂z| ≪ |∂[P _c]/∂z|. The validity of these assumptions is site and survey dependent. For some studies, a single coherent reflector is identified in many traces. If the reflector has spatially uniform reflectivity ((∂[R]/∂x) ~ 0) and has high relief (∂z/∂x), it can justifiably be matched across many traces to calculate attenuation. In other studies, when multiple reflectors in one trace are assumed to have nearly the same reflectivity ((∂[R]/∂z) ~ 0), they are used together to calculate attenuation for that trace individually. Lastly, in the select few surveys where a secondary bed reflection can be observed, it is compared directly to the primary reflection to calculate attenuation. We detail the methodology for each of these three groups below.

Fig. 1. (a) Map of the radar survey at SPL with surface elevation measured by GNSS, the lake outline in blue, and an inset map showing the survey location in Antarctica. The map projection is polar stereographic with the origin at the geographic South Pole (EPSG: 3031). (b) An along-flow radar profile corresponding to A-A’ in (a), with ice flowing from left to right. (c) Corrected power for every sample from colored points in (d). The red lines represents the depth–power regression from Eqn (4), with the slope of the line equal to the average attenuation rate over the depth spanned by the reflector(s). (d) Reflector interpretation for the along-flow profile shown in (b), with gray lines as internal reflectors. The colored points show corrected power for a single reflector over many traces (at the ice–bed interface) and for multiple reflectors from a single trace.

2.2 Method breakdown

2.2.1 Single-reflector methods

Method 1: single-reflector spatially-averaged attenuation. This is the most common method for calculating attenuation using a single reflector. A reflector is selected that can be easily identified across many traces, typically the ice–bed interface, and is used in a regression of reflected bed power as a function of ice thickness (i.e. Eqn (4)) over a spatial domain with minimal basal reflectivity variations (|∂[R]/∂z| ≪ |∂[P _c]/∂z|) (Winebrenner and others, 2003; Jacobel and others, 2009). The permittivity contrast between ice and the glacial substrate is large, so the measured power typically has a high signal-to-noise ratio (SNR). This method is limited by the assumptions that spatial variations in bed roughness and/or permittivity contrast (i.e. wetness) are small and that spatial variations in the attenuation rate are small (|∂N/∂x| ≪ 0) in the domain selected for regression. The resulting attenuation rate from this method is commonly reported as a depth-averaged value for the full ice column; however, using the bed reflector preferentially samples power changes due to more or less deep ice, far from the ice-sheet surface which acts as a fixed (cold) thermal boundary. Assuming that spatial variations in attenuation rate are small, the attenuation rate above the shallowest bed reflection is uniform, and the regression of bed power against changes in range represents the attenuation rate of the deepest ice only. In most real cases, some amount of spatial variation in the depth-temperature profile biases the calculated attenuation rate, violating the method's assumptions and resulting in an ill-posed problem. Importantly, the bias can just as easily be in either direction depending on the glaciological conditions (i.e. toward or away from the depth-averaged interpretation which has commonly been used) (see Supplementary Material).

Method 2: single-reflector depth-resolved attenuation. This method is similar to Method 1 in that it regresses power as a function of depth over many traces, but it uses all the internal reflectors observed throughout the ice column. The value from each reflector represents an interval attenuation rate, so taken together the result is a depth-attenuation profile. This method is analogous to that used for common-midpoint surveys (Holschuh and others, 2016). Unlike the bed, internal reflectors are typically specular, and are less likely to have spatially variable electrical properties over short wavelengths, so the |∂[R]/∂z| ≪ |∂[P _c]/∂z| assumption (also required for Method 1) is more likely valid. Unfortunately, the measured power from internal reflectors generally has a much lower SNR, and specular reflectors exhibit stronger slope-dependent changes in returned power than do rough reflectors (Holschuh and others, 2014). Without corrections (e.g. MacGregor and others, 2015), both the lower SNR and slope effects make it difficult to correctly calculate the attenuation rate, but these effects can easily be overlooked by the formal uncertainty metric.

Method 3: single-reflector spatially-resolved attenuation. This method uses the bed reflector to adaptively constrain attenuation over a spatial window of varying size, optimizing the window so that it is large enough to capture a sufficient depth range but small enough to minimize variations in reflector properties that introduce uncertainty into the calculation (Schroeder and others, 2016a, 2016b). By iteratively moving the window through the survey domain, this method can be used to constrain spatial variation in attenuation. In the description of this methodology, Schroeder and others (2016b) define the optimal attenuation rate as that which most sufficiently decorrelates the attenuation-corrected power from the observed ice thickness. Their derivation is analogous to that which we laid out above; they compute values of N such that (∂([P _c] + 2Nz)/∂z) = 0, again assuming that ∂[S]/∂z and ∂[R]/∂z are negligible. They perform a grid search through physically realistic attenuation rates instead of computing N through inversion as we do here. As in the case of Method 1, we assert that under the given assumptions this method represents an interval result over the span of the bed reflector (see Supplementary Material).

Method 4: known reflectivity. This final single-reflector method uses an interface of known reflectivity (e.g. an ice–seawater interface) to calculate the average attenuation rate through the full ice column by constraining all of the correction terms in Eqn (1) (Shabtaie and others, 1987; Bentley and others, 1998). Implementing this method requires high confidence in the measured, inferred or assumed values of the transmit power as well as the corrections: S, R, q, 4π(2z)² and B. The greatest uncertainty is generally associated with the reflectivity correction; for example, even for an ice–seawater interface, the reflectivity can be greatly reduced by roughness effects (Christianson and others, 2016).

2.2.2 Multiple-reflector methods

Method 5: multiple-reflector depth-averaged attenuation. This method uses many reflectors within one trace for a regression of reflector power and depth, again using Eqn (4). Normally, a majority of the ice column is included in a single regression, excluding only reflectors very near the surface or bed (MacGregor and others, 2015). Using this method, attenuation rates can be calculated for each individual trace, resulting in a sense for spatial variation in depth-averaged attenuation. It should be noted that this method heavily relies on the assumption of constant reflectivity between reflectors. However, previous work has shown that the magnitude of internal reflections can vary substantially, as they represent various dielectric contrasts across interfaces in the ice column (Eisen and others, 2003a; Bingham and Siegert, 2007). Corrections can be made to account for each source of reflectivity variation; for example, Macgregor and others (2015, Appendix A) use dielectric profiles from nearby ice cores to estimate the reflectivity within the ice column and correct the attenuation rate.

Method 6: multiple-reflector depth-resolved attenuation. This method uses many reflectors over a depth-constrained vertical window that does not span the entire ice column (Matsuoka and others, 2010); thus, the result is an interval attenuation rate as opposed to a depth-averaged value. This method can be implemented for each individual trace, as in Method 5, but is more effective for resolving depth variations when many traces are considered together. Theoretically, an adaptive window-fitting technique, as in Method 3, could be used here. Practically though, the window size is effectively constant because it is constrained by the limited thickness of the ice column. Violations in the depth-reflectivity assumption have induced large uncertainties in previous applications of this method. For example, depth-varying attenuation rates computed by Matsuoka and others (2010) showed decreasing attenuation rate with depth, opposite to that expected based on the increasing ice temperature with depth.

2.2.3 Secondary-reflection methods

Method 7: secondary-reflection spatially-averaged attenuation. This method uses the power difference between the primary and secondary bed reflections to calculate attenuation through the full ice column (MacGregor and others, 2011; Christianson and others, 2016). The ‘secondary reflection’, often called the ‘multiple’ reflection, is a repeated instance of the bed in a single trace. As the initial reflected wave is recorded by the receiver, a secondary wave is generated at the ice surface by reflecting off of the firn–air interface, propagating back into the subsurface and producing a second reflection of the bed that is imaged at the surface after twice the travel time. Attenuation is calculated by comparing the measured power between these two reflections. Secondary reflection methods, often called ‘multiple-echo methods’, require knowledge of both the ice–bed and firn–air reflectivities, but, unlike all prior methods, do not require any assumption about the system characteristics. When the reflectivities are assumed to be well known (i.e. a smooth ice–seawater interface) and attenuation calculations are spatially averaged over many traces, high confidence can be assigned to the result.

Method 8: secondary-reflection spatially-resolved attenuation. In exactly the same manner as Method 7, this method uses the primary and secondary bed echoes to calculate attenuation. Here though, attenuation results are reported on a trace-by-trace basis and contrasted between locations.

2.3 Framework for method selection

The necessary assumptions for each method (summarized in Table 1) provide a basis for attenuation method selection. Information about the survey design, instrument characteristics and the glaciological environment can provide an a priori basis for selecting a method, as different assumptions are likely to be violated for surveys in fast-flowing ice versus slow-flowing ice, for surveys with data predominantly collected along-flow versus across-flow, and in surveys that cover a large spatial area versus those that are more spatially confined. In the section below, we outline the survey characteristics that can either enable the use of certain methods in the case of a calibration surface, or inhibit the use of certain methods due to a violation of the central assumptions.

Table 1. A framework for empirical attenuation methods

Note that in all cases where the interface reflectivity is unknown (methods 1, 2, 3, 5 and 6), the critical criterion to be met is |∂[R]/∂z| ≪ |∂[P _c]/∂z|. To meet this criterion, single-reflector methods (1, 2 and 3) require an assumption that reflectivity changes along the reflector are small and that the reflector has sufficient relief (∂z/∂x) for the regression. Multiple-reflector methods require an assumption that reflectivity changes between reflectors are small.

Methods 4, 7 and 8 each utilize a calibration surface of known reflectivity. While these are among the best constrained methods available, they require unusually uniform glacial environments. These cannot be used without data either focused on a region containing floating ice (on an ice shelf or over previously known subglacial lakes with spatially uniform basal reflectivity) (MacGregor and others, 2011; Christianson and others, 2016), or surveys that deliberately connect to one of these regions before going to their true target. Applicable to the narrowest range of environments, Methods 4, 7 and 8 are most useful for small surveys nearby a known reflector, and generally less useful for large surveys over which reflectivity changes are expected or small surveys focused in the ice-sheet interior where no calibration surface can be expected.

Even for surveys without an obvious ice–water interface, the bed reflector is still the most widely used for attenuation calculations (Method 1; Jacobel and others, 2009). Unfortunately, the bed is also the reflector with the highest likelihood of substantial reflectivity variations. This is especially true in cases where heterogeneous water systems are present at the bed, or where unresolvable bed roughness reduces coherent scattering (Schroeder and others, 2014; Christianson and others, 2016; Jordan and others, 2017). If bed roughness is inversely correlated with bed depth, such that shallower sections of a profile reflect less power back toward the instrument receiver, the assumption to minimize reflectivity changes along the reflector ((∂[R]/∂x) ~ 0) is violated, leading to an erroneously low attenuation rate. Where power losses due to roughness are expected, corrections can be made; for example, using an aggregated bed-returned power (Oswald and Gogineni, 2008; Jordan and others, 2016) or by Kirchhoff Theory (Ogilvy, 1991; Peters and others, 2005; MacGregor and others, 2011). With no correction though, single-reflector methods (i.e. 1, 2 and 3) are not well suited for surveys where reflectivity changes are likely along the reflector. These cases will show an attenuation bias toward the reflectivity gradient. Surveys that span a larger area are better suited for Methods 1, 2 and 3, as the effect of heterogeneity in reflector property is reduced with more spatial averaging. Additionally, for these surveys that span a large area, reflectors will generally sample a larger depth range, leading to a more robust regression that represents the attenuation rate over more of the ice column. Both of these advantages are displayed in the correlation metric that underlies Method 3, which often requires a large area to meet the correlation strength minimums.

Dipping internal reflectors have lower reflected power for geometric rather than attenuative reasons (Holschuh and others, 2014). Deeper reflectors conform more closely to the bed topography and will therefore have greater dips (Hindmarsh and others, 2006). In addition, layers imaged along-flow tend to be smoother than layers imaged across-flow, as their shared strain history limits the steepness of layer folds to those that can be produced by englacial vertical velocities (Holschuh and others, 2017). Because dip and power are both depth dependent, geometric losses interfere with attenuation calculations and steeply dipping reflectors cannot be used for attenuation calculations unless corrections are made (MacGregor and others, 2015). Unfortunately, the environments where spatial changes in attenuation can be expected are often those where the layers are most distorted, such as ice-stream shear margins (Holschuh and others, 2019). Thus, in fast flowing ice, Methods 2, 5 and 6 are best applied to data collected along-flow, and in areas with relatively low shear-strain rates.

Paleoclimate transitions can create strong reflectivity changes between reflectors within the ice column. As the atmospheric composition changes, the dust content in depositional snow layers changes which can greatly affect the bulk reflectivity of an interface (Eisen and others, 2003a). For example, in Greenland, the Holocene-glacial transition represents a stark change in reflector character, separating bright reflectors near the surface and dimmer glacial layers near the bed (Karlsson and others, 2013; MacGregor and others, 2015). Using data with such a depth trend will produce attenuation rates that are anomalously high. Methods 5 and 6 are best used in ice deposited during a single climate regime, in ice that preserves many distributed climates with depth, or in regions where changes in dust and salt concentrations do not drive large changes in ice acidity. Using Methods 5 and 6 with radar data that span a single climate transition would require carefully correcting the reflectivity, as is shown by MacGregor and others (2015, Appendix A). The shallowest reflections, from within the firn, are generally associated with density rather than chemical contrasts, so reflectivity variation is expected between those and deeper reflections. Moreover, the magnitude of density contrasts fall as the background density approaches that of solid ice; thus, the reflectivity changes with depth throughout the firn column as well.

Finally, reflectors deep in the ice column are often close to the system noise floor (low SNR) and can therefore show anomalously high variation in measured power. When considered together with other reflectors in the same trace, noisy reflectors can then introduce what looks like an anomalously bright (dim) reflection low in the ice column and lead to a lower (higher) calculated value that is not representative of the true attenuation rate. For Methods 5 and 6, other authors have generally excluded very deep reflectors (MacGregor and others, 2015) and we recommend doing the same.

With this as a framework, we use the site and system characteristics of a new survey from SPL to evaluate the applicability and likely performance of different attenuation estimation methods.

3. South Pole Lake: survey characteristics

SPL is a small subglacial lake ~10 km from the geographic South Pole. The lake was first identified by airborne radar (Carter and others, 2007), but its presence has since been cross-validated by active-source seismic imaging (Peters and others, 2008), which indicated that the water column was up to 30 m deep. The ice-surface manifestation of the lake is a 10-m depression over an area of ~35 km². Our radar survey consists of 15 profiles spanning the surface depression (~150 line-kilometers) that were collected during the 2018–19 austral summer (Fig. 1a). As the surface velocity is only ~9 m a⁻¹ (Lilien and others, 2018) and varies by only a few meters per year on and off the lake, we do not expect strong spatial variation in ice temperature or crystal orientation fabric. Similarly, we do not expect spatial variations in the impurity concentration because the survey extent is small relative to the depositional catchment. Thus, we expect minimal spatial variations in radar attenuation for this survey.

The radar system used here is a ground-based 3 MHz impulse radar (Welch and others, 2009; Christianson and others, 2012, 2014). The transmitter is a commercial pulse generator produced by Kentech Instruments that operates at a 1 kHz pulse-repetition frequency with a pulse amplitude of 4 kV. The receiver is a two-channel, 200 MHz, digitizing PCIe board made by GaGe Applied Electronics that is mounted in a ruggedized computer. The two channels can be used to collect low- and high-sensitivity data simultaneously by distributing the available 16 bits over a different dynamic range. In this case, the high-sensitivity channel was set to ±100 mV and the low-sensitivity channel was set to ±1 V. The antennas are resistively-loaded dipoles that are towed across the snow surface.

Site and system characteristics help winnow the list of applicable methods. This survey covers a relatively narrow region (~20 km by 20 km) over relatively thick ice (~2800 m), it includes profiles both along flow and highly oblique to the flow direction, in a slow-flowing, low-shear region. The lake is up to 30 m deep, but likely shallower over much of its area. Because the imaging wavelength (56 m) is large relative to the resolvable depth of the lake (resulting in power variation over the lake surface), the survey lacks a suitable calibration surface. Additionally, a secondary bed reflection was not observed due to the thickness of the ice column. For those reasons, Methods 7 and 8 cannot be applied here. Likewise, we cannot apply Method 4 because the transmit power for this system is poorly constrained. The remaining methods to be applied include 1, 2, 3, 5 and 6. The limited aerial coverage of the survey, and resulting limited relief in imaged layers, both prevent the use of Method 3 (which required windows of >50 km at Thwaites Glacier to converge) and reduces the quality of estimates from Methods 1 and 2, which require regional averaging and layer relief to produce reliable attenuation estimates. We expect that Methods 5 and 6 will produce the best results for this dataset.

4. South Pole Lake: radar data processing

Data processing is done using ImpDAR, an open-source software package for processing impulse radar data (Lilien and others, this issue). The processing flow follows Christianson and others (2016) including a bandpass filter (fifth-order Butterworth filter from 1 to 5 MHz), time correction for antenna spacing and wave-speed variations through the firn column (permittivity profile from Kravchenko and others (2004)), interpolation to a standard trace spacing of 3 m, a horizontal de-meaning filter and two-dimensional time-wavenumber migration. The de-meaning filter subtracts a moving-averaged trace with a 100-trace window and scaled toward no correction below 500-m depth to emphasize the removal of horizontal artifacts near the surface.

Before interpreting radar power, bright reflections must be extracted from the trace because any low amplitude returns that do not form coherent reflections are more likely to be due to system noise rather than changes in the electrical properties of the ice. Bright reflectors are extracted in two ways (Fig. 2). First, we identify coherent reflectors using the digitizer implemented in ImpDAR (Lilien and others, this issue). The power of an identified wavelet is calculated from the root-mean-square (RMS) amplitude between the two troughs of the wavelet (Fig. 2a). Second, wavelet peaks are automatically selected based on a power threshold, set as the 97–99th percentile of individual sample power for a moving window (Fig. 2b) (as in Matsuoka and others, 2010). The window is 5 wavelengths tall in order to skip past a few potentially dim reflectors. Trace samples selected by this method provide a secondary measure of reflector strength which does not require manual interpretation of internal reflectors. These are referred to below as ‘threshold samples’. After isolating bright reflections, the measured power is corrected for known physical effects. Spherical spreading and refractive focusing in the firn column are corrected as in Eqn (1). Refractive focusing is computed following Bogorodsky and others (1985, Eq. 3.8) with a permittivity profile defined by measurements from the RICE detector at the South Pole (Kravchenko and others, 2004). Below ~150 m, the permittivity profile is constant; thus, focusing corrections are small (as shown by the ray-tracing methods described in Holschuh and others (2016)). Roughness in the reflector reduces the reflectivity (i.e. rough reflectors are diffuse rather than specular (Jordan and others, 2017)). We correct for roughness using Kirchhoff theory (Ogilvy, 1991; Peters and others, 2005; MacGregor and others, 2013). Finally, uncertainty in measured power $\lpar \sigma _{{\rm P}_{\rm c}}\rpar$ is defined layer-by-layer, using the mean difference between crossing profiles within one Fresnel zone of the cross-over locations. Relative uncertainty in reflector depth is defined in the same way (σ_z = 1 m), although this value is admittedly difficult to know precisely. The absolute depth uncertainty is generally higher (~quarter wavelength) and even greater without a careful firn correction. However, the relative uncertainty between measurements is the critical value for the attenuation regression and is determined here by cross-over analysis.

Fig. 2. Power extraction for bright reflectors. (a) An example wavelet selection for determining the power of a coherent reflector. The green-colored area is the picked wavelet with extents indicated by dots at the minimum (or maximum if reversed polarity) amplitude bounds. This area is used to calculate the root-mean-square amplitude. (b) The trace power (black line), interpreted RMS power (green dot), peak power from the isolated wavelet (empty black dot), and power from threshold samples at 97–99th percentile (purple dots), all corresponding to the wavelet in (a).

After processing and power corrections, we apply the appropriate attenuation methods at SPL. First, we implement Methods 1 and 2, calculating an individual attenuation rate for each observed reflector. This is done over the entire survey area by matching reflectors between profiles. Second, we implement Method 5, using all internal reflectors within each individual trace to calculate attenuation for that trace alone. Third, we implement Method 6, using every sample from all reflectors together to calculate vertical variation in attenuation over the entire survey area. These vertical attenuation variations are calculated within a vertically moving window, repeated for windows ranging in size from 5 to 15 wavelengths. One particularly dim packet of reflectors at ~900–1000 m, which corresponds to the glacial transition ~10 ka (Casey and others, 2014), is removed before any multiple-reflector attenuation calculations (Methods 5 and 6).

5. South Pole Lake: results

Methods 1 and 2 give the attenuation rate for each interpreted reflector (Fig. 3). The plotted attenuation rate is representative of ice over the depth range that the reflector spans. The resulting attenuation rates generally increase with depth but with high variability. Curiously, we observe many negative attenuation rates, especially near the ice surface (Fig. 3c); this result is non-physical. Of the reflectors that meet our uncertainty criteria (i.e. below 0.3 dB km⁻¹), the bed (Method 1) shows the highest attenuation rate, 16.7 ± 0.2 dB km⁻¹ (Fig. 3e), and the lowest attenuation rate is −20.9 ± 0.3 dB km⁻¹ from a reflector at 800 m depth. The calculated uncertainty is highest for the shallowest reflectors (up to ~2 dB km⁻¹), but rapidly decreases toward lower uncertainty with depth and is generally <0.3 dB km⁻¹ for reflectors deeper than 1 km.

Fig. 3. Single-reflector depth-resolved attenuation using Methods 1 and 2. (a) Corrected RMS power for every sample from all coherent reflectors. Reflectors are matched between profiles and colored by mean depth. (b) Cross-over error ($\sigma _{{\rm P}_{\rm c}}$), (c) calculated attenuation rates, and (d) corresponding uncertainty calculated for each individual reflector. In each panel (b, c and d) points are faded where the uncertainty exceeds 0.3 dB km⁻¹. (e) The depth and power data (black) and regression (red) for the bed reflector (Method 1). (f) A map of the corrected power from the bed reflector corresponding to (e), plotted here in polar-stereographic coordinates.

Method 5 gives trace-by-trace attenuation rates from both RMS power and threshold samples (Fig. 4). The results using RMS power are normally distributed around 3.5 ± 1.1 dB km⁻¹. The threshold samples give slightly higher attenuation rates and lower uncertainty, 4.4 ± 0.3 dB km⁻¹. The std dev. of the histogram for reported attenuation rates is within the calculated uncertainty for RMS power (0.5 dB km⁻¹) and only slightly greater than the calculated uncertainty for threshold samples (0.6 dB km⁻¹) (Figs 4c and d). As expected, the total range in reported values from Method 5 (i.e. the spatial variation in attenuation over SPL) is minimal (~ ± 2 dB km⁻¹).

Fig. 4. Multiple-reflector depth-averaged attenuation using Method 5 for each individual trace. (a and b) Corrected power for both RMS power from continuous reflectors (a) and threshold power from bright samples (b). Red lines show the regression for each and the black points show power for the entire trace. (c) Histograms of the calculated attenuation rates for each trace throughout the surveyed area with green corresponding to RMS power in (a) and purple corresponding to threshold power in (b). (d) Histograms of the regression uncertainty for each trace as in (c).

Method 6 gives depth-resolved attenuation rates using multiple reflectors from both RMS power and threshold samples (Fig. 5). Low and sometimes negative attenuation is observed near the surface, with the lowest result for RMS power being 0.2 ± 0.01 dB km⁻¹, and for threshold samples −0.5 ± 0.05 dB km⁻¹. The calculated attenuation rates climb as the fitting window moves toward the bed. This continues until the window gets to ~2 km depth, below which there are not coherent reflectors to implement this method. At the deepest window, the calculated attenuation rate by RMS power is 7.9 ± 0.01 dB km⁻¹ and by threshold samples is 9.1 ± 0.05 dB km⁻¹. Both the variation in attenuation and the calculated uncertainty are greatly reduced with increasing window size (values reported in the text are using the 15-wavelength window).

Fig. 5. Multiple-reflector depth-resolved attenuation using Method 6. (a) Attenuation rates calculated from both the RMS power (green) and the threshold samples (purple). This calculation was repeated with the windows of 5, 10 and 15 wavelengths, where the plotted point opacity represents window size, opaquer being a larger window. (b) Corresponding uncertainty for all points in (a). (c) RMS power for all reflectors and all profiles throughout the survey (black). An example regression line (red) is shown for the points over a selected window (green).

We tested two different techniques for isolating ‘bright’ reflectors, one by calculating the RMS amplitude of an isolated wavelet and another by thresholding bright samples within a window of the trace (Fig. 2). Both give similar depth and power for the bright reflectors (Figs 4a and b). The attenuation rates are also similar for each (Figs 4c and 5a), but with slightly higher uncertainty in the RMS case (Figs 4d and 5c).

6. Discussion

Even though we observe some discrepancy between the attenuation methods applied here (Fig. 6), after discarding the results which do not meet our uncertainty criteria, the calculated attenuation rates approximately agree with what is expected based on the predicted attenuation from an Arrhenius model (MacGregor and others, 2007) (Fig. 6). Single-reflector methods (i.e. Methods 1 and 2) result in a highly variable attenuation-depth profile that loosely represents the modeled profile but only for the deepest reflectors. The bed reflector (Method 1) gives the highest reported attenuation rate, which is expected because it is sampling the deepest and warmest ice. Multiple-reflector methods (i.e. Methods 5 and 6) are characteristic of the ice column in a bulk sense, but tend to misrepresent attenuation rates near the surface and bed.

Fig. 6. A comparison of all attenuation calculations and corresponding uncertainty from this study. (a) Attenuation results replotted from Methods 1 and 2 (Fig. 3), Method 5 (Fig. 4) and Method 6 (Fig. 5). The vertical lines represent the mean values from Method 5 over the full depth interval spanned by all internal reflectors. Only results from the 15-wavelength window are shown for Method 6. The black dashed line is the modeled profile from a temperature- and chemistry-dependent Arrhenius model (MacGregor and others, 2007) with uncertainty in the grey shaded region. Temperature input for the model is a Robin (1955) solution with air temperature −51.5°C, accumulation 8 cm a⁻¹, and geothermal flux 75 ± 10 mW m⁻², which generally approximates the nearest measured temperature profile (Price and others, 2002). Approximate ion concentrations are used for the entire ice column: the ice acidity [H ⁺] = 1 ± 0.5 μM; and the sea salt concentration [ss Cl ⁻] = 3 ± 0.5 μM (MacGregor and others, 2009, 2007; Matsuoka and others, 2012). (b) Regression uncertainty plotted against the degrees of freedom (number of points in the regression) for each of the methods implemented in this study (dots) and a profile of the uncertainty coefficient from Eqn (A3) to show that the number of measured points in the regression strongly controls the reported uncertainty.

For the survey presented here, spatial variation in attenuation rates is expected to be exceedingly small. There is no hydrologic or ice-dynamic feature that would cause a spatial gradient in ice temperature (Hills and others, 2017; Holschuh and others, 2019), and the impurity concentration should predominantly vary in the vertical, even at the ice-sheet scale (Siegert and others, 1998; Matsuoka and others, 2012). Hence, it was unsurprising that the trace-by-trace attenuation results were normally distributed with a std dev. close to the calculated uncertainty (Fig. 4). Perhaps more interesting is the agreement to within one std dev. between RMS power results (which require semi-automated interpretation) and threshold sample results (which can be extracted automatically). Manually interpreting reflectors is time intensive. Our result points out that this manual interpretation could potentially be avoided in some studies, as long as the SNR for reflectors is high: that is, above the noise floor for reflectors shallower than ~0.85 of the ice thickness.

Depth-resolved attenuation results (Methods 2 and 6) show increasing attenuation rate with depth which is expected based on the strong temperature dependence. This result is encouraging given that prior attempts to resolve attenuation measurements with depth have shown the opposite result (i.e. decreasing attenuation with depth) (Matsuoka and others, 2010; Holschuh and others, 2016). On the other hand, results from Method 2 are significantly lower than expected (even negative) for most of the reflectors close to the surface, and there is high variance between results for all layers throughout the ice column. Method 2 is an adaptation from the common-midpoint method outlined by Holschuh and others (2016). Admittedly, it did not work well for this common-offset survey; however, results from this method are more encouraging at other locations even with the same radar system (i.e. data from Northeast Greenland Ice Stream (Christianson and others, 2014) and unpublished data from Hercules Dome). One obvious explanation for the counterintuitive results from Method 2 is that the reflectors are too flat near the surface, so using a regression to calculate attenuation is inappropriate because there is insufficient depth variation. This effect can be seen by the increasing uncertainty for near-surface layers (Fig. 3d). However, we also see that the power between reflectors is near-constant for the first ~1000 m below the surface (Fig. 5c), causing slightly negative attenuation rates even from the multiple-reflector methods. An alternative explanation is that there are reflectivity changes between reflectors, which violates the assumption that (∂[R]/∂z) ~ 0. While it would be tempting to make reflectivity corrections using dielectric profiling data from nearby ice cores (MacGregor and others, 2015, Appendix A), prior efforts to do so at other locations have proven challenging (Eisen and others, 2003a, 2003b, 2006). Moreover, these dielectric profiling experiments have not yet been done for the South Pole Ice Core.

The reported uncertainty is much lower than expected. In many cases, the uncertainty for a given method is less than the variance observed between results for that same method; this is especially apparent in Method 2 (Fig. 3). Where this is true, some violation of the assumption, |∂[R]/∂z| ≪ |∂[P _c]/∂z|, must be introducing a hidden bias in the regression that is not represented in the uncertainty calculation. Instead, the uncertainty is most strongly dependent on the number of measurements included in the regression (Fig. 6b). As an alternative uncertainty metric, Schroeder and others (2016b) scale their fitting statistics to the attenuation cross-over error, although this value will still be subject to similar hidden biases.

7. Implications for radio thermometry

The primary advantage of radar sounding is that it can infer subsurface properties without direct sampling. Because radar attenuation in ice is so strongly dependent on temperature, one somewhat elusive goal has been to measure temperature with radar. While depth-averaged temperature results have been reported (MacGregor and others, 2015), published results for depth-resolved attenuation and corresponding ice temperature have given counterintuitive or non-physical results (Matsuoka and others, 2010; Holschuh and others, 2016). Initial attempts with seismic data have been more promising (Peters and others, 2012). Unfortunately, seismic data acquisition requires intensive ground activities, and it is thus not suitable to apply at the ice-sheet scale in the way that radar methods are. Theoretical considerations suggest that radar methods are likely to be as capable as seismic methods with proper experimental design, especially when accounting for ice chemistry (MacGregor and others, 2007). Looking forward, it is important to assess how precisely we can measure attenuation (and associated temperature changes) for a given empirical method.

Prior attenuation studies use variations on the methods described here, both to calculate the attenuation and to constrain uncertainty (Table 2). Published attenuation uncertainty ranges from ~0 to 8 dB km⁻¹, with a few outliers up to 89.9 dB km⁻¹. The uncertainty in the associated ice temperature is up to ~15°C for pure ice at ~−30°C. The reported uncertainty from our study is on the lower end of this range, 0.01–1.6 dB km⁻¹, as are most that report uncertainty from a regression result. It should be noted though that the regression is often calculated using thousands of points (as is done here), and that many studies only report the regression uncertainty with no direct consideration of uncertainty in the regressor variables, both of which lead to lower uncertainty estimates. If any of the sources of uncertainty discussed within our outline of the attenuation framework act to violate the assumption that |∂[R]/∂z| ≪ |∂[P _c]/∂z|, the resulting uncertainty can be significantly greater than predicted by the regression, and even lead to a non-physical results as we show in Fig. 3. Any direct interpretation of ice temperature should come only after correction for any possible hidden biases.

Table 2. Reported attenuation results and associated uncertainty from this and prior studies

~Value taken from the graph.

^a Reported uncertainty does not include any uncertainty estimate for the measured power.

^b High attenuation rates associated with geometric effects of the common-midpoint survey.

^c Many of the reported uncertainty values from this study are lower (in the approximate range 1–10 dB km⁻¹).

8. Conclusions

In this study, we highlighted the site and ice-penetrating radar survey characteristics that affect attenuation estimation, and used those to develop a framework for method selection. To test that framework, we used a recent radar survey from SPL, which shows some discrepancy in attenuation results between methods. While each of the described methods has been implemented before, we provided the first comparison for these methods using a single radar dataset.

Suggestions for method selection based on our framework are confirmed by the results from SPL. When reflectivity changes are expected between reflectors, but the relief of any given reflector is high, single-reflector methods (i.e. Methods 1, 2 and 3) are preferable. At SPL, we observed low reflector relief and correspondingly erratic attenuation results from Method 2, with the deepest reflectors falling in a physically realistic attenuation range but shallower, flatter reflectors showing non-physical, negative attenuation. Alternatively, in cases where reflectivity is expected to be nearly constant between reflectors, many internal reflectors should be used to calculate attenuation as in Methods 5 and 6. At SPL, these multiple-reflector methods showed the greatest fidelity in constraining attenuation, especially for the case of a depth-resolved attenuation profile through the ice column. Regression uncertainty is strongly dependent on the degrees of freedom, and when many thousands of measurements are used in the regression, the reported uncertainty can greatly underestimate the physical uncertainty in the system. In summary, care should be taken to select the method which most effectively minimizes the sources of uncertainty and therefore isolates the depth–power relationship in Eqn (4).