# Acoustic Probing of the Jamming Transition in an Unconsolidated Granular Medium

Xavier Jacob, Vladislav Aleshin, Vincent Tournat, Philippe Leclaire, Walter Lauriks, Vitali Goussev

► To cite this version:

Xavier Jacob, Vladislav Aleshin, Vincent Tournat, Philippe Leclaire, Walter Lauriks, et al., Acoustic Probing of the Jamming Transition in an Unconsolidated Granular Medium. Physical Review Letters, 2008, 100 (15), pp.158003. 10.1103/PhysRevLett.100.158003. hal-01326055v3

HAL Id: hal-01326055

https://hal.science/hal-01326055v3

Submitted on 8 Nov 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Acoustic Probing of the Jamming Transition in an Unconsolidated Granular
Medium

X. Jacob, V. Aleshin, V. Tournat

Laboratoire d'Acoustique, UMR-CNRS 6613, Universite du Maine, Avenue Olivier Messiaen, 72085 Le Mans, France

P. Leclaire,

Laboratoire de Recherche en Mécanique et Acoustique, I.S.A.T - Université de Bourgogne, 49, rue Mademoiselle
Bourgeois, 58000 Nevers, France

W. Lauriks

Laboratorium voor Akoestiek en Thermische Fysica, Katholieke Universiteit Leuven, Celestijnenlaan 200D,
BE-3001 Heverlee, Belgium

V. E. Gusev*

Laboratoire de Physique de l'Etat Condensé, UMR-CNRS 6087, Universite du Maine, Avenue Olivier Messiaen, 72085
Le Mans, France

