## Abstract

We study the dynamics of systems with deterministic trajectories randomly forced by instantaneous discontinuous jumps occurring according to two different compound Poisson processes. One process, with constant frequency, causes instantaneous positive random increments, whereas the second process has a state-dependent frequency and describes negative jumps that force the system to restart from zero (renewal jumps). We obtain the probability distributions of the state variable and the magnitude and intertimes of the jumps to zero. This modelling framework is used to describe snow-depth dynamics on mountain hillsides, where the positive jumps represent snowfall events, whereas the jumps to zero describe avalanches. The probability distributions of snow depth, together with the statistics of avalanche magnitude and occurrence, are used to explain the correlation between avalanche occurrence and snowfall as a function of hydrologic, terrain slope and aspect parameters. This information is synthesized into a ‘prediction entropy’ function that gives the level of confidence of avalanche occurrence prediction in relation to terrain properties.

## 1. Introduction

Many natural systems are characterized by slow dynamics, which can be described as deterministic functions and random perturbations that occur at much shorter time scales (Ridolfi *et al.* 2011). When represented at the time scale of the slow dynamics, such perturbations can be assumed to be instantaneous, causing sudden jumps in the system. For example, daily soil moisture decreases slowly because of seepage and transpiration, and increases rapidly during rainfall events (Rodriguez-Iturbe *et al.* 1999); leaching of soil pollutants and groundwater recharge are similarly affected by rainfall events (Botter *et al.* 2007; Suweis *et al.* 2010; Manzoni *et al.* 2011). The effect of disruptive events, such as natural fires and floods, on above-ground vegetation biomass also occurs at a much shorter time scale than vegetation growth (D’Odorico *et al.* 2006; Perona *et al.* 2009; Crouzy & Perona 2012). Similar type of processes with jumps have also been used in mathematical finance (Cont & Tankov 2006, 2009) and biology (Tuckwell 1988).

Some natural systems can also be affected by a combination of different jumps. Such processes however have been less investigated. An important example of interest here is the dynamics of snow depth on mountain hillsides, which is driven by positive jumps associated with snowfalls (Perona *et al.* 2007), as well as negative jumps caused by avalanches occurring as a function of the actual snow depth (see Schweizer *et al.* 2003 for a review). In this paper, we focus on processes that can be described by a single-state variable, whose trajectory is driven by a deterministic drift and is interrupted by instantaneous random jumps owing to two different mechanisms. We assume that the jumps occur according to two independent compound Poisson processes, one with constant frequency and the other with state-dependent frequency (Daly & Porporato 2006). The first type of jump causes sudden increases of the state variable, whereas the second type of event abruptly restarts the system from a given state. The occurrence of this second type of jump is thus associated with a renewal process (Cox & Miller 1965).

Our theoretical framework is then used to develop a model for snow depth dynamics, which evolves as a function of snow mantle compaction and melting, snowfall and avalanche events. Particular emphasis is devoted to the statistics of avalanche occurrence and magnitude in relation to local terrain features. This is an area of intense research in the geophysical sciences (Harbitz *et al*. 2001; Ancey *et al.* 2003; Schweizer *et al.* 2009), which has traditionally focused on physically based models. Because of the inherent stochasticity driving the avalanche formation and triggering (Eckert *et al.* 2007; Cappabianca *et al.* 2008), such efforts have only been partially successful. In particular, some important yet unexplored questions regard: (i) the probability distribution of avalanche size occurrence in relation to terrain slope and snow compaction processes, (ii) the dependence of the return time of avalanche on the return time of preceding snowfall, and (iii) the identification of those slopes on which prediction of avalanche detachment and size is most uncertain. The model presented here starts addressing some of these issues.

## 2. Model description

