Observation of the process of snow accumulation on the Antarctic Plateau by time lapse laserscanning

Snow accumulation is the main positive component of the mass balance in Antarctica. In contrast to the major efforts deployed to estimate its overall value on a continental scale – to assess the contribution of the ice-sheet to sea-level rise – knowledge about the accumulation process itself is relatively poor, although many complex phenomena occur between snowfall and the definitive settling of the snow particles on the snowpack. Here we exploit a dataset of near-daily surface elevation maps recorded over three years at Dome C using an automatic laserscanner sampling 40–100 m in area. We find that 5 the averaged accumulation is relatively regular over the three years at a rate of +8.7 cmyr−1. Despite this overall regularity, the surface changes very frequently (every 3 days on average) due to snow erosion and heterogeneous snow deposition that we call accumulation by "patch". Most of these patches (60–85%) are ephemeral but can survive a few weeks before being eroded. As a result, the surface is continuously rough (6–8 cm root mean square height) featuring meter-scale dunes aligned along the wind and larger, decameter-scale undulations. Additionally, we deduce the age of the snow present at a given time on 10 the surface from elevation timeseries and find that snow age spans over more than a year. Some of the patches ultimately settle, leading to an heterogeneous internal structure which reflects the surface heterogeneity, with many snowfall events missing at a given point, whilst many others are over represented. These findings have important consequences for several research topics including surface mass balance, surface energy budget, photochemistry, snowpack evolution and the interpretation of the signals archived in ice cores. 15

the 1D Crocus model in parallel. Each simulation represented a decametre-scale cell but there was neither spatial organisation nor notion of neighbourhood between the cells. The simulations were mostly independent except that they received different amount of snow during each snowfall and they exchanged snow during strong wind events (local erosion and redeposition).
These two processes were implemented using stochastic rules, adjusted to mimic some in-situ observations collected at Dome C. The approach was quite successful to reproduce the variability observed in one hundred profiles of snow density and 5 specific surface area collected during a summer campaign. However, their rules lack a physical basis and were empirically parameterised with a limited set of observations.
The aforementioned pioneering modelling studies are limited by the scarcity of observations providing detailed information on the accumulation process. Here, we exploit a new dataset of daily surface elevation maps obtained with an automatic ground-based laserscanner at Dome C (Picard et al., 2016) scanning a surface area of 40-100 m 2 and which operated during 10 three years. The aim is to answer three main questions: 1) what are the spatial and temporal characteristics of the accumulation patterns, 2) how long does snow remain on the surface before being eroded or definitely incorporated in the snowpack, and 3) what is the impact of the accumulation process on the snowpack's upper internal structure ? The study focuses on the meter-scale and daily to annual time scale and follows a statistical approach to describe the accumulation and erosion patterns and their dynamics. The paper is organised as follows: Section 2 presents the data and the algorithms developed to extract 15 information (accumulation, age of the snow on the surface, ...) from the series of elevation maps, Section 3 presents the results and Section 4 provides a discussion.

