2.1 Description of the CA model for snow transport

Cellular automata modelling of snow transport consists of first discretizing the three-dimensional space to be simulated into cells that can only have one of three possible states: air, deposited snow (DS) or mobile snow. At the beginning of the simulation for example, one would typically consider a snow-covered surface with no mobile snow cells. Each of these cells in this study has dimensions of l₀ in the horizontal direction and h₀ in the vertical direction. The ratio of the horizontal and vertical length scales is a model parameter and in our case, $l_{\mathrm{0}}/h_{\mathrm{0}}=\mathrm{5}$. This allows for greater vertical resolution in the model and is required for simulating snow bedforms that are typically much flatter and shallower than sand bedforms (Filhol and Sturm, 2015). Note that this ratio is a free parameter in the model and was chosen as such to have a balance between increased computational expense (due to increased number of cells in the vertical) versus better resolution of the snow surface. In any case, test simulations performed with even higher vertical resolution (not shown) did not modify the results much.

The model proceeds in the form of transitions of the cell states based on certain rules. In the ReSCAL model, these rules are based on nearest-neighbour cell pairings known as doublets. A list of phenomenologically important doublet combinations is made a priori. Each of these doublet combinations can transition to a different doublet combination. These doublet groupings are presented in Fig. 1a, where they are grouped into six different phenomenological mechanisms, namely deposition, erosion, transport, gravity, diffusion and avalanching. Before we elaborate further, it would be beneficial to avoid connecting the above-listed mechanisms directly to physical processes with the same names to avoid confusion. For example, erosion in the CA model is not directly linked to the different grain-scale entrainment mechanisms discussed in Sect. 1 even though at bulk scale, it is intended to produce the same effect.

Each of the transitions (from doublet a to doublet b) described in Fig. 1a is associated with a transition-rate parameter, Λ, with subscripts used to distinguish transition rates of different mechanisms. These transition rates have a dimension of inverse time or frequency. This allows for a timescale to be introduced into the CA model (the length scale being introduced by the discretization of the space). The simulation progresses in a time-dependent stochastic fashion with each transition being regarded as a Poisson process. At each time step, all the doublets eligible for transition are counted and a total transition rate of the system is computed, which is simply a total of all the transition rates weighted by the number of doublets eligible for each transition (note that computing a global transition rate is possible for Poisson processes using the additive property). The global transition of the system is also a Poisson process, and thus the time step can be stochastically chosen. The transition to be performed is chosen next based on the conditional probability of each transition. Finally, one of the doublets eligible for the chosen transition type is randomly picked and a transition is performed. In summary, at each instance of the simulation, the transition type to be performed, the doublet to be transitioned as well as the value of the time step itself are stochastically chosen based on the value of the transition rates and the global state of the cells. A list of transition rates for the CA model used throughout this study is given in Table 1.

It is clear, therefore, that only the relative values of the transition rates are important. These rates are chosen to roughly reflect observations in reality. For example, at bulk scale, transport is far more rapid than erosion or deposition. Erosion itself is typically faster but a more local process compared to deposition, which is slower but occurs over larger length scales (consider a scenario of deposition of effluents from a plume). The processes of gravity and cross-stream diffusion (due to turbulence) are the fastest and slowest respectively. Avalanching is implemented in an extremely simple manner following Bak et al. (1988) – if the local slope is larger than a critical angle, a cell at that location is moved to a random location down the slope. We realized a posteriori that for the scale of the dunes simulated in this study, avalanching is of negligible importance, and thus it is not discussed further.

2.2 Description of the LGCA model for flow over snow surfaces