We study processes that can be described by a single-state variable, *h*, whose dynamics are driven by the Langevin-type equation
2.1
In equation (2.1), the function *ρ*(*h*) is a deterministic drift assumed to be positive and with the linear form *ρ*(*h*)=*ρ*_{0}+*ρ*_{1} *h* (*ρ*_{0}≥0,*ρ*_{1}≥0), such that it causes *h* to decrease in time. The sequences {*t*_{i}} and {*t*_{j}} (with *i* and *j* positive integer numbers) are the instants of occurrence of the *i*th and, respectively, *j*th event of two Poisson processes, the first of which, *N*_{i}(*t*) (*t*≥0), with constant intensity *λ*, and the other, *N*_{j}(*t*) (*t*≥0), with events occurring at a state-dependent rate *ν*(*h*)=*ν*_{0}+*ν*_{1} *h* (*ν*_{0}≥0,*ν*_{1}≥0). *δ*(⋅) denotes the Dirac delta distribution. The random positive jumps {*y*_{i}}, associated with the Poisson process *N*_{i}(*t*), are independent and identically distributed as *g*(*y*)=*γ* *e*^{−γy}. The magnitudes of jumps occurring according to the state-dependent process, *N*_{j}(*t*), are such to reset the state variable *h* to zero (renewal events), and therefore may be considered to be deterministic. This means that the function *h*(*t*) in the last term of (2.1) is computed just before the occurrence of the event *t*_{j}. The ensuing jump distribution of the renewal events is a result of the entire dynamics and will be computed in what follows.

As a result of (2.1), the domain of the state variable is *h*≥0, and the dynamics of *h* are driven by two opposite actions: the positive jumps, which move the process away from zero, and the drift and the state-dependent jumps that move *h* towards the bound at zero.

The probability density function (PDF) of *h*, *p*(*h*,*t*), consists of a continuous part with density *p*^{c}(*h*,*t*), and an atom of finite probability, *p*^{at}_{0}(*t*), at *h*=0. The forward Chapman–Kolmogorov (or master) equation for *p*^{c}(*h*,*t*) can be written as (Cox & Miller 1965; Gardiner 2004)
2.2
while the equation describing the dynamics of the atom of probability *p*^{at}_{0}(*t*) is
2.3
The solution of equations (2.2) and (2.3) needs to be normalized in order to satisfy the unit probability condition
2.4

In steady-state conditions, the value of the atom of probability can be calculated from equation (2.3) as
2.5
The exponential jump distribution, which has often been used in applications (Takács 1960; Weiss 1977; Rodriguez-Iturbe *et al.* 1999), allows us to transform the integro-differential equation (2.2) into a simpler form (Cox & Miller 1965). By multiplying all terms of equation (2.2) in steady state by e^{γh} and then differentiating with respect to *h*, one obtains the second-order ordinary differential equation
2.6
the general solution of which is (Polyanin & Zaitsev 2003)
2.7
where *U*(⋅,⋅,⋅) and *M*(⋅,⋅,⋅) are confluent hypergeometric functions (Abramowitz & Stegun 1965); because *M*(⋅,⋅,⋅) diverges when *h* tends to infinity, *C*_{2}=0. The atom of probability in *h*=0 can be calculated using equation (2.5), imposing the condition in equation (2.4) to calculate the constant *C*_{1}.

In the special case, *ρ*_{1}=0, *p*^{c}(*h*) reads
2.8
When *ρ*_{1}=0 and *ν*_{1}=0, i.e. when the drift is constant and when events of both Poisson processes occur at a constant rate, the model is a modified version of the Takács process (Takács 1960). In this case, the complete solution of equations (2.5) and (2.6) is
2.9
where
A comparison of PDF *p*^{c}(*h*) is shown in figure 1 for different combinations of parameters *ρ*_{0},*ρ*_{1},*ν*_{0},*ν*_{1}.

### (a) Magnitude of the jumps to zero

