RADAR-DERIVED INTERNAL LAYERING AND BASAL ROUGHNESS CHARACTERIZATION ALONG A TRAVERSE FROM ZHONGSHAN STATION TO DOME A, EAST ANTARCTICA

The internal layers of ice sheets from ice-penetrating radar (IPR) investigation preserve critical information about the englacial conditions and ice-flow field. This paper presents a new detailed analysis of the spatial distribution characteristics of internal layers and subglacial topography of East Antarctic ice sheet (EAIS) from Zhongshan station to Dome A. Taking the internal layering continuity index (ILCI) and basal roughness as indicators, it provides an opportunity to evaluate the past internal environment and dynamics of ice sheet. The radar data of 1244 km along a traverse between Zhongshan Station and Dome A of EAIS was collected during the 29th Chinese National Antarctic Research Expedition (CHINARE 29, 2012/2013). Except for the upstream of Lambert Glacier, the patterns of the folds in the internal layers are basically similar to the bed topography. The relatively flat basal topography and the decrease of ILCI with the deepening of the depth provide evidence for identifying previous rapid ice flow areas that the satellite cannot obtain, especially in the upstream of Lambert Glacier. Well continuous internal layers of Dome A almost extend to the bed, with high ILCI and high roughness characteristics. There are three kinds of basal roughness patterns in the whole traverse. The characteristics of the internal layer and basal topography of the traverse between Zhongshan Station and Dome A provide new information for understanding the ancient ice-flow activity and the historical evolution of EAIS. * Corresponding author


INTRODUCTION
The East Antarctic ice sheet (EAIS) is the world's largest ice body, and any minor changes can cause global sea-level changes. The latest research shows that in warmer climates, EAIS may be much more vulnerable than expected, and its contribution to predicting future global sea level may also be underestimated (Aitken et al., 2016). However, the lack of understanding of the structure and long-term evolution of EAIS limits our assessment of its potential changes. Quantifying the internal information and basal characteristics of the ice sheet detected by ice radar is an effective way to achieve this goal (Sime et al., 2014). Ice-penetrating radar (IPR) or radio-echo sounding (RES) systems operate by emitting electromagnetic waves through the air and the ice sheet, which will be reflected at boundaries with different dielectric parameters. The return signal records the formation of the ice sheet. Therefore, the detected internal stratigraphy and basal topography can provide information about the age of the englacial structure ( Winter et al., 2019), ice history (Rippin et al., 2003), and ice conditions (Wilkens et al., 2015), can also provide a basis for building ice sheet models (Young et al., 2017) , which have an irreplaceable role in assessing ice sheet mass balance and ice dynamics. Internal layers in the ice sheet were initially formed by the surface accumulation and were later modified by subsequent burial, compaction and ice flow to form a layered structure. Usually, internal layers geometry is affected by changes in substrate topography, basal conditions or accumulation rates, and ice flow effects . Previous scholars used internal layers to study the relationship between the internal structure of ice sheet and ice flow velocity (Karlsson et al., 2009). Internal layers are divided into two types: one is a continuous and well-defined internal layer which primarily follows surface and/or bed topography and another is a discontinuous layer which is buckled or disrupted layering. The most important reason for the discontinuous layering may be affected by previous or current ice flow at different velocity (Karlsson et al., 2012;Winter et al., 2015). Ice flow behavior and evolution are mainly controlled by subglacial bed conditions (Vaughan et al., 2006). Ice flow depends also strongly on bed conditions (thawed or frozen) and ice rheology (Bell et al., 2011). A significant indicator of bed conditions is basal roughness, which can be used to infer ice dynamics and subglacial geomorphology interpretation . Basal roughness reflects the changing ice dynamic and thermal regimes with history superimposed over of ice sheet, which in turn modify the preglacial topography and the underlying lithology and geological structure (Siegert et al., 2005). The use of ice velocity in combination with internal layers and bed topography to study ice sheet stability and ice dynamics has been widely used, including Greenland (Sime et al., 2014), Institute and Möller Ice streams (Winter et al., 2015), Pine Island Glacier (Karlsson et al., 2012;Wilkens et al., 2015), the Siple Coast ice streams in West Antarctica and the Wilkes Land, Wilkes Subglacial Basin, Aurora Subglacial Basin region (Siegert et al., 2005;Wright et al., 2012) and Dronning Maud Land (Fujita et al. , 2012) in East Antarctica. However, there is still a lack of specific research on the distribution characteristics of the internal layers of Princess Elizabeth Land (PEL) in East Antarctica.
As an important transect of the ITASE (International Trans-Antarctic Scientific Expedition) project, the inland traverse of the CHINARE (Chinese National Antarctic Research Expedition) from the coastal Zhongshan Station to Dome A in East Antarctica has been monitored continuously for 20 years. The traverse, passing through the western side of PEL, eastern upstream of the Lambert Glacier, Gamburtsev Subglacial Mountains (GSM) and Dome A region ( Figure 1) . There are ice thickness, subglacial topography and internal layers by IPR surveys (Tang et al.,2016), surface mass balance (SMB)/snow accumulation by stake measurements and atmospheric observation with automatic weather stations (AWS) (Ding et al.,2015), surface topography and velocity with GPS (Yang et al., 2018), which provide a reliable dataset to analyze age stratigraphy, ice dynamics and mass balance of ice sheet. In this paper, we have a new detailed analysis and interpretation of the subglacial bed topography and the internal layers of EAIS from Zhongshan Station to Dome A. We use the vehicle-based deep IPR data collected during CHINARE29 (2012/2103). The internal structure of the radar profile and bed topography are quantitatively analyzed by internal layering continuous index (ILCI) and basal roughness index, which are applied to initially interpret and assess the ice flow behavior and ice dynamics.

