A 3D map of englacial attenuation rate from radar reflections at Law Dome, East Antarctica

The East Antarctic Ice Sheet (EAIS) is the largest source of potential sea-level rise, containing approximately 52 m of sea level equivalent. To constrain estimates of sea level rise into the future requires knowledge of ice-sheet properties and geometry and ice penetrating radar offers a means to estimate these properties (e.g. ice thickness, englacial temperatures). One of the regions that has been extensively surveyed using ice penetrating radar from the Investigating the Cryospheric Evolution of the Central Antarctic Plate (ICECAP) project in East Antarctica is Law Dome, a small independent ice cap situated to the 5 west of Totten Ice Shelf. The ice cap is slow-moving, has a low melt-rate at the surface and moderate wind speeds, making it a useful study site for estimating the radar attenuation. A new method is proposed for the estimation of attenuation rate from radar data which is mathematically modelled as a constrained regularised l2 minimisation problem. In the proposed method, only radar data is required and the englacial reflectors are automatically detected from the radar data itself. To validate our methodology, attenuation differences at flight crossover points are calculated and statistical analyses performed to assess the 10 accuracy of the results. For spatial analyses, the errors are of the order 22.6%, 15.2%, and 32.8% for mean absolute error, median absolute error, and root mean square error respectively. Also, for the depth analyses, up to the depth of 800m, the errors are under 29.9%, 24.2%, and 38.8% for mean absolute error, median absolute error, and root mean square error respectively. A final product of 3D attenuation rates and uncertainty estimates is provided. The generated dataset is publicly available at https://doi.org/10.25959/5e6851e266f4a (Abdul Salam, 2020). 15

Previously, many studies assumed englacial attenuation rates to be uniform within the ice and proportional to the ice thickness (Bentley et al., 1998;Rippin et al., 2004;Jacobel et al., 2010). However, the attenuation rate varies with location, due to changes in ice properties (Matsuoka et al., 2011;MacGregor et al., 2012). For example, a study in central West Antarctica shows that the one-way depth-averaged attenuation rate varies horizontally by 5 dB km −1 in absolute terms along radar transects (Matsuoka et al., 2010). Matsuoka et al. (2011) argued that the assumption of regionally uniform attenuation rate fails 40 in most cases due to varying ice chemistry and temperature which leads to false attenuation estimates. However, these studies explored attenuation rates only in small regions and the understanding across broader spatial scales needs to be improved. MacGregor et al. (2015) demonstrated the use of radar reflections from englacial layers to constrain the attenuation rates and temperatures as a function of depth in the Greenland ice-sheet. A key advantage of the method derived by MacGregor et al. (2015) is that it does not rely on echoes from the ice-sheet bed to determine attenuation, as the bed is complex and spatially 45 variable (Winebrenner et al., 2003), complicating interpretation of radar data near the bed. The disadvantage of this method is that it requires radiostratigraphy (study of layering by means of radar reflections) data to reliably trace the englacial layers (MacGregor et al., 2015). Many studies considered inferring the englacial attenuation or temperature from airborne radar data (Robin et al., 1969;Barrella et al., 2011;MacGregor et al., 2012MacGregor et al., , 2015Schroeder et al., 2016), but none have done it in a way that only relied on radar data.

