Stochastic Simulations of Bed Topography Constrain Geothermal Heat Flow and Subglacial Drainage Near Dome Fuji, East Antarctica

Topographic variability beneath ice sheets regulates ice flow, basal melting, refreezing processes, and meltwater drainage. The preservation of old ice layers and basal ice stratigraphy is sensitive to these subglacial processes, and Dome Fuji, inland East Antarctica, is one of the few regions where 1.5‐Ma old ice can be preserved for investigating a major climatic change that occurred in the mid‐Pleistocene. We used stochastic simulation methods and radar data to generate an ensemble of simulated bed topography with the continuous and realistic roughness necessary to assess basal conditions. Ensemble analysis reveals the magnitude and spatial distribution of topographic uncertainty, facilitating uncertainty‐constrained assessments of subglacial drainage and topographic adjustments to geothermal heat flow (GHF). We find that topographic variability can lead to widespread local GHF variations of ±20% of the background value, which aggregate to raise the regional value and suggest previously underestimated distributions and rates of basal melting. We also find that survey profile spacing has an increasing influence on topographic uncertainty for rougher bed, deriving an empirical relationship that could guide future survey planning based on uncertainty tolerance.

Dome Fuji, inland East Antarctica, is a candidate site for the existence of ice over 1.5 Ma that could provide continuous climate records spanning the mid-Pleistocene Transition (Fischer et al., 2013).The preservation of old ice and undisturbed basal stratigraphy are sensitive to basal melt, ice flow, and accreted ice from meltwater refreeze (Van Liefferinge & Pattyn, 2013), so constraints on GHF and subglacial hydrology are crucial in this region for characterizing 1.5 Ma ice and understanding catchment-scale ice dynamics.At the interior of ice sheets where ice surfaces are relatively flat, bed topography mostly controls meltwater drainage (Shreve, 1972).Basal meltwater flows in extensive drainage networks (Dow et al., 2022), is stored in and released from subglacial lakes (Livingstone et al., 2022), and accretes onto the base of the ice sheet (Bell et al., 2011).Here, we compiled multiple ice-penetrating radar data sets to constrain ice thicknesses and investigate the bed beneath Dome Fuji (Figure 1a).The region has complex topography with mountains and valleys (Karlsson et al., 2018;Tsutaki et al., 2022), evidence for subglacial lakes (Karlsson et al., 2018;Popov & Masolov, 2007;Siegert et al., 2005), and estimated near-basal temperatures at the pressure melting point (Talalay et al., 2020).This indicates the potential for sensitive basal conditions and active hydrological systems, emphasizing the need for constraining the impacts of topographic variability and uncertainty.
The study region is from 596 to 1,020 km Easting, 816 to 1,240 km Northing, on a Polar Stereographic projection parallel to 71°S (EPSG: 3031).We added a 100 km buffer to this region before clipping available data to prevent too few data being available close to the region boundary when running the stochastic simulations.Two-way travel time to depth conversions for ice-bed picks were standardized (Text S1.2 in Supporting Information S1) between data sets using a radio-wave propagation speed of 1.69 × 10 8 m s −1 , a depth-averaged assumption applied previously in this region (Tsutaki et al., 2022).The ice-thickness data sets were compared at 783 crossover points, yielding median differences −1, 3, 18 m between CReSIS data and AWI, JARE59, and JARE60 respectively (Text S1.3 in Supporting Information S1).