IPR data
In order to accurately measure the ice thickness, high-resolution internal structure, and subglacial topography, a vehicle-based deep IPR was used along the 1244km traverse ( Figure 1) from Zhongshan (76.281°E ， 69.574°S) to Dome A (77.133°E ， 80.409°S) in EAIS during the CHINARE 29 (2012/2013). Zhongshan here refers to the starting point of the surveying line, which is 22.7km away from the Zhongshan Station. The IPR, also called High-Resolution Ice-Sounding Radar (HRISR), is designed to operate as a linear frequency modulated (LFM) pulse system by the Institute of Electronics, Chinese Academy of Science, with a center frequency of 150MHz and a bandwidth of 100MHz. The antenna system has two logperiodic antennas on a specialized versatile mobile observation cabin (SVMOC), which is 2m above the snow surface. The main technical parameters of deep IPR are given in Table 1. The data acquisition system digitizes the return signal by 12-bit A/D converter at a sampling frequency of 500MHz and simultaneously records GPS signal. The GPS horizontal error is <7m, and the vertical error is <10m .  (Fretwell, et al.,2013). The red stars represent the location of Zhongshan station, Taishan station and Kunlun station, respectively. The white point marks the location of Dome A.
IPR data is simply processed to get a clearer radar profile. Firstly, the first 2000 samples of the raw data are eliminated, which are the blind area caused by the receiver-blanking switch for avoiding receiver saturation. Secondly, pulse compression and coherent stacking 60 traces are carried out. Then, DC filters and mean filters processing is carried out to reduce blurring noise. Finally, conversion of time to depth is achieved by h=v*t/2, where v is the electromagnetic wave velocity in the ice sheet, and t is the time delay between the air-ice interface and the target. Here, we select a speed equal to 0.168m ns -1 , which produces an ice thickness error of less than 20m . In the field, ice radar detection was not performed before the snowmobile was still running smoothly, which made the acquired radar data geographically discontinuous. Matching each trace of the radar data with the GPS latitude and longitude recorded at the same time, the calculated gap of radar trace is 4m with a relative error of 5cm. The geographic location of the radar line as shown in Figure 1 is obtained. Based on the surface and subglacial characteristics, we divide the transect from Zhongshan to Dome A into five zones for IPR profiles analysis. Respectively, zone A (0.0-262.0km away from Zhongshan), whose surface elevation is rapidly rising, has an average slope along the transect of 6.3m km -1 . Zone B (207.2-431.7km) is a slightly rising area with an average slope of 1.8m km -1 . Zone C (461.7-752.1 km) is a relatively flat area with an average slope of 0.76m km-1. Zone D (752.1-972.2 km) is an upstream area of the Lambert Glacier and zone E (972.2-1243.7km) is the Dome A area.

