A non-Markovian model of rill erosion

We introduce a new model for rill erosion. We start with a network similar to that in the Discrete Web and instantiate a dynamics which makes the process highly non-Markovian. The behavior of nodes in the streams is similar to the behavior of Polya urns with time-dependent input. In this paper we use a combination of rigorous arguments and simulation results to show that the model exhibits many properties of rill erosion; in particular, nodes which are deeper in the network tend to switch less quickly.


Reinforcement and Rill Erosion
Stochastic processes with reinforcement are inherently non-Markovian and therefore may model some real phenomena more accurately than can their Markovian counterparts. Reinforcement is a mechanism that provides a bias to a system, making it more likely to occupy states the more often those states are visited. Some well-studied examples include variations on the urn of Pólya (the original introduced in [4] and this and subsequent models studied, for example, in [1] and [9]) and reinforced random walks ( [3,15]). The infinite memory exhibited in these examples can force a system to spend most (or almost all) of its time in a small subset of its state space. Many natural phenomena exhibit similar behavior; for instance, the overall pattern of erosion on a hillslope is relatively stable once it is established, although small details of the pattern may change frequently and catastrophes that permanently alter it may occasionally occur.
We investigate a discrete time, infinite-memory random process defined on the nodes and edges of an oriented diagonal lattice ( Figure 1) that we propose as a simple model of hillslope erosion. The lattice starts out smooth in the sense that it has no edges initially, but it sprouts edges everywhere the instant the process starts, much as rain can start soil erosion everywhere on a hillslope at once. Edges may connect an interior node to two, one, or neither of the two nodes directly above it. Exactly one edge descends from each interior node, and it points either left or right. At every node and at every time step a simple two parameter reinforcing law, based on the entire history of the network above a given interior node, randomly determines the direction of the nodes descending edge and then is updated. Obvious modifications of these statements apply to nodes at the top or bottom (if one exists) of the lattice.
The current pattern of connections among nodes represents the present state of the process, and the patterns stability -measured by the tendency of the same state, or one similar to it, to occur on subsequent iterations of the processrepresents the patterns strength as a memory. The degree of reinforcement is set by tuning two parameters, r and α. At any given moment the current pattern is a collection of dendritic networks that appears similar to drainage networks found in nature; indeed, lattice models have often been used to investigate the morphology of natural drainage networks (e.g. [17]). We focus on the surficial dynamics of rill networks [10], rather than their morphology. Put in terms of erosion, we are more interested in the process of erosion than we are in the result.
The analogy between our model and erosion, specifically rill erosion, is straightforward: r can be interpreted as a rainfall rate (or equivalently, as the rate of sediment generation) and α −1 as the resistance of soil to erosion, while the reinforcement dynamics correspond to the overland flow of water and sediment down a hill. Rills are small, ephemeral channels that transport sediment down hillslopes when it rains [18]. They form when rainfall and runoff dislodge particles from the soil surface and transport them along flow paths governed by variations in the surface roughness of soils and the soil's ability to resist erosion. Flow depths in rills are typically on the order of a few centimeters or less, while the longest channels in rill networks can be several meters long. Processes affecting rill erosion take place over timescales ranging from milliseconds to hours.
The topology of rill networks is relatively unstable when compared to larger scale natural drainage systems (of which rills may be a part) like gulley systems and river basins. Rill networks are most unstable at their tops where boundaries between rills and inter-rill areas are not well defined and shift often, but connectivity can change downhill as well, usually at a slower rate than uphill. Some rills grow throughout a rainfall event, others are filled by sediment and disappear, still others alternate. A detailed description of rill erosion 1) must account for complicated interactions among rainfall, soil properties, and topography [19], and 2) often depends on obtaining a set of physical parameters that are difficult to measure.
Despite the high degree of complexity of rill erosion at small scales, at macro-scopic scales it is principally determined by particle detachment, overland flow, and sediment transport ( [16]). In turn, each of flow, detachment, and transport depends critically on the rate of rainfall and the soils resistance to erosion. It is not completely surprising that our simple two parameter model exhibits some important elements of the macroscopic behavior of rills formation. In fact, similar to rill erosion, each node in the model network switches direction infinitely many times but the switching rate depends on position up or down hill. Furthermore, floods that carry unusually large amounts of water and catastrophes that significantly alter the flow pattern occur occasionally in the model, as they do in nature.

