The probabilistic nature of dune collisions in 2D

Dunes are bedforms of different size and shape, appearing throughout aeolian, subaqueous and extraterrestrial environments. Collisions between dunes drive dune field evolution, and are a direct result of interacting dunes of different heights, travelling at different speeds. We perform 2D cellular automaton simulations of collisions between dune pairs migrating in a steady flow. Modelled collisions can result in either ejection, where dunes exchange mass before separating, or downstream- or upstream-dominant coalescence (merging of dunes). For each of these three elementary types of interaction, we identify the mass exchange mechanism and the distinctive intermediate morphologies. Surprisingly, we show that the collision outcome depends probabilistically on the initial dune area ratio r and can be described by a narrow sigmoidal function centred on $r=\mathrm{1}/\mathrm{2}$. Finally, we compare our simulations with laboratory experiments of dune collisions, finding good agreement concerning the intermediate morphology and the collision outcome. Our results can motivate further observational or experimental studies that validate our probabilistic collision predictions and fully determine the controls on the coalescence–ejection transition.

1 Introduction

Dunes are self-organising structures that form spontaneously on a particle-laden surface overlain by a sufficiently strong flowing medium. In nature, they can be found in aeolian landscapes, e.g. deserts (Elbelrhiti et al., 2008; Lü et al., 2021) or coastal beaches (Parteli et al., 2006), aqueous environments such as river beds (de Almeida et al., 2016) or subjected to extraterrestrial atmospheres (Bishop, 2007). In these settings, dunes have long been thought to grow to attain a maximum height which, in settings with sufficient sediment, could be controlled by the flow depth (in aqueous environments) or the thickness of the planetary boundary layer (for aeolian systems; Andreotti et al., 2009; Andreotti and Claudin, 2013). Alternatively, sediment supply may be an additional limiting factor (Gunn et al., 2022; Jarvis et al., 2022). Regardless, the migration velocity c of a dune depends inversely on the dune size. Based on mass conservation, it has frequently been proposed that $c\sim \mathrm{1}/\mathcal{H}$ (Bagnold, 1941; Southard, 1991), where ℋ is the dune height, but other relations, including $c\propto \mathrm{1}/\mathcal{L}$ (Kroy et al., 2002), where ℒ is the dune length, or $c\propto \mathrm{1}/(\mathcal{H}+{\mathcal{H}}_{\mathrm{0}})$ (Andreotti et al., 2002a), where ℋ₀ is a constant, have also been suggested.

The inverse relationship between dune size and velocity means that faster, smaller dunes can approach and collide with slower, larger dunes. Such collisions have been frequently observed in subaqueous experiments (Coleman and Melville, 1994; Bradley and Venditti, 2019; Jarvis et al., 2022) and numerical simulations (Gao et al., 2015), as well as for aeolian dunes (Lü et al., 2021), even though the latter evolve over much longer timescales. These collisions have been shown to play an important role in pattern coarsening, whereby a larger number of smaller dunes transition to become a smaller number of larger dunes (Coleman and Melville, 1994; Gao et al., 2015; Bradley and Venditti, 2019; Jarvis et al., 2022). This coarsening process occurs during the development of a dune field from a flat bed, whereas more mature dune fields can attain a steady state, where their mean wavelength and amplitude remain relatively constant. In 2D, coarsening primarily occurs through dune coalescence, where a pair of colliding dunes merge. Two distinct types of coalescence have been observed, downstream and upstream dominant, defined according to which peak survives the coalescence process (Coleman and Melville, 1994; Gao et al., 2015; Jarvis et al., 2022). Another possible collision outcome that has been recognized experimentally is ejection, where mass is transferred from the downstream larger, slower dune to the upstream smaller, faster dune until the leading dune becomes sufficiently small that it can accelerate away (Diniega et al., 2010; Gao et al., 2015; Jarvis et al., 2022). The balance between the different outcomes ultimately means that dune collisions are a strong and dominant control on dune field evolution (Elbelrhiti et al., 2008; Durán et al., 2009; Kocurek et al., 2010; Hugenholtz and Barchyn, 2012), regulating both the size and the spacing of dunes (Hersen and Douady, 2005; Génois et al., 2013a, b).

