Stochastic analysis of micro-cone penetration tests in snow

Cone penetration tests have long been used to characterize snowpack stratigraphy. With the development of sophisticated digital penetrometers such as the SnowMicroPen, vertical profiles of snow hardness can now be measured at a spatial resolution of a few micrometers. By using small penetrometer tips at this high vertical resolution, further details of the penetration process are resolved, leading to many more stochastic signals. An accurate interpretation of these signals regarding snow characteristics requires advanced data analysis. Here, the failure of ice connections and the pushing aside of separated snow grains during cone penetration lead to a combination of (a) diffusive noise, as in Brownian motion, and (b) jumpy noise, as proposed by previous dedicated inversion methods. The determination of the Kramers–Moyal coefficients enables differentiating between diffusive and jumpy behaviors and determining the functional resistance dependencies of these stochastic contributions. We show how different snow types can be characterized by this combination of highly resolved measurements and data analysis methods. In particular, we show that denser snow structures exhibited a more collective diffusive behavior supposedly related to the pushing aside of separated snow grains. On less dense structures with larger pore space, the measured hardness profile appeared to be characterized by stronger jump noise, which we interpret as related to breaking of a single cohesive bond. The proposed methodology provides new insights into the characterization of the snowpack stratigraphy with micro-cone penetration tests.

1 Introduction

Snow is an essential component of our environment and can significantly impact our lives: from the wishful dream of a white Christmas to the misfortune of avalanche accidents. Having a closer look at snow, one discovers many microstructural patterns and realizes that snow on the ground undergoes constant evolution (Colbeck et al., 1990). The snow microstructure fully controls its physical and mechanical properties, which are essential for diverse applications, such as avalanche forecasting (Schweizer et al., 2003). A snowpack is typically structured in numerous layers composed of different snow types, where such stratigraphy will determine the snowpack stability. Cone penetration tests have long been used to characterize the snowpack stratigraphy (Bader and Niggli, 1939). The SnowMicroPen (SMP) can perform cone penetration tests of snow in the field (Schneebeli et al., 1999). It measures the force needed to drive a cylinder with a millimetric conic tip into the snowpack. With its high resolution (250 measurements mm⁻¹), the measured force or hardness is supposedly linked to the snow microstructure (Johnson and Schneebeli, 1999). A typical consequence of such a high-precision measurement is that more and more details of the penetration process can be resolved, leading to many more stochastic signals. In this context, it is of particular interest to employ advanced data analysis to find out how different kinds of stochastic signals are related to different snow types. There are also other penetrometers used for cone penetration tests in the field of snow; for example McCallum (2013) used a large-scale penetrometer with a tip diameter of 36.7 mm in polar snow. For our analysis, we specifically used the measurement data of micro-cone penetration tests from Schneebeli et al. (1999).

Johnson and Schneebeli (1999) developed the first model to estimate micromechanical properties of snow from measured penetration profiles. They assumed that the material compaction is negligible and that the penetration resistance is only composed of friction between the penetrometer and snow grains and a superposition of spatially uncorrelated and identical brittle failures of individual snow microstructural elements (e.g., the bonds between the snow grains). Marshall and Johnson (2009) extended the theory of Johnson and Schneebeli (1999) to account for simultaneous ruptures by Monte Carlo simulations. They precisely resolved the snow micromechanical parameters, such as the deflection length, rupture force and rupture intensity. Löwe and van Herwijnen (2012) re-stated the pioneering idea of Johnson and Schneebeli (1999) and described the fluctuating penetration hardness as a Poisson shot noise process. In their model, the micromechanical parameters can be simply derived from the cumulants and the co-variance of the penetration signal. Peinke et al. (2019) further extended the homogeneous Poisson process of Löwe and van Herwijnen (2012) so that the scale of variation in the rupture intensity could be decoupled from the scale of variations in the deflection length and the rupture force. These models are now commonly used to characterize the snowpack stratigraphy from SMP measurements. Indeed, Proksch et al. (2015) related the micromechanical parameters derived from SMP measurements to some of the most critical snow characteristics, namely density, specific surface area and structural correlation length. These relations are now routinely used to quantify the snowpack stratigraphy (e.g., Calonne et al., 2020). Besides, Reuter et al. (2019) estimated the elastic modulus and fracture energy from the micromechanical parameters, which can then be used to compute point snow stability for avalanche hazard assessment (e.g., Reuter et al., 2015).

Here, we consider the measured fluctuating hardness as a consequence of summing up the interactions between the penetrometer tip and individual snow particles. We describe this penetration process in analogy to the well-known Brownian motion (Einstein, 1905), where a microscopically visible particle suspended in fluid is moving randomly. Due to the sum of several collisions with the molecules in the fluid as illustrated in Fig. 1, the large red particle undergoes a motion described by a stochastic process. Such a stochastic process is driven by white noise and is known as a Langevin process (see Gardiner, 1985, and Sect. 2). While very similar elementary collision events are assumed for classical Brownian motion, we also need to consider brittle failures of individual snow microstructural elements (bonds between snow grains, crushing of grain clusters), which cause sharp declines in penetration hardness. These brittle failures were modeled as a Poisson shot noise process by Löwe and van Herwijnen (2012) and Peinke et al. (2019). Jump noise acts as discontinuous paths inside the diffusion process, and low-jump events can be considered Poisson distributed noise (which corresponds to the abovementioned shot noise). The idea of this work is to employ a method that allows estimating the underlying stochastic differential equation from empirical data and differentiating between a Langevin (pure diffusive) and a jump-diffusion process (Anvari et al., 2016). Via this advanced analysis, we seek more detailed snow characterization from micro-cone penetration test resistance data.¹