Rugged LaserScan (RLS)
The rugged laserscan (RLS) is composed of a lasermeter mounted on a two-axis rotation stage to perform the elevation and 20 azimuthal rotations, enabling 2-D scanning of the surface. It is described in detail in Picard et al. (2016) and only limited information is recalled here. The lasermeter (DIMETIX FLS-CH 10) measures the radial distance to the snow surface with an intrinsic accuracy of ±2 mm (statistical confidence level of 95.4%) which proved to be effective for a wide range of illumination and temperature conditions (Picard et al., 2016). To achieve this constant accuracy however, the rate of measurements, of 20 Hz in optimal conditions, is automatically reduced when the conditions are unfavourable. The consequence for our particular setup 25 where the lasermeter is rotated at a constant speed, is a reduced spatial resolution. We also found that despite a design for outdoor operations, the lasermeter performance was greatly improved in high light illumination by adding a band-pass optical filter at the laser operating wavelength (650 nm) on the optical window. The lasermeter is heated with internal components and regulated at 0°C for maximal stability of the internal time reference. We fitted an additional heating patch (20 W) on the external box, enabling operations up to -80 • C, which also contributes to the removal of the frost and snow (by sublimation) 30 that occasionally builds up on the device.
The two-axis stage is composed of two identical reduced motors controlled by a feedback loop on the position (servomotor).
The precision and accuracy are of the order of 0.03 • and 0.1 • respectively, which is small but nonetheless is the main limiting 3 The Cryosphere Discuss., https://doi.org /10.5194/tc-2019-4 Manuscript under review for journal The Cryosphere Discussion started: 14 January 2019 c Author(s) 2019. CC BY 4.0 License. factor in terms of accuracy. The scan is performed by moving the azimuth (horizontal) stage at constant speed from nearly -90 • to +90 • , then by increasing the zenith angle (angle from the vertical axis) by a small increment, and finally by moving the stage back from +90 • to -90 • . The process is repeated many times for zenith angles from 17 • up to 62 • . The increment in zenith angle and the speed in azimuth are not constant, they are calculated to obtain a uniform measurement sampling over the whole area. Nevertheless the resolution effectively obtained also depends on the actual lasermeter rate which can vary a lot. A 5 normal scan contains about 200,000 points and takes a total of 4 hours to complete. In Picard et al. (2016) we found that the vertical accuracy is better than 1 cm on the long term (a year). The reproducibility within a scan is better than 5 mm.
Scanning is scheduled at 9 pm local time (GMT+8) every day to avoid high illumination. However, not all the scans are completed with sufficient points to produce a useful map. The main reasons include 1) downtime due to major failures of the RLS or for maintenance, 2) snow or frost deposition on the laser window, and 3) blowing snow crystals intercepting the laser 10 beam, which greatly reduces the acquisition rate. The two latter causes of failure obviously occur during or after snowfalls and blowing snow events. This results in fewer valid scans in the periods of greater surface change. This correlation between the observation quality and the observed phenomenon represents a potential source of statistical bias which must be kept in mind for the analysis. Other, more occasional, reasons of scan failure include power supply shutdown and interruption of the scanning process due to undocumented errors raised by the lasermeter. 15 The RLS was operating at Dome C ( Fig. 1) over 2 periods with different configurations (Table 1). From 1 January 2015 to 17 January 2016, it was set up at a height of 2.8 m (period 1). A major failure occurred during this period between 17 October 2105 and 5 December (49 days). After the maintenance during the summer campaign, it was reinstalled at a height of 4.5 m (period 2) on 1 February 2016. Because of the difference of height during period 1 and 2, the scanned area is different: 40 m 2 and 110 m 2 respectively, and the effective spatial resolution was accordingly adjusted to 2 cm and 3 cm respectively. Although 20 the area scanned during period 1 was located inside the area scanned during period 2, no attempt to co-register the two was done. The two time-series are interpreted as independent datasets. Only the mean surface elevation at the end of period 1 is taken as reference for shifting the starting elevation of the period 2 in order to get a continuous mean elevation time series. Furthermore, after a year of operation in period 2, the outer sheath of the lasermeter cable got damaged by the recurring friction during the azimuthal rotations. The electric contacts became less and less reliable. Nonetheless, the lasermeter continued to 25 operate but returned a decreasing amount of valid data until the dismounting on 28 January 2018. Because of the random nature of the problem, a few scans are useful, at least to estimate the mean elevation of the surface. Hence, we split the period 2 into a first high quality part, named 2a, and a second part named 2b where only the few best scans are selected (Table 1). Despite a much lower temporal resolution, the period 2b time-series still provides a useful extension of nearly a year. In the following, we mainly use the period 2a for its finer temporal resolution and wider scanned area, while the periods 1 and 2b are only exploited 30 to investigate the general trends over the 3 years.