Definition of the Model
Consider the vertices in the even sub-lattice of Z 2 which have second coordinate non-positive. That is, the set x + y is even and y ≤ 0} and edges be a node with left parent w 1 = (x − 1, y + 1), right parent w 2 = (x + 1, y + 1) (for those nodes with second coordinate 0, parents will not exist), and with left child (x − 1, y − 1), right child (x + 1, y − 1). We will use the term depth k to refer to those nodes with y coordinate equal to 1 − k. Conversely, for any node v the term depth(v) will denote the numerical value of the depth of v.
First we describe the algorithm for the behavior of v heuristically. At the end of the 0 th second, v receives I v (0) = r units of rain (but does nothing else). During the n th second, for n ≥ 1, the following sequence occurs: 1. v flips a coin, heads-biased with probability P L v (n), which reflects T L v (n), the total "sediment" load v has sent to its left child by time n.
2. If this coin shows heads (tails), v sends its current input of sediment I v (n − 1) to its left (right) child. v adds this number to the total sediment, , it has sent to the left (right) for all time.
3. v receives sediment load from 0,1, or 2 parents and receives r units of rain. Call the sum of these two I v (n). Increment time and return to step 1.
The evolution of the node's behavior depends on two parameters: the rainfall rate r > 0 and a term α −1 > 0 that resists change. To make this precise, we make several definitions. We start by initializing variables. For each v ∈ Z 2 even , let even , we define a Bernoulli variable, D L v (n), the biased coin, that is conditionally independent from vertex to vertex given the variables We create the bias for the next coin: where η = α/r compares the effect of the rain to the system's inherent resistance to change. In this paper, we shall always take r = 1 so that η = α. Last we define the input w1 (n)) + I w2 (n − 1)D L w2 (n) + r otherwise and the filtration Figure 1 for an illustration of the process at node v. (2), ...) the sequence of directions that node v chooses (for example (L,R,L,...)). At the end of time t, after all nodes have sent loads to their children, we may update certain edge variables. Define a sequence of edge configurations {ω n } n≥0 , where for each n, ω n is a map from E 2 even → {0, 1}, using the following rule. If the node v = (x, y) has d v (n) = L then let If, on the other hand, d v (n) = R then let Figure 1 for a realization of ω n .
We say that nodes v and w are connected at time n if there exists a path of distinct adjacent edges e 1 , .., e m with ω n (e i ) = 1 for all i so that e 1 connects v

Regimes for η
The parameter η plays an important role in the behavior of the model. For a fixed node v (at depth k) we have that for all n ≥ 1 This indicates that when η is small the node v chooses a direction at time 0 and has a high probability of sticking to this direction for most values of n ≥ 1.
Since this is true for each node v, the evolution of {ω n } is somewhat simple. In the limit as η → 0, each node picks a direction and stays with that direction for all time. That is, for each n ≥ 1, and for each finite subset E ⊂ E 2 even , and the dynamics has no effect on the configuration in any finite subset of Z 2 even . The configurations at any moment are the same as those in the discrete web ([2], [7]).
In the other direction, as η → ∞ each node "forgets" its history. That is, for each node v, the conditional probability given F n that it chooses left at time n + 1 is given in (1), and the limit of this quantity is 1/2. By symmetry, Therefore the variables in any finite subset of {d v (n) : n ≥ 0} converge in distribution to i.i.d. Bernoulli(1/2) variables. Furthermore, the variables in any finite subset of {d v (n) : n ≥ 0, v ∈ Z 2 even } converge in distribution to i.i.d. Bernoulli(1/2) variables. Intuitively this holds because distinct nodes only interact with each other through their input and output loads, and both of these are eventually dominated by large η. These statements indicate that when η is large, the dynamics of our erosion model are similar to those in a network in which each node flips a fair coin at each time n, independently from site to site and from time to time, to determine in which direction to send sediment. Thus the configurations {ω n } resemble those taken from the dynamical discrete web ( [11], [7]).
Given the relation both these extreme cases have to the variables {ω n }, it is natural to view the present model (with 0 < η < ∞) as an interpolation between the discrete web and the dynamical discrete web. Indeed, for each fixed n, the distribution of ω n is the same as that in the case of η = 0 at time n = 0 or that in the case of η → ∞ at any time n. (In both cases, all directions are chosen by independent fair coin flips).
As we shall see in section 3.2.2, the model with 0 < η < ∞ can be likened to the case η → ∞ in the following way. Each level k is associated to a measure θ k (defined in eq. (8)). Each node at level k samples (non-independently) from this measure a value p v . For any sequence n 1 , n 2 , ..., n m of times and for any sequence x 1 , ..., x m of elements from the set {L, R}, Because of this fact, we may view the model (for large time) as one in which each node fixes a Bernoulli parameter p v and flips a p v -biased coin independently each second n (but not independently from site to site) to determine the direction in which to move sediment.

