1. Introduction

Saltation was divided up by Reference BagnoldBagnold (1941) into four distinct subprocesses: aerodynamic entrainment, the grain trajectory, the splash process and wind modification. The splash process is when a grain impacts the bed and may ricochet and eject additional particles. This process is critical in determining the mass flux and response times of a saltating system. The main contribution of the wind is to accelerate the grains horizontally, at least at friction velocities where suspension is negligible, and without the conversion of this to vertical motion in the splash process saltation would rapidly die out. Indeed, in equilibrium saltation, aerodynamic entrainment may not occur at all (Reference Doorschot and LehningDoorschot and Lehning, 2002) and the population of saltating grains is only maintained by the splash process. An early analytic investigation of the splash function was performed by Reference RumpelRumpel (1985), who considered idealized collisions between regularly packed discs in two dimensions. Reference Willetts and RiceWilletts and Rice (1985) have performed detailed experiments on the splash function for sand. Reference Nalpanis, Hunt and BarrettNalpanis and others (1993) summarized many experiments and showed that they were consistent and that there was agreement with experiments on snow (Reference Araoka and MaenoAraoka and Maeno, 1981).

Despite the extensive work that has been carried out, published splash functions suffer from many difficulties both theoretical and practical. In this paper we will use the term splash function to refer to the joint probability density function (jpdf) of the velocities after impact given the impact velocity; thus

(1)

where u is the velocity of the impacting grain, v the velocity of the same grain after the collision and v_i are the velocities of any grains which reach a certain height above the bed after the collision. This definition ensures that the function is well defined, and highlights the fact that a certain separation is necessary if one wishes to treat the collision as a discrete process rather than a continuous process in which energy is radiated into the bed as a compression wave. Many published splash functions use unbounded distributions for splash velocities so that the energy of the outgoing particles can be greater than that of the incoming ones. A physically sensible splash function should satisfy

(2)

so that energy is never created. The splash function should respect appropriate symmetries and should be defined for all angles and velocities. For example, the dependence on angles should only enter through trigonometric functions so that it is periodic. Of particular importance is correctly accounting for the dependences and correlations between the variables. For example, it is not sufficient to calculate separate probability density functions (pdf’s) for ejection angle and ejection velocity, unless it is shown that they are independent, which they are not. The correlations between the variables are of course critical in determining the efficiency with which horizontal momentum is converted into vertical movement. The splash function, as defined in this paper, should have no dependence on the wind properties. This is because collisions in the bed occur over such small time periods that aerodynamic forces are unimportant. The observed splash velocities will, however, depend on the wind properties since these determine the distribution of the impact velocities u. If f_i(u) is the distribution of impacting particles, the observed splash distribution will be

(3)

This is a significant complicating factor in estimating the splash function, which needs to be properly accounted for.

All published splash functions (to the authors’ knowledge) separate the splash function into a binary probability, for whether or not a particle is observed, and then a conditional distribution function for its velocity. This is a useful empirical approach, but it is objectionable theoretically. It is not at all clear that the distribution of ejection velocities should be bimodal. A much more natural approach is to have a unimodal ejection velocity pdf, but then regard particles that have negative vertical velocities as being absorbed by the bed.

An important practical difficulty is due to censoring of the data. Most experiments have observed collisions in a two-dimensional slice illuminated by a light sheet. If the particles have a large velocity component normal to this plane they will not be observed. Another very important censoring mechanism is that particles can only be observed if they move clear of other particles in the bed. The length of a particle streak in a video frame is proportional to its speed, so that slow particles will leave short streaks and also be harder to observe. For example, suppose that, due to roughness in the bed, the grain-size and the pixel size, only particles that rise up h = 2 mmare observed. This corresponds to a speed of so in this case no particles that are splashed with a speed less than this would be observed. Even if a particle is observed, its measured velocity will be reduced by

(4)

A carefully designed experiment should include an analysis of the importance of these effects.