RLS data processing
Processing the raw data to produce elevation maps on a common and regular grid is performed in several steps. For each single acquisition, raw data comprise the radial distance and two angles (azimuth and zenith). After filtering the obviously erroneous 4 The Cryosphere Discuss., https://doi.org /10.5194/tc-2019-4 Manuscript under review for journal The Cryosphere Discussion started: 14 January 2019 c Author(s) 2019. CC BY 4.0 License. distances (less than 3 m or more than 17 m), the data are projected into Cartesian coordinates (x, y) with z the vertical axis.
z(x, y) is the surface elevation. A second filter is applied to remove the points (x, y) for which |z(x, y) −z| > 5 cm, withz the averaged elevation of the points in a 5 cm-radius circle centred at (x, y). This operation removes outliers and small objects like blowing crystals or the RLS mast guy lines. This set of curated but irregularly spaced points is then interpolated onto a regular grid using the bilinear method interpolation provided by the matplotlib.mlab.griddata Python function (version 2.2). The grid 5 spacing is set to 2 and 3 cm respectively for the periods 1 and 2, in relationship with the different setup heights. To avoid filling large gaps with the bilinear interpolation, a grid point was attributed a valid z value only if at least one measurement was taken within 6 cm around it (independently of the period and grid space). If the final map contains less than 100,000 valid grid points, it is completely discarded. All these operations yield a time-series of elevation maps on a common grid. The grid orientation was determined by observing the shadow of the RLS mast in the scanned area. We found that the x-axis lies towards 10 116 • (East-South-East), which allowed us to draw the South direction on the maps.
Various additional datasets are then derived from the generated maps. The accumulation between every successive scan is calculated as the difference between the elevation for each point of the grid. In most cases, both scans are acquired one day apart (90% over the period 2a) so the computed accumulation is representative of the daily variation. However, due to scan failure or rejection during the processing, longer time intervals are present in the time-series (2 days in 7% of the cases, and 15 3-9 days in the remaining 3%). We nevertheless interpret the accumulation time series as if it was daily accumulation only.
The age of the snow on the surface is another dataset derived from the elevation maps. The algorithm is applied independently for each point as follows. For each point, the age is initialised to 0. For each date i the accumulation a i since the last detected deposition or erosion event (or the beginning of the time-series) is calculated. If a i is larger than a given positive threshold δ a , there was deposition of new snow, and the age is reset to 0. If this amount a i is null or positive but lower than δ a , the variation 20 is considered insignificant and the age is incremented by the time elapsed since the previous available date. If this amount is negative a i < 0, snow was removed by erosion. In this latter case, the algorithm then searches back in time the first date j for which the elevation was equal (or closest below) to the present elevation at the point. The snow present at the surface at the date j is now emerging again. The age of this snow is calculated as its age at date j plus the time elapsed since j, that is i − j. We choose a threshold value δ a = 1 cm which tends to avoid sporadic age change for small variations (possibly due to 25 residual noise). This algorithm is relatively robust because even if δ a is chosen too close to the noise level, age errors are not accumulated over time, and only the statistics of the age at a given date may be affected. A low threshold tends to bias the age distribution towards younger snow due to the over-estimation of new snow accumulation. In addition to the age, we derive the time of residence on the surface, which is the same as the age except that when erosion is detected, the time on the surface at time i is taken as equal to the time on the surface at time j, not counting the time spent while buried (i−j). From it, we also 30 derive the time of residence before erosion, which is equal to the time spent on the surface just before erosion is detected.

