Supplement of Nunataks as barriers to ice flow: implications for palaeo ice sheet reconstructions

Abstract. Numerical models predict that discharge from the polar ice sheets will become the largest contributor to sea-level rise over the coming centuries. However, the predicted amount of ice discharge and associated thinning depends on how well ice sheet models reproduce glaciological processes, such as ice flow in regions of large topographic relief, where ice flows around bedrock summits (i.e. nunataks) and through outlet glaciers. The ability of ice sheet models to capture long-term ice loss is best tested by comparing model simulations against geological data. A benchmark for such models is ice surface elevation change, which has been constrained empirically at nunataks and along margins of outlet glaciers using cosmogenic exposure dating. However, the usefulness of this approach in quantifying ice sheet thinning relies on how well such records represent changes in regional ice surface elevation. Here we examine how ice surface elevations respond to the presence of strong topographic relief that acts as an obstacle by modelling ice flow around and between idealised nunataks during periods of imposed ice sheet thinning. We find that, for realistic Antarctic conditions, a single nunatak can exert an impact on ice thickness over 20 km away from its summit, with its most prominent effect being a local increase (decrease) of the ice surface elevation of hundreds of metres upstream (downstream) of the obstacle. A direct consequence of this differential surface response for cosmogenic exposure dating is a delay in the time of bedrock exposure upstream relative to downstream of a nunatak summit. A nunatak elongated transversely to ice flow is able to increase ice retention and therefore impose steeper ice surface gradients, while efficient ice drainage through outlet glaciers produces gentler gradients. Such differences, however, are not typically captured by continent-wide ice sheet models due to their coarse grid resolutions. Their inability to capture site-specific surface elevation changes appears to be a key reason for the observed mismatches between the timing of ice-free conditions from cosmogenic exposure dating and model simulations. We conclude that a model grid refinement over complex topography and information about sample position relative to ice flow near the nunatak are necessary to improve data–model comparisons of ice surface elevation and therefore the ability of models to simulate ice discharge in regions of large topographic relief.