The most serious problem in developing a sound splash function is the lack of data. How many data would be necessary to accurately map the splash function? Suppose we consider just the impacting particle so that there are six degrees of freedom (three for the impact velocity and three for the ejection velocity). If we divided up each variable into ten ranges and wanted on average ten samples per region, we would need 10⁷ samples to be able to perform a chi-squared test of distribution. If we assume that the bed is isotropic, that is only the speed and angle to the vertical of the incoming particle are important, then this reduces the free variables by one, so only 10⁶ events would be necessary. If we only wanted to consider the behaviour in a plane −and this is all we can do with two-dimensional data−another degree of freedom would go and 10⁵ events would be sufficient. A fully automated experiment that could measure this number of events would be extremely interesting, but in this paper we are restricted to datasets containing O(103) elements, so a full mapping and hypothesis testing on the splash function is not possible and will not be attempted. Instead we will examine the correlations between different variables and attempt to make a choice that simplifies the problem. Given the small size of the datasets, we will emphasize developing a theory with very few parameters that fits the data approximately.

2. Method

This paper uses data from experiments described in Reference Sugiura and MaenoSugiura and Maeno (2000). Experiments were conducted in two wind tunnels in Japan. Compact snow, which had been kept for a year at –15°C, was used in the wind tunnel of the Institute of Low Temperature Science, Hokkaido University, whereas fresh snow was used in the wind tunnel of the Shinjo Branch of Snow and Ice Studies, National Research Institute for Earth Science and Disaster Prevention. The compact snow had mean diameter 0.35 mm with standard deviation 0.28 mm. The fresh snow was dendritic, with dimensions of 0.5–5.0mm. The experiments with the fresh snow were carried out with a friction velocity of 0.17 ms^–1, and the experiments with compact snow were carried out with friction velocities of 0.19 and 0.25 ms^–1. At the downwind end of the wind tunnel, impacting and ejecting grains were illuminated by a 10 mmwide laser sheet cut with a rotary shutter at intervals of 1 ms for natural snow and 0.5 ms for fresh snow. The illuminated particles were recorded using a video camera, and the velocity before and after impact calculated from the trajectories.

3. Results

The software for analyzing the trajectories and calculating the velocities calculates values quantized along the vertical and horizontal axes of the camera. For compact snow, there were 1366 usable* collisions (out of 1385) and the velocities were quantized at an angle rotated by 0.83° spacing 0.079 ms^–1corresponding to pixel size of 0.079mm. For fresh snow, there were 1629 usable collisions (out of 1634) and the velocities were quantized at an angle rotated by 0.12° spacing 0.063 ms^–1 corresponding to pixel size of 0.13mm. Errors are presumed to be limited to the quantization error. Velocity is calculated from the trajectories extrapolated to the bed surface, so no sampling bias is expected in the data as would occur if the velocity was calculated at a height above the bed. When no grain was observed to be ejected from the bed, the speed was taken to be 0.

The observed data are shown in Figure 1–4. The two figures of ejection velocities show that there are very few observed non-zero velocities within the ellipse

Fig. 1. Impact velocity for compact snow. Area is proportional to number. Total is 1366.

Fig. 2. Ejection velocity for compact snow. Area is proportional to number. Total is 1366.

Fig. 3. Impact velocity for fresh snow. Area is proportional to number. Total is 1629.

Fig. 4. Ejection velocity for fresh snow. Area is proportional to number. Total is 1629.

(5)

This may be correct, but it is much more likely that in fact there were particles with these velocities but that they were unobserved, that is they were censored. This is for the reasons discussed in the introduction, which are that if a particle is moving slowly it only leaves a short streak on a video frame, which is less likely to be observed, and also that a particle must move up high enough from the bed to no longer be obscured by other grains. This assumption is made for the rest of the paper.

4. Theory

In this section, we consider the velocity pdf of only the ejected particle with the largest velocity, and we call this the ricochet function f_r(v; u). This can be obtained from the complete splash function by integrating out the velocities v_i of the ejected particles. In nearly all cases, we expect this to be the impacting particle but this is not necessary. We choose this approach to describe the splash function since we imagine a cascade of collisions between particles in the bed, with the energy and velocities decreasing after each one. Thus it is natural to regard the smaller ejection velocities as being dependent on the larger ones rather than vice versa.