Meteorological data
Meteorological data are used to relate surface changes and weather. We use precipitation forecast provided by the ERA Interim reanalysis (ERA-I) (Dee et al., 2011) near Dome C. Due to the absence of reliable methods to record in-situ precipitation in Antarctic conditions, ERA-I is one of the most reliable sources of information to date (Bromwich et al., 2011;Wang et al., 2016). For wind speed and direction, we use data from the Weather station at Concordia station ( 75.1 • S, 123.3 • E, 3230 m a.s.l., http://www.climantartide.it/). Over the three years, the mean wind speed at 3 m is 7 ms −1 which is higher than the ERA-I prediction of 5.2 ms −1 at 10 m. Nevertheless, both sources agree on the temporal variations, with a correlation of 0.8, which is the most important point for our comparison. The wind regime at Dome C is characterised by a prevailing direction from the 5 South (74% of the time), the most likely direction being 190 • . The distribution and maximum remain identical when selecting only strong winds (e.g. over the mean, 7 ms −1 ) which are relevant for this study.

Results
The RLS dataset is studied first by addressing the general characteristics of the elevation changes (Sec. 3.1) and annual accumulation over the whole area (Sec. 3.2), to then focus on the smaller temporal and spatial scales (Sec. 3.3, 3.4 and 3.5). Lastly, 10 we investigate the internal structure of the snowpack deduced from the elevation changes (Sec. 3.6).

Surface elevation changes
The elevation of the surface averaged over the scanned area generally increases with time at a mean rate of 8.7 cmyr −1 (Fig.   2). However, different dynamics are observed for the different years. As noted by Picard et al. (2016) for the first year, a unique and marked accumulation event occurred in the period 4-17 July 2015 which accounted for most of the annual accumulation 15 during that year. In contrast, the second and third years of the dataset feature a fairly regular increase of the surface elevation suggesting the absence of dramatic accumulation events. Interestingly, this regularity in the cumulated precipitation is also found in the ERA-I forecast (Fig. 2, green curve) with a remarkable absence of seasonal signal despite the huge variations of temperature between summer (typ. -25 • C) and winter (typ. -70 • C). However, the precipitation in ERA-I is overall too weak with only 16 kgm −2 yr −1 over the three years (Petit et al., 1982). This mass is equivalent to only 5 cmyr −1 of snow taking a 20 typical surface density value of 320 kgm −3 for the conversion (Genthon et al., 2015;Picard et al., 2014;Leduc-Leballeur et al., 2017). The missing mass in ERA-I could be due to underestimation of the synoptic precipitation or of the condensation that seems indeed very low (Genthon et al., 2015) compared to that reported by other sources (Krinner et al., 2006;Agosta et al., 2018).
The standard deviation of the surface elevation (a.k.a root mean square height of the surface) is shown in Figure 2 as a 25 blue shaded area. Physically, this metric measures the amplitude of the elevation variations over the scanned area including all the spatial scales of variations, that is those due to the meter-scale roughnesses, the local slope and potential instrumental artefacts (such as the repeatability error of the RLS). Overall, the standard deviation has a large magnitude compared to the annual accumulation. Over the three years, the standard deviation was smaller (4 cm) at the beginning of the first season, before the event in July 2015, where it nearly doubled and remained fairly constant after that to 7 cm and then 8 cm from June respective role of the meter-scale roughness and the local slope over the scanned area, we fitted a plane on each scan using the least square method. The local slope is fairly constant, of the order of 1-1.5 • over the time-series, indicating that a large-scale terrain undulation was present. In comparison, the Dome C area has a very small overall slope, less than 1 m per kilometer (0.06 • ). Once the plane is subtracted from the elevation map, the standard deviation is reduced by about half, implying that half of the elevation variations are attributed to meter-scale roughness (typically sastrugi and small dunes) whilst the other half are 5 caused by terrain undulation (e.g. decameter-scale dune Picard et al. (2014)). It is finally worth noting the absence of seasonal signals in the roughness evolution, which differs from reports for other locations on the ice-sheet (Gow, 1969;Adodo et al., 2018).

Spatial and temporal characteristics of the daily accumulation
The mean and standard deviation of accumulation (notedĀ and σ A respectively) between every two consecutive dates (called 25 daily accumulation here despite the irregular sampling in time) were calculated and compared to precipitation (from ERA-I) and wind speed and direction (from the Concordia Automatic Weather Station). Figure 4 shows scatter plots between these variables (except wind direction as no interesting signal was found). Overall, there is little correlation, if any, between the variables, which is unexpected because in many regions of the world the daily accumulation is strongly related to precipitation (e.g. Alpine regions), and in Antarctica the role of the wind has been emphasised (Groot Zwaaftink et al., 2013). This role plot showing the mean and standard deviation of the accumulation. This pattern indicates that significant (positive or negative) mean accumulations are always associated with high standard deviations over the area, which can be expressed mathematical by σ A > |Ā|. We interpret this relationship by the fact that 1) many events change the surface but do not affect the overall mass in the area (σ A >> |Ā| ≈ 0 cm) and 2) that the events leading to significant accumulations or erosions cause highly heterogeneous changes over the area (with both negative and positive changes), and conversely that the process of accumulation 5 or erosion by "layer" (which would correspond to σ A << |Ā|) does not exists at Dome C. However this relationship does not link accumulation to meteorological conditions. This result suggests that further investigation concerning the positive and negative changes of surface elevation in every point is necessary. Figure 5 shows the cumulated sum (average over the scanned area) of daily accumulation considering either only the increments (deposition above 0, 0.5 or 2 cm) or the decrements (erosion higher than 0, 0.5 and 2 cm). The figure also shows 10 the net cumulative accumulation for comparison, which is exactly the same as the mean surface elevation changes in Figure   2. Over the one-year period, each pixel has received a total of 73 cm of snow on average. Most of it has been removed (63 cm, or 87%) leaving a final net accumulation of about 10 cm (Section 3.1). The ratio between deposition and net accumulation suggests that a snow particle is remobilised at least 73/10 = 7.3 times before its definitive settling. This value is to be taken with caution, firstly because of the limited frequency of the scans that only allow to probe the "long term" deposition/removal 15 cycles, which excludes the rapid rebounds that occur during saltation and blowing snow, and secondly because of the limited precision of the laserscanner that leads to count any small measurement errors as positive or negative accumulation. For this latter reason, we additionally computed the increments above 0.5 cm (this value is twice larger than the RLS precision). These "significant" events of deposition bring 56 cm of snow in a year which is still much larger than the net annual accumulation and again suggest that snow is subject to many deposition/removal cycles. Interestingly, the increments above 2 cm (call hereinafter 20 "major events" as they represent more than a 5th of the annual net accumulation) amount to 22 cm, meaning that every pixel receives many times a year at least 2 cm of snow in a single day (or a few days for the few cases the scanner was not working).
The process of erosion shows similar behaviour, with 18.4 cm removed by events larger than 2 cm on average over the area.
Overall, the results of this section depict a surface subject to frequent changes that remains invisible to the observers who only access long time-period averages or spatial averages of accumulation. To further explore the accumulation process, in the 25 next section we focus on the "major events" which account for 30% of the total deposition. Figure 6 shows the accumulation map between 4 and 17 July 2015 corresponding to the exceptional event rising the surface by 17 cm on average and clearly visible in Figure 2. Unfortunately no data is available between these two dates, preventing to describe the precise sequence of this event. This RLS failure was certainly caused by snow jamming the lasermeter window,  Figure 7 shows very different accumulation maps, corresponding to the accumulation between 28 and 29 May 2016 and between 29 and 30 May 2016. It is worth noting that the x-and y-axis scales and the colorscale are different from those in Figure 6. The first day, the net accumulation is -0.8 cm, with 80% of the area in erosion (and 75% in large erosion (>2 cm)).