Figure 1 (a) Penetration resistance caused by the interactions of the snow particles (gray) with the penetrometer tip (red); (b) Brownian motion where a microscopically visible particle (red) suspended in the fluid is moving randomly due to the collisions with the molecules (gray) in the fluid.

The article is organized as follows. In Sect. 2, we summarize the stochastic analysis method and show how it is possible to distinguish between diffusive and jump noise. In Sect. 3, the method is applied to centimetric snow samples whose microstructure is also captured by tomography. In Sect. 4, as a proof of concept, the SMP profile of a natural snowpack is analyzed with this technique.

2 Stochastic method

A stochastic process x(t) can be described through stochastic differential equations. This section explains the equations used to model cone penetration tests in snow. Since the SMP is driven by a motor with constant speed $u=\frac{\mathrm{d}z}{\mathrm{d}t}$ (z is depth; t is time) and samples the measurement every 4 µm, the measured penetration force or snow hardness R is considered the depth dynamics R(z(t)) and handled like a stochastic variable x(t).

2.1 Langevin equation

A diffusive process x(t), which is a continuous stochastic process, follows the Langevin equation, where for a small step size dt has the following expression (Risken, 1984; Friedrich et al., 2011; Tabar, 2019):

$$\begin{array}{}\text{(1)}& \mathrm{d}x\left(t\right)=D^{\left(\mathrm{1}\right)}(x,t)\mathrm{d}t+\sqrt{D^{\left(\mathrm{2}\right)}(x,t)}\mathrm{d}{W}_t,\end{array}$$

where D⁽¹⁾(x,t) and D⁽²⁾(x,t) are the drift and the diffusion coefficients, respectively, and W_t is the Wiener process. The drift term D⁽¹⁾ describes how fluctuations relax to the local mean values of x, defined by $D^{\left(\mathrm{1}\right)}(x,t)=\mathrm{0}$. The diffusion term D⁽²⁾ represents the amplitude of the noise. The coefficients D⁽¹⁾ and D⁽²⁾ are also known as the first- and second-order Kramers–Moyal (KM) coefficients, respectively. In general, KM coefficients can be directly determined from the given data x(t) using their definitions of the conditional incremental average (Friedrich et al., 2011; Tabar, 2019); i.e.,

$$\begin{array}{}\text{(2)}& K^{\left(j\right)}(x,t)=\underset{\mathrm{\Delta }t\to \mathrm{0}}{lim}{\displaystyle \frac{〈{\left(x(t+\mathrm{\Delta }t)-x\left(t\right)\right)}^j{\mathrm{|}}_{x\left(t\right)=x}〉}{\mathrm{\Delta }t}}.\end{array}$$

Further details on methods of this estimation can be found in Friedrich et al. (2011) and Rinn et al. (2016)². The Langevin equation describes a continuous diffusion process where $K^{\left(j\right)}(x,t)=\mathrm{0}$ for j≥3 and $D^{\left(j\right)}(x,t)=K^{\left(j\right)}(x,t)$. According to Pawula's theorem, all KM coefficients K^(j) are vanishing for j≥3 if $K^{\left(\mathrm{4}\right)}(x,t)=\mathrm{0}$ (Risken, 1984; Pawula, 1967). In our case, however, the higher-order coefficients are not always vanishing; hence we extend the discussion to the jump-diffusion process (see Sect. 2.2).

For the SMP data considered, the drift term D⁽¹⁾ in Eq. (1) describes how hardness fluctuations relax to the local mean values of hardness R, defined by $D^{\left(\mathrm{1}\right)}(R,z)=\mathrm{0}$. The diffusion D⁽²⁾ term represents the amplitude of the hardness fluctuations. The coefficients D⁽¹⁾ and D⁽²⁾ are z-dependent for non-stationary (inhomogeneous) processes. Here, we assume that for a chosen small depth interval (z), D⁽¹⁾(R,z) and D⁽²⁾(R,z) only depend on R (similarly to Löwe and van Herwijnen, 2012, for their shot noise model).

2.2 Jump-diffusion dynamics

Typically, when the signal of a stochastic process presents sharp changes at some instant (discontinuity events), higher-order Kramers–Moyal coefficients (especially K⁽⁴⁾(x,t)) become non-negligible. In this case, an extension to Langevin-type modeling with additional jump noise is needed (see Tankov, 2003; Stanton, 1997; Johannes, 2004; Bandi and Nguyen, 2003; Anvari et al., 2016; Tabar, 2019). Such a jump-diffusion dynamic is given by the following stochastic differential equation:

$$\begin{array}{}\text{(3)}& \mathrm{d}x\left(t\right)=D^{\left(\mathrm{1}\right)}(x,t)\mathrm{d}t+\sqrt{D^{\left(\mathrm{2}\right)}(x,t)}\mathrm{d}{W}_t+\mathit{\xi }\mathrm{d}{J}_t,\end{array}$$

where, again, D⁽¹⁾ and D⁽²⁾ are the drift and the diffusion coefficients, respectively, and W_t is the Wiener process. The quantity ξ is the size of the jump noise and is assumed to be normally distributed; i.e., $\mathit{\xi }\sim N(\mathrm{0},{\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}})$, where ${\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}(x,t)$ is the so-called jump amplitude. The term J_t is the Poisson jump process, which is the zero–one jump process with a jump rate (or intensity) λ(x,t) (Hanson, 2007; Tabar, 2019).

For jump-diffusion processes, the drift and diffusion coefficients (D⁽¹⁾, D⁽²⁾), the jump rate ${\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}$ and amplitude λ are related to the KM coefficients as (Anvari et al., 2016):

$$\begin{array}{}\text{(4)}& {\displaystyle }{D}^{\left(\mathrm{1}\right)}(x,t)=K^{\left(\mathrm{1}\right)}(x,t),\text{(5)}& {\displaystyle }{D}^{\left(\mathrm{2}\right)}(x,t)+\mathit{\lambda }(x,t)\langle {\mathit{\xi }}^{\mathrm{2}}\rangle =K^{\left(\mathrm{2}\right)}(x,t),\text{(6)}& {\displaystyle }\mathit{\lambda }(x,t)\langle {\mathit{\xi }}^j\rangle =K^{\left(j\right)}(x,t),\text{for}j>\mathrm{2}.\end{array}$$

The estimate of the drift coefficient is the same for the diffusion process (Eq. 1) and the jump-diffusion process (Eq. 4). Jump amplitude ${\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}$ and rate λ can be estimated by using Eq. (6) with j=4 and j=6 and Wick's theorem (Isserlis, 1916; Wick, 1950) for Gaussian random variables, which states that $\langle {\mathit{\xi }}^{\mathrm{2}n}\rangle =\frac{\left(\mathrm{2}n\right)\mathrm{!}}{{\mathrm{2}}^{n}n\mathrm{!}}\langle {\mathit{\xi }}^{\mathrm{2}}{\rangle }^n$:

$$\begin{array}{}\text{(7)}& {\displaystyle }{\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}(x,t)={\displaystyle \frac{K^{\left(\mathrm{6}\right)}(x,t)}{\mathrm{5}{K}^{\left(\mathrm{4}\right)}(x,t)}},\text{(8)}& {\displaystyle }\mathit{\lambda }(x,t)={\displaystyle \frac{K^{\left(\mathrm{4}\right)}(x,t)}{\mathrm{3}{\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{4}}(x,t)}}.\end{array}$$

To improve the estimation of KM coefficients K^(j)(x,t) and in particular of high-order coefficients, the Nadaraya–Watson estimator, which is a kernel estimator, can be used (Nadaraya, 1964; Watson, 1964):

$$\begin{array}{}\text{(9)}& K^{\left(j\right)}(x,t)=\underset{\mathrm{\Delta }t\to \mathrm{0}}{lim}{\displaystyle \frac{{\sum }_{i}k\left(\frac{x_{i\mathrm{\Delta }t}-x}{h}\right)(x_{(i+\mathrm{1})\mathrm{\Delta }t}-x_{i\mathrm{\Delta }t}{)}^j}{{\sum }_{i}k\left(\frac{x_{i\mathrm{\Delta }t}-x}{h}\right)\mathrm{\Delta }t}},\end{array}$$

where here we use a Gaussian kernel k(u) (Tabar, 2019). With the kernel-based method the conditional moments can be calculated more smoothly by controlling the kernel bandwidth h (Lamouroux and Lehnertz, 2009), where here we use the kernel bandwidth h=0.3.

For the SMP measurements considered, we also assume that for a chosen small depth interval (z), D⁽¹⁾(R,z), D⁽²⁾(R,z), ${\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}(R,z)$ and λ(R,z) only depend on R. The stochastic differential equation for the jump-diffusion process thus reads, at small depth intervals, as

$$\begin{array}{}\text{(10)}& \mathrm{d}R\left(z\right)=D^{\left(\mathrm{1}\right)}\left(R\right)\mathrm{d}z+\sqrt{D^{\left(\mathrm{2}\right)}\left(R\right)}\mathrm{d}{W}_z+\mathit{\xi }\mathrm{d}{J}_z,\end{array}$$

with an interpretation of the drift and diffusion terms (D⁽¹⁾, D⁽²⁾) analogous to the purely diffusive process (Eq. 1) but now extended by a jump noise term. This is the same as Eq. (3) in terms of depth z instead of time t. The jump rate λ has the dimension of $\frac{\mathrm{1}}{\left[Z\right]}$ and can be related to the shot noise intensity described by Löwe and van Herwijnen (2012) and Peinke et al. (2019). Typically, λ dz corresponds to the stochastic average of the number of jumps for a penetration increment of depth dz (Anvari et al., 2016; Tabar, 2019). The jump amplitude ${\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}$ represents the square of the typical size of a jump. Note that the jump can be negative (failure of a microstructural element) or positive (loading of a microstructural element). Here, we do not consider any progressive loading of a microstructural element as described by Löwe and van Herwijnen (2012) and Peinke et al. (2019) with the microstructural deflection length δ. Here, the loading of a microstructural element and its contribution to penetration hardness are considered instantaneous.