The flow overlying the evolving surface is simulated using the lattice-gas cellular automata (LGCA) approach. In this numerical technique, the fluid is modelled as a set of particles lying on the nodes of a square (or cubical) mesh that is called a lattice. A particle must lie on one of the nodes of the lattice at all times and each node of the lattice can hold only one particle at any given time. Furthermore, a particle can move only to the nearest or next-nearest neighbours in one time step of the simulation. In other words, the velocities that a particle may have are extremely limited, in both magnitude and direction. This is illustrated in Fig. 1b. The top panel in this figure shows for example, that a particle may move only to nearest or next-nearest lattice points. Thus, a particle may move in one of eight different directions and have two possible speeds. A time step in the LGCA consists of two sub-steps. The first is known as the propagation step, shown in the lower left panel, where all particles with non-zero velocities move to their destination lattice nodes. The incoming particles at each node are represented by the green arrows that all lead into lattice points. Note that different lattice points have different numbers of incoming particles. After this sub-step, it may happen that multiple particles (temporarily) lie on the same node. To impose the constraint that a node may have only one particle, N-body collision calculations are performed between the incoming particles, and the particles obtain new velocities. These new velocities are represented by outgoing green arrows from many lattice points in the lower-right panel of Fig. 1b. This step is known as the collision step. The lattice velocities as well as the N-body collision rules are adopted from D’Humières et al. (1986). As boundary conditions for the fluid particles, the collision of a particle with a solid object or a wall is modelled as a simple elastic collision (similar to the model of the ideal gas). The LGCA methodology further provides a simple way to convert lattice-based velocities to “real” velocity of the flow. Typically it is simply the average of the velocities of particles in a given neighbourhood.

The LGCA approach is in some sense a reduced-order model of the full Navier–Stokes equations achieved by imposing strict constraints on directions and velocity values. Its development began with the pioneering work of Frisch et al. (1986) and was the precursor to more advanced lattice–Boltzmann methods. The LGCA model in ReSCAL is adopted from D’Humières et al. (1986). This modelling technique has the advantage that it is an extremely rapid method to simulate flow over complex surfaces, which is typically quite challenging for more traditional fluid simulation techniques such as large-eddy simulations or even Reynolds'-averaged Navier–Stokes (RANS) models. In the context of its use in this study, the LGCA, by simulating flow over complex bedforms on the surface, provides values of the surface-shear stress at every location of the surface. The surface-shear stress is essentially computed as a gradient of the velocity (computed using LGCA) in the direction normal to the local surface. The surface (as well as the normal) is of course the result of the CA model for snow transport.

Of all the transition types, erosion is the only one directly linked to the morphology of the surface while all other transitions are independent of their location in 3-D space. This link between erosion and morphology is established by modifying the transition rate for erosion according to location and making it a function of surface-shear stress, which due to surface morphology is heterogeneous. Thus, areas in the domain with larger shear stress experience more erosion. In the ReSCAL model, the erosion rate is linearly dependent on excess shear stress (τ_s−τ₁) as

where τ_s is the local surface-shear stress and τ₁ and τ₂ are the lower and upper limits of the linear regime. For τ_s<τ₁, the erosion rate is set to zero while for τ_s>τ₂, the erosion rate is at the maximum possible value of Λ₀. The values of Λ₀ and τ₂−τ₁ are constant and act only as scale and slope parameters respectively. τ₁ is equivalent to a threshold shear stress and is kept as a free parameter. As is explained in Sect. 2.4, it is used to specify the wind speed. A subtle point to note is that the shear stress quantities are scaled with respect to the shear stress scale, τ₀. This scale is not discussed or identified in this study as it is not directly relevant to finding the kinematic scales of length and time.

Table 1CA model scales and parameters using in this study. See Fig. 1a for more information about the transitions and the doublets involved.

^a Ratio of vertical to horizontal rates for both erosion and transport mechanisms. ^b Enhancement factor for deposition due to DS-type cells.

Download Print Version | Download XLSX

