Transport capacity and saturation mechanism in a real-space cellular automaton dune model

. In a real-space cellular automaton dune model, individual physical processes such as erosion, deposition and transport are implemented by nearest neighbor interactions and a time-dependent stochastic process. Hence, the transport capacity, the saturation mechanism and the characteristic wavelength for the formation of dunes are emergent properties that can only be determined a posteriori from the output of the numerical simulations. Here we propose a simpliﬁed version of the model to establish asymptotic relations be-tween the microscopic erosion-deposition-transport rate parameters and the characteristic length and time scales of the ﬂux saturation mechanism. In particular, we show that, in the cellular automaton, the saturation length is a mean transport distance controlled by the deposition of mobile sedimentary cells. Then, we discuss how these results can be used to determine the sediment ﬂux within dune ﬁelds and the rate parameters of a new class of discrete models that concentrate on the effect of heterogeneities in grain-size on dune mor-phodynamics.


Introduction
Given the apparent complexity of dune fields, it is often necessary to rely on numerical methods to estimate the overall sediment flux and the migration history of individual grains (Zhang et al., 2013).With this purpose in mind, there are actually two efficient numerical approaches for the modeling of dunes and transport phenomena at the interface between a granular bed and a flow: -Continuous models are based on a set of partial differential equations that combine the principle of mass conservation with an analytical description of the flow and relations between the transport capacity and topography (Sauermann et al., 2001;Kroy et al., 2002;Andreotti et al., 2002;Fischer et al., 2008).Their main advantage is that all these relations have well-defined scaling properties.Nevertheless, using such a continuoushomogeneous framework, it remains difficult to evaluate how the details of the formalism (e.g.parametrization of the relation between the curvature and the basal shear stress, the properties of recirculation flows) affect the overall dynamics of the system, in particular in presence of superimposed bed forms and dune-dune interactions.In addition, large domains are numerically difficult to achieve with continuous models.
-Cellular automaton models are discrete models that implement a set of rules to mimic sediment transport and the effect of the flow on topography (Nishimori and Ouchi, 1993;Werner and Gillespie, 1993;Werner, 1995).Their main advantage is that they reproduce a huge variety of bed forms at a reasonable computational cost.Nevertheless, they are constructed from an arbitrary elementary length-scale which is always difficult to define.
The real-space cellular automaton dune model developed by Narteau et al. (2009) is an attempt to link these two lines of research.It belongs to the class of discrete models and can be described as a hybrid approach that combines a cellular automaton of sediment transport with a lattice-gas cellular automaton for high Reynolds-flow simulation (Frisch et al., 1986;d'Humières et al., 1986).At the elementary length scale l 0 of the cellular automaton of sediment transport, a 3-D square lattice is a real-space representation of the physical environment under investigation (Rozier and Narteau, 2013).Each cell can be in a finite number of states and all the physical mechanisms are implemented by nearest neighbor interactions and a time-dependent stochastic process.In practice, there is a finite number of transitions associated with erosion, transport and deposition.Each transition is characterized by a Poisson rate parameter Λ expressed in units of frequency using the characteristic time scale t 0 of the model.Thus, rate parameters can be varied at will to modify the relative amplitude of the physical mechanisms under consideration (Narteau et al., 2001).Meanwhile, a lattice-gas cellular automaton is implemented in the fluid state of the model of sediment transport to simulate a basal shear stress, which is the main control for the erosion rate.Then, because the rebound dynamics of the lattice-gas model depends on both the topography and the top boundary condition, we end-up with a fully coupled system for the fluid and granular phases.As these couplings result only from the configuration of cells and the subsequent patterns of interaction, they are responsible for a wide range of emergent dune features such as superimposed bed forms (Zhang et al., 2010) or arm growth (Zhang et al., 2012).
The characteristic length and time scales {l 0 , t 0 } of the real-space cellular automaton dune model have been entirely defined with respect to a linear stability analysis (Narteau et al., 2009).Used in many other branches of physics, this standard technique identifies stable/unstable regimes across the entire parameter space of the model to reveal the mechanisms that, starting from a flat sand bed, spontaneously generate bed forms.Then, in all natural environments where such a dune instability is observed, the model can be used by setting the l 0 and t 0 -values to the appropriate values of the relevant physical parameters (Claudin and Andreotti, 2006).In this way, we can quantitatively compare the numerical results to observations in nature and laboratory experiments, as well as to the predictions of other numerical models.In particular, we can set up the model to reproduce aeolian or subaqueous dune features on Earth and other planetary bodies.For example, the barchan dune presented in Fig. 1 has a height of 25 l 0 and a propagation speed of 0.028 l 0 /t 0 in units of the model.In arid desert on Earth, it corresponds to a dune height of 12.5 m and a propagation speed of 17.5 m yr −1 for a grain size d =200 µm; in water, it corresponds to a dune height of 1.3 cm and a propagation speed of 10 m h −1 ; on Mars, it corresponds to a dune height of 375 m and a propagation speed of 9 cm yr −1 for a grain size d =87 µm (see Zhang et al., 2013, Table 2).
Here our goal is to quantify the saturated flux and the saturation length of our model with respect to the transition rates that characterize the mechanisms of erosion, transport and deposition at the elementary length scale of the cellular automaton of sediment transport.Such analytical relations are critical to predict sand flux on natural dune features.As described in the discussion, they can also help to constrain the rate parameter values of a new class of models that accounts for dune pattern formation in heterogeneous sand bed (i.e. grain size, shape and density).

A simplified real space cellular automaton model
To measure sediment flux in the real-space cellular automaton model, we reproduce the experimental protocol proposed by Bagnold (1941) using cells in a neutral state (i.e.no interaction) to construct the solid boundaries of our virtual physical environment.As shown in Fig. 2a, the initial condition is a thick flat layer of sediment preceded by a non-erodible flat bed in the direction of the flow.When the flow arrived on the erodible part of the bed, the sediment load increases.Then, the principle of the experiment consists in measuring sediment flux as a function of distance and time in the direction of flow.As shown by Andreotti et al. (2010), such a saturation experiment can be performed in wind tunnels although it is impossible to maintain a flat sediment bed over long times.
The two main reasons are (a) the formation/migration of an erosion front at the discontinuity from non-erodible to erodible bed and (b) the formation of sedimentary structures such as ripples or dunes.The same problem exists in our simulations.It is even more pronounced due to the discrete nature of the model and the relatively high ratio between the elementary scale of the model and the characteristic grain size of the sediment layer (i.e.l 0 /d > 10 3 ).Fortunately, our numerical approach offers the opportunity to overcome this obstacle.
To eliminate the influence of bed topography on the flow, we simplify the cellular automaton of sediment transport using the same three states (fluid, immobile sediment, sediment in motion) but only three transitions (Fig. 3).This new set of transitions has the advantage to preserve the initial bed topography composed of stable cells.Nevertheless, each transition of erosion generates a new mobile sedimentary cells and therefore violate the principle of mass conservation if this cell have a non-null density.That is why we consider here that the mass of mobile sedimentary cells is null.Conceptually, we should admit that, when a sediment cell becomes mobile, it is immediately replaced by a stable sediment cell with exactly the same properties.By symmetry, when a mobile sediment cell becomes stable, it is immediately removed from the surface of the bed.Then, the flow is only disturbed by mobile sedimentary cells and is not affected by changes in bed topography (Fig. 2b, c).In this case, we can focus on the saturation mechanism by studying the evolution of shear stress as well as the screening/armoring effect produced by the mobile sedimentary cells.Overall, the simplicity of this new cellular automaton of sediment transport also allows to derive (and numerically verify) a number of analytical relations involving the rate parameters of the model (i.e.Λ e , Λ t , Λ d ).

Saturation length and saturated flux
In the model, the basal shear stress result from the dynamical coupling between the topography and the flow.To account for the initial topography (Fig. 2), we first stabilize the flow at the beginning of each simulation by running the  lattice-gas cellular automaton model alone.Thus, we reach a quasi-stationary equilibrium state and the basal shear stress can be described by the same probability density function over space and time.Then, we activate the transitions of doublets (Fig. 3) incrementing the time for each of them (see Narteau et al., 2009, Appendix B).This corresponds to a sudden change in flow speed from 0 to the mean value at the end of the flow stabilization period.
In vertical planes perpendicular to the flow, we systematically count the total number of transport transitions to average the sediment flux over the entire width of the domain.Figure 4a shows the sediment flux with respect to distance downstream for Λ e t 0 = 4 × 10 −5 , Λ t t 0 = 50 and Λ d t 0 = 2. Above the erodible area (i.e. x > x 0 = 50 l 0 ), we observe that the flow relaxes to an equilibrium value Q sat over a characteristic distance l sat .The quality of the approximation by a law of the form with .45 l 2 0 /t 0 , x c = x 0 and l sat = 27 l 0 shows that this relaxation process is exponential.We discuss below the Q sat and l sat -values according to the rate parameters of the model of sediment transport.
Before, to better understand the origin of the saturation mechanism in the model, we take the same initial condition as in the previous numerical experiment (Fig. 4a), but add a non-erodible zone downstream.Fig. 4b shows how the sediment load returns to zero when the flow reaches the nonerodible zone at x 1 = 500 l 0 .Once again, Eq. ( 1) with Q 0 = Q sat = 3.45 l 2 0 /t 0 , Q 1 = 0, x c = x 1 and l sat = 27 l 0 shows an excellent fit to the numerical data.Then, the sediment load relaxes towards its new equilibrium value over the same characteristic distance l sat and the charge and discharge mechanisms follow the same exponential process.This indicates that, in our simplified cellular automaton model, the saturation mechanism does not depend on specific sand bed properties (e.g.porosity, cohesion) or particle-particle interactions (e.g. grain ejection by saltating grains).Instead, it is entirely controlled by deposition as in suspended load (Ashida and Okabe, 1982).Then, the saturation length can be simply described as the mean distance over which a mobile sedimentary cell is transported.
Let us consider that, in our model, the number X of transport transition of a mobile sedimentary cell is a discrete random variable with probability density function f .Then, we can write (2) In the stochastic dynamical system implemented at the elementary length scale of the model of sediment transport, a mobile sedimentary cell moves if it participates in a transition of transport (Fig. 2).It becomes stable if it participates in a transition of deposition.For a given mobile sedimentary cell, these two transitions occur with probabilities Then, f can be described as the probability for a cell to operate X transitions of transport before a transition of deposition: By injecting these last three relations in the expression of l sat , we obtain (5) Thus, we show that the saturation length is directly the ratio between the characteristic time for deposition and the characteristic time for transport.For a given flow strength, this dependence of the l sat /l 0 -value on deposition and transport rates is exactly the same as in many continuous models that estimate transport from a number of moving particles (Charru and Hinch, 2006).Not surprisingly, Fig. 5 shows that the prediction of the numerical model are in perfect agreement with the analytical solution for a wide ranges of Λ t and Λ d -values.Within the 3-D cubic lattice of our cellular automaton, the normalized saturated sediment flux Q sat /l 2 0 is the number n t of transition of transport per unit of time ∆t in a vertical plane perpendicular to the flow.Then, the probability for a mobile sedimentary cell to be observed somewhere depends only on a transition of erosion upstream followed by a sufficient number of transitions of transport.In all respects, the saturated flow can be written where P e = Λ e ∆t is the probability to have a transition of erosion in the limit of ∆t → 0. This formula shows that the saturated flow is the direct product of the erosion rate by the average distance traveled by a mobile sedimentary cell.Figure 6 shows that this analytical solution is verified by the numerical experiments.

The saturation time and the density of mobile sedimentary cells
The sediment flux does not adjust immediately to its local Q sat -value.It always takes some time to reach this dynamic equilibrium state.To estimate this time, we can count the number N m of mobile cells on the erodible bed.These mobile cells protect the bed from further erosion.More exactly, because they cannot move vertically, they prevent erosion of cells that compose the flat erodible bed just below them.Then, if N s is the number of erodible cells exposed to the flow, we can write (7) Such a system of coupled linear equations have a solution of the form exp(λt).The eigenvalues of the coefficient matrix are the solutions of The only non-zero solution is equal to −(Λ e + Λ d ).As a consequence the number of mobile sedimentary cells in the system relaxes exponentially to its equilibrium value over a characteristic time Over short time just after the stabilization of the flow and the triggering of the transitions of the model of sediment transport, we verify numerically this exponential relaxation process at all location downstream of the transition from a non-erodible to an erodible bed (Fig. 7).Interestingly, using Eqs. 5 and 9, we can also show that l sat /t sat ≈ l 0 Λ t when Λ d Λ e .Then, the ratio between the saturation length and the saturation time is the mean velocity of the mobile sedimentary cells.
Finally, the equilibrium solution of our system of coupled linear equations is simply Since the density of mobile cell on the bed can be written we have With Fig. 8, we verify numerically this relation within the simplified version of the cellular automaton dune model.Such a dependence of the density of mobile sedimentary cells on Λ e and Λ d -values explains why we systematically take Λ e Λ d in the numerical experiments.Indeed, when d s 1, there are almost no interaction between mobile sedimentary cells and their proportion on the bed do not alter the magnitude of the overall sediment transport.

Concluding remarks
The saturation length l sat is of primary importance in determining the characteristic length scale for the formation of dunes (Hersen et al., 2002;Elbelrhiti et al., 2005).Then, the most relevant time scale of the underlying instability may be derived from the saturation flux Q sat .Here we address these scaling issues by showing how we can estimate the l sat and Q sat -values in a real-space cellular automaton dune model (Narteau et al., 2009;Zhang et al., 2010Zhang et al., , 2012)).More exactly, we use a simplified version of this discrete approach to derive analytical solutions that only depend on the rate parameters of the model of sediment transport .Our results confirm that the saturation mechanism and the transport properties may be directly related to the magnitude of the {Λ e , Λ t , Λ d }-values.Obviously, as the bed topography and the variability of the flow increase, the proposed relations may be not correct to the nearest degree under natural conditions as in the full version of the cellular automaton dune model.However, the dependence on the rate parameter values can be generalized to all models that incorporate erosion and deposition rates to predict transport from the motion of independent particles.For example, we find that the saturation length is a deposition length as for suspended load (Ashida and Okabe, 1982) and in other models dedicated to sediment transport in water flows (Charru, 2006;Lajeunesse et al., 2010).Although the saturation length seems rather controlled by grain inertia in aeolian settings (Andreotti et al., 2010), the same scaling can be performed to derive, from the l sat -value, the most unstable wavelength of the instability mechanism responsible for the formation of dunes.Then, the inverse of the growth rate 1/σ max of this most unstable mode has to be much smaller than the saturation time to ensure a quasi-stationary sediment flux in the neighborhood of the crests.This is also the case in our model when Λ d Λ e .
The relation between Q sat and the {Λ e , Λ t , Λ d }-values offers the opportunity to estimate variations of transport properties across the synthetic dune fields produced by the realspace cellular automaton dune model.Locally, it is necessary to know the probability distribution function of the basal shear stress g(τ s ) and the dependence of the erosion rate on this basal shear stress Λ e (τ s ).If all the other parameters remain constant, the local saturated flux is In addition, as it take some time and distance to adapt to this flux, the model will offer the opportunity to estimate the impact of the saturation mechanism in more realistic bed form configurations.Nevertheless, given the changes in topography, the formation of secondary flow and the permanent evolution of the dune features, it is likely that the local fluxes of sediment will have a significant dispersion around the estimates provided by Eq. ( 13).
In order to study polydisperse granular beds with the realspace cellular automaton dune model, different transition rates can be assigned to different sedimentary states.Then, we can use the results presented here to constrain the transition rate values with respect to the physical properties of the granular medium (density, grain size).With this purpose in mind, we can take advantage on empirical or theoretical relations derived from observations (Lettau and Lettau, 1978;Sørensen, 1991;Iversen and Rasmussen, 1999;Charru, 2006;Lajeunesse et al., 2010;Houssais and Lajeunesse, 2012).For example, considering that the deposition rate is proportional to the settling velocity, we can take where d is the grain diameter, ∆ρ = (ρ g −ρ f )/ρ f and ρ g and ρ f are the grain and the fluid density, respectively.Similarly, for a given and uniform basal shear stress τ s , we can define from the inertial balance between the weight of a particle and the drag force.However, in this case, the shear stress also incorporates a dependence on grain size as well as a threshold value τ c below which Λ e = 0 (see Narteau et al., 2009, Eq. 2).
In our cellular automaton, when the density of mobile sedimentary cells increases (Fig. 8), it may happen more frequently that two of them are aligned in the direction of the flow.Such a configuration prevents the occurrence of a transport transition for the upstream mobile sedimentary cell.So, there is a feedback of the number of mobile sedimentary cells on transport properties.If this feedback takes different forms (see the paragraph below), it is important to note that our analytical solutions of Q sat and l sat are only valid when the occurrence of transport transitions is not affected by the number of mobile sedimentary cells.To solve this problem, it is possible to add a dependence between vertical motions of mobile sedimentary cells and their number.Nevertheless, this is a new level of complexity which is unnecessary for dune morphodynamics since d s < 0.1 across the entire domain.
Even in this simplified model, there is a retroaction of the sediment in motion on the flow and transport properties.First, mobile sedimentary cells create a dynamical armoring of the bed that prevents for further of erosion the underlying sedimentary layer.Second, they also reduce the basal shear stress by increasing the roughness of the sedimentary layer.This second effect has not been studied here because it requires more realistic trajectories for mobile sedimentary cells.Interestingly, the real-space cellular automaton model possesses all the ingredients to be applied at the microscopic scale of a grain in order to concentrate on dynamic equilibrium of the transport layer in the near-bed region.