The continuous part of the distribution of the variable *h* at steady state, *p*^{c}(*h*), can be used to derive the PDF, *p*_{r}(*h*), of the magnitudes of the jumps associated with the state-dependent renewal events that reset *h* to zero. As discussed by Daly & Porporato (2007), at steady state, *p*_{r}(*h*) can be obtained as the conditional probability of being at the state *h* and of having an independent renewal event from that state in the time interval d*t*, renormalized to the total probability of having a renewal event in the time interval d*t*,
2.10
When *ν*_{1}=0, jumps to zero occur according to a Poisson process with constant rate *ν*_{0} and thus jumps occur from any value of *h* with the same probability in any time interval. Examples of *p*^{c}(*h*) and *p*_{r}(*h*) are shown in figure 2. Figure 2*a* shows that the distribution of *h* can have a mode larger than zero, depending on the role of the state-dependent dynamics. As *ν*_{1} increases, i.e. the frequency of jumping to zero becomes larger, becomes larger and *p*^{c}(*h*) presents an increasing peak at *h*=0. As the influence of the state-dependent term increases, the descending part of *p*_{r}(*h*) gradually approximates the exponential PDF of jumps (figure 2*b*). This suggests that jumps to zero tend to occur right after the positive jumps.

### (b) Times between jumps to zero

In some applications, it is important to be able to quantify the statistics of the intertime, *τ*, between renewal events. The probability distribution function of *τ*, *p*_{τ}(*τ*), can be computed by defining an ancillary function indicated with the subscript *ν*, *p*_{ν}(*h*,*τ*), accounting for a continuous part *p*^{c}_{ν}(*h*,*τ*) and an atom of finite probability *p*^{at}_{ν}(*τ*), and driven by the equations (Daly & Porporato 2006, 2007)
2.11
and
2.12

In order to obtain the PDF of the renewal event intertimes, we consider a process where at *τ*=0, the system is subject to a renewal and find the distribution of the time of occurrence of the next renewal. This corresponds to the boundary conditions
2.13
and
2.14

Equations (2.11) and (2.12) are similar to equations (2.2) and (2.3) with the exception that equation (2.12) does not contain the integral term, , related to the input of trajectories into the atom due to renewal events. Therefore, as shown in figure 3, *p*_{ν}(*h*,*τ*) can be associated with a process identical to that described by equation (2.1) with the only difference that, whenever a renewal event occurs, the trajectory is not re-injected at *h*=0 but the process stops. Accordingly, *p*(*h*,*t*) describes the probability distribution of an ensemble of trajectories each driven by equation (2.1) (figure 3*a*). In contrast, in the case of equations (2.11) and (2.12), each trajectory stops when the first renewal event occurs (figure 3*b*), so that the area under *p*_{ν}(*h*,*τ*) is no longer constant and equal to one, but decreases as *τ* increases. This area thus represents the fraction of trajectories that started at *h*=0 and did not experience a jump to zero up to the time *τ*. In other words, the area under *p*_{ν}(*h*,*τ*) defines the survivor function of the process, which decreases in time because of the loss of probability when a trajectory is lost as a result of a renewal (Daly & Porporato 2007). The PDF of the time of the first jump to zero when starting at *h*=0 can hence be calculated as
2.15

The term containing the atom accounts for transitions from the atom to the continuous part and vice versa.

Because the system of equations (2.11) and (2.12) is difficult to solve, the solution of the problem can be alternatively obtained by finding the Laplace transform of *p*^{c}_{ν}(*h*,*τ*), i.e.
2.16
and calculating its derivative with respect to *τ* in *s*=0, i.e.
2.17
This function will be analysed in detail for the specific case of the avalanche model.

## 3. Stochastic dynamics of snow avalanche occurrence

The stochastic dynamics previously analysed can be used to model the process of snow avalanche occurrence described at a time scale that is longer than the actual duration of snowfall or avalanche events. In this case, the snow mantle depth decreases slowly under a quasi-deterministic dynamics during no-snow days and increases suddenly during snow events until an avalanche resets the process.

The return-time statistics of avalanche occurrence and size are known to differ from those of preceding precipitation (Ancey *et al.* 2003; Eckert *et al.* 2007). In particular, the slope of the hillside determines the equilibrium between the weight of the snow mass and the friction forces, while snowfall events are responsible for the increased local load. Thus, on high-slope terrains, one may expect the risk of detachment to be inherently high because of mere stability reasons, and avalanche occurrence to relate more closely to preceding precipitation events.

