Julien Brondex
email: julien.brondex@gadz.org

Olivier Gagliardini

Fabien Gillet-Chaulet

Mondher Chekki

Comparing the long-term fate of a snow cave and a rigid container buried at Dome C, Antarctica

Keywords: Ice Memory, Polar subsurface constructions, Rigid container, Snow cave, Dome C

Ice Memory is an international project aiming at creating a global ice archive sanctuary in Antarctica. The design of a perennial subsurface storage space for the cores is a cornerstone of this project. Here, we use an ice/firn flow model to investigate possible storage solutions that would meet the specific requirements of the project. To this end, we consider two extreme cases in terms of rigidity of the facility: an ice cave dug :::::::: excavated : into the firn and a perfectly rigid container buried within it. We focus on the rate of sinking of the facility as well as on the rate of closure of the cave and the evolution of the normal stresses supported by the container. Our results show that the lifetime of a cave is highly affected by the initial density of snow in its surrounding. On the other hand, the presence of the rigid container within the domain perturbs the flow of snow, creating patches of high density in its surrounding and leading to significant normal stresses on its walls. In particular, strong stress concentrations are obtained at the container angles. These results prove that unreinforced shipping containers are unsuited for this task.

Introduction

Ice Memory is an international project, which aims to create a global ice archive sanctuary in Antarctica gathering ice cores collected all over the world on glaciers that will likely have melted away in the coming decades due to climate change. The design of a perennial facility guaranteeing a safe storage of the ice cores over coming decades to centuries is the cornerstone of this project. Since it is the coldest place on Earth and because several scientific bases already stand at its surface, the interior of the Antarctic ice sheet appears to be the most suitable location for this undertaking. Moreover, burying the storage facility within the polar firn is a simple and efficient way to get a constant temperature that is not affected by seasonal variability, which is essential to preserve the quality of the cores.

