Simulating the internal structure of the Antarctic Ice Sheet-towards a spatio-temporal calibration for ice-sheet modelling

Ice Sheet Models are a powerful tool to project the evolution of the Greenland and Antarctic Ice Sheets, and thus their future contribution to global sea-level changes. Probing the fitness of ice-sheet models to reproduce ongoing and past changes of the Greenland and Antarctic ice cover is a fundamental part of every modelling effort. However, benchmarking ice-sheet model data against real-world observations is a non-trivial process, as observational data comes with spatio-temporal gaps in coverage. Here, we present a new approach to assess the ability of ice-sheet models which makes use of the internal layering 5 of the Antarctic Ice Sheet. We simulate observed isochrone elevations within the Antarctic Ice Sheet via passive Lagrangian tracers, highlighting that a good fit of the model to two dimensional datasets does not guarantee a good match against the three dimensional architecture of the ice-sheet and thus correct evolution over time. We show, that paleoclimate forcing schemes commonly used to drive ice-sheet models work well in the interior of the Antarctic Ice Sheet and especially along ice divides, but fail towards the ice-sheet margin. The comparison to isochronal horizons attempted here reveals, that simple heuristics of 10 basal drag can lead to an overestimation of the vertical interior ice sheet flow especially over subglacial basins. Our modelobservation intercomparison approach opens a new avenue to the improvement and tuning of current ice-sheet models via a more rigid constraint on model parameterisations and climate forcing which will benefit model-based estimates of future and past ice-sheet changes.

. Antarctic topography from BedMachine Antarctica (Morlighem et al., 2020) overlain by radar transects analysed in this study. The enlarged sections highlight three regional clusters which we focus on in the following sections. number of published and available traced and dated Antarctic isochrones is limited at the current stage, a situation which initiatives like AntArchitecture aspire to improve in the coming years. For this intercomparison we make use of three excellent compilations of dated isochrones (Leysinger-Vieli et al., 2011;Cavitte et al., 2016;Winter et al., 2019) which provide about 10, 000 km of analysed radar transects covering the past 38 ka to 170 ka BP (see Table 1). With this wealth of information 95 we are able to not only analyse the performance of our model simulations close to ice core locations but also along and across ice divides spanning the whole distance from Dronning Maud Land to the George IV Coast (see Figure 1). The area connected by the radar-transects covers a wide range in the spectrum of bedrock relief, geothermal heat flux as well as climate regimes and ice-sheet elevations. This allows us to analyse the impact of uncertainties in the paleo mass balance, basal melt, ice dynamics and bedrock on modelled isochrone elevations. Table 1 gives an overwiew of the radar-transects used in this 100 work. The Dome C region is included in all three cited studies and provides an excellent region to test the validity of the climate forcing scheme employed by Sutter et al. (2019). The EDC ice core itself provides tight constraints on the age of the isochrones via the AICC2012 chronology (Veres et al., 2013;Bazin et al., 2013) as well as the accumulation history during the last glacial-interglacial cycle (Cavitte et al., 2018). It is important to note that, as the ISM's climate forcing was modulated  via the EDC deuterium record, we expect the best fit of our computed isochrones to be in the region encompassing this ice 105 core. The areas around EDC, Dome Fuji and EDML are all well sampled by reliable radar observations, therefore providing very accurate estimates of bedrock elevation (Karlsson et al., 2018). This makes it possible to investigate the impact of bedrock uncertainties on the mismatch between observed and computed isochrone elevations. Finally, the transects in the vicinity to faster flowing marine drainage regions (e.g. Wilkes Subglacial Basin) allow for an assessment of model parameterisations of ice flow on the computed isochrone elevation.