Patchy accumulation
The second day, the accumulation is +0.3 cm with 62% of the area in accumulation. In both cases, elongated patterns are clearly visible. The second day it is possible to observe that most erosion patterns correspond to the accumulation patterns 5 of the previous day (e.g. the two large blue areas on the left the first day are green the second day) meaning that most of the accumulation features in the first day have been blown away and replaced by new ones with similar shapes (size, elongation, orientation).
To explore the geometrical characteristics of these deposition patterns, that we called patches, we applied a threshold (2 cm) to the daily accumulation maps and segmented the map. Each individual patch received then a unique label and its geometrical 10 properties were computed using the Python function skimage.measure.regionprops. The properties include area A, eccentricity e and major axis length of the best-fitting ellipse (0 indicates a circle and 1 an infinitely-thin segment), and the orientation of the major axis (with respect to North, increasing Eastward). Figure 8 reports the orientation and eccentricity of the elongated patterns only (e > 0.8). On 29 May 2016, it is clear that the patches were aligned with orientations ranging from 140 • to 192 • with respect to the North. The average orientation is 174 • (we applied weighting proportional to Ae 2 to account for the size and 15 eccentricity). This orientation is nearly South which unsurprisingly corresponds to mean wind direction of the day (169 • ) and to the prevailing wind direction (Champollion et al., 2013). The same computation is repeated for the 280 days of the period 2a, yielding 1103 patches that appeared for 93 different days (33% of the time). As found before for the mean accumulation, there is no clear relationship between the patch properties and wind speed (not shown). The most marked relationship is between the daily-mean orientation and the wind direction on the same day ( Fig. 9) with a correlation of 0.33 (the correlation is weighted 20 by the number of patches as in Fig. 9). This simply means that the patches form longitudinal dunes or sastrugi approximatively aligned with the wind direction.
As a consequence of the frequent deposition/removal cycles, many patches studied here are in fact ephemeral. They can be rapidly removed by the next wind event following their formation, as highlighted on the sequence of 29 and 30 May 2016.
These ephemeral patches do not contribute to the snowpack on the long term. To capture and investigate more specifically the 25 few patches that remain and settle and that in the end are the important ones for the snow mass balance and the snowpack internal structure, we estimated the time of survival as follows: for each patch detected at time i, we compute its initial volume (between the surface at time i and the surface at time i − 1) and how this volume evolves at all times i > i (volume between the surface at time i and the surface at time i − 1). More precisely we seek for the minimum value of this volume, and the date at which this minimum is reached. The patch is fully eroded if the minimum volume is 0 cm 3 or negative, which happens for 30 652 of the patches among 1103 detected in the period 2a (59% of the cases). Partial erosion occurs when the minimal volume is between 0 and the initial volume, which concerns all the other patches. We found none being fully and continuously buried after deposition in our dataset. Full or partial erosion of at least half of a patch volume occurs for a large proportion of them (930, i.e. 84%). The mean time of survival before full erosion is 17 days, and for partial erosion this time increases to 46 days, meaning that partial erosion is typically possible after a longer period. In other terms, if a patch is not removed soon after This result confirms that deposition/removal of snow is a very common process and that only a small part of the deposited 5 snow remains on the surface and contributes to the snowpack. Moreover, even when the snow is removed, the number of days spent on the surface can be significant, which lets metamorphism operate and modifies the physical and chemical properties of the snow before it is removed and blown away.