After World War II, several subsurface structures aiming at providing yearround accommodation to militaries and/or scientists have been constructed on permanent snow fields (e.g. [START_REF] Mellor | Undersnow structures: Byrd station, antarctica[END_REF][START_REF] Clark | Camp century evolution of concept and history of design construction and performance[END_REF][START_REF] Kohlberg | Georg von neumayer station (gvn) and neumayer station ii (nm-ii) german research stations on ekström ice shelf, antarctica[END_REF]. To this end, various construction techniques have been developed and experimented. They range from the direct burying of hard buildings into the firn to the excavation of unreinforced snow caves and tunnels, including a combination of both with hard buildings placed inside a network of snow tunnels that can be lined or roofed with some hard materials or left unrestrained (e.g. [START_REF] Abele | Production analysis of cut-and-cover trench construction[END_REF][START_REF] Mellor | Methods of building on permanent snowfields[END_REF][START_REF] Steffensen | Report on the neem 2012 balloon trench experiment[END_REF]. However, snow behaves as a viscous fluid flowing under the influence of gravity and the natural fate of any cavity dug ::::::::: excavated into the firn is to close-off. In addition, in regions where the annual surface mass balance is positive, any rigid structure buried into the firn has to withstand an ever-increasing pressure. As a result, although highly variable, the lifespan of the subsurface structures mentioned in the available literature seldom reaches the decade. Beyond this time window, severe distresses are usually reported, including failure of roof trusses, crushing of lined tunnels and severe closure of unrestrained tunnels [START_REF] Mellor | Methods of building on permanent snowfields[END_REF][START_REF] Kovacs | Camp century revisited, a pictorial view -june 1969[END_REF].

Beside the construction technique used, the lifespan of a subsurface structure depends on the snow accumulation rate, the initial depth of the construction, the initial size of the excavations, the vertical profiles of temperature and density, the potential occurrence of melting event, the flow regime at the considered location, and possibly other parameters not identified yet. In this context, Dome C is probably one of the most favourable locations for this type of construction because of the very ::::::::::: combination ::: of : specific conditions prevailing there: very low snow accumulation, very dry atmosphere, very low temperature making the occurrence of melting events unlikely in the coming decades, and almost purely vertical motion of snow making shear stresses negligible. In addition, since the facility aims at storing ice cores only, it will not contain any internal heat source, which has been shown to strongly increase the closure rate of snow cavities [START_REF] Clark | Camp century evolution of concept and history of design construction and performance[END_REF]. On the other hand, regular maintenance of the installation should not be considered as a possible mean to increase its lifetime. Indeed, we must consider the possibility that a continuous deployment of staff at Dome C could be hampered in the future. All these particularities strengthen the need for an ad hoc numerical study aiming at investigating optimal solutions to meet the specific needs of this project. This numerical study is the subject of the present paper.

Snow is a complex material and mathematical laws implemented in numerical models to account for its mechanical behaviour cannot claim to capture the full complexity of underlying physical processes [START_REF] Arthern | In situ measurements of antarctic snow compaction compared with predictions of models[END_REF]. In addition, these laws involve many parameters, some of which are poorly constrained or affected by significant spatial or temporal variability. As a consequence, the first step of the present study consists in constructing an initial steady state of the considered domain, which is then compared to available field measurements.

Starting from this initial state, we investigate the mechanical interactions between the storage structure and the surrounding flowing firn, considering two end-member cases in terms of rigidity of the structure: an unreinforced snow cave dug :::::::: burrowed : into the firn and a perfectly rigid container buried within it.

We put a particular focus on the rate of sinking of the storage facility as well as on the rate of closure of the cave and the temporal evolution of the normal stresses supported by the container. In Sect. 2, we introduce the model used to conduct this study and present the experimental setup. The results obtained for each experiment are presented in Sect. 3 and discussed in the last section.

Methods

Model description

The mechanical interactions between the storage solution and the surrounding compressible firn are modelled using the code Elmer/Ice [START_REF] Gagliardini | Capabilities and performance of elmer/ice, a new-generation ice sheet model[END_REF]. Elmer/Ice is an open-source finite-element software for ice sheets, glaciers and ice flow modelling, which also includes a module for the rheology of compressible firn [START_REF] Gagliardini | Flow simulation of a firn-covered cold glacier[END_REF] that has already been applied in several studies [START_REF] Zwinger | A full Stokes-flow thermo-mechanical model for firn and ice applied to the Gorshkov crater glacier, Kamchatka[END_REF][START_REF] Gilbert | A 3-D thermal regime model suitable for cold accumulation zones of polythermal mountain glaciers[END_REF][START_REF] Licciulli | A full stokes ice-flow model to assist the interpretation of millennial-scale ice cores at the high-alpine drilling site colle gnifetti, swiss/italian alps[END_REF].

All variables and parameters used in this study are summarized in Table 1.

Flow law

We adopt the flow law first proposed by [START_REF] Gagliardini | Flow simulation of a firn-covered cold glacier[END_REF] and later corrected by [START_REF] Zwinger | A full Stokes-flow thermo-mechanical model for firn and ice applied to the Gorshkov crater glacier, Kamchatka[END_REF] to model the flow of firn, which behaves as a non-linear viscous compressible fluid. It ::::

This ::: law :: is :::::::: intented ::: to ::::::::

represent :::: the ::::::::: secondary :::::: creep :: of ::::: firn. ::: As :::::: such, ::::::::: processes :::::: related ::: to :::::::: primary :::::

creep :::: and ::::: snow ::::::::::::: metamorphism :::: are ::: not :::::::::: accounted :::: for. ::: In :::::::: addition, :: it : must be stressed that this approach does not capture brittle fracture of snow, which involves other mechanical processes. In particular, the present model is limited to the description of the continuous deformation of a snow cave related to the flow of firn and cannot provide any information on a potential collapse of the latter. Nevertheless, none of the reports regarding subsurface constructions on permanent snow fields available in the literature mention such an event (e.g. [START_REF] Mellor | Methods of building on permanent snowfields[END_REF][START_REF] Kohlberg | Georg von neumayer station (gvn) and neumayer station ii (nm-ii) german research stations on ekström ice shelf, antarctica[END_REF][START_REF] Steffensen | Report on the neem 2012 balloon trench experiment[END_REF].

The Cauchy stress tensor σ can be decomposed into an isotropic part -pI, where p = -(tr σ)/3 is the isotropic pressure (tr denotes the trace operator)

positive for compression and I the unit tensor, and a traceless, deviatoric stress tensor, i.e.

τ = σ - tr σ 3 I = σ + pI . (1)
Similarly, the deviatoric part ė of the strain-rate tensor

˙ = 1/2[grad u + (grad u) T ],
where u is the velocity vector, is obtained by

ė = ˙ - tr ˙ 3 I = ˙ - div u 3 I . ( 2
)
Invariants for the strain-rate can then be defined as

γ 2 e = 2 tr ( ė) 2 = 2 ėij ėij , ˙ 2 D = γ 2 e a + (div u) 2 b , (3)
and for the stress as

τ 2 = 1 2 tr (τ ) 2 = 1 2 τ ij τ ij , σ2 D = aτ 2 + bp 2 . ( 4
)
The two functions a = a(D) and b = b(D) depend only on the relative density D = ρ/ρ i , where ρ is the snow density and ρ i is the ice density. Following [START_REF] Gagliardini | Flow simulation of a firn-covered cold glacier[END_REF], the relationships between the deviatoric and isotropic parts of the stress and strain-rate tensors write, respectively,

τ = 2 a B -1/n ˙ (1-n)/n D ė , (5)
and

p = - 1 b B -1/n ˙ (1-n)/n D div u . ( 6
)
The parameter n used in Eqs. ( 5) and ( 6) is the flow law exponent, usually set to n = 3, while B is a fluidity parameter which depends mostly on the temperature of the snow/ice. The firn temperature at Dome C is affected by the seasonnal variability of surface temperatures only within a 10 m-thick layer. Below this surface layer, the firn has a constant temperature equal to the annual mean surface air temperature, which is estimated at -55 • C (Leduc- [START_REF] Leduc-Leballeur | Modeling L-Band Brightness Temperature at Dome C in Antarctica and Comparison With SMOS Observations[END_REF]. For this reason, following recommendations of [START_REF] Cuffey | The physics of glaciers[END_REF], the fluidity parameter B is set to a uniform value of B = 0.078 MPa -3 a -1 . Regarding the functions a and b, we use the analytical solution proposed by [START_REF] Duva | Analysis of consolidation of reinforced materials by power-law creep[END_REF] for high relative densities (0.81 < D < 1.0): Relations ( 8) and ( 9) have been shown to produce acceptable results for simulated density profiles and surface velocities at Col du Dôme, Mont-Blanc, France [START_REF] Gilbert | A 3-D thermal regime model suitable for cold accumulation zones of polythermal mountain glaciers[END_REF]. It is important to note that the parametrizations adopted for these two functions is a source of uncertainty on the simulated density profiles as well as on flow velocities and, therefore, on the computed stresses. Ideally, these functions should be re-calibrated on a case-by-case basis for each particular applications. However, this would be a time-consuming task requiring a lot of in-situ measurements for a limited benefit since it would anyway not allow to capture the high spatio-temporal variability of the density as measured over the first few tens of meters of the polar firn (e.g [START_REF] Leduc-Leballeur | Modeling L-Band Brightness Temperature at Dome C in Antarctica and Comparison With SMOS Observations[END_REF]. For this reason, we have used for the present work the relations ( 8) and ( 9) without further modifications.

a 0 (D) = 1 + 2(1 -D)/3 D 2n/(n+1) , b 0 (D) = 3 4 (1 -D) 1/n n 1 -(1 -D) 1/n 2n/(n+1) (

Field equations

To keep the computational cost affordable, we solve the full Stokes equations, with the constitutive law of firn described above, on a two-dimensional rectangular domain in a (x, z) vertical plan, the z axis being the vertical pointing upwards. The momentum balance, in which acceleration terms are neglected,

reads div σ + ρg = 0 , ( 10
)
where g = (0, -g) is the gravity vector. The spatio-temporal evolution of density is governed by the mass conservation equation:

∂ρ ∂t + div ρu = 0 . ( 11
)
Equations ( 10) and ( 11) are coupled to an advection equation governing the evolution of the domain top surface, which is a free surface. This equation reads:

∂ t z s + u∂ x z s = w + a s , (12)
where ∂ i z s denotes the partial derivative of the function z s relative to the variable i, z s is the top surface altitude, w the vertical component of the velocity vector and a s the annual surface mass balance. Although snow precipitations in Antarctica are expected to increase slightly in the future [START_REF] Palerme | Evaluation of current and projected antarctic precipitation in cmip5 models[END_REF], uncertainties regarding the local evolution of surface mass balance are high and the value of a s is kept constant in all simulations of the present study.

Boundary conditions

The top surface being a stress-free surface, the following Neumann condition applies:

(σ • n)| zs = 0 , ( 13
)
where n is the normal unit vector to the considered boundary. The domain is assumed to be a perfect dome, which implies purely vertical flow in the absence of any obstacle and yields the following Dirichlet condition on both sides:

(u • n)| x l = (u • n)| xr = 0 , (14)
where x l and x r are, respectively, the left and right horizontal boundaries of the domain. The bottom boundary condition is based on the assumption that, on the time scales at stake in our simulations, the altitude of the surface at Dome C is constant. As a consequence, all the snow accumulated at the top surface over an arbitrary period of time must be balanced by an equivalent mass of snow/ice flowing out of the domain through its bottom boundary over the same time period, which yields the condition:

(u • n)| z b = a s ρ(z = z s ) ρ(z = z b ) , ( 15
)
where z b is the altitude of the bottom boundary, which is kept fixed in time.

Finally, the resolution of Eq. ( 11) requires a Dirichlet condition on the top boundary, which is the only inflow boundary of the domain. This condition is simply given by:

ρ(z = z s ) = ρ s , ( 16
)
where ρ s is the density of fresh snow. initialisation step is shown in Fig. 1a. Its horizontal extension is arbitrarily set from x l = 0 m to x r = 22.44 m. The sensitivity of the results to this choice is discussed in the Supplement. Regarding the vertical dimension, it must be stressed that the ice thickness at Dome C is estimated at :: of 3273±5 m (Parrenin et al., 2007), while, on the considered time scales, the storage facility is expected to evolve within the first few tens of meters below the surface only. The choice of the vertical extension is then a trade-off between a restricted domain size limiting the computational cost of simulations and a domain deep enough so that the sensitivity of the results to the bottom boundary condition (15), which is poorly constrained, is limited. As a consequence, the vertical limits of the initial domain are set to z b = -250 m and z s = 0 m.

Description

The resolution of Eq. ( 11) requires the prescription of an initial density field.

We use the parametrized function derived by [START_REF] Leduc-Leballeur | Modeling L-Band Brightness Temperature at Dome C in Antarctica and Comparison With SMOS Observations[END_REF], which corresponds to the best fit in the least squares sense to density measurements performed on two 80-m-long ice cores drilled at Dome C in the austral summer 2012-2013. This parametrization is also used to calculate the values of densities at the top and bottom boundaries of the domain as required by conditions (15) and ( 16).

For this simulation, we use an unstructured mesh made of 14656 threenode triangular elements of uniform size. The typical spacing between two neighbouring nodes is of ∼ 1 m. To be sure that a steady solution is reached at the end of the simulation, the prescribed total simulation time is of 10 ka and the timestep size of 1 a.

Snow cave experiment

A well-proven technique to build a snow cave into the polar firn relies on the use of an inflatable balloon. First, the balloon is placed in a trench excavated in the snow. Then, the balloon is inflated and snow is blown back on the top of it until the trench is filled up. After sufficient age hardening of snow has occurred (typically after a few days), the balloon is deflated and removed from the resulting snow cave. This technique is already mentioned in the report of [START_REF] Mellor | Methods of building on permanent snowfields[END_REF], but is not detailed very much. A recent implementation of this concept was performed in 2012 by members of the Center for Ice and Climate from the University of Copenhagen for the NEEM ice core drilling project in North-West Greenland [START_REF] Steffensen | Report on the neem 2012 balloon trench experiment[END_REF]. Careful monitoring of the deformations of this construction showed encouraging results, leading the Ice Memory steering committee to initiate its own field test at Dome C [START_REF] Ascione | Balloon cave construction test site at concordia station, 18/19 summer campaign[END_REF].

The numerical setup that we are presenting here, and in particular the initial depth of the modelled cave, is based on this field test (Fig. 1b). The transverse section of the cave corresponds to the combination of a square that is 5 m on each side and of the incircle of this square: the top half of the section is thus a half circle of diameter 5 m while the bottom half is a rectangle. The highest point of the cave is located 2.2 m below the firn surface, which is perfectly flat at the beginning of the simulations. As for the initialisation, the mesh is made of three-node triangular elements. However, for this experiment, the mesh is not uniform. High refinement is prescribed within a rectangle surrounding the cave where the typical spacing between two neighbouring nodes is of ∼ 5 cm.

Outside of this highly refined rectangle, the elements size increases following a geometric progression towards the boundaries of the domain, such that the typical spacing between two neighbouring nodes at the bottom boundary of the

domain is of ∼ 1.2 m.
The cave/firn interface is treated as a stress-free surface, which will closeoff over time due to the flow of firn. From a mechanical point of view, this is equivalent to state that the cave/firn interface has no rigidity at all. At each time step, the displacement δx k of a node k belonging to this interface is given by:

δx k = u k δt x + w k δt z , ( 17
)
where u k and w k are, respectively, the horizontal and vertical components of the flow velocity computed at node k, and δt is the timestep size. To deal with these displacements, we use a module available in Elmer/Ice which displaces the nodes through a deformation of the whole mesh, with the constrain of keeping the nodes located on the lateral and bottom boundaries of the domain at fixed positions. However, a problem arises after a few tens of years of simulation, when the mesh is too deformed with elements overlapping each others. For this reason, an automatic remeshing procedure was developed. Every 20 a, the run is stopped and a new mesh of the geometry in its current state is created. For this new mesh, the elements size is a function of the distance relative to the firn/cave interface, which has moved and deformed over the course of the run.

This function is parametrized in order to get small elements of uniform size, with a distance between two neighbouring nodes of about 5 cm, within the first 10 cm surrounding the firn/cave interface. As going further away, the size of elements increases linearly. All the physical fields computed at the end of the 20 a run preceding the remeshing step, in particular the velocity and density fields, are then linearly interpolated on the new mesh, which constitutes the initial state

for a new 20 a transient run. This procedure is repeated until the end of the simulation. The total simulation time is set to 150 a as the deformation of the cave becomes too high beyond this time window. The consecutive remeshing steps necessarily induce a loss of information through interpolation diffusion between two consecutive runs, especially in places where the considered fields are strongly non-linear. However, since the mesh is always highly refined in the area of interest, it turns out that the final results are not significantly affected by this procedure.

The snow cave experiment includes four simulations which only differ by the initial density prescribed in the close vicinity of the cave. For the reference simulation, the initial density field corresponds to the steady density field obtained at the end of the initialisation run. For the three other simulations, this initial density field is slightly modified to account for the trench involved in the construction process of the cave, which is backfilled with snow of higher density by a snow blower. In order to assess the importance of the trench size and shape on the lifetime of the cave, we investigate three cases: (1) a single trench with a width corresponding exactly to the balloon diameter (narrow trench case),

(2) a single trench which is two meters larger than the balloon on both sides (wide trench case), and (3) a T-shape trench with 1 m wide ledges on both sides of the cave mimicking the on-going Dome C field test mentioned above. The exact dimensions and shapes of the trenches involved in these three cases are reported in Fig. 1b. The sensitivity to the width of the top horizontal branch of the T-shape is investigated by running two additional simulations for which this branch is, respectively, 2 and 3 m wider than the cave on both sides. Inside the trench, the initial density of blown snow is forced to ρ(x, z, t = 0) = 550 kg m -3 , which corresponds to the density measured after backfilling at Dome C. Outside the trench, the initial density is again the steady density field obtained at the end of the initialisation. Note that the sensitivity of the lifetime expectancy of the cave to its initial geometry has also been evaluated by running four similar simulations as the ones presented here, but considering a perfect circle for the initial shape of the cave (Supplement).

Rigid container experiment

In this experiment, we consider the opposite end-member in terms of rigidity of the storage solution: a perfectly rigid container buried within the firn. The chosen dimensions for the container width and height, i.e. w c = 2.44 m and h c = 2.9 m respectively, are based on the dimensions of the 20ft shipping containers commonly operated in Antarctica by the French Polar Institute (Institut Paul-Emile Victor, IPEV). The container floor is initially located at z = -10 m. This initial depth is motivated by the fact that below 10 m, the firn temperature is constant, which is preferable for the long-term quality of the cores.

Imposing a perfect rigidity of the container is equivalent to impose a no flux condition at the firn/container interface. However, because the container is advected with the flowing firn, the velocity of firn at this interface is part of the solution to the flow problem. In other words, the resistance opposed by the container to the flow of firn translates into an implicit Dirichlet boundary condition for the flow problem, defined on a moving boundary which corresponds to the firn/container interface. As commonly done when dealing with flows around moving rigid bodies (e.g. [START_REF] Glowinski | A distributed lagrange multiplier/fictitious domain method for flows around moving rigid bodies: application to particulate flow[END_REF], the vertical velocities of all nodes belonging to the container roof (i.e. normal velocity to the roof) are forced to be equal to each other through the use of Lagrangian multipliers.

In addition, we impose periodic boundary conditions for the firn normal and tangential velocities between the container roof and floor, while the horizontal velocity (i.e. normal velocity) is forced to zero on the two container sides.

There are no imposed conditions regarding the tangential velocity of firn on the container sides, which means that a relative motion between the snow and the container is allowed over these two interfaces. However, another possibility is that the snow sticks to the container sides. In that case, the boundary condition which needs to be applied to the container sides is a no slip condition. In between these two end-members, it could also be assumed that there is some kind of friction between the container sides and the snow flowing around. In that case a friction law should be implemented as the boundary condition. Since such a law is difficult to constrain, we have made the choice to run all the simulations of this experiment with a free slip boundary condition on the container sides.

The sensitivity of the results to this choice turns out to be low, as discussed in the Supplement.

As for the snow cave experiment, the initial mesh of the domain is made of three-node triangular elements, and is characterized by a high refinement around the container which decreases linearly as going further away. The container being perfectly rigid, the container/firn interface undergoes a rigid body displacement. Concretely, this is achieved by extracting the vertical flow velocity of firn at the nodes belonging to the roof of the container. The vertical displacement of all nodes belonging to the container/firn interface over one timestep is given by the product between the vertical flow velocity computed at this particular node and the size of the timestep, while the horizontal displacement of these nodes is forced to zero. As in the snow cave case, the whole mesh deforms as the container sinks, and needs to be re-designed every 20 a of simulation through the same procedure as the one described previously.

The rigid container experiment gathers a total of six simulations, corresponding to a combination of three different initial densities and two different container weights. Because the container is perfectly rigid, the total simulation time is not limited by its deformation and is set to 200 a for all the simulations.

A reference simulation is run for which the initial density field is the steady density field obtained at the end of the initialisation run. The second and the third simulations include a trench of initial density ρ(x, z, t = 0) = 550 kg m -3 , which has the same width as the container for the former and which is two meters larger than the container on both sides for the latter. These first three simulations do not consider the weight of the container itself or of the payload it contains. Therefore, we run three similar experiments which differ only by the fact that the weight of the container and its payload is accounted for. To this end, we consider a total weight of 26 tons, which corresponds to a container tare weight of about 2 tons and a payload weight of about 24 tons. This total weight is divided by the total length of the container, i.e. l c = 6 m, since only a vertical section of the container is considered. Finally, the weight is applied as a uniform pressure spread on the container floor.

Results

Initialisation

Figure 2 shows the evolution of the computed relative density and vertical velocity profiles over the course of the initialisation run. Although the vertical density profile used to initialize the model is derived from measurements, the model is initially out of equilibrium. This is first the consequence of the poor constrain prevailing on several model parameters, and most notably, on the parameters at stake in the two functions a and b used in the constitutive laws ( 5) and ( 6). Second, the prescribed initial density profile is itself poorly constrained.

Indeed, the density measurements performed on cores drilled at Dome C show very high variability over the first 20 to 30 meters below the firn surface. As a consequence, the prescribed initial density profile is poorly representative over this surface layer. Moreover, this profile is derived from measurements performed on two 80-m-long ice cores and it turns out that, beyond this depth, it underestimates densities compared to deeper measurements performed on other cores drilled at Dome C (Fig. 2a).

Initial density profile @10a @100a @1ka @10ka @0a @10a @100a @1ka @10ka (
After 10 a of simulation, the vertical velocity profile shows an inflection point at a depth of about 50 m below the firn surface. Above this depth, vertical velocities increase rapidly towards the firn surface to reach about 41 cm a -1 at the surface. Below this depth, the vertical velocities are of a few cm a -1 only and decrease slowly as going deeper to reach the value imposed by boundary condition (15), i.e. 2.87 cm a -1 , at the bottom of the domain. This flow regime induces significant compressive stresses which ease the densification of the low density surface layers as they are advected downward. Thus, while the modelled densities are still very close to the initial profile at depth (as well as at the surface because of boundary condition ( 16)), in a depth range comprised between 5 m and 45 m below the firn surface they are getting slightly higher than measured densities. After 100 a of simulation, the vertical velocity has considerably decreased over the whole firn column, and particularly within the top 50 m below the firn surface. The vertical velocity at the surface is of about 16 cm a -1 at 100 a. The firn density at a given depth keeps increasing. In particular, by that time, the modelled density of the deepest layers of the domain is getting higher than that of the initial profile. Beyond 1000 a of simulation, the modelled density at a given depth is systematically in the higher range or higher than the measured ones. Between the penultimate and the last timestep of the initialisation run, the density field has evolved by less than 0.001%, and we consider that the density profile obtained at 10 ka corresponds to a steady state.

Given the high spatial variability of measured densities over the top 20 to 30 meters, the initial density profile is hardly more consistent with measurements over this surface layer than this steady density profile. In addition, the first computed surface vertical velocity is of about 41 cm a -1 , which is largely above the observations. In contrast, as expected from the implemented boundary condition (15), the computed vertical velocity at the surface for the steady state corresponds to the observed mean annual surface accumulation, i.e. slightly less than 8 cm a -1 (Parrenin et al., 2007).

Snow cave experiment

Figure 3 shows the evolution of the shape and the position relative to the surface of the cave over the course of the experiment for each of the considered initial density fields. Although not represented in Fig. 3, a few centimer-high bump tends to form at the firn surface in all cases, except for the reference run for which a few centimeter-deep cavity appears from the first years of the simulation.

These surface deformations can be seen on the animation of the full simulations included in the Supplement. The deformation of the cave is strongly sensitive to the initial density in its immediate surrounding. First, the deformation patterns of the roof are different when the trench is included than when it is not. Indeed, for the reference simulation, the shape of the roof evolves progressively from concave to convex, inducing a close-off of the cave that is much faster than for the three other cases. In contrast, for all the considered cases, the lateral walls and floor show similar deformation patterns with progressive curvature toward the interior of the cave. Furthermore, the magnitude of the deformation also depends on the initial density field: the wider the high density trench, the less readily the cave closes-off. These different magnitudes of deformation can be explained by the fact that a patch of high density snow surrounding the cave enables stress transfer around the cave to the underlying firn layers. Thus, apart from the reference simulation for which the cave rapidly becomes unusable, the close-off is slow, with a height reduction of the order of ∼ 1 m after 50 a of simulation for the T-shape trench case and even slightly less for the large trench case.

The evolution over time of the position of the cave pseudo-center, which is defined as the middle of the segment joining the highest and lowest points of the cave belonging to its central axis, is represented in Fig. 4a. After 139 a of simulation for the reference run, the cave is too deformed and the pseudo-center of the cave is no longer defined. Note also that the fact that the rate of sinking seems much higher for the reference simulation than for the other runs is actually due to the strong asymmetrical deformation of the cave which tends to shift the pseudo-center of the cave downward. Similarly, although deformations are much smaller for all the other considered cases, the magnitude of the roof deformation is systematically a bit more important than that of the cave floor, which again tends to shift slightly the position of the cave pseudo-center downward. With this in mind, Fig. 4a shows that, unsurprisingly, the initial density field has little influence on the rate of sinking of the cave, which is of the order of ∼ 5 cm a -1 for all the considered cases. This rate of sinking is similar to the rate of sinking of a snow particle in the unperturbed initial steady density field (orange line in Fig. 4a which is almost superimposed on the blue line).

From results reported here, it appears that a perfectly rectangular trench larger than the cave is the best option to maximize the lifetime of the latter.

However, for practical reasons related to the construction process, the T-shape is the most convenient shape for the trench as it enables to keep the balloon in place while snow is blown back on it. For this reason, two additional simulations considering a T-shape trench with a horizontal branch, respectively, 2 and 3 m wider than the cave on each sides have been run. It turns out that the wider the horizontal branch of the T-shape, the longer the lifetime of the cave. However, the gain in terms of lifetime between a branch that is 1 m larger than the balloon diameter on each side and one that is 2 m larger is more important than the gain obtained when the branch is 3 m larger instead of 2 m larger. In other words, beyond a certain width of the horizontal branch of the T-shape, the gain in terms of lifetime of the cave is not sufficient to justify the amount of work required to make the branch larger. Our results suggest that a branch being 2 m larger than the balloon diameter on each side is a good compromise for a balloon with a 5 m diameter.

Rigid container experiment

Figure 4b shows the evolution of the position of the roof of the container over time for the six main simulations. As for the snow cave experiment, the rate of sinking of the container is not very sensitive to the initial density field in its surrounding. It also shows little sensitivity to weight consideration: for each case of initial density, the version of the container for which the weight is accounted for sinks hardly faster than the unloaded version. Overall, the rate of sinking of the container is similar to that of the snow cave, i.e. of the order of ∼ 5 cm a -1 . As the container sinks, the density of snow in its immediate surrounding increases. Figure 5 shows the evolution over the course of the simulation (one plot every 20 years) of the snow relative density along the container roof for the six considered cases. The same figure for the container floor is included in the Supplement (Fig. S1). The progressive densification of snow observed on the container roof is obviously partly due to the fact that snow is denser at higher depths, but also to the apparition of a patch of higher density on the container roof as snow encounters the container while flowing downward.

Within this patch of higher density, firn velocities are slightly lower than at the same depth further away from the container, in places where the density field is not perturbed by the presence of the latter. This explains why, when with weight) large and narrow trenches simulations, we report the results obtained for the no weight (resp. with weight) reference simulation to ease comparison (black dotted lines).

compared to the sinking of a snow particle flowing in the unperturbed initial steady density field (orange line in Fig. 4b), the sinking of the container is systematically very slightly slower than the latter. This patch of higher density is not very sensitive to the weight of the container and is very similar for the narrow and large trench cases. In contrast, the reference simulation (for both the unloaded and loaded cases) shows lower snow density on the container roof over the whole simulation, which is not surprising since, by construction, it starts with a much lower density at the firn/container interface than when a high density trench is included. The patch of higher density also appears in the vicinity of the container floor (Fig. S1). At this location, the sensitivity to the initial density field around the container as well as to whether the weight is considered or not is stronger. The relative density along the container floor in the case of the unloaded reference run does not depart very much from the unperturbed relative density. In contrast, when the weight of the container and its payload is accounted for and/or when the high density trench is included, the modelled relative density along the container floor turns out to be higher than the corresponding unperturbed relative density at the same depth. This is because the weight of the container and/or the weight of the dense snow of the trench above the container is transferred to the underlying snow layers (the container being perfectly rigid) and tends to compress them, enhancing their densification. A remarkable feature concerning every considered cases -although less marked on the container floor for the unloaded reference simulation -is the occurrence of very localised peaks of high relative density at each corner of the container roof and floor. At these places, the relative density is significantly higher than that at the middle of the container roof/floor. Figure 6 shows the evolution over time (one plot every 20 years) of the normal stresses exerted by snow along the container roof for the six considered cases.

The same figure for the container floor is included in the Supplement (Fig. S2).

Note that all these stresses are negative because they are compressive stresses.

The patterns of these normal stresses are correlated to the previously described patterns of relative density. In particular, the peaks of high density obtained at corners at both container floor and roof are correlated to strong stresses concentrations in these same areas. Further away from the corners, normal stresses are up to four times lower. Apart from the two reference simulations which systematically produce lower normal stresses, the order of magnitude of the normal stresses supported by the middle parts of the container roof/floor are mostly independent of the initial density around the container and whether

or not the weight is considered. The stress concentrations at corners are more sensitive to these two parameters, but these values must be considered with caution since they are much more sensitive to the numerics, including mesh refinement and linear interpolation following remeshing steps. At the end of the simulation, the order of magnitude of the normal stresses experienced by the middle parts of the container roof/floor is of ∼ 120 kPa, whereas it reaches up to ∼ 450 kPa at floor corners where the stress concentrates, and slightly less, i.e. up to ∼ 400 kPa, at roof corners.

In order to quantify the gap between these normal stresses and ::::::: evaluate :::

the :::::::::: importance ::: of ::::: using :: a ::::::: complex ::::: flow :::::: model :: to :::::::: estimate ::: the :::::: loads ::::::::: supported ::

by :::: the ::::::::: container ::::::: instead ::: of :::::: simply ::::::::: assuming ::::: that ::::: these :::::: loads ::: are ::::::: limited ::: to the hydrostatic pressure prevailing at the considered depths, ::::::::: calculated :: at :::: any shows : the temporal evolution (one plot every 20 years) of the ratio between these two quantities is represented in Figure 7 :::::::: modelled ::::::: normal ::::::: stresses :::: and ::::::::::

hydrostatic :::::::: pressure : for the six considered cases. The hydrostatic pressure is calculated as:

P static = -ρ mean gH , ( 18
)
where g is the gravity, H is the height of snow above the container roof at the considered time of the simulation, and ρ mean is the mean density integrated from the firn surface to remaining :::::: ∼ 35% :: is :::: due :: to ::: the ::::::: gradual :::::::::::: densification :: of ::::: snow :: in ::: the ::::: close ::::::: vicinity ::

of ::: the ::::: roof. :::: For :::: the :::: two ::::: cases ::::::::: including :: a ::::::: trench, :::::::::: dynamical ::::::: stresses ::::::: explain ::::

only :::::: ∼ 30% ::: of :::: this :::::::::: difference, ::: the :::::::::: remaining :::::: ∼ 70% :::::: being ::::::: mostly :::: due :: to :::: the :::::::

presence ::: of :::: the :::::: trench ::: of :::: high :::::: initial :::::::: density :::: but :::: also :::::: partly ::: to ::: the :::::::: gradual ::::

snow ::::::::::::: densification. ::: It :: is ::::::::::: interresting ::: to ::::: note :::: that :::::::::: significant ::::::::: deviation ::::: from ::::::::: lithostatic :::::::: pressure :: at ::: the :::::::: interface :::::::: between :::: two :::::: bodies :::: with :: a :::::::: viscosity ::::: ratio :: of ::

10 :: or ::::: more :::: has :::: also ::::: been :::::::: reported :: in :::: the ::::::: context :: of ::::::::::: geodynamic ::::::::::::: reconstruction :::::::::::::::::::::::::::::::::::: [START_REF] Moulas | Stress field associated with elliptical inclusions in a deforming matrix: Mathematical model and implications for tectonic overpressure in the lithosphere[END_REF][START_REF] Luisier | Metamorphic pressure variation in a coherent alpine nappe challenges lithostatic pressure paradigm[END_REF]. Regarding the floor center part, the ratio between σ nn and P static is slightly lower, with values around 1 when the trench is not included and around 1.3 otherwise (Fig. S3) : , :::: this :::::: latter ::::: value :::::

being ::::: again :::::: partly ::: due ::: to ::: the :::::::::::::::: non-consideration :: of ::: the ::::: high :::::: density ::::::: trench :: in ::: the :::::::::: hydrostatic :::::::: pressure :::::::::: calculation. ::: In ::::: other ::::::: words, ::: the ::::::: normal ::::::: stresses :::::: along ::: the ::::::::

container ::::: floor ::: are :::: less ::::::: affected ::: by ::: the :::: flow ::::: than ::::: those ::::: along :::: the :::::::: container :::: roof.

Note that this ratio does not evolve very much over the course of the simulation, especially beyond the first 20 years. As stated above, values obtained at corners must be handle with care as they are more sensitive to the numerics. Yet, results

presented here tend to show that normal stresses concentrated at roof and floor corners are of the order of 4 to 6 times higher than corresponding hydrostatic pressures in an unperturbed density field, :::::::::::: emphasizing ::: the :::::::::: imperative ::::: need ::: for : a :::: flow ::::::: model :: in :::::: order ::: to ::: get :::::::: realistic ::::::::: estimates ::: of ::::: snow :::::: loads, :::::: which :::: are :: a ::::::::::

prerequisite ::: to ::: the :::::: design :: of :: a ::::::: reliable ::::::: storage :::::::: solution.

Figure 8 shows the evolution over time (one plot every 50 years) of normal stresses exerted by snow along one container side for the six considered cases.

Note that the experimental setup being perfectly symmetric by construction, the results obtained on the right and left sides of the container are almost identical (slight differences occur over time, especially in the narrow trench case, because of the remeshing procedure and induced linear interpolation but are negligible).

The normal stresses on container sides turns out to have very low sensitivity to the consideration of weight. Indeed, if tiny differences obtained between the loaded and unloaded versions for each of the three initial density cases exist, they are mostly attributable to the small differences in the sinking rate, with the loaded cases inducing slightly deeper containers than the unloaded cases at each of the considered times. Another remarkable feature is the low, or even null, normal stress localised on the top corners of the container sides. This is because the snow tends to detach from the wall of the container in this area. However, the contact between snow and the container wall is firmly restored just a few centimetres lower. The maximal normal stress supported by the container wall occurs just a few centimetres above the bottom corners of the container sides.

It is interesting to note that :::: even ::: for :::: the :::::::: reference :::::::::: simulation : these maximal stresses are lower than the ones supported by the container roof which confirms that normal stresses depart from purely hydrostatic loads, in which case deeper snow layers would necessarily induce higher pressure on the container.

Discussion

::::::

Among :::: the ::: few ::::: polar ::::::::: subsurface :::::::::::: constructions :::::::: reported :: in :::::::: available ::::::::: literature, entirely :::::: buried ::: by ::::::::: ∼ 180 cm :: of ::::: snow ::: by ::::: 1960 ::: due ::: to ::::: snow :::: drift :::::::::::::: [START_REF] Mellor | Methods of building on permanent snowfields[END_REF]. after ::: its ::::::::::: construction, :::::::: pictures :::: were :::::: taken ::::: inside ::: the ::::::: station :::::::::::::::::: [START_REF] Barna | Remediation of old south pole stationphase i: ground-penetrating-radar surveys[END_REF]. Although the initial position of the snow cave, which has been fixed to match the conditions of the current field test on-going at Dome C, is a few meters shallower than that of the container, it turns out that both solutions show similar sinking rates, of the order of ∼ 5 cm a -1 , independently of the initial density in their immediate vicinity. Such slow sinking rates are valuable from a logistical point of view as they imply that the cores should remain easily accessible in the future: accessing the storage facility after several decades should not require much more work than that needed for its initial implementation. However, one must keep in mind that this result relies strongly on the assumption that the annual surface mass balance will not evolve significantly over the targeted lifespan of the storage facility.

Provided that a sufficiently large trench of high initial density is dig :::::::: excavated to welcome the balloon during the construction process, the rate of closure of a snow cave is relatively slow. Without in situ measurements, it is difficult to say whether such slow closure rates are related to the uncommon :::: fully ::::::::::: attributable ::

to ::: the ::::::::: favorable : conditions prevailing at Dome C(very low accumulation, very cold snow, ...), or if they are reflecting model biases, such as an overestimated initial density profile leading to underestimated flow velocities. First measurements of the test cave currently monitored at Dome C should be available soon and will shed light on this issue. In any case, despite these slow deformation rates, the cave should not be considered as perennial for timescales beyond ∼ 100 a. As a consequence, if the snow cave is the adopted solution for the Ice Memory project and if targeted lifetime goes beyond the century, a strategy regarding the long-term maintenance of the storage solution should be planned.

The choice of the time interval at which maintenance of the cave should be performed is a trade-off between the amount of work required to reach the cave (which increases as the height of snow above the cave increases), and the cost and logistic required to perform this work, which tends to prohibit a too short return time.

At the other end of the rigidity spectrum, we have considered the case of a perfectly rigid container. In fact, traditional shipping containers such as those used for logistics in Antarctica are designed in such a way that it is possible to stack up several of them on top of each other, the loads being fully supported by the four corner posts. The floor is composed of a reinforced base structure usually made of two bottom side rails and a number of crossmembers which are welded together as a sub-assembly on which plywood boards are longitudinally laid and fixed to support and transfer the load of the freight. In contrast, the roof and sides are usually constituted of several die-stamp corrugated steel sheets that are butt jointed together to form large panels by automatic welding.

In the traditional use of freight containers, these steel panels are not supposed to bear any significant load and their thicknesses are usually of a few millimetres only. As a consequence, given its overall dimensions, an unreinforced shipping container is obviously not able to bear the normal stresses that have been highlighted in the present study.

An option is then to strengthen the shipping containers operated at Dome C through ad hoc structures, which could either be set up inside the containers or mounted outside. Such structures must be designed in order to tackle both the mechanical and functional issues, i.e. to take over the snow loads while preserving sufficient space for core storage. This state of fact underlines the importance of the present study: a good knowledge of the magnitude, spatial repartition and temporal evolution of the normal stresses supported by the container is a prerequisite for the design of an optimal reinforcement structure.

Conclusion

One of the major challenges raised in the frame of the Ice Memory project is the design of a permanent storage facility, which would ensure safe subsurface storage over coming decades to centuries at affordable cost. To tackle this challenge, we have used a numerical model describing the viscous flow of firn to assess the temporal evolution of the two opposite end-member cases in terms of rigidity of the structure, i.e. an unreinforced snow cave and a perfectly rigid by :::: the ::::::::: container. :

The particular conditions prevailing at Dome C leads to slow rate of sinking of the storage facility and tend to increase its lifetime expectancy in comparison to what has been reported in available literature regarding subsurface constructions tried in the past. In particular, provided that the trench of high initial density surrounding the inflatable balloon involved in its construction process is sufficiently large, the closure rate of a snow cave dug :::::::: burrowed into the firn turns out to be relatively slow. Yet, if targeted lifespan of the storage solution exceeds the century, the use of rigid structure buried within the firn will have to be considered.

As the annual mass balance at Dome C is expected to remain positive in the long term, such a construction will be submitted to ever-increasing snow loads.

For a structure buried at an initial depth of 10 m, these loads are significant from the beginning of the simulations. In particular, strong stress concentrations occur at the container angles. Because their roof and sides are not designed to provide any significant mechanical resistance, traditional shipping containers such as those operated in Antarctica are not suited for this purpose, unless specifically reinforced. A possibility is then to take advantage of the results presented here in order to design optimal ad hoc reinforcement structures.

stitut Paul-Emile Victor (IPEV) and the italian Agenzia Nazionale per le Nuove Tecnologie, l'Energia e lo Sviluppo Economico Sostenibile (ENEA). All the computations presented in this paper were performed using the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is supported by Grenoble research communities : . :::: We :::::: thank ::::::::::: ENEA-UTA ::: for :::: the :::::::::::: coordination :: of :::: the ::::: snow Supplement 2: Sensitivity to the domain width

Here, we assess the sensitivity of the results to the total width of the domain.

In results presented in the main paper, the total width of the domain is of 22.44 m, which corresponds to a domain that is 10 m wider than the obstacle (i.e. the container or the cave) on both sides. Given the experimental set up, for which none of the model parameters have any dependency on the horizontal coordinate x, the flow of firn should not be perturbed by the presence of the obstacle beyond a certain distance to the latter. The question is whether 10 m is a sufficient distance so that the flow regime at the right and left boundaries of the domain is not affected by the presence of the obstacle. Otherwise, it can be expected that the results obtained in the immediate vicinity of the obstacle will be affected by the total width of the domain. This sensitivity analysis was performed for the container case only. Indeed, the container being a rigid obstacle that the firn has to bypass in its flow, we expect that the latter will be perturbed over a larger distance than in the case of the cave, which is a free surface. In other words, should the modelled domain be sufficiently wide for the case of the container, it will also be for the case of the cave.

To confirm that a horizontal dimension of 22.44 m is sufficient, we run an additional simulation, which corresponds to the reference simulation of the container case without weight as presented in Section 2.2.3, except that the domain total width is extended to be comprised between x l = 0 m and x r = 102.44 m.

The container is then placed in the middle of the domain, i.e. between x = 50 m and x = 52.44 m, so that the domain is 50 m wider than the container on both sides. Figure S4 shows a comparison of the normal stresses obtained along the container roof, the container floor and the container sides for the two considered domain width. As shown by this figure, results turn out to be almost insensitive to the domain width. This means that a total width of 22.44 m is sufficient so that the flow regime at the left and right boundaries of the domain is not affected by the presence of the container, and it is then useless to make the domain wider. Indeed, the computation time is close to 3 times larger for the shape and the position relative to the surface of the cave over the course of the simulation for each of the considered initial density fields.

@0a @50a @100a For the two cases for which the trench is larger than the balloon, the cave is sort of pinched and similar pattern of deformations occur on the roof and on the floor of the cave, so that the points corresponding to the right and left extremities of the firn/cave interface almost define a horizontal symmetry axis. In contrast, the initial horizontal symmetry of the cave rapidly vanishes for the case of the narrow trench and even more so for the reference run. For each one of the four considered cases, the deformation of the upper half of the circle is very similar to the deformation of the roof of the cave when the latter is the combination of a square and half circle over the first 40 to 50 years of simulation. In the same time, the lower half of the circular cave deforms much slower than that of the half circle/square combination cave, for which the sides of the square tend to curve inward. These differences in the deformations of the lower halves of the two shapes end up affecting the upper halves, and after the first 40 to 50 years of simulation the differences in terms of deformation between the upper half of the circular cave and the upper half of the circle/square combination cave become perceptible, with the upper half of the circle/square combination cave shrinking faster than that of the circular cave. At the end of the simulation, although the initial volume of the cave was higher for the circle/square combination than for the circular cave, the relative volume loss of the former is higher than that of the latter.

Despite the slight differences in terms of deformation patterns described above, the initial shape of the cave turns out to have minor influence on the lifetime expectancy of the cave, at least for the two tested shapes. In particular, this initial shape is much less critical than the initial shape of the trench in which the balloon is placed during the construction phase. Therefore, we suggest to opt for the circle/square combination as the initial shape of the cave, which is obviously more convenient for storage due to its flat floor, and to concentrate the efforts on having a large trench of high initial density surrounding it. Ideally this trench should have a rectangular shape or, if not possible, a T-shape with a sufficiently large top branch.

Supplement 4: Sensitivity to the firn/container boundary condition

Here, we assess the sensitivity of the results obtained for the rigid container case to the implemented boundary condition at the firn/container interface.

As stated in Section 2.2.3, all simulations regarding the rigid container case presented in the main text rely on a free slip condition at the firn/container interface. However, it is also possible that, instead of flowing freely around the container, the snow sticks to its walls. In terms of modelling, this latter case corresponds to a no slip boundary condition. Therefore, we run an additional simulation, which corresponds to the reference simulation without weight as the firn/container interface turns out to be low, with only 4% (resp. 5%) of difference for the normal stress on the middle part of the container roof (resp.

container floor) at the end of the simulation. This difference is slightly higher on the container sides, with a difference of 7% on the maximum normal stress exerted on the lower part of the container sides at the end of the simulation. In addition, the use of a free slip boundary condition at the firn/container interface put us on the safe side since it produces the highest normal stresses, which are the ones that need to be considered when dimensioning a reinforcement structure.

  7) Thus, in the limiting case for which D = 1 (pure ice), a = 1 and b = 0, and the snow/firn law reduces to the classical incompressible Glen's flow law traditionally used to model glacier flow. For smaller relative densities, we use parameterized forms of functions a and b that were first proposed by Gagliardini and Meyssonnier (1997) by fitting cold room experiments and densification measured at Site 2 (Greenland), before being slightly corrected by Zwinger et al. (2007) for a study of the Gorshkov crater glacier (Kamchatka, Russia) as: .09371 -20.46489D), D ≤ 0.81 b 0 (D), 0.81 < D ≤ 1.0 (9)

Figure 1 :

 1 Figure 1: ::::::: Modelled : domain and associated boundary conditions implemented for (a) the initialisation, (b) the snow cave experiment and (c) the rigid container experiment. Top and lateral boundary conditions represented on (a) are also valid for (b) and (c) but were not reported for the sake of readability. The hatched areas in (b) and (c) represent the various considered shapes for the trench of high initial density surrounding the cave/container as mentioned in the text: the narrow trench (green), the wide trench (blue), and the T-shape trench (red) for (b) only.

Figure 2 :

 2 Figure 2: Evolution of (a) the vertical relative density profile and of (b) the vertical velocity profile over the course of the initialisation run. In (a), red crosses correspond to measurements of densities on the 1999 FireTrack firn core (Augustin et al., 2004), green triangles are density measurements on the 2010-2011 Volsol core (Burr et al., 2018), light and dark red points correspond to measurements carried out on the two 80-m-long 2012/2013 ice cores by Leduc-Leballeur et al. (2015) and used to derive the initial density profile (black).

Figure 3 :

 3 Figure 3: Shape and position of the cave at t = 0 a, t = 50 a and t = 100 a for the various prescribed initial density field: reference simulation (black), narrow trench case (green), large trench case (blue) and T-shape trench case (red). Note that the four considered cases have been positioned arbitrarily on the horizontal axis and that, for the sake of readability, the slight deformations of the firn surface occurring over time are not represented.

Figure 4 :

 4 Figure 4: Evolution over time of the position of (a) the pseudo-center of the cave and (b) the roof of the container for the various considered cases: reference simulation (black), narrow trench case (green), large trench case (blue), and T-shape trench case (red) for (a) only.The position over time of a snow particle in the unperturbed initial steady density field is represented in orange (superimposed with or very close to blue line in both cases). In (b), results obtained without considering the weight of the container (dashed lines) and accounting for the weight of the container and its payload (continuous lines) are both represented.

Figure 5 :

 5 Figure 5: Evolution over time of the snow relative density along container roof for the six considered cases. For the two reference simulations (i.e. with and without weight), we report the value of the unperturbed initial relative density at the depth at which the container roof finds at the considered time of the simulation (black dashed lines). For the no weight (resp.

Figure 6 :

 6 Figure 6: Evolution over time of the normal stresses along container roof for the six considered cases. For the no weight (resp. with weight) large and narrow trenches simulations, we report the results obtained for the no weight (resp. with weight) reference simulation to ease comparison (black dotted lines).

:

  :::::::: considered :::::: depth ::::: from ::::::::: available ::::::: density ::::::: profiles, :::: we ::::: make :: a ::::::: further :::::::: analysis ::::: which :::::::: consists :: in ::::::::::: quantifying ::: the :::: gap :::::::: between ::::: these :::: two :::::::::: quantities. ::::::: Figure :: 7 :::::

Figure 7 :

 7 Figure 7: Evolution over time of the ratio between calculated ::::::: modelled normal stresses and hydrostatic pressure along container roof for the six considered cases.

Figure 8 :

 8 Figure 8: Evolution over time of the normal stresses along container side for the six considered cases. The vertical axis is the normalized distance along container side, with z = 0 m corresponding to the roof position and z = -2.9 m corresponding to the floor position. For the no weight (resp. with weight) large and narrow trenches simulations, we report the results obtained for the no weight (resp. with weight) reference simulation to ease comparison (black dotted lines).

  the :::: Old ::::: South :::: Pole ::::::: Station :: is :::::::: probably :::: the ::: one ::: for :::::: which ::::::::::: surrounding ::::::::: conditions ::: are ::: the :::::: closest ::: to ::::: those ::::::::: prevailing :: at :::::: Dome :: C ::::: (very :::: low ::::::::::::: accumulation, :::: very :::: cold ::::: snow, ::::::: almost :::::: purely :::::::: vertical :::: firn ::::::: motion, :::: ...). :::::: The :::: Old :::::: South ::::: Pole ::::::: Station ::: was ::::: built :::::: during :::: the ::::::: 1956-57 ::::::: austral ::::::: summer ::::: and, :::::: despite ::::::::: relatively :::: low :::::: annual :::::::::::: accumulation ::::::::::::: (∼ 20 cm a -1 , ::::::: mostly :::: due ::: to ::::: wind ::::::::::: transport), :::: the ::::: roofs ::::: were :::::::

  :: the ::::::: station :::: kept :::::: being :::::::: operated ::::: until ::: the :::::: 1970s. ::: In ::::: 1997, ::: i.e. ::: 40 ::::: years ::::

  :::::: pictures ::::: show :::::: some :::::::: localised ::::::::: distresses ::: of :::: the ::::::::: structure, ::::: such ::: as :::: the :::::: failure :: of ::::: some ::::::: timbers :: in :: a ::::::: corridor :::: and ::: the :::::::: crushing ::: of : a :::::::::: corrugated ::::: steel :::: arch ::::::(which ::::: must :::: have ::::: been ::::::::: enhanced :: by ::::::::::::: asymmetrical ::::: loads :::::: related ::: to :::::: drifted ::::: snow :::::::::::: accumulation :: in ::: the ::::: main ::::: wind :::::::: direction ::: as ::: the ::::::::::: construction :::: was ::::::: initially :::::: raised ::on :::: the :::::::: surface). :::::::::: However, ::::: these :::::::: failures ::::: seem ::::: quite ::::::: limited :::::: when ::::::::: compared ::to :: a ::::::: similar :::::: photo :::::: report :::::::::: performed :: at ::::::: Camp :::::::: Century :::::::::: (northwest :::::: corner ::: of ::: the :::::::::: Greenland ::: Ice :::::: Sheet) ::: in ::::: 1969, :::: only :::: ten ::::: years ::::: after :::: the ::::::::::: construction :: of :::: the :::::: station ::::::::::::::[START_REF] Kovacs | Camp century revisited, a pictorial view -june 1969[END_REF]. ::::: This ::::: tends ::: to ::::::: support ::::: that ::: (1) :::::::: locations ::::: with :::: very :::: low :::::::::::: accumulation ::::: rates :::: and :::::: almost :::::: purely ::::::: vertical ::::: flow, ::::: such :: as :::::: Dome ::: C, ::: are ::::: ideal ::for ::::::::::: undertaking ::: of ::::::::: long-term ::::::::: subsurface ::::::::::::: constructions ::: and :::: (2) : it :: is :::::::::: preferable :: to :::: have ::: the ::::::: storage ::::::: facility :::::: buried :::::: below :::: the ::::::: surface :::: from :::: the ::::::::: beginning ::: to ::::: avoid :::::::: problems ::::::: related ::: to ::::: snow ::::: drift, ::::: such ::: as :::: fast :::::::::::: accumulation :::: and ::::::::::::: asymmetrical ::::: loads. :

  container. :::: The ::: use :: of ::::: such :: a ::::: model ::::::: proved ::::::::: necessary ::: not :::: only ::: to ::::: assess :::: the :::: cave ::::::::::: deformation :::: over :::: time :::: but :::: also :: to :::: get ::::::: realistic ::::::::: estimates :: of :::: the ::::: loads ::::::::: supported ::

Figure S1 : 1 TimeFigure S2 :Figure S3 :

 S11S2S3 Figure S1: Evolution over time of the snow relative density along container floor for the six considered cases. For the two reference simulations (i.e. with and without weight), we report the value of the unperturbed initial relative density at the depth at which the container roof finds at the considered time of the simulation (black dashed lines). For the no weight (resp.with weight) large and narrow trenches simulations, we report the results obtained for the no weight (resp. with weight) reference simulation to ease comparison (black dotted lines).

Figure S5 :

 S5 FigureS5: Shape and position of an initially perfectly circular cave at t = 0 a, t = 50 a and t = 100 a for the various prescribed initial density field: reference simulation (black), narrow trench case (green), large trench case (blue) and T-shape trench case (red). The results regarding the half circle/square combination cave presented in the main paper are also reported for ease of comparison. Note that the four considered cases have been positioned arbitrarily on the horizontal axis and that, for the sake of readability, the slight deformations of the firn surface occurring over time are not represented.

Table 1 :

 1 List of parameter values used in this study.

	of the experiments
	2.2.1. Initialisation
	Using the model described above, we run a first transient experiment in
	order to produce an initial steady density field. The domain considered for this

Acknowledgements

This work was funded by Fondation Université Grenoble Alpes. This study was performed as a part of the Ice Memory project, which is led by the french In-

Code and data availability

Elmer/Ice code is publicly available through GitHub (https://github.com/ElmerCSC/elmerfem).

All the simulations were performed with the version 8.4 (Rev: 6b6e592f) of Elmer/Ice. All scripts used for simulations and post-treatment as well as model output are available upon request from authors.

Author contribution

Olivier Gagliardini and Julien Brondex designed the study. Julien Brondex conducted the simulations with support of all co-authors. Julien Brondex wrote the paper with contributions of all co-authors.

Declaration of Competing Interest

All authors declare that they do not have any conflict of interest Supplement 3: Sensitivity to the initial shape of the cave Here, we want to assess the sensitivity of the results obtained regarding the snow cave to its initial shape. To this end, we performed similar simulations as the ones presented in section 2.2.2, the only difference being that the initial shape of the cave is a perfect circle. Figure S5 shows the evolution of the presented in Section 2.2.3, but with a no slip condition instead of the free slip condition at the firn/container interface. Figure S6 shows a comparison of the normal stresses obtained along the container roof, the container floor and the container sides for the two considered boundary conditions.

The sensitivity of these stresses to the implemented boundary condition at


Maud Alligier

Anne-Laure Borel

Véronique Savey

Claire Rives-Lange

Marie-Claude Brindisi

Xavier Piguel

David Nocca

Maud Monsaingeon-Henry

Emilie Montastier

Sophia Beliard

Borel Al

Brindisi Mc

Bossu Estour

Lemoine Ay

Cécile Bétry
email: cbetry@chu-grenoble.fr

A series of severe neurologic complications after bariatric surgery in France: the NEUROBAR Study

Keywords: Neurological complications after bariatric surgery Bariatric surgery, Nutritional Status, Wernicke Encephalopathy, Thiamine Deficiency, Vitamin B Deficiency

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Obesity is a major health problem with increasing prevalence [START_REF] Gbd ; Afshin | Health Effects of Overweight and Obesity in 195 Countries over 25 Years[END_REF]. Inducing significant weight loss with medical management remains difficult. In this context, bariatric surgery (BS) is recognized as being the most efficient treatment of severe obesity [START_REF] Angrisani | IFSO Worldwide Survey 2016: Primary, Endoluminal, and Revisional Procedures[END_REF][START_REF] Colquitt | Surgery for weight loss in adults[END_REF]. Most studies in this field have focused on weight loss and resolution of comorbidities. In contrast, the specific area of long-term nutritional consequences of BS has been overlooked [START_REF] Montastier | Long-term nutritional follow-up post bariatric surgery[END_REF].

In 2016, more than 50 000 patients underwent a BS in both public and private centers in France with a gradual increase of BS performed every year since 1997 [START_REF] Galvin | EFNS guidelines for diagnosis, therapy and prevention of Wernicke encephalopathy[END_REF]. The increasing number of bariatric procedures in France raises the question of the nutritional long-term followup [START_REF] Oudman | Preventing Wernicke Encephalopathy After Bariatric Surgery[END_REF]. In this context, French Centers of Obesity Care Management were labeled by the French Health Ministry. They are referral centers for obesity management, as well as for BS and their complications.

Neurological diseases are one of the most severe complications of BS related to nutritional impairment. They can have dramatic consequences, including death and permanent disabilities, as we previously reported [START_REF] Czernichow | Current challenges in providing bariatric surgery in France: A nationwide study[END_REF]. Several case reports, but only few case series, with such complications have been published yet [START_REF] Thereaux | Longterm adverse events after sleeve gastrectomy or gastric bypass: a 7-year nationwide, observational, population-based, cohort study[END_REF][START_REF] Bétry | Need for Intensive Nutrition Care After Bariatric Surgery[END_REF][START_REF] Algahtani | Neurological complications of bariatric surgery[END_REF]. The most documented neurological complication after BS is Wernicke encephalopathy. This disease is due to thiamine deficiency and was initially described in subjects with alcoholism [START_REF] Fragoso | Neurological complications following bariatric surgery[END_REF][START_REF] Juhasz-Pocsine | Neurologic complications of gastric bypass surgery for morbid obesity[END_REF]. The factors associated with such complications remain to be determined.

Materials and Methods

Study description

This was a French national multicentric observational retrospective study. To limit the heterogeneity of practices, only patients who underwent BS after January 2010 and up to November 2018 have been included. Every patient with a neurological complication including encephalopathy or peripheral neuropathy and with a history of BS could be included.

Data collection

Each medical coordinator of the 37 French Centers of Obesity Management received an email explaining the purpose of the NEUROBAR study. Public or private hospitals regularly refer patients with complications after BS to these referral centers. Medical coordinators were asking, if they were aware of patients with neurological complications following BS in their territory. Patients may have been managed inside or outside of their centers. Inclusion of cases was left to the discretion of the medical coordinators. Then, they or their collaborators were provided with a link to an online questionnaire, and, were asked to fulfill one questionnaire for each case. This online questionnaire was defined by 50 items related to 7 categories, namely characteristics, of the case, of the BS, and of the neurological complications, differential diagnosis, factors contributing to the neurological complications, case management, current status of the case (Supplemental data 1). We collected data about diabetes, which is a frequent cause of peripheral neuropathy. We also collected data about psychiatric history including psychotropic drug intake to determine the active psychiatric comorbidities. Indeed, psychiatric diseases could be a factor for poor vitamin adherence after bariatric surgery [START_REF] Sunil | Predictors of Vitamin Adherence After Bariatric Surgery[END_REF]. We used an online survey management system, hosted on a secured and authorized system for personal health data collect (Grenoble Hospital). The questionnaire was closed at the end of November 2018. The FORCE network has been strongly involved within the NEUROBAR study.

Ethics and data protection:

The Hospices Civils de Lyon are the sponsor of this retrospective study. None personally identifying information (such as name, day of birth, or social security number...) were collected. Data were collected during routine medical examinations. Patients did not undergo any treatment or examination specifically devised to collect data for this study. Thus, in accordance with the French regulation, the study was approved by the CCTIRS (Comité Consultatif sur le Traitement de l'Information en matière de Recherche dans le domaine de la Santé), a French committee dedicated to data protection for biomedical research. Patients were also individually informed that their data would be used for research (MR-003 reference methodology). This study has been registered by the data protection officer (DPO) of the Hospices Civils de Lyon.

Data analysis:

Data analysis was performed using RStudio software version 1.0.136 with occasional reliance on 'prettyR', 'RVAideMemoire' [START_REF]RStudio: Integrated Development Environment for R[END_REF]. Owing to the small number of patients, data are presented as median and interquartile range (IQR).

Results

Twenty centers participated in this study and seventeen did not. Three centers had no patient corresponding to the targeted phenotype. Forty-six patients were eligible. Eight patients from six centers declined to take part in the study. Overall, 38 patients from 13 centers were included (Supplemental Table 1).

Subjects characteristics

There were 34 females and 4 males. At the time of the surgery, the median age was 39 [IQR: 26-51years old and the median body mass index (BMI) was 43 [IQR: 39-47] kg.m 2 . Eight subjects were on psychotropic drugs at the time of surgery and five subjects had diabetes history.

Surgery characteristics

Sleeve gastrectomy and gastric bypass were the two main surgeries (Figure 1A). Among the gastric bypass goup, one-third was anastomosis gastric bypass (OAGB) and two-third were Roux-en-Y gastric bypass (RYGB). Patients were operated in university hospital, private local hospitals or public local hospitals with a fairly similar proportion (Figure 1B). Patients from private or public local hospitals were secondarily referred to French Centers of Obesity Management given the neurological complications. Five patients had a history of two or more previous bariatric surgeries before the neurological event. Only 22 cases (58%) received vitamin and trace element supplementation with uncertain treatment adherence in ten patients before the neurological complication.

Surgery and gastrointestinal complications

Twenty cases (53%) had both neurological and surgical complications, including six cases of fistula, four cases of anastomotic ulcer of gastric bypass, three cases of stenosis of sleeve gastrectomy and two cases of gastric twist after sleeve gastrectomy.

At least 26 cases had gastrointestinal symptoms (68%). Among them, 77 percent of cases had vomiting, which was the most frequent symptom (Figure 2). In 15 out of the 26 cases (58%), gastrointestinal symptoms were associated with a surgical complication. For 7 out of the 12 patients without gastrointestinal symptoms, limited oral energy intake was reported in the online questionnaire.

Neurological complications after BS

The median duration between the surgical procedure and the onset of the neurological complication event was 6 [IQR: 3-12] months. Ten cases had an encephalopathy, fifteen had a peripheral neuropathy, twelve had both and one case had a pyramidal syndrome.

Every case of encephalopathy had at least one of the three signs of Gayet-Wernicke encephalopathy namely eye signs, cerebellar signs, and either mild memory impairment or an altered mental state in a context of dietary deficiencies. Eighteen out of twenty-two had an imaging examination (Magnetic Resonance Imaging, MRI or computed tomography scan, CTscan). In the majority of cases, the imaging examination was considered normal. Vitamin B1 deficiency was proven in 14 out of 38 cases (Table 1).

Nutritional status of subjects with neurological complications after BS At the time of neurological care, nutritional status assessed with BMI, weight loss or albumin level was extremely variable (Figure 3). Twenty-eight cases (74%) had at least one B vitamin deficiency: vitamin B1 deficiency was the most frequent. Deficiencies in trace elements were not rare including deficiency in selenium in 14 cases (Figure 4).

Nutritional care of neurological complications and prognosis

Twenty patients required artificial nutrition, enteral and/or parenteral nutrition. In the majority of cases, patients were supplemented with B vitamins (vitamin B1: 76%, vitamin B6: 68%, vitamin B12 42% vitamin B3 32%). Twelve patients (32 %) underwent a reoperation including three cases of conversion of OAGB to RYGB and two cases of reversal gastric bypass to normal anatomy.

The neurological symptoms were completely resolved in only nine cases. Neurological sequelae included cognitive impairments and symptoms associated with peripheral neuropathy (neuropathic pains, hypoesthesia, cognitive disorders, etc.). Two cases with encephalopathy died with one case directly imputable to the neurological complication of BS.

Discussion

With the increasing development of BS worldwide, nutritionists are facing new challenges related to medical complications of BS [START_REF] Angrisani | IFSO Worldwide Survey 2016: Primary, Endoluminal, and Revisional Procedures[END_REF]. In a recent French series of patients requiring intensive nutrition care after BS, half of the 12 cases displayed neurological disorders, namely encephalopathy or neuropathy, following BS with possible dramatic consequences [START_REF] Czernichow | Current challenges in providing bariatric surgery in France: A nationwide study[END_REF].

Despite the importance of this topic, few studies reported more than 10 cases of neurological complications after BS [START_REF] Thereaux | Longterm adverse events after sleeve gastrectomy or gastric bypass: a 7-year nationwide, observational, population-based, cohort study[END_REF][START_REF] Bétry | Need for Intensive Nutrition Care After Bariatric Surgery[END_REF][START_REF] Algahtani | Neurological complications of bariatric surgery[END_REF]. With the multicentric NEUROBAR study, we aimed to report a French series of patients with severe neurological complications after BS.

The most striking result emerging from the data is that gastrointestinal symptoms (68%) and especially frequent vomiting was a common factor associated with neurological complications in our case series. In contrast, only 7.4 % of patients had often or daily vomiting after BS in a previous study [START_REF] Cano-Valderrama | Long-Term Food Tolerance After Bariatric Surgery: Comparison of Three Different Surgical Techniques[END_REF]. Our results are consistent with the Oudman et al.'s review showing that vomiting is the most common symptom among patients with Wernicke encephalopathy after BS [START_REF] Juhasz-Pocsine | Neurologic complications of gastric bypass surgery for morbid obesity[END_REF]. Surgical complications are also present in more than half of the patients, which is higher than what is reported [START_REF] Gloy | Bariatric surgery versus non-surgical treatment for obesity: a systematic review and meta-analysis of randomised controlled trials[END_REF] and may be involved in decreased food intake. In the NEUROBAR study, 32% of patients underwent a reoperation versus 8% in the systematic review of Gloy et al. [START_REF] Gloy | Bariatric surgery versus non-surgical treatment for obesity: a systematic review and meta-analysis of randomised controlled trials[END_REF]. In contrast, some characteristics (sex, age, type of surgery or comorbidities, etc.) are relatively similar in our study compared to a French nationwide study in 2013 [START_REF] Galvin | EFNS guidelines for diagnosis, therapy and prevention of Wernicke encephalopathy[END_REF]. Moreover, the delay between BS and the neurological complication is variable, as well as the nutritional status. Some patients can develop early neurological complications without large weight loss or decreased albumin level. It suggests that protein-energy malnutrition is not necessarily associated with neurological complications after BS.

Another remarkable result of our series is the high rate of B vitamin deficiencies (74%) and especially vitamin B1 deficiency. Most cases of encephalopathy, described here, displayed symptoms compatible with Wernicke encephalopathy according to the definition in the European Federation of Neurological Societies (EFNS) guidelines [START_REF] Fragoso | Neurological complications following bariatric surgery[END_REF]. Despite this, brain imaging is often normal. It highlights the difficulty to diagnose Wernicke encephalopathy with atypical presentation in non-alcoholic patients [START_REF] Fragoso | Neurological complications following bariatric surgery[END_REF]. Moreover, in our study, no deficiency in thiamine was reported in some cases, despite compatible clinical features. It can be explained by a supply in thiamine before any blood analysis or no vitamin B1 blood test, which is not always indicated in the medical records.

We are aware that our research may have some limitations. There are several sources of selection bias: 1) We collected cases only from the 37 French Centers of Obesity Management.

2) Inclusions were left to the discretion of the medical coordinator of every center 3) Sixteen centers did not reply, and some of the patients declined to participate in the study. Thus, our findings might not be transferable to the French population. On a wider level, a study based on Health Insurance Fund for Salaried Workers, which covers a large majority of the French population, would be useful to estimate the prevalence of neurological complication after BS in France. Secondly, it was not possible to determinate odds ratio for risk factors of neurological complications after BS, given the design of our study. We compared our results with previous published results to suggest that some factors can be risk factors. Neurological complications after BS remain rare, estimated below 5% [START_REF] Algahtani | Neurological complications of bariatric surgery[END_REF][START_REF] Tabbara | Rare Neurological Complications After Sleeve Gastrectomy[END_REF], and thus, it is difficult to design prospective studies. Further work based on matched-case control studies could confirm the overrepresentation of surgical complications and gastrointestinal symptoms in patients with neurological complications after BS. Thirdly, we are aware that our research may have limitations related to its retrospective design. To limit the potential bias and to ensure homogeneity, we used a standardized electronic questionnaire with a limited number of questions and we collected data only from 2010. Eventually, we cannot state that all neurological pathologies are directly related to BS. To limit this bias, we included some items in the questionnaire to rule out differential diagnoses, e.g. infections.

We hope that our research will help prevent neurological complications after BS.

Regarding our case series, we suggest that every patient with gastrointestinal symptoms, in particular, frequent vomiting, should benefit from a prompt consultation with the surgeon and the nutritionist. It could help identify surgical complications to and may limit neurological complications. Indeed, we found that gastrointestinal symptoms and surgical complications were associated in more than 50% of the cases of severe neurological complications of BS.

Moreover, most of the neurological complications were associated with deficiencies in B vitamin and especially in thiamine. We suggest to follow the EFNS guidelines for prevention of Wernicke encephalopathy, namely, parenteral thiamine supplementation for every patient at risk of such complications [START_REF] Fragoso | Neurological complications following bariatric surgery[END_REF]. Furthermore, deficiencies in trace element were not rare. Our study supports the importance of micronutrient deficiency screening and systematic supplementation.

Conclusions 260

In conclusion, gastrointestinal symptoms and surgical complications were common in our series. Being aware that such circumstances could be risk factors of neurological complications after BS might help to early set up prevention strategies. B vitamins deficiencies, especially vitamin B1 deficiency, are also frequent in our study. 265

65The

  NEUROlogical complications following BARiatric surgery (NEUROBAR) study was conducted in France. It aimed to widen current knowledge on severe neurological diseases after BS. It was coordinated by the French Obesity Research Centre of Excellence (FORCE) network, which is a clinical investigation network, labelled by French Clinical Research Infrastructure Network (F-CRIN) in 2014. This network has developed a well-structured 70 national coordination for harmonization and standardization of clinical investigations and data collection.

Figure 1 .Figure 3 .Figure 4 :

 134 Figure 1. Surgery characteristics of cases with neurological complications after bariatric

Figure 2 .

 2 Figure 2.: Gastrointestinal symptoms among subjects with neurological complications after

Figure 3 . 20 Figure 4 :

 3204 Figure 3. Nutritional status ((A) BMI, (B) percentage of body weight loss defined as total

Table 1 .

 1 Characteristics of patients with encephalopathy 345 MN: motor neuropathy; SN: sensory neuropathy; UE: upper extremities; LE: lower extremities; GWE: Gayet Wernicke encephalopathy N/A: data not available Figure 1. Surgery characteristics of cases with neurological complications after bariatric surgery 1A. Type of bariatric surgery (data available in 34 cases) 1B. Type of hospital for the place of surgery (no available data for 4 cases) 5 OAGB, one anastomosis gastric bypass