2.3 Synthetic examples

In this section, we illustrate how diffusive and jump noises affect the stochastic fluctuations on a synthetic example. An Ornstein–Uhlenbeck (OU) process x is described by a stochastic differential equation (SDE) with linear relaxation dynamics and an additive uncorrelated noise. With an additional jump term, it is defined as

$$\begin{array}{}\text{(11)}& \mathrm{d}x=-\mathit{\gamma }x\mathrm{d}t+\sqrt{D}\mathrm{d}{W}_t+\mathit{\xi }\mathrm{d}{J}_t.\end{array}$$

where γ is the relaxation rate, D the constant diffusion coefficient and W_t a scalar Wiener process. The noise $\mathit{\xi }\sim N(\mathrm{0},{\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}})$ is assumed to be normally distributed with the constant variance or jump amplitude ${\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}$. J_t∼P(λt) is the Poisson jump process, which is a zero–one jump process with constant jump rate λ. Here, we have triply stochastic processes W_t, J_t and ξ, which are assumed to be independent of each other.

Three synthetic time series of the OU jump-diffusion process were generated for $\mathrm{\Delta }t={\mathrm{10}}^{-\mathrm{3}}\mathrm{s}$ with $\mathit{\gamma }=\mathrm{100}{\mathrm{s}}^{-\mathrm{1}}$, $D=\mathrm{10}{\mathrm{s}}^{-\mathrm{1}}$ and for the additional jump terms with $\mathit{\lambda }=\mathrm{100}{\mathrm{s}}^{-\mathrm{1}}$ and ${\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}=\mathrm{1}$. The generated data were normalized with their respective standard deviation. In Fig. 2 the normalized time series of the OU jump-diffusion processes are shown. A pure diffusion process (left), a drift-jump process (middle) and the combined jump-diffusion process (right) are shown. The fourth-order KM coefficients K⁽⁴⁾(x) of each time series are also plotted in Fig. 2, bottom row. K⁽⁴⁾ of the diffusion process is negligibly small compared to drift-jump and jump-diffusion processes.

Figure 2 Normalized time series of Ornstein–Uhlenbeck processes with only diffusion, only jump (drift-jump) and jump-diffusion terms (a, b, c) and their fourth-order KM coefficients K⁽⁴⁾ (d, e, f). K⁽⁴⁾ of the diffusion process is negligible compared to drift-jump and jump-diffusion processes. For all three examples, we used the same noises in the stochastic part of the stochastic differential equation.

For a jump-diffusion process, another parameter that we considered was the ratio of diffusion and jump noise $\frac{D^{\left(\mathrm{2}\right)}}{\mathit{\lambda }{\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}}$, which becomes here for our OU jump-diffusion process $\frac{D}{\mathit{\lambda }{\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}}$. To validate our method, based on the KM coefficients of Eq. (5) to Eq. (8), three pairs of parameters were chosen: (i) $D=\mathrm{5}{\mathrm{s}}^{-\mathrm{1}}$, ${\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}=\mathrm{1}$; (ii) $D=\mathrm{5}{\mathrm{s}}^{-\mathrm{1}}$, ${\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}=\mathrm{0.5}$; and (iii) $D=\mathrm{20}{\mathrm{s}}^{-\mathrm{1}}$, ${\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}=\mathrm{1}$. The other parameters are the same as in the previous example. The normalized time series of these examples are plotted in Fig. 3, top row. In Fig. 3, bottom row, the ratios of diffusion and jump noise $\frac{D}{\mathit{\lambda }{\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}}$ estimated from each time series were compared to the expected values (blue lines). As we used normalization and the same noises in simulation, all time series are very similar; however, one can observe clearly the much noisier fine structure in case (iii) where the diffusion coefficient is larger.

Figure 3 Normalized time series of OU jump-diffusion processes with $\mathrm{\Delta }t={\mathrm{10}}^{-\mathrm{3}}\mathrm{s}$, $\mathit{\gamma }=\mathrm{100}{\mathrm{s}}^{-\mathrm{1}}$ and $\mathit{\lambda }=\mathrm{100}{\mathrm{s}}^{-\mathrm{1}}$ for (i) $D=\mathrm{5}{\mathrm{s}}^{-\mathrm{1}}$, ${\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}=\mathrm{1}$; (ii) $D=\mathrm{5}{\mathrm{s}}^{-\mathrm{1}}$, ${\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}=\mathrm{0.5}$; and (iii) $D=\mathrm{20}{\mathrm{s}}^{-\mathrm{1}}$, ${\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}=\mathrm{1}$ (a, b, c) and the corresponding ratio of diffusion and jump noise $\frac{D^{\left(\mathrm{2}\right)}}{\mathit{\lambda }{\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}}=\frac{D}{\mathit{\lambda }{\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}}$ (d, e, f). Dots are results from the KM coefficients, and the blue line denotes the theoretical values given by the constants. For all three examples, we used the same noises in the stochastic part of the SDE.

