Accumulation patterns around Dome C, East Antarctica, in the last 73 kyr

. We reconstruct the pattern of surface accumulation in the region around Dome C, East Antarctica, since the last glacial. We use a set of 18 isochrones spanning all observable depths of the ice column, interpreted from various ice-penetrating radar surveys and a 1-D ice ﬂow model to invert for accumulation rates in the region. The shallowest four isochrones are then used to calculate paleoaccumulation rates between isochrone pairs using a 1-D assumption where horizontal advection is negligible in the time interval of each layer. We observe that the large-scale (100s km) surface accumulation gradient is spatially stable through the last 73 kyr, which reﬂects current modeled and observed precipitation gradients in the region. We also observe small-scale (10 s km) accumulation variations linked to snow redistribution at the surface, due to changes in its slope and curvature in the prevailing wind direction that remain spatially stationary since the last glacial.

Abstract. We reconstruct the pattern of surface accumulation in the region around Dome C, East Antarctica, since the last glacial. We use a set of 18 isochrones spanning all observable depths of the ice column, interpreted from various ice-penetrating radar surveys and a 1-D ice flow model to invert for accumulation rates in the region. The shallowest four isochrones are then used to calculate paleoaccumulation rates between isochrone pairs using a 1-D assumption where horizontal advection is negligible in the time interval of each layer. We observe that the large-scale (100s km) surface accumulation gradient is spatially stable through the last 73 kyr, which reflects current modeled and observed precipitation gradients in the region. We also observe small-scale (10 s km) accumulation variations linked to snow redistribution at the surface, due to changes in its slope and curvature in the prevailing wind direction that remain spatially stationary since the last glacial.

