High-resolution subglacial topography around Dome Fuji, Antarctica, based on ground-based radar surveys conducted over 30 years

. The retrieval of continuous ice core records of more than 1 Myr is an important challenge in palaeo-climatology. For identifying suitable sites for drilling such ice, the knowledge of the subglacial topography and englacial layering is crucial. For this purpose, extensive ground-based ice radar surveys were done over Dome Fuji in the East Antarctic plateau during the 2017–2018 and 2018–2019 austral summers by the Japanese Antarctic Research Expedition, on the basis of ground-based 20 radar surveys conducted over the previous ~30 years. High-gain Yagi antennae were used to improve the antenna beam directivity and thus attain a significant decrease in hyperbolic features of unfocussed along-track diffraction hyperbolae in the echoes from mountainous ice-bedrock interfaces. We combined the new ice thickness data with the previous ground-based data, recorded since the 1980s, to generate an accurate high-spatial-resolution (up to 0.5 km between survey lines) ice thickness map. This map revealed a complex landscape composed of networks of subglacial valleys and highlands, which sets. Based 25 on the new map, we examined the roughness of the ice-bed interface, bed surface slope, the driving stress of ice, and the subglacial hydrological condition. These new products and those analyses set substantial constraints for identifying possible locations for new drilling. In addition, our map was compared with a few bed maps compiled by earlier independent efforts based on airborne radar data to examine the difference in features between sets of the data. Recognizing the improvement in data quality and comparing existing topographic products, we suggest that widely available bed topography products should be validated with in-situ observations where it is possible. conventionalreal aperture radar soundings. The present study uses conventionalreal aperture radar sounders with high-gain Yagi antennae, resulting in improved antenna beam directivity and thus a decrease in hyperbolic features in echoes from the mountainous ice-bedrock interfaces. KS, and TS conducted the field survey. SF analysed the ice radar data with support from TO. ST carried out the gridded ice thickness data generation with input from AA, TO, and 515 FS. ST and SF wrote the manuscript with input from all co-authors.

1 Introduction   (Fig. S1a). The bed topography map compiled by Karlsson et al. (2018) cannot identify the bedrock 65 undulations of a horizontal scale less than approximately 10 km due to the spatially sparse (typical survey line spacing of 10 km) observational data and using interpolation with solidintense smoothing to fill in blank areas of data. These authors suggested primary areas for potential drilling sites on the subglacial mountains based on the distinction of melted/frozen bottom by ice sheet modelling under the assumption of surface mass balance and geothermal heat flux. Detailed basal topography data with higher spatial resolution are essential to provide 70
New ice-penetrating radar surveys with high spatial resolution and the compilation with existing data have been necessary to improve the ice thickness maps, which would facilitate the identification of candidate sites for the oldest ice drilling. 80 However, many earlier radar sounding data showed unfocussed along-track diffraction hyperbolae at the ice-bedrock boundary, resulting in potential biases in the estimated ice thickness. The hyperbolic effects mask sharp peaks and steep valleys near such peaks, resulting in underestimation of ice thickness. Focused synthetic aperture radar (SAR) processing can be used to correct the errors (e.g., Young et al., 2017;Rodriguez-Morales et al., 2020). Another approach is to make the footprint of the radar beam smaller with conventionalreal aperture radar soundings. The present study uses conventionalreal aperture radar sounders 85 with high-gain Yagi antennae, resulting in improved antenna beam directivity and thus a decrease in hyperbolic features in echoes from the mountainous ice-bedrock interfaces.
For identifying suitable sites for drilling very old ice, gaining knowledge of the subglacial topography and englacial layering is crucial. Here, we present an ice thickness dataset based on ground-based measurements with a high spatial resolution of up to 0.5 km in terms of survey line spacing for the Dome Fuji region. We combined the latest data from ground-based surveys 90 during the 2017-2018 and 2018-2019 austral summers with earlier data obtained from ground-based surveys carried out by the JARE from the 1980s until 2013. We demonstrate how the selection of antennae affects ice thickness assessment in the southern region of Dome Fuji. We constructed new 0.5 km gridded ice thickness data with the same scale as the line spacing (up to 0.5 km) of measurements during the 2017-2018 and 2018-2019 austral summers. We compare our compilation with a few bed maps compiled by earlier independent efforts, which are mainly composed of airborne radar data, to examine the 95 difference in features between the datasets. Understanding this difference will facilitate the merging of JARE ground-based data and other airborne data in the future. We use the new dataset to examine subglacial conditions to gain further knowledge for identifying the presence and distribution of old ice. Finally, we suggest that the candidate area for the deep coring of old ice is much more limited than previously estimated.