Fig. 1 .
Fig. 1.Morphodynamics of barchan dunes in a wind tunnel using the real-space cellular automaton dune model.(a) Morphology of a barchan dune in a quasi-stationary equilibrium state at t/t0 = 3962 and t/t0 = 9890.The flow is from the bottom to the top (black arrow) and the initial condition at t = 0 is a conical sand pile upstream.An equilibrium state can be reached because all the sediment ejected downstream is randomly reinjected upstream.Note the presence of superimposed bed forms.(b-c) Sedimentary fluxes and basal shear stress on the barchan dune (see Narteau et al., 2009, for details about the model).

Fig. 2 .
Fig.2.The real-space virtual environment: (a) Initial condition to reproduce the experiments ofBagnold (1941).Fluids cells (cyan) and neutral cells (yellow) that make up the ceiling are transparent for the sake of readability.Immobile sedimentary cells are shown in black.(b) Distribution of mobile sedimentary cells (red) when the system reaches a dynamic equilibrium state.Note the increasing number of mobile sedimentary cells in the direction of the flow from the limit between the non-erodible and the erodible bed.(c) A vertical slice of cells in the direction of the flow.Note the presence of the ceiling that confines the flow.

Fig. 3 .
Fig. 3. States and active transitions of doublet in the simplified cellular automaton model for sediment transport.Three transitions rates Λe, Λt and Λ d with units of frequency (1/t0) are associated with erosion, transport and deposition, respectively.This simplified set of transition ensure that the sediment bed remains flat.When it is not the case, the cellular automaton model has to include horizontal transitions for erosion and diffusion as well as vertical transitions for deposition and transport (see the 10 transitions of the complete model in Narteau et al., 2009, Fig. 2).