50
Another radar-based method to constrain attenuation rates was proposed by MacGregor et al. (2007). In this method it assumed constant reflectivity values for the internal layers and the method also exploits the dependence of radar attenuation on ice temperature. However, this approach requires undisturbed englacial layers (i.e., having clear boundaries in englacial layers), which cannot be applied to many important regions because of the lack of englacial layer information. This depthaveraged method is applicable at the ice-sheet scale, which requires contiguous and undisturbed englacial layers, and cannot be 55 applied to areas which include highly-crevassed and fast-flowing regions near grounding zones and shear margins (MacGregor 2 https://doi. org/10.5194/essd-2020-146 Open Access Earth System Science Data Discussions Preprint. Discussion started: 27 October 2020 c Author(s) 2020. CC BY 4.0 License. et al., 2015). However, as these regions often influence the dynamics of the ice-sheet to a greater extent layer-based approaches cannot be used to obtain englacial temperatures in many critical areas for ice-sheet modelling.
As an alternative to the methods described above, Schroeder et al. (2016) developed an adaptive approach for estimating englacial attenuation rates for the entire ice column in Thwaites Glacier, West Antarctica. In this method, the unfocused radar 60 bed echoes are fitted based on the correlation of ice thickness and the corrected bed power echo. The method performs weakly across the catchment near steeply sloping bed topography. In addition, another concern is that this method returns a 2D spatial map; however, deriving englacial temperatures requires estimates of attenuation throughout the ice column, requiring a 3D attenuation map. The temperature gradient can then be used to constrain the geothermal heat flux in the ice-sheet.
The key motivation for this work is to enable the mapping to englacial temperature profiles possible, which can help in 65 estimation of geothermal heat flux. Accurate estimates of englacial temperature and geothermal heat flux are incredibly important for constraining model simulations of ice dynamics (e.g. viscosity is temperature-dependent) and sliding. However, we currently have few direct measurements of vertical temperature (i.e. only at boreholes/ice domes) and geothermal heat flux in Antarctica. This method derives attenuation rates, that can then be mapped directly to englacial temperatures and geothermal heat flux. Several earlier studies have derived Arrhenius based functions to relate attenuation and ice temperatures in Green-70 land (MacGregor et al., 2015;Jordan et al., 2016) and West Antarctica (MacGregor et al., 2007). A mapping function can also be derived using a temperature profile as output and attenuation rates (for a certain radius around the borehole) as input. A temperature profile from the ice-borehole at Dome Summit South is available (Dahl-Jensen et al., 1999) which can help in generating mapping function. The attenuation rates derived here can then be transformed to temperature profiles by utilising the mapping function. Geothermal heat flux can be estimated from the gradient of these vertical temperature profiles (working 75 paper).
Our proposed method overcomes some of the limitations of the previous methods. Initially, the englacial reflectors are detected in the reflection data and then the mathematical technique that minimises the mismatch is used to obtain the estimates of englacial attenuation. This method does not require any additional data-sets (such as englacial stratigraphy). Therefore, this method can be applied to regions which only have radar reflection data. A large amount of radar data has already been acquired 80 for a number of purposes across vast areas of the ice-sheet. In East Antarctica, stratigraphy is not available except for limited regions typically around ice coring sites (e.g. at Dome C, where some internal layers are traced (Cavitte et al., 2018)), and while initiatives such as AntAchitecture aim to address this issue in the longer term, the majority of radar data is currently without stratigraphy. This method can utilise the existing acquired data and will extract the attenuation rates, which can further be used to estimate temperature profiles; consequently, it can improve the geothermal heat flux maps and ultimately can help 85 improve ice-sheet modelling predictions. The 3D map of englacial attenuation rate is not available from any previous study and it can be a good data set for testing different novel algorithms. Here, the radar data from the Investigating the Cryospheric Evolution of the Central Antarctic Plate (ICECAP)  mission for Law Dome East Antarctica are used. In the next section, the importance of the Law Dome region is discussed. Section 3 describes the methodology used for estimating the englacial attenuation rates. An assessment of the quality of the results, crossover and uncertainty analyses are presented in   Roberts et al., 2015) which can be helpful in mapping attenuation to temperature profiles; and 3) the ice on Law Dome is slow-moving and stable, and there is negligible melting at the ice surface (Etheridge et al., 1998) and relatively moderate wind 100 speed (8.3 ms −1 ) (Morgan et al., 1997), which means that the englacial layers are relatively undisturbed.
The ICECAP geophysical data that we use in this study were collected over the period 2008-2012, covering over 14800 linekm of the Law Dome region Young et al., 2011). The survey aircraft was fitted with a High Capability Radar Sounder (HiCARS) instrument, with a central frequency of 60 MHz and bandwidth of 15 MHz. In order to retain full energy of the radar reflections, the radar data was processed after pulse compression (Wright et al., 2012). Young et al. (2011) 105 determined the ice thickness from ICECAP airborne radar data using an ice speed of 169 m/µs for electromagnetic wave propagation and no firn correction was applied for the reduced velocity of radio-waves. In this work the bed-echo and aircraft height data, ice thickness, and radar reflections will be used in order to extract the attenuation rate. Figure 1 shows the coverage of these data sets from ICECAP project for the entire Law Dome region.
The spatial resolution of the ICECAP radar data is approximately 20 m along the flight line while the spacing between the 110 flight lines is in order of kilometres and is not evenly sampled as shown in Fig. 1.