Study area
The JARE's study area around Dome Fuji is situated on the East Antarctic plateau (Fig. 1a). The region shown in Fig. 1b  which covers the southern region of the highest dome summit, has about 12000 km 2 and an elevation range of about 3700 to 3810 m a.s.l. An annual mean air temperature of −54.4°C was observed at Dome Fuji station (Kameda et al., 2009). The annual 110 accumulation rate ranges from about 24 mm w.e. a −1 to 10 % below this value (Fujita et al., 2011) (Fig. S2), as estimated from snow stakes and microwave radiometry.
In the 2017-2018 austral summer, the JARE conducted ground-based radar surveys over a 12000 km 2 area with a spacing of 5 km or less. The total length travelled was approximately 2950 km, covering the vicinity of Dome Fuji station and the southern subglacial mountains (Fig. 1b). Based on the results of these surveys, a more detailed radar survey was conducted on 115 the south side of Dome Fuji in the 2018-2019 season as a collaborative research project by the University of Alabama, the University of Kansas, the National Institute of Polar Research (NIPR), and the Norwegian Polar Institute (NPI). Radar data for a total of 2700 km were acquired in a 1000 km 2 area, covering the subglacial mountain range around NDF site (77°47′19″.789°S, 39°03′10″.053°E, 3763 m a.s.l.). The final spacing between the survey lines was in the range of 0.25-0.5 km. In the present study, we focus on the performance of the conventionalincoherent JARE pulse-modulated radar. Ice 120 thickness detected with the wideband radar sounder provided by the Center for Remote Sensing of Ice Sheets (CReSIS), the University of Kansas, is discussed elsewhere (Rodriguez-Morales et al., 2020).