It is through the erosion mechanism that the flow simulation performed using LGCA and the snow transport performed using CA are tightly coupled (see Fig. 1c). This robust coupling between the flow, its modification due to the undulating surface and the evolution of the surface itself is the defining feature of the ReSCAL model and thus makes it perhaps one of the best cellular automata models for simulating surface morphology that exists. This can be evidenced from the successful application of the ReSCAL model to simulate various complex sand bedforms on Earth (Zhang et al., 2010; Ping et al., 2014; Lü et al., 2017) and Mars (Zhang et al., 2012) and also gives us the confidence of introducing this model to the cryospheric science community.

2.3 Implementation of a simple sintering mechanism for deposited snow

In this study, we intend to focus on the effect of sintering on snow bedform dynamics. For this purpose we implement a simple sintering mechanism that suits the CA modelling approach. Every time a cell transitions to the “deposited snow” (DS) type, the time step of the simulation is noted. This provides a way of measuring the age of the snow cell, i.e, the period of time a snow parcel rests in one place. All transition rates for transitions consisting of DS-type cells are then multiplied by an erodibility factor, f_E, which is simply chosen to be a linearly decreasing function of the age of the cell with the most important transition being erosion. Thus, Eq. (1) describing the erosion mechanism is modified as

with the erodibility factor f_E being defined as

where t is the current time of the simulation, t_dep is the time when the cell transitioned to a DS-type cell and t_s is the sintering timescale, i.e, the time after which a DS-type cell cannot perform any further transitions. This is a new parameter that must be chosen for such a model and we have chosen it to be 24 h. Thus, after 24 h, a DS-type cell will permanently remain in the same location for the rest of the simulation.

Cells that remain erodible, i.e., have ages less than 24 h are henceforth denoted as eDS-type cells while immobilized; sintered cells with ages greater than 24 h are denoted as neDS-type cells. This model is directly inspired by the approach of Filhol and Sturm (2015) and circumvents the requirement of accurate modelling of the highly complex and poorly understood phenomenon of sintering. The erodibility factor, f_E, of snow cells as a function of their age, t−t_dep, is shown in Fig. 1d. The coloured boxes in Fig. 1d are shown only to visually represent transition of eDS-type cells, coloured light blue, to neDS-type cells, which are coloured dark blue. A similar colour scheme is adopted in the following sections. Note that in the context of this study, erosion and erodibility refer only to the action of the wind. Erosion due to the impact of airborne snow grains on a moderately sintered snow bedform is not simulated and is kept for future work.

2.4 Finding the length and timescales of the CA model

A crucial contribution of Narteau et al. (2009), in addition to the development of the ReSCAL model, was the development of a methodology of translating results of CA models, in which length (l₀) and timescales (t₀) are unknown, to real units, thus allowing for intercomparison between simulations and data collected from field experiments. We present their approach and the related calculations of length and timescales used in this study below.

Consider a system consisting of air blowing over a flat granular bed. If the flow velocity is faster than a threshold velocity, the system is mechanically unstable, resulting in aeolian transport and leading to the formation of the bedforms. For a small period of time after aeolian transport commences, the system can be analysed using linear stability analysis, which identifies the fastest growing mode of the evolving surface. Past theoretical analyses (Hersen et al., 2002; Elbelrhiti et al., 2005; Claudin and Andreotti, 2006) have established this length scale, λ_max, as being equal to 50 ${\mathit{\rho }}_{\mathrm{s}}\stackrel{\mathrm{‾}}{d}/{\mathit{\rho }}_{\mathrm{f}}$, where ρ_s is the density of the grains, ρ_f is the density of the overlying fluid (air in our case) and $\stackrel{\mathrm{‾}}{d}$ is the mean grain size diameter.

We performed a series of numerical experiments using the ReSCAL model where we simulate flow over a wavy surface with a different wavelength for each experiment and with an amplitude (A) of 2 h₀, which is the smallest amplitude we can have in our discrete system. We allow the simulation to proceed for a very short time and measure the growth of the natural logarithm of amplitude with time (dln (A)∕dt). The results of these numerical experiments are shown in Fig. 2a. It can be clearly seen that the fastest-growing wavelength is λ=28l₀. Thus,