Observation of isochrones
Radar surveys of ice-sheets observe numerous internal layers, i.e., changes in the dielectric properties of the ice, which lead to a partial reflection of the downward propagating radar wave back to the surface. Most of these layers are caused by changes in the conductivity content of the ice by deposition of acids from volcanic eruptions at the surface of the ice-sheet. After volcanic eruptions these depositions usually occur over a few years and are rather homogeneous over spatial scales of 10s to 100s of 115 km. The absence of melting on the Antarctic plateau leads to continuous burrial and, thus, submergence of these conductivity layers, which are typically spread over less than one meter (Eisen et al., 2006). Lateral variations in accumulation rates cause differences in submergence velocities and, thus, to deformation of these initially surface-parallel layers. With increasing depth, the influence of ice dynamics on the (vertical) submergence velocity increases and deforms the layers further.
Radar systems used to map internal layers and the thickness of ice-sheets operate in the range of 10s to 100s of MHz, 120 corresponding to a wavelengths on the order of meters. The width of the observed reflections, caused by the interaction of the electromagnetic radar waves with the conductivity peaks, is determined by the bandwidth or pulse length of the radar systems.
They are typically in the range of 0.5 to several tens of meters (Winter et al., 2017).
With increasing depth the internal layers often follow the bed topography in an often smoothed fashion. Thus, the lateral variation of the depth of these internal reflection horizons is usually sufficiently small to be continuously detected with the 125 spatial sampling distance of airborne radar systems. They operate with a horizontal sampling distance in the range of meters 5 https://doi.org/10.5194/tc-2020-349 Preprint. Discussion started: 2 December 2020 c Author(s) 2020. CC BY 4.0 License.
to tens of meters (Winter et al., 2017) and can be traced over hundreds to thousands of kilometres. As each of these layers originates from the same physical change in conductivity, which the ice inherited during deposition as snow at the surface, they are considered isochrones and have the same age along the reflection horizon.
Dating can most reliably be achieved by tying internal layers to ice core sites, where age-depth estimates are available from 130 ice-core proxy analysis. For the comparison of observed and modelled isochrone elevations, the uncertainty of the observed isochrones is an important parameter. Apart from the radar-inherent uncertainties coming along with identifying and continuously tracing the same radar signal for an horizon and the conversion of two-way traveltime to depth, the uncertainty in the ice-core ages and the actual matching procedure between horizons and ice core are the largest sources of error. We use various internal layers in this study. Details on the radar resolution, dating methods and resulting uncertainty can be found in the depth towards the bed. This range covers most of the horizons considered in this study. At places where the deepest ice is younger than around Dome C, the age gradient with depth will be less steep towards the bed. Thus, the age uncertainties of the 140 horizons will be lower. We consider this uncertainty always to be smaller than the proposed age derived from the ISM data.

Large-scale ice-sheet modelling
In order to compute Antarctic isochrones, time resolved 3D velocity data as well as the transient ice-sheet geometry (ice thickness and bedrock topography) are necessary. We therefore ran a paleoclimate model ensemble covering the last 220 ka and make use of previously published model ensemble results covering the last 2 Ma ). An overview over 145 the experimental setup is provided in Figure 2. The model forcing is described in detail in Sutter et al. (2019) and provides spatio-temporal information on ocean temperatures, surface temperature and precipitation derived from climate snapshots during the LIG (Pfeiffer and Lohmann, 2016) and the LGM as anomaly forcing with respect to the modelled pre-industrial climate state. In between the climate snapshots, surface temperature and ocean temperatures are interpolated on the basis of the EDC deuterium data (Jouzel et al., 2007)

Lagrangian tracer advection 160
We employ the open source python software OceanParcels v0.9 (Lange and van Sebille, 2017) which was originally designed to analyse ocean data but can be easily adopted for transient or steady state ISM data. At its core, OceanParcels provides a Lagrangian particle tracking scheme which is computationally efficient and adaptable to a variety of input grids ( Sebille (2017)): where X is position of the particle in space, v the three-dimensional velocity at the particles' position and ∆X b (t) is the change in position due to "behaviour", which does not play a role in our application as we assume that every ice parcel is passively advected with the velocity field without any additional advection terms. For further details regarding OceanParcels please refer to Lange and van Sebille (2017). For OceanParcels to work we need to provide three dimensional velocity data 170 covering the time frame we are interested in (i.e., the last 160 ka as this is the age of the oldest isochrone we take into account) slow ice flow in East Antarctica. We set the advection time step (∆t) to one year and for each experiment we employ 200, 000 individual tracers which we seed at the ice surface at the time of the age of the isochrone we want to compute (e.g. at 38 ka BP for the 38 ka isochrone). Misfits of the ice-sheet model state in terms of elevation and velocity field relative to the true (unknown) ice sheet state at that point in time in the past, will lead to deviations of the modeled isochrone as observed in the ice sheet today. This information can then be used to identify such misfits and improve the model representation of the ice 180 sheet. We find that for less than 10 5 tracers, gaps in coverage due to the dispersive nature of ice flow are too large to resolve the internal layers away from ice divides. The seeding mask (i.e., the region where tracers are initiated) consists of the area with an ice thickness of more than 2200 m and surface flow slower than 7 m/a (see Figure