Changes in snow depth due to snow metamorphism and consequent compaction play a less obvious role. Snow layers are generally influenced by post precipitation meteorological conditions, as well as the existing snow height, which causes snow metamorphism, melting and refreezing processes (see Schweizer *et al.* 2003 for a review). Whether changes of the microscopic structure of snow favour, or not, snow layer stability is not yet well understood, and this remains a source of uncertainty in current models of snow avalanche prediction.

### (a) Avalanche model and solution

In the following, we assume snowfall to occur as an instantaneous Poisson process of constant rate *λ*; snowfall depths are assumed to be exponentially distributed with average magnitude *α*=1/*γ*. Avalanches are assumed to take place as a Poisson process with state-dependent rate *ν*(*h*)=*ν*_{1}*h*, where *ν*_{1} is a parameter accounting for topographic factors such as the terrain slope. For a given snow depth, higher terrain slope can be assigned a higher probability of avalanche occurrence, and will therefore be characterized by a larger value of the parameter *ν*_{1}. Snow compaction due to snow metamorphism and partial melting processes determine a reduction of *h* that we approximate as a deterministic decay of the type d*h*/d*t*=−*ρ*_{1}*h*. Although this assumption is clearly a simplification of the real process, where snow depth can actually vanish at a certain time, it shows a certain similarity with the snow cover evolution of physically based simulated processes (Schweizer *et al.* 2009). Accordingly, the state-dependent model of equation (2.1) with *ρ*_{0}=*ν*_{0}=0 can be considered a minimalist representation of the process of snow avalanche occurrence at a point.

With these specifications, the PDF of the snow depths can be derived from equation (2.7), obtaining
3.1
with a normalization constant
3.2
where *A*=*γλ*/(*γρ*_{1}+*ν*_{1}), *B*=*γλ*/(*γρ*_{1}+*ν*_{1}), , *D*=*ν*_{1}/(*γρ*_{1}+*ν*_{1}). The functions *Γ*(⋅) and _{2}*F*_{1}(⋅,⋅,⋅,⋅) are the gamma and the hypergeometric _{2}*F*_{1} functions, respectively (Abramowitz & Stegun 1965).

Once *p*^{c}(*h*) is determined, the atom at zero can be calculated by putting the rate of jumps to zero, , equal to the rate of jumps from zero, , yielding
3.3
The distribution of the avalanche events is readily obtained from (2.10) with *ν*(*h*)=*ν*_{1}*h*.

The exact PDF of the intertimes between avalanche events, which can be computed as in appendix A, reads 3.4

We show in figure 4*a* some examples of PDFs *p*_{τ}(*τ*) for the intertimes between renewal events. In figure 4*b*, we plot how the mode of *p*_{τ}(*τ*) changes as a function of the two parameters *ρ*_{1} and *ν*_{1}. Given the high skewness of such a distribution, the mode is more indicative of the typical frequency of renewal events, and more relevant also for practical applications. The mean of the density distribution *p*_{τ}(*τ*) can be analytically computed from the moment-generating function obtained by Laplace transforming equation (3.4). Because the corresponding expression for is rather cumbersome, we omit giving its explicit expression here; however, it will later be used to compute the mean frequency of avalanches, .

### (b) Avalanche versus snowfall return period

The PDFs obtained before can be used to study how snow depth and times between avalanche occurrences depend on the statistics of heavy precipitation, thereby complementing the existing studies on the statistical link between precipitation events and avalanche detachment (Harbitz *et al*. 2001; Schweizer *et al.* 2009). We start by computing the return period (*R*(*h*)) for snowfalls and avalanche events above a given snow depth *h*. The return period, defined as the mean number of days between an event *e* (either snowfall *s* or avalanche *a*) of depth at least *h*, can be computed as *R*_{e}(*h*)=1/(1−*P*_{e}(*h*)) (Brutsaert 2005), where *P*_{e}(*h*) is the non-exceedance probability of the daily maximum of *e*. The non-exceedance probability can be computed exactly for snowfalls, using the properties of the related homogeneous Poisson process with a constant frequency *λ* and the peak over threshold theory (Cramer & Leadbetter 1967). Accordingly, starting from the non-exceedance probability of a single snowfall of depth *h*, *G*_{h}(*h*)=1−*e*^{−γh}, the non-exceedance probability of the maximum of snowfalls over a period of 1 year, *T*=1 year, is , so that the return period for snowfalls of depth *h* is *R*_{s}(*h*)=1/(1−*P*_{s}(*h*)).