Stochastic Simulation of Bed Topography
Sequential Gaussian Simulation (SGS) is a stochastic method (Deutsch & Journel, 1997) that can be applied to simulate values between measurements based on the statistical characteristics of nearby data (e.g., Graham et al., 2017;Law et al., 2023;MacKie et al., 2020MacKie et al., , 2021MacKie et al., , 2023)).Ice-thickness data were decimated to 100 m spacing using a median reduction filter, and an SGS algorithm from the python package GStatSim (Mackie et al., 2022) was implemented to simulate values for cells without measurements using 40 conditioning data points selected using a nearest-neighbor octant function (Mackie et al., 2023) with a search radius of 30 km.We empirically determined the optimum parameters for an exponential statistical model based on the experimental semivariogram (Text S2 in Supporting Information S1).The SGS algorithm moves sequentially over a random path on a 2D grid, selecting values for unsurveyed grid cells using a Gaussian probability function.The probability distribution is generated using a local mean and semivariogram estimated from the statistical model, and conditioned by surrounding data points which are sequentially updated with previously simulated values.
We used SGS to generate an ensemble of 100 ice-thickness grids, and bed elevation was estimated by subtracting ice thicknesses from the Reference Elevation Model of Antarctica (REMA) ice-surface elevation, which should have vertical uncertainty below 1 m in this region (Howat et al., 2019).To examine variability in short-scale bed roughness, we also calculated the Terrain Ruggedness Index (TRI), which takes the square root of the sum of squared differences between grid cells and eight surrounding cells (Riley et al., 1999).

Topographic Modification of Geothermal Heat Flow
To simulate the impacts of topographic relief on GHF distribution, an adjustment factor for large-scale GHF fields was calculated following an empirical geostatistical approach (Colgan et al., 2021).A function for the unitless local topographic modification to GHF (ΔG/G) was calculated: where G = modeled large-scale GHF estimate, α = empirically determined local topographic deviation from regional mean required to induce 100% change in local GHF, z = local bed elevation,   = mean elevation averaged over an empirically determined horizontal radius r, and i, j = two-dimensional horizontal grid indices.In reality α will change with variations in thermal conductivity and radiogenic heat production between different geological settings, however, we used parameterizations α = 950 m and r = 5 km based on a synthesis of parameter combinations against both observed and modeled local GHF anomalies that arise from topographic relief (Colgan et al., 2021).

Subglacial Water Flow and Ponding
Subglacial hydraulic pressure potential (ɸ) was estimated for each simulated bed following Shreve (1972).We assumed that water pressure equals ice overburden pressure, a steady-state assumption appropriate for inland Antarctica where high-pressure distributed water systems are generally maintained between channels (Dow et al., 2022).Water routing was determined for gradients in ɸ using a D∞ algorithm (Tarboton, 1997).Locations of hydraulic sinks were used to predict potential sites of subglacial lakes, and fill depths required to reach the sink spill points provided a proxy for lake depths (Arnold, 2010;Willis et al., 2016).Hydraulic sinks generated through stochastic simulation can be artificial at grid cells far from measurements, so median lake depth over the ensemble was used to screen for persistently predicted sites of water ponding.