Internal layer continuity index
The internal layer of the ice sheet not only records the stratigraphy history, but also holds the key information of changing the ideal fallout state of ice and snow by the daily wind, ice flow and base melting. Continuous internal layers indicate that these processes are weak or absent. Previous studies have used internal layers to qualitatively interpret the relationship between the internal structure of the ice sheet and ice-flow velocity (Karlsson et al., 2009;Rippin et al., 2003). Recently, the continuity index (CI) method proposed by Karlsson et al. (2012) has been used to quantitatively characterize the apparently continuous and fractured internal layer regions from large-volume IPR datasets, and has been widely applied to assess ice-flow history and ice dynamics Sime et al., 2014). The internal layer continuity index (ILCI) is a calculated value of each A-scope that increases with the number and strength of the internal layer detected. Discontinuities in the number or strength of the internal layer make it easy to identify spatial changes detected. Depending on the rapidity of the signal change between extremely low and high values, the gradient in each A-scope trace can be used to quantify the continuity of internal layers, which is an internal layer continuous index (ILCI) (Karlsson et al., 2012). Specifically, we calculate the mean of the absolute value of the derivative of reflected relative power in a time window sub-interval [L1: L2] between the ice surface and the basal interface, i.e.
where Pi is the reflected relative power (dB) at each time sample i and is the depth (m). If the total number of time samples between the surface and basal is M, we choose a subinterval N as [(1+M/5) : (M-M/5)] =[ L1: L2] to avoid surface noise and the reduced signal from the deepest part of the ice. In essence, a high continuity index corresponds to clear internal layers, while sections with discontinuous or absent layers return relative low continuity index. We use the middle three fifths of the ice column for analysis, and average over windows containing 500 traces (2km) to assess overall layer continuity along the traverse. Smoothing over 500 traces is useful for identifying anomalies, while the smoothing over 2500 traces gives a clearer regional picture, hence both plots are included here.

Basal roughness
Generally, basal roughness is defined as the degree of the vertical variation of topography within a given horizontal distance, which can be used to quantitatively describe the irregularity in a basal interface. The basal roughness of ice sheet as an indicator of subglacial conditions and as a potential control on ice-flow has been applied to studies of ice dynamics and subglacial geomorphological interpretation Siegert et al.,2005;Wright et al., 2012;Wilkens et al., 2015). In most of the recent glaciological studies, a single-parameter roughness produced by Fast Fourier Transform (FFT) analysis of bed topography is defined as the integral of the spectrum within a specified wavelength range Siegert et al.,2005). This parameter as total roughness ( ξ t) represents the amplitude of the undulations, which is used to quantify the vertical variation of the bedrock, but the information about the horizontal frequency is lost (Wilkens et al., 2015). Li et al. (2010) introduced a second parameter as frequency roughness (η) to explore the frequency variation of fluctuations. This measure is calculated as the total roughness index (ξt) divided by the bed slope index (ξsl). The parameter (η) corresponds to the horizontal irregularity of the section and invariant of the vertical stretching of the profile .
Here, we use the same method of a two-parameter roughness index (ξt and η) to quantitatively describe the basal roughness of both vertical and horizontal irregularities in the subglacial basal interface. The basal roughness is calculated starting in a profile of length, l. The length (l) is a moving window with 2 N sampling points, where N=5 (32 sampling points =128m), which is the suggested minimum value (Taylor et al., 2004). The length of the moving window is also in the meter-scale range to explore the smallest-scale roughness patterns possibly. First, the spectral power density, S(k), is calculated by FFT on the bed elevation (Z(x)) of which mean bed elevation (<Z(x)>) is subtracted. The total roughness (ξt), indicating the amplitude of the elevation fluctuation, can be obtained by integrating S(k).