Modelled Antarctic isochrone elevations
Our main goal is the identification of systematic mismatches between predicted and observed isochrone geometry. To achieve 195 this, we focus on three major sources of model-observation mismatches: (i) climate forcing, (ii) model parameterisation, (iii) bedrock and geothermal heat flux. As computed isochrone elevations will be affected by a superposition of uncertainties from all three sources we separate our analysis into regions where we expect only one factor to dominate. To isolate the effect of (i) climate forcing we turn to areas characterised by limited ice flow (e.g. ice domes or close to ice divides) and by a good representation of the bedrock elevation in the model (difference between input data set and local high resolution radar data).

200
Here, we will primarily focus on Dome C as the model ensemble was tuned to match the regional ice-sheet configuration. To assess the impact of model performance and ice flow parameterisation we discuss the model-fitness to compute isochrones in regions of elevated surface velocities and in subglacial basins (bedrock below sea level) as well as the effect of different parameter choices under the same boundary conditions (climate, geothermal heat flux, bedrock topography). To test (iii) we analyse the impact on bedrock elevation uncertainty around Dome C and use different geothermal heat flux fields to investigate 205 their influence on the modelled internal layer architecture.

Dome C -evaluating the paleoclimate forcing
We now turn to the first (i) of the aforementioned three main factors determining simulated isochrone elevations, the climate forcing. Due to the near cessation of surface flow at ice domes the local internal stratigraphy is by and large dominated by the regional surface mass balance history and potentially affected by basal melting due to a large ice thickness, low surface  Above, we identified which paleo temperature-precipitation scaling leads to the best fit of modelled and observed isochrone elevations. Making use of this paleo-accumulation estimate we can plot the transient accumulation history at EDC during the last glacial interglacial cycle as well as the regional precipitation forcing during the LGM and the LIG. Figure 5 A) RACMO (van Wessem et al., 2014). It is important to note, that our approach does not strive to achieve a detailed reconstruction of paleo-accumulation as done in Cavitte et al. (2018). They are fitting englacial layer horizons via a 1-D pseudo-steady ice flow model  leading to more accurate estimates of the local surface mass balance along the ice divide, while we are merely trying to establish the validity of our paleoclimate forcing approach. Nevertheless, Figure 5 illustrates  the missing imprint of the paleo-spinup. Figure 6 shows the difference in isochronal elevation between pal and pd-pal. It is 265 relatively small in comparison to the difference between pal and pd-pd despite the fact that surface elevations for pd-pal and pd-pd are very similar. This illustrates that tuning an ISM to present day observations and mass balance trends (pd-pd) might yield a suitable fit to ice-sheet surface observations but fails to capture the integrated internal layering which is a product of the paleo surface-mass-balance as well as of paleo ice-dynamics. This might have important repercussions for paleo ice-sheet studies as well as for projections of the future sea-level contribution of the Antarctic Ice Sheet, as the ice-sheet's initial state 270 can affect its future behaviour over centennial or even decadal time scales (Seroussi et al., 2019) due to misrepresentations e.g. in the parameterisation of basal friction, the internal flow fields (and thus mass transport) and the thermal state as a consequence of overfitting. Ideally, to assess the impact of isochrone elevation calibration on ice sheet model projections, one would compare three cases where the model parameters are tuned a) via thickness/surface velocity inversion, b) against present day surface observables and paleo-proxies (ideally including a paleo spinup) and c) isochrone elevation calibration combined with 275 b). This is work in progress and beyond the scope of this manuscript.

