A high-resolution synthetic bed elevation grid of the Antarctic continent

Digital elevation models of Antarctic bed topography are smoothed and interpolated onto lowresolution ( > 1 km) grids as current observed topography data are generally sparsely and unevenly sampled. This issue has potential implications for numerical simulations of ice-sheet dynamics, especially in regions prone to instability where detailed knowledge of the topography, including fine-scale roughness, is required. Here, we present a high-resolution (100 m) synthetic bed elevation terrain for Antarctica, encompassing the continent, continental shelf, and seas south of 60 S. Although not identically matching observations, the synthetic bed surface – denoted as HRES – preserves topographic roughness characteristics of airborne and ground-based ice-penetrating radar data measured by the ICECAP (Investigating the Cryospheric Evolution of the Central Antarctic Plate) consortium or used to create the Bedmap1 compilation. Broad-scale (> 5 km resolution) features of the Antarctic landscape are incorporated using a low-pass filter of the Bedmap2 bed elevation data. HRES has applicability in high-resolution ice-sheet modelling studies, including investigations of the interaction between topography, ice-sheet dynamics, and hydrology, where processes are highly sensitive to bed elevations and fine-scale roughness. The data are available for download from the Australian Antarctic Data Centre (doi:10.4225/15/57464ADE22F50).


Introduction
The largest source of uncertainty in projections of sea-level rise to the end of the 21st century is derived from poorly constrained estimates of mass loss from the Antarctic and Greenland ice sheets (Church et al., 2013).As the most vulnerable regions of the Antarctic Ice Sheet are grounded below sea level, the ice-sheet response to climate warming will be determined by dynamics operating at the grounding line (Schoof, 2007b;Drouet et al., 2013).Where an ice sheet rests on a bed topography that is below sea level and deepens towards the ice-sheet interior, marine ice-sheet instability (MISI) could occur, leading to increased ice flow, thinning, and rapid grounding line retreat (Weertman, 1974;Thomas et al., 2004;Schoof, 2007a;Durand et al., 2009;Goldberg et al., 2009;Favier et al., 2014;Joughin et al., 2014).It follows that bed elevation is one of the most important controls in modelling ice-sheet dynamics and constraining estimates of future sea-level rise.

F. S. Graham et al.: Synthetic Antarctic bed elevation
Coordinated international efforts over recent decades have vastly increased the coverage and density of bed elevation measurements in Antarctica (Lythe et al., 2001;Le Brocq et al., 2010;Fretwell et al., 2013).These data have been used to improve the fidelity of gridded digital elevation models (DEMs) spanning the whole Antarctic continent.Building on a 5 km gridded bed elevation DEM (Lythe et al., 2001;Le Brocq et al., 2010), the most recently compiled Antarctic bed topography dataset, Bedmap2, is available at 1 km resolution, having been generated from over 25 million measurements (Fretwell et al., 2013).Nevertheless, much of the Antarctic continent is difficult to access and remains poorly sampled.In such regions, bed elevation DEMs rely on interpolation, resulting in geometric inconsistencies that adversely impact numerical simulations of ice dynamics (Warner and Budd, 2000;Fürst et al., 2015;Gasson et al., 2015).Uncertainties in bed elevation are particularly problematic given that, for much of the Antarctic Ice Sheet, the simulated large-scale velocity field depends only on the local-scale details of the geometry and boundary conditions despite the elliptic nature of the governing equations for ice flow.
Recent effort has focussed on understanding the impact of low-resolution bed elevation data on ice mass flux.Durand et al. (2011) performed a sensitivity analysis of an outlet glacier susceptible to MISI, demonstrating that at least 1 km spatial resolution in bed topography is required for accurate estimates of ice mass flux.However, bed elevation data of a higher resolution than 1 km may be necessary in some applications to capture both the channelised landscape that guides glacier flow and the fine-scale roughness that impacts basal sliding (Goff et al., 2014;Ritz et al., 2015).Importantly, a question that has yet to be addressed for the Antarctic continent is what minimum resolution in bed elevation is required to accurately simulate ice-sheet dynamics.
The purpose of this study is to generate a high-resolution synthetic bed topography dataset for Antarctica (HRES) for investigating the sensitivity of ice-sheet dynamics to bed elevation resolution, including the interaction with subglacial hydrology (Fricker et al., 2007;Goff et al., 2014).We emphasise that this dataset is intended to be synthetic (i.e.HRES is not intended to be a substitute for other bed elevation datasets that preserve the observations) but has covariance properties that are consistent with those of the measured bed elevations from available radar transects.The generation of HRES relies on bed elevation data used to create the Bedmap1 compilation and from the ICECAP (Investigating the Cryospheric Evolution of the Central Antarctic Plate) airborne radar survey where they are available at high resolution.The lowfrequency (> 5 km) component of HRES is identical to a similarly low-pass filtered Bedmap2.HRES covers the same domain as Bedmap2 and is available at a spatial resolution of 100 m.The length scale of the topographic roughness used in this study is limited to 200 m.3).Numbers 1-27 correspond to drainage basins defined by the Goddard Ice Altimetry Group from ICESat data (Zwally et al., 2012).