Methodology
We use radar reflectors from within the ice column to estimate radar attenuation. Radar attenuation is the loss (dB) of signal strength due to dielectric absorption, scattering and geometrical spreading. The attenuation rate given by following Equation 1 describes the loss of energy from source to receiver per unit distance (Reynolds, 2011).
The attenuation rate (N ) strongly depends on electrical conductivity (σ) , permittivity ( ) and angular frequency (ω). As some of these properties are themselves strongly dependent on temperature, variations in the attenuation coefficient can be used to estimate englacial temperature, subject to making assumptions, such as constant conductivity (same chemical properties).
Radar attenuation is proportional to conductivity and depends on the ice chemistry/acidity but it varies slowly spatially. In the 120 case of Law Dome, the concentrations of impurities in the ice column are low (Etheridge et al., 1998).  In our analysis of radar reflectors from within the ice column the first crucial step is to locate the englacial layers with an assumption of having the same reflection strength, and all the reflectors are assumed to be specular (mirror-like) (Schroeder  et al., 2015). Once the englacial layers are detected down the ice column, the located depth and strength of radar can be used to infer the attenuation rate. The resulting attenuation rate can be validated by using crossovers of the flight lines to check the 125 robustness of the method and estimate the uncertainty in the obtained attenuation. After obtaining the attenuation profiles, the data is post-processed for visualisation.