Where Z0(x) is the detrended elevation profile, Z0(x)= Z(x)-<Z(x)>. x is the horizontal location. Z(x)is bed elevation. <Z(x)> is the mean bed elevation of a profile of length, l.
Then, the slope of the detrended bed profile, is used to obtain a frequency roughness of the second parameter, η, indicating the horizontal irregularity of the bed. That is, Finally, if the integral is between 0 and ∞ , the total roughness ( ξ t) is equal to half the mean square of the bed elevation profile with units of m 2 Wright et al., 2012). So, we use the same basal roughness indexes ( 2 t  and 2 ), which were adopted by Wright et al. (2012) in their analysis, to quantitatively analysis the basal roughness of subglacial bed topography.

RESULTS AND ANALYSIS
Figure 2 shows a regional view of the IPR profile along the traverse and of the overall internal layering properties calculated using the ILCI method. We have stitched the profile of the subregion, where the breaks indicate that the profile interval is more than 5 km. The color bar above the profile represents ILCI of the middle three fifths of the full ice column, calculated from 500 trace moving windows (~2km). The result of internal layering properties by ILCI can also be qualitatively seen from the radar profile.
ILCI can quantitatively and efficiently analyze the continuity of internal layers of large-volume IPR datasets. As shown in Fig. 3, areas with clear and consecutive internal layers return a high ILCI; while areas with disrupt or absent internal layers return a low ILCI. When comparing the results of ILCI directly between different areas, it should be ensured that the data is completely acquired and processed with the same parameters. Here, the overall IPR data was processed in the same way. However, the radar profiles at 71-185km and 267-320km show obvious quality differences. These regional data are interfered with by the field acquisition, which makes the radar data quality decline, and the ILCI value also decreases. In addition, Zone A, Zone B (except 400km) and Zone D (except 880km) also have relatively low ILCI, where the internal layers are discontinuous or disrupted. In particular, Zone D is the upstream Lambert Glacier region characterized by a wide valley with deep depression at the bottom (Fig. 3d), where most of the bedrock is not detected. It may be caused by the scattering and absorption of radar signals in the ice. The ice thickness in this area is more than 3000m, and the undulation amplitude of the internal layers is small. It is speculated that the bed undulation range may be small and the base groove may be shallow. The parallel and continuous internal layers near the surface are located on the deep disrupted layers, and most of the deep internal layer structures resemble "whirlwinds" or "tornados" (Karlsson et al., 2009), such as the white rectangle of T area in Figure 2d. The structure in T is a complex response of the radar signal to the high-amplitude buckling of layers within the ice sheet. Therefore, the whirlwind can be used to identify and map the areas with buckled and disrupted layering (Winter et al., 2015). Except for the Zone D area, the internal layers along the traverse are similar to the undulation of the bed topography. The closer to the internal layer of the surface, the weaker the undulation. Zone C has a median ILCI with clear continuous internal layers distributed parallel to the bed and discontinuous or absent internal layers in the deep ice. Zone E is the area near Dome A in the inland center of the Antarctic ice sheet, with a relatively high ILCI and significant fluctuations ( Figure  2e, Figure 3a). In Zone E, the clear internal layers extend almost to the bed, with the largest proportion of clear internal layer thickness, exceeding 90% of the ice thickness. However, Zone C and Zone E have low ILCI at local locations (such as near 735km, 740km, 1010km, 1100km, 1166km and 1207km), which may be the result of destructive interference caused by stacking adjacent traces with asynchronous reflected signal at steep peaks and deep valleys.
From the distribution of the bed elevation in Fig. 3a, it can be seen that there is a variety of mountainous topography in PEL. The subglacial topography is mainly mountainous and valley glaciers. Within 0-100 km near the coast, the bed undulation is relatively flat and the bed elevation is near sea level. The average bed elevation between 120km and 390km is 380m and the fluctuation range is obvious. At 390-420km, the average elevation of bedrock is about 1000km, and the frequency of undulation is faster. From 420km to 750km near east upstream of Lambert Glacier, the bed is undulating with many deep valleys and high peaks. The entire of the bed topography in upstream of the Lambert Glacier was not detected due to insufficient ice-penetrating capability of IPR, and the basal roughness could not be accurately assessed. However, there are no significant mutations of amplitude in the distribution of the internal layers, and the overall distribution is similar. Based on the detected partial bed topography, it is speculated that the area may have similar topographical features. Therefore, we assume that the bed topography depicted in the upstream of the Lambert Glacier is established. There is the largest deep valley about 165 km wide, and the bed elevation may be below sea level. Due to the existence of the GSM, there are deep valleys and peaks at Dome A, where the average bed elevation is above 2000m and the changes relatively frequently and sharply. At the edge of the Lambert Glacier and GSM, the relative bed elevation difference is approximately 1200 m. Due to the significant uplift of the bedrock in the GSM region and the thinning of the ice sheet, change of basal thermal mechanisms has caused an important impact on bedrock roughness .   3b is shown a basal roughness distribution except for the area where the bed is not detected in upstream of Lambert Glacier and where the radar is not detected. To better characterize the region, we performed a smoothing average of 500 (2 km) and 2500 (10 km) traces for the total roughness ξ t and frequency roughness η (Fig. 3b). There is a high total roughness in the 390-410km and Zone E regions, and the others have a low total roughness value along the traverse; while most of the frequency roughness has a relatively high value, only in 220-390km and Zone D have low frequency roughness. Thus, the basal roughness along the profile can be divided into three types, which are: low ξ t low η area (220-390km and Zone D), low ξt high η area (0-220km and 410-750km), and high ξt high η area (390-410km and Zone E).