3 Application to snow measurements

In this section, our main aim is to show how the jump-diffusion model can be used to distinguish snow types from hardness profiles measured with the SMP. Firstly, small snow samples whose microstructure was also fully characterized by tomography before being measured by the SMP were used to test the developed methodology. Secondly, as a proof of concept, we analyzed one penetration profile of a snowpack measured in the field and we provided the subsequent profile of microstructural parameters. Last, the results were discussed.

Figure 4 Three-dimensional view of the microstructure of some representative samples: precipitation particles (PP1), depth hoar (DH1), rounded grains (RG1) and large rounded grains (RGlr1). The ice matrix is shown in gray; the pore space is transparent. Sub-samples shown are cubic with a side length of 3 mm. Details on the data acquisition can be found in Peinke et al. (2020).

3.1 Laboratory samples

3.1.1 Measurement data

We tested several snow samples composed of four different natural snow types, namely precipitation particles (PP), depth hoar (DH), rounded grains (RG) and large rounded grains (RGlr) as classified in Fierz et al. (2009). The samples were prepared by sieving snow into small sample holders (diameter and height of 20 mm) and letting them sinter for a couple of days at −10 ^∘C. Their microstructure was captured with X-ray tomography at a nominal resolution of 15 µm (Fig. 4). The cone penetration test was conducted with a modified version of the SMP, as shown in Fig. 5. More information on sample preparation and the SMP measurement can be found in the study of Peinke et al. (2020). The main sample properties are summarized in Table 1, and the measured hardness profiles (one example for each snow type) are plotted in Fig. 6. The first 4 mm is affected by the progressive penetration of the conic tip and was not considered in the stochastic analysis (Peinke et al., 2019). The remaining profiles were divided into smaller segments of a depth of 10 mm.

Table 1 Overview of the detailed properties of the snow samples used. Snow types are classified according to the international classification of snow on the ground (Fierz et al., 2009). The density and specific surface area (SSA) are derived from the tomographic images from Peinke et al. (2020). Additionally, the standard deviations σ_R of the detrended profiles are calculated.

Figure 5 Setup of micro-cone penetration test (Peinke et al., 2020). The cone penetration tests (CPTs) were conducted by inserting a cylinder of diameter 5 mm, with a conical tip of an apex angle of 60^∘, into the snow samples. The samples were placed in the cylinder sample holder with a diameter of 20 mm. The cone was inserted vertically at a constant speed of 20 mm s⁻¹. The SMP force sensor (Kistler 9207) measures forces in the range of [0, 40] N with a resolution of 0.01 N.

Figure 6 Segments of snow hardness profiles of PP1, DH1, RG1 and RGlr1. These four different types occur as natural snow types. Precipitation particles (PP) have the smallest trend and fluctuation force; large rounded grains (RGlr) are the largest, while depth hoar (DH) and rounded grains (RG) have similar trends and fluctuation forces between those of PP and RGlr. The first 4 mm is affected by the progressive penetration of the conic tip and was not considered in the stochastic analysis.

To work out the significance of advanced stochastic features for snow, we focused on the fluctuations in the hardness profiles. Each profile was first detrended. The trend $\stackrel{\mathrm{‾}}{R}$ was computed as the convolution of the original signal with a Gaussian kernel with a standard deviation of 0.6 mm.³ The fluctuation amplitude σ_R was computed as the standard deviation of $R-\stackrel{\mathrm{‾}}{R}$. The detrended profiles are defined as $R^{\prime }=\frac{R-\stackrel{\mathrm{‾}}{R}}{{\mathit{\sigma }}_R}$. The detrended profiles R^′ are characterized by a zero mean and a standard deviation of 1. The average value of $\stackrel{\mathrm{‾}}{R}$ on the segment and the value of the standard deviation σ_R are shown for each segment in Table 1. The detrended profiles R^′ are shown in Fig. 7, for all four snow types.

To estimate errors, we divided the detrended and normalized data into different sub-samples. Given are two PP, three DH, three RG and six RGlr measurement profiles. These profiles were separated into smaller segments, which finally gives four PP, five DH, five RG and six RGlr samples. We estimated the KM coefficients of each sample using Eq. (2) and averaged them over the sub-samples of each snow type. Thus, drift, diffusion functions and jump parameters and their errors were estimated. The errors were reported as the standard error of the means.

Figure 7 Detrended snow hardness profiles of four different snow types, PP, DH, RG and RGlr. For the snow type RGlr, we only show the samples RGlr4, RGlr5 and RGlr6, as described in Table 1. The detrended profiles R^′(z) of each snow type are shifted vertically for better visualization.

3.1.2 Results