Data synthesis
A two-step approach was used to generate the highresolution synthetic bed elevation terrain, HRES.First, we simulated a non-conditional "roughness" terrain (i.e. a stochastic realisation of roughness that does not necessarily honour the exact values of the original data) using highspatial-resolution radar data obtained from the 2009-2012 ICECAP campaigns (Roberts et al., 2011;Young et al., 2011;Wright et al., 2012;Blankenship et al., 2011Blankenship et al., , 2012) ) that were used to create the Bedmap1 compilation (hereafter BC1; Lythe et al., 2001).The locations of the data included in this step are shown in Fig. 1.The ICECAP bed elevation data are measured using a High-Capability Radar Sounder (Hi-CARS) high-bandwidth airborne ice-penetrating radar (Peters et al., 2005); BC1 combines data from multiple airborne and ground-based radar sounding campaigns, from a variety of systems.Our method for the generation of the roughness terrain can easily incorporate additional bed elevation data as they become available.Once generated, the roughness terrain was high-pass filtered using a Gaussian kernel with a 5 km half-power cutoff.
Second, the Bedmap2 bed topography DEM was low-pass filtered, using a low-pass Gaussian kernel with a 5 km halfpower cutoff.The two filtered terrains were combined (preserving all wavelengths of the original datasets), resulting in the high-resolution bed topography, HRES.
An alternative method for the production of highresolution (250 m) bed elevation data has recently been applied to the Thwaites Glacier region (Goff et al., 2014).Goff et al. (2014) combine both conditional and non-conditional simulations of a range of data with the intent to avoid the inconsistencies and artefacts introduced through interpolation techniques such as kriging.The resulting terrain is of sufficiently high resolution to facilitate characterisation of the subglacial landforms and landscape of the Thwaites Glacier, which will lead to improved understanding of ice flow and its sensitivities to external forcing.However, the methods used to produce this terrain rely on a higher data density than is available for most of Antarctica.Our methodology was chosen because of its ability to handle spatially sparse data with highly inhomogeneous sampling resolutions while being computationally tractable for all of Antarctica at 100 m resolution.
In the following sections, we provide a detailed outline of the methods used to generate the roughness terrain and to compile the final synthetic bed topography dataset.The pseudo-algorithm for the generation of HRES is provided in Appendix A.

Roughness terrain synthesis
Ideally, the spatial covariance characteristics of the nonconditional roughness terrain (the high-frequency component of the synthetic topography dataset) should match those of ICECAP and BC1.The method of Cholesky decomposition of the observed covariance matrix can be used to produce such correlated data (Davis, 1987).Specifically, the positive definite covariance matrix C calculated from the ICECAP and BC1 datasets can be decomposed into lower L and upper U triangular matrices: where L has real and positive diagonal entries and U is the transpose of L. This method results in a unique decomposition for positive definite matrix C. Now, given a vector z of uniformly distributed random numbers with zero mean, we find The product Lz can be used to construct a non-conditional realisation of bed topography, of which the covariance structure is consistent with ICECAP and BC1.We next describe how the covariance matrix C and resulting simulated roughness terrain are calculated.

Covariance structure
In order to perform the Cholesky decomposition, Eq. ( 1), we first calculated the covariance distribution for the ICECAP and BC1 datasets.The along-track covariance distribution for each flight or traverse line was estimated using 16 km sliding windows with 8 km offsets.For each window, the along-track data were averaged into 100 m bins and the following exponential decay model was fitted to obtain the covariance C(d) (Goff and Jordan, 1989): In Eq. ( 3), d is the along-track distance (i.e. between two points in our bins), A is the topographic variance, and 1/B is the e-folding distance.Both A and B are free parameters that were fit to minimise the RMSE between C(d) and the observations.It is C(d), and not A or B, that is required in the generation of HRES.To ensure that the data density was approximately consistent for each calculation of the covariance distribution, windows were included only if data were present in at least 90 % of the 100 m bins.Additionally, data were omitted from the calculation if A < 0, A > 500, or the ratio of A to the maximum covariance in any 100 m bin was outside the range [0.33, 3] (the latter condition ensured a reasonable fit to the exponential model).A total of 9272 points satisfied the criteria for inclusion in the non-conditional roughness terrain (Fig. 1).