Using values of ρ_f=1.00 kg m⁻³, ρ_s=910 kg m⁻³ and $\stackrel{\mathrm{‾}}{d}=\mathrm{200}$ µm, we find the length scale of our model to be l₀=0.325 m.

Once the length scale has been identified, we can now proceed to identify the timescale of the system. Returning to the system of air being blown over a flat surface, it has already been established that a long time after the initiation of the aeolian transport and with the wind speed being constant, the flux of grains in the air achieves a steady-state value known as the saturated flux Q_sat. Past work, beginning already with Bagnold (1941) and refined over successive studies using both field and wind tunnel data, has resulted in a semi-empirical formulation of Q_sat as a function of material and flow properties (Bagnold, 1936; Iversen and Rasmussen, 1999; Ungar and Haff, 1987). In effect, the saturated flux can be computed as

with

where u_* and u_c are the friction velocity and the threshold friction velocity for aeolian transport respectively.

If we consider an idealized scenario where the threshold velocity is zero with the resultant saturated flux value being $Q_{\mathrm{sat}}^{\mathrm{0}}$, a relationship relating saturated fluxes only as a function of u_* and u_c can be found as

This relationship is quite useful for CA-based modelling as the free parameter of τ₁ is essentially equivalent to $u_{\mathrm{c}}^{\mathrm{2}}$ and for modelling purposes can be chosen to be equal to zero. Thus, for different values of τ₁, the CA model provides the left-hand side (LHS) of Eq. (6). We performed a series of experiments beginning with the idealized scenario of setting τ₁=0 and gradually increasing the value of τ₁ for each individual experiment. In each experiment we allowed the system to reach steady state and calculated the steady-state flux of snow in the air. The resulting values of $Q_{\mathrm{sat}}^{\mathrm{ratio}}$ as a function of τ₁ are shown in Fig. 2b. Once $Q_{\mathrm{sat}}^{\mathrm{ratio}}$ ratios for different values of τ₁ are found, we can compute the u_* in real units for a given value of τ₁. Note that u_c is computed using Eq. (5b). Using the log law, a given u_* value can be converted to wind speed above a chosen height over the surface. Using a roughness length of $z_{\mathrm{0}}={\mathrm{10}}^{-\mathrm{4}}$ m, wind speed at a height of 1 m above the surface, denoted as U_1 m, is computed and shown in Fig. 2c. Once the u_* value for each τ₁ is found, we can compute the real saturated flux using Eq. (5a) and equate it to the model saturated flux. Thus $Q_{\mathrm{sat}}^{\mathrm{model}}{l}_{\mathrm{0}}{h}_{\mathrm{0}}/t_{\mathrm{0}}=Q_{\mathrm{sat}}^{\mathrm{real}}$ m⁻² s⁻¹. As the length scales l₀ and h₀ are known, t₀ can be found as a function of τ₁. Values of the timescale t₀ for different values of τ₁ are shown in Fig. 2d. Note that the stress scale, τ₀, is not identified explicitly due to the fact that in addition to length and timescales, it would require a mass scale as well. Furthermore, since it is important only in terms of ratios and is directly related to the wind speed, it is not of major importance to the current study.

Having established the length and timescales of the CA model, we choose three particular values of τ₁ that are typically used in this study. These are ${\mathit{\tau }}_{\mathrm{1}}\in \left\{\mathrm{5},\mathrm{20},\mathrm{60}\right\}$ and are equivalent to $U_{\mathrm{1}\mathrm{m}}\in \left\{\mathrm{20.5},\mathrm{12.5},\mathrm{7.0}\right\}$ m s⁻¹ and denoted further as high-wind (UH), medium-wind (UM) and low-wind (UL) cases respectively. Values of different quantities involved in calculating the timescales for these three τ₁ values are provided in Table 2 with the full dataset of computed values for all the values of τ₁ provided in Table S1 in the Supplement.