Locating reflectors
In locating the reflectors, we process each vertical ice column of radar data (also known as a trace) independently to produce high resolution spatial data at approximately every 20 m along the flight lines. For a column with M reflectors (M −1 englacial 130 layers), where the index i shows the i th reflector from the surface, i.e., i = 0 corresponds to the ice surface and i = M − 1 corresponds to the lowest englacial reflector, the general equation for specular internal reflections from the i th reflector is given by (MacGregor et al., 2009(MacGregor et al., , 2015: Here, P t is the transmitted power, λ air is radar signal wavelength in the air, G a is the antenna gain, T is the transmission 135 loss at the interface of air and ice, L a is the one-way attenuation loss, L vs is the total loss due to volume scattering, L b is the loss due to birefringence, L sys is the total system loss, G p is the processing gain, h is the aircraft height above the ice surface, d is depth of the ice and is real permittivity of the ice. In Equation 2, P t , λ air , G a , T, L sys , G p , are assumed invariant for any given trace of radar data, which remains the same vertically down the trace, but can vary horizontally along the flight line. Similar to Matsuoka et al. (2010), we assume that L vs and L b are negligible. The received signal power P i r from each reflector 140 needs to be corrected first for the geometrical spreading so that the attenuation due to spreading of radiowave travelling down the ice column is compensated and the main attenuation factor should be only dependent on the ice properties. In Equation 3, the received raw radar signal is added with the geometrical correction factor (h distance travelled in the air and d within the ice), and then converted to decibels ([x] dB = 10log 10 x) for easier calculations. The radar sounding geometrically-corrected power is given by (MacGregor et al., 2015;Schroeder et al., 2016): where [P i r ] dB is the received power in decibels. After the power correction, the next step is to locate the englacial reflectors. Englacial layers have unique electrical properties, realised as peaks in englacial reflectivity. Sometimes there are false reflectors which means that two consecutive peaks may occur very close together, implying that the englacial layers are very thin. To avoid this situation, an assumption is made that 150 englacial layers are no thinner than 60 m. These potential false reflectors make the process unstable and sometimes result in too many reflectors if the condition of minimum thinner layer is ignored. The value less than 60 m fails in certain regions and this value is minimum value which worked for the whole region. Any reflection peak occurs closer than 60 m is considered as one englacial layer which makes our method capable of capturing vertical resolution of 60 m. In an example profile in the detected reflectors are shown by small red circles. It can be seen that out of 14 reflectors, three reflectors are discarded 155 because the attenuation rate is always positive and the signal is always attenuating over the depth. In other words, for constant strength reflectors, the resulting reflections will always lead to decreasing strength. However, the deeper layers can have higher reflectivity if they contain higher conductivity materials. The assumption of uniform layer reflectivity is a simplistic approach that we adopt at this current stage but lack of conductivity profiles at Law Dome. It will be worthwhile in the future to look into other approaches if englacial reflectivity data become available.

Retrieval of Attenuation
After detecting the englacial layers, the thicknesses of each englacial layer are calculated. Now, the attenuation rate estimation problem is formulated as a Ridge regularised minimisation problem (Warton, 2008). In this problem, the required components are englacial layer thicknesses and strength of power from reflectors, while the attenuation rates are calculated. A penalty term is added to the minimisation problem to suppress unphysical high frequency oscillations due to overfitting but which provides the problem becomes a minimisation of the least squares solution which means relying entirely on the data, but the solution may contain physically unrealistic oscillations. The given minimisation problem is as follows: where Z is the englacial layer thickness square matrix (M − 1 x M − 1), which contains the information of the thicknesses of all the englacial layers, ∆z i is thickness of i th englacial layer, N is the radar signal two-way attenuation rate vector (M − 1 x 1) within englacial layers, P is the power differences (in dBs) vector (M − 1 x 1) of the strongest reflection with each i th layer and λ is the regularisation parameter to tune the model. Here, Z resulted in a matrix because we formulated for englacial layers attenuation (i.e. exactly within the layer) rather than total attenuation rate (MacGregor et al., 2015). As we are resulting attenuation must be a non-negative value. Our formulated problem is not a simple Ridge regression, it has additional constraint of N > 0 which will not behave exactly simple Ridge regression, but the general behaviour is of Ridge regression.
The optimisation problem of Equation 4 is shown graphically in Fig. 3, which shows how it is different from least square 180 solutions and conventional Ridge regression.
In Ridge regression (Hoerl and Kennard, 1970), selection of λ is main parameter. A higher value of λ will force the solution to be smoother but at the cost of diminishing the effect of data on the modelled solution. In addition, higher values of λ will also suppress attenuation towards zero (not exactly zero like lasso regression (Tibshirani, 1996)). With higher value of λ, the optimisation cost function rely more on λ N than the traditional least square error component. Our aim here is to find a 185 value of λ which leads to a more robust and balanced solution. We explored number of values for λ ranging from λ = 0 to λ = 10. When the value exceeds λ = 0.85, the solution is physically unrealistic (extremely small values for the attenuation rate  which is using a similar survey system of ours (HiCARS with frequencies ranging from 52.5 M hz to 67.5 M Hz). The reported value of attenuation rates helped us to determine acceptable range of attenuation rates for a similar system. Attenuation rates varies with frequencies and survey systems, so it is important to consider radar survey system when comparing typical values of attenuation rates. The reported attenuation rates for the same system at Thwaites Glacier (Schroeder et al., 2016) helps us in 200 determining a physically realistic attenuation rates at Law Dome. As we change the value of λ, the value of attenuation rate is also affected, because the amount of bias in solution depends on λ in Ridge regression. The higher the value of λ, the values of attenuation rates are pushed more towards zero. So, the reported attenuation rate in Schroeder et al. (2016) helped us to constrain λ, which leads to acceptable physically realistic range of attenuation rates.
The primary disadvantage of this method is that numerous englacial layers must be correctly detected to reliably constrain 205 the attenuation rate. Furthermore, in order to avoid dependence on false reflectors, several conditions are required to be met for processing the ice column, i.e., the ice thickness must be at least 200 m in depth and the number of internal reflectors must be greater than or equal to four.

Data quality
In this section, the compilation of 3D attenuation data from the processed traces is discussed in detail. As each trace is processed 210 independently, the resultant attenuation rates are not well-distributed spatially over Law Dome. The reasons for this are: firstly, the flights do not follow uniformly spaced gridded paths; and secondly, the proposed method sometimes discards the traces with less than four internal reflectors or the ice thickness at that location is less than 200 m. For this reason, a 3D grid is generated which can help in visualising and binning the attenuation rates to the corresponding 3D cell.
Next, the results are further processed to calculate 3D distribution of the attenuation rates. In the 3D distribution, the 2D 215 spatial coverage is 1 km x 1 km bins while each bin point is a trace down the ice column depth, at a spacing of 60 m. Each 3D cell point may be influenced by many independent processed radar traces. Many radar trace values may fall within a 1 km x 1 km bin, be a single value in a bin or have no radar data. If there are many values falling in that 3D cell, the median is selected.
The median values are selected for filling the grid because it performs better due to the reduced number and magnitude of outliers in the data. If no points fall within that cell, it will result in missing values (given as NaNs).

Depth-averaged attenuation and uncertainties
In this section, the processed data is visualised in terms of depth-averaged attenuation rate along with the uncertainties and the number of points spatially contributing to the plot. In the 3D grid, the spatial map over Law Dome is obtained by weighted average over all the depths in the traces. In the post processing of the attenuation results the following steps are done: 1. A 3D grid is created where each cell represents 1 km x 1 km x 60 m.

225
2. For each flight line, each trace is mapped to the closest cell in the 3D grid. We also have an additional 3D grid which retains information about how many values are contributing to the corresponding cell. Another spatial 2D version is added which represents how many traces from the flight lines are contributing to each bin as shown in Fig. 6.
3. The depth-averaged attenuation rate is calculated as shown in Fig. 5(a).
4. Standard deviation within each bin is calculated in order to find the uncertainty in results as shown in Fig. 5(b). Figure 5(a) shows the estimated englacial attenuation rate and the corresponding standard deviation in each bin is shown in Fig. 5(b). As the flight spacing is non-uniform, there are different numbers of points contributing spatially, as indicated in Fig. 6.

Crossover Analysis
In this section, the attenuation rate results are analysed at the crossover points in order to see the robustness of the method. A 235 crossover is a location where two or more flight lines cross each other (here crossover is considered within 35 m). There are 582 crossovers. Here, we consider results averaged both over the depth and spatially.
In the first spatial analysis, three errors are computed: mean absolute deviation (MAD), median absolute deviation (MEDAD) and root mean square error (RMSE). In MAD, the deviation is greatly affected by outliers; in contrast, MEDAD values are much less sensitive to outliers. In the second analysis, again the same three crossover differences of MAD, MEDAD and RMSE are calculated at the crossover as a function of depth. Here, the attenuation differences are seen down the depth which shows that the differences are 250 low at the top englacial layers and high at the bottom. This suggests that the deeper the englacier layer, the harder it will be to detect. Figure 8 shows the attenuation rate over the depth for all the crossovers. The number of points contributing to the error calculations is also shown as a function of depth. The error at the ice base is high because detecting deep reflectors are difficult, the number of samples at these depths for crossovers are limited, and the absolute values of attenuation rates increases with depth. In terms of percentage depth errors the error increases with depth. Up to the depth of 800 m, the errors are under 29.9% 255 for MAD, 24.2% for MEDAD, and 38.8% for RMSE relative to the absolute values of attenuation rates. As we go deeper the number of samples are also reduced and the deeper reflectors are hard to detect correctly so the error increases. Because of less number of samples at deeper depths, the results uncertainty increases.
As discussed above, the minimum MSE occurs for the values range 0.3 ≥ λ ≤ 0.4, which is shown in Fig. 4. We further explored which values will provide a more robust solution in context of crossover differences. The normalised median absolute 260 deviation is plotted for all the crossovers as shown in Fig. 9 with different values of λ. In addition, these errors are averaged and is shown in Fig. 10. Both figures show the solution is more robust at λ = 0.3 and the attenuation rates result in values closer to typical reported values (Schroeder et al., 2016). There are several reasons for the crossover differences. The crossovers can be from different field seasons and hence the physical situation of the ice-sheet may differ and also small differences in equipment, result in varying radio echos. Another 265 important reason is that ice is an anisotropic material (Thorsteinsson, 2001) and the acquisition of radar data is from different flight path directions. For majority of crossovers, the error values are low in comparison with the estimated attenuation rate values which implies the robustness of our proposed method. Also, to each cell in the 3D grid, there are many samples contributing to each cell which further reduces the uncertainty in the estimated attenuation rates.

270
In this paper, a new method for attenuation rate estimation is developed and used at the Law Dome region. The importance of this method is that it only relies on radar data and can be applied to regions where there is no information about the stratigraphy of englacial layers. We show that our method can be applied where only radar data is acquired. This is an important advancement because from attenuation rates, the temperature profiles can be mapped, which could provide further insight into heat sources at the base of the ice-sheet (e.g. constraints on geothermal heat flux).

275
A high resolution attenuation rate can be obtained by this method with acceptable uncertainty and errors; however, it can be further improved by treating the assumptions made. Firstly, the quality of the estimated attenuation rates can be further improved if the actual strengths of the reflectors are used rather than assuming constant strength reflectors. This is only possible when we have reflectivity of englacial layers or other proxy to get this information. Secondly, the ice anisotropy should be taken into account which is responsible for the different behaviour of radio waves when incident on the ice from different directions.

280
Finally, the automatic detection of internal reflectors can be further improved if the reflectors detected are traced in the adjacent ice columns or using prior information about the englacial reflectors using stratigraphy. However, in reality, the availability of such reflector information is very rare compared to the radar coverage. Northing (  Northing ( (Schroeder et al., 2016). As the values of λ increases, the normalised deviation for crossover increases, i.e., c) λ = 0.4 and λ = 0.6.
Also note that the higher values of λ will force the resulting attenuation rates closer to zero and hence will be smaller than the typical values of attenuation rates.