Although this current study only considers the collision of 2D dunes, for which only coalescence and ejection have been observed, it is important to acknowledge that collisions between 3D barchan dunes can produce a much wider range of outcomes owing to turbulence and lateral sediment transport (Endo et al., 2004; Durán et al., 2005; Elbelrhiti et al., 2005; Hersen, 2005; Katsuki et al., 2011; Assis and Franklin, 2020, 2021). Additionally, it has also been shown that a turbulent wake shed by an upstream dune (Bristow et al., 2018, 2019, 2020) can enhance the migration velocity of a downstream neighbour or even cause it to split into two (Elbelrhiti et al., 2005; Worman et al., 2013), further complicating collision dynamics (Bacik et al., 2020, 2021; Assis and Franklin, 2021). Although recirculation zones downstream of dune crests occur in both 2D and 3D configurations (Hermann et al., 2005; Araùjo et al., 2013; Michelsen et al., 2015), turbulence itself is an inherently 3D phenomenon, even present in quasi-2D experiments (Bacik et al., 2020, 2021). Thus, in a pure 2D domain and for dunes sufficiently large to be scale invariant (Kroy et al., 2002; Andreotti et al., 2002b), only the size ratio between dunes controls the collision outcome. Katsuki et al. (2005) developed a simplified model of 2D dune collision and showed that for $r=A_{\mathrm{u}}/A_{\mathrm{d}}\gtrsim \mathrm{0.5}$, where A_u and A_d are the areas of the upstream and downstream dunes respectively, ejection should occur. Conversely, coalescence is predicted for r≲0.5. A more sophisticated continuum model from Diniega et al. (2010) which accounts for the dependence of sediment flux on shear stress, however, predicted the transition to occur for $r\approx \mathrm{1}/\mathrm{3}$. For 3D barchans, both Katsuki et al. (2011) and Bo and Zheng (2013) found the transition also depended on the perpendicular distance between the dune axes (axis offset distance). Later, Zhou et al. (2019), motivated by the results of large eddy simulations of on-axis collisions between 3D barchans, suggested ejection occurs if the sand flux received by the downstream dune exceeds what is lost from the barchan horns. These predictions and hypotheses remain largely unverified, either through field observations or laboratory experiments. However, Assis and Franklin (2020) showed experimentally that, for subaqueous 3D barchans, all collision outcomes can be mapped in a regime space defined by the Shields number (proxy for flow strength), the axis offset distance and the dune size ratio. Such a regime diagram has not been experimentally created in a 2D system though.

In this paper, we use the cellular automaton model ReSCAL (Narteau et al., 2009; Zhang et al., 2010; Rozier and Narteau, 2013) to simulate the collision of 2D dunes and quantitatively constrain the coalescence–ejection transition. We are able to produce both downstream- and upstream-dominant coalescence outcomes, as well as instances of ejection. Such interactions have previously been generated in ReSCAL simulations of bedform coarsening (Gao et al., 2015), although this study represents a first attempt to constrain the transition between regimes. Limiting our domain to 2D results in a narrower range of outcomes than can be observed in natural 3D systems. However, this means we can more completely study the phenomena, something that would be computationally very expensive in the much larger parameter space of 3D systems. By performing a large number (1600) of simulations for dune pairs with the same size ratio, we show that the collision outcome can be modelled probabilistically rather than deterministically, and we find an empirical relationship for the coalescence–ejection transition. Additionally, we note that intermediate structures during collisions are associated with distinct morphologies. Finally, we compare these results with qualitative and quantitative observations of colliding dunes in subaqueous experiments reported by Jarvis et al. (2022). Although our study is strictly only valid for 2D systems, our results should motivate further research to test if the fundamental characteristics of collisions that we observe can be preserved in 3D environments, where turbulence and flow perturbations, as well as lateral sediment transport during interactions, also have an influence.

2 Methods

We simulate interactions between discrete 2D dunes using the dune model ReSCAL (Narteau et al., 2009), which couples a cellular automaton model of sediment transport and a lattice gas model of turbulent fluid flow. We summarise the numerical method below, whilst full details can be found in Narteau et al. (2009) and Rozier and Narteau (2013).

Simulations consider a two-dimensional domain of square cells of length l₀ in one of the following states: fluid, neutral, immobile sediment or mobile sediment. Neutral cells denote the upper and lower boundaries. Sediment transport is modelled using transitions of pairs of nearest-neighbour cells corresponding to physical processes (erosion, deposition, transport). Additionally, a repose angle of the sediment is imposed to account for avalanching. Whilst direct numerical simulations (Lefebvre and Winter, 2016) and experiments (Kwoll et al., 2016) on subaqueous dunes have shown that the leeside slope angle strongly controls the flow downstream of the dune, with low-angle dunes significantly reducing turbulence and flow separation, we fix the angle of repose at 35^∘. We briefly discuss the implications of this choice in Sect. 4. Each process has a characteristic timescale expressed in units of t₀ (the model time step), which determines the probability of a particular transition to occur. This probabilistic approach means the model is not entirely deterministic, but takes a random seed as an input.

The fluid flow is simulated using a lattice gas model. Particles move according to their velocity vectors that can change according to collisions. Particle fluxes are averaged to produce a velocity field. Particles colliding with the upper boundary conserve their horizontal velocity but experience a change of sign to their vertical velocity, creating a free-slip boundary condition. At the lower boundary, all fluid particles colliding with the surface rebound in their incident direction, creating a no-slip condition. This reproduces the expected logarithmic velocity profile for turbulent flow over a flat wall (Narteau et al., 2009), with surface shear stress τ_s defined as the normal derivative of the velocity field with respect to topography. The erosion rate is defined for three distinct regions as