Bed Elevation, Variability and Uncertainty
We generated an ensemble of n = 100 ice thickness and bed elevation grids (Shackleton et al., 2023), which have elevations between −500 and 2,100 m above the WGS84 ellipsoid.An example result (#001) is presented in Figure 1b and all results are presented in Movie S1.Steep, mountainous terrain with valleys and ridges occurs at a subglacial massif south of Dome Fuji Station, with smoother terrain directly beneath and north of the ice dome and in lowland valleys converging to the southwest.Regional subglacial topography described in previous work (Fujita et al., 2012;Karlsson et al., 2018;Tsutaki et al., 2022) are observed in the results.
The median absolute deviation (MAD) for the bed topography ensemble stack quantifies the variability between simulated topography ensemble members (Figure 1c) and has a regional average of 85 ± 41 m (mean [μ] ± 1 standard deviation [σ]).Areas with radar data are distinct with low MAD, with close to zero variability near dense survey profiles over the subglacial massif, and less than 25 m in the vicinity of other survey profiles.Less than 50 m variability is typically observed within 2 km of measurements, increasing to over 250 m between the most widely spaced survey profiles (10 km spacing).The topographic roughness of nearby measurements is replicated in simulated regions, and the MAD is highest close to the rough bed and generally below 100 m in smoother regions.
Individual results were validated by calculating the difference between ice-thickness measurements (before median filter decimation) and simulated ice-thickness grid cell values extracted at measurement locations (Text S3.1 in Supporting Information S1).Probability density functions for all validation results show little variability (Figure 1d), indicating that error distributions remain consistent across the ensemble.The ensemble mean difference is close to zero but distributed widely (0.9 ± 50.7 m).However, the range of observed ice-thickness within each grid cell is comparable at 37.2 ± 44 m (Figure 1d), suggesting that the uncertainty of simulated topography is similar to the measured variability within individual grid cells.

Topography-Adjusted Geothermal Heat Flow
Bed topography grids were used to estimate local topographic adjustments to GHF (Shackleton et al., 2023) by calculating a unitless adjustment factor ΔG/G for large-scale GHF fields (Figure 2a; Movie S2).The adjustment factor can be ±40% in regions with high topographic relief, with little to no adjustment in regions of smooth bed.At 10s of km-scale, adjustments show that GHF is reduced along ridges and increased along valley floors (Figure 2b).The probability density functions for adjustment factors over ensemble results show negligible variation (Figure 2c), reflecting consistent regional roughness patterns.The regional impact of topographic adjustments therefore remains consistent between results, and ensemble median of 0.6% suggests elevated GHF is expected, with local adjustments mostly ±20% (2σ).

Meltwater Drainage and Subglacial Lakes
Our analysis predicts dendritic networks of streams routing water away from central sectors and the subglacial massif, broadly conforming to ice drainage divides at the regional scale and mostly following bed topography at <10 km scale (Figure 3).Some dense grouping of predicted streams are observed, which indicate consistent water-flow predictions, more commonly associated with regions with dense measurements, for example, at lower elevations in ice-drainage basin-3 (B3 in Figure 3).Spatially distributed stream predictions are ubiquitous and occur between otherwise stable stream networks, for example, at Easting 700 km and Northing 1,050 km, where inconsistent water routing coincides with a larger gap between survey profiles (Figure 1a).Other regions contain inconsistent water routing despite regular spacing between profiles, coinciding with high topographic uncertainty where streams diverge around an obstacle (e.g., northern B10) and low topographic uncertainty within lowlands (e.g., southern divide of B3/B10).Ensemble analysis indicates potential sites for subglacial lakes at the bed which have ensemble median depths up to 27 m and extents ranging from 0.25-395 km 2 .Using ensemble median values excludes infrequently predicted lakes in interpolated regions, resulting in more lakes near survey profiles.

Topographic Variability and Subglacial Processes
All topographic features with scales of 10s km are consistently observed across the ensemble, and regional patterns of GHF adjustments and subglacial water drainage remain consistent.The exact geometry and location of <10 km-scale topography can vary between simulations (Movie S1) and local-scale features are affected by topographic uncertainty (Figures 2a and 3), which is highest where the bed is rougher and/or measurement density is lower (Figure 1c).Region-wide small-scale bed roughness has consistent distributions between simulated beds (Figure 4a), with ensemble median TRI of 213 ± 117 m (1σ).Using a single variogram for the entire domain implies that inter-region variability in roughness is likely underestimated.The regional mean TRI for the same domain in Bedmachine v3 is 34 ± 28 m (1σ) and Bedmap2 is 46 ± 36 m (1σ) due to smoothing in interpolated regions, suggesting that bed topography in this region could be an order-of-magnitude rougher with larger GHF variability than shown in other interpolated bed topography maps.
Subglacial drainage analyses predict potential locations and extents of subglacial lakes and show drainage sensitivities to topographic uncertainty (Figure 3).There are 20 subglacial lakes inferred from radar data in this region (Karlsson et al., 2018;Livingstone et al., 2022;Wright & Siegert, 2012), of which 11 directly coincide with our ensemble median lake predictions.Radar-derived bed elevation inherently masks the bathymetry of radar-detected lakes, so the topography analysis in this study is unlikely to predict all radar-detected lakes in bankfull condition.Of the nine radar-inferred lakes that are not predicted in this study, four are located upstream of consistently predicted streams which delineate drainage pathways.Uncertainty-constrained subglacial water flow directions also outline potential connectivity between drainage systems and across catchments that otherwise might not be detected using a single bed interpolation.We identified several drainage pathways across the continental ice-flow divides which reflect regionally low ice-surface slopes, and predicted lakes close to grid-south Radar-inferred subglacial lakes as yellow triangles at lake center coordinates (Siegert et al., 2005) or lines along-profile (Karlsson et al., 2018).Ice sheet drainage basins delineated and labeled in white from Zwally et al. (2012).
10.1029/2023JF007269 7 of 11 have uncertain drainage toward either Basin-3 or Basin-10.The topographic uncertainty analyses reveal regions where an ensemble approach to predicting subglacial processes could be necessary or more data are needed to reduce uncertainty, especially important in regions where basal conditions are sensitive to topographic variability.

Basal Conditions and Presence of 1.5 Ma Old Ice
The preservation and stratigraphy of deep ice layers are sensitive to basal melting or refreezing processes (Fischer et al., 2013), so diagnosing basal conditions is crucial for deep ice core site selection.Our topographic adjustments to GHF suggest that local basal conditions could be sensitive to km-scale physiographic settings and aggregated effects lead to regionally elevated GHF, indicating that the use of large-scale GHF models may lead to inaccurate delineation of thawed bed and estimates of basal melt rates.Thermodynamical modeling suggests 3.1% of our study area had a frozen bed for the past 1.5 Ma which could preserve old ice (Van Liefferinge et al., 2018).The subglacial massif is well-constrained by measurements, and we consistently estimate local GHF adjustments reaching ±30%.The DF2 core was drilled in a locally elevated GHF region ca.3% higher than background values (Figure 2e), equating to 1.5-2.2mW m −2 when applied to the GHF range from the three borehole temperature estimates.Recent 1D ice flow modeling by Obase et al. (2023) showed that GHF differences of 2 mW m −2 or 100 m change in ice thickness can lead to basal temperature variations of 2°C.Many localities even within 10 km of the drill site could have topographically adjusted GHF over 5% lower than the regional value (Figure 2a), a crucial difference here where basal temperatures are close to or at the pressure melting point.Obase et al. (2023) found that a GHF of 60 mW m −2 most accurately reproduced the age distribution of layers observed in the ice core.Previous borehole-based GHF estimates using an ice thickness of 3,090 m reported earlier (Ageta et al., 1998) were lower at 59 mW m −2 (Hondoh et al., 2002) and 50.4 mW m −2 (Mony et al., 2020), though a revised ice thickness of 3,028 ± 15 m (Fujita et al., 2006) was used by Talalay et al. (2020) to yield an estimate of 78.9 ± 5.0 mW m −2 .Our data compilation suggests that radar-derived ice thicknesses within 500 m of the DF1 core site are 3,018 ± 16 m (1σ), based on averaged values from four independent surveys (Table S2 in Supporting Information S1).Over the subglacial massif, Van Liefferinge et al. ( 2018) estimated the threshold GHF required to keep basal ice frozen over the past 1.5 Ma is ca.65 mW m −2 , which can be lowered to 55 mW m −2 considering new surface mass balance data averaged over the past 300 years as a long-term regional value (Van Liefferinge et al., 2021).Some pan-Antarctic GHF models show that the revised threshold is met over the subglacial massif, although 2 out of 3 local GHF estimates derived from borehole-temperatures suggest regional GHF could be over the threshold despite our topographic adjustment factor (Figure 2d).
No large relief features are observed in radar data or simulated beneath the current ice dome, suggesting that local relief may not be a significant control on the current ice dome location and that the location could have been different in the past.Thermomechanical modeling suggests that ice thicknesses could have fluctuated up to 250 m at Dome Fuji Station over the last 800 ka (Parrenin et al., 2007).Re-routing of water between catchments is plausible with shifting overlying ice, and historical water storage and drainage patterns could vary significantly from present-day conditions, which could affect the delineation and characterization of 1.5 Ma old ice.Meltwater refreezing onto basal ice also has significant potential to disturb basal ice stratigraphy, as observed in the nearby Gamburtsev Mountains (Bell et al., 2011).Our bed topography ensemble could be utilized further to constrain potential basal ice flow disturbances to layer stratigraphy through ice flow modeling.

Survey Spacing Required for a Given Bed Uncertainty Tolerance
The wide-ranging radar survey spacings in this region (0.25-10 km) allowed us to explore the relationship between distance d from radar data, bed roughness (TRI), and topographic uncertainty (MAD).Based on ensemble-averaged TRI and survey locations, we manually identified regions (Figure 1c: white boxes) that have characteristically rough bed with dense (grid-center) and sparse (grid-southeast) measurements, and smooth bed with both dense (grid-west), and sparse (grid-south) measurements.TRI in these regions is distinct across the ensemble (Figure 4a).MAD increases with increasing d but plateaus for d > 6 km, at ca. 200 m for rough bed and ca.75 m for smooth bed (Figure 4b).Both examples with distinct data density show similar distance dependencies for MAD, which implies a bed-roughness-dependent relationship between topographic uncertainty and survey spacing.Regression analysis between TRI and MAD for data grouped into 0.5 km d intervals for 0 < d < 10 km (Figure 4c) demonstrates that MAD increases with TRI for each d interval, and the rate of increase is enhanced for larger d.Linear models were fit to binned data with y-intercept fixed at zero assuming MAD approaches zero for smooth surfaces, showing that MAD ≈ βTRI, where β represents the dependence of topographic uncertainty on bed roughness.
Figure 4d shows the distance dependence for β, and the empirical relationship could be used to guide survey planning given an initial estimate for regional roughness based on existing data or preliminary surveys.Interpolated grids based on sparse data lack roughness between survey profiles and are not suitable for preliminary roughness analyses.Using a derived TRI, the separation between survey profiles (i.e., twice the distance to data points, 2d) can be prescribed for a given uncertainty tolerance (Figure 4d).If bed roughness is similar to the Dome Fuji region (median TRI = 213 m), survey profiles should be separated by 1.2, 5.7, and 27 km for uncertainty tolerance of 50, 100, and 150 m, respectively.However, accounting for the regional standard deviation of 117 m, the TRI could be 330 m, which would reduce the survey spacings to 0.7, 1.9, and 14 km.Although it is difficult to accurately estimate bed roughness prior to extensive surveys, this relationship could be used to empirically implement scientific requirements into survey planning.Bed topography mapping guided by uncertainty tolerance might be particularly useful in coastal regions of Antarctica, where uncertainty could impact estimated ice 10.1029/2023JF007269 9 of 11 discharge and projected grounding-line retreat of fast-flowing glaciers (Matsuoka et al., 2022).To further develop this approach, the measurement scale and bed roughness metric should be investigated as well as the impacts of grid cell size.New pan-Antarctic data sets (Frémand et al., 2023) could facilitate analyses of this relationship for different regions and glaciological settings.

Conclusions
Stochastic simulations provided detailed bed topography with realistic roughness an order-of-magnitude larger than previous interpolations.Topographic uncertainty is significantly larger in regions with rougher bed, and we quantified the dependence of topographic uncertainty on distance to radar measurements and bed roughness.An empirical relationship can guide the spacings of future surveys in the Dome Fuji region and could be applied and developed elsewhere to plan surveys within a framework of uncertainty tolerance.Constraining the bed to within 50 m requires radar profile spacings of 0.7-1.2km in this region.The impacts of topographic variability on the distribution of GHF and subglacial drainage were investigated through ensemble analysis of simulated beds.Topographically adjusted GHF is elevated in valleys and reduced over ridges, with km-scale variability leading to widespread local GHF adjustments ±20% the regional value, aggregating to enhance overall GHF.These effects could strongly regulate the distribution of basal ice at the pressure melting point and should be included in future estimates.Hydraulic potential analysis reveals meltwater drainage pathways and their relative probability, which maps water discharge from subglacial lakes and basal melt from the upper regions of ice-flow catchments, as well as highlighting regions potentially sensitive to ice dome migration.Our analysis indicates widespread potential for subglacial lakes in the region, which coincide with just over half of the lake locations proposed in previous work and identify potential locations of undetected lakes.The uncertainty-constrained drainage analysis and GHF adjustment approach we apply here in the Dome Fuji region could be applied elsewhere to diagnose basal conditions and delineate 1.5 Ma old ice.This work was supported by funding from the Research Council of Norway (FRINATEK 315246) and is a contribution to EU Horizon 2020 Beyond EPICA Oldest Ice Core (BE-OIC) project (Grant 815384).CReSIS radar data used in this study were collected with financial contributions from the Norwegian Polar Institute and National Institute of Polar Research, Research Organisation of Information and Systems, as a collaborative project with The University of Alabama and the University of Kansas.Japanese Antarctic Research Expeditions, as well as the British Antarctic Survey, the Belgian Science Policy Office, and the International Polar Foundation contributed to fieldworking to collect this data set.We thank E. Mackie, A. Robel, and two other anonymous reviewers for their constructive reviews of this work.We also acknowledge NSF Grant 2126503 and the use of the CReSIS toolbox from CReSIS generated with support from the University of Kansas, NASA Operation IceBridge Grant NNX16AH54G, and NSF Grants ACI-1443054, OPP-1739003, and IIS-1838230.This is BE-OI contribution number 33.

Figure 1 .
Figure 1.Bed topography simulations and analysis.(a) Radar survey tracks, Dome Fuji Station (red diamond) and Reference Elevation Model of Antarctica ice surface contours (Howat et al., 2019).Survey patterns for Center for Remote Sensing of Ice Sheets (CReSIS) and Japanese Antarctic Research Expedition (JARE) 60 are similar, with mostly 0.5 km separation.(b) Hillshaded bed simulation #001 from the ensemble.(c) Median absolute deviation (MAD) between 100 results with vertical and horizontal profiles plotted.Four boxes indicate sample locations (Section 4.3).(d) Simulated bed minus measured ice thickness (100 curves) and measurement range within each map grid cell (bars).All maps projected to EPSG: 3031 with elevations in meters referenced to the WGS84 ellipsoid.

Figure 2 .
Figure 2. Spatial variations in geothermal heat flow (GHF).(a) Topographic adjustments to GHF (result #001) and inset a 10 × 15 km region (red box) around Dome Fuji Station.(b) The central subglacial massif (result #001).(c) Probability density curves for 100 topographic adjustment results, with ensemble-averaged standard deviation (gray), mean (black), and median (green).(d) GHF distributions in the study region from pan-Antarctic models with corresponding value at Dome Fuji Station (circles).Local GHF estimates from borehole temperatures and reported uncertainty are plotted (diamonds), with a box plot showing potential regional GHF range after ensemble adjustments to local value.(e) Topographic adjustment factor over the ensemble at Dome Fuji Station, with ensemble-averaged mean (black), median (green), and standard deviation (gray).

Figure 4 .
Figure 4. Topographic uncertainty dependence on bed roughness and survey spacing.(a) Topographic Roughness Index (TRI) in the Dome Fuji region, for ensemble results (black curves), the sample regions shown in Figure 1c, Bedmachine v3, and Bedmap2 at native resolution.(b) Median Absolute Deviation (MAD) at increasing distances from measurements.Characteristic regions plotted in the same colors as panel (a), with all other locations in gray.(c) MAD for increasing TRI binned at 0.5 km d intervals between 0 and 10 km.Only 4 of the 19 intervals are plotted to show overall trends, with full model statistics in TableS3in Supporting Information S1.(d) Slopes of fitted regression lines in panel (c) at increasing distance to measurements.