Outline of the Paper
In section 2 we discuss the (relatively simple) behavior of nodes at depth 1.
Since these nodes receive constant input load over time, we can use the well known model of Pólya's Urn to analyze their output. In section 3 we discuss the more complicated behavior of nodes at depth at least 2. Here we make use of results of Pemantle [14] for the time-dependent Pólya Urn. We look more closely at properties of the input load, of the output load, and of the dynamics of these lower-depth nodes.

Top Level
Since our top level nodes are equivalent to the model of Pólya's urn, we recall basic facts of Pólya's model Start with an urn containing R 0 red balls and B 0 black balls and draw one ball from the urn. Return this ball to the urn, along with another ball of the same color. After this round there are R 1 red balls in the urn and B 1 black balls in the urn, with either R 0 = R 1 or B 0 = B 1 . Repeat this process infinitely many times, creating sequences {R n } n≥0 and {B n } n≥0 so that for each n, It is well known that the fraction F R n = Rn Rn+Bn has an almost sure limit and that this limit is distributed as β(R 0 , B 0 ) (see e.g. [8]).
Let v be a node at depth 1. At the beginning of each second, v receives an amount of sediment equal to 1 and this input load amount does not change with time. The node sends this load either to the right or left, depending on the bias rule in (1). We are interested in the fraction of total load the node sends left (right) up to time n. To this end, define the load fraction Theorem 2.1. The quantities LF L v (n) and LF R v (n) have limits as n → ∞. These limits are random: they are distributed as β(η, η).
Proof. We will indicate the proof only for the case LF L v (n). An easy calculation shows that P L v (n) is a martingale w.r.t. {F n } and, since it is bounded for all n, it has an almost sure limit. Solving for the limiting distribution is similar to solving for the related quantity in the standard Pólya urn model. See, for instance, [5]. This gives Note that the limiting distribution in Theorem 2.1 is supported on [0,1] and has no atoms. This implies that with probability 1, the node v switches states (L,R) infinitely often and that neither of these states is transient. This is quite unlike the "sticking" associated to the dynamics in the η → 0 limit (refer to (2)). The distribution from the above theorem for different values of η is pictured in Figure 2. For 0 < η < 1 the limiting load fraction has a bimodal distribution, and for η > 1 the distribution is unimodal, symmetric about 1 2 . This means that when η is small, each node is likely to have a relatively strong preference for one direction and that when η is large, each node is likely to favor L and R somewhat equally. The case η = 1 gives a uniform distribution.
Here v is equally likely to have a strong directional preference as it is not to.