Fig. 4 .
Fig. 4. Saturation length and saturated flux in the simplified cellular automaton model.(a) Downstream of the transition from a non-erodible to an erodible bed, the sediment flux adjusts to the basal shear stress and reaches a stable value Qsat over a characteristic distance lsat.(b) An erodible bed between two non-erodible zones: charge and discharge mechanisms follow the same exponential relaxation regime with the same characteristic length lsat (see blue/black lines and Eq. 1).

Fig. 5 .
Fig. 5. Relation between lsat and the rate parameters of cellular automaton model.(a) Increase of the saturation length with respect to an increasing Λt-value for Λ d t0 = 2. (b) Decrease of the saturation length with respect to an increasing Λ d -value for Λtt0 = 0.4.(c) lsat with respect to Λt/Λ d , the mean distance traveled by a mobile sedimentary cell: Λ d t0 = 2 (blue squares), Λtt0 = 0.4 (red crosses).

Fig. 6 .
Fig. 6.Relation between Qsat and the rate parameters of the cellular automaton model.(a) Increase of Qsat with respect to an increasing Λt value for Λ d /t0 = 2. (b) Increase of Qsat with respect to a decreasing Λ d -value for Λt/t0 = 1.(c) Qsat with respect to Λt/Λ d , the mean distance traveled by a sedimentary cell: Λ d t0 = 2 (blue squares), Λtt0 = 0.4 (red crosses).The slope of the linear relation is proportional to Λe.

Fig. 7 .
Fig. 7. Relation between the saturation time tsat and the rate parameters of the cellular automaton model.(a, b) Evolution with respect to time of the density of mobile sedimentary cells in vertical plane perpendicular to the flow at 2 l0 and 20 l0 from the transition from a non-erodible to an erodible bed.Despite different stationary values, note a similar relaxation time at these two positions.(c) Density of mobile sedimentary cells in vertical plane perpendicular to the flow with respect to the distance from the transition from a non-erodible to an erodible bed.