According to the description provided in Sect. 2, the drift D⁽¹⁾(R^′) and diffusion D⁽²⁾(R^′), as well as the fourth-order KM coefficients K⁽⁴⁾(R^′), for normalized data were determined for the four different snow types, PP, DH, RG and RGlr as shown in Fig. 8. In addition, the autocorrelation function (ACF) was determined from the signals. If the fluctuations in R^′(z) belong to the diffusion processes, one would expect that $K^{\left(\mathrm{4}\right)}\left(R^{\prime }\right)=\mathrm{0}$. However, we find that in general, $K^{\left(\mathrm{4}\right)}\left(R^{\prime }\right)\ne \mathrm{0}$ and the higher-order KM coefficients are not negligible. This indicates the presence of discontinuities in the snow hardness profile, so the jump-diffusion model is considered. Therefore, we estimated jump parameters, i.e., jump amplitudes ${\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}\left(R^{\prime }\right)$ and jump rates λ(R^′) from data of R^′(z) in Fig. 9.

Figure 8 State-dependent drift D⁽¹⁾(R^′), diffusion D⁽²⁾(R^′), fourth-order KM coefficients K⁽⁴⁾(R^′) and their respective autocorrelation functions (ACFs) of four different snow types, PP, DH, RG and RGlr (left to right). The errors are shown as gray-shaded background. The red lines in the ACF plots indicate the correlation length scales determined from ACFs. Comparing the correlation length scales $L_{\mathrm{C}}=\frac{\mathrm{1}}{\mathit{\gamma }}$ where $D^{\left(\mathrm{1}\right)}=-\mathit{\gamma }{R}^{\prime }$ with those of the autocorrelation functions (ACFs), we find that both length scales have the same ordering of their values for all snow types.

As shown in Fig. 8, the drift coefficients, $D^{\left(\mathrm{1}\right)}(R^{\prime },z)$, are mostly linear functions with negative slopes, which describe how fast the system tends back toward the stable fixed point. Due to our normalization, the fixed point of dynamics is located at the origin, i.e., $R^{\prime }=\mathrm{0}$. Taking $D^{\left(\mathrm{1}\right)}=-\mathit{\gamma }{R}^{\prime }$, the correlation length scale is given by $L_{\mathrm{C}}=\frac{\mathrm{1}}{\mathit{\gamma }}$. For each snow type, L_C was determined for $-\mathrm{2}<R^{\prime }<\mathrm{2}$: $\mathit{\left\{}\text{PP, DH, RG, RGlr}\mathit{\right\}}=\mathit{\{}\mathrm{0.01},\mathrm{0.04},\mathrm{0.02},\mathrm{0.08}\mathit{\}}\mathrm{mm}$. Comparing L_C with the correlation lengths of the autocorrelation functions – L_ACF, $\mathit{\left\{}\text{PP, DH, RG, RGlr}\mathit{\right\}}=\mathit{\{}\mathrm{0.006},\mathrm{0.025},\mathrm{0.016},\mathrm{0.038}\mathit{\}}\mathrm{mm}$ – we find that both length scales have the same ordering of their values for all snow types. The snow types PP and RG have the shorter correlation length scale, in comparison to the snow types DH and RGlr.

Figure 9 Jump amplitudes ${\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}\left(R^{\prime }\right)$, jump probabilities λ(R^′) Δz, and diffusion and jump ratios $\frac{D^{\left(\mathrm{2}\right)}}{\mathit{\lambda }{\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}}$ of four different snow types, PP, DH, RG and RGlr (left to right). More data are present in the range of $-\mathrm{2}<R^{\prime }<\mathrm{2}$; we focus our statistical analysis in this range with fewer uncertainties. The errors are shown as gray-shaded background. The ratios of diffusion and jump noise $\frac{D^{\left(\mathrm{2}\right)}}{\mathit{\lambda }{\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}}$ for PP and RG are minimum near zero and maximum for DH and RGlr, which means that the larger the ice structure, the stronger the diffusion noise, and vice versa. The horizontal blue lines show the mean values of the respective parameters in the range of $-\mathrm{2}<R^{\prime }<\mathrm{2}$.

The jump amplitudes, ${\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}\left(R^{\prime }\right)$, and the jump probabilities, λ(R^′) Δz, are shown in Fig. 9. The jump amplitude, ${\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}\left(R^{\prime }\right)$, indicates how large the jump noise for different R^′ values is. The jump probability describes how probable jumps or discontinuities in forces can occur. To analyze whether diffusion or jump noise is dominating, the dimensionless ratio of diffusion and jump noise $\frac{D^{\left(\mathrm{2}\right)}}{\mathit{\lambda }{\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}}$ (Fig. 9, bottom row) was calculated. For a rough estimation, the mean values were determined in the range of $-\mathrm{2}<R^{\prime }<\mathrm{2}$ and are plotted as horizontal blue lines. The mean of $\stackrel{\mathrm{‾}}{\mathit{\lambda }}$ can be used to define the second characteristic length scale $L_{\mathrm{J}}=\frac{\mathrm{1}}{\mathit{\lambda }}$ (apart from $\frac{\mathrm{1}}{\mathit{\gamma }}$). For the abovementioned range of R^′ we obtained $\mathit{\left\{}\text{PP, DH, RG, RGlr}\mathit{\right\}}=\mathit{\{}\mathrm{0.006},\mathrm{0.01},\mathrm{0.007},\mathrm{0.02}\mathit{\}}\mathrm{mm}$. Results are summarized in Table 2, and we discuss these in Sect. 4.