Lower Levels
The simplicity of behavior at the top level comes from the fact that each node has an input load which is constant w.r.t. time. This is not true at lower levels. Each node has an input load whose magnitude is non-trivially time dependent.
To make this more apparent, isolate an arbitrary node v at depth k. If at time t = n, v is not connected to either of its parents in ω n , then its input load is 1 unit (coming only from rain). If, on the other hand, v is connected to at least one of its parents, then its input load will be strictly greater than 1 unit. Therefore, the geometry of the connected components of ω n determines the behavior of each node. This relationship is complex for at least two reasons. First, not only does the geometry of the network influence node behavior, the node behavior in turn determines the future geometry of the network. In this sense, our system generates its own randomness. Second, the method by which this randomness arises involves propagation. The geometry of nodes at depth k − l at time m affects the behavior of nodes at depth k at time n if and only if m = n − l. In other words, it takes l seconds for the output load from depth k − l to reach nodes at depth k. In spite of these complications, we set out to analyze these lower level nodes.
The node v has an input load sequence . We are interested in analyzing the nature of this input sequence, the nature of the output sequence, and the relationship between the two. of η does not matter, as a consequence of Theorem 3.1). The simulation was conducted with periodic boundary conditions, with 10 6 nodes per row, and with 10 rows. Therefore, the histogram for depth k at time n = 300s should closely approximate the probability mass function of the distribution of the input load for depth k at time n = 300s. One notices a few things. First, the support of the distribution at depth k is integers in the interval [1, 1 2 k(k + 1)]. Next, the mass function appears to decrease from load value 1 to a local minimum at k − 1, to increase for a bit to a local maximum, and then to decrease to the edge of its support. About 1/4 of nodes are at the heads of rills, while the fraction of rills starting short of the top increases with depth. The "bump" in the load distribution to the right of the value k − 1 appears to travel to the right as depth increases. Looking at Figure 3, it is tempting to guess that the load distribution at a given level is a mixture of a distribution for loads that start at the top and one for loads that do not. Last, the different mass functions have several common values. For example, the probabilities for load values 1 to 4 are the same in each figure, and the probabilities for load values 1 to 6 are the same in the center and right figures.