Ice radar systems
The JARE uses conventionalincoherent pulse-modulated VHF radar sounders with a peak transmission power of 1 kW (Table   S1). Radar systems with various centre frequencies (179, 60, 30 MHz and other values) have been used to investigate the 125 mechanisms of radio wave reflection and the frequency dependence of birefringence in the ice sheet in addition to ice thickness measurements (e.g., Maeno et al., 1994Maeno et al., , 1995Maeno et al., , 1996Maeno et al., , 1997Fujita et al., 1999Fujita et al., , 2002Fujita et al., , 2003Fujita et al., , 2006Fujita et al., , 2011Fujita et al., , 2012Matsuoka et al., 2002Matsuoka et al., , 2003. A transmitter pulse width of 60, 250, 500, or 1000 ns was chosen depending on the scientific objectives (measurements of ice thickness or internal layers) and logistics limitations (i.e., antenna size). The JARE often chooses a wider pulse for ice thickness measurements to detect the target bed topography with higher energy electromagnetic waves. The radar 130 systems were mounted on a snow-tracked vehicle (Ohara Co. SM100 S-type) (Fig. S3). The transmitting and receiving antennae were attached to either side of the vehicle. The polarization plane of the antennae was arranged either parallel or perpendicular to the vehicle movement direction. The JARE uses Yagi antennae and thus the radio waves are linearly polarized, with the polarization plane matching the antenna plane. In the field observations from 1992 to 2013, as shown in the radar survey lines in Figure S1b, the antennae consisted of one or two antenna stacks and three or eight antenna elements per stack 135 (either for transmitting antennae (Tx) or receiving antennae (Rx)), resulting in a total antenna gain (Tx + Rx) that ranges from c.a. 15 to c.a. 33 dBi. For the two antenna stacks with eight antenna elements each, the half-power beamwidth was typically ±20° and ±10° for the E-and H-planes, respectively. For the one antenna stack with three antenna elements, the half-power beamwidth was typically ±35° and ±45° for the E-and H-planes, respectively. The number of antenna stacks and elements determines the antenna gain and directivity. With radar, we need to detect significant echoes from the bed and thus higher 140 antenna gain allows thicker ice sheets to be detected. In addition, higher directivity is necessary for resolving complex subglacial topography, such as steep mountains, slopes, or valley basins. To meet these requirements, for the 2017-2018 and 2018-2019 field campaigns, we increased the number of antenna stacks and/or elements of the Yagi antennae. In the surveys, the number of antenna stacks was 4, 2, or 1 either for Tx or Rx. The number of antenna elements was 8 to 16 per stack. These settings resulted in an antenna gain (Tx + Rx) of c.a. 35 to c.a. 27 dBi. The directivity of the radar beam was also improved by 145 the increase in antenna gain; the half-power beamwidth was typically from ±15° to ±20° for the E-plane and from ±5° to ±22° for the H-plane. In addition, our radar platform was situated on the ice sheet surface. Thus, the geometry of the wave spreading effect is expected to be much smaller than that for airborne radar sounding, which typically has a radar height of a few hundred metres above the ice sheet surface. Therefore, the footprint of the radar beam for our ground-based radar platform after 2017 is much narrower than that of our one antenna stack with three antenna elements used before 2014. Similarly, the footprint of 150 the radar beam for our ground-based radar platform should be narrower than that of the airborne radar antenna with a total antenna gain (Tx + Rx) of about 28 dBi (Nixdorf et al., 1999). The data were recorded on a digital oscilloscope connected to a laptop PC situated inside the tracked vehicle. The three-dimensional coordinates of the radar-sampled points were measured by a global navigation satellite system (GNSS) device attached to the vehicle's roof, which was almost at the same height as the radar antennae. 155

Initial data processing
The ice-bed interface was determined by extracting the peak power of echoes from the bottom of the radargram with semiautomatic detection routines and manual modification. A horizontal smoothing filter was applied to the data using a moving average over 7 s, which corresponds to a horizontal distance of about 20 m at a vehicle speed of 10 km h −1 , to increase the signal-to-noise ratio. The ice thickness at the time of observation was converted from the two-way travel time from the surface 160 to the ice-bed interface under the assumption of a propagation velocity of 1.690 × 10 8 m s −1 determined based on an approximation using the relative permittivity of ice (e.g., Fujita et al., 2000;Saruya et al., 20212022). Here, we estimated this velocity using the relative permittivity value of ice (~3.147±0.010) considering depth-dependent variations of crystal orientation fabrics (Azuma et al., 2000;Saruya et al., 2021) and temperature within the ice sheet. (Motoyama, 2007). Timing errors of the initial trigger in the oscilloscope are possible sources of systematic errors. As an initial step, we calibrated the 165 radar echoes based on a downhole radar target experiment at the Dome Fuji ice coring site. This calibration inherently included corrections for the thickness of both firn and bubbly ice. We estimated the systematic error of our ice thickness values to be within ±15 m at a thickness of ~3000 m. The vertical resolution depends on pulse widths; it was 5 and 21 m for our pulse widths of 60 and 250 ns, respectively (Table S1). The bed elevation was obtained from the difference between the GNSSderived surface elevation and the measured ice thickness. Radar systems used before 2014 had relatively low antenna gains, 170 so the detected bed topography can contain more significant uncertainty than that in the data acquired after 2017. We used the data obtained after 2017 (JARE59 pol179, JARE59 vhf179, and JARE60) as a reference to calibrate the ice thickness obtained from radar measurements before 2014 (JARE33, JARE37, JARE40, JARE49 pol179, and JARE49 vhf60). We considered data within 50 m in horizontal distance from each other to be crossover points. We examined the difference in the ice thickness at the nearest points within the area. A potential bias in the ice thickness was examined for JARE59 and JARE60 data through 175 crossover analysis. For example, JARE59 pol179 data were examined using JARE59 vhf179 and JARE60 data. The ice thickness data showed biases ranging from −37 to 75 m. The ice thickness data acquired from each measurement were corrected for these biases. We subsequently combined the ice thickness point data from each observation into a single dataset.
The combined ice thickness data were interpolated to a 0.5 km resolution grid using an ordinary kriging interpolation method in the open-source GIS software SAGA GIS (http://www.saga-gis.org, last access: 13 August 202130 March 2022). The 180 method is based on the experimental variogram with a lag distance of 1 km. The experimental variogram is fitted to a linear model whose parameters are determined by minimizing the average squared difference between the observational variogram and the model. We set a search distance of 6 km and a maximum data number of 600, which were determined to be the minimum values required to generate weakly smoothed gridded data over the measured area without data gaps.

Assessment of uncertainties in ice thickness 185
We examined the uncertainties of our gridded ice thickness data. The uncertainties were assessed in terms of three error components, namely (1) the vertical resolution of the radar system, (2) the standard deviations of the ice thickness difference, and (3) the standard deviations derived from the kriging interpolation scheme used for generating the gridded data. We estimated the total measurement error at sampled points from components (1) and (2)    The uncertainties of the JARE data were also examined by comparing the gridded ice thickness map with measurement 200 points from radar survey lines. The ice thicknesses of the gridded data were interpolated at measurement points using a linear interpolation scheme. The differences in ice thickness were then calculated by subtracting the measurement data from the gridded data.   (Fig. 1b). The radargram acquired from the JARE54 radar has many hyperbolic shapes at the ice-bed interface (Fig. 2a). The electromagnetic waves that bounce offreflect the convex terrain mask the electromagnetic wave signals that bounce offreflect the valley terrain and slopes, yielding considerable uncertainty when 215 capturing the convex terrain. The radargram acquired from the 11-element Yagi antenna used in JARE59 shows that the hyperbolic effects decreased dramatically (Fig. 2b) compared with those for JARE54. The concave-convex terrain is more distinct in the bedrock topography. In addition, the shapes of subglacial mountains are sharper. The pre-improved antenna with the hyperbolic effects underestimated the ice thickness. For example, in the range of 35-55 km of Dome Fuji, where the hyperbolic effect was observed at the ice-bed interface, the ice thickness in JARE54 was thinner than that of JARE59 with a 220 maximum and an average of 263 and 21 m, respectively.
The higher return power obtained with the improved antenna gain enhanced the visibility of the internal layer stratigraphy.
For example, we detected several continuous horizontal layers near the ice bottom in the mountain range near the NDF site ( Fig. 2b), where the ice thickness is approximately 2500 m. As imaged by the radargram, the internal ice stratigraphy provides valuable information on the integrity of ice layering and the undisturbed ice column, where coherent scattering appears as 225 layered reflections. The reflection of the layer stratigraphy near the bottom becomes weaker as the ice thickness increases. The continuous horizontal layer from NDF is obscured around Dome Fuji station, where the ice thickness exceeds 3000 m.

Ice thickness and bed topography
The gridded ice thickness constructed from the compiled JARE radar survey data (hereafter, referred to as JARE data) reveal the presence of a complex landscape around Dome Fuji, with ice thicknesses ranging from 1800 to 3400 m and an average 230 thickness of ~2670±70 m for the investigated area (Fig. 3a). Because the surface topography in the study area is relatively flat (Fig. 1b), the ice thickness closely reflects the bedrock topography. In other words, areas with thick and shallow ice indicate valleys and mountain ranges underlying the ice sheet, respectively. A valley terrain with an ice thickness of 2700-3500 m extends from the north to the west of Dome Fuji, and a bedrock bump with the ice thickness of 2300-2700 m exists further to the northwest. The southern mountains of Dome Fuji consist of complex subglacial ridges and valleys with the ice thickness 235 of 1800-2600 m.

Uncertainties in ice thickness
The total uncertainty of ice thickness on the interpolated map consists of three factors. Two of them are associated with individual measurements, and another one is associated with interpolation. The error components (1) the vertical resolution of the radar system and (2) the standard deviations of the ice thickness difference were estimated to be 5-85 m (Table S1) and 240

43-80 m (Fig.3.3 Uncertainties in ice thickness
The total uncertainty of ice thickness on the interpolated map consists of three factors. Two of them are associated with individual measurements, and another one is associated with interpolation. The error components (1) the vertical resolution of the radar system and (2) the standard deviations of the ice thickness difference were estimated to be 5-85 m (Table S1) and 43-80 m (Fig. 4), respectively. The error component (3), the standard deviations derived from the kriging interpolation, was 255 estimated to be 1-34 % (relative standard deviations). Accordingly, the JARE data uncertainties are estimated to be from ±60 to ±80 m with an average of ±71 m (Fig. 3b), which corresponds to within 4 % of the ice thickness in the study area. The uncertainty is typically small near the data points and increases with distance from the data points.
For the 107138 data points along the radar lines, the mean difference between the interpolated and measured ice thicknesses is less than −1 m, with a standard deviation of 35 m (Fig. 5a). The absolute differences in ice thickness are < 50 m for 88 % of 260 the points, and > 100 m for only < 2 % of the points. We do not find clear systematic bias in the gridded ice thickness along radar measurement lines (Fig. 5b). We also examined the overall errors in the ice thickness on gridedgridded data generated by kriging interpolation. A comparison of the gridded data with the measurement data showed an anomaly in ice thickness of < 1 % (Fig. 5c), suggesting that the effect of kriging interpolation is constant for all ranges of ice thickness. Note that there is a significant variation in the measurement data relative to the gridded data in the scatter plot because multiple measurement 265 data are compared with the same gridded data in each cell. The gridded ice thickness is > 50 m thinner than the measurement data are situated along the radar lines of JARE37 and JARE49 that run from Dome Fuji to the east and west, respectively (Fig.   5b). The bias along the JARE49 line can be attributable to the bias in the ice thickness in the JARE49 vhf60 measurement data with respect to the JARE59/60 data, with a mean positive anomaly of 75 m (Fig. 4) (note that the gridded data is corrected for the bias). For the JARE37 line, it is unclear why the gridded data are thinner than the measurement data, because there is no 270 appreciable bias in the JARE37 ice thickness with respect to the JARE59/60 data (Fig. 4).

Improvements in ice thickness estimates with high-gain and high-directivity antennae
The use of high-gain and high-directivity antennae resulted in improved accuracy of bedrock topography detection. A crossover analysis of the ice thickness showed that the standard deviations of the ice thickness difference based on data 275 acquired after 2017 (43-56 m) are smaller than those before 2013 (58-80 m) (Fig. 4). Our data suggest that high-gain and high-directivity antennae provide a significant improvement for ice thickness measurement in mountainous terrains. Our 11and 16-element Yagi antennae had lengths of 5.0 and 5.4 m, respectively, for a frequency of 179 MHz. Even though such vertically long antennae would have difficulty in mounting on ground-based platforms, they provide an important opportunity for improving conventionalincoherent radar systems. Modern coherent radar sounders with phase analysis and SAR processing 280 (e.g., Rodriguez-Morales et al., 2020;Lilien et al., 2021) are able to remove hyperbolic effects. For improving the accuracy of bedrock topography measurements, we can apply radar with high-gain and high-directivity antennae, state-of-the-art modern radar with focused SAR processing, or both.

Comparison with other ice thickness data
We compare our gridded ice thickness with the point ice thickness data acquired with the JARE59 system to assess the extent 285 bedrock topography measurements, we can apply radar with high-gain and high-directivity antennae, state-of-the-art modern radar with focused SAR processing, or both. 295

Comparison with other ice thickness data
We compare our gridded ice thickness with the point ice thickness data acquired with the JARE59 system to assess the extent to which it reveals the fine undulations of the bedrock topography. We also evaluated the JARE gridded ice thickness with Bedmap2, BedMachine Antarctica v1, and AWI ice thickness along with the survey route between Dome Fuji and NDF (Fig.   1bS4). We used a 500 m0.5 km horizontal resolution of BedMachine and AWI data, and a 1 km resolution of Bedmap2 with 300 resampled to a 500 m0.5 km grid. The ice thickness data acquired by JARE along the route was combined into respective ice thickness maps to assess the impact of the smoothing effect in the interpolation scheme used for each ice thickness map.
Bedmap2, BedMachine, and AWI gridded data have compiled point ice thickness data acquired in JARE33 and JARE37 along the survey route between Dome Fuji and NDF (Fig. S1b). This comparison demonstrates differences in ice thickness between gridded data that contain common point data and created by different data compilation, smoothing methods, and parameter settings. Figure 6a shows that our gridded ice thickness data accurately resolve the undulations in the basal topography with a horizontal scale of > 500 310 m0.5 km. The absolute difference in ice thickness between the gridded and radar-sampled data along the survey line is ~158 m, with an average value and a standard deviation of −6 and 44 m, respectively (Fig. 6b), indicating that our data slightly underestimate the ice thickness on this spatial scale. The most significant differences appear in the sharp peaks at basal topographic highs, where the horizontal scale is smaller than the grid size. The mean differences between the radar measurement points and ice thicknesses from Bedmap2, BedMachine Antarctica, and AWI gridded data are 1, −36, and −8 m, 315 with standard deviations of 120, 108, and 148 m, respectively (Fig. 6b). These gridded data are smoother than the radar data and thus tend to overestimate the ice thickness in convex areas of the subglacial terrain and underestimate it in concave areas. For example, BedMachine Antarctica resolves a convex terrain on a horizontal scale of about 5 km at a horizontal distance of 7 and 12 km from Dome Fuji, whereas Bedmap2 and AWI data show a flat terrain (Fig. 6a). BedMachine Antarctica shows a large discrepancy in the vicinity of Dome Fuji, implying that the effectiveness of mass conservation methods of bed interpolation is limited by low ice velocities near the ice divide. 325 We also examined the differences in ice thickness between the JARE gridded data and the data from Bedmap2, BedMachine Antarctica, and AWI over the entire study area. Significant variance was found in all comparisons. The results of a statistical analysis of ice thickness differences are summarized in Table S2. For Bedmap2 and AWI gridded data, the most significant differences appear in the subglacial mountainous terrain south of Dome Fuji (Figs. 7a and 7b7c). A comparison revealed overestimation of the ice thickness in the area by ~510 m, resulted from the JARE60 surveys with a fine line spacing (~0.5 330 km). The ice thickness is underestimated by ~810 m in an area further south of the subglacial mountains, where the subglacial valleys extend from west to east with an ice thickness of > 3000 m. The differences in ice thickness are smaller in the vicinity of the Dome Fuji station and further north because Bedmap2 includes JARE data for this area. Although the AWI data also include the JARE measurement data in this area, the ice thickness obtained from AWI data is > 200 m thinner northwest of the station. The JARE49 pol179 data acquired in this area show a negative bias in ice thickness, with an average value of −37 335 m relative to the data from JARE59 and JARE60 (Fig. 4). Presumably, one reason is that the JARE49 data were included in the AWI data without correction of this bias, resulting in a thinner ice thickness. For BedMachine Antarctica data, the differences in ice thickness are smaller than those for the other two datasets (Fig. 7c7b). Small areas with thicker or thinner ice are distributed in a mosaic pattern. The differences in ice thickness from the three maps relative to the JARE data have small mean values and high standard deviations, suggesting that an interpolation scheme inevitably introduces a significant 340 smoothing effect. Accordingly, our data better delineate the critical features of the subglacial mountain ranges south of Dome Fuji compared with the previously published bedrock topography. The high spatial resolution of the ice thickness data is attributed to the data interpolation smoothing effect being much smaller than that in previous data compilations, which is due to the JARE radar data having a horizontal spacing that is equivalent to the grid size.

Subglacial conditions 345
To gain further knowledge for identifying the optimum drill site, we examined subglacial topographical, glaciological and hydrological conditions with our ice thickness data. Here, we present the roughness of the ice-bed interface, bed surface slope, the stress state of the ice and the subglacial hydrological condition.

Basal roughness and slope
Bed surface roughness aids the ice to develop a complex deformation pattern in the lowermost part of the ice sheet. The 355 complex deformation field could disturb the stratigraphic continuity of the ice core record. We calculated the bed surface roughness along the ice radar surveyed route as the root mean square deviation (RMSD) in detrended bedrock elevation on an 800 m baseline (Shepard et al., 2001;Young et al., 2011Young et al., , 2017. The RMSDs at the 800-m length scale typically ranges from 0-80 m with an average of 26 m in the study area (Fig. 8a).
The RMSD obtained in this study is larger than the results by Van Liefferinge et al. (2018), who used the AWI compilation 360 ice thickness data , probably due to the spatially dense basal topography data. The RMSD is lower in the summit regions and surrounding subglacial basins of mountain ranges. For example, much of the mountain summits around the NDF, the subglacial basin at 20 km south of the NDF, and a midway region between Dome Fuji and NDF exhibit smoother basal topography. Significantly rougher bedrock surfaces are observed on the slopes, valleys and saddles surrounding the mountain range, where the basal topography is highly variable. Basal surface slope shows a similar spatial pattern with the 365 basal roughness (Fig. 8b). Steep basal surfaces are distributed on the ridges and valleys of the mountain range south of Dome Fuji. Basal topography is relatively flat in the basal trough around the Dome Fuji station. The bed slope ranges from 0-14° with an average of 4° in the study area.
Our results show that the basal roughness is generally high in the subglacial mountainous terrain and lower in the subglacial basins. These spatial patterns are consistent with previously reported in Dome Fuji, Dome A, Ridge B, and Dome C in the 370 inland of the EAIS (Bingham and Siegert, 2009;Young et al., 2017;Van Liefferinge et al., 2018). Based on chemical analysis of lowermost part of the Dome C ice core, Tison et al. (2015) proposed that high basal roughness may change the vertical stress field, changing the thickness of the internal layer in the basal ice. However, the low basal roughness may also be an unfavourable condition for preserving old ice. Eisen et al. (2020) found that in areas in East Antarctica with low flow velocities, such as areas near the inland ice divide, there is a weak spatial relationship that the bedrock roughness is systematically smaller 375 where the basal conditions are temperate. Indeed, a spatial correlation between smooth basal topography and subglacial lakes has been investigated in the basal temperate area around Dome Fuji . In the Dome Fuji area, the bed is estimated to be temperate at ice thicknesses more than 2500 m (Fujita et al., 2012), and the presence of subglacial lakes has also been found in these locations (Popov and Masolov, 2007;Karlsson et al., 2018;Van Liefferinge et al., 2018). Even in areas less than 2500 m of ice thickness, where the bed was estimated to be frozen by Fujita et al. (2012), there is a complex 380 distribution of smooth and undulating basal surfaces (see areas within red contours in Fig. 8a). Further discussion on basal roughness combined with discriminant analysis of basal thermal conditions based on new radar data will be required to narrow down the area where old ice may exist.

Driving stress
Basal ice disturbance is likely to occur when the ice is subjected to horizontal stress due to the surface slope. However, the 385 reduction in horizontal stress may also be due to the basal melting, which would eliminate the old climatic record. To examine the influence of the ice stress field on the disturbance of old ice at the base of the ice sheet, we calculated the driving stresses from new gridded ice thickness data in conjunction with satellite-derived ice sheet surface elevations. The driving stress ! (kPa) is defined as where " is the density of ice (910 kg m −3 ), g is the gravitational acceleration rate (9.8 m s −2 ), h is the ice thickness, and is the slope of the ice sheet surface. Surface slopes were derived from the 1 km resolution Cryosat-2 surface DEM (Helm et al., 2014). The surface slope was then interpolated to a 0.5 km resolution to combine with the gridded ice thickness data.
Driving stresses spatially correlates with surface slope (Figs 1b and 8c). We observe that driving stresses are low around the Dome Fuji station and increased with distance. Basal topography around the Dome Fuji station is less undulating, resulting in 395 less spatial variation in ice thickness. Elevated driving stresses were observed over the subglacial mountain ridge around NDF.
These high driving stresses may disturb the basal ice in these locations. Suppose ice flows over beds with heterogeneous basal topography and resistance. In that case, it will cause undulations on the ice surface, affecting the spatial variation of the surface slope, which affects the local driving stresses (e.g., Gudmundsson, 2003;Sergienko et al, 2014). The relationship between undulating basal topography and locally elevated driving stresses can be confirmed by the spatial correlation between basal 400 roughness, basal slope and driving stresses. Meanwhile, regions around the subglacial mountain with low driving stress correspond to local subglacial basins. It implies a condition of basal lubrication, and thus melting may have modified the basal ice.

Hydrological conditions
In the Dome Fuji region, the presence of subglacial water, including subglacial lakes and water-saturated sediments, has been 405 identified based on reflected radio echoes from bedrock surface (e.g., Popov and Masolov, 2007;Fujita et al., 2012;Karlsson et al., 2018). Subglacial water flows from areas of high hydropotential to areas of low hydropotential following to the steepest gradient (Shreve, 1972). Hydraulic potential analysis reveals where subglacial water can migrate and provides essential insights into selecting old ice sites. In the Dome Fuji area, hydraulic potential and the potential subglacial drainage routes have been investigated using airborne radar measurement data . We calculated the hydraulic potential and the 410 potential drainage routes of the subglacial meltwater using new ice thickness data. The hydraulic potential ϕ (m) is defined as where % is the density of water (1000 kg m −3 ) and $ is the bedrock elevation. The potential water drainage routes were calculated by a multiple-flow-direction algorithm using the topographic analysis software package TopoToolbox developed in MATLAB (Schwanghart and Scherler, 2014). We set an upslope area threshold for defining the subglacial stream network to 415 100 connected cells (50 km 2 ).
The spatial distribution of hydraulic potential is generally consistent with the basal topography (Fig. 8d). The water flow routing in the hydrological network initiates over the subglacial mountains. It flows along the valleys to the low subglacial basins, indicating that the geometry of the subglacial hydrological network is dendritic. In the study area, we identified four main drainage systems; flow from the subglacial mountain beneath NDF to the basins located southwest of Dome Fuji (easting 420 of 868 km, northing of 1070 km), west of NDF (825 km, 1070 km), south of NDF (810 km, 1025 km), and east of NDF (860  (Popov and Masolov, 2007;Fujita et al., 2012;Karlsson et al., 2018). Our subglacial hydraulic analysis confirms the features of the basal topography on which subglacial lakes exist, suggesting that subglacial lakes may be present at two other basins and other local topographic lows. Basal ice in the subglacial valleys is considered to be melted, as drainage streams flow even where the ice thickness is less than 2500 m. Subglacial hydraulic analysis using the new ice thickness data suggested 430 the possible locations of old ice to the summits and ridges of subglacial mountains in the southern area of Dome Fuji. Further clarification of the thermal conditions of the ice-bed interface and identification of subglacial lakes is required in this area, using the reflective properties of radar radio signals or ice flow modellings.

Prospects for finding oldest ice around Dome Fuji
The subglacial mountain ranges that extend from the south to southeast of Dome Fuji, where we conducted dense radar 435 measurements, are frozen at the ice-bed interfaces where the ice thickness is approximately < 2500 m (Fujita et al., 2012). We could delineate areas with an ice thickness of < 2500 m with the JARE data, indicating that a frozen ice-bottom area can be restricted above the ridges and summits of the subglacial mountain range (Fig. 8a9a). The subglacial terrain exhibits dense dendritic patterns of ridges and summits (Fig. 8b9b). These areas, where the ice-bed interface should be frozen, appear only slightly eroded by ice sheet flow. On the other hand, selective linear erosion occurred in valleys with large ice thickness, 440 suggesting the further erosion of troughs. Fluvial erosion on the East Antarctic continent before 34 MaMyr, when the ice sheet was formed, is also thought to be a factor in evolving basal troughs (e.g., Jamieson et al., 2005). Regarding the requirements for Note that even in areas more than 2500 m of ice thickness, old ice may also exist where the ice-bed interface is frozen due to low basal geothermal heat flux, and the ice stratigraphy near the base is undisturbed due to low horizontal flow velocity.
Regarding the requirements for potential oldest ice sites, smaller ice thickness is preferred to avoid the possibility of bottom 450 ice melting due to an increase in the insulation effect. In addition to this general view for frozen/melt distinction, our map revealed the presence of complex and steep terrain in the area with reduced errors from hyperbolic effects and high spatial density. This new landscape suggests that finding a potentially suitable drilling site in this area may be more challenging than previously thought.
Even if the bottom of the ice sheet is frozen, the ice sheet is either above bedrock ridges, summits, steep slopes, or troughs. 455 Finding undisturbed ice layers just above the bed is a prerequisite for oldest ice drilling. We expect undisturbed ice layers just above the bed to be more likely located above ridges or summits than steep slopes or troughs because the ice over slopes are under high shear stress and that over troughs are subjected to the effects of complex ice flow and bottom melt.. Our analyses of driving stresses and hydraulic conditions suggested that the ice over slopes is under high shear stress, and that over troughs are subjected to the effects of complex ice flow and bottom melt. Horizontal flow velocity in our study area is < 1 m a -1 (see 460 Karlsson et al., 2018), suggesting that basal ice rheology is dominated primarily by vertical normal stress and horizontal shear stress is relatively small. Under the dominance of the vertical normal stress, horizontal shear appears mainly on subglacial slopes rather than ridges or troughs. Basal troughs are often influenced by basal melt or connected to deeper troughs of more basal melt. Then, troughs tend to be fast pathways for ice flow. Therefore, we suggest that subglacial ridges in our study area are under simple ice flow conditions, compared to slopes or troughs to preserve layered conditions. 465 Low snow accumulation rates are one of the requirements for the presence of old ice in the ice sheet base, as they inhibit vertical ice advection and increase the time of existence of the inner layer. The large-scale spatial pattern of the accumulation rate in the Dome Fuji area decreases towards the inland ice sheet (Fig. S2). On the windward (low latitude) side of an ice divide around the Dome Fuji station, the wind releases precipitation as it upslopes over the ice sheet surface. Then it brings relatively dry air on the leeward (NDF) side due to the orographic effect. For the local scale, previous studies have indicated 470 that the accumulation rate is negatively correlated with surface slope, and there is also a relationship between the surface and the basal slopes (e.g., Black and Budd, 1964;Furukawa et al., 1996;Fujita et al., 2002Fujita et al., , 2011. For example, ground-based shallow ice radar measurements around NDF revealed that the accumulation rate is locally variable ~30 % and negatively correlated with surface slope at a length scale of ~2 km (Van Liefferinge et al., 2021). Further analysis between the detailed basal topography revealed by our observations and the accumulation rate will be required to provide new insights into the 475 complex spatially varying surface mass balance distribution, contributing to selecting possible locations for old ice.
Finding suitable locations for the very old ice drilling in mountainous areas requires a pinpointing assessment to avoid sites with irregular bed topography. This is the basis as to why we needed the data with a significant decrease in hyperbolic features in the echoes from mountainous ice-bedrock interfaces. After locating possible areas for deep drilling based on the present 0.5 km resolution bed topography, further radar sounding with a higher spatial resolution (i.e., 0.1 km or less) is needed to find 480 the best candidate points in the identified areas.

Conclusions
To better understand the detailed bedrock topography for finding potential sites that contain ice that extends to > 1 MaMyr, we conducted ground-based radar measurements with a high spatial resolution across the Dome Fuji region, East Antarctica, in the 2017-2018 and 2018-2019 austral summer seasons. The antenna performance of the VHF radar systems was greatly 485 improved by high gain and high directivity. We constructed an ice thickness map from the improved radar data and previous data collected since the late 1980s. The differences in ice thicknesses of our gridded data and previously published gridded data were investigated in relation to the radar measurement data to assess the impact of the smoothing effect in the interpolation scheme on the bedrock topography. We also examined the spatial distribution of the ice thickness differences between our data and previous gridded data. 490 The data acquired using the improved radar systems had a significant decrease in the hyperbolic features at the ice-bedrock interface and allowed basal topography to be identified with higher accuracy. The improved radar systems clearly resolved the internal layer stratigraphy at the lower part of the ice sheet. The new ice thickness data show the bedrock topography, particularly the complex terrain of subglacial valleys and highlands south of Dome Fuji, with substantially higher detail than that in previously published data. The high spatial resolution of the ice thickness data is attributed to the data interpolation 495 smoothing effect being much smaller than that in previous data, which is due to the JARE radar data having a horizontal spacing that is equivalent to the grid size. Our topographic analysis indicates that steep and rough bedrock surfaces extend at the slopes of the basal mountains, while relatively flat surfaces exist at the mountain summit and subglacial basins. Hydraulic analysis implies the presence of subglacial water flow networks and subglacial lakes in subglacial valleys and basins, where these were not explicitly detected before. The new ice thickness map setsand subglacial environments findings are set 500 substantial constraints for identifying possible locations for oldest ice drilling areas. Finally, recognizing the improvement in data quality and comparing existing topographic products, we suggest that widely available bed topography products should be validated with in-situ observations where it is possible. Financial support. This study was supported by KAKENHI from the Japan Society for the Promotion of Science and MEXT (grant numbers 17H06320, 17H06104, 18H05294, 18K18176, and 20H04978JP17H06320, JP17H06104, JP18H05294, JP18K18176, and JP20H04978). 530