Kévin Fourteau
email: kfourteau@protonmail.com

Pascal Hagenmuller

Jacques Roulle

Florent Domine

On the use of heated needle probes for measuring snow thermal conductivity

Keywords: Snow, snow physics, glaciological instruments and methods

Heated needle probes provide the most convenient method to measure snow thermal conductivity. Recent studies have suggested that this method underestimates snow thermal conductivity; however the reasons for this discrepancy have not been elucidated. We show that it originates from the fact that, while the theory behind the method assumes that the measurements reach a logarithmic regime, this regime is not reached within the standard measurement procedure. Using the needle probe without this logarithmic regime leads to thermal conductivity underestimations of tens of percents. Moreover, we show that the poor thermal contact between the probe and the snow due to insertion damages results in a further underestimation. Thus, we encourage the use of fixed needle probes, set up before the snow season and buried under snowfalls, rather than hand-inserted probes. Finally, we propose a method to correct the measurements performed with such fixed needle probes buried in snow. This correction is based on a lookup table, derived specifically for the Hukseflux TP02 needle probe model, frequently used in snow studies. Comparison between corrected measurements and independent estimations of snow thermal conductivity obtained with numerical simulations shows an overall improvement of the needle probe values after application of the correction.

Introduction

Understanding heat transfer through a snowpack is essential to model the heat budget of the Earth surface as well as snow metamorphism [START_REF] Colbeck | Theory of metamorphism of dry snow[END_REF]. Quantifying the snow thermal conductivity, which links the heat flux to the temperature gradient through Fourier's law is thus of primary importance, as it plays an integral role in the energy budget and thermal state of snowpacks, as well as in the energy budget of the ground (e.g. Zhang and others, 1996), ice caps and glaciers (e.g. Gilbert and others, 2012), or even sea ice (e.g. [START_REF] Lecomte | On the formulation of snow thermal conductivity in large-scale sea ice models[END_REF]. Consequently, the thermal conductivity of snow has been investigated in detail for more than half a century (Yosida and others, 1955;[START_REF] Jaafar | Thermal conductivity of snow by a transient state probe method[END_REF][START_REF] Sturm | Thermal conductivity measurements of depth hoar[END_REF]Morin and others, 2010;[START_REF] Calonne | Numerical and experimental investigations of the effective thermal conductivity of snow[END_REF][START_REF] Riche | Thermal conductivity of snow measured by three independent methods and anisotropy considerations[END_REF]Domine and others, 2015).