4.1. Normal impulse model

There is a wide choice of different variables that can be used to describe the velocities of a particle before and after a collision: Cartesian or polar coordinates can be used; various restitution coefficients and rotations can be defined. They are all equivalent since a theory defined using one set of variables can be transformed into another. However, some choices of variables may have fewer correlations and have a simple physical interpretation. For example, if the velocities after a collision are roughly proportional to the impact speed, then restitution coefficients will be independent of the impact speed and an important simplification is obtained.

Figure 5 shows two collisions of grains off a surface. If the grain is spherical, the surface hard and smooth and the surface properties can be described by a restitution coefficient e and Coulomb friction coefficient μ then it is straightforward to calculate the ejection linear and angular momentum from the impact linear and angular momentum. Of course these criteria are not satisfied for a grain impacting a bed, but we will show that regarding a collision in these terms defines two useful variables, which we call impulse normal n̂ and impulse restitution e_n. We define these as the surface normal and restitution that would result in the observed velocities for an ideal collision with zero friction. Thus

Fig. 5. Schematic of two ricochets.

(6)

This can be inverted to give

(7)

(8)

Note that the e lies in the range [-1; 1]. The presence of surface friction will act to rotate the impulse normal away from the surface normal towards the incoming particle. The impulse restitution has a simple interpretation as a ratio of masses. If the impacting particle has mass m and it collided with a particle of mass M with restitution coefficient e then

(9)

Thus we would expect that the upper limit on the observed e_n would be the true material restitution e. Equation (9) also shows that if the part of the bed the grain impacts consists of a group of grains joined by stiff contacts, due to sintering for example, then M/m will be large and the impulse restitution will be close to e. If, however, it hits a smaller grain on the top of a ripple so that M/mis small, the grain will move through and e_n will be close to –1. This also suggests how to apply the theory to grains of different sizes. Since no information about the size of impacting grains was collected in these experiments, we will not pursue this idea. But if such data were to be collected they could easily be included by working with the mass ratio rather than e_n.

4.2. Correlations

To see how the variables are correlated, we calculated Kendall’s tau coefficient for all variable pairs. This is a non-parametric test that is sensitive to all monotonic correlations (Reference Hollander and WolfeHollander and Wolfe, 1973). Table 1–Table 3 show the number of standard deviations from 0 for Kendall’s tau normalized by the diagonals. Positive values indicate positive correlations, and negative values negative ones. Values small in absolute magnitude indicate lack of a monotonic correlation and suggest that the variables are independent. There is a very important caveat in interpreting these data. No allowance has been made for separating out the correlations in the impact velocities because of the lack of data, and this problem will not be further discussed. Also correlation induced by censoring of the data will be ignored. A more detailed study should carefully address these issues.

Table 1. Kendall’s normalized tau coefficient of correlation for different variable pairs for compact snow

Table 2. Kendall’s normalized tau coefficient of correlation for different variable pairs for fresh snow

Table 3. Definition of variables used in Tables 1 and 2

One can see for example that there is a strong correlation between the ejection vertical velocity υ_y and the impact horizontal velocity for u_x. It is also shown that there is large negative correlation between the ejection angle and the total restitution and speed, which would have to be considered if we were to construct a jpdf for them. The horizontal and vertical restitution coefficients are also highly correlated. Rather surprisingly, ejection horizontal and vertical speeds do not appear to be strongly correlated, at least for the compact snow. This is odd because there must be a strong negative correlation, at least in collisions when little energy is dissipated, in order that energy is not created. A very pleasing result is that the impulse angle and restitution are only very weakly correlated so we might hope to model them as independent. There is also only a small correlation between them and the impact velocity, and the impulse restitution does not depend strongly on the impact angle. These statements are all much less true for the fresh snow than the compact snow. We hope that this is due to greater problems with censored data in the more dissipative fresh-snow case.

This analysis of the correlations suggests that we base a theory on the impulse angle and impulse restitution and that they are treated as independent. Since we are trying to find a simple approximate theory and the correlations are low, we assume that both are independent of the impact velocity and that only the impulse angle depends on the impact angle.