Age of the snow on the surface
We estimate the age of the snow on the surface with the algorithm described in Section 2.2 for all the scans. Figures 10 and 11 10 show the age of snow near the end of the period 2a (18 and 27 January 2017) with respect to the beginning. We have chosen these dates because they are less affected by the fault that increasingly affected RLS from February 2017.
The map and the distribution of age show that snows with very different ages are present on the surface, from 0 (accumulation of the day) to the maximum possible given the duration of the time-series. About 55% of the surface is younger than 100 days (which is thus the approximate median age), while about 18% is older than 200 days and 5% older than 300 days. These 15 proportions change from day to day because of the ephemeral deposition/erosion as illustrated on 27 January 2017 (Fig. 10) where a few patches of recent accumulation covered 11% of the surface. These patches disappeared the day after (data not shown), forming the surface as shown in the 18 January map. It is also clear from the histograms that almost all ages are present at the same time, which is in agreement with the apparent regularity of the accumulation (Fig. 3).
Note that period 2a had a relatively high net accumulation and narrow spatial distribution of accumulation, so it is likely that 20 the distribution of age is even wider for years with less accumulation.
The age maps provide another interesting piece of information. We can indeed see clear patterns with marked alignments on these maps (Fig. 10). This confirms the patchy nature of the accumulation process already noted in the previous section (Sec. 3.4). We indeed showed that a significant part of the daily accumulation over 2 cm is organised in patches but that most of them (59%) do not contribute at all on the long term to the surface mass balance. Here, we complement this by showing that some 25 of the partially eroded patches remain clearly visible on the surface even after a year.

Structure of the snowpack
The accumulation process depicted in the previous section may affect the internal structure of the snowpack and its spatial variability. The burial of the snow can be tracked using the RLS, at least over a short depth given the limited length of the time-series. Figure 12 shows the snowpack internal structure along the x-axis transect at y = 4 m (see e.g. Fig. 10) at the end of (1 Feb 2016). The z origin corresponds to the mean elevation in the whole scanned area at the beginning of period 1. The two grey lines represent the mean at the beginning and end of the period 2a respectively.
It is remarkable that the figure shows distinct coherent patterns and fairly little noise. The figure confirms that the heterogeneity on the surface transfers into the snowpack. For instance, around x = 0 m old snow is present at 11 cm under the surface (yellow) whilst the surface is young (black). This area presents the greatest diversity in terms of age and the largest number Another remarkable feature is around x = 7.5 m where we can see a 10 cm high dune deposited in July 2016 (orange) where the surface was already higher than the surrounding as marked by the grey lines. Then, two successive events (in September and November, violet) accumulated more snow on the sides of this dune which was then about 15 cm higher than the surrounding surface. This may be the result of interaction with the pre-existing dune leading to preferential deposition on a side, here the 15 windward side. A similar behaviour is observed around x = −4 m where the space between the two small dunes (in light orange) has been filled by some subsequent events.

Discussion
The laserscanner that was operating at Dome C for about three years provides very rich and new information on the accumulation process. We analyse these results successively from a temporal, spatial and snowpack perspective.