Table 2 Summary of the results for the correlation length scales $L_{\mathrm{C}}=\frac{\mathrm{1}}{\mathit{\gamma }}$, $L_{\mathrm{J}}=\frac{\mathrm{1}}{\mathit{\lambda }}$, jump amplitude ${\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}$, and diffusion and jump ratio $\frac{D^{\left(\mathrm{2}\right)}}{\mathit{\lambda }{\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}}$ for all snow types analyzed. The results were evaluated in the range of $-\mathrm{2}<R^{\prime }<\mathrm{2}$.

3.2 Application to field snow data

3.2.1 Measurement data

Next, measurements from a field campaign are presented (Hagenmuller and Pilloix, 2016). The measurements were also performed with an SMP, but the tip had a slightly different shape corresponding to the standard version of the SMP (Johnson and Schneebeli, 1999). The spatial sampling is again 4 µm. This difference in the measurement methods was irrelevant, as we subsequently show with these preliminary results that in principle, the stochastic methodology can also be applied to real snow data and that qualitatively comparable results are obtained.

3.2.2 Results

The snow hardness profile of a field measurement is shown in the top left of Fig. 10. The measurement profile is strongly inhomogeneous; therefore, we used the Nadaraya–Watson estimator to determine the local characteristics of the profile. Using the moving-window technique, the profile was separated into non-overlapping windows of 500 data points (2 mm), and the detrending was performed on each window by the Gaussian kernel with a kernel size of 0.6 mm, normalized with its standard deviation as in our previous analysis of laboratory data. For each depth value z and the corresponding value R^′(z), the local values of the fourth-order KM coefficient K⁽⁴⁾(z) and the jump amplitude ${\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}\left(z\right)$ can be determined, as shown in Fig. 10. The local characteristic of each snow type from the previous section is plotted in the right column of Fig. 10 for better comparison with the field measurement data. Interpretation of these results will be discussed next.

Figure 10 The snow hardness profile of a field measurement together with its fourth-order KM coefficients K⁽⁴⁾ and the jump amplitude ${\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}\left(z\right)$ is shown in (a), (c) and (e). These parameters were determined using non-overlapping moving window with 500 data points (2 mm) by means of the Nadaraya–Watson estimator. Snow hardness profiles of the laboratory-prepared snow types and their local parameters are also plotted to enable better comparison (b, d, f); they are shifted horizontally for better visualization. The results of field measurements are shown at the depth of $\mathrm{700}\mathrm{mm}<z<\mathrm{740}\mathrm{mm}$. With reference to the local characteristic snow types from laboratory measurements, we can see in the dynamics that mixtures of different snow types are present in this section of measurement. For $\mathrm{732}\mathrm{mm}<z<\mathrm{740}\mathrm{mm}$, the high K⁽⁴⁾ and ${\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}$ indicate the presence of small and less dense structures of snow which resemble the RG-like snow types.

4 Discussion

Our work is based on the proposed analogy of Brownian motion and the SMP penetration process, as illustrated in Fig. 1. The events of bond breaking or of collision with molecules are summed up in a mean force and noise. For continuous Brownian noise, we need an integration over sufficient micro-scale events, as discussed by Einstein (1905) in his original paper. In our interpretation, it is found that sufficiently large particles lead to this integration; see Fig. 1b. To our interpretation, this integration over discrete single events of bond breaking in the immediate surroundings of the SMP, in addition to the pushing aside of loose snow grains during the penetration process, forms continuous Brownian noise. The jump noise may represent the bond-breaking events occurring directly at the tip of the SMP, and the amplitude of the jump noise should depend on the strength of the ice bonds and void sizes. From this interpretation, it is clear that snow type morphology, shown in Fig. 4, is essential for effective stochastic analysis as outlined herein.

We start the discussion with the mean values of R and the standard deviation σ_R (see Fig. 6 and Table 1). The less dense PP and the dense RGlr snow can be well separated, whereas the differences are less clear for DH and RG. In the following, we discuss the measured SMP penetration profiles based on our stochastic results. Since we now focus on a stochastic investigation of the fluctuations in the penetration profiles, the detrended and normalized data are R^′ used. Furthermore, we can note that the normalization of the snow profiles affects neither the correlation length scales $L_{\mathrm{C}}=\frac{\mathrm{1}}{\mathit{\gamma }}$ from the drift coefficients nor the jump characteristic length scales $L_{\mathrm{J}}=\frac{\mathrm{1}}{\mathit{\lambda }}$. Because this analysis depends on the number of available data, our discussion of the estimated KM coefficients is limited to the range $-\mathrm{2}<R^{\prime }<\mathrm{2}$.

The drift terms D⁽¹⁾, Fig. 8, are all monotonously decaying with increasing R^′ and can be approximated by a linear decay, $D^{\left(\mathrm{1}\right)}=-\mathit{\gamma }{R}^{\prime }$. The slope indicates how fast the signals relax to a fixed point located at $R^{\prime }=\mathrm{0}$. The magnitudes of the slopes are PP > RG > DH > RGlr; thus the PP snow has the fastest relaxation or shortest correlation length scale L_C. We find that the larger the ice structures, the slower (or longer) the relaxation. If we compare this result with the snow structures shown in Fig. 4, we conclude that L_C or γ is clearly related to the size of the snow structures.