(received

Experiments with acoustic waves guided along the mechanically free surface of an unconsolidated granular packed structure provide information on the elasticity of granular media at very low pressures that are naturally controlled by the gravitational acceleration and the depth beneath the surface. Comparison of the determined dispersion relations for guided surface acoustic modes with a theoretical model reveals the dependencies of the elastic moduli of the granular medium on pressure. The experiments confirm recent theoretical predictions that relaxation of the disordered granular packing through non-affine motion leads to a peculiar scaling of shear rigidity with pressure near the jamming transition corresponding to zero pressure. Unexpectedly, and in disagreement with the most of the available theories, the bulk modulus depends on pressure in a very similar way to the shear modulus.

PACS numbers: 43.20.+g, 45.70.-n, 46.40.Cd, 43.20.Mv

Transition from fluidity to rigidity in disordered
media is a poorly understood phenomena in mechanics.
This problem keeps attracting the attention of researchers
studying its particular manifestations such as the
transition to glassiness in liquids [1,2] and to jamming in
colloidal dispersions or granular packings [2,3]. General
interest in the unjamming transition in unconsolidated
granular media is supported by its macroscopic
manifestation in nature in the form of avalanches of sand
and snow [4].

Recent theoretical investigations [1,3] as well as
numerical experiments with packing of elastic grains [5-
10] have predicted that in the vicinity of the jamming
transition in the solid phase the elasticity of a disordered
material is controlled not just by the elasticity of the
individual contacts between the grains but also by
material relaxation via non-affine motion of the grains. In
particular, the dependence of the shear rigidity G on
pressure p deviates from the law G ∝ p¹/³ that follows
from the Hertz-Mindlin theory of contacts between
spheres [11,12]. The theory [3] indicates that the shear
modulus G should also be proportional to the excess
coordination number of contacts δz with respect to the
critical isostatic number zc at the jamming transition,
where the number of the contacts in the system is equal to
the number of force balance equations [1]. This results,
due to the relation δz ∝ p¹/³ [3], in G ∝ p²/³, while the
bulk elastic modulus B does not contain additional
scaling with pressure: B ∝ p¹/³. The scaling exponent 2/3
had earlier been found numerically for random packings
of frictionless particles at zero temperature interacting
through the repulsive Hertzian potential [5,6]; the
corresponding range of pressures p evaluated numerically
can be estimated for glass beads to be 100 Pa ≤ p ≤ 10 MPa.
These results have been confirmed for p ≥ 100 kPa in
Ref. [7,8]. Moreover very recently numerical experiments
[9,10] have demonstrated that the scaling is not modified

by friction at least for the 2D packing of spheres. The
deviation from the scalings G ∝ p^{α_G} (with α_G = 2/3) and
B ∝ p^{α_B} (with α_B = 1/3) is predicted [7,10] only at
higher pressures p ≥ 10 MPa, where the condition
δz << z_c (i.e. the vicinity of the jamming transition) no
longer holds.

To the best of our knowledge there has only been
a single attempt to experimentally verify these theoretical
simulations, using photoelastic disks near the 2D
jamming transition [13], in which reasonable agreement
with theory was found. However the features predicted
for the packing of spheres have not been tested
experimentally. The measurement of the dependence on
pressure of the shear and longitudinal sound velocities c<sub>s,L</sub>,
which under the condition of negligible variations of
density scale as the square root of the elastic moduli, is a
classical method for the evaluation of the state of an
unconsolidated granular medium [14-16]. That the
behavior c<sub>s</sub> ∝ G<sup>1/2</sup> ∝ p<sup>α<sub>s</sub></sup> with α<sub>s</sub> = 1/3 has not been revealed
in most of the earlier experiments could be attributed to
an important deviation from the jamming threshold
pressure p = 0 (for example, p ≥ 2 MPa in [16],
p ≥ 5 MPa in [9]). However, the exponent α<sub>s</sub> = 1/3 has
not been documented even in experiments conducted at
much lower pressures ( 5 kPa ≤ p ≤ 50 kPa )[17] although
there are indications that α<sub>s</sub> > 1/4 holds over a limited
range of pressures 5 kPa ≤ p ≤ 20 kPa . In a very limited
interval 50 kPa ≤ p ≤ 100 kPa the exponent α<sub>s</sub> = 0.3, which is
closer to 1/3 than to 1/4, has been proposed to fit the
experimental data in Refs. [14,15]. It could be concluded
that there are indications of fast scaling (α<sub>s</sub> > 1/4) in
some acoustic experiments, but this phenomena has never
been observed over a wide enough interval of sufficiently
low pressures for the precise evaluation of the scaling
exponent.
To access the elastic moduli of granular media at much lower pressures than usual we use waveguide acoustic modes propagating along the mechanically free surface of an unconsolidated granular medium [18,19]. The localization of these modes near the surface is due to the bending of the acoustic rays towards the surface caused by an increase with depth of the rigidity of the gravity-loaded granular packing (the acoustic mirage effect²⁰). Through experiments at ultrasonic frequencies (500 Hz ≤ f ≤ 5 kHz, an order of magnitude higher than in earlier experiments [21,22]), we are able to test near-surface layers at low pressures corresponding to the vicinity of the jamming transition. The lowest-order waveguide mode at its highest frequency f = 3 kHz propagates at a velocity ~15 m/s and penetrates only ~0.5 cm beneath the surface, probing the material at pressures less than 75 Pa. This is below the lowest pressure level accessed in numerical simulations [7,8]. However, even at the bottom of the experimental container (Fig. 1) at a depth of 20 cm, which can be accessed through the low-frequency higher-order modes, the pressure does not exceed 7.5 kPa, and is lower than in the most of the earlier experiments with bulk acoustic waves. The new experimental method proposed here to study such properties of granular materials consists in measuring at the surface of a granular layer the motion associated with the propagation of guided surface acoustic modes. The waves are excited using a 1 mm thick aluminium plate attached to a shaker and partly buried in the granular material (Fig. 1). The vertical surface velocity as a function of time and distance is measured with the help of a laser Doppler vibrometer with a surface focused beam (sensitivity 1mm/s/V within the frequency band 80 Hz – 20 kHz). The particle velocity is recorded every millimeter along the propagation direction from the source over a 300 mm distance. Typically, a maximum vertical velocity of 5.10⁻⁴ m/s at 1 kHz was observed close to the source. As the mechanical nonlinearity of a granular medium increases on approaching the jamming point, [6,17] we checked that the data are not influenced by the wave amplitude in the chosen range of excitation strengths. Using the measured signal in the time-space domain it is straightforward to perform the time and space double Fourier transform of the signals and to obtain frequency-wavenumber signals and plots of the dispersion curves for different guided acoustic modes. The experiments are performed on a granular material consisting of glass beads 150 µm in diameter in a large tank of dimensions 80 cm × 50 cm × 20 cm (see Fig. 1). The sample is prepared before each experiment as follows: the container is filled up with the grains and is shaken gently. Then acoustic excitations of different amplitudes and frequency ranges are applied for several hours until a stable response of the medium is obtained.

A typical experimentally obtained pattern in the $(\omega, k)$-domain is presented in Figure 2. The range $0<k<50-100 \text{ m}^{-1}$ is probably affected by the reflection of the waves from the container walls, but the short-wave part corresponds to the modes guided in the horizontal channel [18,19,22].

The dispersion relations in Fig. 2 are plotted using an intensity scale corresponding to the logarithm of the normalized Fourier transform of the particle velocity. In order to further increase the contrast,

threshold filtering was performed. The measurements were repeated a number of times with different excitation functions, namely short pulses with Gaussian spectra centered at different frequencies. For each dataset, maxima of the smoothed intensity function for different k-points were found. In 80% of the different experimental conditions these maxima were in good agreement with each other. We selected only these points and calculated the power-law approximation curves together with the standard deviation for each mode. In Fig. 2 the solid lines represent the averages ± standard deviations.

<figure><img src="image_2.png" /><figcaption>FIG. 1. Experimental setup.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIG. 2. A typical dispersion curve. Thin lines bounding the “tubes” indicate average + standard deviation and average - standard deviation curves obtained as a result of a statistical treatment of a number of patterns.</figcaption></figure>

For the interpretation of the experimental data, guided surface acoustic modes (GSAMs) [18,19] polarized in the sagittal (vertical) plane are modeled with the continuous medium approximation assuming local isotropy of the elastic properties. The horizontal $u_x$ and the vertical $u_y$ components of the mechanical displacement vector $\bar{U} = \{u_x(y), u_y(y)\}e^{i\omega t-ikx}$ are controlled in microscopically a inhomogeneous medium with vertical stratification $c_{S,L} = c_{S,L}(y)$ by a system of coupled Helmholtz equations [18,19]

$$ (c_s^2 u_x')' + (\omega^2 - c_L^2 k^2) u_x = ik \left[ (c_L^2 - 2c_s^2) u_y' + (c_s^2 u_y')' \right], $$

$$ (c_L^2 u_y')' + (\omega^2 - c_s^2 k^2) u_y = ik \left[ c_s^2 u_x' + ((c_L^2 - 2c_s^2) u_x')' \right]. $$

Here the prime denotes the y-derivative. In addition to conditions of zero stress tensor at the mechanically free boundary (y=0) the eigen-problem for GSAMs includes the conditions of localization ($u_{x,y}(y \to \infty) \to 0$). In unconsolidated granular media, the stratification of sound velocities is controlled by the dependence of the elastic moduli on pressure, which increases linearly with depth because of gravity ($p = \rho gy \propto y$), where $\rho$ is the density and g is the gravitational acceleration).
Assuming $c_{S,L} = \gamma_{S,L}(\rho_g y)^{\alpha_{S,L}}$, where $\gamma_{S,L}$ are depth-independent coefficients, normalizing the depth-coordinate to the inverse wave number $k^{-1}$ ($ky \Rightarrow y$) and using the substitution $iu_y \Rightarrow u_y$, we transform the governing equations in a form convenient for numerical evaluation:

$$ (y^{2\alpha_s} u_x')' + (\Omega^2 - \delta y^{2\alpha_L}) u_x = (\delta y^{2\alpha_L} - 2y^{2\alpha_s}) u_y' + (y^{2\alpha_s} u_y)' , \quad (1) \\ \delta(y^{2\alpha_l} u_y')' + (\Omega^2 - y^{2\alpha_s}) u_y = -y^{2\alpha_s} u_x' - (\delta y^{2\alpha_L} - 2y^{2\alpha_s}) u_x' . $$

Here $\Omega = \omega/(\gamma_S k)(k/\rho_g)^{\alpha_S}$ is the normalized frequency and $\delta = (\gamma_L/\gamma_S)^2 (\rho_g/k)^{2\alpha_L-2\alpha_S}$ is a dimensionless parameter equal to the square of the ratio of the longitudinal to shear velocities at a depth of about the GSAM penetration ($y \propto k^{-1}$). In the previously studied case, [18,19] $\alpha_S = \alpha_L = \alpha$, there was no characteristic spatial scale in the system; the eigen frequencies $\Omega_n$ ($n=1, 2, ...$) of GSAMs and the acoustic field spatial distributions (eigen modes) obtained through the numerical solution [19,20] did not depend on $k$ (because for $\alpha_S = \alpha_L$ the parameter $\delta$ is $k$-independent). The power-law dispersion relation for GSAMs then follows from the definition of $\Omega$: $\omega_n = \Omega_n \gamma_S (\rho_g)^\alpha k^{1-\alpha} \propto k^{1-\alpha}$. The solutions lose their self-similarity [18,19] as soon as $\alpha_S \neq \alpha_L$; in this case a characteristic spatial scale $y_0$ might be formally defined by $c_S(y_0) = c_L(y_0)$, and the eigenfunctions of the problem (1) would start to depend on $k$. However, theoretically the self-similarity should be recovered in the asymptotic limits $\delta \ll 1$ and $\delta \gg 1$. The former case is impossible in physical reality because $\delta > 1$ due to $c_L \propto \sqrt{B+(4/3)G} > c_S \propto \sqrt{G}$. For the theoretically predicted relation $\alpha_S > \alpha_L$ the inequality $\delta \gg 1$ could be expected for sufficiently short-wavelength GSAMs due to $\delta \propto k^{2(\alpha_S - \alpha_L)}$. In the asymptotic case $\delta \gg 1$, similar to the situation in which acoustic waves propagate at the interface of the ocean with unconsolidated water-saturated sediments [23], the system (1) can be transformed into equations describing two types of waveguide modes: fast and slow.

The modes, which we call here fast modes (in accordance with the notation of fast seismo-acoustic waves [23]), are purely compressive; their description follows from Eq. (1) when shear rigidity is completely neglected ($c_S = 0$). Correspondingly, the fast modes are controlled by the distribution of the velocity $c_L(y)$, and the dispersion relation scales as $\omega \propto k^{1-\alpha_L}$. These modes can be easily detected in experiments with normal laser incidence since they have a significant vertical component of the surface displacement.

However in the same limiting case $c_S \to 0$, $\delta \to \infty$ there exists a second type of waves termed here the slow modes. They are controlled only by the profile $c_S(y)$ and have a dispersion relation scaling as $\omega \propto k^{1-\alpha_S}$, but they are not pure shear and contain a compression contribution ($\text{div}\vec{U} \propto \delta^{-1}$). However, as it has been shown by our numerical analysis of Eq. (1), the corresponding vertical

component of the surface displacement is typically very weak except for very low order modes.

<figure><img src="image_4.png" /><figcaption>FIG. 3. Experimental dispersion curves (solid lines) represented as $\pm$ standard deviation intervals and theoretical points fitted to them. The numbers indicate mode orders. The bunch of solid grey and black lines provides for comparison the theoretical slopes $5/6=1-\alpha_L^{\text{theor}}$ and $2/3=1-\alpha_S^{\text{theor}}$. The inset shows a plot of the theoretical normalized ratios of the vertical component of the particle displacement vector at the surface to the modulus of this vector, as a function of mode number.</figcaption></figure>

In the real experimental situation the condition $\delta \gg 1$ is fulfilled with a limited precision: our fits (Fig. 3) of the experimental data (Fig. 2) using the numerical solution of Eq. (1) provide $\delta \approx 4.6-4.9$ only. Thus, although the experimentally detected modes are neither purely fast nor slow, some of the theoretical dispersion curves for $\delta \approx 4.6-4.9$ are mainly determined by $\alpha_S$ and $\gamma_S$ (quasi-slow modes), while the others are mainly determined by the values of $\alpha_L$ and $\gamma_L$ (quasi-fast modes). In this sense, the modes 1-10, 13-22 and 25-30 in our analysis covering the first 30 modes are quasi-slow and thus poorly visible except for some low order modes, whereas the modes 11-12 and 23-24 are quasi-fast and easily detectable. Variations of about 2% in the parameters $\alpha_{S,L}, \gamma_{S,L}$ shift these mode numbers by 1-2. The insert in Fig. 3 shows the theoretical vertical projections of the polarization vector at the surface, i.e. the vertical components of the particle displacement normalized on the displacement modulus, for the modes of different order. This could be used to qualitatively estimate the visibility as measured by the laser vibrometer. We identify the first four lowest “tubes” in Fig. 2 with the first four lowest order ($n=1-4$) GSAMs, and the two other tubes with the four quasi-fast modes 11, 12, 23, 24. The other modes are invisible in experiment (Fig. 2) and are not plotted in Fig. 3.

Matching of the curves is typically realized as the minimization of an objective function, penalizing only points coming out of the corresponding experimental intervals and leaving complete freedom to the points within intervals. For modes higher than the fourth, only points with visibility greater than 0.5 were considered. The minimization was performed by the principal axis method (PRAXIS) developed by Brent [24] and also known as the modified Powell algorithm, which is one of the most effective amongst methods of function optimization without analytical knowledge of derivatives.
This procedure has provided the optimized parameters of the rigidity profiles in Eq. (1): $\gamma_s^{\text{exp}} = 6.43 \pm 0.13$ [SI], $\gamma_L^{\text{exp}} = 14.8 \pm 0.4$ [SI], $\alpha_s^{\text{exp}} = 0.320 \pm 0.006$, and $\alpha_L^{\text{exp}} = 0.305 \pm 0.01$ with the residual value of the objective function around $200 \text{ s}^{-1}$. The errors given above have the following sense: if the fitted value is perturbed by the value of the error, the objective function doubles.

Our fitted scaling exponents $\alpha_s^{\text{exp}} \approx 0.96\alpha_s^{\text{theor}}$, $\alpha_L^{\text{exp}} = 0.64 \pm 0.01 \approx 0.96\alpha_L^{\text{theor}}$ provide the first experimental confirmation for the theoretically predicted behavior of shear rigidity of the 3D disordered packing of spherical elastic beads near the jamming transition $\alpha_G^{\text{theor}} = 0.666...$. At the same time, the experimentally revealed scaling exponent of the bulk modulus $\alpha_B^{\text{exp}} = 0.61 \pm 0.02 \approx 1.83\alpha_B^{\text{theor}}$ is much closer to the scaling exponent of the shear modulus than to theoretical value $\alpha_B^{\text{theor}} = 0.333...$. Our attempt to impose the theoretical value $\alpha_L$ and compensate this misfit by matching the other parameters $\gamma_{s,L}$ and $\alpha_s$ led to an important deviation from the experimental data and to a drastic increase in the residual objective function. We have estimated that the wave velocity dispersion due to medium microinhomogeneity at the scale of the grain size and the diffraction of GSAMs caused by the finite dimensions of our acoustic source [both are not taken into account by the model (1)] give a negligible contribution to the experimentally determined dispersion relations in Fig. 2. We also verified by examining the depth profiles of all eigenmodes visible in experiment that their contact with the bottom of the container can be neglected and that the results obtained for the half-space are self-consistent. The experiment therefore indicates that in the vicinity of the jamming transition near the mechanically free surface, relaxation of the granular media via non-affine motion in response to compression could be very similar to that in response to shear loading. This is in contradiction to the theoretical predictions obtained earlier [3,5-10] for macroscopically homogeneous and isotropic media.

The key to understanding of this discrepancy is the realisation that the granular medium in our experiments is both macroscopically inhomogeneous (not only because of the vertical stratification, but also just because of the presence of the surface [25]) and macroscopically anisotropic [17,26] (because of the preferential direction of its loading by the gravity field). Even in continuum mechanics, inhomogeneity — including the surface as a particular case — and anisotropy induce, through the breaking some symmetry constraints, the coupling of compression and shear motion. Acoustic mode conversion $S \leftrightarrow L$ takes place at each point of the inhomogeneous media. Pure compression and shear waves can be the eigenmodes of homogeneous anisotropic media only along particular directions, otherwise the eigenmodes are quasi-longitudinal and quasi-shear, combining both compressional and shear motion. We suggest that it could be macroscopic inhomogeneity and anisotropy that through breaking some symmetry constraints couple non-affine motions caused by compression and shear loading on the microscale, allowing disordered granular media under local compression and shear to relax in a similar way. To verify this hypothesis the numerical simulations of the type

described in Refs. [5-10] need to include the effects of macroscopic inhomogeneity and anisotropy, and the model in Eq. (1) should be extended to the anisotropic case. The rattlers, particles that do not overlap with the other particles [6] and commonly excluded from the analysis of scaling exponents in granular media [6,13], could also (as any other type of defects) couple compression and shear, contributing an extra relaxation pathway for the granular medium under compressional loading. Equal exponents $\alpha_G$ and $\alpha_B$ are in fact expected according to disordered lattice models for rigidity percolation where all elastic moduli obey the same scaling behavior (see [27] and the references therein).

The success in fitting our experimental results using theoretical predictions for the static (equilibrium) shear modulus indicates that the characteristic time scale $\tau_{n-a}$ of the disordered granular packing non-affine relaxation near the jamming transition ($p \le 7.5 \text{ kPa}$) does not exceed few microseconds, as is evident from the quasi-equilibrium condition $\omega\tau_{n-a} \ll 1$. To determine the relaxation time $\tau_{n-a}$ or a distribution of relaxation times, higher acoustic frequencies than those reported here should be used.

This study was supported by ANR project No. NT05-341989.

\*Electronic address: vitali.goussev@univ-lemans.fr

[1] S. Alexander, Phys. Rep. **296**, 65-236 (1998).

[2] L. E. Silbert, A. J. Liu, S. R. Nagel, Phys. Rev. Lett. **95**, 098301 (2005).

[3] M. Wyart et al. Phys. Rev. E **72**, 051306 (2005); M. Wyart, Ann. Phys. Fr. **30**, 1-96, (2005).

[4] R. A. Bagnold, *The Physics of Blown Sand and Desert Dunes* (Methuen, London, 1941).

[5] C. S. O'Hern et al. Phys. Rev. Lett. **88**, 0755507 (2002).

[6] C. S. O'Hern et al. Phys. Rev. E **68**, 011306 (2003).

[7] H.A. Makse et al. Phys. Rev. E **70**, 061302 (2004).

[8] W.G. Ellenbroek et al. Phys. Rev. Lett. **97**, 258001 (2006).

[9] Shundyak, M. van Hecke, W. van Saarloos, Phys. Rev. E **75**, 010301(R) (2007).

[10] E. Somfai et al. Phys. Rev. E **75**, 020301(R) (2007).

[11] K. L. Johnson, *Contact Mechanics* (Cambridge University Press, Cambridge, 1985).

[12] R.D. Mindlin and H. Deresiewicz, J. Appl. Mech. **20** 327-344, (1953).

[13] T.S. Majmudar et al. Phys. Rev. Lett. **98**, 058001 (2007).

[14] B.O. Hardin and F.E. Richart, J. Soil Mechanics and Foundations Division, Proc. ASCE, **89** (SM1), 33-65 (1963).

[15] B.O. Hardin and W.L. Black, J. Soil Mechanics and Foundations Division, Proc. ASCE, **92** (SM2), 27-42 (1966).

[16] S.N. Domenico, Geophysics **42**, 1339-1368 (1977).

[17] V. Tournat et al. Phys. Rev. Lett. **92**, 085502 (2004).

[18] V. E. Gusev, V. Aleshin, V. Tournat, Phys. Rev. Lett. **96** 214301 (2006).

[19] V. Aleshin, V. Gusev, V. Tournat, J. Acoust. Soc. Am. **121**, 2600-2611 (2007).

[20] C-h. Liu and S. R. Nagel, Phys. Rev. Lett. **68**, 2301-2304 (1992).

[21] B. Andreotti, Phys. Rev. Lett. **93**, 238001 (2004).

[22] L. Bonneau, B. Andreotti, E. Clement, Phys. Rev. E **75**, 016602 (2007).

[23] O. A. Godin and D. M. F. Chapman, J. Acoust. Soc. Am. **110**, 1890-1907 (2001).

[24] R. P. Brent, *Algorithms for Minimization without Derivatives*, (Prentice-Hall, Englewood Cliffs, 1973) Chapter 7.

[25] V. V. Krylov, Progress in Surf. Science **32**, 39-110 (1989).

[26] X. Zheng and B. Ni, J. Geotech Geoenviron. Eng. **125**, 741-749 (1999).

[27] D. A. Head et al. Phys. Rev. E **68**, 025101(R) (2003).