Temporal perspective
The RLS provides contrasted results regarding the surface evolution depending on the spatial scale of interest. On the one hand, the time-series of surface elevation averaged over the whole scanned area depicts a slow and relatively regular accumulation, without seasonality or rapid events, except a single event which raised the surface by 8 cm in a few days during the first winter of observation. The accumulation rate measured by RLS ranges between 8 and 10 cmyr −1 for the three years of observation, 25 which agrees with values reported by previous studies for Dome C (e.g. Petit et al., 1982;Urbini et al., 2008). The time-series of standard deviation of the surface elevation (RMS height) also is relatively steady, which again suggests that only slow changes occur over the scanned area. This overall stability is further illustrated in the sequence of photographies in Figure   13. The photographies were taken from 20 m above ground at Dome C. They depict a landscape dominated by barchan dunes with little differences between the pictures of January 2017 and December 2017, 11 months apart. This steadiness is surprising 30 because the accumulation over this period has been of the order of 8 cm, the typical annual accumulation expected in the area (Petit et al., 1982;Genthon et al., 2015).
On the other hand, a very different picture is obtained by investigating the daily changes of surface elevation (i.e. the daily accumulation) and the small spatial scales. Indeed, many local changes of elevation affect the surface almost everyday. Most of these changes are however ephemeral (a few tens of days) so that the overall statistics of the surface are relatively constant or slowly changing. These changes are typically caused by migrating patches on their way windwards or remobilisation of loose snow deposited by a recent storm. The picture of September 2017 in Figure 13 illustrates how major these ephemeral 5 changes of the landscape can be. Nevertheless, these temporarly accumulated snow masses have little consequence as far as the surface mass balance and the snowpack internal structure are concerned. They however do have consequences on other aspects as they shield the snow surface for a few days, weeks or even months from solar and infrared radiation and suppress the photochemistry activity and the exchanges of heat, water and chemical components between the consolidated surface and the atmosphere. During the period of residence, these ephemeral snow masses are subject to transformations resulting 10 in changes of their physical, isotopic and chemical properties . When these masses are transported in a downwind region and deposited, they can be very different compared to fresh snow coming from direct local snowfall. This phenomenon of deposition/erosion/transport repeats itself several times. We provide a rough estimate of about 7 cycles before settling, considering only the timescales longer than a day, thus excluding the many rebounds occurring during the saltation.
We however are unable to estimate the distance traveled over these cycles. 15 The meteorological conditions triggering the surface changes is an important question for modelling, but has not been elucidated here. Snowfalls seem to be frequent at Dome C according to the ERA-I reanalysis, which is in agreement with the slow and regular accumulation observed with RLS, but these events are not related to the observed changes of the surface. The occurrence of snowfalls in ERA-I has been compared to satellite data (Palerme et al., 2014;Lemonnier et al., 2018) in a coastal region and found to be reliable, but on the other hand, it is well known that ERA-I misses a large part of the accumulation 20 amount around Dome C (Genthon et al., 2015). Apart from snowfall, common sense suggests that wind speed should be a driver of change. Groot Zwaaftink et al. (2013) considered in a modelling study that snow deposition occurred after sustained wind over 3 ms −1 for 100h. Libois et al. (2014) similarly considered that wind speed over 7 ms −1 is required to initiate drift which then continues as long as this speed limit is reached at least once within 24h. This criteria yields a drift rate of 70 events per year, which compares well with our estimates of 93 days with patch formation over a year. Nevertheless, our analysis has 25 not revealed any robust simple relationship with the wind speed. The accuracy of the daily accumulation derived from RLS is limited by noise and some artefacts, which could be an explanation. However, a complex interaction between wind and the surface is also not to be excluded. For instance, wind with a direction perpendicular to the prevailing surface roughness exerts a stronger drag on the surface than when parallel, potentially resulting in stronger erosion. This effect was at least once shown on the removal of surface hoar at Dome C (Champollion et al., 2013). Similarly (Amory et al., 2016) showed how drag decreased 30 during a snowfall episode (though not at Dome C, but in a coastal region) as the surface roughness direction was adjusting to the wind direction. Another key parameter, here missing in our analysis, is the cohesion of the snow (Sommer et al., 2017) that is able to modulate to a very large extent the snow mobility (Vionnet et al., 2012). The cohesion is due to sintering and increases over time at a rate increasing with temperature. Hence, time elapsed and temperature conditions since the initial deposition could be relevant parameters to be included in a relationship between wind speed and surface changes. However, we found that many patches of accumulation are removed within a few days only, which is probably too short for sintering to be effective (Sommer et al., 2018).