Suppose that is the distribution of surface normals over the bed. Then the distribution of impulse normals will be approximately proportional to this but weighted by a factor where is the unit vector in the direction of the impacting particle. Friction during the collision will also weight the distribution of towards as will shadowing effects by ripples. That is to say, a particle sees more of the surface that is normal to its motion. This suggests that the distribution of impulse normals should be of the form

(10)

No direct measurements were taken of the surface properties, so f_n is unknown, but it would be very interesting to try to measure this in future experiments. Instead we assume that it the very simple and isotropic form Thus we have

(11)

This form is valid for three dimensions, but since velocities normal to the light plane were not observed we restrict the analysis to two dimensions which, ignoring censoring, is equivalent to integrating out this direction (this changes β). The distribution then becomes

(12)

with θ_i ϵ [0; π/2] being the impact angle to the horizontal and θ_n the impulse angle from the horizontal measured clockwise. This distribution is similar to a beta distribution but differs because of the dependence on θ_i.

For the impulse restitution, the natural choice of distribution is beta since it is in the range [–1, 1]. For a less elastic material, the upper limit should be taken as the coefficient of restitution or the theory developed in terms of the apparent mass ratio (Equation (9)). Since the correlations with impact velocity and impulse angle are small, we assume that the parameters are constant for all collisions and only depend on the bed properties. Thus we take

(13)

Combining these distributions and including the Jacobian factor (which would be slightly different in three dimensions) de dθ = dv/[u ⋅ (u - v)], the pdf for f(v) is

(14)

where N is the normalization factor, which does not have a simple closed form. If there were not the censoring issue to contend with, it would be straightforward to calculate the parameters independently for angle and restitution. Instead Equation (14) is modified so that all negative vertical velocities and those falling within the ellipse defined by Equation (5) are regarded as having velocity 0. The parameters were then found using maximum likelihood estimation and are shown inTable 4.

Table 4. Splash function parameters found using maximum likelihood estimation and assuming elliptical censorship

Figure 6–9 compare the predicted distributions and the observations. In each of these figures the solid curve shows the predicted distribution if there were no censoring of the data, the lefthand histogram is the observed data, and the righthand histogram is the predicted observations with censoring. The substantial effect of the censoring on the data is clear. In particular, no impulse angles less than half the impact angle will be observed, as these result in negative vertical velocities after the impact, so that a particle would be absorbed. It is of course possible to get much better fits to these particular histograms by fitting directly to them, but these figures show the results of fitting the complete distribution for each collision so that the dependence on the impact velocities and the correlations will be correct.

Fig. 6. Observed impulse-angle distribution compared to fitted distribution for fresh snow.

Fig. 7. Observed impulse restitution distribution compared to fitted distribution for fresh snow.

Fig. 8. Observed impulse-angle distribution compared to fitted distribution for compact snow.

Fig. 9. Observed impulse restitution distribution compared to predicted distribution for compact snow.

Because of the construction of this splash function, it can be extrapolated to angles outside of the range of observation. Figure 10 shows the pdf of the impulse angle for a range of impact angles. Note that this is symmetrical for a vertical impact angle. Figure 11 shows the probability of observing a ricocheting particle as a function of impact speed and angle. At low speeds, the elliptical censoring hypothesis dominates, that is even if a particle bounces without losing energy it will not bounce high enough to be observed. As the impact speed increases for shallow angles, the probability tends towards one. For steeper angles the probability tends towards the probability of a particle being scattered with a negative vertical velocity so that it is captured by the bed. Due to the small range of impact angles observed, it is not possible to say whether these predicted probabilities of absorption are correct. Experiments at steep impact angles would be very interesting for testing splash functions and could answer this question. If it turns out that these predicted probabilities are too low, it may be necessary to allow the impulse restitution to depend on the impact angle. This could be done by considering the effective mass ratio, which one would expect to be larger for steep collisions. Another possible modification would be to allow multiple collisions. That is, a downwardly scattered particle would be allowed to collide repeatedly until it scattered upwards or its energy had dropped below the observable level.