Caveats to modelling isochrones with large scale ice-sheet models
Differences between modelled isochrones (based on paleo-velocity fields) and their observed elevation above bedrock are small for all isochrones sampled in this study which are younger than the LIG (< 120 ka, root mean square deviation of RMSD < 5%). It is striking however, that for isochrones older than 130 ka the gap between model results and observation widens (see 280 Figure 7). Due to the lack of climate model data for glacials and interglacials preceding the LIG, we estimated the climate forcing before 130 ka BP to be a linear combination between the climate state of the LGM and the LIG using the EPICA Dome C deuterium record as an index to interpolate between the two climate states (see section 2.1 in Sutter et al. (2019)).
Naturally, this approach is just an approximation of the actual transient climate conditions and therefore will necessarily lead to discrepancies in the surface mass balance as well as the ocean forcing. Another potential cause of the relatively poor 285 representation of isochrones before 130 ka BP is the high circum-antarctic subsurface ocean temperature peak assumed for isochrone elevations was only 60 and 90 ka, respectively. This could be another aspect influencing the particularly poor match between modelled and observed isochrone elevations for these ages.
Next to the uncertainties associated with paleoclimate forcing fields and ice-sheet model parameterisations the uncertain-295 ties in the basal boundary conditions (iii) applied in a large-scale paleo-ice-sheet model further complicate the computation of isochrone elevations. Uncertainties in the bedrock elevation in areas with only sparse radar observations can feature differences of several hundred meters if compared to newly acquired high-resolution radar data (Karlsson et al., 2018;Morlighem et al., 2020). This can affect the basal flow regime and corresponding thermal state of the local ice column significantly. We would expect, that line transects characterised by a bedrock profile in relatively good agreement with reality should yield a good 300 agreement between modelled and observed isochrones as well. While this is true e.g. for some parts of transect DC-X57 in the Dome C region (see Figure 8), the assumption cannot be generalised to all East Antarctic transects discussed in this work as uncertainties in the parameterisations of ice flow, paleoclimate forcing and geothermal heat flux can have a dominating or confounding effect over differences in the observed and modelled bedrock relief. Additionally, the impact of bedrock undulations on isochrone elevation will be stronger close to bedrock and might not be observable in relatively young isochrones where 305 the surface mass balance is the dominating factor. A further caveat is the largely unknown distribution of geothermal heat at the base of the ice which is only inferred indirectly and exhibits large differences in between data-sets Shapiro and Ritzwoller  The difference between the model bedrock and the 1 km BedMachine data (Morlighem et al., 2020) is highlighted by the grey shading. The labeled circles denote the surface flow and magnitude in m/a relative to the direction of the transect.
The coloured bars on the bottom of depict the geothermal heat flux and the computed basal melt rate along the transect as in Figure 6.
surface boundary conditions as well as model parameterisations on East Antarctic isochrone elevations across the whole ice column is beyond the scope of this study but could be the focus of an upcoming paleo-ice-sheet model intercomparison.  of the whole ice-sheet. The focus on the EDC ice core region in the previous section allowed us to appreciate the impact of uncertain forcing fields and bedrock elevation on a relatively small length scale (100 km). We now turn to larger distances, 315 where the Dome C-derived transient paleoclimate forcing and model optimisation might lead to more pronounced diversions between the modelled isochrone elevations and the radar data. We first assess the 90 ka isochrone along transect CEA-10 (see Figure 9). It starts close to TALDICE crossing the upper reaches of the Wilkes Subglacial Basin after which it traverses the Concordia Subglacial Trench, passing EDC and LDC and finally arching upwards to end at Lake Vostok (total distance ca.
1500 km). 320 Figure 9 illustrates nicely, that there are two main regimes covered by profile CEA-10. The northern half of the transect bottom part in Figure 9) is dominated by the imprint of the Wilkes Subglacial Basin with a large misfit between the modelled and observed isochrone. The striking excursions of the modelled isochrone elevation in Figure 9 point towards a too fast internal vertical ice motion. This could possibly be due to the heuristics involved in determining the yield stress in subglacial basins. While a spatially varying temperature-precipitation relationship could also affect the isochrone elevation in this region, 325 this seems not to be the case here. We show this by computing the 90-ka isochrone elevation in a paleo-simulation with a temperature-precipitation relationship of 8%K −1 which cannot mitigate the drop in elevation in the first 300 km of the transect (see Figure 9). We did not test temperature precipitation relationships higher than 8%K −1 , as they would be beyond what is reconstructed from East Antarctic ice cores. Geothermal heat flux can also be ruled out from being the culprit, as for the most part the thermal regime in the data set from Shapiro and Ritzwoller (2004) at Dome C is similar or even lower (of course 330 the actual Geothermal heat flux could be much higher affecting local isochrone elevations, this could be tested in a more