Because avalanches do not occur as a homogeneous Poisson process, their return period is more complicated to compute. It suffices here to compute it approximately referring to a Poisson process of mean rate (or equivalently ), where *p*^{c}(*h*) is given by equation (3.1). Thus, from the non-exceedance probability of a single avalanche *P*_{r}(*h*), which can be calculated by integrating equation (2.10), the approximate non-exceedance probability of the annual maximum of an avalanche of depth *h* is . Finally, the (approximate) return period for avalanches is again evaluated as *R*_{a}(*h*)=1/(1−*P*_{a}(*h*)). Figure 5 shows the results of this analysis for two different terrain orientations and characteristic slopes. Terrain orientation is assumed to influence the drift coefficient *ρ*_{1}, whereas the slope affects *ν*_{1}. High slope parameter values mask the effect of the compaction parameter, and the return period for snowfall and avalanche processes coincide (figure 5*a* and figure 5*b*), suggesting also a high correlation between the events. As the slope parameter decreases, the return period becomes clearly different and governed by the parameter *ρ*_{1}, which is the rate of compaction and shows that terrain orientation has an effect on the dynamics. If the rate of compaction is high enough (figure 5*c*), then enough memory is introduced in the stochastic process to determine an avalanche return period greater than that of an equal snowfall, and the two processes will be less correlated. As the effect of the drift diminishes (figure 5*d*)), the return period of the two curves may eventually cross, with that of avalanches becoming shorter than that of an equal precipitation. Although surprising, this result is indeed logical and related to the nature of the stochastic process, which may allow avalanches that are larger or smaller than the snowfall events, depending on the sequence of multiple accumulation and compaction of snow before an avalanche.

### (c) Avalanche prediction entropy

Depending on the process parameters, avalanches may be more or less correlated with the preceding snow event. The degree of de-correlation existing between the magnitude of the last snowfall and the successive avalanche event decreases with slope, and the avalanche size distribution tends to that of snowfalls as slope increases (figure 2*b*). Figure 6*a* show the increase in the correlation of numerically generated avalanche size versus last snowfall events. The decorrelation function, *χ*, defined here as 1 minus the crosscorrelation between snowfall depth and the preceding snow depth, is plotted in figure 6*a* as a function of the slope parameter *ν*_{1}, showing how avalanche events tend to become decorrelated from snowfalls as the slope of the hillside decreases.

To quantify the prediction uncertainty, which has both descriptive and practical purposes (Barbolini & Savi 2001; Chernous & Fedorenko 2001), we introduce also a prediction entropy (*H*) function, which measures the degree of confidence with which avalanche occurrence may be predicted depending on process parameters. This measure could, for instance, be integrated in geographic information system platforms, similar to Cappabianca *et al.* (2008), once the parameters are linked to local terrain properties. Thus, the avalanche predictability will be investigated as a function of mechanical (e.g. slope coefficient *ν*_{1}) and metamorphism processes (e.g. the compaction coefficient *ρ*_{1}). We first define the avalanche versus snowfall frequency ratio, which is the ratio between the average number of avalanches and that of snowfalls in a given period of time. Because at most one avalanche event can take place between snowfalls, such a ratio is also the probability of observing an avalanche between two successive precipitation events,
3.5

In other words, *P*_{a} provides a measure of ‘synchronization’ between avalanche and preceding snowfall sizes. Moreover, because of the form of the state-dependent frequency, *ν*(*h*)=*ν*_{1}*h* increases with the slope and tends to 1 when the distribution of avalanche intertimes tends to that of snowfalls. It emerges that predicting avalanche occurrence following snowfall events has high reliability, both on steep and gentle (though bigger than the critical slope, approx. 20^{°}÷25^{°}) terrain because of the high correlation and the low occurrence probability, respectively, while at intermediate slopes, these two factors concur to increase the uncertainty.