Figure 2Establishing the length and timescales of the CA model. (a) Identifying the more unstable wavelength (see Eq. 4). λ_max=28l₀. (b) Variation in the $Q_{\mathrm{sat}}^{\mathrm{ratio}}$ with τ₁. (c) Variation in the velocity at 1 m above the surface U_1 m as a function of τ₁. (d) Identifying the timescale t₀ of the model for different values of τ₁

Download

Table 2Details of the chosen wind scenarios for further analyses along with calculation of different quantities leading up to finding the relevant timescales for each wind speed.

^a Computed using $\frac{Q_{\mathrm{sat}}\left({\mathit{\tau }}_{\mathrm{1}}\right)}{Q_{\mathrm{sat}}^{\mathrm{0}}}$, where $Q_{\mathrm{sat}}^{\mathrm{0}}=\mathrm{2.273}\times {\mathrm{10}}^{-\mathrm{2}}{l}_{\mathrm{0}}{h}_{\mathrm{0}}/t_{\mathrm{0}}$. ^b Using Eq. (6), where u_c=0.134 m s⁻¹, using Eq. (5b). ^c Through the log law, $u=\frac{u_{*}}{\mathit{\kappa }}log\left(\frac{z}{z_{\mathrm{0}}}\right)$, where κ is the von Kármán constant (=0.4) and z₀ is the roughness length $(={\mathrm{10}}^{-\mathrm{4}})$. ^d Computed using Eq. (5a) using material properties as described in the text. ^e By equating model and real saturated flux.

Download Print Version | Download XLSX

3.1 Steady-state barchan motion

Under the influence of constantly blowing wind, a conical pile morphs into a barchan dune that starts to move downstream. In Fig. 3a, we show the top and side views of the evolution of the conical pile (case C20, with a constant wind speed of U_1 m=20.5 m s⁻¹), into a barchan dune. Note that the time is scaled using the sintering timescale (t_s) of 24 h even though the sintering mechanism itself is switched off. The cone essentially flattens and elongates in both the streamwise (along-wind) and crosswind directions, thus increasing the length and the width of the snow deposit. This evolution can be quantitatively seen in Fig. 3b. The quantities of length (L) and width (W) are shown on the left axis, whereas the height (H) of the cone is shown on the right axis. The length and width of the dune is initially the same as the diameter of the cone, in this case, approximately 15 m, while the height of the cone is initially 1.1 m. After approximately 3 t_s, the morphology of the barchan, particularly its height, is approximately constant, and thus we consider that steady state has been achieved. At steady state, the L, W and H dimensions of the C20 dune are respectively 26.4 m, 22.75 m and 0.59 m.

We additionally also show the evolution of the length of the longest streamwise section of the barchan. This quantity is termed as the maximum streamwise length or MSL of the bedform (L_s). It has been shown that this is the slowest moving part of a barchan and thus is representative of the speed of the dune (Zhang et al., 2014). At steady state, we find L_s to be equal to 12.18 m. The relevance of L_s will be more apparent in the coming paragraphs. It is interesting to note that it takes approximately 3 t_s or 72 h for this dune to reach steady state. This is quite a long time considering that we have the wind blowing constantly for this period at U_1 m=20.5 m s⁻¹. Thus morphodynamics of the barchans are much slower than typical sintering timescales. This point is elaborated upon further in the rest of the paper.

In Fig. 3c, we show four different measures of the speed of the dune. In the field, one would typically track the location of the crest of the dune as a function of time. Other measures are also possible, such as tracking the displacement of the horn or the tail of the dune. Locations of the crest, tail and the horn of the barchan are tracked and plotted as a function of time. We find that all the speed measures are quite noisy, with the horn-based measure being most noisy and thus one to be avoided. To find the speeds, we thus resort to finding the slope of the best-fit line, values of which are provided in the legend of Fig. 3c. The horn-based speed measure is found to be the fastest, whereas the crest-based measure is found to be the lowest. Admittedly, the differences are minor, especially compared to the fluctuations itself. We use an additional measure based on tracking the centre of mass (COM) of the dune. In the field, this measure may be obtained by assuming constant snowpack density and a laser scanner. This measure may be preferred as a more physically based and unbiased measure of dune movement and thus in the rest of this study dune speed is meant to be the speed of the COM of a dune.