Cholesky decomposition
The covariance matrix C defined by the exponential model in Eq. ( 3) is necessarily symmetric and positive definite.As such, Cholesky decomposition can be applied to C to obtain the lower triangular matrix and its transpose.For a box with 16 km side lengths and easting and northing coordinates centred on each of the 9272 valid data points from Sect.2.1.1 and divided into 100 m bins, a covariance matrix C was calculated using coefficients A and B from Eq. (3), and with d varying appropriately.Cholesky decomposition was applied to each of these covariance matrices C, yielding lower triangular matrices L l (here, we have used the superscript "l" to denote the fact that these lower triangular matrices are local -i.e.defined over a box with 16 km side lengths centred on each of the 9272 valid data points and comprising 161 × 161 cells).
Next, each L l was multiplied by a uniform random matrix with zero mean -denoted as z l -defined over the same spatial domain as L l .Each local matrix z l was a subset of a uniform random matrix -z -spanning the same spatial domain as Bedmap2 but at 100 m resolution (i.e. a matrix of 66 661 × 66 661 grid cells).We generated the Cholesky decomposition roughness terrain (CDRT) on the same spatial domain as z by calculating the average of the squared and weighted inverse distance of the 20 closest values from the product L l z l .The choice of 20 artefacts minimised by squared and weighted inverse distance points was associated with sparse data.In this way, CDRT has a spatial covariance structure consistent with that of the original ICECAP and BC1 data.Note that although CDRT is one realisation of a near-infinite number of unique roughness terrains, this realisation suffices for our purposes (namely, to generate a synthetic bed topography dataset suitable for investigating the impact of resolution and roughness on ice-sheet dynamics).The calculation of CDRT was spatially independent for each data point; thus, computational parallelism through the use of OpenMP directives was utilised to reduce computational time (which was on the order of 2000 CPU hours).

