Review of Mas e Braga and others, "Nunataks as barriers to ice flow: implications for palaeo ice-sheet reconstructions"

In this manuscript, the authors apply the ice sheet model Ua to an idealized, rectangular geometry including Gaussian nunataks, and investigate how the ice flow and geometry react to these barriers. They report that, for conditions mimicking those typical for Antarctica, nunataks can influence the geometry significantly, such that the ice thickens upstream and thins downstream of the obstacle. Details vary depending on scenarios and mesh resolution. The results are relevant because they point out that the interpretation of cosmogenic exposure dating in the vicinity of nunataks in terms of past ice-thickness changes is far from being straightforward and can be grossly wrong if the flow disturbance caused by the nunataks is ignored.


Introduction
Ongoing changes in climate are already causing significant mass loss and ice-margin retreat of both the Antarctic and Greenland ice sheets (Garbe et al., 2020;King et al., 2020). Near-future (2100 CE) projections of sea level rise point to ocean thermal expansion as the main cause, but the wide uncertainty (0.43-2 m) is attributed to uncertainties in ice mass loss from the 25 Antarctic Ice Sheet (Oppenheimer et al., 2019). Over multi-centennial timescales, sea level contribution from Antarctica is expected to become dominant , especially through ice loss from outlets with retrograde bed slopes . Numerical ice sheet modelling efforts aimed at reducing uncertainty by better understanding the processes that lead to sea level rise focus on both the near future (Goelzer et al., 2020;Seroussi et al., 2020), and inferences from the geological past (Pollard and DeConto, 2009;Albrecht et al., 2020). Recent efforts include improvements in some key model 30 components such as grounding line dynamics (e.g. Gladstone et al., 2017;Seroussi and Morlighem, 2018), model coupling to solid Earth and sea level models (e.g. Gomez et al., 2020), and improved treatment of ice-ocean interaction processes (e.g. Reese et al., 2018;Kreuzer et al., 2020). The importance of bedrock topography  and grid resolution (Durand et al., 2011) have been acknowledged previously, and studied particularly for marginal regions of the ice sheet (e.g. Sun et al., 2014;Robel et al., 2016;Favier et al., 2016). 35 Bedrock topography plays a strong role in regulating mass loss of ice sheets. Topography that deepens towards an ice sheet interior can accelerate ice sheet retreat due to a positive feedback on ice flux across the grounding line -referred to as Marine Ice Sheet Instability (Schoof, 2007). Conversely, spatial variations of bedrock topography (Robel et al., In Review) and the resulting drag exerted at the ice base and on the sides of ice streams (Jamieson et al., 2012(Jamieson et al., , 2014, outlet glaciers (Jones et al., 2021), and fjords (Åkesson et al., 2018), can slow down or even stabilise grounding line retreat. When evaluating ice loss 40 beyond this century, inland regions also become important, especially those with large subglacial topographic relief, such as the overridden mountain ranges that fringe the glaciated cratons of Greenland and Antarctica (Howat et al., 2014;Burton-Johnson et al., 2016).
The accuracy of ice sheet models is limited by grid resolution (e.g. Cuzzone et al., 2019), simplifications in model physics (Cuffey and Paterson, 2010;Hindmarsh, 2004) and uncertainties in the climate forcing (e.g. Seguinot et al., 2014;Alder and  predominant setting from where samples are acquired (e.g. Ackert et al., 1999;Lilly et al., 2010;Bentley et al., 2014;Suganuma et al., 2014;Jones et al., 2017). Cosmogenic nuclide exposure ages are determined by the concentrations of cosmogenic nuclides in erratic boulders or cobbles (i.e. glacially transported clasts of a different lithology than the underlying bedrock, deposited during periods of ice cover) or bedrock surfaces. The concentration of a cosmogenic nuclide increases the longer a rock surface is exposed to cosmic rays from space (Gosse and Phillips, 2001). Assuming no cosmogenic nuclides have been 60 inherited from a prior period of exposure, and by means of known production rates, the cosmogenic nuclide concentration allows inference of the last time the ice sheet covered that specific location, and consequently the ice surface elevation at the time of sample exposure. Comparing the sample's elevation to the present-day ice-surface elevation then enables determination of the magnitude of ice thinning between the present and the time of exposure. Yet, the gradient in ice surface elevation up and downstream of a nunatak range often does not have a clear established relationship with the regional ice surface, a deficiency a.s.l.) in both types of experiments, which is representative of continental shelf depths of Antarctica .
The grounding line in Úa is defined by the flotation condition, and although its migration is not the focus of our study, we allow it to move freely, refining the mesh elements around the grounding line as it evolves (as described at the end of this section).
Nunataks in this domain are represented as three-dimensional Gaussian surfaces, which are superimposed on B with their centre (i.e. the nunatak summit) at |x| = 50 km and y = 0 km. All generated nunataks have an outcrop size (i.e. exposed area above the ice surface) after spin up of approximately 12 × 5 km along its main axes. Depending on the experiment, the nunatak 120 was elongated transverse to flow (i.e. along the y axis) or parallel to flow (i.e. along the x axis). The adopted nunatak dimensions and the value of β are within the interval constrained by 33 profiles across nunataks in Antarctica (Figs S1-S3).
The effects of ice rheology (including temperature) in Úa are accounted for by the strain-rate factor A in Glen's flow law. To compute A, which we treat as spatially uniform and constant over time for the purpose of our experiments, we assume an ice temperature (T) of −20 • C. This yields a value similar to that found for regions surrounding nunatak escarpments in Antarctica, 125 as obtained in Gudmundsson et al. (2019) when inverting for A and basal slipperiness (C) based on satellite-derived ice surface velocities. Following the same reasoning for A, we assume basal sliding to follow Weertman's sliding law (Weertman, 1957), and use a constant value for C of log 10 (C) ≈ −4.5 (C = 2.9 · 10 −5 m kPa −1 a −1 ). The main model parameters are summarised in Table 1, and sensitivity analyses of values for A and C are presented in the supplementary material (Figs. S4, S5).
The surface mass balance (SMB; in metres per year, ma −1 ) parameterisation applied to the spin up and subsequent experi-130 ments is given by Eq. (2): where SM B 0 is the SMB at the ice divide (x = 0), SM B e is the SMB at the edge of the domain (x = 150 km), Lx = 150 km, and b(t) is a time-dependent SMB term, through which perturbations to the total SMB are applied. We prescribe SM B 0 = 1.3 ma −1 , and SM B e = −0.3 ma −1 , which results in no ablation occurring over our region of interest (black line in Fig. 3a), 135 as usual for Antarctic settings.
An unstructured mesh was generated for the domain, which is refined during simulation time (including during spin up) based on a series of glaciological refinement criteria. Element size is refined according to ice thickness, from a maximum of 8 km down to 205 m around the interface between ice and nunatak, where ice thickness approaches the minimum (which we set as 1 m). The mesh at the grounding line is refined to 500 m within a 2 km buffer zone on both sides. The mean element size 140 after spin up was 740 m (414 m median).

Model spin up
Model spin up starts from a 1200 m-thick uniform ice distribution, to which a constant (but spatially variable) SMB is applied (i.e. b(t) = 0 ma −1 in Eq. 2) for 20 thousand years (kyr). This period is long enough for the system to reach equilibrium with the SMB forcing. The ice surface slope from the divide to the grounding line, which after spin up was located at x = 136 km 145 (Fig. 3b), is ∼ 1.3 %. This inclination is representative of measured profiles along nunataks for various regions of the Antarctic Ice-divide bedrock elevation (B0) 750 and −750 m a.s.l.

Ice thinning experiments
In order to understand whether ice thinning occurs uniformly up and downstream of a nunatak, we impose three different degrees of ice thinning uniformly over the domain. The perturbations to SMB that induce ice thinning are applied through b(t) in Eq.
(2). The evolution of b(t) is based on a smoothed step curve that applies a weight, evolving from 8 % to 100 %, mimicking the deglaciation progression recorded in ocean sediment and ice cores for the past 20 kyr (e.g. Lisiecki and 160 Raymo, 2005;Jouzel et al., 2007). This "weight curve" is then multiplied by a constant chosen in order to match Last Glacial Maximum-to-present ice thinning at three contrasting regions as inferred from a continent-wide transient modelling experiment . These regions are Thwaites Glacier, in West Antarctica ('thw', −0.75 ma −1 ), and the Transantarctic Mountains ('tam', −0.60 ma −1 ) and Dronning Maud Land ('dml', −0.45 ma −1 ) in East Antarctica (Fig. 3c). These three different thinning scenarios were applied to nunataks that were elongated along and transverse to ice flow, and to reference 165 experiments without a nunatak for comparison purposes.

Three-nunatak experiments
Across much of Antarctica and Greenland, nunataks are more common in groups than in isolated cases. The aim of this set of experiments is therefore to test whether three nunataks, separated by narrow glaciers, yield ice surface elevations and gradients that differ from the control run due to a combined effect of all nunataks on ice flow. For this test, we perform sets of 170 experiments with the same SMB as applied in the 'thw' experiment with one nunatak (control run), but using three nunataks

Mesh-resolution experiments 175
In a final series of sensitivity experiments, we assess how well regular-spaced grids of coarser resolutions typically used in ice sheet models (5, 10, and 20 km) resolve the ice surface elevation pattern around nunataks compared to the solution using an unstructured, locally refined mesh. We do so by repeating the full set of three-nunatak experiments and the one-nunatak 'thw' scenario (as control) for these regular meshes (i.e. without refinement).
180 All sets of experiments, their respective surface mass balance (SMB) perturbations, and the number of nunataks in each set are summarised in Table 2. 3 Results

Ice thinning experiments 185
Our experiments clearly demonstrate that the presence of a nunatak impacts ice surface elevations, and that the response magnitude depends on its orientation relative to ice flow. Using the 'thw' scenario as an example, the ice surface directly upstream of the nunatak is higher than its surroundings for the entire simulation period, while downstream it is lower (Figs. 4a,b). In the experiment where the nunatak is elongated transverse to flow (Fig. 4a), this effect extends up to 30 km upstream and downstream along the centreline, and 15 km perpendicular to the centreline (of which the first 6 km are exposed bedrock). At the 190 start of the simulation (i.e. 0 kyr after spin up), the ice sheet surface intersects the nunatak on its upstream face 360 m above the general ice sheet surface 15 km away from it (transverse to the centreline). Downstream, the ice sheet surface intersects the nunatak 360 m above the lowest ice surface elevation recorded in the nunatak vicinity (located 3 km downstream of the summit). The lowest ice surface elevation is also 100 m lower than the ice surface elevation 15 km away from the centreline.
Along the nunatak flanks, the ice surface elevation falls about 300 m. In summary, while the general ice surface slope over 195 the grounded part of the domain is 1.3 %, the slope is about 2-3 times steeper around the nunataks. For experiments where the nunatak is elongated parallel to ice flow, a similar, but less pronounced response is observed, including a lower ice surface elevation gradient, and relative elevation-change effects that can be discerned 20 km up and downstream along the centreline, and up to 5 km transverse to the centreline (Fig. 4b). In both cases, relative changes in ice surface elevation across the centreline due to the nunatak presence can be discerned from the general ice surface elevation until where the subglacial continuation of 200 the nunatak can be distinguished from the smoothly sloping bedrock (Fig. 4a,b, and insets).
Both ice thinning experiments detailed in Fig. 4a,b reveal an overall steepening of the ice surface gradient during the 20 kyr thinning period (profiles from red to green). The steepening of the profile through time is illustrated in Fig. 4c,d, where points equidistant from the nunatak summit (upstream and downstream) indicate a steepening profile 30, 10 and 3 km from the summit for nunataks elongated transverse to ice flow, and 30, 10 and 8 km from the summit for nunataks elongated parallel to ice flow, 205 respectively. For equidistant locations closer to the nunatak summit (< 2.5 km and < 7.5 km for Figs. 4c,d, respectively) this pattern becomes disrupted by the ice behaviour around the nunatak. This is because, relative to the case without a nunatak (the grey profile that accompanies each colored profile set), the ice surface profile steepening over 20 kyr increases with the introduction of the nunatak. The ice surface profiles are the steepest in immediate conjunction with the nunatak (see differences between solid lines and dashed and dotted lines), which has important potential implications for bedrock exposure patterns.

210
Indeed, for locations 7.5 km and 2.5 km away from the nunatak summits in Figs. 4c,d respectively, the initial pattern of steepening is reversed as the downstream face of the nunatak becomes ice free (and so elevation does not further decline) whereas the upstream ice surface slowly declines throughout the deglaciation experiment. This pattern of increased steepening around nunataks is potentially important for cosmogenic nuclide studies to consider because from these modeling experiments, significant differences in exposure of up and downstream faces of nunataks become apparent (Fig. 5).

215
To analyse the difference in timing of ice surface evolution and nunatak surface exposure up and downstream of the nunatak, we select five pairs of points equidistant from the nunatak summit. They span a distance between 1.5 km, where the downstream side is already exposed at the start of the experiment, and 2.5 km, the closest element to the nunatak where the ice continues to thin normally (i.e. it does not become exposed or stops thinning) in any ice-thinning scenario. The chosen points (Fig. 5) are also separated in the model mesh by one element size from one another. In our analysis, we consider that nunatak surface 220 exposure commences when ice thickness falls below 10 m. A thickness threshold larger than the minimum thickness in the model is used for a consistent identification of the timing of surface exposure between the different experiments and different points analysed, and 10 m yielded the best results among the thresholds tested. For the purpose of cosmogenic exposure dating, most cosmogenic nuclide production occurs when ice is thinner than~10 m (Gosse and Phillips, 2001). At the upstream side, if the previous criterion fails, we then determine when the ice surface stabilises and thus reaches its minimum elevation.

225
The implications of comparing time of exposure (thickness < 10 m) with time of stabilisation (thinning < 5 cm/100 yrs) are discussed further in Sect. 4.1. For visualisation purposes, only the experiments with nunataks elongated transverse to ice flow are shown in Fig. 5, but the patterns shown also hold true for nunataks elongated parallel to ice flow. Because of the elongation of the nunatak, its initial area exposed along the centreline is larger, and thus the pair of points compared for the latter are located further from the nunatak summit (4.5, 4.8, 5.0, 5.5, and 6.0 km). only the downstream side is exposed, while its upstream counterpart is still thinning at the end of the simulation time. Within 2 km from the nunatak summit (or 6 km for the nunatak elongated parallel to flow), and for those scenarios where both up 235 and downstream sides are exposed (or meet the stabilisation criterion), the upstream side lags its downstream counterpart by 2 to 14 kyr (Fig. 5b). In a retrograde-slope setting, however, this effect is not observed. Rapid and uniform thinning happens once accelerated grounding line retreat is triggered, akin to Marine Ice Sheet Instability, yielding a similar adjustment up and downstream of the nunatak (Fig. S7). for nunataks elongated transverse to ice flow under the different thinning scenarios (see Table 2

Three-nunatak experiments 240
The experiments with three nunataks show how different glacier widths produce dissimilar responses of ice surface elevation ( Fig. 6), mimmicking situations where multiple nunataks are separated by narrow glaciers. After 20 kyr of ice thinning, the 0 km-width experiment (i.e. where the nunataks merge into a single, 36 km-wide nunatak; Fig. 6a) yields the highest ice retention upstream (further delaying surface lowering at the upstream side), with its surface up to 250 m higher compared to the control one-nunatak case (Fig. 6e). For the experiments with the smallest glacier widths (0, 5 km; Figs. 6a, b), the 245 constructive interference between the nunataks results in an ice retention that is strong enough to instigate a strong deficit response downstream. Compared to the control scenario, this results in an even lower ice surface elevation up to 50 km  The differences in glacier width also result in different organisations of ice flow around the nunataks. The wider barriers formed by the 0 and 5 km-wide glaciers yield a different pattern of ice flux (Fig. 7), where the peak is much stronger and concentrated 13-30 km away from the centreline. In the case of the two experiments that have wider glaciers (10 and 15 km), ice flux also peaks further away from the centreline than in the control experiment, but is more uniformly distributed across the domain. Their largest peak values are closer to those of the control experiment, and distinct smaller peaks occur exactly where 255 the glaciers are located. The difference between the flux downstream and upstream is inversely proportional to glacier width (i.e. the narrower the glacier, the larger the difference), which points to an increased retention of ice upstream as the cause for the increase in relative heightening/lowering of the ice surface.

Mesh-resolution experiments
The coarser regular-mesh experiments (5, 10, and 20 km horizontal resolution) applied to the series of three-nunatak configura-260 tions show important differences from the refined-mesh experiments (Fig. 8, cf. Fig. 6). The 5 km-mesh experiments deviate the least from the refined-mesh experiments, and capture the relationship between glacier width and ice surface elevation gradient.
Still, a 5 km mesh does not capture the same magnitude of ice surface lowering downstream, underestimating it by as much as 80 % (100-200 m) relative to the refined mesh, an effect that is particularly pronounced close to the nunataks (cf. Figs. 8a,b and Figs. 6a,b).

265
The 10 km mesh captures, to some degree, the original nunatak shape and consequent surface elevation increases upstream and decreases downstream of the nunataks (Figs. 8f-j). However, there are problems with both signal strength (under/overestimations of at most 250 m) and location (not coinciding with the pattern seen in the refined-mesh experiments). In the 20 km mesh, the original nunatak shape is not properly captured (Fig. 8k-o), and the large element sizes cause increased heightening/lowering of ice surface elevation to occur much farther away up and downstream of the obstacles compared to their 270 respective refined-mesh counterparts (Fig. 6). Finally, the 20 km mesh yields similar results for the 0, 5, and 10 km glacierwidth experiments (Figs. 8k-m), but yields a much different response for the 15 km-width experiment (Fig. 8n), which is much more in tune with the results of the refined mesh.
The differences between the regular-mesh and the refined-mesh experiments also result in different thinning rates between them. For example, in the control simulations used as reference for the elevation differences shown in Figs. 6 and 8, a point 30 275 km upstream (downstream) of the nunatak in the refined mesh thins to 1500 (750) m a.s.l. after 13.3 (11.2) kyr ( Fig. 9). In the 5 km mesh, the same pair of points thins to these levels after 15.2 kyr and 11.9 kyr, respectively, and in the 20 km mesh, after 7.2 kyr and 10.8 kyr respectively. In other words, the 20 km resolution model run overestimates the rate of thinning under the same forcing (even though, in this particular case, the point 30 km downstream already starts from a lower elevation; Fig. 9c), while the 5 km resolution model run shows much closer values to the refined-mesh experiment, despite the underestimated thinning.

280
The finest regular mesh tested (5 km) performs the best among the regular meshes, since it partially captures the widest glaciers tested (15 km wide), and best represents the effect of glacier width on ice flow constriction (Fig. 10). The differences between ice flow downstream and upstream shown in Fig. 10 evolve similarly in the refined and 5 km meshes, but are not well represented in the 10 or 20 km meshes. The coarsest mesh (20 km) shows similar results between the 5 and 15 km-wide glacier experiments, and between the wider nunatak and the 10 km-wide glacier experiments. This is because the coarse mesh 285 resolution misses two nunataks summits, which lie between two nodes. As a result, the lower topography that is captured by this coarse mesh becomes a subglacial continuation of the single central nunatak. This also explains why the 15 km-with experiment showed an ice surface elevation pattern similar to the control run. Further tests indicate that the 20 km mesh only Figure 8. Difference in surface elevation (m, as in Fig. 6) between each 'thw' three-nunatak configuration (varying glacier widths) and the 'thw' control with one nunatak (right-most panels), for three regular mesh resolutions of 5, 10, and 20 km. The reference experiments were also performed on a regular mesh. These figures illustrate how differential elevation changes up and downstream of the nunataks (digitised in white based on their outcrop size, linearly interpolating the ice surface between the nodes and vertices) differ from the experiments using a refined mesh (Fig. 6), and how this mismatch decreases with increasing glacier width.
captures the existence of three nunatak summits when they are spaced by glaciers that are at least 20 km wide (not shown).
Still, in these tests the glaciers are not wide enough for the mesh to capture their existence, and thus the results only reflect the 290 effect of a wider obstacle.   Minchew et al., 2018). However, our sensitivity tests show that the surface response around the nunatak is insensitive to variations in the rheology factor within typical values found for our target regions in Antarctica (Fig. S4). The experiments indicate that, across ice flow (i.e. along the y direction), the extent of the ice surface that is impacted is more likely related to the 305 interaction of ice flow with the subglacial extension of the nunatak, as commonly reported for different spatial and temporal scales, and modelling setups of different complexities (e.g. Siegert et al., 2005;Durand et al., 2011;Cuzzone et al., 2019;Paxman et al., 2020).
The ice surface steepening and consequent mismatch between the up and downstream sides increases as the ice thins, up to the point when the downstream side becomes exposed. Exposure happens earlier downstream, as expected due to lower ice 310 surface elevation, and an equidistant point upstream becomes exposed (or has its thinning stabilised) up to 14 kyr later than its downstream counterpart. While we used stabilisation of thinning upstream to determine whether the ice surface attained its minimum elevation before reaching the minimum thickness criterion, such stabilisation does not happen because of an equilibrium of the modelled upstream surface with the applied SMB. If the upstream surface had attained equilibrium, stabilisation would have happened earlier when the imposed thinning was lower, which is the opposite of what was observed when com-315 paring the three thinning scenarios. The fact that equidistant points up and downstream of the nunatak were in some cases not both exposed also implies that important changes in ice surface elevation might not be recorded upstream of a nunatak. The difference in time of exposure up and downstream is also higher in the experiments when ice flow is further constricted by the narrow glaciers (Fig. S8).
Although our setup is largely based on Antarctic settings, the idealised SMB forcing deviates from what is more commonly 320 observed in Antarctica in two ways, which could have influenced our results. First, we prescribe ice thinning through surface melting, which is more characteristic of Greenland (e.g. Kjeldsen et al., 2015), while Antarctica's main source of mass loss is through dynamic thinning (Pritchard et al., 2009). Since mass loss occurs only downstream of the nunataks and is highest at the downstream end of the model domain, we believe that the interaction between ice flow and nunataks, and the consequent ice surface elevation response, were not significantly impacted by how ice thinning was imposed. Second, the temporal pertur-325 and how the differential response between the upstream and downstream surfaces introduce a bias on the timing of bedrock surface exposure around the nunatak. These effects are important for the interpretation of cosmogenic-nuclide ages and for comparing such ages with results from ice sheet models, and are discussed in more detail next.

Implications for the interpretation of past ice sheet reconstructions
The steepening of the ice surface around nunataks has important implications for the interpretation of in-situ constraints on 335 past ice thickness changes from surface exposure dating. A commonly adopted practice is to assume that regional ice surface elevations are directly reflected by the absolute elevations of samples, but this is likely to yield inaccurate results for past ice sheet reconstructions. Our study demonstrates that this assumption would often yield an error of up to almost 400 m, which is of the same magnitude as many reported thickness-change estimates in regions of significant ice surface relief (e.g. Ackert et al., 2007;Suganuma et al., 2014;Kawamata et al., 2020). A different practice, which allows for improved comparisons between 340 sites, uses sample elevations relative to the modern ice surface elevation to infer past ice sheet thickness changes (e.g. Johnson et al., 2008;Jones et al., 2015). A key assumption in this approach is that the ice surface gradient and organisation of ice flow around a nunatak remained the same through time. This assumption contradicts our modelling results, which show increasing ice surface gradients around nunataks during ice sheet thinning. Our experiments further indicate that when samples are taken upstream (downstream) of a nunatak, estimates of past regional ice surface elevation will be overestimated (underestimated). In 345 directions transverse to ice flow, the ice surface is typically at a lower (higher) elevation than directly upstream (downstream) of the nunatak summit. The interpretation of thickness evolution is further complicated considering that the direction of ice flow could have changed as the ice thinned (e.g. Suganuma et al., 2014;Fogwill et al., 2014). This means that the present-day ice surface used as reference elevation could have behaved differently at the time of sample exposure, resulting in an elevation mismatch different from that of the modern ice sheet.

350
An alternative practice for inferring ice thickness changes is to determine minimum and maximum estimates. When constraining the ice sheet thickening necessary for the ice surface to reach sample elevations in Dronning Maud Land, Andersen et al. (2020) use two approaches to determine their reference present-day surface elevation. For a minimum estimate of required thickening, the reference point for ice surface elevation is chosen where there is a major break-in-slope (treating as the "regional" reference) and the point at a local depression (treating as the "local" reference). For a maximum estimate, the 355 difference between the lowest point in a 100 km swath zone along the ice stream is chosen as reference. Compared to our modelling results, this minimum thickening estimate accounts for the effects observed downstream should the sample also be taken downstream. Conversely, the maximum estimate could yield overestimated thickness changes of several hundreds of metres, as it is comparable to the case where a sample would have been collected at the upstream side of the nunatak in our experiments.

360
Our modelling results can also be used to guide the collection of samples for cosmogenic dating. The areas of nunataks located perpendicular to ice flow, where the subglacial signature of a nunatak is not as pronounced, are likely to provide more accurate estimates of regional ice sheet thickness change, as these areas are less impacted by the differences in ice surface steepening. Also, sampling at the nunatak flanks should diminish the ice surface gradient effect on exposure ages, while deviations toward the up or downstream faces would increase age differences (see Sect. 4.1). While nunatak flanks 365 are commonly sampled, there has been no strong preference for such locations (Fig. 2). Taking as an example cosmogenic ages younger than the Last Glacial Maximum (21 kyr before present) to minimise the effect of inherited concentrations from prior exposure, and considering exposure ages computed from 10 Be (due to its wider availability) and 14 C (due to its much lower level of inheritance), samples taken up and downstream at the same elevation interval show significant age differences regardless of the cosmogenic isotope chosen (Fig. S9;Heyman, 2021;Balco, 2021). However, complex exposure histories and 370 different sampling strategies among studies, which target specific thinning histories, introduce a spatial bias to the data, and to our knowledge, no sampling strategy has been designed to test for lags in exposure up and downstream of nunataks. Such experiment would aid the validation and interpretation of our findings. Nevertheless, reporting sample placement relative to ice flow near the nunatak should help in the interpretation of its exposure age, and when comparing with other sites and modelling results.

375
The recommendations for sample collection based on our results also apply to subglacial bedrock locations that are targeted to test for past ice sheet collapse (e.g. Spector et al., 2019). In particular, sampled subglacial bedrock ridges on the downstream side of a nunatak may record past exposure indicative of a thinner-than-present ice sheet, but an equivalent subglacial ridge on the upstream side may record no such exposure. High-resolution ice flow modelling around a sampled nunatak is therefore necessary to understand how representative the sample location is of regional-scale ice loss. Irrespective of sampling strategy 380 and application, it is important to keep in mind that not only modern ice flow direction should be considered, but also past variations in flow patterns. Local signs of palaeo ice flow (e.g. glacial striations) can help with such interpretations.

Implications for modelling ice flow in areas of large topographic relief
Model resolution plays a substantial role in how modelled ice flow interacts with nunataks. Regular-mesh experiments with resolutions typical for ice sheet models that simulate multi-millennial changes (5-20 km), show a stronger elevation gradient 385 between the ice surface up and downstream of the nunataks. These resolutions do not properly capture the prescribed nunatak shape, patterns of ice flow or thinning rates, which are different in each experiment despite the same applied forcing. Palaeo ice sheet models are run at relatively coarse horizontal resolutions (5-40 km; e.g. Golledge et al., 2012;Whitehouse et al., 2012;De Boer et al., 2014;Kingslake et al., 2018;Tigchelaar et al., 2018;Gomez et al., 2020) to keep computational times reasonable, and often use cosmogenic exposure dates as constraints (i.e. to rule out choices of uncertain parameters that yield 390 unrealistic results) or benchmarks (i.e. to assess the model's ability to reproduce the geological record). In light of our results, it is understandable that experiments using a regular grid at these resolutions struggle to closely match ice surface elevation over mountainous regions reconstructed from cosmogenic exposure dating, and do not fully capture their recorded timing and magnitude of ice thinning (Spector et al., 2019;Stutz et al., 2020). In our experiments, a grid-cell size smaller than the glacier width manages to capture the drainage effect to some degree. Still, a higher number of grid cells is needed to properly resolve ice flow through the narrow glaciers between nunataks and the large topographic relief, thus resolving the observed variations in deglaciation age up and downstream. A better representation of ice flow diminishes the overestimation of the ice-surface elevation gradient. To overcome this limitation in palaeo ice sheet models, the model grid could be refined around nunataks so that it properly resolves this pattern of ice flow. The use of adaptive meshes in ice sheet models (e.g. Berends et al., 2021), nested regional models, or downscaled setups, which take a lower-resolution ice sheet model state as 400 boundary/initial conditions, could be potential solutions. Furthermore, due to their appropriate representation of regional ice flow patterns, higher resolution simulations could also help when designing sampling strategies and reconstructing regional thickness changes, thus diminishing the mismatch between modelled and reconstructed ice surface elevation.

Summary and conclusions
Ice flow in regions of complex terrain, where mountain ranges create steep ice-surface profiles, provide challenges for recon-405 structing and modelling past ice sheet changes. The choice of a reference present-day ice surface elevation to be used when determining regional estimates of thickness change from surface exposure dating is not straightforward, and models struggle to match the available elevation reconstructions and their timing of surface exposure. In order to improve our understanding of ice flow over these regions of large topographic relief, we used an ice flow model that represented an idealised portion of an ice sheet. Five experiment ensembles were carried out in order to better understand how the ice surface responds to the presence of 410 nunataks under thinning scenarios, and how the local response compares to the regional response. The first ensemble comprised simulations where different degrees of ice thinning and different shapes (elongated transverse and parallel to ice flow) were tested for a single nunatak. The other four sets were performed to assess the interaction of three obstacles aligned transverse to flow with the ice surface, and to which extent grid resolutions commonly employed by ice sheet models capture smaller-scale ice flow patterns between these obstacles.

415
Overall, we found that the interaction of ice flow with a nunatak results in a steepening of the ice surface, caused by increasing elevation upstream and lowering downstream. Locally, this ice surface mismatch results in an earlier bedrock exposure downstream, and a delayed exposure upstream during ice sheet thinning. At a regional scale, a single nunatak is able to impact the ice surface up to 30 km in both directions along the centre line, while a surface steepening is present transverse to ice flow over the entire extent to which the subglacial continuation of a nunatak stands out from the otherwise linearly sloping bedrock 420 elevation. As a result of this surface mismatch, a difference in the time of exposure (or minimum elevation) of equidistant points from the nunatak summit exists, which was found to be up to 14 kyr. Although a compilation of cosmogenic nuclide ages indicates that samples taken upstream and downstream yield different ages, spatial biases in the sample distribution could have introduced biases in the resulting age distribution. A positive interference between closely spaced nunataks, or a more extensive nunatak perpendicular to flow can further increase the ice surface elevation gradient, while efficient drainage through 425 glaciers formed between the nunataks is able to alleviate it. This mismatch and its consequences should be taken into account when sampling for surface exposure dating and when inferring past ice sheet thickness change, since they directly influence the interpretation of a sample's elevation relative to the regional ice surface.
Finally, we found that the current grid resolutions employed by ice sheet models cannot adequately resolve the flow through glaciers formed in nunatak ranges. The inability to capture smaller scale interactions between ice flow and bed topography 430 results in an overestimated ice surface gradient across these obstacles, and the models miss important variations in the time of bed exposure up and downstream. A resolution of 5 km can to some degree resolve a glacier width of 10 km or more, reducing the overestimation but not entirely solving the problem. Although such high-resolution simulations are currently too computationally expensive to be carried out at full ice-sheet scale, especially over millennia, accurate data-model comparisons with cosmogenic-nuclide dated surfaces require a proper regional refinement around the nunataks. This could be achieved with 435 the aid of adaptive meshes, nested grids, or higher resolution regional models, which should perform better at reproducing the timing and magnitude of ice loss in regions of complex and large topographic relief. Better understanding the relationship between sample location and regional patterns in ice flow and ice surface elevation, combined with improved model simulations over sampled sites, will ultimately allow to account for potential biases in exposure ages and improve the comparison between in-situ data and ice sheet models.