Figure 3Morphodynamics of a solitary dune (case C20, high-wind scenario U_1 m=20.5 m s⁻¹). (a) Visual representation of evolution from a cone pile to a barchan. The final image is annotated for identifying different descriptors of a barchan. (b) Evolution of the length (L), width (W) and the maximum streamwise length (L_s) of the barchan. (c) Different versions of calculating the dune speed using the displacement of the centre of mass (COM), tail, horns or crest of the dune.

Download

Figure 4Influence of (a) dune length and (b) wind speed on dune velocity. The instantaneous speeds are represented by symbols while the time-averaged speed is shown by thick lines. The wind speed in (a) is the UH scenario whereas the dune in (b) is the cone C24.

Download

An example of the cone-based experiments presented in the previous paragraphs and shown in Fig. 3 is repeated for 16 different cone (and thus barchan) sizes and two additional wind speeds, namely the medium-wind UM (U_1 m=12.5 m s⁻¹) and low-wind UL(U_1 m=7.0 m s⁻¹) scenarios. In Fig. 4a, we show the variation in dune speed with dune length for four different dunes for the UH wind scenario. We show both the instantaneous dune speed (computed every 50 t₀ steps, represented by circular symbols) and the time-averaged velocity trends. This is to contrast the large instantaneous dune speed fluctuations with a comparatively constrained time-averaged value of dune speed. As expected, dune speed decreases with the length (and thus size) of the dune. In Fig. 4b we show the variation in the dune speed of the same dune for the three different wind scenarios. Similar to the previous figure, there are large fluctuations of dune speed while the time-averaged value is more constrained. While differences between UH and UM cases are not significantly different, the dune in the UL case is almost 20 % slower than that in the UH and the UM cases. It must be noted that the effect of wind speed is not as significant as the dune size. Increasing the wind speed by a factor of 3 between the UL and the UH cases does not seem to induce a proportional response from the dune speed.

Thus, we have shown that the dune speed is inversely proportional to its height and directly proportional to wind speed. This is in fact a well-known property of barchans first recognized by Bagnold (1941) and quantitatively explored using field data by Elbelrhiti et al. (2005) for sand and by Kobayashi and Ishida (1979) for snow and in CA-type numerical experiments by Zhang et al. (2014). All these studies roughly proposed that

where c is the dune speed, Q is the saturated snow flux, c^* is the normalized dune speed, H is the height of the dune and f is a linear function. We explore this scaling for snow barchans in Fig. 5.

Long time averaging of both the height and velocity of the dune in each of the 48 simulations is carried out and the velocity-versus-height data points for all these simulations are shown in Fig. 5a. Simulations are classified based on wind speed alone and coloured accordingly. Note that velocity is normalized by the flux of snow in each simulation. We find that all the simulation results, for different barchan dimensions as well as wind speeds, lie on a hyperbolic function of H. The least-squares fit is found to be

where a, b and d are parameters with values of 1.7 (dimensionless), −0.1 m and 0.94 m⁻¹ respectively. In the same figure, we also show the instantaneous values of normalized velocity as a function of height for each of the simulations. To further quantify the fluctuations in instantaneous velocity, we show histograms of this quantity for four different barchans for U_1 m=20.5 m s⁻¹ in Fig. 5b. We find that dune speeds become more constrained with increasing dune size with the smallest (largest) dune having the most (least) broad histogram.

Figure 5(a) Dune speed normalized by snow flux (c∕Q) as a function of height (H) for all the solitary dune simulations performed. Note that the instantaneous speeds are presented by small lightly coloured symbols. Time-averaged dune speeds for each simulation are shown by large symbols. Note that all symbols are coloured according to the wind speed of the simulation. (b) Probability distribution function of dune speeds for barchans of four different heights.