Input Load
We present three structural theorems regarding the load distribution. The first gives basic information needed to make calculations, and the second gives us the value of the first moment of the distribution. The third discusses a limiting measure for the family of loads that do not originate at the top. Because of the simplicity of the first theorem, we state it without proof. c. The distribution of I v (n 0 ) is equal to the distribution of |C + w,n0 | for any node w with depth equal to min(n 0 , k). Therefore Theorem 3.2. Let v be a node at depth k. The mean of the load distribution is Proof. We prove by induction on k. For k = 1, the statement is trivial, so consider k > 1. Since the distribution of I v (n) is constant for n ≥ k, we assume n ≤ k. Let N v,k−1 be the number of nodes at level k−1 which send sediment to v at the end of time n − 1. This variable takes values in {0, 1, 2} with probabilities {1/4, 1/2, 1/4}, respectively. Call w 1 (w 2 ) the left (right) parent of v.
where to go from the second line to the third line, we use the fact that the variables I w1 (n − 1) and I w2 (n − 1) have the same distribution (see b. under Theorem 3.1). Theorem 3.1 lets us use geometric properties of clusters of a static network (ω n0 ) to study something which is dynamic: the load at time n at node v. That load may have come from a pathway that no longer even exists at time n. We further exploit this relationship, but to do this we must consider the concept of the dual web, defined in, for example, [7], and of whose definition we remind the reader. Consider the odd sublattice x + y odd and y ≤ 1} For any node v * = (x * , y * ) ∈ Z 2 odd we call the node (x * + 1, y * + 1) the right child of v * and we call the node (x * − 1, y * + 1) the left child of v * . Similarly, we call the node (x * + 1, y * − 1) the right parent of v * and we call the node (x * − 1, y * − 1) the left parent of v * . We define the set E 2 odd in the obvious way. The set of configurations {ω n : n ≥ 0} induces a set of configurations {ω * n : n ≥ 0} ⊂ {0, 1} E 2 odd by the following rule. If, in the configuration ω n , a node v = (x, y) is connected to its left child, we form a connection between the node v * = (x, y − 1) and its right child in the configuration ω * n by setting the image under ω * n of the edge in E 2 odd between v * and its right child to 1, and the image of the edge between v * and its left child to 0. If, on the other hand, v is connected to its right child in ω n then we set the image of the edge from v * to its left child under ω * n to 1 and the image of the edge from v * to its right child to 0. See Figure 4 and notice that we construct clusters in ω * n so that no occupied edges in ω * n cross any occupied edges in ω n . The upward paths in ω * n now resemble the downward paths in ω n . That is, the upward path starting at the node v * is a simple symmetric random walk which is killed at depth 1. Random walks starting at different nodes are independent until they meet, at which point they coalesce into one random walk. (This is similar to the coalescing random walks picture of the discrete web, described in [2], [11], [7].) There is an obvious physical interpretation for the paths in the dual web. For any two adjacent paths in the configuration ω n , there is a path in ω * n separating them. If the paths in ω n represent rills or drains, the paths in ω * n represent the divides or ridges between them. Just as divides between rills do not cross rills, paths in ω * n do not cross paths in ω n . We now characterize the load distributions for our model. For any node v = (x, y) (with depth k), let v * L = (x − 1, y) and let v * R = (x + 1, y). Consider the set of edges in the dual lattice contained in the paths emanating from the vertices v * R and v * L in ω * n until either (a) they meet at some vertex w * or (b) they reach a depth of 1. The set of nodes in Z 2 even in the interior of this set of edges is exactly the backward cluster of v in the configuration ω n .
We now make some definitions so that we can work with this load distribution. Let {X L i : i ≥ 2} and {X R i : i ≥ 2} be independent sets of random variables (also independent of each other) which take the values 1 and -1 each with probability 1 2 . For i ≥ 2 let Y i = 1 2 (X R i − X L i ) and for i ≥ 1, let W i = 1 + Y 2 + ... + Y i . Consider the stopping time τ = min{n : W n = 0} Up until the stopping time τ , the random variable W i represents the width of the backward cluster of the node v in the real lattice (only valleys and not separating ridges), where we only consider nodes in this cluster whose depths are between k − i + 1 and k. Therefore the total number of nodes in this partial cluster should be Now we can make an equivalent definition of the distribution of the load I v (n) by saying that for each fixed n, it is the same as the distribution of the random variable L k (n) := L min(τ,n,k) This variable is essentially a discrete integral of the symmetric random walk Let v be a node at depth k v and let w be a node at depth k w ≥ k v . For any n ≥ k v and for any l < k v we have P(I v (n) = l) = P(I w (n) = l) Therefore the limit exists in distribution. This limit is a.s. finite but has infinite mean.
Proof. On the event τ ≥ k v , i.e. the load originated from the top, I v (n) = L min(τ,n,kv) = L kv ≥ k v > l and I w (n) = L min(τ,n,kw) ≥ L min(τ,n,kv) > l Hence, we need only consider τ < k v .
The random variable L kv (n) is constant for n ≥ k v , so where in the last equality we use (6). Consequently, for any fixed l, the limit exists. By the definition (4), a random variable with this limiting distribution is Since τ ≤ L τ ≤ τ (τ + 1) 2 the third statement of the theorem will follow if we show that τ is a.s. finite and has infinite mean. But since the increments {Y i } of the random walk {W i } have mean zero, the walk is recurrent. In addition, it is a standard result that the entrance time of the set {0} has infinite mean. This completes the proof.

Dynamics
Now we investigate some aspects of the effect of η on the stability of configurations over time. As noted, the dynamics creates an interpolation between the discrete web and the discrete dynamical web. The evolution of the system mirrors some aspects of rill erosion, one being that nodes through which a large amount of water passes at time n 0 have a non-trivial probability to channel a large amount of water at any time n 1 > n 0 . The degree to which this is true depends on the parameter η, as we will see.

Load Correlation
where I v (n) = I v (n)−E(I v (n)). This quantity is only defined for n ≥ N because two load vectors for a box of depth N are in some sense incomparable if they are taken at times n 0 , n 1 with n 0 < N ≤ n 1 . For example, a node at depth n only has a maximum possible load of n(n+1) We can see in Figure 5 that the coefficient approaches 1 as η → 0. This makes sense because, as remarked in section 1.3, the η → 0 limit of the dynamics (in any fixed box) is the same as the dynamics (or rather non-dynamics) of the discrete web. Therefore the load vector for this box should be similar (if not the same) at any two times. As η increases, the correlation coefficient decreases and appears to approach 0. Indeed, additional simulations give the following data: for η = 10, 100, 1000, 10000, the coefficients were .2391, .0896, .0189, and .0101,. From the discussion of the η → ∞ limit given in section 1.3, the correlation coefficient should approach that computed from two load vectors from independent realizations of the discrete web.