Compilation of HRES
Due to the statistical properties of large samples of distributions, the bed elevation extrema from CDRT were −72 897 and 70 838 m, well outside the observed range from the ICE-CAP/BC1 measurements of −3373 to 3380 m.To address this, we defined a scaling factor based on a comparison of roughness values from the original and simulated datasets.Roughness G is defined (Shepard et al., 2001) as the root mean squared deviation between points of detrended bed elevation y separated by a lag ( x), We used a lag x = 1600 m, consistent with that used by Gooch et al. (2016).Roughness values were calculated using Eq. ( 4) for each of the ICECAP/BC1 flight lines separately and for the points in CDRT that overlaid these data.A linear fit to the spread of roughness values from ICECAP/BC1 and CDRT yielded a slope of 14.42 and an R 2 value of 0.49 (significant at the 95 % confidence interval using a two-sided Student's t test) for the correlation between the observed and predicted values (Fig. 2a).We used the median value of the ratio between ICECAP/BC1 and CDRT roughnesses -approximately 22.87 -to uniformly scale CDRT.This median value was close to the ratio of CDRT to ICECAP/BC1 extrema of 21.29.Once scaled, CDRT was high-pass filtered using a Gaussian kernel with a 5 km half-power cutoff.The corresponding Gaussian kernel was used to low-pass filter the Bedmap2 DEM, which was first interpolated to the same 100 m grid as CDRT.The two filtered datasets were then added to produce HRES.

Results
The HRES terrain is plotted below Bedmap2 in Fig. HRES was generated from a non-conditional simulation of the ICECAP/BC1 data that is unlikely to honour the exact values of the underlying data.For this reason, the magnitude of the differences between HRES and Bedmap2 is not necessarily the most robust measure of the quality of HRES.Instead, the extent to which the distribution of D differs from a normal distribution provides an indication of the fidelity of HRES to the original ICECAP/BC1 data.We calculate the deviation of the distribution of D from the normal distribution using the D'Agostino-Pearson K 2 test (D'Agostino et al., 1990).The test statistic K 2 is approximately chi-square distributed with 2 degrees of freedom.K 2 calculates deviation from normality as a result of skewness and/or kurtosis and is defined as where Z( , and Z(b 2 ) is a test of kurtosis (b 2 ).The test statistic is calculated for each of the Antarctic drainage basins defined using Ice, Cloud, and Land Elevation Satellite (ICESat) altimetry (Table 1; Zwally et al., 2012), and the corresponding distributions of D are compared with the normal distribution in Fig. 4.
The deviation of the distribution D from the normal distribution N is a result of the covariance structure of the underlying observations used to construct HRES.In particular, we note marked differences between the distribution D and the normal distribution N, where (1) more ICECAP/BC1 data are available and/or meet the criteria for inclusion in the simulation of CDRT and (2) a basin covers terrains of highly contrasting roughnesses (e.g.basin 17, which spans part of the relatively smooth Ross Ice Shelf as well as the Transantarctic Mountains).In East Antarctica, the distributions of D and N differ the most in ICESat basins 12-17, which include areas of Wilkes Land and the northern tail of the Transantarctic Mountains, and an area within Palmer Land in the Antarctic Peninsula (basin 24).Basin 21 is the only basin that is not statistically significantly different from the normal distribution at the 95 % confidence interval.Nevertheless, the distribution of D is generally closer to the normal distribution in regions with the poorest ICECAP/BC1 data coverage, including much of West Antarctica (basins 1, 20-23, and 25, which encompass Marie Byrd Land and the Siple Coast, Ellsworth Land, and the Filchner-Ronne Ice Shelf) and basins 5-9 in Queen Maud Land, East Antarctica.Sharply peaked distributions of D in basins 17-19 delineate smooth terrain over the Ross Ice Shelf from rougher, continental terrain.
Differences between ICECAP/BC1 data points and the corresponding overlay points in HRES and Bedmap2 along 18 selected ICECAP/BC1 flight or traverse lines are compared in Fig. 5.These flight-traverse lines encompass a  5) for each of the drainage basins 1-27 in Fig. 4. The √ b 1 and b 2 statistics are the bases for tests of skewness and kurtosis, respectively.For a normal distribution, K 2 is approximately chidistributed with 2 degrees of freedom.range of landscapes, from smooth subglacial basins to highelevation highlands.For over half of the selected flighttraverse lines, the along-track roughness values from HRES are within 20 % of the corresponding roughness values from ICECAP/BC1.Flight lines O, P, and R show the poorest agreement in roughness between HRES and ICECAP/BC1, with roughness values deviating by more than 50 % of the higher value in each case.However, flight-traverse lines O, P, and R are derived from regions with a paucity of highresolution data available for inclusion in the generation of CDRT (flight lines O, P, and R were themselves not included in the generation of CDRT for this reason).As expected, where Bedmap2 data are in better agreement with ICECAP/BC1, the normalised along-track RMSE between HRES and ICECAP/BC1 is minimised (Table 2).This relationship holds independent of the underlying terrain roughness.

Errors and uncertainties
Sources of uncertainty exist in the datasets, methods, and processes used to generate HRES.We do not quantify these errors explicitly because HRES is a synthetic terrain that has been generated predominantly for investigating the impact of resolution and roughness on ice-sheet dynamics, rather than as a realistic, specific representation of Antarctic bed topography.Nevertheless, the following caveats should be considered in the generation of HRES: i.The roughness terrain incorporated in HRES is a non-conditional simulation of high-resolution flighttraverse line data from the ICECAP and BC1 compilations, which are themselves sparsely available over the Antarctic continent (Fig. 1).The flight-traverse line data have associated errors from instrumentation and processing (e.g.Peters et al., 2005) -these errors will propagate through the simulation of HRES.
ii.In order to generate the high-frequency roughness terrain, we assume that topographic variance is a smoothly spatially varying function.In reality, we have too few high-resolution data points to adequately assess the rigour of such an assumption.
iii.Only ICECAP/BC1 data of sufficiently high resolution (i.e.greater than 100 m resolution, chosen as it is twice the Nyquist frequency of the observations) were included in the simulation of HRES.This limits how well the final HRES dataset matches the observations, especially in regions of West Antarctica.The roughness terrain will be updated to incorporate additional highresolution bed elevation data as they become available.
iv.The Bedmap2 DEM, of which the low-pass component is included in the generation of HRES, suffers from artefacts through the particular gridding and interpolating methods used compared with other ice thickness interpolation methods, especially in regions with no nearby measurements (Roberts et al., 2011).
v. The non-conditional simulation technique based on the Cholesky decomposition of ICECAP/BC1 covariances makes a number of assumptions that influence the outcome bed elevations (notably, that the original data are isotropic and that high-frequency noise is normally distributed).
vi. HRES is simulated using data that are not independent.
vii.It is possible that even finer-scale topographic features than those captured in HRES play a role in modulating ice dynamic processes (e.g.Ritz et al., 2015).This has implications for the degree to which future modelling will ascertain what resolution in bed topography is enough for consistent and accurate simulations of ice dynamics (i.e.we can only assess the impact of bed topography features of a scale greater than 100 m).We will explore this further in subsequent studies.