Download

We needed to perform time averaging for at least 18 h of instantaneous velocities to converge the time-averaged dune speed with the smallest dunes requiring up to 36 h of averaging. While this is fairly simple to do in an idealized environment of numerical computations, it is highly unlikely that such time averaging will be possible in the field, primarily due to fluctuations of mean wind speeds, effects of topography and the fact that the sintering process has similar timescales. On the other hand, due to large fluctuations of instantaneous velocity, a limited time series, even of a few hours, is unlikely to show any systematic trend. A possible solution could be to sample multiple mobile dunes, hopefully of different sizes, at the same time using a laser scanner or photogrammetry.

We would like to note that values of time-averaged dune speeds found using the CA model are fairly similar to those reported in (a few) published studies (see for example, Table 2 in Filhol and Sturm, 2015). The values for the UH and the UM cases in particular are close to the studies by Doumani (1967), Kuznetsov (1960) and Kotlyakov (1961). It must be noted that these studies are quite old, and in future work we shall further compare our results with the latest dataset from Kochanski (2018) and Kochanski et al. (2018). Dune dimensions and speeds for all 48 simulations performed are provided for reference in Table S2.

3.2 Effect of sintering on barchan motion

We now turn our attention to understanding the effect of sintering on dune morphodynamics. In this section, we focus on the effect of sintering on barchans that are already in steady motion. This is achieved by “switching on” the effect of sintering only once the morphology and the mean dune speed has reached a steady state.

Sintering is activated in all 48 simulations with different barchan shapes and wind speeds. In general, three types of behaviours are observed. Firstly, the fastest-moving barchans continue their motion without much difference. On the other hand, the slowest-moving barchans seem to sinter “in place”; i.e, the barchan immediately ceases to move, without any significant change in morphology. The intermediate behaviour is found for a range of dune speeds in between the extreme cases where a small part of the dune is deposited as a non-erodible layer while the dune continues to move, albeit with slightly reduced dimensions on account of loss of mass due to sintering.

Figure 6Effect of sintering on the morphodynamics of a mobile barchan (case C20): (a) medium-wind (U_1 m=12.5 m s⁻¹) and (b) low-wind (U_1 m=7.0 m s⁻¹) scenarios. Note that the colour scheme is such that light blue colours represent mobile (eDS-type) cells while darker shades of blue are sintered (neDS-type) cells.

Download

To illustrate these three different types of behaviours, we show in Fig. 6 the morphodynamics of the same barchan dune (C20: L, W, H, L_s=26.4, 22.75, 0.59 and 12.3 m) for the UM and UL cases. The UH case is not shown because no perceptible difference in the morphodynamics is detected. For the UM case, shown in Fig. 6a, we find that the dune begins to leave behind a mass of snow that can no longer be eroded. This mass originates at the tail end of the dune, which is the oldest part of the dune as we shall see further. The dune velocity is such that there is a continuous ejection of mass from the tail as the dune continues its downwind motion. Ultimately, there is a split, where a sintered non-erodible mass of snow is left behind while a smaller barchan remains intact and mobile. In the UL case, the barchan sinters in place, and very quickly comes to a standstill. Notice that the shape and the dimensions of the barchan in this case do not change much once the sintering is activated. The morphodynamics of the two cases presented in Fig. 6 along with two additional cases are provided as Movies M1–M4 in the Supplement.