Fig. 10. Impulse-angle pdf for a range of impact angles.

Fig. 11. Observed ricochet probability as a function of impact angle and speed.

5. Ejected Particles

After a bed collision, additional particles may be ejected, and these were counted, though their velocities were not measured. However, the number of ejected particles is only well defined if it is understood to mean the number of particles with an energy greater than some observable critical value E_c (or some other function of velocity). Because the velocity of the ejected particles was not measured, only a limited analysis can be performed and it is hard to separate the limitations of the measurement process from the underlying probabilistic process.

Let f_e(v1; v₂; . . . ; u; v) be the pdf for the ejection velocity given the impact and ricochet velocities, where as before u is the impact velocity, v the ricochet velocity and v_i the velocities of the ejected particles. Since the particles are not distinguished, we take υ₁ > v₂ > ⋅ ⋅ ⋅. If the probability of observing an ejected particle is given by p_o(v) then the probability of observing at least n ejected particles is

(15)

The primary constraint that f must satisfy is of course that the sum of the energies of the ejected particles and the ricocheted particle is less than or equal to the energy of the impacting particle. Thus

(16)

We propose a very simple model constructed to satisfy this constraint that contains only one parameter λ, the fraction of energy lost during each collision. We consider the entire splash process as a series of binary collisions. Suppose that after the kth collision the remaining energy is E_k. Then during a collision (1 = λ)E_k is assumed to be dissipated or radiated into the bed. Some amount λX_iE_k, with X_i uniformly distributed between 0 and 1, is given to particle k + 1 so that and the remaining energy is passed forward to the next collision so that Any particle with an energy greater than is assumed to be observed. Particles with less energy are ignored.

The mean energy for the kth collision will be (λ/2)_k E₀, so the mean velocity of the kth ejected particle will be (λ/2)_k E₀/2. The mean number of ejected particles will thus be roughly log(2E_c/E₀)/log(λ/2). Unfortunately this approach is ill conditioned when the velocity of the ejected particles is not observed. Instead we consider an exponential distribution designed to have similar properties that combines the observation process with the ejection process. Let the probability that at least n ejected particles are observed be

(17)

where the ejection energy scale E_e combines material properties and the observation accuracy E_c, is a dimensionless parameter and E = m(u² - υ²)/2. With this distribution, the number of observed particles should scale as 1/γlog(E/E_e), that is it increases logarithmically not linearly. The parameters were found using maximum likelihood estimation and were:

where υ_e is the ejection velocity scale defined by The excellent agreement between the observed distributions and the data is shown in Figure 12 and Figure 13. Such good agreement is of course not surprising when we are fitting two parameters and only comparing a few categories. Another test, shown in Figure 14, is to compare the probability of at least one ejection with the theory. This shows that energy in all the observed collisions is sufficiently low that the probability increases approximately linearly over the observed range and it is not possible to test its exponential dependence. Though the data are not sufficient to validate the theory, they are at least not in contradiction.

Fig. 12. Predicted and observed number of ejected particles for compact snow.

Fig. 13. Predicted and observed number of ejected particles for fresh snow.

Fig. 14. Probability of at least one ejection as a function of (available energy per unit mass converted to velocity). The line is the prediction and the circles are observations.

6. Conclusions

This paper has developed a splash function that is based on a physical model of the grain–bed collision process. It is well defined for all impact velocities, respects appropriate symmetries, and ascribes zero probability to states that create energy. The theory is based on the impulse-angle distributionwhich is related to the roughness of the surface. It would be interesting to measure this directly in future experiments and to compare this with the implied impulse angle pdf. The impulse restitution is related to the stiffness of the contacts between the grains in the bed and the material restitution properties. Though there are not enough data for a rigorous validation of the splash function, the available data are in agreement. Clearly there is much scope for future work. Such experiments should take great care to understand the limitations of the measurement process if bias in the results is to be avoided. A fully automated experiment that could detect of the order O(105) collisions in three dimensions would be immensely useful. Until such data are available, so as to accurately study the correlations between impact and ejection velocities, no splash function can be considered established.