The results for D⁽²⁾ show that about the same diffusive noise amplitude is found for all snow types. In contrast, we see clear differences for the fourth-order KM coefficients K⁽⁴⁾. Although K⁽⁴⁾≠0 is always the case, clear differences in the magnitude of this KM coefficient are found. K⁽⁴⁾ is the largest in PP, followed by RG, DH and RGlr.

The amplitudes of jump noise, ${\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}$, are the highest for PP, followed by RG, DH and RGlr. For the jump probabilities λ Δz, we distinguish a group composed of PP and RG, with higher values, and one composed of DH and RGlr, with lower values. One can interpret this finding such that for the precipitation particles (PP, recent snow) with very small ice structures and high porosity, the breaking occurs easily and frequently, which explains that the jump probability has the largest contribution here. Similarly, RG is also less dense with smaller ice structures than DH and RGlr. Thus, we can also interpret this finding such that the smaller the ice structures and the less densely the snow is packed, the stronger the jump noise. For the densely packed snow with larger ice structures, the breaking of the ice structures is less frequent, which explains a lower jump probability. The mean of the jump rate $\stackrel{\mathrm{‾}}{\mathit{\lambda }}$ in the range of $-\mathrm{2}<R^{\prime }<\mathrm{2}$ can be used to define the second characteristic length scale $L_{\mathrm{J}}=\frac{\mathrm{1}}{\mathit{\lambda }}$. Similarly to the correlation length scale $L_{\mathrm{C}}=\frac{\mathrm{1}}{\mathit{\gamma }}$, PP has the smallest length, followed by RG, DH and RGlr.

Besides the features of the different terms in the stochastic processes, the contributions of the diffusive and the jump noise can be compared by the dimensionless quotient of $\frac{D^{\left(\mathrm{2}\right)}}{\mathit{\lambda }{\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}}$, i.e., the relation between the two noise contributions. Consistent with our earlier discussion, the smallest values for $\frac{D^{\left(\mathrm{2}\right)}}{\mathit{\lambda }{\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}}$ are obtained for PP; i.e., the jump noise dominates due to the frequent fracture of small (soft) ice structures. For the other snow types, we see that within the range of $-\mathrm{2}<R^{\prime }<\mathrm{2}$, the values of $\frac{D^{\left(\mathrm{2}\right)}}{\mathit{\lambda }{\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}}$ increase with larger ice structures, in accordance with Fig. 4. The quotient $\frac{D^{\left(\mathrm{2}\right)}}{\mathit{\lambda }{\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}}$ is smaller for PP and RG and larger for DH and RGlr. For RGlr, the diffusive noise dominates in a broad range of R^′ values. For the densely packed snow with larger grain sizes, it also takes much force to push the ice grains on the side but not necessarily to break the cohesive bonds close to the tip. Therefore the penetration signal is dominated by Brownian noise. This result is also consistent with the value of K⁽⁴⁾, which is relatively small for RGlr. It is interesting to see that $\frac{D^{\left(\mathrm{2}\right)}}{\mathit{\lambda }{\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}}$ and K⁽⁴⁾ enable differentiation between RG and DH. In contrast, according to the classical statistical features of the snow signals shown in Table 1, the differences for DH and RG are less clear.

In Sect. 3.2, we analyze field measurement data which are highly inhomogeneous. With reference to the laboratory measurements, we see dynamics that suggest mixtures of different snow types within this depth segment. For $\mathrm{732}\mathrm{mm}<z<\mathrm{740}\mathrm{mm}$, the high K⁽⁴⁾ and ${\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}$ indicate the presence of small and less dense structures of snow which resemble the RG-like snow types. Based on these preliminary results of real field data, the developed methodology appears promising for interpreting cone penetration tests in the field, but further quantitative evaluation is required.

5 Conclusions

In conclusion, we observe that the advanced stochastic analysis of SMP measurements of snow layers allows differentiation of snow types. The diffusive and jump-noise contribution can be quantified and gives new insights into the stochastic behaviors of the cone penetration test in snow. For different snow types, we find an interesting mixture of diffusive- and jump-like noise. We propose the interpretation that the dominant diffusive process is due to the pushing aside of many snow grains, whereas the breaking of ice structures leads to dominant jump noise. Our results show that the denser structures typical of DH and RGlr lead to a more collective diffusive behavior, whereas for the highly porous snow structures of PP and RG, the single breaking events lead to a relatively strong jump noise. For this interpretation, all R values were detrended and normalized; thus the absolute values of the snow hardness R are not essential but more the resulting collective behavior of the snow types.

Finally, we would like to point out that our characterization of a complex material, snow, by a penetration process should have the potential to be generalized to, for example, biological tissue or ground layers. Last but not least, we would like to point out that our work provides additional insight into analyzing and modeling the complex nature of snow types, complementing existing methods.