The prediction uncertainty can be precisely defined for the binary random variable representing the occurrence of an avalanche between snow events, with probability *P*_{a} as the Shannon entropy of a binary variable (Cover & Thomas 2006, p. 15)
3.6

Because *P*_{a} depends on the parameter *ν*_{1}, *H* measures the level on uncertainty in predicting avalanche occurrence at different slopes. This function is plotted in figure 6*b* for different values of the drift parameter *ρ*_{1}, which shows the effect associated with the metamorphism dynamics that may occur in changing terrain aspects such as elevation and exposition. Notice that *H* tends to a limiting function as , i.e. when no drift is present. In general, the fact that *H*=1 at intermediate slopes indicates a minimum in the reliability of predicting avalanche occurrence. This is caused by the comparable role played by deterministic compaction (ascribable to snow metamorphism) and the state-dependent stochastic term (ascribable to mechanical stability) at intermediate slope, whereas either one or other process prevails at low and high slopes.

In summary, the measure *H* shows that this simple model explains why statistics of avalanche occurrence coincide with those of intense precipitation only in the limit of very steep slopes. Within this picture, events occurring on high slopes are well determined by preceding intense precipitation and are likely powder-type avalanches. Snow metamorphism occurring between successive precipitation events, although increasing the stability of the snow mantle, also contributes to decorrelating the occurrence of avalanches from preceding precipitation. Thus, at low slopes (yet bigger than the critical angle), the higher stability of the snow mantle makes the prediction reliability increase again, while avalanches become more slab type.

## 4. Conclusions

In this paper, we studied a family of stochastic processes characterized by the interaction of compound Poisson processes, and applied it to describe the statistics of snow avalanche occurrence at a point over long time scales. The interplay between the drift and the state-dependent jumps was found to strongly affect the resulting PDF of the state variable (figures 1 and 2). The main role of the state-dependent events is to modify the probability distribution function of both the sizes and the intertimes between these events, which thus acquire a non-exponential distribution.

Despite its simplicity, the snow avalanche occurrence model contains the main basic ingredients of the real process, allowing us to explain the deviation of observed snow avalanche statistics from those of preceding precipitation (Schweizer *et al.* 2009). As a result, avalanche sizes in the detachment zone were shown to correlate to snowfalls only on high slopes, where instability due to load increase dominates and snow compaction plays a minor role. On the contrary, snow metamorphism was found to become important at lower slopes (depending on orientation), while snowfall and avalanches tend to decorrelate because of snowfall accumulation before avalanche events. At intermediate slopes, snow accumulation by multiple snow events and compaction by metamorphisms have comparable importance and this results in a maximum entropy prediction uncertainty of avalanche, given the occurrence on a snowfall event.

## Acknowledgements

E.D. thanks the Institute of Environmental Engineering (IfU), ETH Zurich, CH, for hosting him in July–August 2010, and the Australian Academy of Science for funding this visit within the IRSES programme. P.P. wishes to thank the European Commission for financing the Programme ‘PEOPLE’ Call ID FP7-PEOPLE-IRSES-2008 Proposal N230845 (WARECALC), and the Swiss National Science Foundation for financially supporting the Group AHEAD (grant no. PP00P2-128545/1).

## Appendix A. Probability density function of the time between avalanche occurrence

In this appendix, we solve equations (2.11) and (2.12). When *ρ*_{0}=*ν*_{0}=0, equation (2.12) can be simplified to
A 1
Because after an avalanche *h*=0, one obtains
A 2
By inserting (A2) in (2.11), we obtain, after Laplace transformation,
A 3

In the Laplace representation, the boundary condition of the equation becomes simply *p**(*s*,*τ*=0)=0. Now, by introducing the parametrization
A 4
and
A 5
One can transform equation (A3) into an ordinary differential equation, which can easily be solved. Reverting to the original variables, one then uses equation (2.17) in order to obtain the distribution of the renewal intertimes. The final result reads as equation (3.4).

- Received July 4, 2012.
- Accepted September 3, 2012.

- This journal is © 2012 The Royal Society