The reason for this behaviour becomes clear when we analyse the distribution of ages of the constituent cells of a dune. Recall that “age” is defined as the time since a cell transitioned to the deposited snow (DS) type and has remained as such. In Fig. 7a, the three panels show the distribution of ages on the central slice of a C24 dune for three different wind speeds with the UH and the UL cases consisting of the youngest and oldest barchans respectively. The age increases when moving from the leeward to the windward face of the barchan. The distribution in the flanks of the dune interior is also shown and it is found to be an extension of the distribution found in the central slice with no abrupt variations. It is interesting to note the stratigraphy of the dune and the distinct bands of each new layer added on the leeward face of the barchan. While the distribution indeed changes as a function of the wind speed, the banded pattern is consistent in all three panels. Furthermore, the banded pattern also highlights the fact that deposition on the leeward side of the dune occurs in pulses rather than in a continuous fashion. This is correlated to the large fluctuations to dune speeds at short timescales as described earlier in Figs. 4 and 5.

The adjoining Fig. 7b quantifies the age distribution of the three panels in the form of cumulative distribution functions (CDFs) of ages with respect to wind speeds. For the UH case, we find that all of the barchan is younger than 0.3 t_s, and thus sintering does not have a perceptible impact on its dynamics. On the other hand, for the UM case, even though most of the barchan is younger than the sintering timescale of 24 h, approximately 50 % of the barchan is older than 0.3 t_s. Thus, there is an increasing influence of sintering on the dynamics of this barchan, evidenced previously in Fig. 6a. Finally, in the UL case, almost 50 % of the barchan has ages greater than the sintering time! Thus when the effect of sintering is activated, the barchan almost immediately ceases its motion and comes to a halt, thereby sintering in place.

Figure 7Distribution of age within a mobile barchan (case C24) prior to sintering. (a) Three panels show the distribution of age in the central slice as well as the arms of the barchan for the three different wind scenarios. (b) Cumulative distribution function of the age of the entire barchan for the three different wind scenarios.

Download

Recall that for each wind speed (UL,UM and UH), we performed simulations for 16 different barchans with increasing dimensions. These simulations allow us to identify the size of the barchan at which the effect of sintering is strongly felt. To do so, in Fig. 8a, we show the location of the tail of the different barchans once the sintering is activated for the UH case. The slope of each of these lines would correspond to the tail-based speed measure shown previously in Fig. 3c. With increasing size, the barchan speed decreases as previously discussed. Lines are coloured green to indicate that the barchan speed does not change due to sintering. However, at a certain barchan size, there is a transition where, upon reaching sintering time t_s, the tail no longer moves; i.e, the motion of barchan has been perturbed by sintering. We identify that for the UH case, this occurs for barchan with MSL of L_s=21.5 m (line coloured black); all barchans with sizes greater than this limit are also affected by sintering and show behaviour similar to that shown in Fig. 6b. A similar analysis for the UM case shows the limit to be at a barchan size of L_s=12.3 m. In the UL case (not shown) all the chosen barchan shapes sinter in place and thus the MSL is less than 8.8 m for this case.

These numerical experiments highlight the fact that the interplay between barchan dimensions, wind speed and the sintering rate impose a maximum length scale of any snow bedform that can remain erodible and thus mobile. For example, in the UH case (U_1 m=20.5 m s⁻¹), the largest streamwise dimension of any eDS-type bedform is limited to 21.5 m. Any bedform with dimensions greater than this limit (which perhaps arose earlier due to even higher wind speeds) will be strongly affected by sintering and will most likely sinter in place, thus converting to a neDS-type bedform. This limit decreases with wind speed and increases with the sintering rate.

In the following section, we move to a more realistic case where instead of beginning with a cone pile or even a barchan in steady motion, we directly simulate the effect of wind blowing over an initially flat fresh snow layer. As described earlier, the flat snow surface is mechanically unstable and rapidly evolves into various bedforms. Interestingly, we find that the results presented in this section remain valid there as well.

Figure 8Identification of the maximum streamwise length (MSL, L_s) for (a) high-wind and (b) medium-wind scenarios. In each figure, the location of the tail of the barchan is plotted as a function of time. Lines from top to bottom represent barchans with increasing size. The lines are coloured to identify mobile (green) and immobilized (red) behaviours. The black line identifies the barchan at which this transition occurs.

Download