de Finetti Measures
Whereas we can compare nodes at the top level to standard Pólya urns, we can compare lower level nodes to time-dependent input [14] or random input [13] Pólya urns. We start with an urn with R 0 red balls and B 0 black balls, as before, but we also have a time-dependent (or random) input sequence I = (I 0 , I 1 , ...). At time t = n we draw a ball from the urn and we return it to the urn along with I n balls of the same color. Notice that this process with I = (1, 1, ...) is just the standard Pólya urn.
To analyze these lower level nodes, we will also make use of a fundamental result in the theory of exchangeable variables.
Theorem 3.5 (de Finetti). Let (Ω, F, P) be a probability space and suppose that {X n } n≥0 are exchangeable {0, 1}-valued random variables defined on Ω. Then there exists a random variable F on Ω so that conditioned on F , the random variables X n are independent Bernoulli with parameter F .
It is easy to verify that if depth(v) = 1, then the {0, 1}-valued variables {D L v (n)} n≥0 are exchangeable. In our case, the variable F from Theorem 3.5 is actually Thus if we know the asymptotic fraction of left choices for a node, then our node is just flipping independent coins each second with the same bias. At lower levels, the variables {D L v (n)} are not exchangeable. However, they are asymptotically exchangeable. We use the definition of Kingman [12]. In the language of Theorem 3.5, let F be the random variable associated with the exchangeable variables {X n }. We call F the de Finetti measure for the sequence {Y n }.
Let v be a node with depth k ≥ 1.
Because of lateral translation invariance, the de Finetti measure for v depends only on the depth k. In light of this, we define θ k = de Finetti measure for row k With this framework we will be able to study the switching rate of each node v once we have the following lemma.
Lemma 3.10. For any node v, almost surely, Proof. Similar to the proof of [13, Theorem 2.3].
We now define the switching function s v for n ≥ 2 by . Define the switching rate S v (n) to be the time average of s v , that is Theorem 3.11. The n → ∞ limit of S v (n) exists a.s.
Proof. The proof is similar to the proof of [13, Theorem 2.3]. Let d n = D L v (n) and p n = P L v (n). A straightforward calculation gives Now, by Lemma 3.10. We must show that the last limit above is a.s. equal to where we use both equations in (11). Note also . We leave the reader to verify that M n is a martingale. Using L 2 -orthogonality of martingale differences, Therefore M n is an L 2 bounded martingale and converges a.s. This completes the proof.
If p v ∈ (0, 1) then neither of the choices L or R are transient for v. This prompts the question of whether or not the de Finetti measures θ k have atoms at 0 or 1. For any fixed k, the answer is no.
Theorem 3.12. For each k ≥ 1, the measure θ k has no atoms.
Proof. In [14,Theorem 4] it is shown that a time-dependent input Pólya urn's de Finetti measure cannot have atoms if there is a C so that I v (n) ≤ C for all n. For each realization of the dynamics and for each v, we have I v (n) ≤ kv(kv+1) 2 for all n. The result follows.
Corollary 3.13. Each node v has a nonzero asymptotic switching rate. Therefore, for each v, the states L and R are recurrent.
Proof. This is a direct consequence of Lemma 3.10 and Theorem 3.12. Here an interesting picture of our network emerges. On the one hand we may view the system as an infinite lattice (the lower half plane), where each node is a random input Pólya urn. The output of the urns at depth k at time n becomes the input of the urns at depth k+1 at time n+1. On the other hand, as remarked in section 1.3, we may first sample (non-independently) values {p v : v ∈ Z 2 even } from the de Finetti measures {θ k : k ≥ 1} to create an infinite array. As time n approaches infinity, the behavior of the system approaches the behavior of the same network in which each node v chooses to send its current load left with probability p v and right with probability 1 − p v , independently at each second. Therefore this picture is of a network of two variables, a realization of values p v from the de Finetti measures, and realization of dynamics which coincides with the dynamics of a much simpler network. This second network is an obvious generalization of the Dynamical Discrete Web. Figure 6 shows histograms for the de Finetti measures θ k for k = 2, 5, 9 and for values of η = .5, 1, 2. One sees that the measures become more biased as k increases (for fixed η). In other words, the mass of θ k is concentrated on domains closer to 0 and 1 than is the mass of θ k−1 . This would seem to imply that the η = 1 2 , row 2 η = 1 2 , row 5 η = 1 2 , row 9 η = 1, row 2 η = 1, row 5 η = 1, row 9 η = 2, row 2 η = 2, row 5 η = 2, row 9 expected asymptotic switching rate of a node at level k (which is 2p v (1 − p v )) must decrease with k. Similarly, if k is fixed and η decreases to 0, it seems that the expected switching rate should decrease. Figure 7 represents data given by simulations conducted with an erosion network with width 10 5 , depth 50, and η = .1, 1, or 10. The simulation ran for n = 1000 steps and at the end, switch rates for each node in the network were computed. In each row, each node's rate was averaged. Since two nodes v 1 and v 2 with the same depth k have independent behavior as long as they are at least a distance of 2k apart, the ergodic theorem gives that, as the network size approaches infinity, the resulting average should resemble the expected switch rate for a row. The above results results were averaged by row over 3 independent trials. Finally, the data were plotted by row. Not only do the average switch rates appear to decrease as k increases, there appears to be a non-trivial (i.e. non-zero and η dependent) limit for the switch rate. This indicates that for at least some values of η, the limit of the measures θ k (if it exists) is most likely not equal to 1 2 (δ 0 + δ 1 ).