Introduction
The Dome C region, located on the East Antarctic interior plateau, has long been the focus of extensive research: it is the site of the oldest continuous ice core as yet retrieved , the EPICA Dome C ice core, dating to ∼ 800 ka . Modern surface precipitation on the Dome C plateau is extremely low (∼ 25 mm yr −1 , Stenni et al., 2016), with infrequent storm events representing more than 50 % of the total annual precipitation (Frezzotti et al., 2005). Coastal air masses lose moisture as they are driven inland to higher elevation, resulting in a characteristic precipitation gradient with higher measured and modeled precipitation on the north side of Dome C (Arthern et al., 2006;Genthon et al., 2016;Kållberg et al., 2004;Gallée et al., 2013;Palerme et al., 2014;Van Wessem et al., 2014). Present-day moisture-bearing air mass trajectories (Scarchilli et al., 2011;Genthon et al., 2016) point to a western Indian Ocean provenance for the snow precipitation at Dome C (85 % of the precipitation), and suggest this could have persisted through glacial-interglacial cycles.
Snow precipitation is homogeneous on a large-scale, whereas local variations in snow accumulation are controlled by local surface topography as a function of wind direction. Black and Budd (1964) and Budd (1971) were the first to observe the close relationship between bedrock relief, surface slope and accumulation rates in Wilkes Land. Whillans (1975) details how wind speed and direction can affect total mass balance in Marie Byrd Land. Frezzotti et al. Published by Copernicus Publications on behalf of the European Geosciences Union. (2007) show that surface slope in the prevailing wind direction (SPWD) is a key constraint in determining spatial and temporal variability of precipitation; a higher SPWD can lead to significant ablation and redeposition of snow (Frezzotti et al., 2002a(Frezzotti et al., , b, 2005(Frezzotti et al., , 2007. Das et al. (2013) show that SPWD is a strong threshold for the formation of wind scour or megadune fields. Evidence for a persistent westerly wind circulation pattern comes from mineral dust measured at EPICA Dome C which shows a uniform geographic provenance from South America and Australia to the East Antarctic plateau during glacial-interglacial cycles (Delmonte et al., 2010;Albani et al., 2012).
Airborne and ground-based ice-penetrating radar data have long been used to constrain the surface and bedrock topography over large parts of the Antarctic Ice Sheet (Gudmandsen, 1971;Drewry et al., 1980;Millar, 1981;Siegert, 2003;Bingham and Siegert, 2009, and many others), as well its internal stratigraphy (e.g., Siegert, 1999;MacGregor et al., 2012;Cavitte et al., 2016). As the internal stratigraphy represents isochronal surfaces throughout much of the ice sheet, dated internal radar reflectors can be used to directly constrain the surface mass balance of the ice sheet in the highest part of the ice column (Medley et al., 2013;Verfaillie et al., 2012). Reconstructing accumulation history from deeper isochrones is an ill-posed inverse problem as both accumulation variations and changes in ice flow can affect isochrone geometries (e.g., Koutnik et al., 2016;Neumann et al., 2008;Parrenin and Hindmarsh, 2007;Parrenin et al., 2006;Waddington et al., 2007;Nereson and Waddington, 2002), and assumptions have to be made about one or the other (Martin et al., 2009;Leysinger Vieli et al., 2011;Siegert, 2003;Morse et al., 1998, see companion paper for more discussion). Assumptions on the vertical strain rate will also affect reconstructed paleoaccumulation rates (e.g., MacGregor et al., 2015;Waddington et al., 2007). Several radar isochrone studies have shown the existence of a coastto-dome precipitation gradient: Verfaillie et al. (2012) show a continuous existence of a precipitation gradient through the last 300 years, while Siegert (2003) shows the persistence of a strong accumulation gradient between Dome C and Ridge B (a topographic high upstream of Lake Vostok) over glacial-interglacial timescales.
Better constraining accumulation rates through time is important for several reasons: 1. The spatial distribution of snow accumulation affects the position of topographic domes, which ultimately affects the geometry of the ice sheet (with its resulting sea-level implications) through time (Scarchilli et al., 2011;Fujita et al., 2011;Morse et al., 1998).
2. In addition to other parameters, accumulation rates are required for accurate dating and interpretation of ice cores. Constraints on accumulation and flowline geometries of ice particles through time are necessary to reconstruct ice core chronologies and correct for the effects associated with deposition at a different location and elevation than the ice coring site (Koutnik et al., 2016;. Especially in the context of the search for 1.5 million-year-old ice, such constraints will have a significant influence on the choice of an ice core site. 3. The temporal evolution of accumulation rates provides important constraints on ice sheet mass balance through time for modeling experiments (Koutnik et al., 2016;Fischer et al., 2013).
Here, we reconstruct paleoaccumulation rates for the Dome C region using a 1-D pseudo-steady ice flow model (described in the companion paper: Parrenin et al., 2017) for the last 73 kyr, using the isochronal constraints obtained from radar surveys. We discuss the large-scale accumulation and small-scale variations in accumulation calculated around Dome C.

Dome C region
The Dome C region represents a topographic high in the middle of the EAIS and is at the confluence of several ice divides, the largest of which separates the Byrd Glacier catchment from the Totten Glacier catchment. The topography of the Dome C region is gentle: the change in elevation is ∼ 10 m across 50 km (Genthon et al., 2016), reaching a maximum elevation at Dome C of ∼ 3266 m a.s.l. (geoid height). A saddle connects Dome C to Lake Vostok along the ice divide, with a secondary dome referred to as "Little Dome C" (LDC) just south of the Dome C ice core site. The bedrock is characterized by a large subglacial massif ∼ 40 km to the south of the Dome C ice core site and ∼ 10 km south of the LDC, easily identifiable on Fig. 1, where the radar survey grid is tightest. For ease of description, we refer to it as the "Little Dome C Massif" (LDCM) to differentiate from the surface topographic high. The deep Concordia Subglacial Trench (CST) runs along its eastern edge and is followed by a steep ridge, ∼ 2000 m high , which we will refer to as the Concordia Ridge (CR). Both the LDCM and the CR (see Fig. 1) have been identified as promising targets for retrieving 1.5 million-year old ice (Van Liefferinge and Pattyn, 2013).

Radar data
We use several airborne ice-penetrating radar surveys collected in the Dome C region by the University of Texas at Austin Institute for Geophysics (UTIG) and the Australian Antarctic Division (AAD) as part of the ICECAP project (International Collaborative Exploration of the Cryosphere through Airborne Profiling, Cavitte et al., 2016) and the Oldest Ice candidate A (OIA) survey flown by ICECAP in January 2016 (Fig. 1, . All surveys use the same center frequency of 60 MHz and the same bandwidth of 15 MHz; radar isochrones can therefore be easily matched from one season to the next. A set of 18 isochrones is traced throughout the region, using multiple crossovers, thus ensuring the reliability of the tracing as outlined in Cavitte et al. (2016). All 18 isochrones consistently match at intersections, by virtue of both the spatial density of the radar survey and the relatively smooth geometry of the 18 isochrones in the region (see Supplement 5 for an example radargram). The colocation of the EPICA Dome C ice core in the survey region enables the dating of the isochrones using the AICC2012 chronology (Bazin et al., 2013;Veres et al., 2013). Obtaining ages and associated uncertainties for each isochrone is described in Cavitte et al. (2016). We extend the same isochrones to the newly acquired OIA survey and add a num-  (Veres et al., 2013;Bazin et al., 2013). b Age uncertainties are a combination of the AICC2012 chronology age uncertainties and the radar isochrones' depth uncertainties (Veres et al., 2013;Bazin et al., 2013). ber of isochrones in the OIA region. We use all 18 isochrones for the 1-D model inversion but only use the youngest four isochrones going back into the last glacial (10-73 ka) for paleoaccumulation reconstructions, explained below. All four isochrone depths, ages and uncertainties at the Dome C ice core site are given in Table 1.

Modeling
We use 18 radar isochrones, dated from 10 ka (before 1950) to 366 ka, and the 1-D pseudo-steady ice flow model described in the companion paper (Parrenin et al., 2017). The model inverts for time-averaged geothermal heat flux G 0 , time-averaged accumulation rateā and time-averaged vertical strain rate profile parameter p every kilometer along a radar line. Pseudo-steady-state means that all parameters in the model are considered steady except for R(t), a temporal factor applied to both basal melting and accumulation (see companion paper). R(t) is obtained from accumulation variations inferred from the AICC2012 chronology (Veres et al., 2013;Bazin et al., 2013) and represents the ratio of the accumulation at time t to the average accumulation over the last 800 kyr. When inverting the radar isochrones using the pseudosteady ice flow model, ages and accumulations are all used in steady-state form, with glacial-interglacial accumulation variations normalized. The calculated time-averaged accumulation rateā, p and G 0 result from the best fit of all the radar isochrone depths (dropping x for simpler notation). However, some differences between modeled and observed isochrones remain as all isochrones have to be simultaneously fitted for each point x. The 18 isochrones have to be used in the inversion as the deepest isochrones provide the strongest constraints on p and G 0 .
To reconstruct paleoaccumulation rates through time, we use the G 0 and p values calculated and assume they remain unchanged over each time so that the remaining misfit between modeled and observed isochrones is entirely a result of the uncertainty in a. The calculated a χ represents the paleoaccumulation rate for a layer with an age interval χ ,  bounded above and below by a radar isochrone of AICC2012 age. We refer to these as isochrone-bounded layers. To calculate a χ values for each layer, we adjust the value of a such that modeled and observed isochrone-bounded layer age intervals χ are fitted perfectly for each layer.
In mathematical form, if z is the depth of the isochrone and χ the age of the isochrone, we can write a o, χ , the "observed" paleoaccumulation rate for a certain age interval χ as follows: where τ m is the vertical thinning function obtained from the model, i.e., the ratio of the vertical thickness of a layer to its initial vertical thickness at the surface and χ is the age interval. This is similar to the "shallow-layer approximation" used by Waddington et al. (2007).
Using Eq.
(1), we calculate the best fit paleoaccumulation rates through time in one iteration after the model inversion. This gives the spatial variations of the paleoaccumulation rates through time.
Care must be taken in not over-interpreting the paleoaccumulation maps obtained. We do not argue that we have reconstructed absolute paleoaccumulations for the past 73 kyr. The 1-D pseudo-steady ice flow model used here (see companion paper, Parrenin et al., 2017) does not take horizontal advection into account. Paleoaccumulation rates calculated are valid at the ice divide and the dome where horizontal ice flow speeds are negligible. Farther away, horizontal advection has a larger influence. A full 3-D model is required to reconstruct accumulation rates more extensively in space and further in time.
To respect our assumption that τ is modeled perfectly, we only calculate paleoaccumulation rates a o, χ for the first four isochrone-bounded layers. Our fourth and deepest layer The Cryosphere, 12, 1401-1414, 2018 www.the-cryosphere.net/12/1401/2018/ Figure 3. Time-averaged accumulation ratesā along the radar lines over the Dome C region over the last 73 kyr. Accumulation rates are given in mm of water equivalent per year. There is a clear large-scale north-south accumulation gradient, with accumulation decreasing with distance from the Indian Ocean coast, the main pathway of snow precipitation. Black lines outline areas of small-scale high accumulation: they correlate to areas where surface contours (in gray) become further apart, i.e., where surface slope is reduced. Background is the same as in Fig. 1. used reaches an average depth of 30 % of the ice thickness, with calculated thinning never reaching below 0.6. Furthermore, to avoid ill-posed conditions for our 1-D paleoaccumulation reconstructions, we only retain data points that have experienced a maximum of 5 km of horizontal advection. We do this for each layer, using Van Liefferinge and Pattyn (2013) ice surface balance velocities, corrected for temporal velocity variations using R(t) (Parrenin et al., 2017), and the age interval spanned by the layer considered. Any point that has traveled more than 5 km horizontally is masked.
We show the paleoaccumulation rates calculated for the four youngest age intervals spanning 0-73 ka in Figs. 2 and 4.
The Metropolis-Hastings (MH) algorithm (described in the companion paper, Parrenin et al., 2017) enables the calculation of an accumulation rate uncertainty (see S2) which takes into account the age uncertainty of the radar isochrones. The age uncertainty of the radar isochrones is a combination of the radar depth uncertainties translated to age uncertainties (Cavitte et al., 2016) and the AICC2012 ice core chronology uncertainties (Veres et al., 2013;Bazin et al., 2013). Cavitte et al. (2016) describe the various sources of radar depth un-certainty and how they are calculated. The radar isochrone depth and age uncertainties are given in Table 1. We plot the time-averaged accumulation rate and the paleoaccumulation rates for each isochrone-bounded layer over the survey region (see Figs. 3 and 4). The accumulation rate uncertainties are given in Fig. S2 and further discussed in Supplement 2.

ECMWF ERA40 snow precipitation rate
We compare our calculated large-scale patterns of precipitation to present-day measurements, ECMWF (European Center for Medium-Range Weather Forecasts) ERA40 reanalysis data  is used to obtain a map of present-day estimated precipitation rates over the survey region. The snow accumulation rates in the Dome C region result from precipitation in the form of snow (snowfall and diamond dust) and are then modified by wind-driven processes. The wind erosion, wind redistribution and sublimation, as well as other processes during or after a precipitation event, leads to spatial deposition at the surface that is much less homogeneous than the original precipitation (e.g., Frezzotti et al., 2004). The ECMWF ERA40 model seems to correctly reproduce the observed precipitation's spatial and temporal variability at Dome C, but systematically underestimates the precipitation magnitudes (Genthon et al., 2016;Stenni et al., 2016), probably because clear-sky precipitation is not adequately parameterized (Bromwich et al., 2004;Van de Berg et al., 2006) and blowing snow transport and sublimation processes are not considered. However, since the Dome C site is not influenced by strong winds, this is expected to have a minor effect within the summit area, but cannot be completely neglected farther than 25 km from the dome and ice divide. ECMWF ERA40 data have been normalized using the surface accumulation average of the last centuries from existing ground-penetrating radar (GPR) within 25 km from Dome C summit (Urbini et al., 2008).
A number of steps went into adjusting the ECMWF ERA40 modeled precipitation rates to field measurements, to calculate the "ECMWF ERA40 estimated present-day surface accumulation rates", shown in Fig. 5. These steps are as follows: 1. ECMWF ERA40 monthly average precipitation rates were used to calculate a long term precipitation average over the 1989-2011 period; 2. precipitations were then interpolated over the region of interest as a 1 km grid; 3. precipitation values were increased by 12.9 mm yr −1 to match GPR measurements in the area (Urbini et al., 2008) as ECMWF ERA40 precipitation values are systematically too low compared to ground-based measurements.

Spatial trends of paleoaccumulation
To look at small-scale paleoaccumulation variations more closely, we remove large-scale precipitation gradients (see Sect. 5). For this, we calculate a quadratic fit of the ECMWF ERA40 present-day surface accumulation values (calculated as described above) with each isochrone-bounded layer's paleoaccumulation and subtract the calculated fit from the layer's paleoaccumulation values. The result is a map of The Cryosphere, 12, 1401-1414, 2018 www.the-cryosphere.net/12/1401/2018/ Figure 5. Historical average accumulation rates a 100 yrs along the radar lines superimposed on ECMWF ERA40 estimated presentday surface accumulation rates (see Sect. 3.3). There is a very good agreement in the magnitude of accumulation values between the two datasets and in their north-south accumulation gradient on large-scales (100s km), with accumulation decreasing with distance from the coast. White lines outline the same areas of small-scale high accumulation as in Fig. 3. Background is the same as in Fig. 1. detrended paleoaccumulations for each isochrone-bounded layer (Fig. 6).

Slope and curvature in the prevailing wind direction (SPWD and CPWD)
In Sect. 5, we discuss the importance of surface slope in the prevailing wind direction (SPWD) and curvature in the prevailing wind direction (CPWD). We use ECMWF 5 year average wind directions  and Bamber et al. (2009) surface elevations to calculate SPWD and CPWD values over a 3 km radius in the survey region (Fig. 6). A positive value of surface curvature indicates a surface trough, while a negative value of surface curvature indicates a surface bump.

Results
The reconstructed paleoaccumulations a o, χ are shown in the top panel of Fig. 2 along the A-A' radar transect (VCD/JKB2g/DVD01a) marked on Fig. 1. The A-A' radar line runs along the ice divide, and a marked decreasing accu- mulation rate can be seen going from the north-east side towards the south-west consistently over all age intervals. The lower panel of Fig. 2 displays a o, χ along the B-B' radar transect (OIA/JKB2n/Y77a) marked on Fig. 1. This transect runs across the divide and again there is a visible spatial accumulation gradient, clearly observable for the youngest layer. We also show reconstructed accumulation rates in map view in Figs. 3 and 4. Figure 3 displays the time-averaged accumulation rateā and Fig. 4 displays the paleoaccumulation rate per isochrone-bounded layer a o, χ . We show all four age intervals calculated. We observe that the time-averaged accumulation (Fig. 3) has a clear north-south gradient, decreasing from > 21 mm w.e. yr −1 (water equivalent per year) in the north to 15 mm w.e. yr −1 in the south. Superimposed, we observe a number of regions ∼ 20 km wide that show a ∼ 25 % accumulation increase over the LDCM, to ∼ 50 km wide or more east of the CR with a ∼ 75 % increase. These are outlined by black lines on Fig. 3. Around the CR, we also note that the extended area of high accumulation is adjacent to an area of very low accumulation, parallel to it and just east of www.the-cryosphere.net/12/1401/2018/ The Cryosphere, 12, 1401-1414, 2018 the CST. This corresponds to an area of drastic surface slope and curvature change (see also Figs. 6, S3 and S4). The spatial pattern of paleoaccumulation rates per isochrone-bounded layer (Fig. 4) is similar to that of the time-averaged accumulation: a large-scale gradient northsouth with superimposed areas of higher accumulation in the same locations as for the time-averaged accumulation reconstruction. We note a striking similarity between the timeaveraged accumulation rate (Fig. 3) and the paleoaccumulation rates for the ages 0-10 ka (Fig. 4). We also note that accumulation rates are higher for the interglacial age interval (0-10 ka) than for the glacial age intervals (see Figs. 4 and  S4). The small-scale accumulation patterns are visible in the 0-10 ka age interval, we see the same three areas of high accumulation as outlined in Fig. 3. For older layers, the smaller spatial extent of the paleoaccumulation data makes it difficult to conclude on the persistence of these small-scale high accumulation areas.
Bedrock elevations from Bedmap2 (Fretwell et al., 2013) augmented with new OIA survey data outlined with a dashed rectangle , as well as Bamber et al. (2009) surface elevation contours, are plotted in the background of Figs. 3 and 4. The areas of higher accumulation are co-located with areas of low surface slopes, visible from the surface contours. The accumulation variations we observe are also co-located with significant bedrock relief changes, which reach e.g., ∼ 2000 m for the CR escarpment, and ∼ 500 m for the south side of the LDCM (see Fig. S1).
We use the time-averaged accumulation,ā, obtained from the model and R(t) to plot historical average accumulation rates so as to compare with the ECMWF ERA 40 model results directly. For this, we take the ratio of the average accumulation rate of the last 100 years to that of the last 800 kyr using the AICC2012 chronology, which has a value of 0.65. The time-averagedā is scaled by this factor of 0.65 to obtain historical average accumulation rates, which we call a 100 yrs here. We plot a 100 yrs together with ECMWF ERA40-derived surface accumulation data in Fig. 5 (see Sect. 3.3 for details). We observe that the large-scale north-south accumulation gradient in a 100 yrs closely resembles that of the ECMWF ERA40-derived surface accumulation rate in Fig. 5: high accumulation in the north nearer to the coast, and lower accumulation in the south as you move towards the interior. The magnitude of the accumulation rates also matches surprisingly well.
The calculated accumulation rate uncertainties from the model, with an average value of 0.16 mm w.e. yr −1 (see Fig. S2), are an order of magnitude (or more) smaller than the values of reconstructed time-averaged accumulation rate, providing confidence in the time-averaged accumulation rates calculated. However, errors have been treated as uncorrelated so we cannot apply these uncertainties to the paleoaccumulations. We hope to improve this in the future.
To focus on the small-scale variations in paleoaccumulations, we plot detrended paleoaccumulations (see Sect. 3.4) for the region on top of SPWD and CPWD values (see Sect. 3.5), as shown on Fig. 6. We only show this relationship for the first layer, spanning the past 10 kyr, as older layers are not as extensive. Looking at the spatial distribution of these detrended paleoaccumulations in relation to SPWD, we observe that areas with high accumulation are co-located with areas of markedly reduced SPWD values with respect to the surrounding values (∼ 0.5-1.2×10 −3 of absolute SPWD decrease, see Fig. S3). But more striking is the relationship between the magnitude of the curvature (and polarity) and the magnitude of the residual paleoaccumulation (Fig. 6). The areas of high accumulation in Fig. 3 are outlined in black. They correspond to areas of high positive detrended paleoaccumulation, > 1.2 mm w.e. yr −1 , and are well correlated with areas of strongly positive curvature values (> 2×10 −7 m −1 ). This is evident in the LDCM area. Areas of high negative detrended accumulation, < −1.6 mm w.e. yr −1 , are also well correlated with areas of strongly negative curvature. This is best seen east of the CR. The correlation holds particularly well for the youngest layer (0-10 ka) over the entire region. We plot detrended paleoaccumulation for layers older than 10 ka and observe that this relationship holds over the LDCM with a slightly increasingly offset with increased ages (see S4).

Reconstruction uncertainties
The 1-D assumption to calculate paleoaccumulation rates is clearly the largest source of uncertainty in our reconstructions. In the 1-D pseudo-steady ice flow model described in the companion paper (Parrenin et al., 2017), the goal is to constrain the age of the deep ice. For that work, trade-offs in the strain thinning (i.e., p and G 0 ) and accumulation rates do not matter, as their combined effects dictate the age of the ice. However, to calculate the layer-by-layer paleoaccumulation rates, we have to assume that τ (Eq. 1) is modeled perfectly, which breaks down as horizontal advection increases. We reckon that for the first layer whose average depth is ∼ 150 m, that is ∼ 5 % of the total depth, the error in the thinning is small enough to not affect significantly our accumulation results (after strain thinning, the layer is always at least at 90 % of its original thickness). For the other layers, it is difficult to imagine an error in the thinning function that would produce, by chance, a similar accumulation pattern to that of the first layer. In addition, by setting a limit on the maximum horizontal advection allowed for each age interval, the described accumulation patterns and variations are reasonably unaffected by the 1-D assumption. The threshold of 5 km is chosen such that horizontal advection is negligible compared to the scale of the observed accumulation rate variability. The small-scale areas of high accumulation are at least 20 km wide in the region, therefore the 5 km thresh-The Cryosphere, 12, 1401-1414, 2018 www.the-cryosphere.net/12/1401/2018/ old on horizontal movement does not affect our conclusions. We are only able to reconstruct paleoaccumulation rates back through 73 ka, therefore a 3-D model is required to look at paleoaccumulation rates further back in time. Furthermore, the model assumes a constant ice thickness through time. Even though small variations in the ice thickness through time will affect the absolute value of the reconstructed accumulation rates, the assumption of constant ice thickness is fair for the center of the EAIS where modeled ice thickness variations have been reported below 200 m (Bentley, 1999;Ritz et al., 2001; (representing a 5 % error on the ice thickness) and little is known of the spatial distribution of these ice thickness variations in the center of the ice sheet. A 5 % error on the ice thickness will produce a 5 % error on the thinning function τ  and therefore a 5 % error on the accumulation rates calculated. This error can be ignored for two reasons. First, it is small compared to the accumulation variations that we observe (larger than 10 %), and second, it only affects the absolute value of the accumulation rates reconstructed but not the relative differences in accumulation rates from one location to the next in the Dome C region. Since we focus exclusively on changes in gradients and patterns in accumulation rates, this additional source of error does not affect our conclusions. Despite this error, we observe a clear reduction in the magnitude of the accumulation rates as we go back in time and enter the last glacial maximum, as expected and measured in ice cores (Bazin et al., 2013;Veres et al., 2013;Parrenin et al., 2017).

Large-scale spatial stability of paleoaccumulation rates
The observed patterns of paleoaccumulation agree well with previous studies of surface snow accumulation variability in the Dome C region. Considering first the large-scale patterns in the accumulation reconstructions, we observe a consistent large-scale gradient (large-scale here refers to 100s of km) for each age interval, with accumulation values decreasing from the north side of Dome C to the south side. Scarchilli et al. (2011) suggest moisture provenance from the Indian Ocean sector is the most consistent with the clear north-south gradient in precipitation observed as we near Dome C. The fact that our paleoaccumulation reconstructions reproduce the present-day large-scale surface accumulation gradient and that this remains true back to 73 ka suggests persistence of the source of moisture for this part of the East Antarctic plateau through the last glacial and deglaciation. Transects A-A' and B-B' in Fig. 2 clearly show the north-south orientation of the accumulation gradient. This large-scale accumulation gradient is also clearly seen in the ECMWF ERA40 data for the region (Fig. 5), as well as in other large-scale accumulation models of the region (e.g., Genthon et al., 2016) or Regional Climate Model (MAR) (Gallée et al., 2013(Gallée et al., , 2015. GPR data collected dur-ing traverses across Dome C and along the divide also show a strong north-south gradient in accumulation (Urbini et al., 2008;Verfaillie et al., 2012;E. Le Meur, personal communication, 2016). We note a good agreement between our accumulation values and trends along a-a' going from Dome C along the ice divide towards Vostok (top panel of Fig. 2) and the GPR transect measured by Verfaillie et al. (2012) on the other side of the Dome C divide. A SPRI airborne transect collected over Dome C shows a strong accumulation gradient of 10s of mm yr −1 over a spatial scale of 100s of km (Siegert, 2003). Urbini et al. (2008) show a small component of counter-clockwise rotation of the accumulation pattern in historical times centered on Dome C, but the general north-south gradient difference in accumulation across the dome remains. Measurements made in other areas of the ice sheet, e.g., across Talos Dome (Frezzotti et al., 2007), point to similar patterns: accumulation is highest near the moisture source and decreases with distance from the coast. Fujita et al. (2011) point to the same patterns of reduced accumulation inland across Dronning Maud Land.

Small-scale spatial stability of paleoaccumulation rates
Considering the small-scale (10s of km or a few ice thicknesses) patterns of accumulation shown earlier, we described several regions of locally increased accumulation. The colocation of the areas of higher accumulation with areas where surface slope is reduced, as seen from the surface contours or the markedly reduced SPWD values with respect to the surroundings (Fig. S3), fits well with the model put forward by Frezzotti et al. (2007). Frezzotti et al. (2007) show that accumulation increases when SPWD decreases over Talos Dome and attribute the correlation between the absolute magnitude of SPWD and accumulation rates to katabatic wind-driven ablation. Note that the prevailing wind direction over the area is more or less along the long axis of Dome C flowing from higher up the ice divide towards Dome C (see Fig. 6, Frezzotti et al., 2005;Urbini et al., 2008). The spatial correlation we obtain between the detrended paleoaccumulations and the CPWD can be explained by the same mechanisms as for SPWD, since SPWD and CPWD are directly related. Layer 0-10 ka shows high detrended paleoaccumulation values where the surface curvature is strongly positive (i.e., surface trough), and low values where the surface curvature is strongly negative (i.e., surface bump). This is true for both the LDCM and CR regions. The proximity of the isochrone-bounded layer to the surface influences how well the correlation holds, particularly visible in the CR region which is furthest from the divide. For any deeper layer (Fig. S4), this relationship is slightly offset in space; a likely cause is the increased amount of horizontal advection with depth, up to the set maximum of 5 km.
Even though the absolute magnitudes of slope and curvature changes we observe are relatively small (on the orwww.the-cryosphere.net/12/1401/2018/ The Cryosphere, 12, 1401-1414, 2018 der of 10 −3 and 10 −7 m −1 , respectively), other studies have shown that even very small slope changes can have a strong influence on wind-borne redistribution of snow (Grima et al., 2014;King et al., 2004;Whillans, 1975). However, a single mechanism has yet to be described that would explain the relationship between CPWD (and therefore SPWD) and small-scale accumulation variations. Grima et al. (2014) observe strong surface density variations linked to surface slope breaks; however, some increases in accumulation occur over steeper surface slopes, which is surprising when steep slopes are usually associated with reduced accumulation (Hamilton, 2004;Frezzotti et al., 2004). King et al. (2004) show that local slope changes of 0.01 can create up to 30 % variations in accumulation, and invoke a highly non-linear relationship between wind speed and snow transport to explain the type of accumulation variability they observe. Whillans (1975) also shows that slope changes as small as 0.001 over a distance of 3 km can affect snow deposition, and argues for a relationship between slope, wind strength and mass drift transport. The extreme pattern of high and low accumulation parallel to the CST and east of the CR seems to be the ideal example of how surface topography variations affect accumulation rates. The ice flowing radially away from Dome C has to flow over the CST and the prominent bedrock CR. CPWD shows strongly negative values over the subglacial CR; it creates a surface which is concave down perpendicular to the wind direction. We can imagine a scenario in which snow is strongly plucked away on this steepest surface slope, but further down-wind, as slope reduces and reaches contrastingly strongly positive CPWD, the snow can then be redeposited directly down-wind as suggested in Frezzotti et al. (2004).
We noted in the results that the small-scale accumulation variations were co-located with bedrock relief variations (see Fig. S1). Frezzotti et al. (2007) explain that bedrock topography can be the underlying influence on the variability of snow accumulation at scales of 1-20 km, corresponding to the length-scales of the accumulation variations we calculate here. Bedrock topography will have a stronger influence on the overlying ice in the presence of subglacial lubrication (Rémy et al., 2003). Rémy et al. (2003) show that for the Dome C region, the most positive surface curvatures are directly linked to the largest ice thicknesses and the presence of subglacial lakes. It is interesting to note that areas of higher detrended paleoaccumulation correlated to high positive CPWD outlined in Fig. 3 are above deep bedrock valleys dotted with many observed subglacial lakes .
We attempted a series of low order (linear and quadratic) fits between CPWD and our detrended paleoaccumulations but none explain all the variability. The data is suggestive of threshold behaviors between low and high CPWD magnitudes.. For the youngest layer, we calculated a correlation coefficient of 0.255 for high CPWD magnitudes and 0.103 for low CPWD magnitudes, using ±2 × 10 −7 as cutoffs (Fig. 7). This is perhaps not surprising as we correlate  Figure 7. Relationship between CPWD and detrended accumulation rates for the youngest layer (0-10 ka). Magenta and red lines are linear best fits for regions of high (>2 × 10 −7 or < −2 × 10 −7 ) and low (between −2 × 10 −7 and 2 × 10 −7 ) CPWD magnitudes, respectively. The correlation coefficient is higher for regions of high CPWD magnitudes. present-day surface topographic features with 10 kyr of surface accumulation history. We therefore expect small bumps and troughs in the surface topography to be prone to spatial migration, while more substantial troughs co-located with subglacial topography (e.g., Fig. 6) might remain relatively fixed. ECMWF wind speed magnitudes over the LDCM and CR areas  are below the 5 m s −1 threshold for dune processes to be active in the region, and the radar data used do not show any buried dune structures. The accumulation patterns observed are more suggestive of the preferential infill of surface troughs by winds. These troughs might not fill-up easily because of the very low surface precipitation rates in the region (Genthon et al., 2016;Urbini et al., 2008) combined with the presence of areas of subglacial melting in the region , creating additional draw-down of the surface.
Although we cannot yet explain the mechanisms causing the small-scale paleoaccumulation variability we observe in the Dome C region, which is beyond the scope of this manuscript, our observations have important ramifications for better understanding the region's stability through time.
In the future, we hope to improve our paleoaccumulation rate reconstructions, and in particular go back further into the last glacial cycle with a full 3-D model. Further GPR data was recently collected over the LDCM, and strain nets and various other instruments have been deployed. These new measurements will add to the existing data set and provide important constraints if we hope to develop 3-D inversions.

Conclusions
We reconstructed accumulation rates for the last 73 kyr. Looking at both large-and small-scale accumulation gradients, we show that these have not changed significantly since the last glacial. Large-scale accumulation gradients will remain constant if moisture-bearing air mass trajectory interactions with surface topography do not vary. Small-scale accumulation variations are strongly controlled by SPWD and CPWD and therefore, if the pattern of high and low accumulations remains fixed over a long period of time, this requires consistent interactions between local surface slopes and prevailing winds over the last 73 kyr, independent of whether the control comes from the bedrock topography and/or potential basal melting. This points to a spatially stationary and persistent accumulation pattern in the Dome C region over the last glacial, an important constraint for modeling efforts in the area, both for dating existing ice cores as well as for the prospecting of a 1.5 million-year-old ice core site.
Data availability. The radar isochrones used in this manuscript will be made publicly available in 2018 as a separate publication. Code for the model is available publicly under https://github.com/ parrenin/IsoInv.
Author contributions. MGPC interpreted and analysed the radar isochrones, FP developed the model and ran experiments with MGPC with CR input, DAY, JLR and DDB were involved in survey design and data acquisition, BVL and MF provided data and discussion material. MGPC prepared the manuscript with contributions from all co-authors. The authors declare that they have no conflicts of interest.