F. S. Graham et al.: Synthetic Antarctic bed elevation
We refer to the original dataset and method papers for a more detailed discussion of errors inherited by the HRES dataset from the underlying terrains (Alabert, 1987;Davis, 1987;Bourgault, 1997;Lythe et al., 2001;Le Brocq et al., 2010;Young et al., 2011;Fretwell et al., 2013).

Data availability
A detailed description of the data used can be found in the first paragraph of the conclusions.

Conclusions
The result of this study is a 100 m resolution gridded synthetic Antarctic bed elevation terrain -referred to as HRESthat has been made available for download at the Australian Antarctic Data Centre (doi:10.4225/15/57464ADE22F50). HRES combines a high-frequency non-conditional simulation of bed elevation with the low-frequency component of the Bedmap2 bed elevation terrain.This dataset is available in NetCDF standard format on a 100 m resolution grid in a polar stereographic projection (central meridian 0 • , standard parallel 71 • S) with respect to the WGS84 geoid.The 100 m grid has 66 661 rows by 66 661 columns, where the corner of the lower left cell is located at a polar stereographic easting and northing of −3 333 000 and −3 333 000 m, respectively.The value for missing data is −9999.The file size is approximately 17 GB.
HRES is not intended as a realistic depiction of highresolution Antarctic bed topography and is, therefore, not meant as a substitute for datasets such as Bedmap2 (although, the low-frequency component of HRES is identical to the Bedmap2 bed elevation dataset).Instead, HRES is a synthetic terrain generated for the specific purpose of assessing the sensitivity of ice-sheet dynamic processes to the resolution of the underlying bed topography.The sufficiency of the resolution of HRES for addressing the sensitivity of ice-sheet dynamic processes to bed elevation resolution will be addressed in a subsequent numerical modelling study.The results of the modelling study will also emphasise regions where high-resolution bed elevation data are needed, which will facilitate targeted efforts in data collection.The Cholesky decomposition method used to simulate HRES may be extended to isotropic fields in other areas of research where observations are sparse, such as in the mapping of bathymetry in oceanographic studies or of roughness in the topography under ice shelves.

Figure 1 .
Figure 1.Locations of ICECAP/BC1 bed elevation data included in the synthesis of the Cholesky decomposition roughness terrain (CDRT; Sect.2.1.2).Data are coloured by the natural log of the amplitude coefficient in the covariance data fit, namely log A (m 2 ) in Eq. (3).Numbers 1-27 correspond to drainage basins defined by the Goddard Ice Altimetry Group from ICESat data(Zwally et al., 2012).

FFigure 2 .
Figure 2. (a) Roughness values for each of the ICECAP/BC1 compilations (x axis) and the corresponding overlay points in the Cholesky decomposition roughness terrain (CDRT; y axis) calculated from Eq. (4).The fitted line is calculated using linear least squares.(b) Binned distribution of bed elevation points from the ICECAP/BC1 compilations and the corresponding overlay points in HRES.(c) Cumulative probability density function of bed elevations.
3. HRES bed elevations range from −8848 to 4008 m: within 25 and 10 % of the corresponding bounds in Bedmap2, which are −7054 and 3972 m, respectively.The very low bed elevations in both HRES and Bedmap2 are in the deep ocean.The low-frequency components of the two datasets are es-Earth Syst.Sci.Data, 9, 267-279, 2017www.earth-syst-sci-data.net/9/267/2017/ sentially identical; thus, the difference between them (D = Bedmap2 − HRES; Fig.3c) is essentially a measure of the CDRT roughness introduced in HRES.

Figure 5 .
Figure 5. (a) Locations of selected flight lines from the ICECAP/BC1 compilations.Tracks are read from cyan (start) to purple (end) to provide reference for data in panel (b).(b) Along-track bed elevations from ICECAP/BC1 (black) and corresponding overlay points from Bedmap2 (blue) and HRES (green).The x axis shows the along-track distance (km) from the first point in the flight line.

Table 2 .
Along-track roughness (m) from the ICECAP/BC1 flight lines in Fig.5and the corresponding roughness values from the HRES and Bedmap2 overlay points.RMSE (m) between the ICECAP/BC1 data and the corresponding HRES and Bedmap2 data was normalised by the square root of the number of points in each track.For the ICECAP flight lines, the unique PST (project-season-track) identifier is reported; for the BC1 flight lines, the mission number is reported.