Catastrophes
Next, we study the following situation. Suppose a node v at a large depth k (for this section we assume the depth is at least 2) starts with a small input load and keeps a relatively small input load until a much later time. Then v's load changes dramatically. If this new load is sufficiently large, it could bring v's de Finetti measure much closer to 1 2 (δ 0 + δ 1 ). This analysis is from the point of view of the node v, whereas the analysis of the last half of the section will be from the point of view of the parent. Let Definition 3.15. For any n ≥ 1, define the flood ratio F v (n) = Iv(n) Av(n) . For c ≥ 1, we say that a flood of order c occurs at time n if F v (n) ≥ c.
Proposition 3.17. For any v, A v (n) has a limit a.s. Therefore, Proof. We show the first statement by induction on depth(k). Clearly this is true if depth(v) = 1. Otherwise, let w 1 be the left parent of v and assume that for all nodes w with depth equal to that of w 1 , A w (n) has a limit. Let N R w1 (n) be the number of i ≤ n such that D R w1 (i) = 0. By Lemmas 3.7 and 3.10, exists. The same argument shows that, if w 2 is the right parent of v, lim n→∞ exists. Therefore, For the second statement of the proposition, we use Lemma 3.14 to see that Remark 3.18. The above proof shows also that lim n→∞ Remark 3.19. A simple extension of the above proof using Theorem 3.12 shows that for any level k, the distribution of the time-average of the load for level k has no atoms.
From the previous proposition we know that if v has an asymptotic average input load A v , then v will have infinitely many floods of order k(k+1) 2Av+ for any > 0. This number can be quite large if A v is small. We study numerically the rates of these large floods as they relate to both depth(v) and η.
For a fixed node v, define the random measure Let {v i } i≥1 be an enumeration of the nodes with the same depth as that of v. By the ergodic theorem, the average as M → ∞. See Figure 8. The graphs come from a simulation run for 10 5 seconds on a network with a width of 10 4 nodes and a depth of 50 nodes. We graphed the density function for the measure 1 η = 1 10 , row 5 η = 1 10 , row 20 η = 1 10 , row 50 η = 1, row 5 η = 1, row 20 η = 1, row 50 η = 10, row 5 η = 10, row 20 η = 10, row 50 Figure 8: Histograms approximating the measure E(µ v,n ) for n large.
As η → 0 with a fixed row or as the row increases with fixed η, each measure seems to concentrate its mass at 0 and 1. In other words, the measure of any interval which does not include either of these two points appears to approach 0. In addition, Table 1 shows that as η → 0 with a fixed row or as the row increases with fixed η, the expected fraction of time during which a large flood occurs (ratio above 5) increases. Since the time average of flood ratios approaches 1 as n → ∞, the above facts indicate a trend that as η → 0 or as the row increases, the time variance of measures increases, giving more possible variability of the flood ratios. Although this variance seems to increase, values near 1 (on the x-axis) show that the fraction of time flood ratios spend near 1 increases. In spirit, this is in accordance with previous results, as we explain. Figure 7 shows that as the row increases, the expected switching rate of a node decreases. It is reasonable to believe that the same conclusion holds if the depth is fixed but η decreases. Therefore the network prefers to be more static in these circumstances and we would expect a node to receive a relatively constant load, forcing flood ratios to be near 1.
We now change focus to the parent. Make the following definition. .03711 Table 1: Expected fraction of time during which a large flood (ratio at least 5) occurs.
We investigate the relationship between floods and catastrophes. Figure 9 shows the expected fraction of right catastrophes from the left parent which result in a flood of at least the same magnitude. The simulation was performed with a network of width 1000, depth 50, η equal to either .1, 1, or 10, and for a duration of 10 5 seconds. The calculation of fractions was only made between 9000s and 10000s and we only consider catastrophes with ratio at least 1. It is clear from the figure that as the row increases, this expectation decreases. As η decreases for a fixed row, the expectation also decreases. As n → ∞, by Remark 3.18 and Lemma 3.10. Therefore, we may use Remark 3.19 to show that almost surely for large n a right catastrophe of order c occurs for the node v at time n if v has a flood of order c at time n − 1. In other words, whenever a node receives a flood of order c at a large time, it has either a left or right catastrophe of the same order at the next second. Now we may interpret the probability that a node has a flood given that its left parent has a right catastrophe as the probability that a parent's catastrophe incites a catastrophe in the child. This would be a step of a possible catastrophe cascade. The simulation results indicate that cascades become less present at lower levels (on average) but that they should never cease to exist. Two questions naturally arise. Given that a node has a right catastrophe of order c, how far does its catastrophe cascade travel? At each step of the cascade, the relevant (right or left) catastrophe ratio will generally increase. Indeed, a child node may even receive a catastrophe from both parents. How large does this ratio become in a typical cascade?