$$\begin{array}{}\text{(1)}& {\mathrm{\Lambda }}_{\mathrm{e}}=\left\{\begin{array}{ll}\mathrm{0},& \mathrm{for}{\mathit{\tau }}_{\mathrm{s}}<{\mathit{\tau }}_{\mathrm{1}},\\ {\mathrm{\Lambda }}_{\mathrm{0}}\left({\mathit{\tau }}_{\mathrm{s}}-{\mathit{\tau }}_{\mathrm{1}}\right)/\left({\mathit{\tau }}_{\mathrm{2}}-{\mathit{\tau }}_{\mathrm{1}}\right),& \mathrm{for}{\mathit{\tau }}_{\mathrm{1}}\le {\mathit{\tau }}_{\mathrm{s}}\le {\mathit{\tau }}_{\mathrm{2}},\\ {\mathrm{\Lambda }}_{\mathrm{0}},& \mathrm{for}{\mathit{\tau }}_{\mathrm{s}}>{\mathit{\tau }}_{\mathrm{2}},\end{array}\right.\end{array}$$

where Λ₀ is a constant, τ₁ the critical shear stress for sediment motion and $({\mathit{\tau }}_{\mathrm{2}}-{\mathit{\tau }}_{\mathrm{1}}{)}^{-\mathrm{1}}$ the linear coefficient between Λ_e and τ (τ₂ is an adjustable parameter). This erosion, along with transport and deposition, modifies the lower boundary of the fluid domain, providing feedback on the fluid flow. This two-way coupling generates acceleration of the flow on the upstream side of dunes, shear layers and a recirculation zone on the downstream side (Narteau et al., 2009; Zhang et al., 2010), as observed naturally (Sweet and Kocurek, 1990) or imposed by other numerical models (Hermann et al., 2005).

We use a two-dimensional periodic domain of height H=300l₀ and length L=2000l₀ (see Appendix A for justification that the domain is sufficiently large). Two triangular sand piles with slope angle 35^∘ (consistent with the angle of repose used in previous cellular automaton simulations, e.g. Zhang et al., 2014) were placed 1000l₀ apart. The slope angle was equal to the repose angle. In all simulations, the downstream dune was initially larger than that upstream, ensuring interaction. We set τ₁=0 and ${\mathit{\tau }}_{\mathrm{2}}-{\mathit{\tau }}_{\mathrm{1}}=\mathrm{1000}{\mathit{\tau }}_{\mathrm{0}}$ but found that the observed interactions did not depend on the absolute values for τ₂≫τ₁, i.e. as long as the flow regime is far from the sediment transport threshold. This is to be expected since we are assuming flow conditions are far above the threshold for sediment transport. Table B1 lists the parameter values used for the suite of simulations.

2.1 Physical scaling

Later in this article (Sect. 4), we will compare our simulated collisions with those observed in the subaqueous experiments of Jarvis et al. (2022). However, since cellular automatons are defined on a discrete domain, and given the arbitrary nature of the transition rules, there is no prior relationship between l₀ and the time step t₀ on one hand and physical length and timescales on the other. Instead, we define these by comparing model results with experimental and natural observations, following Narteau et al. (2009) and Zhang et al. (2014). For l₀, this was previously achieved (Narteau et al., 2009; Zhang et al., 2014) by simulating fluid flow over an initially flat sediment bed and determining the most unstable initial wavelength λ_max, which is predicted to be λ_max=50 $({\mathit{\rho }}_{\mathrm{p}}-{\mathit{\rho }}_{\mathrm{f}})d/{\mathit{\rho }}_{\mathrm{f}}$ (Hersen et al., 2002; Elbelrhiti et al., 2005; Andreotti and Claudin, 2013), where ρ_p(f) is the sediment (fluid) density and d the particle diameter. Within ReSCAL, λ_max≈40l₀ (Narteau et al., 2009; Zhang et al., 2014), whilst in the experiments of Jarvis et al. (2022), ρ_p=2500 kg m⁻³, ρ_f=1000 kg m⁻³ and $d=\mathrm{1.21}\times {\mathrm{10}}^{-\mathrm{3}}$ m, resulting in $l_{\mathrm{0}}=\mathrm{2.27}\times {\mathrm{10}}^{-\mathrm{3}}$ m. We then determine the timescale t₀ by comparing the modelled equilibrium sediment flux with field observations (Narteau et al., 2009; Zhang et al., 2014). For the saturated flux, we take (Meyer-Peter and Müller, 1948)

$$\begin{array}{}\text{(2)}& Q_{\mathrm{sat}}={\displaystyle \frac{\mathrm{8}{\mathit{\rho }}_{\mathrm{f}}{\left(u_{*}^{\mathrm{2}}-u_{*\mathrm{t}}^{\mathrm{2}}\right)}^{\mathrm{3}/\mathrm{2}}}{\left({\mathit{\rho }}_{\mathrm{p}}-{\mathit{\rho }}_{\mathrm{f}}\right)g}}.\end{array}$$

We note that many similar transport laws exist for bedload transport in turbulent flow (Paintal, 1971; Bagnold, 1973; Engelund and Fredsoe, 1976; Fernandez-Luque and Beek, 1976; Bridge and Dominic, 1984; Lajeunesse et al., 2010; Pähtz and Durán, 2020), whose mathematical form only varies slightly. Since there is little experimental evidence distinguishing between these transport laws (Lajeunesse et al., 2010), we choose for simplicity Eq. (2). We also take (Bagnold, 1936; Iverson and Rasmussen, 1999)

$$\begin{array}{}\text{(3)}& u_{*\mathrm{t}}={\displaystyle \frac{\mathrm{1}}{\mathrm{10}}}\sqrt{{\displaystyle \frac{\left({\mathit{\rho }}_{\mathrm{p}}-{\mathit{\rho }}_{\mathrm{f}}\right)gd}{{\mathit{\rho }}_{\mathrm{f}}}}}\end{array}$$

for the threshold friction velocity. Substituting Eq. (3) into Eq. (2) gives

$$\begin{array}{}\text{(4)}& Q_{\mathrm{sat}}={\displaystyle \frac{\mathrm{8}{\mathit{\rho }}_{\mathrm{f}}}{\left({\mathit{\rho }}_{\mathrm{p}}-{\mathit{\rho }}_{\mathrm{f}}\right)g}}{\left(u_{*}^{\mathrm{2}}-{\displaystyle \frac{\left({\mathit{\rho }}_{\mathrm{p}}-{\mathit{\rho }}_{\mathrm{f}}\right)gd}{\mathrm{100}{\mathit{\rho }}_{\mathrm{f}}}}\right)}^{\mathrm{3}/\mathrm{2}}.\end{array}$$

In the experiments of Jarvis et al. (2022), u_* varies, so we take an intermediate value of 0.1 m s⁻¹. Equation (4) then gives $Q_{\mathrm{sat}}=\mathrm{5.29}\times {\mathrm{10}}^{-\mathrm{4}}$ m² s⁻¹, which is matched with the modelled saturated flux (Narteau et al., 2009; Zhang et al., 2014) to determine that $t_{\mathrm{0}}=\mathrm{9.74}\times {\mathrm{10}}^{-\mathrm{3}}$ s. The computed values of l₀ and t₀ enable comparisons between our simulations and the experiments of Jarvis et al. (2022).

3 Results and discussion

3.1 Three elementary types of dune interaction

The simulations reproduce both coalescence and ejection behaviour, as shown in Fig. 1, where the initially upstream dune has been coloured red and the initially downstream dune blue. Two distinct types of coalescence (Coleman and Melville, 1994; Gao et al., 2015), downstream- (Fig. 1a) and upstream-dominant (Fig. 1b) coalescence, feature distinctly different intermediate mechanics despite both resulting in only a single dune surviving the interaction. Figure 2 shows more details regarding the evolution of the positions and heights of the peaks and troughs in the intermediate structures.

Figure 1 The three elementary types of dune interaction observed in the simulations: (a) downstream-dominant coalescence ($A_{\mathrm{u}}=\mathrm{300}{l}_{\mathrm{0}}^{\mathrm{2}}$, $A_{\mathrm{d}}=\mathrm{3000}{l}_{\mathrm{0}}^{\mathrm{2}}$, r=0.1), (b) upstream-dominant coalescence ($A_{\mathrm{u}}=\mathrm{1200}{l}_{\mathrm{0}}^{\mathrm{2}}$, $A_{\mathrm{d}}=\mathrm{3000}{l}_{\mathrm{0}}^{\mathrm{2}}$, r=0.3) and (c) ejection ($A_{\mathrm{u}}=\mathrm{1800}{l}_{\mathrm{0}}^{\mathrm{2}}$, $A_{\mathrm{d}}=\mathrm{3000}{l}_{\mathrm{0}}^{\mathrm{2}}$, r=0.6).

Figure 2 Peak and trough positions and heights in the intermediate bedform during (a, d) downstream-dominant and (b, e) upstream-dominant coalescence and (c, f) ejection. t=0 is defined as when the dunes first touch. The trough is defined as the range of points which form the minimum height in the intermediate bedform, and thus the extent is finite (see sketch). In (a, d), downstream-dominant coalescence completes once the trough and upstream peak coincide and is represented by dotted line segments.

Downstream-dominant coalescence (Fig. 1a) prevails if the upstream dune is much smaller than the downstream dune ($r=A_{\mathrm{u}}/A_{\mathrm{d}}\ll \mathrm{1}$). The upstream dune climbs the stoss side of the downstream dune (raising the trough), but the increased shear stress on the stoss slope of the larger dune spreads the upstream material, some of which is transported over the downstream crest and deposited on the lee side. The trough increases in height faster than the upstream peak until they reach the same height and the upstream peak vanishes. The sediment which initially constituted the upstream dune (red in Fig. 1) ends up in a layer parallel to the slip face of the resultant dune. Migrating dunes typically contain many slip face-parallel layers (Bagnold, 1941; Allen, 1970) and, since the upstream material becomes incorporated in such a layer, evidence of the interaction is lost over time.

Upstream-dominant coalescence (Fig. 1b) preserves the peak of the upstream dune and occurs for larger values of ratio r than downstream-dominant events. Upon touching, the upstream peak starts to climb the stoss side of the downstream dune while retaining its avalanche face and recirculation zone. Sediment at the foot of the stoss slope is trapped by the upstream dune and the gradient in shear stress between the new trough and the downstream peak increases. This means the downstream dune shrinks as it erodes faster, whilst the upstream dune grows due to the extra mass it has assimilated during its migration as a superimposed bedform. Meanwhile, the trough height increases until a flat plateau forms as the now smaller downstream peak accelerates away (Fig. 2b and e). The critical characteristic is that the downstream peak shrinks faster than it migrates; it is unable to escape and it sinks into the plateau which then disappears as the upstream peak migrates forward.

The defining distinction between downstream- and upstream-dominant coalescence is the formation of a bounding surface within the internal structure of the merged dune, which extends into an open plateau located between two avalanche slopes. However, it is complicated to quantify the transition between the behaviours. Firstly, the simulated outcome is probabilistic, as will be demonstrated in the following section, so many simulations would be required to gain meaningful outcomes. Additionally, plateaus always form for all superimposed bedforms reaching the crest (Appendix A), even though they may be very small and short-lived. This means that, close to the transition, superimposed bedforms on the dune surface make it challenging to track the separate peaks and, consequently, difficult to distinguish the two types of behaviour. Following the coalescence, the upper part of the resultant dune is well mixed with sediment from both the upstream and downstream dunes, whilst the lower part contains only initially downstream material. The contact surface between the two sediment source transitions from the original slope of the impacted dune to a horizontal plane. However, the bounding surface is ultimately eradicated as the dune migrates forward and material is transported from the stoss to the lee side.

For still larger values of r, ejection (Fig. 1c) occurs. As with upstream-dominant coalescence, the upstream peak grows at the expense of the downstream peak and the trough becomes a plateau. However, this time the plateau is continually eroded and shrinks, providing sediment to the downstream peak, slowing its decrease in height. Eventually, the plateau in the trough is eroded, the upstream dune disconnects and the two peaks become distinct bedforms (Fig. 2c and f). After the interaction, the downstream dune still only contains initially downstream sediment, whilst the upstream dune has a structure similar to that seen at the end of upstream-dominant coalescence. Again, this intermediate structure is ultimately lost. It needs noting that, close to the transition to coalescence, the plateau between the two peaks can be completely eroded in multiple locations, resulting in the ejection of multiple small bedforms. Nevertheless, since the interaction mechanism is the same regardless of whether one or more downstream dunes are ejected, we classify all of these events the same.

The exchange of mass resulting from the three different types of interaction produces distinctive internal structures within the resultant dunes, as can be seen in Fig. 1. For downstream-dominant coalescence, the upstream dune breaks up as it climbs the stoss slope of the downstream dune, with the associated sediment transport up to and over the crest. Avalanching and deposition on the lee side creates a mixed layer, inclined at the angle of repose, within which the sediment from the initially upstream dune is confined. Assuming there is no difference between the sediment in the two initial dunes, this layer will just appear as one of many that form in the final dune as it propagates (Bagnold, 1941; Allen, 1970) and, thus, no evidence of the interaction will be preserved. Conversely, in the case of upstream-dominant coalescence, it is the downstream dune which is lost. In this case, whilst some of the downstream material becomes mixed into the upstream dune as it propagates over the downstream stoss slope, much of the material becomes contained within a basal horizontal layer. This is preserved only for the time it takes for the final dune to propagate a distance equal to its length. During this time, the material from the basal layer is eroded at the foot of the stoss slope, and transported up and over the dune, resulting in the final dune becoming well-mixed. Finally, in the ejection case, whilst the downstream dune only loses sediment, the upstream dune initially gains material in a similar fashion to the upstream-dominant coalescence scenario. A basal layer of initially downstream material forms in an entirely analogous way before ultimately becoming mixed into the upstream dune. Thus, for all three scenarios, although the collision creates distinct sedimentary structures, these are lost in the time it takes each dune to migrate a distance equal to its length. Therefore, these surfaces are unlikely to be preserved in the geological record except in zones of high deposition rates where dune interactions are synchronous with sediment accumulation.

3.2 Coalescence or ejection: randomness in the interaction outcome

All transitions in the cellular automaton are governed by stochastic processes as the model takes a random seed as an input (Narteau et al., 2009). We find that the transition from coalescence to ejection can depend on the seed; i.e. the phenomenon can be described probabilistically as opposed to deterministically. Figure 3a–c demonstrate this variability for a case close to the transition, showing three different outcomes (respectively coalescence and ejection with two and one ejected dunes) for simulations with $A_{\mathrm{u}}=\mathrm{2100}{l}_{\mathrm{0}}^{\mathrm{2}}$ and $A_{\mathrm{d}}=\mathrm{4000}{l}_{\mathrm{0}}^{\mathrm{2}}$ (r=0.525) but different seeds. Interactions that lead to an increase in the number of dunes only occur for values of r close to the transition.

Figure 3 (a–c) Interactions for $A_{\mathrm{u}}=\mathrm{2100}{l}_{\mathrm{0}}^{\mathrm{2}}$ and $A_{\mathrm{d}}=\mathrm{4000}{l}_{\mathrm{0}}^{\mathrm{2}}$ (r=0.525) but different random seeds. The colour shows the height h and time runs vertically. In (a) coalescence occurs whilst in (b) and (c) ejection is observed, with (b) showing multiple ejected bedforms. (d) Fraction of coalescence events as a function of $r=A_{\mathrm{u}}/A_{\mathrm{d}}$ for different values of A_d. Variation in the data is due to only 50 simulations being used to produce each point. The curve is the fitted function $f=\mathit{\left\{}\tanh \right[a(r-b)]+\mathrm{1}\mathit{\}}/\mathrm{2}$, with $a=\mathrm{16}\pm \mathrm{3}$ and $b=\mathrm{0.506}\pm \mathrm{0.006}$. (e) Contours of probability of ejection as a function of A_u and A_d.

To estimate the outcome probability for a given r, we performed 50 simulations with different seeds for 32 different pairs of dune areas, resulting in a total of 1600 simulations. Performing further simulations would be prohibitively computationally expensive. Then, for each set, we quantified the proportion of simulations resulting in either outcome. Figure 3d shows the fraction of ejection events f_eject as a function of r. We observe a narrow transition from no ejection events to all ejection events around $r=\mathrm{1}/\mathrm{2}$, with no apparent dependency on the absolute dune size. We note that this transition value is identical to that found by Katsuki et al. (2005) but greater than that found by Diniega et al. (2010). Unfortunately, we have no theoretical explanation for this result. Nonetheless, we fit the data to

$$\begin{array}{}\text{(5)}& f_{\mathrm{eject}}={\displaystyle \frac{\tanh \left[a\right(r-b\left)\right]+\mathrm{1}}{\mathrm{2}}},\end{array}$$

finding $a=\mathrm{14}\pm \mathrm{2}$ and $b=\mathrm{0.509}\pm \mathrm{0.005}$. Although any sigmoidal function could be used, the hyperbolic tangent provides the best fitting according to a sum of squared residuals criterion. We use Eq. (5) to define contours of probability for interaction outcome in the space {A_u,A_d} (Fig. 3e). It is important to note that the precise fitted values of a and b, and therefore the precise locations of the regime domains in Fig. 3e, are strictly only valid in 2D systems.

4 Comparison between dune interactions in the experiments and the simulations

We compare both the qualitative morphology of interacting bedforms in the experiments of Jarvis et al. (2022) and the simulations, as well as the quantitative regime transition for ejection and coalescence. The experiments involved the creation and evolution of discrete, quasi-2D dunes from a thin layer of glass beads in a narrow, counter-rotating, water-filled annular flume. Full details can be found in Jarvis et al. (2022). However, during the experiment, dunes were observed to undergo both coalescence and ejection interactions with their neighbours. Given the experiments enabled measurements of the heights of the dunes undergoing interactions, we are able to compare the observed behaviour with the simulation results presented here. However, there are caveats to consider. Firstly, we only simulated interactions between two discrete dunes, whilst, in the experiments, there is a train of interacting dunes. Secondly, the number of simulations (50) for each pair of dune areas is relatively small. However, since we have performed simulations for 28 different pairs of dune areas, we expect further simulations to just reduce the uncertainties on a and b. Thirdly, dunes in the experiments have finite width equal to that of the channel (9 cm) and some three-dimensional effects may impact the results (dunes are typically longer at the outer wall than the inner wall by ∼6 %). In particular, turbulent eddies will be shed by the dunes despite the narrow channel. Such 3D turbulent flow fields and their effect on dune collision dynamics cannot be reproduced in our 2D model. Finally, we observe experimental leeside slope angles $\mathit{\theta }\approx (\mathrm{18}\pm \mathrm{2})$^∘, whilst in our simulations we set θ=35^∘. Although we have performed some additional simulations to verify that ejection and coalescence occur for θ=18^∘ as they do for θ=35^∘, we currently neglect any influence θ may have on the coalescence–ejection transition. Given these caveats, we identified experimental interactions (see Table 1) for comparison according to specific criteria: (1) events took place between discrete dunes and (2) during the interaction there was no physical overlap with other neighbouring dunes. Nevertheless, it is necessary to state that the existence of a dune train will cause flow disturbances that could still affect collision outcomes. These events were then labelled as coalescence or ejection. Jarvis et al. (2022) reported that downstream-dominant coalescence was not observed in any experiment as the majority of bedforms were of comparable size due to the initial conditions. This is consistent with the prediction from the numerical results, which show that a large initial size difference is required to feature downstream-dominant coalescence. Jarvis et al. (2022) also did not report any observations of dune repulsion or collision suppression, as seen by Bacik et al. (2020) in their experiments on interactions between isolated dune pairs in periodic domains. This is likely because dune repulsion acts to push a two-dune system towards an antipodal configuration (Bacik et al., 2021). However, at any given time, the dunes in the experiments of Jarvis et al. (2022) were always relatively evenly spaced. Consequently, any dune-repulsion effects would have been very small and difficult to observe.

Table 1 Upstream and downstream heights, H_u and H_d respectively, along with the corresponding area ratio r and collision outcomes of the selected experiments from Jarvis et al. (2022). Also shown is the predicted probability of ejection f_eject as given by Eq. (5). Where f_eject>0.99, Eq. (5) gives a probability of ejection that, within the uncertainty, is equal to 1.

Morphologically, interacting bedforms in the experiments and simulations are similar. A distinctive feature of upstream-dominant coalescence seen in both datasets is the formation of a plateau towards the end of the interaction. Figure 4a and b compare the morphology of an experimental and a numerical example respectively. The plateau is seen in both cases and was ubiquitous in all upstream-dominant coalescence events. The same intermediate morphology was observed experimentally by Groh et al. (2009) and in active dune fields (Gao et al., 2015). This plateau could have implications for boundaries in the sedimentary record of previous dune-forming environments.

Figure 4 The morphology of intermediate bedforms in (a) experiments of Jarvis et al. (2022) and (b) simulations. Both show the distinctive plateau. (c) The location of experimental ejection and coalescence events on the numerically determined regime diagram. Most interactions result in their most likely outcome, with just four having a probability less than 0.5.

Figure 4c shows the location of the selected events in the regime diagram (as determined from simulations), which has been transformed to be in terms of dune heights as opposed to areas, since height is easier to measure experimentally. Given the assumption of dune self-similarity (Diniega et al., 2010), the height and area measurements will be equivalent. There is reasonable agreement between experiments and numerical predictions, with most interactions resulting in their most likely outcome; only 1 in 9 coalescence events and 3 in 10 ejection events had a probability less than 0.5. One ejection event, though, only had a probability of 0.002 according to our empirical model, suggesting that there may be other factors at play. This is unsurprising since in the experiments, the dunes interact with multiple neighbours, whilst the simulations consider only a pair of dunes. This can lead to changes in the sand flux entering and leaving an interacting pair of dunes which may affect the outcome. Appendix C presents a visual means of comparing the experimental observations with the empirical probability distribution determined by the simulations.

5 Concluding remarks

Dune collisions are an important part of dune field evolution and pattern coarsening and, in 2D, can result in two possible outcomes: coalescence and ejection (Coleman and Melville, 1994; Endo et al., 2004; Durán et al., 2005; Katsuki et al., 2005; Diniega et al., 2010; Bo and Zheng, 2013; Gao et al., 2015). Here we use a cellular automaton model (Narteau et al., 2009) to simulate collisions between discrete two-dimensional dunes and show that the collision outcome can be modelled as probabilistically depending on the dune size ratio. The observation that this is a probabilistic rather than a deterministic dependence is an interesting result which, to the authors' knowledge, has not been seen before. Whilst this requires experimental validation, it does not seem unreasonable given that sediment transport is, by nature, stochastic (Pähtz et al., 2020). Furthermore, we determine an empirical relationship for the probability of ejection as a function of the dune area ratio r (Eq. 5) with the transition centred on $r\approx \mathrm{1}/\mathrm{2}$. This is in contrast to the work of Diniega et al. (2010), who previously used a continuum model to predict the transition at $r=\mathrm{1}/\mathrm{3}$. However, it is in agreement with the prediction of $r=\mathrm{1}/\mathrm{2}$ from Katsuki et al. (2005). Nonetheless, the models of both Katsuki et al. (2005) and Diniega et al. (2010) are deterministic, so they cannot reproduce our probabilistic results. We also found no reason for the difference in the predicted transitional value between that found by Katsuki et al. (2005) and this study and that found by Diniega et al. (2010). This requires further investigation.

Our numerical simulations also show that coalescence occurs in two varieties, upstream or downstream dominant, which, together with ejection, result in three elementary types of dune interaction in 2D. Despite their probabilistic nature, these elementary types of interaction arise in specific dune size ranges and can be recognised in 3D even in the presence of turbulence and multi-directional flow regimes. Thus, they potentially provide an efficient means of decomposing coarsening phases and associated timescales (Lü et al., 2021).

Furthermore, we compare the numerical observations with ideal experimental interaction events from the study of Jarvis et al. (2022). Morphologically, the interactions appeared very similar, as evidenced by the existence of a flat plateau that appears during upstream-dominant coalescence. Additionally, the interaction outcomes agree well with the numerically determined regime diagram. Although these comparisons are favourable, restricting our simulations to a 2D domain means we are unable to reproduce the full range of dune interaction phenomena that can be observed in 3D systems (Endo et al., 2004; Durán et al., 2005; Katsuki et al., 2005; Assis and Franklin, 2020, 2021). In particular, the model cannot capture the role of transverse recirculation and the channelisation of an oblique flow between dunes when they are laterally offset. A fully 3D fluid velocity field, and possibly a variable lee slope angle, would be required for the cellular automaton model to reproduce these behaviours and the more general phenomenology of dune interactions observed in laboratory experiments, including the collision-suppression and dune-repulsion phenomena observed by Bacik et al. (2020). Our results highlight a need for experiments on interactions between two discrete dunes (a) to verify if the outcomes are indeed probabilistic and (b) to quantify the transition between the different regimes. Such experiments would also enable further quantitative comparison with numerical simulations, including constraints on mass exchange and velocity evolution with time.

Appendix A: Single-dune simulations

To determine the required domain size and initial dune separation, simulations of a single dune were performed. These simulations were initiated in a periodic domain of variable length L=600–19 200l₀ with a triangular pile of sand of variable area (400–$\mathrm{4000})l_{\mathrm{0}}^{\mathrm{2}}$ and a fixed slope angle of 35^∘ in the centre. During each simulation, the mound morphed into a steady shape with a curved stoss slope and straight slip face at the angle of repose. Regardless of the size of the initial mound, this equilibrium shape was attained within 5000 time steps (Fig. A1a) and a migration distance of 350l₀, consistent with models for 3D barchans (Zhang et al., 2014). Hence, an initial separation > 1000l₀ in the simulation is sufficient. Once equilibrated, the position of the dune crest with time was used to calculate the migration velocity c. Our results for the dependency of c on A agree with Zhang et al. (2010), although for L≲1200l₀, c depends on L due to long-range interactions between the dune and its repeated images in the periodic domain (Fig. A2). We are therefore justified in our choice of L=2000l₀. We also note that the dune aspect ratio $\mathcal{H}/\mathcal{L}$ is independent of A over the full range of dune sizes we used, thus demonstrating scale invariance.

Superimposed bedforms were observed on the stoss slopes of these dunes (Fig. A1b). These disturbances to the dune surface appeared at the upstream foot of the dune and propagated to the crest, growing as they climbed. Once a superimposed peak approached the dune top, it would sometimes be higher than the crest. In this case, the crest would shrink in height and a short plateau would form at the top of the slip face. The new peak migrated across this plateau whilst migration of the dune itself paused until the peak reached the slip slope and avalanching recommenced.

Figure A1 (a) Time evolution of the dune aspect ratio $\mathcal{H}/\mathcal{L}$ for different dune areas A. The bedform is initiated as an isosceles triangle with slope angle 35^∘, and an equilibrium dune shape is seen to form within 5000 t₀ for all areas considered. (b) Time-averaged RMS amplitude of superimposed bedforms ${\stackrel{\mathrm{‾}}{\mathcal{H}}}_{\sup }$ as a function of the fractional distance along the stoss slope f for equilibrated dunes. It can be seen that the superimposed bedforms grew as they climbed the stoss slope.

Figure A2 Migration velocity c of an individual dune as a function of dune area A and domain length L. We see that the dune velocity decreases with dune area. For L≲2000l₀, c also increases with L, particularly at large A. This is due to long-range interactions between the dune and its repeated images in the periodic domain.

Appendix B: Parameter values

Table B1 shows values of parameters used in the simulation.

Table B1 Values of parameters used in the cellular automaton in model and physical units. Conversion to physical units was performed using the scales determined in Zhang et al. (2014).

Appendix C: Separation plot

An alternative means of comparing the experimental data with the probability model (Eq. 5) determined from the numerical simulations is presented in Fig. C1. This shows a separation plot, a visual tool that can be used to test the reliability of models with binary outcomes (Greenhill et al., 2011). To create this plot, the 19 experimental data points shown in Fig. 4c are ordered from the smallest $H_{\mathrm{u}}/H_{\mathrm{d}}$ on the left to the largest on the right. Events which result in ejection are represented as red columns whilst those represented as coalescence remain white. If the model given by Eq. (5) were to be entirely uncorrelated with the experimental results, we would expect the red columns to be uniformly dispersed along the axis. However, there is an increasing concentration of red columns towards the right hand side of the plot, showing that ejection events indeed occur more frequently when Eq. (5) predicts them to be more likely. Overlain on the plot is a line representing the calculated probability of ejection given by Eq. (5) for each observation. This allows us to see that the only true outlier is the ejection event that occurs for H_u=26.0 mm, H_d=46.3 mm and has a probability of 0.002. As suggested by Fig. 4c, it seems that Eq. (5) does a reasonably good, but not perfect, job of describing the data.

Figure C1 A separation plot, generated by ordering the 19 experimental data points from Fig. 4c from the smallest $H_{\mathrm{u}}/H_{\mathrm{d}}$ on the left to the largest on the right. Instances resulting in ejection are represented as red columns. Overlain is a line chart showing the probability of each interaction resulting in ejection, as given by Eq. (5).