For the last 50 years, snow thermal conductivity has generally been measured via the heated needle probe method, and this technique still remains the only one easily applicable in the field [START_REF] Jaafar | Thermal conductivity of snow by a transient state probe method[END_REF][START_REF] Sturm | Thermal conductivity measurements of depth hoar[END_REF]Domine and others, 2015). This method is based on the transient heating of a wire embedded in a snow sample and the monitoring of the subsequent increase of the wire temperature. The thermal conductivity of the tested medium can be directly related to this temperature increase. For instance, a highly conductive medium will dissipate heat away from the probe and the wire temperature will increase slowly. Thanks to its simplicity, the needle probe method can be used for both measurements performed in the field (e.g. Morin and others, 2010) and in the laboratory (e.g. [START_REF] Sturm | Thermal conductivity measurements of depth hoar[END_REF][START_REF] Calonne | Numerical and experimental investigations of the effective thermal conductivity of snow[END_REF][START_REF] Riche | Thermal conductivity of snow measured by three independent methods and anisotropy considerations[END_REF]. Moreover, needle probes can be left in the snowpack and their use automated, allowing the measurement of the temporal evolution of thermal conductivity over an entire snow season (Morin and others, 2010;[START_REF] Domine | Soil moisture, wind speed and depth hoar formation in the arctic snowpack[END_REF]others, 2015, 2018).

However, recent studies have highlighted that the needle probe method leads to systematic underestimations of the thermal conductivity by up to about 50% when compared to other techniques [START_REF] Calonne | Numerical and experimental investigations of the effective thermal conductivity of snow[END_REF][START_REF] Riche | Thermal conductivity of snow measured by three independent methods and anisotropy considerations[END_REF]. While different hypotheses have been proposed for the origin of this underestimation [START_REF] Riche | Thermal conductivity of snow measured by three independent methods and anisotropy considerations[END_REF], the actual reasons remain unclear. The goal of this paper is to explore these reasons, and to propose a correction for the use of the needle probe method in snow when possible.

Background

This section aims to provide the theoretical background underlying the needle probe method and to detail its specific application to snow.

Theory

The needle probe (hereafter NP) is an established method, used for media other than snow, such as soil (e.g. He and others, 2018) and food (e.g. [START_REF] Murakami | Recommended design parameters for thermal conductivity probes for nonfrozen food materials[END_REF]. As explained in the Introduction, it relies on the measurement of the internal temperature of a cylindrical probe, heated at constant power. In order to quantitatively relate the increase of temperature of the probe to the thermal conductivity of the investigated medium, the equation governing the temperature increase of an infinitely long and perfectly conducting cylindrical probe heated with a constant power and embedded in a homogeneous and isotropic medium have been investigated in detail, notably by [START_REF] Blackwell | A transient-flow method for determination of thermal constants of insulating materials in bulk part I-Theory[END_REF] or [START_REF] Jaeger | Conduction of heat in an infinite region bounded internally by a circular cylinder of a perfect conductor[END_REF]. The heat equation has been solved taking into account a potential thermal contact resistance between the probe and the surrounding medium. As the surrounding medium is considered infinite, no boundary condition needs to be prescribed at its boundary. For a constant heating rate Q (expressed in Watts per meter of probe length) the temperature elevation ΔT of the probe is given by Eqn (18) of [START_REF] Jaeger | Conduction of heat in an infinite region bounded internally by a circular cylinder of a perfect conductor[END_REF],

DT(t) = Q k 2a 2 p 3 1 0 (1 -e (-(t/t)u 2 ) ) u 3 D(u) du, ( 1
)
where t is the elapsed time since the start of the heating cycle expressed in s, k is the thermal conductivity of the isotropic homogeneous medium investigated expressed in W K -1 m -1 , α is a parameter defined as α = 2 c medium /c probe , with c medium and c probe the volumetric thermal capacities of the investigated medium and of the probe, respectively, and both expressed in J m -3 K -1 , τ is a characteristic time expressed in s and defined as τ = r 2 /χ, with r the radius of the probe in m and χ the thermal diffusivity of the investigated medium expressed in m 2 s -1 . Further Δ(u) is a function defined as

D(u) = [uJ 0 (u) -(a -hu 2 )J 1 (u)] 2 + [uY 0 (u) -(a -hu 2 )Y 1 (u)] 2 , ( 2
)
where J n and Y n are the Bessel functions of order n of the first and second kinds, and h = H k/r with H the thermal contact resistance between the probe and the external medium, expressed in m 2 K W -1 , and u is an integration variable without physical significance.

Equation (1) contains potential unknowns beside thermal conductivity, such as the thermal contact resistance between the probe and the medium, and is therefore seldom used in practice. For long times t, it can be simplified to its asymptotic development [START_REF] Jaeger | Conduction of heat in an infinite region bounded internally by a circular cylinder of a perfect conductor[END_REF],

DT = Q 4pk 2h + ln 4 C t t - (4h -a) 2a t t + a -2 2a t t ln 4 C t t , (3)
where C = e γ ≃ 1.7811 (γ being Euler's constant). For even larger times, the two high-order terms are negligible compared to the first logarithm term, and the asymptotic development can be further simplified to

DT = Q 4pk 2h + ln 4 C t t . (4)
In this case, the rate of increase of ΔT is related to the logarithm of time through the thermal conductivity k and the heating power Q.

One can therefore estimate the thermal conductivity k with

k = Q 4p dDT d ln (t) -1 . ( 5
)
This equation is the one used in practice to measure the thermal conductivity of a medium with the NP method (e.g. [START_REF] Sturm | Thermal conductivity measurements of depth hoar[END_REF], as it does not require to know the physical properties of the probe (its radius and volumetric thermal capacity), the potential thermal contact resistance at the probe/ snow interface, or the volumetric thermal capacity of the snow.

Here, only the heating power Q needs to be known without any additional calibration. To estimate the thermal conductivity of the medium, the measured temperature increase ΔT is first plotted against ln(t). Then, in the portion where the curve appears linear, the thermal conductivity k is deduced from the inverse of the slope and Eqn (5).

The needle probe method in snow sciences

Due to its ease of use, the NP method has been employed for several decades to measure snow thermal conductivity [START_REF] Jaafar | Thermal conductivity of snow by a transient state probe method[END_REF][START_REF] Sturm | Thermal conductivity measurements of depth hoar[END_REF]Domine and others, 2015). As explained previously, Eqn (5) only applies if the heating time is long enough for Eqn (1) to reduce to a simple logarithm (Eqn (4)). However, several practical factors limit the heating time usable in the case of snow. Most importantly, as reported by [START_REF] Sturm | Thermal conductivity measurements of depth hoar[END_REF], the heating of snow near the probe can trigger air convection in the pore space of the snow, which increases heat transport away from the needle, and slows down the temperature increase of the probe. Thus the presence of convection may lead to an overestimation of the thermal conductivity of the snow sample. The onset of convection can be identified as a decrease of the slope of the ΔT against ln(t) curve (e.g. Figure 3 of [START_REF] Sturm | Thermal conductivity measurements of depth hoar[END_REF]. The occurrence and timing of convection events depends on the snow type and the supplied heating power. They are therefore not systematic, but some have been reported to occur after about 50 s of heating (e.g. [START_REF] Sturm | Thermal conductivity measurements of depth hoar[END_REF]Domine and others, 2015). Besides convection, a too long heating cycle would also increase the snow temperature sufficiently to induce snow metamorphism or even melt, compromising the integrity of the snow sample. For these reasons, the heating cycle for measuring snow with the NP method is usually kept within 100 s (Morin and others, 2010;[START_REF] Riche | Thermal conductivity of snow measured by three independent methods and anisotropy considerations[END_REF]Domine and others, 2015). The heating curve data are thus usually analyzed starting around 30 s and until the onset of convection or the end of the 100 s heating cycle [START_REF] Sturm | Thermal conductivity measurements of depth hoar[END_REF]Morin and others, 2010;Domine and others, 2015). However, discrepancies between thermal conductivity values estimated with the NP method and values estimated with other techniques have been reported in the last decade [START_REF] Calonne | Numerical and experimental investigations of the effective thermal conductivity of snow[END_REF][START_REF] Riche | Thermal conductivity of snow measured by three independent methods and anisotropy considerations[END_REF]. One of these independent techniques is based on the numerical simulation of heat conduction through a measured 3D snow microstructure. The effective thermal conductivity can be obtained as the ratio of the modeled heat flux over the temperature gradient across the sample [START_REF] Riche | Thermal conductivity of snow measured by three independent methods and anisotropy considerations[END_REF]. Note that with this methodology, the energy transported as latent heat by the water vapor is usually not taken into account. Using this type of method, [START_REF] Calonne | Numerical and experimental investigations of the effective thermal conductivity of snow[END_REF] report thermal conductivity values that are greater for a given density than that derived from the density-thermal conductivity curve obtained by Sturm and others (1997) and based on over 500 values obtained using the NP method. Moreover, Calonne and others (2011) characterized two samples with both the NP method and numerical simulations, resulting in the same discrepancy, namely that the NP method yields thermal conductivity values underestimated by about 25%. In accordance with [START_REF] Calonne | Numerical and experimental investigations of the effective thermal conductivity of snow[END_REF], [START_REF] Riche | Thermal conductivity of snow measured by three independent methods and anisotropy considerations[END_REF] confirmed that NP measurements yielded lower values compared to numerical simulations, even when applied to the same snow sample. They report NP values that are on average 35% lower than numerical simulation results. However, when used under proper condition the TP02 needle probe used in these studies is expected to have an accuracy of 3% and a repeatability of 1% (TP02 Manual, 2010). These deviations are relatively small and cannot explain the systematic large underestimations reported by [START_REF] Calonne | Numerical and experimental investigations of the effective thermal conductivity of snow[END_REF] and [START_REF] Riche | Thermal conductivity of snow measured by three independent methods and anisotropy considerations[END_REF].

The reasons for these biases have not been clarified so far. [START_REF] Riche | Microstructural change around a needle probe to measure thermal conductivity of snow[END_REF] argued that the damage caused to the snow microstructure by the insertion of a needle may create an air gap between the probe and the snow, potentially affecting the estimation of the thermal conductivity. However, based on 2D numerical simulations, Morin and others (2010) suggested that an air gap around the needle does not affect the measurement as long as at least 30 s of heating time have elapsed. [START_REF] Riche | Thermal conductivity of snow measured by three independent methods and anisotropy considerations[END_REF] also argued that the homogeneous nature of snow assumed to derive Eqn (1) is questionable and might cause some deviations from the NP theoretical background. Indeed, [START_REF] Riche | Thermal conductivity of snow measured by three independent methods and anisotropy considerations[END_REF] pointed out that the typical size of a snow grain (from 0.2 to 2 mm; [START_REF] Fierz | The International Classificationi for Seasonal Snow on the Ground[END_REF]) is close to the typical radius of the probe (0.75 mm), potentially invalidating the homogeneous medium assumption. Finally, the theoretical work of [START_REF] Jaeger | Conduction of heat in an infinite region bounded internally by a circular cylinder of a perfect conductor[END_REF] assumes an isotropic medium, while it is known that snow thermal conductivity exhibits some degree of anisotropy [START_REF] Calonne | Numerical and experimental investigations of the effective thermal conductivity of snow[END_REF][START_REF] Riche | Thermal conductivity of snow measured by three independent methods and anisotropy considerations[END_REF]. Following the work of Grubbe and others (1983), [START_REF] Riche | Thermal conductivity of snow measured by three independent methods and anisotropy considerations[END_REF] argued that for horizontally inserted probes, Eqn (5) applies, but with k being the geometric mean of the horizontal and vertical thermal conductivities k x y and k z , i.e. k = k x yk z . However, accounting for anisotropy does not explain the systematic discrepancy between NP measurements and numerical simulations [START_REF] Riche | Thermal conductivity of snow measured by three independent methods and anisotropy considerations[END_REF].

The goal of this work is to systematically evaluate the potential sources of bias in the needle probe technique, such as those suggested above. In particular, we evaluate the over-simplification of the heating curve to a simple logarithm, the validity of the homogeneity assumption for snow, the impact of structural damages done during the probe insertion, and the impact of ice crystallization below a needle probe buried in a snowpack for season-long thermal conductivity monitoring. We finally propose a correction of NP measurements, applicable to needle probes buried, rather than inserted, in snowpacks. To these ends, we performed a set of NP measurements in a cold room using both hand-inserted and fixed probes, buried in the snow for an extended period. These experimental results, and their potential underestimations, are discussed in the light of computed microtomography images complemented with numerical simulation of energy transport, as described in the following methods section.

Methods

Needle probes used

The experimental measurements done for this study were performed using two TP02 probes manufactured by Hukseflux, the same model as in Morin and others (2010), [START_REF] Riche | Thermal conductivity of snow measured by three independent methods and anisotropy considerations[END_REF] or Domine and others (2015). The TP02 model respects the ASTM D5334-14, ASTM D5930-97 and IEEE 442-1981 standards on the NP method. This model consists of a 15 cm long needle, heated on the 10 cm part adjacent to the handle of the needle. As seen in Figure 1a, two thermocouples are inserted in the needle, one at the tip in the non-heated zone of the needle, where the temperature does not increase during measurements, and one in the middle of the heated zone (also displayed in Figure 0.1 of TP02 Manual, 2010). With these thermocouples, the temperature increase due to the heating can be monitored. The TP02 probe is comprised of a stainless tube (1.5 mm external diameter and 0.2 mm thickness), encapsulating two heating constanstan wires (0.127 mm diameter) with PTFE jackets (0.0865 mm thickness). The remaining space of the probe is filled with glass beads (information provided by Hukseflux in January 2020). Based on this information, we estimated that the TP02 probe has a density of 4500 kg m -3 and a specific thermal capacity of 500 J kg -1 . We evaluated that during a measurement the TP02 probe is equivalent to a homogeneous and highly conducting cylinder with the given density and specific thermal capacity, and that the temperature increase during the measurement of extruded polystyrene foam is very well represented by Eqn (1) (Appendix A). Indeed, the TP02 needle presents a very favorable length to diameter ratio of 66 (Kömle and others, 2011), which ensures that the heated zone of the probe can be treated as infinitely long, and thus that Eqn (1) applies [START_REF] Blackwell | The axial-flox error in the thermal-conductivity probe[END_REF].

During a measurement, and similarly to previous snow studies (e.g. Morin and others, 2010;[START_REF] Riche | Thermal conductivity of snow measured by three independent methods and anisotropy considerations[END_REF], the probe was heated with a constant power Q for a total duration of 100 s. The precise heating power was monitored, and set to be near 0.4 W m -1 when measuring snow samples. This typically led to a temperature increase of about 1 K in the center of the heated zone. The temperature was also recorded during 100 s prior to the heating cycle, in order to establish a baseline. This baseline was used to estimate the initial temperature and to compute a potential trend. Indeed, even though care was taken to ensure the isothermal condition of the probe before the measurements, some datasets exhibited a non-negligible trend during the 100 s leading to the heating cycle, as seen in Figure 1b. Similarly to Kömle and others (2008), the data during the heating cycle were detrended. Note that overcorrecting the temperature trend will tend to under or overestimate the measured thermal conductivities depending on the sign of the temperature trend. If possible, users of the NP should therefore use data not affected by a temperature trend. However, the potential errors produced by overcorrection in this article remain mostly below 5%, and are therefore small compared to the other errors discussed in the article. Thus, performing the data analysis without detrending does not fundamentally change the results and conclusions of the article.

The resulting temperature increase was then expressed as a function of the logarithm of time, and the heating curve data were selected over a given time window. The thermal conductivity value was estimated using the slope of the selected window and Eqn (5), as exemplified in Figure 1c. In this article, the selected data systematically correspond to a time window of 30 s, centered around a given time referred to as the measurement time. In principle, if the ΔT against ln(t) curve is linear, the estimated slope should be independent of the width of this time window. We found that 30 s of data is long enough to limit the influence of experimental noise on the slope estimation, and short enough so that the slope of the ΔT against ln(t) curve does not vary much over this interval. Increasing or decreasing the width of this time window has little effect on the slope estimation of noisefree data (produced using the analytical Eqn (1)), if the measurement time is larger than 30 s. Again, processing the data with a different time window does not fundamentally change the results of the article.

Cold-room gradient temperature-controlled experiment

For the purpose of this study a gradient box experiment was set-up in a cold room. The gradient box is a thermally insulated box that can be filled with snow, and whose bottom and top temperatures can be imposed to create a thermal gradient. In our case, 15 cm of light snow comprised of decomposing and fragmented precipitation particles (density of 125 kg m -3 and specific surface area of 40 m 2 kg -1 ) was sieved into the box, burying a fixed needle placed at mid-height, and let to evolve for 3 weeks under a thermal gradient of about 50 K m -1 and with a mean temperature close to 268 K, as displayed in Figure 2. We observed the snow microstructure and its evolution during the 3 weeks of experiment and saw its metamorphism to depth hoar. The classification of the snow types was done according to [START_REF] Fierz | The International Classificationi for Seasonal Snow on the Ground[END_REF].

During the first 2 weeks of metamorphism, manual NP measurements were also performed in the gradient box. The probe was either inserted through a pre-drilled hole in one of the walls, or freely inserted in the snow after the removal of one of the walls. The fixed probe buried during the sieving of the snow was meant to reproduce the use of NPs that are set-up in the field prior to snow falls and subsequently stay embedded in the snow for the rest of the season (Morin and others, 2010; 1) with a zero contact thermal resistance (in orange). Domine and others, 2015). We recall that our measurements were performed following the typical methodology used in previous snow studies (i.e. a heating power of 0.4 W m -1 and a heating cycle of 100 s).

As seen in Figure 3, the temperature increase curves obtained with the fixed needle are reproduced by the application of Eqn (1) with a zero contact thermal resistance and reasonable values of density and thermal conductivity.

X-ray computed microtomography

Snow was also regularly sampled from the gradient box with a cylindrical corer made out of polymethylmethacrylate (PMMA, also known under the commercial name Plexiglas), designed specifically to be placed in an X-ray tomograph. Every time snow was sampled from the gradient box, a computed microtomography (μCT) scan was performed, allowing us to access the snow microstructure and its evolution during the 3 weeks of experiment. The μCT scans were performed with a 15 μm resolution for all samples, and were subsequently segmented following the methodology of Hagenmuller and others (2016) to produce binary images of ice and air. Note that after about a week of metamorphism the snow became very brittle in the gradient box and even though care was taken during the sampling process, we sometimes observed collapses with a clear settlement of the snow in the PMMA core.

In addition to NP measurements in the gradient box, NP measurements were also performed in the snow sampled in the PMMA core, after some of the μCT scans. An insertion guide was used to perform the measurements in this case, in an effort to minimize the structural damages done to the snow. The insertion guide consists of a rail at the end of which the needle probe can be held in place by an operator. The PMMA core can then be glided on the rail to insert the probe in the snow. Schematics of the gradient box, with the position of the measurements performed in it, are given in Figure 2.

Finite element simulations

The μCT binary images were exploited with numerical simulations performed using the Finite Element Method (FEM) and the open-source software ElmerFEM [START_REF] Malinen | Elmer finite element solver for multiphysics and multiscale problems[END_REF]. These numerical simulations were used for two main purposes: (i) to estimate the thermal conductivity of a snow sample based on its microstructure and (ii) to simulate the heat curves of NP experiments based either on a full snow microstructure or on an equivalent homogeneous medium. All the simulations were performed at a temperature of 268 K, similarly to the coldroom measurements. In total, three distinct types of numerical simulations, schematized in Figure 4, were carried out: 3D homogenization simulation A first type of simulation was used to estimate the effective thermal conductivity of a snow sample using its μCT-based microstructure, similar to [START_REF] Calonne | Numerical and experimental investigations of the effective thermal conductivity of snow[END_REF] or [START_REF] Riche | Thermal conductivity of snow measured by three independent methods and anisotropy considerations[END_REF]. For each sample, a tetrahedral mesh was first obtained using a sub-cube of the segmented μCT scan and the CGAL meshing library. A steady-state simulation of heat conduction under a thermal gradient was performed, by imposing the top and bottom temperatures of the snow sample and enforcing zero-flux conditions on the remaining sides [START_REF] Riche | Thermal conductivity of snow measured by three independent methods and anisotropy considerations[END_REF]Fourteau and others, 2021). As in [START_REF] Riche | Thermal conductivity of snow measured by three independent methods and anisotropy considerations[END_REF], the effective thermal conductivity was then obtained as the ratio of the resulting heat flux to the imposed gradient. In order to capture the anisotropy of snow, these steady-state simulations were performed in the vertical and horizontal directions for each snow sample. These simulations are referred to as 3D homogenization simulations in the rest of the article, as they yield the effective thermal conductivity of a snow sample when viewed as a homogeneous medium with a given thermal conductivity relating the macroscopic heat flux to the macroscopic thermal gradient. Simulations were performed on cubes with side length between 6 and 9 mm, in order to use samples larger than the Representative Elementary Volumes (REV; [START_REF] Calonne | Numerical and experimental investigations of the effective thermal conductivity of snow[END_REF].

3D heterogeneous NP simulation

The second type of simulation was used to simulate the results of a transient NP measurement in snow seen as a heterogeneous mixture of ice and air. For this, the probe was heated with a constant power Q, and its simulated temperature increase recorded every second. These simulated data were then analyzed similarly to actual NP measurements data. In this case, it was necessary to perform a complete 3D simulation to properly represent the heterogeneous nature of the medium. The probe was numerically inserted in the segmented images of snow samples, thus ensuring a perfect thermal contact between the needle and the snow. Similarly to the first type of simulation, the 3D meshes were produced from segmented images using the CGAL meshing library. The transient heating of the needle was thus simulated, accounting for heat transport through both the ice and the air.

2D homogeneous NP simulation

A third type of simulation was finally used to simulate the results of a transient NP measurement, this time performed in a homogeneous medium. As with the 3D heterogeneous NP simulation, the temperature increase of the probe was recorded every second, and the data analyzed as for actual NP measurements. As the problem exhibits a translational invariance along the probe axis, it can be simplified to a 2D version in the plane perpendicular to the probe axis. This type of simulation allowed us to easily include potential defects at the probe/snow interface, such as air gaps created during the needle insertion. Finally, in the case where the thermal contact between the probe and the investigated medium is perfect, the 2D simulations yielded the exact same results as the analytical formula of Eqn (1) with h = 0. The 2D meshes used in the simulations were produced using the gmsh software.

In snow heat is transported by conduction through the ice and air phases. However, the motion of water vapor carrying latent heat, also leads to heat transport (Yosida and others, 1955;[START_REF] Vries | The theory of heat and moisture transfer in porous media revisited[END_REF][START_REF] Hansen | A macroscale mixture theory analysis of deposition and sublimation rates during heat and mass transfer in dry snow[END_REF]. The importance of these latent heat processes depends on the surface kinetics of water molecule sublimation and condensation (Fourteau and others, 2021). If water molecules are rapidly accommodated and sublimated, then heat transport occurs as if the air had an extra thermal conductivity of magnitude D 0 L(dc sat /dT ) (Yosida and others, 1955;[START_REF] Vries | Simultaneous transfer of heat and moisture in porous media[END_REF][START_REF] Vries | The theory of heat and moisture transfer in porous media revisited[END_REF]Moyne and others, 1988), where D 0 = 2 × 10 -5 m 2 s -1 is the diffusion coefficient of water vapor in the air (Calonne and others, 2014), L = 28 × 10 5 J kg -1 is the latent heat of sublimation of water [START_REF] Lide | CRC handbook of Chemistry and Physics, Properties of Ice and Supercooled Water[END_REF], and dc sat /dT is the derivative of the saturation concentration of water vapor with respect to temperature, estimated to be 2.6 × 10 -4 kg m -3 K -1 at 268 K. On the contrary, if the surface kinetics of water molecules deposition is very slow, then latent effects can be neglected and snow behaves as an inert bi-phasic medium. Based on measured variations of the thermal conductivity of snow with temperature, Fourteau and others (2021) suggested that the fast kinetics hypothesis is better suited to model the thermal conductivity of snow, rather than the slow kinetics one. Therefore, in this work we consider that snow is characterized by the infinitely fast kinetics hypothesis and we included latent heat effects under the form of an increased air thermal conductivity. In particular, this means that for the 3D homogenization simulations, the thermal conductivity of air is increased to take into account latent heat effects, in contrast to [START_REF] Calonne | Numerical and experimental investigations of the effective thermal conductivity of snow[END_REF] and [START_REF] Riche | Thermal conductivity of snow measured by three independent methods and anisotropy considerations[END_REF]. The fast kinetics assumption has no influence on the 2D homogeneous NP simulations, as in these we directly consider snow as a homogeneous medium with a given thermal conductivity.

Results

Potential sources of error for the NP method in snow

Deviation of NP measurements from the logarithmic regime As previously mentioned, to be applicable, the NP method requires the heating time to be long enough for Eqn (1) to simplify to the simple logarithmic Eqn (4). With snow, it is thus necessary that the convergence toward the logarithm regime occurs before the onset of convection or the end of the 100 s heating cycle. To our knowledge the convergence toward a logarithm regime has been justified with two reasons. First, the selected portion of the ΔT against ln(t) curve is visually linear. Second, the characteristic time τ = r 2 /χ appearing in the analytical Eqn (1) and which partly governs the rate of convergence to the logarithmic regime, can be estimated to be of the order of a few seconds only. One may therefore expect the heating curve to be close to a logarithm after about 30 s, since t ≫ τ (Morin and others, 2010).

To assess whether such a convergence actually occurs within 100 s in snow, one can compare the complete analytical formula of Eqn (1) with the simplified versions of Eqns (3) and (4). To apply these formulas, we considered that the NP measurements were made with a TP02 probe heated at 0.4 W m -1 in a snow with a thermal conductivity of 0.15 W K -1 m -1 , a density of 200 kg m -3 , a specific thermal capacity of 2000 J kg -1 (Morin and others, 2010), and without any thermal contact resistance between the needle and the snow. The parameters to use in Eqns (1), ( 3) and ( 4) thus have values of τ = 1.5 s, α = 0.356, and h = 0.

Figure 5 displays the temperature increase given by Eqns (1), ( 3) and (4) over a 300 s period. It indicates that while the exact formula and the asymptotic development converge within the first 30 s, they are still significantly different from a simple logarithm in the first 100 s. Unsurprisingly, this affects the slope of the ΔT against ln(t) curves. In particular, Figure 6 displays the thermal conductivity values deduced if one were to apply the standard Eqn (5) to the blue curve, which represents the true theoretical heating curve. As explained earlier, the slope of the ΔT versus ln(t) curve is obtained by fitting a straight line over a 30 s time window centered around a time t mes , referred to as the measurement time. Results show that for measurement times between 30 and 85 s, the slope of ΔT versus ln(t) given by the exact formula is larger than that expected with the pure logarithmic form, which leads to an underestimation of the thermal conductivity, by about 15% in this case. In agreement with the idea of a gradual convergence toward a pure logarithm, the more time is elapsed the better is the estimation of the thermal conductivity.

Theory therefore indicates that in the case of snow measured with a TP02, a heating cycle of 100 s is by far not long enough to reach the logarithmic regime and for Eqn (5) to apply. This runs against common practice in snow studies and against the reasons provided at the beginning of the section for the validity of this assumption, namely that in the case of snow measurements (i) the ΔT against ln(t) curve appears to have reached linearity after ∼ 30 s, and (ii) t ≫ τ after 30 s. Our understanding is that these reasonings are not valid as (i) visually judging the linearity of a curve is subjective, especially with real measurement data that are subjected to noise and potential convection starting before 100 s, (ii) Eqn (3) suggests that the decrease of high order terms in the asymptotic development follows 1/t, which is relatively slow to decrease (as seen in Fig. 5). In other words, heating curve data obtained before 100 s should be considered as shorttime, and consequently Eqn (5) does not directly apply.

As highlighted in Figure 7, the degree of error made by applying Eqn (5) to short-time measurements strongly depends on the parameter α, directly related to the ratio of the volumetric thermal capacity of the investigated medium over that of the probe. Media with low thermal capacities compared to that of the probe (small α) are characterized by large underestimations, while media with high thermal capacities compared to that of the probe (high α) are rather characterized by overestimations. For snow densities ranging from 100 to 500 kg m -3 and measured with a TP02 probe, the parameter α ranges from 0.18 to 0.89. This indicates that while the snow thermal conductivity tends to be systematically underestimated with the NP method, the degree of underestimation is quite sensitive to snow density.

Impact of the heterogeneous nature of snow

We have seen that the limited duration of NP measurements in snow could explain the underestimations reported by [START_REF] Calonne | Numerical and experimental investigations of the effective thermal conductivity of snow[END_REF] and [START_REF] Riche | Thermal conductivity of snow measured by three independent methods and anisotropy considerations[END_REF]. Yet, other mechanisms may also contribute to thermal conductivity NP measurement errors. Notably, [START_REF] Riche | Thermal conductivity of snow measured by three independent methods and anisotropy considerations[END_REF] proposed that the discrepancy between the NP and other methods could be due to the heterogeneous nature of snow. Indeed, the theoretical development of [START_REF] Jaeger | Conduction of heat in an infinite region bounded internally by a circular cylinder of a perfect conductor[END_REF] is derived assuming that the external medium is homogeneous, and it is not clear whether the calculations apply to a medium such as snow.

To test this hypothesis, we compared simulations of NP measurements done in an actual snow sample, heterogeneous by nature with its complex microstructure, and in its equivalent homogeneous medium. For the simulation of the heterogeneous medium, we performed a 3D FEM transient simulation of NP heating within a snow microstructure obtained by μCT. The selected snow sample corresponds to a depth hoar snow type [START_REF] Fierz | The International Classificationi for Seasonal Snow on the Ground[END_REF], with a density of 175 kg m -3 and a specific surface area of 18 m 2 kg -1 . We chose this sample as it has grains with a size similar to that of the needle radius, and is thus a good candidate to quantify the potential impact of heterogeneities on the NP measurement. In order to have an as large as possible simulation domain, and thus mitigate boundary effects, the original μCT image (1.35 cm side-length) has been spatially extended thanks to plane symmetries, as shown in Figure 8. For the simulation of the NP measurement in the homogeneous medium, one needs to prescribe an effective thermal conductivity to the equivalent homogeneous medium. The effective thermal conductivity of the equivalent homogeneous medium has been estimated to 0.125 W K -1 m -1 , thanks to a 3D homogenization simulation performed on the central, non-duplicated part of Figure 8. We recall that in our case, all simulations include the effect of latent heat transport, as described in the methods section.

Outputs of the simulations in the heterogeneous and homogeneous media are displayed in Figure 9. Results indicate that the measurements performed on heterogeneous and homogeneous media lead to mostly similar estimations of the effective thermal conductivity, with a maximal difference between the two cases of about 5% for a measurement time of 30 s. This difference appears to diminish with larger measurement times. In order to confirm this result, we performed the same comparison of heterogeneous and homogeneous simulations for a NP measurement performed in a rounded grains snow sample, characterized by a density of 280 kg m -3 and a specific surface area of 23 m 2 kg -1 . As with the depth hoar snow, we observe that the heterogeneous simulation results in an overestimation of <5% compared to the homogeneous simulation.

Thus, while there exist discrepancies between our two types of simulation, these differences are small compared to the overall underestimation caused by the insufficient heating time during NP measurements. The systematic underestimation of the NP method thus cannot be attributed to the heterogeneous nature of snow, and to a first approximation it appears reasonable to assume that during a NP measurement snow behaves as a homogeneous medium, with a thermal conductivity consistent with Fig. 6. Thermal conductivity values obtained by fitting a straight line between ΔT and ln(t) over a 30 s time window centered around a measurement time t mes . The temperature increase is computed using the analytical formula of Eqn (1) and the thermal conductivity is deduced by applying Eqn (5) to the temperature increase. The snow sample has a thermal conductivity of 0.15 W K -1 m -1 , represented as a dashed line. The green zone corresponds to the typical measurement times used in snow studies.

steady-state numerical simulations. Furthermore, this implies that simulations of NP measurements in snow can be performed with a reasonable accuracy using a simple homogeneous medium instead of a complex 3D microstructure.

Impact of microstructural damages: theory and experiment

When inserting a NP into a snow layer, the initial snow microstructure is partially damaged. Using μCT, Riche and Schneebeli (2010) report that damage zones of a few tenths of a millimeter were formed when inserting graphite rods with a 0.6 mm radius in snow, creating a poor thermal contact between the two. If the heating cycle is long enough to reach the logarithmic regime, the presence of a thermal gap between the NP and the investigated medium should not influence the result. However, as seen above the logarithmic regime is not reached for NP measurements lasting 100 s or less. Equation (3) indicates that before the logarithmic regime, the presence of a heat gap affects the increase of temperature with time, and therefore affects the estimation of the effective thermal conductivity. In particular, a non-zero h parameter (i.e. a non-zero thermal contact resistance between the probe and the snow) in Eqn (3) increases the derivative of ΔT with respect to time, and thus leads to a further underestimation of the effective thermal conductivity when Eqn ( 5) is applied.

Yet, Morin and others (2010) report, using numerical modeling, that the presence of an air gap around the needle does not affect the estimation of the thermal conductivity, in apparent contradiction with the theoretical work of Jaeger ( 1956). An analysis of the computer program used by Morin and others (2010) (S. Morin, personal communication, December 2019) reveals that in this numerical model the entirety of the heat produced in the needle is transferred to the surrounding medium, without any heat storage in the needle. This is equivalent to setting a zero heat capacity for the needle probe, i.e. α → ∞. Letting α go to ∞ in Eqn (3) leads to the following asymptotic development for the temperature increase over time:

DT = Q 4pk 2h + ln 4 C t t - 1 2 t t + 1 2 t t ln 4 C t t . ( 6
)
The parameter h now only appears as an additive constant, and does not affect the increase of ΔT over time. In this peculiar case, the presence of an air gap does not influence the estimation of the thermal conductivity, consistent with the modeling of Morin and others (2010). However, this case does not apply for NP measurements in snow, as the thermal capacity of snow is small relatively to that of the probe.

To assess the influence of the presence of an air gap around a TP02 needle probe, we performed 2D transient simulations of NP measurements, with the presence of air gaps around the probe. The simulations cover three types of snow and four air gap sizes around the needle. The results are displayed in Figure 10, and indicate that the presence of an air gap greatly influences the estimation of the effective thermal conductivity, especially for the denser and more conductive snow. Air gaps of 0.2 mm, in line with the values reported by [START_REF] Riche | Microstructural change around a needle probe to measure thermal conductivity of snow[END_REF], can lead to an underestimation much larger than in the case of a perfect thermal contact. The idea of larger underestimations due to the presence of air gaps is consistent with the measurements we performed during the gradient box experiment. We recall that during this experiment we performed NP measurements (i) with a fixed needle probe, (ii) with manually inserted probes with the help of an insertion guide, and (iii) with manually inserted probes without the help of a guide. The NP data were analyzed by fitting a straight line in the ΔT against ln(t) curve, in a 30 s window centered around a 45 s measurement time. We selected a 45 s measurement time as this is a time typically chosen by previous snow studies (e.g. Morin and others, 2010;[START_REF] Riche | Thermal conductivity of snow measured by three independent methods and anisotropy considerations[END_REF]. This allows an easier comparison of our results with the published literature. Nonetheless, we have tested that using a different measurement time of 65 s or 85 s leads to equivalent results. The different experimental estimations of the snow thermal conductivity are displayed in Figure 11. Due to hand movements while manually holding the needle probe during measurements without the insertion guide, the air gaps around the needle were usually much larger than the 0.2mm gap used in Figure 10. Therefore, the expected underestimation due to the presence of large air gaps are larger than those presented in Figure 10. The results show that the data with manually inserted probes, prone to microstructural damages, are systematically lower than those obtained with the fixed needle data. Moreover, the use of a guide to minimize structural damages during the needle insertion tends to improve the measurements, even though it does not entirely solve the problem. Note that these results are consistent with the recent review of the use of NP in soil by He and others (2018), who report that 'data for the first 5 to 100 s, depending on the duration of the [heating cycle] and the probe dimensions, are affected by the probe characteristics and the contact resistance between the probe and the soil'.

Impact of ice recrystallization below fixed needle probes: theory and observation

When a fixed NP is buried within a snowpack under a macroscopic thermal gradient, some of the water vapor diffusing in the pores recrystallizes below the needle. It is therefore possible that the microstructure near the needle is not representative of the microstructure of the rest of the snowpack. In such conditions, one might wonder how the effective thermal conductivity estimated with the NP relates to that of the undisturbed snow, far from the needle.

To test this hypothesis we performed a 2D numerical simulation of a NP measurement in homogeneous snow, assuming that of a 1 mm thick block of pure ice crystallized below the needle. The presence of a 1 mm thick block of ice below the needle is the worst case scenario and maximize the potential effect of ice accretion below the needle. The comparison with a simulation not including ice accretion, displayed in Figure 12, indicates that the estimation of the effective thermal conductivity is only slightly influenced by the presence of such a block of ice.

Moreover as seen in Figure 13, at the end of our gradient box experiment we carefully excavated the fixed needle probe in order to quantify the nature of the crystals formed around the needle. After visual inspection, no clear differences were observed between these crystals and the rest of the snowpack. A μCT scan of the needle with ice crystals still attached to it was also performed and confirmed that the crystals attached to the needle were depth hoar crystals, similar to those in the rest of the snowpack.

For all these reasons, we believe that the snow microstructure surrounding a fixed NP is not impacted enough to significantly affect the effective thermal conductivity estimation.

Correcting the underestimation of NP measurements in snow

We have seen that the limited duration of the heating cycle leads to a systematic underestimation of the effective thermal conductivities of snow if one directly applies Eqn (5), even in the case of a perfect thermal contact between the probe and the snow. Yet, the NP technique remains of particular interest in snow sciences as it is currently the only available experimental method to easily estimate snow thermal conductivity in the field, as well as the only method that can be automated to produce in-situ measurements over an entire snow season. This is why we propose here a methodology to correct the results given by the application of Eqn (5) at short times. This correction, with results presented here for T = 268 K and subsequently extended to other temperatures, is based on the following assumptions:

(1) There is no thermal contact resistance between the snow and the needle (i.e. h = 0). (2) The needle probe is a TP02 model, with the material properties validated in Appendix A.

(3) The snow effective thermal conductivity is related to its density, following the formula given by Fourteau and others (2021) at 268 K.

The protocol to perform a corrected measurement is simply to:

(1) Run a heating cycle at constant power per meter Q.

(2) Estimate a first value of the effective thermal conductivity k raw using Eqn (5), as with standard NP measurements (e.g. [START_REF] Sturm | Thermal conductivity measurements of depth hoar[END_REF][START_REF] Riche | Thermal conductivity of snow measured by three independent methods and anisotropy considerations[END_REF]. This estimation of k raw is performed by fitting a slope between ΔT and ln (t) in a time window centered around a measurement time denoted t mes . This values k raw thus corresponds to an underestimated value.

(3) Correct k raw with a correction factor C(t mes , k raw , T ), giving the corrected value k corrected = C k raw . Note that the correcting factor C depends both on the measurement time t mes and on k raw .

The correcting factor C has been empirically estimated by applying the analytical formula of Eqn (1). Practically, heating curves of NP measurements in homogeneous snow have been computed for effective thermal conductivities ranging from 0.054 to 0.50 W K -1 m -1 (corresponding densities ranging from 50 to 450 kg m -3 ). For each simulation, the correcting factor at time t mes has been estimated as the ratio between the true thermal conductivity, and the estimation k raw resulting from the application of Eqn (5) at measurement time t mes . We thus obtained a lookup table giving C(t mes , k raw , T ) for an array of t mes and k raw . Values of C can then be interpolated over an entire continuous range of t mes and k raw . A visualization of the correcting factors depending on k raw and t mes is given in Figure 14. It indicates that the correcting factor exceeds 1.6 for small t mes and low k raw . For large t mes and k raw it is <1.05. The correction factor, and thus the degree of underestimation, shows a general decrease with larger t mes and larger k raw . This is to be expected as larger t mes results in a measurement performed closer to the logarithmic regime where Eqn (5) applies, and larger k raw implies a denser snow, and thus a larger alpha parameter in Eqn (1). The correction table is available as a Supplementary Material. We also propose correction tables for measurements performed at other temperatures, ranging from 223 to 273 K (based on the densitythermal conductivity relationships proposed by Fourteau and others (2021) at various temperatures).

We tested our correction protocol on the needle probe measurements done with the fixed needle probe in the gradient box. The validity of the correction is assessed by comparing the corrected values with effective thermal conductivities deduced from 3D homogenization simulations performed on snow microstructures measured on the same day as the NP data. As seen in Figure 15a, application of the correction to the NP data leads to an increase of the estimated thermal conductivities of about 30 %. During the first 10 days the corrected NP values are well in line with the μCT-based estimations made under the fast kinetics hypothesis. Figure 15 also displays the μCT-based thermal conductivity estimations under the slow kinetics hypothesis, which are discussed in the discussion section, but we recall that the fast kinetics hypothesis appears more suited to model snow thermal conductivity (Fourteau and others, 2021). This suggests that our correction protocol provides adequate thermal conductivity values.

After 10 days of experiment, the μCT-based estimations of the thermal conductivity display an overall increase accompanied by large variations. Such an increase is not reflected in the NP measurements, which remain stable with and without correction. Our understanding is that the increase and variations of the μCT-based estimations reflect the partial damaging of the snow during sampling, rather than real variations of the snow thermal conductivity in the gradient box. Indeed, as discussed in the methods section, we observed collapses while coring welldeveloped depth hoar snow. This led to an increased snow density in the core compared to that of the gradient box, which explains the high thermal conductivities deduced from μCT scanning. This hypothesis is confirmed by the fact that the variations of the μCT-based thermal conductivities are highly correlated with the density variations of the snow in the PMMA core, as shown in Figure 15b.

For this article, the needle probe measurements in the gradient box were performed with a measurement time of 45 s as it is close to the value chosen by previous snow studies. Nonetheless, we have also tested our correction protocol with different measurement times of 65 and 85 s. Both cases lead to results similar to those obtained with a measurement time of 45 s.

Finally, for our correction protocol we made the choice to use a methodology as close as possible to standard NP measurements, which permits a direct comparison with previous published NP data in snow. An alternative correction method would be to fit the whole temperature increase curve with the thermal conductivity as the sole free parameter (thus still assuming that thermal conductivity and density are related by the formula of Fourteau and others, 2021). This method has the advantage of utilizing the whole temperature curve, which could potentially yield more robust results. As a test, we implemented such a method, and observed results in agreement with those of the correction protocol proposed here for temperature curves measured both before and after 10 days.

Discussion

Hand-used and fixed needle probes

As discussed above, non-quantifiable microstructural damages due to the insertion of the needle in snow might lead to a significant underestimation of the effective thermal conductivity. While it is possible to limit those damages with careful manipulation, it appears unlikely to totally avoid them, especially in brittle snow such as depth hoar. On the other hand, the use of a fixed NP does not lead to such microstructural damages. Moreover, we have seen that (i) the ice crystals formed below a fixed NP have characteristics similar to those of the rest of the snowpack and (ii) even the worst case scenario of a block of pure ice forming below the needle does not strongly perturb the effective thermal conductivity measurements. We therefore encourage the use of fixed needle probes rather than hand-used ones whenever possible. That being said, users should be aware that fixed NP data will still systematically underestimate the effective thermal conductivity of snow, unless a correction is specifically applied to mitigate this problem.

Limitations of the proposed correction

As the NP is currently the only easily applicable and easily automated method for thermal conductivity measurements in the field, a correction for NP measurements in snow has been [START_REF] Fierz | The International Classificationi for Seasonal Snow on the Ground[END_REF], with a transition from decomposing and fragmented precipitation particles in green to depth hoar in blue. Note that we did not observe the faceted crystals transition, as the snow had already metamorphized to depth hoar for our first sampling after 3 days.

proposed in this article. Here, we discuss the limitations of this correction, as well as those of three hypotheses used to derive the correction.

We first had to assume a perfect thermal contact between the snow and the probe. It therefore follows that the proposed correction only applies for fixed needle probes, and cannot be used to correct hand-done measurements. Indeed, in the latter case the damages done to the microstructure greatly influence the temperature increase of the probe, as seen in the results section, and our computations no longer apply. Furthermore, the precise damage done to the snow is essentially non-quantifiable during a real NP measurement. This means that the degree of underestimation due to these damages is unknown and therefore cannot be corrected.

We then assumed that the probe had the physical properties validated for the TP02 in Appendix A, and our correction thus does not apply to other NP models. Still, similar correction tables for other probe models could easily be obtained, as long as the physical properties of the probe used (i.e. its radius and volumetric thermal capacity) are known.

Finally, we assumed that the density and the thermal conductivity of snow were related by the formula proposed by Fourteau and others (2021). This hypothesis is needed to restrain the number of independent snow parameters governing the temperature increase of the probe. Indeed, in the general case the precise temperature of the probe is governed by both the thermal conductivity of snow as well as its density (controlling the snow volumetric thermal capacity). A same slope of the ΔT versus ln(t) curve at a measurement time t mes could therefore be attributed to various snow samples, if thermal conductivity and density are allowed to vary independently. Assuming that density and thermal conductivity are dependent ensures the uniqueness of the snow sample yielding a given ΔT versus ln(t) slope at a given measurement time t mes , a necessary assumption for our correction method. One should therefore be aware that snow samples that deviate from the assumed relationship between density and thermal conductivity would not be perfectly corrected for. The snow thermal conductivity versus density parametrization used in the article has been derived for alpine snow with density mostly ranging between 100 and 350 kg m -3 (Fourteau and others, 2021) and the provided correction table thus primarily applies to these snow types. The thermal conductivity of dense or arctic snow may deviate from the used density parameterization and a specific correction table should be therefore adapted with the described methodology.

As an example of this last limitation, let us assume that we are using the NP method to measure a snow sample characterized by a density of 100 kg m -3 and a given thermal conductivity of 0.055 W K -1 m -1 . This snow deviates from the density-thermal conductivity relationship assumed in our correction protocol, which predicts a thermal conductivity of 0.079 W K -1 m -1 for this density. Performing a standard NP measurement without correction would lead to a value of about 0.04 W K -1 m -1 , i.e. an underestimation of about 30% compared to the true thermal conductivity value. Applying the correction derived in this article leads to a final result of 0.061 W K -1 m -1 , about 11% larger than the true value of 0.055 W K -1 m -1 .

A potential way to relax the necessity assuming a densitythermal conductivity relationship would be to fit the entire ΔT curve against the analytical formula of [START_REF] Jaeger | Conduction of heat in an infinite region bounded internally by a circular cylinder of a perfect conductor[END_REF] and to retrieve both the snow thermal conductivity and the snow density. Yet, it is known that obtaining density from NP measurement is particularly difficult and leads to results with large uncertainties (Waite and others, 2006). As a test, we nonetheless tried to fit the whole temperature curves with density and thermal conductivity as free parameters. Results show that overestimation of density is compensated by underestimation of the thermal conductivity, resulting in both unrealistic values of density and thermal conductivity despite a good fit to the temperature data.

Therefore, users of the proposed correction should be aware of its limitations. While the correction tends to improve the estimated thermal conductivity by accounting for the systematic underestimation of the NP, the correction will poorly work for snow samples that strongly deviate from the assumed relationship between density and thermal conductivity and for experiments with a poor thermal contact between the snow and the probe.

Validity of the fast kinetics hypothesis for energy transport in snow

In this article, we made the assumption that at the microscopic level snow is characterized by a fast adsorption and desorption of water molecules onto ice. In this case, snow behaves as if the thermal conductivity of air were increased compared to its normal value (Yosida and others, 1955;[START_REF] Vries | Simultaneous transfer of heat and moisture in porous media[END_REF][START_REF] Vries | The theory of heat and moisture transfer in porous media revisited[END_REF]Moyne and others, 1988). Justifications for this assumption have notably been proposed by Fourteau and others (2021). We stress that in this article the fast kinetics hypothesis was only required to assume a relationship between snow density and thermal conductivity, and to compare our corrected results to 3D homogenization values. Therefore, the part of the results section about the sources of underestimation of the NP method in snowwould remain valid even if the fast kinetics hypothesis did not apply.

Yet, a further justification for the better suitability of the fast kinetics hypothesis over the slow kinetics one can be provided thanks to our fixed needle probe measurements. A comparison between these fixed NP data and thermal conductivity estimations based on 3D homogenization simulations under the slow kinetics assumption can be made in Figure 15a. This indicates that the uncorrected NP data are in good agreement with the 3D homogenization values under slow kinetics. However, we have seen that even in the case of a perfect thermal contact between the probe and the snow, the uncorrected NP method still yields an underestimated value. We can thus conclude that the 3D homogenization values under slow kinetics are also underestimated. This underestimation likely comes from the absence of latent heat transport under the slow kinetics hypothesis.

Conclusions

This work investigated the origins of the systematic underestimation of needle probe measurements of snow effective thermal conductivity suggested in previous studies. We have shown that the discrepancy originates from the fact that the heating duration used for NP snow measurements is not long enough to reach the logarithmic regime, while this logarithmic regime is at the core of the computations yielding the effective thermal conductivity of snow from the rate of temperature increase of the needle. Moreover, as the logarithmic regime is not reached, the effective thermal conductivity estimation is impacted by the density of snow, the potential presence of a thermal gap around the needle, and the material properties of the needle. In particular, the presence of microstructural damages done to the snow during the insertion of the needle leads to a further underestimation of the effective thermal conductivity. This implies that hand measurements, which are prone to microstructural damages difficult to quantify, should be treated with caution. For this reason, we recommend whenever possible the use of fixed needles, which are set-up before the snow season, buried during snowfalls, and not subjected to microstructural damages. Yet, even with such probes underestimation of the snow thermal conductivity is still present, and can reach up to 45% for low-density snow.

In order to correct the underestimation of needle probe measurements, we proposed a simple measurement protocol, complemented with a correction lookup table. This correction table is specifically designed for the use of a fixed TP02 needle probe in snow, a model frequently employed in snow studies. Similar tables could be developed for other probe models if required. We applied this new correction protocol to measurements performed in the cold-room during a gradient box experiment, and compared the results with thermal conductivity estimations derived from 3D numerical simulations based on μCT microstructures. Application of the correction shows a net improvement of the NP values with a good agreement between the NP and μCT-based thermal conductivity estimations, when the snow was not damaged during sampling.

Fig. 1 .

 1 Fig. 1. General use of the NP method. (a) Schematics of a TP02 probe, with the two thermocouples (blue dots) and the heated wire (orange line). (b) Example of a temperature increase recorded during a NP measurement, with a non-negligible trend during the heating cycle. Note that the measurement shown here has been selected for illustration and presents one of the strongest trend observed during the gradient box experiment. Most measurements exhibit relatively small temperature trends of 0.02 K min -1 or less. (c) Estimation of the thermal conductivity from a detrended ΔT vs ln(t) curve, using the inverse of the slope of the dashed orange line, computed over a time interval centered around t mes .

Fig. 2 .

 2 Fig. 2. Schematics of the gradient box experiment with side and top views. The top view corresponds to a horizontal cut half-way up the box, and shows the position of the NP measurements and of the snow coring.

Fig. 3 .

 3 Fig. 3. Comparison of temperature increases measured in snow with a NP (in blue) and predicted by Eqn (1) with a zero contact thermal resistance (in orange).

Fig. 4 .

 4 Fig. 4. Schematics of the three types of numerical simulations performed in this work. (a) 3D homogenization simulation used to estimate the thermal conductivity of a snow sample based on its microstructure (ice phase shown in gray). (b) 3D heterogeneous NP simulation, modeling the temperature increase of a probe (in blue) taking into account heat transfer through both phases composing a snow sample (ice phase shown in gray). (c) 2D homogeneous NP simulation, modeling the temperature increase of a probe (in blue), in which the external medium is considered homogeneous.

Fig. 5 .

 5 Fig. 5. Temperature increase of a needle probe over a 300 s heating cycle, according to Eqns (1), (3) and (4). Time is given both on a linear (left) and a log (right) scale.

Fig. 7 .

 7 Fig.7. Impact of the parameter α on the error made by applying Eqn (5) to short-times NP data and for a reference thermal conductivity of 0.15 W K -1 m -1 , represented as a dashed line. The different thermal conductivity estimations are done at constant τ parameter, set to 1.5 s. Higher α values characterize snow of higher density and probes of lower volumetric thermal capacity.

Fig. 8 .

 8 Fig. 8. Extension of the simulation domain for the 3D heterogeneous NP simulation. The symmetry planes used to extend the simulation domain are shown as red lines. The probe is displayed in blue in the middle of the sample, the ice phase is shown in black and the air in light gray.

Fig. 9 .

 9 Fig. 9. Simulations of NP measurements considering snow as a heterogeneous medium (blue) or an equivalent homogeneous medium (orange). The dashed line indicates the known thermal conductivity of the snow sample.

Fig. 10 .

 10 Fig. 10. Simulations of NP measurements of snow, depending on the snow physical properties and the size of the air gap around the needle probe. The dashed lines mark the true thermal conductivity value.

Fig. 11 .

 11 Fig. 11. Needle probe measurements performed during the gradient box experiment. The thermal conductivity values are obtained directly by applying Eqn (5) to heating curve data centered around a 45 s measurement time. In blue: Measurements done with the fixed needle probe. In green: Measurements performed by manually inserting the needle probe with the help of an insertion guide. In orange: Measurements performed by manually inserting the needle probe without an insertion guide.

Fig. 12 .

 12 Fig. 12. (a) Simulations of NP measurements with and without a block of ice accreted below the needle. (b) Geometry of the accreted block of ice used for the simulations shown in light blue.

Fig. 13 .

 13 Fig. 13. Picture of recrystallized ice below a fixed needle probe after 3 weeks of temperature gradient metamorphism.

Fig. 14 .

 14 Fig. 14. Correction factor C as a function of k raw for measurement times t mes ranging from 40 to 100 s and at 268 K.

Fig. 15 .

 15 Fig. 15. (a) Snow thermal conductivity measurements performed during the gradient box experiment. In blue: uncorrected values obtained with the fixed NP and the direct application of Eqn (5). In orange: Corrected values obtained with the fixed NP and the methodology proposed in this article. In dashed black: thermal conductivity values obtained with 3D homogenization simulations, using μCT-based microstructures and under the fast surface kinetics hypothesis. In dotted black: thermal conductivity values obtained with 3D homogenization simulations, but under the slow kinetics hypothesis. (b) Density of the snow samples deduced from μCT over time. The colored zones indicate the snow type with their associated standardized symbols[START_REF] Fierz | The International Classificationi for Seasonal Snow on the Ground[END_REF], with a transition from decomposing and fragmented precipitation particles in green to depth hoar in blue. Note that we did not observe the faceted crystals transition, as the snow had already metamorphized to depth hoar for our first sampling after 3 days.

Fig. 16 .

 16 Fig. 16. Simulated and measured temperature increase of a TP02 probe during the measurement of XPS foam.

Fig. 17 .

 17 Fig. 17. Simulated and measured determination of thermal conductivity of XPS foam using a TP02 NP and Eqn (5). The dashed line marks the known thermal conductivity of XPS foam.

Downloaded from https://www.cambridge.org/core. 25 Feb 2022 at 16:48:14, subject to the Cambridge Core terms of use.

Kévin Fourteau and others Downloaded from https://www.cambridge.org/core. 25 Feb 2022 at 16:48:14, subject to the Cambridge Core terms of use.

Journal of Glaciology

Acknowledgements. This work contributes to the APT project (Acceleration of Permafrost Thaw), funded by the Climate Initiative program of the BNP-Paribas Foundation. We thank Laurent Arnaud for developing the original program used to control the NP, Yannick Deliot for his help fixing our NP electronics, and Laurent Pezard for setting up the automated control of the NP in the cold room. We acknowledge Antoine Bernard and Amaury de Nacquard for their help in the cold room. Neige Calonne and Frédéric Flin provided helpful inputs on the use of the gradient box. We are thankful to the two anonymous reviewers for reviewing our manuscript and to Argha Banerjee for editing it.

Supplementary Material. The supplementary material for this article can be found at https://doi.org/10.1017/jog.2021.127. The correction table specifically derived for the use of a TP02 probe in snow is provided as a Supplementary Material.

Author contributions.

The research was designed by F.D., K.F. and P.H. F.D. obtained funding. K.F., P.H. and J.R. performed the research. K.F. wrote the paper with inputs from F.D., P.H. and J.R. The authors declare having no competing interests.

Data availability statement. The codes used for the numerical simulations and the analysis of the heating curves will be provided upon direct request to the corresponding author.

Appendix A: Validation of the TP02 probe's physical properties

Throughout this work, we use the analytical formula of Eqn (1) and FEM simulations of NP measurements in order to quantify the degree of underestimation made with the NP method in snow. As seen in the main article, the degree of underestimation is notably governed by the physical properties of the probe. It is therefore necessary to ensure that the physical properties used to represent a TP02 probe in the article, i.e. a highly conducting material with a density of 4500 kg m -3 and a thermal capacity of 500 J kg -1 , are adequate.

To assess this validity we compared the results of a 2D FEM simulation of NP measurement with actual measurements performed in extruded polystyrene foam (XPS foam). We chose XPS foam as the test material, as it has thermal properties relatively close to light snow and which are well documented. Furthermore, XPS foam is a non-brittle material, and we can thus expect a good thermal contact with the probe. It has a density of 50 kg m -3 [START_REF] Lide | CRC handbook of Chemistry and Physics, Properties of Ice and Supercooled Water[END_REF], a thermal conductivity 0.033 W K -1 m -1 [START_REF] Lide | CRC handbook of Chemistry and Physics, Properties of Ice and Supercooled Water[END_REF], and value provided by the manufacturer), and a thermal capacity of 1270 J K -1 kg -1 (Engineering ToolBox, 2003). Finally, in order to improve the signal to noise ratio of the temperature data, the heating power Q was increased to 1.3 W m -1 , leading to a temperature increase of about 10 K after 100 s.

Results indicate a good agreement between our 2D NP simulation and the actual measurement data, both for the temperature increase of the probe shown in Figure 16, and the estimation of thermal conductivity given by application of Eqn (5) shown in Figure 17. This confirms that for NP measurements, a TP02 probe is equivalent to a highly conducting cylinder with a density of 4500 kg m -3 and a specific thermal capacity of 500 J kg -1 .


Dawin Baden

Pierre Henry

Ginette Saracco

Lionel Marié

Investigation of anisotropy in non-saturated and heterogeneous carbonate series using crosshole acoustic technics

Physical properties of carbonate rocks cannot be fully captured from laboratory-sized samples. Indeed, heterogeneous facies distribution and/or diagenetic alterations may lead to significant variations in petrophysical properties within few meters. In carbonates, diagenetic transformations are tightly related to nature of fluids flowing through the formations, e.g. via fractures network. Consequently, reservoir properties may have patchy distribution, and may not be correlatable (e.g. using facies distribution or wells-logs correlations) within few meters. Our works aim at characterizing carbonates anisotropy at different scales, and are subject of two presentations at SEG's 87 th Annual Meeting. This abstract deals with the first part of our approach, that's to say characterizing impact of diagenetic alteration on reservoir properties and seismic anisotropy, from meter to multimeter scale. Crosshole acoustic-survey has been carried out using ultrasonic frequencies (50kHz), and vertical-resolution of 10cm. The main conclusions are: (1) multi-meter crosshole surveys, enable measuring elastic properties of both matrix, and fracturedmatrix. Wave velocities measured at multi-meter scale are comparable to that from laboratory measurements, and this analogy is attributed to poor likelihood of finding dense fractures-sets -which would drastically change elastic moduli-at reduced scale. (2) Weak anisotropy is detected, varying between +10% (slow axis perpendicular to strata) and -2% (slow axis parallel to strata), and likely results from combination of matrix anisotropy (mineral orientation and pore geometry) and fractures effects. (3) In heterogeneous layered formations, variation of mean velocity with source and receivers locations significantly impact anisotropy parameters estimated by curve fitting. One should be aware of this potential bias, while interpreting crosshole surveys with sparse data, especially in carbonate formations.

Introduction:

This study is consecutive to a series of published works aiming at characterizing both reservoir and elastic properties of Lower Cretaceous limestone from the Urgonian platform of Provence (SE, France). References describing the characteristics of this platform in terms of facies, and biostratigraphy are available e.g. [START_REF] Leonide | Relation between stratigraphic architecture and multi-scale heterogeneities in carbonate platforms: The Barremian-lower Aptian of the Monts de Vaucluse[END_REF]). These microporous bioclastic limestones are analogues of different hydrocarbon reservoirs in Middle East, e.g. Kharaib and Shuaiba formations [START_REF] Borgomano | Petrophysics of Lower Cretaceous Platform Carbonate Outcrops in Provence (SE France): Implications for carbonate reservoir characterisation[END_REF]). In addition, several published databases document coupled porosity, and elastic waves velocities measurements from the Urgonian platform [START_REF] Fournier | Critical porosity and elastic properties of microporous mixed carbonate-siliciclastic rocks[END_REF][START_REF] Fournier | Elastic properties of microporous cemented grainstones[END_REF][START_REF] Fournier | Pore space evolution and elastic properties of platform carbonates (Urgonian Limestone, Barremian-Aptian[END_REF][START_REF] Borgomano | Petrophysics of Lower Cretaceous Platform Carbonate Outcrops in Provence (SE France): Implications for carbonate reservoir characterisation[END_REF]. In these works, the authors used conventional laboratory approach, i.e. velocity measurements carried out on 1" to 1.5" mini-cores, under effective stress; and they discussed "effective media theory", and the concept of "pore aspect ratio" extensively. Nevertheless, a known limitation of these approaches is the implicit assumption of medium's isotropy. Despite this assumption is generally true at the plug scale, it is rarely the case at larger scale. Carbonate rocks are often fractured and may have either "facies-related" or "patchy" porosity distribution that may lead to seismic anisotropy.

The works presented here focus on the impact of up-scaling onto elastic properties of carbonates rocks. Crosshole ultrasonic measurements have been carried out in situ, and the results are compared to published works, and to our multi-scale laboratory approach (cf. our second presentation).

Geological context:

The study took place in SE France, in an underground facility, namely LSBB (Low Noise Underground Laboratory). The LSBB is a network of sub-horizontal galleries dug in thick Barremian bioclastic limestone, within the vadose zone of "Fontaine de Vaucluse" aquifer system. The study site was located in the GAS-gallery (Figure 1), ~250m below the surface topography, which enables to set up experiments in subsurface conditions, at least in terms of hydrostatic stress, [START_REF] Guglielmi | Seismicity triggered by fluid injection-induced aseismic slip[END_REF] estimated the main stress to 6MPa. The gallery gives direct access to the 30° southward-dipping layers, and five vertical fully-cored boreholes (~20m deep) named P1 to P5 were drilled in 2009, and were studied as part of previous works [START_REF] Jeanne | Architectural, petrophysical and hydromechanical properties of fault zones in fracturedporous rocks: Compared studies of a moderate and a mature fault zones[END_REF][START_REF] Jeanne | Dissimilar properties within a carbonate-reservoir's small fault zone, and their impact on the pressurization and leakage associated with CO2 injection[END_REF]Jeanne, Guglielmi, and Cappa 2012;Jeanne et al. 2012). The geological background over the studied interval is summarized in Figure 2. The south-dipping layers intersected by the boreholes are inner-platform facies i.e. low to moderate energy, muddy or peloidal carbonate sands with abundant rudists. The cores' analysis emphasized two main facies types:

(1) calcarenites (grainstone to rudstone), and (2) calcilutites (wackestone to packestone). Since the reservoir properties are not related to facies, but instead to diagenesis [START_REF] Léonide | Diagenetic patterns and pore space distribution along a platform to outer-shelf transect (Urgonian Limestone, Barremian-Aptian[END_REF], we do not describe them any further in this abstract. The studied interval can be divided in three mechanical units: two are densely fractured (from 1 to 6m, and from 14 to 21m), and one is moderately fractured (from 6 to 14m). The formation is affected by two sets of sub-vertical fractures (dip ~80°), and oriented N30 and N120. The fractures density is controlled by both porosity and the fault zone nearness (Jeanne, Guglielmi, and Cappa 2012;[START_REF] Jeanne | Dissimilar properties within a carbonate-reservoir's small fault zone, and their impact on the pressurization and leakage associated with CO2 injection[END_REF]. The open-fractures sets are located mainly in area with gradual porosity changes close to the calcilutites. Under pressuresolution conditions, the calcilutites may have produced carbonate rich fluids (cf. stylolites) that migrated and precipitated nearby, occluding some of the fractures and decreasing locally the overall porosity of the porous layer.

Material and methods:

The acoustic logging has been carried out between P4 (Emission Borehole), and P2 (Reception Borehole) located two meters apart, along azimuth N20°. The study was conducted from 1 to 14.5m depth. The acquisition was divided into sections for technical reasons, each section is 1.5m high, and the spacing of the sensor array was 10cm. For each section 256 full-waveforms acoustic signals were recorded. The P-and Swaves travel times were manually picked and were processed, section by section, according to raypaths angles. This approach emphasized angular dependency of elastic wave velocities, and a correlation with beddings. Data were investigated using a model of anisotropic media, the transverse isotropy (TI). TI media have same properties in one plane (e.g. the x-y plane) and different properties in direction normal to this symmetry plane (e.g. z-axis). If the z-axis is the vertical axis, then the material is a Vertical Transversely Isotropic medium (VTI). For such material, the stiffness tensor (CIJ) written in Voigt's notation, requires only five independent components to fully describe its elasticity.

                     (1)
After equation ( 2), the components of CIJ can be obtained from 5 velocity measurements [START_REF] Mavko | The rock physics handbook[END_REF]. (2)

                                
The model of VTI media enables to compute phase velocities for quasi-longitudinal (Vp), quasi-shear (Vsv), and pure-shear (Vsh) modes, see equation below. (3)

                                                        
Where ρ is the density, and θ is the phase's angle. Finally, the phase velocities were converted into group velocities for direct comparison with the field data, and the stiffness components were compared to values documented in literature on Urgonian limestones.

Results and Discussion:

The velocity distribution is consistent with porosity variations within the studied area (Figure 2), and the dipping of the layers (~30°), is well figured. The highest velocities [~6.5 -6.3km/s] are between 1 and 5.5m depth, which is a low porous interval. The lowest velocities [~4.6-4.8km/s] are from 9 to 13m, where porosity reaches up to 20%. Two intervals with intermediate velocities ~5.2km/s located between ~7 and 8m and from 13 to 14.5m, show porosity variation. Although the tomography have fair consistency with the porosity, yet the anisotropy visible in data (e.g. Figure 3) seems smoothed by such the process. One can see that transverse isotropy approximation was relevant in this study, since phase velocity models are consistent with experimental data (Figure 3). Taking the example of S2, which is a high velocity interval with P-wave velocity anisotropy ~10%. The direction of maximum velocity is 30° which means it follows the stratification. One would expect having slow axis normal to beddings, but in our context, fractures-sets are normal to layers, since they are dry the reduction of elastic modulus is maximum perpendicular to the their orientation [START_REF] Hudson | Overall properties of a cracked solid[END_REF][START_REF] Hudson | Wave speeds and attenuation of elastic waves in material containing cracks[END_REF][START_REF] Rüger | P-wave reflection coefficients for transversely isotropic media with vertical and horizontal axis of symmetry[END_REF][START_REF] Guéguen | Crack-induced anisotropy in crustal rocks: Predicted dry and fluidsaturated Thomsen's parameters[END_REF][START_REF] Tsvankin | Reflection moveout and parameter estimation for horizontal transverse isotropy[END_REF]. So that, the anisotropy observed in S2 could be explained either by stress relief (and thus opening of fractures-sets), or presence of stylolithes [START_REF] Baud | Impact of stylolites on the mechanical strength of limestone[END_REF].

The second example, S8, shows despite detected anisotropy is relatively small, the value is reliable. Data from this interval showed great consistency for successive shots, and anisotropy measurements made on the matrix (cf. our 2 nd presentation at SEG's 87 th Annual Meeting), showed that the higher the intercrystalline porosity (i.e. between micrite crystals), the more the matrix was anisotropic. Moreover, this interval have few macroscopic fractures, which leaves porosity fabric and mineral orientation as the most likely explanations for this anisotropy. The particularity of S6, which is our last example, lies into the fact that there is an obvious influence of the source location on the anisotropy. This interval has increasing porosity values bottomward, plus, it is an area were fractures attenuate as they reach the porous layer. The change in the anisotropy simply show that two anisotropies (matrix and fractures) are convoluted there.

In Table 1, our results are compared to other works conducted on the same carbonates, at seismic-scale (≤180m) in [START_REF] Bereš | Seismic anisotropy analysis at the Low-Noise Underground Laboratory (LSBB) of Rustrel (France)[END_REF], at field-scale (≤13.5m) in this study, and at laboratory-scale (=2.54cm) in [START_REF] Fournier | Pore space evolution and elastic properties of platform carbonates (Urgonian Limestone, Barremian-Aptian[END_REF]. The smaller the sample the higher the coefficients variability, this is explained, by averaging effect of the scale. Indeed, different factors, i.e. facies, grains sorting, and diagenesis cause heterogeneity in carbonate rocks. Hence the smaller the samples the greater the probability of sampling particular features. Accordingly, elastic properties are averaged with the increasing scale. Interestingly, even if major fractures are present in the sections, our stiffness coefficients are closer to laboratory values than seismic-scale ones. This suggests that at meter-scale, fractures start to be sampled but do not dominate the signal, which enables to capture both signatures of the matrix and the fractures. In addition, the faults present in the LSBB area [START_REF] Jeanne | Architectural, petrophysical and hydromechanical properties of fault zones in fracturedporous rocks: Compared studies of a moderate and a mature fault zones[END_REF], and which did not affect our measurements (no fault between the tested boreholes), but which were unavoidably measured by [START_REF] Bereš | Seismic anisotropy analysis at the Low-Noise Underground Laboratory (LSBB) of Rustrel (France)[END_REF], may have played a role even bigger in lowering the velocities than fractures-sets alone. For example [START_REF] Al-Harrasi | Seismic anisotropy in a hydrocarbon field estimated from microseismic data[END_REF] related the anisotropy in the Natih carbonate reservoir (Oman) the rock type and proximity to major faults. So that, the coefficients presented in this study may be representative of the rock matrix and of the fractured rock, whereas at seismic-scale, major faults may obliterate the effects of the matrix, hence, seismic-scale coefficients may be more representative of faulted rock mass.

Conclusion:

A crosshole acoustic-survey has been carried-out using ultrasonic frequencies (50kHz), and with vertical-resolution of 10cm. This works were conducted in order to characterize, in situ, a microporous carbonate series from meter to decameter scale. The main conclusions can be summarized as follow:

(1) Crosshole ultrasonic surveys enable measuring elastic properties of both matrix, and fractured-matrix. The wave velocities measured at multi-meter scale are comparable to that from conventional laboratory measurements. This similarity can be attributed to the poor likelihood of finding multiple fractures as the dimensions of investigated area decreases.

(2) In most cases, a weak anisotropy is detected, varying between +10% (slow axis perpendicular to strata) and -2% (slow axis parallel to strata) and likely results from a combination of matrix anisotropy (mineral orientation and pore geometry) and fractures effects.

(3) In heterogeneous layered formations, variation of mean velocity with source and receivers locations significantly impact anisotropy parameter estimated by curve fitting. One should be aware of this potential bias, while interpreting crosshole surveys with sparse data, especially in carbonate formations.

Figure 1 .

 1 Figure 1. Localization GAS-gallery in the LSBB, and disposition of the boreholes.

Figure 2 .

 2 Figure2. Tomography of inter-borehole space, put in perspective with a schematic view of the geology, and the porosity (adapted fromJeanne et al. (2012). The tomography of P-wave velocity was computed using the eikonal equation(D. Baden et al. 2017, submitted).

Figure 3 .

 3 Figure3. Velocity models computed using the computed stiffness coefficients. The P-and S-wave velocities are figured by open circles and dots, respectively. The solid, the dotted and the dashed models stand for: Vp, Vsh, and Vsv, respectively. The cosine models (grey) are also plotted for comparison but are not discussed in this abstract.

Table 1 .

 1 CIJ components for Urgonian limestones (Baden et al., @submitted article).

	(×10 9 ) C11	Bereš et al., [2013] 44,9 -36,9	Baden et al., (submitted) 108,0 -41,4	Fournier et al., [2014] 110,4 -26,8
	C13	19,7 -18,3	52,5 -6,1	52,7 -9,2
	C33	69,2 -65,9	99,6 -43,4	C11
	C44	19,7 -18,3	30,3 --13,7	34 -7,9
	C66	15,4 -14,8	30,4 -12,5	C44

Acknowledgements:

This works were supported by two Research programs: H-CUBE project (ANR-12-SEED-0006), D. Baden's PhD grant; and ALBION-HPMSCa project, funded by TOTAL .

EDITED REFERENCES

Note: This reference list is a copyedited version of the reference list submitted by the author. Reference lists for the 2017 SEG Technical Program Expanded Abstracts have been copyedited so that references provided with the online metadata for each paper will achieve a high degree of linking to cited sources that appear on the Web.