Conclusion
We have shown that the erosion model exhibits many properties of rill erosion. Each node chooses a random initial direction (right or left) in which to send sediment and further such choices become biased at a rate largely determined by the parameter eta. This is similar to the method by which rills are cut into a hillslope. As more water and sediment flows through a rill, a channel is cut deeper, giving reinforcement to the path, making it more likely to carry sediment in the future. Though the dynamics manifests itself through reinforcement, no fixed node can become fully biased (i.e. have a de Finetti measure equal to a sum of two delta masses). That is, since each node has a non-trivial asymptotic switching rate, sediment flow emerging from it will take both a left and right path a positive fraction of time. This rate of switching appears to decrease as we move further down the hill.
There are a number of questions which deserve careful analysis. Do the measures θ k have a limit? If so, one would expect the limit to be the de Finetti measure associated with the "infinity process". To define this process, we start with a lattice of nodes which extends infinitely far in both positive and negative y directions. Since the behavior of a node v at time n in the present model depends only on the nodes in the n − 1 levels above it we may consider the input to the node v at time n in the infinity model to be a function of the output of this finite number of ancestors. In the same way we have analyzed in this paper, it is possible to show that a de Finetti measure θ ∞ for this process exists and that where the term inside the limit is the measure given by θ n (n)([a, b]) = P(P L v (n) ∈ [a, b]) for depth(v) = n Does the measure θ ∞ have atoms for some values of η? If so, is there a critical η * so that for 0 < η < η * , θ ∞ has atoms? If the limit of θ k (assuming it exists) is not the same as θ ∞ , does this limit have atoms and is there a critical η associated with it?