CONCLUSIONS
We use the IPR data along the traverse between Zhongshan Station and Dome A of EAIS, to analyze in detail the internal stratigraphy and the basal roughness in the PEL. Except for the upstream of the Lambert Glacier, the large number of folds of the internal layers in the remaining areas are mainly the dominant factors affecting the bed topography. There are decreased ILCI by depth deepen, and the basal roughness of three types in the region of 0-430km (Zone A and Zone B). 460-750km (Zone C) region has medium ILCI, where the internal layers are distributed parallel to the bed with low-amplitude and long-wavelength bed topography. The buckled and disrupted internal layers of the upstream of Lambert Glacier (Zone D) indicate that there may be a fast or enhanced ice flow before. Dome A (Zone E) has clear internal layers with high ILCI, and the ratio of the continuous and visible layer thickness is more than 90%., and the rough bed topography is largely unaffected by slow ice flow.
There are many uncertainties in the results of the discussion with a single radar route in this paper, such as ice flow direction, basal conditions, and surface accumulation, and the understanding of the internal structural properties and the basal roughness is not yet complete. However, our results provide an opportunity to discuss the internal mechanisms and ice flow history in PEL of EAIS. It also provides a data foundation for ice sheet modeling and age calculations. In future research, airborne IPR data may cover this rote as well as the entire PEL. By comparison, it will provide a wide range of basal roughness and ice flow history of EAIS for further discussion. We will select higher-dimensional geophysical observations and modeling to help explain the ice dynamics of the PEL, East Antarctica. .