Spatial perspective
All our results show great variability at the meter scale on the surface (e.g. daily accumulation, age of snow on the surface) and within the snowpack (snow layers). Most of this variability results from the heterogeneous accumulation of snow and 5 the selective erosion. We propose to call this former process "patchy" accumulation, in contrast to the even accumulation "in layers" which is typical of the alpine regions. Representing this process in one-dimensional snow models is a challenge. (Groot Zwaaftink et al., 2013) managed to represent the fact that most patches are ephemeral by summing snowfalls and delaying the effective deposition based on wind speed based criteria as aforementioned. Even though this results in thicker deposition per event than when all snowfalls are deposited independently, it is impossible with one-dimensional models to 10 account for the fact that the patches often cover a very small fraction of the surface and therefore can be much thicker than if a snowfall is evenly deposited. This is fundamental to produce layer thickness of a few centimetres, as observed in the field, when most snowfalls bring less than a millimeter per day (4 mm is the maximum of snow over the three year period in ERA-I).
Only using a distributed or three-dimensional approach enables the representation of this feature. Libois et al. (2014) attempts (already mentioned) were a relative success with a good representation of the spatial distribution of annual accumulation 15 compared to the GLACIOCLIM stake network. However, some of their hypotheses need to be reevaluated in the light of our new results. For instance, they assumed that snow deposition only occurs in the 20% lowest part of the surface. This tends to smooth the surface, yet we have shown an example of accumulation near the highest dunes, which conversely tends to enhance the spatial variability. This may be the reason why they obtained an under-estimation of the variability of the density and specific surface area profiles in the snowpack. Further work is needed to implement a process able to produce dunes and more 20 generally to transfer the observational findings of the present study to concrete processing, adequate for numerical modelling.
It results from the patchy accumulation and the strong erosion that the surface is continuously rough at Dome C, as evidenced by the standard deviation of surface elevation. This is usually not the case in alpine regions where roughness increases after snowfalls (Naaim-Bouvet et al., 2016). Surface roughness can amplify itself, first because it plays the role of aerodynamic obstacles that promote heterogeneous deposition and the formation of new rough features overlying the old ones. Moreover, on 25 a longer term, snow on the different faces of the roughness features is exposed to different radiation and wind shear conditions, likely leading to different evolutions of the microstructure (different sintering, sublimation, condensation and metamorphism).
The RLS is limited on this aspect. An avenue is to exploit the laser backscatter signal available from some lasermeters, to retrieve the specific surface area, micro-roughness (hoar), cohesion and potentially other properties.
It is also worth mentioning that the patches studied throughout this paper are smaller in general than the gradient of accu-30 mulation that appeared in July 2015 (Fig. 6) and the barchan dunes visible in the photographies in Figure 13. The former are also well aligned with the wind direction (longitudinal bedform) while the latter have tails elongated at ≈45 • with respect to the wind direction (Filhol and Sturm, 2015). These clearly are different objects, with different size (meter versus decameter) and different dynamics (daily versus yearly).

Snowpack perspective
The heterogeneity of the surface eventually transfers to the snowpack in depth. Figure 14 shows a vertical section of snow extracted from 5 m depth with evidence of past windpacks. Several layers appear distinctively despite the age (over 50 years old).
It is not excluded that transformations of the snow after burial may have amplified the initial differences of snow properties, but in any case, these layers are thick, and were certainly thick when deposited, compared to the annual accumulation. These 5 layers are maybe even thicker than most patches identified with RLS or in the snowpack reconstruction (Fig. 12). An important consequence of this internal heterogeneity concerns the interpretation of ice cores at high resolution or measurements along profiles. With snow age on the surface spanning at least one year, it is clear that some precipitation events, volcanic eruptions, or deposition of nuclear substances are not recorded everywhere, at least not with the same intensity. Gautier et al. (2016) explored in detail this aspect for volcano traces using five cores extracted one meter apart at Dome C. They found that vol-10 canic events were missing in 30% of the cores on average and that the flux uncertainty reached 65% when a single core was used. Another issue is the variable age of the snow at a given depth. We showed variations up to one year, but it was likely under-estimated due to the limited duration of our time-series. This age spread not only hinders the analysis of any annual and sub-annual signals, but also may explain some apparent multi-year signals as discussed by Laepple et al. (2018). This also suggests that depth synchronisation should be applied between different cores. Gautier et al. (2016) found a maximal offset of 15 40 cm between two of their five cores, which is the high-end of what can be explained with our estimates of standard deviation of surface elevation (up to 8 cm) or annual accumulation range (up to 30 cm).

Conclusions
The laserscan dataset collected at Dome C over three years provides for the first time quantitative information on the snow surface dynamics in a site typical of the ridge area on the East Antarctic Plateau. The main results demonstrate that i) the 20 surface elevation increases on average with an apparent regularity, without seasonality, iii) the variations at meter-scales are in contrast large and highly dynamical, iii) the surface is continuously rough, with a standard deviation of up to 8 cm, iv) the accumulation and erosion events are frequent, spatially uneven and significant with respect to the annual net accumulation which implies that snow is remobilised several times before settling, v) the age distribution of snow on the surface spans over more than a year, vi) the snowpack internal structure reflects the surface heterogeneity. These results are useful and 25 have great consequences for several research topics including surface mass balance, surface energy budget and thus climate, photochemistry, snowpack evolution and signals archived in ice cores, which requires further work. We also plan to improve the RLS to capture smaller time scales (hourly) and attempt to increase its robustness in order to collect longer time-series as this proved to be important to assess the age of surface snow and capture the inter-annual climate variability. The present results can also be exploited to build stochastic or physical modelling of the accumulation and erosion process. Investigating 30 other locations on the Antarctic plateau, with different annual accumulation and wind speed is necessary.