Compilation of sample position relative to ice flow
For Figs. 2 (main text), S12, and S13, we compile 10 Be and 14 C exposure ages in Antarctica, and show how their position relative to the nunatak summit differs from the direction of ice flow at the nunatak surroundings. Boulder samples containing 10 Be ages were extracted from the expage compilation, by Jakob Heyman (https://expage.github.io/data.html last accessed: March 19 th , 2021), considering the "published" ages. Samples containing 14 C ages were extracted from the Ice-D database 10 (http://antarctica.ice-d.org/ last accessed: April 30 th , 2021), considering those calculated with the LSDn scaling method. The 14 C ages dataset includes the following publications: Balco et al. (2016), Johnson et al. (2019Johnson et al. ( , 2020, Goehring et al. (2019), . For both datasets, a sample's position vector (P) was calculated relative to the closest nunatak summit. The summits were extracted from BedMachine Antarctica (Morlighem et al., 2020) by calculating the Topographic Position Index (TPI, 15 Wilson et al, 2007), which yields only the most prominent features, removes noise, and allows a faster computation time of the raster calculations. A morphological feature map (Wood, 1996) was created from the TPI map, considering slopes steeper than 5 o within a processing window of 5 grid points. The orientation of P was compared to the orientation of ice-flow in the closest vicinity of the nunatak (F). The ice-flow vectors were extracted from the Measures 20 km dataset in Quantarctica v3.2 (Matsuoka 20 et al., 2021), and they indicate regional ice flow directions, rather than local. We use polar histograms to plot the offset between P and F (i.e., O = P -F).
In this compilation, we refrained from filtering for sample distance from the summit. We first considered samples within 30 km of the summit, but found that the number of statistically significant 25 results was considerably smaller. We divided the data into four sectors that span a 90 o arc each: downstream face, upstream face, and two flanks. We tested the use of narrower sectors (45 o as opposed to 90 o ) for classifying a sample, but this also yielded too few samples. To minimise the effect of inheritance, we restricted 10 Be and 14 C samples to those with exposure ages younger than 21 ka. In the 14 C dataset, samples with multiple measurements have their age and uncertainties represented as the 30 mean of all measurements.

Compilation of the nunatak dataset
Thirty-three surface elevation profiles across nunataks parallel to ice flow were analysed to constrain 35 the obstacle dimensions used in our experiments (Fig. S1). This dataset includes examples from different areas of the Antarctic continent (Fig. S2), and were selected so that they would be at least 50 km distanced from other major obstacles in the direction of ice flow flow. Most examples were taken from Dronning Maud Land due to the higher relief and nunatak availability. To determine the axis lengths of the Gaussian surface used to represent nunataks, nunatak axis lengths were determined by 40 fitting an ellipse around the exposed part of the nunatak, and taking the lengths of its major and minor axes. Similarly, outlet-glacier widths were compiled by analysing the spacing between adjacent nunataks in nunatak ranges. Figure S3 shows the frequency distribution for major and minor axes, and the distribution of outlet-glacier width. Our major axis size lies at the upper end of the most common nunataks observed, and the remaining values used for these three variables in our model experiments 45 cover 80 -90 % of the typical sizes. Finally, we consider the three-nunatak 0 km-spacing experiment to satisfactorily cover the largest nunatak cases.

Sensitivity of the ice surface evolution to end-member parameter choices
In order to assess the robustness of our results relative to the choice of model parameters, namely C, 75 which controls basal sliding in Weertman's flow law, and A, which controls the ice rheology in Glen's flow law, we repeated the 'thw' experiment using other possible choices of model parameters. These were also based on the distribution maps for A and C in Gudmundsson et al. (2019), focusing on areas of strong topographic relief around nunataks. The parameter space covered in these sensitivity tests are T = -5 o C and T = -10 o C for the ice temperature influence on A, and log10(C) = -3.5 and log10(C) = -5.5 80 for basal sliding. We chose the nunatak elongated transverse to flow for these experiments because its shape was shown to yield the strongest ice-surface response compared with the nunatak elongated parallel to flow.
Overall, the setup was rather insensitive to the adopted temperature, since results using -5 and -10 o C 85 yield similar results close to the nunataks (Fig. S4). Temperature played a much stronger role in the position of the grounding line, which is beyond the region of interest of our study. Conversely, the ice sheet configuration was more sensitive to the choice of basal slipperiness (Fig. S5). Choosing a higher sliding parameter yielded a much thinner ice sheet, and surface velocities were twice as high as in the control run. When prescribing higher friction (i.e., less sliding) at the bed, the resulting ice cap was 90 much thicker, and did not have a floating-ice terminus. The increased driving stress caused by the thicker ice resulted in mean (median) surface velocities similar to the control run: 32 (30) ma -1 , compared to 33 (29) ma -1 , but slower maximum velocities at the nunatak flanks (39 compared to 53 ma -1 ).
The choice of parameters used in the experiments presented in the main text (treated here as 'control') 95 were shown to yield a profile that is a better generalisation of the profiles in Fig. S1. In the higher sliding case (Fig. S5a), the ice becomes too thin after the spin up, and consequently would yield an almost ice-free domain after the thinning experiments. In the lower sliding case (Fig. S5b), flow velocity at the nunatak flanks s still within the range of values observed along the nunatak profiles (Fig.  S9), and the absence of a floating terminus would require a higher ablation to perform marine ice sheet 100 instability experiments. While these are still plausible conditions for our setup, the resulting surface slope of 1.7% is almost 50% steeper than the control run (1.3%), which is already at the upper end of the range of values observed for Antarctic nunataks (Fig. S1).

Sensitivity of the time of bedrock surface exposure to model parameters
Certain choices of model parameters are likely to impact the time of bedrock surface exposure obtained in the experiments presented in the main text. We performed two additional experiments with a single nunatak under the 'thw' thinning scenario, but using a minimum thickness of 0.5 and 5 m, and an 125 experiment where the ice rheology is softer around the nunatak (i.e. as soft as an ice at -5 o C; Fig. S6), as a representation of a crevassed region around the obstacle. The differences in the time of exposure are shown in Fig. S7, where the standard one-nunatak experiment under the 'thw' scenario (and 1 m minimum ice thickness) is used as "control". 130 Figure S6. Spatial distribution of the ice rheology factor A for the "crevassed" experiment. The value for A away from the nunatak is that used in all other experiments (log10(A) ~ -8.5, or T= -5 o C). Fig. 6 of the main text, but for the minimum ice thickness experiments where the minimum thicknesses allowed are 0.5 and 5 m, and a "crevassed" experiment, where ice is softer around the nunatak compared with the rest of the domain (see Fig. S6).

145
In order to evaluate the effect of having prescribed a linearly varying SMB profile, as opposed to a more realistic representation, we have performed two additional simulations. In both tests, the SMB profile is smoother at the upstream side of the domain, and decays rapidly downstream, through the following equation (cf. Eq. 2 of the main text): ( , ) = SMB 0 − (SMB 0 − SMB e ) * tan ( |x| Lx ) * sin 2 ( |x| Lx ) + ( ) We prescribe two different values for SMB0, -0.92 ma -1 and -0.78 ma -1 , and call the respective experiments 'margin1' and 'margin2'. While 'margin1' ensures that the integrated SMB at t=0 kyr is the same as in the 'thw' scenario (referred here as 'control'), 'margin2' is designed so that the SMB 155 upstream is close to zero. Figure S8 shows how ice thinning is impacted by this different SMB curve, as in Fig. 5a of the main text. Overall, the smoother profile at the nunatak vicinity caused the surface elevation differences up and downstream to be smaller than in the 'control' run, although the general surface elevation is lower. Thus, the differences in ice surface elevation seem to be more sensitive to the gradient in SMB rather than the absolute values themselves. 160    Relationship between sample elevation above sea level (m a.s.l.) and published 10 Be/calculated 14 C exposure ages (in thousand years ago, ka) for the samples that compose panels a and d, respectively. (c,f) kernel density function of the ages shown in panels b and e, respectively. In these plots, the dashed line shows the median age, and the shading shows the uncertainty interval based on the median uncertainty in 10 Be ages, and median external error in 14 C ages. The median ages for each sample group are also displayed, along with the p value for the two-sided Mann-Whitney U test. Figure S13. Sample location relative to flow and age distribution for 10 Be (upper row) and 14 C (lower row) as in Fig. S12, but considering only samples located at the nunatak flanks