Barchan swarm dynamics from a Two-Flank Agent-Based Model

In this work, we simulate barchan swarms using the Two-Flank Agent-Based Model and investigate how changes to model parameters and environmental drivers lead to different swarm dynamics. In particular, we explore how the parameter q_shift, which controls the rate of equilibration between the flanks of an asymmetric dune and thereby the stability of an asymmetric morphology, influences the frequencies of different collision types in the swarm and can be adjusted to produce swarms with dune sizes that are longitudinally homogeneous. Such size stability has been observed in real-world swarms but has not been obtained by previous agent-based models. We also find that, for certain densities of newly added barchans, the dune number density remains constant with downwind distance, something which has also been reported in nature but not in previous agent-based models. We also investigated how boundary condition changes propagate through swarms, something which will become increasingly important in the face of human infrastructure expansion and anthropogenic climate change. Finally, we are able to investigate how asymmetry and spatially patterning in the swarms are influenced by bimodal wind regimes. Since these simulations produce many more realistic phenomena than previous models, this work represents a significant step forward in the our understanding of the interplay between environmental conditions and dune interactions in shaping the dynamics of barchan swarms.

1 Introduction

Crescent-shaped barchan dunes are one of the most striking aeolian bedforms and are found on Mars (Bourke, 2010; Rubanenko et al., 2022; Bourke and Goudie, 2009; Robson et al., 2022) and many different locations on Earth (Goudie, 2020; Elbelrhiti et al., 2008; Hersen et al., 2004; Robson et al., 2022). Compared to other types of bedform, barchans migrate quickly, often up to hundreds of metres per year (Sparavigna, 2013; Corbett, 2018; Barnes, 2001). The rate at which the dunes move is controlled primarily by their size (Long and Sharp, 1964; Bagnold, 1941), which means that collisions can occur between the bedforms.

Individual barchans are typically tens of metres wide (Elbelrhiti et al., 2008; Durán et al., 2009) and evolve over timescales on the order of years. The evolution of small numbers of such bedforms can, therefore, be studied using numerical models of the microscopic behaviours of the sediment and fluid flow (Durán et al., 2011; Xiao et al., 2023). It is also possible to study proxies for aeolian barchans in a laboratory setting using a fluid, e.g. water, which is much denser than air (Hersen, 2005; Hersen et al., 2002; Assis and Franklin, 2020, 2021). In nature, however, we do not typically observe barchans forming in small numbers but rather vast populations of interacting dunes (Hersen et al., 2004). These huge systems, known as swarms (Hesse, 2009; Robson et al., 2022), can span tens or hundreds of square kilometres and contain hundreds, thousands, or even tens of thousands of bedforms (King, 1918; Durán et al., 2009; Elbelrhiti et al., 2008).

Swarms display many emergent phenomena such as size selection (Durán et al., 2009; Robson et al., 2022; Génois et al., 2013a), longitudinal homogeneity (Durán et al., 2011; Elbelrhiti et al., 2008), and horn-to-toe alignment (Elbelrhiti et al., 2008), all of which must result from interactions between the bedforms. Yet the exact relationships between barchan interactions and the observed properties of swarms are not well understood. This is an important issue, since it has been shown that anthropogenic climate change will affect many active aeolian systems (Engel et al., 2018; Gunn et al., 2022; Baas and Delobel, 2022), yet it will be difficult to understand the effect these changes will have on swarms without a sufficiently realistic model. This goes beyond merely academic curiosity, since the migration of barchans can pose a direct threat to human infrastructure, heritage sites, and settlements (Boulghobra, 2016; Hermas et al., 2012; Effat et al., 2011; Ding et al., 2020; Finkel, 1959; Lorenz et al., 2013); thus, understanding how dunes in swarms will evolve under contemporary climate change is vital for future planning.

Not only are the length scales of swarms orders of magnitude larger than those of the dunes themselves, but the timescales over which these systems evolve are also much longer than for individual barchans (Hesse, 2009; Xiao et al., 2023). Such length and timescales prohibit use of microscopic simulations due to computational costs (Durán et al., 2011) and laboratory experiments which cannot reach a sufficient size. Instead, one must turn to an alternative approach, agent-based models (ABMs), which treat the dunes themselves as the fundamental objects and distil the complexities of the microscopic motion of sediment into rules for the size evolution and migration of barchans and interactions between the bedforms (Lima et al., 2002).

Several ABMs have been implemented to study barchan swarms (Lima et al., 2002; Durán et al., 2011; Worman et al., 2013; Génois et al., 2013b), each of which included a different subset of known interaction types. Alongside these, several one-dimensional ABMs have also been applied to aeolian systems (Parteli and Herrmann, 2003; Lee et al., 2005; Diniega et al., 2010); however, since barchan swarms are two-dimensional, the results can be very different.

The earliest barchan ABM (Lima et al., 2002) demonstrated that sand flux exchange could lead to spatial structuring, while, since only merging interactions were simulated, the overall size of dunes in the system increased with downwind distance. Allowing exchange collisions to dominate was shown to produce more realistic global size distributions (Durán et al., 2011) but did not eliminate the problem that dune size increased with downwind distance (Durán et al., 2011). When fragmentation collisions, instead, dominate the collision dynamics, clusters form, which exhibit a runaway increase in density which is mitigated only by the destruction of small dunes and an overall loss of sand that has to be compensated for through the spontaneous introduction of new dunes throughout the swarm (Génois et al., 2013b). Furthermore, when fragmentation and sand loss completely dominate the behaviour of a swarm; the largest dunes in the system are the newly injected dunes (Génois et al., 2013b), which contradicts the intuition that newly formed dunes should be the smallest in a swarm (Elbelrhiti, 2012). Finally, when collisions result only in merging, while dunes above a certain size are allowed to spontaneously fragment (calve), the distributions are found to display a large peak at the size at which calving begins to occur and a smaller peak at the size at which newly calved dunes form with very few barchans at other sizes (Worman et al., 2013). The swarms generated with all these model versions contrast with real-world swarms which are found to have typical distributions which are similar to a log-normal distribution (Durán et al., 2009). Thus, none of the existing barchan swarm ABMs has successfully reproduced realistic swarms. Furthermore, each of the previous works has considered the idealised case of symmetric barchans subject to an idealised unidirectional wind.

We have developed a new barchan ABM, the Two-Flank Agent-Based Model (TFABM), which has been shown to not only reproduce all of the observed collision behaviour of barchans but which can also be used to simulate the formation of asymmetric bedforms under the influence of variable wind regimes (Robson and Baas, 2023), thus successfully reproducing the angular dependence of asymmetry growth due to bimodal winds as suggested in previous studies (Parteli et al., 2014; Lv et al., 2016). In this paper, we report the characteristics and behaviours of barchan swarms under unidirectional, as well as bimodal winds, demonstrating that the new model is able to more realistically reproduce certain aspects of barchan swarms than previous ABMs and thus furthers our understanding of the connection between barchan interactions and swarm properties.

2 Model description

The TFABM and its capabilities to replicate all known barchan collision dynamics were introduced in Robson and Baas (2023), and extensive details of the model structure can be found there. However, in this section, we review the key features and provide definitions of the model parameters which are explored in this study.

2.1 Dune morphology

The major difference between the TFABM and previous barchan ABMs is the structure of the dunes (agents) themselves. In the TFABM, each dune is represented by its two flanks, which are able to change size semi-independently of one another. The position of the dunes is recorded as the coordinates of the upwind toe. The central axis of the dune extends in the direction of the dominant wind, which we define to be the positive x direction. The size of the dune is represented by the widths W_p,s of the port and starboard flanks. The widest point of the dune occurs at a distance ${\mathit{\lambda }}_{\mathrm{1}}(W_{\mathrm{p}}+W_{\mathrm{s}})/\mathrm{2}$ downwind of the toe, where we use λ₁=1, based on comparison with an existing dune dataset (Sherman et al., 2021). The total lengths of the flanks, from the toe to the tip of the horn, are assumed to be proportional to their widths $L_{\mathrm{p},\mathrm{s}}={\mathit{\lambda }}_{\mathrm{2}}{W}_{\mathrm{p},\mathrm{s}}$, where we use λ₂=1.8, based on real-world dunes (Sherman et al., 2021). Note that because these relationships between length and width are fixed, the TFABM cannot produce asymmetry where a single limb becomes highly extended longitudinally without the bulk of the barchan changing significantly. However, such asymmetry was reported to primarily occur due to bimodal winds with obtuse angular separation (Parteli et al., 2014), and in the later sections, we focus only on bimodality with acute angular separation for which the model performs well.

The width of each horn is a linear function of the flank width $H_{\mathrm{p},\mathrm{s}}=\mathit{\alpha }{W}_{\mathrm{p},\mathrm{s}}+\mathrm{\Delta }/\mathrm{2}$, with α=0.05 and Δ=4.6 m (Hersen et al., 2004; Worman et al., 2013). Finally, the volume of a flank is given by

$$\begin{array}{}\text{(1)}& V_{\mathrm{p},\mathrm{s}}={\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{2}}{\mathit{\lambda }}_{\mathrm{3}}{W}_{\mathrm{p},\mathrm{s}}^{\mathrm{3}}}{\mathrm{6}}},\end{array}$$

where λ₃ is interpreted as the proportionality coefficient of flank height/width, we set this to one-third, such that a symmetric dune has volume $V_{\mathrm{tot}}=\mathrm{1}/\mathrm{40}$ $W_{\mathrm{tot}}^{\mathrm{3}}$, which is in agreement with Elbelrhiti et al. (2008). For further details of the morphology of dunes in the TFABM, see Robson and Baas (2023).

2.2 Sand flux

Unlike previous barchan ABMs, the more realistic morphology of dunes in the TFABM allows for variation in the direction of wind. Incoming sand flux is absorbed by each flank across the windward projection of its width. For oblique winds, this means that the windward flank will absorb proportionally more of the flux, starving the leeward flank of material.

As well as absorbing material across their wind-facing width, dune flanks also lose material across the width of their horns. In the initial study in which the TFABM was introduced, flux was emitted from the horns at a saturated rate, regardless of the influx the horn received, in line with findings from continuum simulations (Hersen et al., 2004). Such a flux rule has also previously been applied in an earlier ABM (Worman et al., 2013), and we will refer to it as “unscaled” flux. However, others have suggested that the outflux from the horns should scale with the influx such that a dune which receives no incoming sediment will also emit comparatively less (Durán et al., 2010). Such a rule, which we will call “scaled” flux, has also previously been applied in an ABM (Durán et al., 2011).

Finally, in asymmetric dunes, there is an effective transfer of material from the larger to the smaller flank. If this were not the case, then any asymmetry in a barchan would persist, which contradicts with the view of symmetric barchan morphology as an attractor (Werner, 1995; Groh et al., 2008), as well as with the results of continuum simulations (Parteli et al., 2014) and water tank experiments (Assis and Franklin, 2020). In the TFABM, the rate of this lateral transfer of material is governed by the phenomenological model parameter q_shift and occurs at a rate proportional to the difference in the streamwise projection of the flank lengths (note that this is not the same as the full length during periods of oblique flow). Combining the absorption, emission, and lateral equilibration flux, the evolution of the flank volume is governed by

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\mathrm{d}{V}_{\mathrm{p},\mathrm{s}}}{\mathrm{d}t}}=q_{\mathrm{in}}{\stackrel{\mathrm{̃}}{W}}_{\mathrm{p},\mathrm{s}}-q_{\mathrm{out}}{H}_{\mathrm{p},\mathrm{s}}+q_{\mathrm{shift}}\sin |\mathrm{90}-\mathit{\theta }|\left(L_{\mathrm{s},\mathrm{p}}-L_{\mathrm{p},\mathrm{s}}\right),\end{array}$$

where q_in is the received influx, including both ambient flux and that streaming off of the horns of upwind dunes which overlap the flank; θ is the wind direction; and ${\stackrel{\mathrm{̃}}{W}}_{\mathrm{l},\mathrm{r}}$ is the flank width perpendicular to the wind. The outflux term q_out is simply q_sat, the saturated sand flux, set to 79 m² yr⁻¹, which is the average value of those reported for Tarfaya in Elbelrhiti et al. (2008) in the case of unscaled flux but for scaled flux is a linear function of the influx as follows:

$$\begin{array}{}\text{(3)}& q_{\mathrm{out}}=\left(a{q}_{\mathrm{in}}+b{q}_{\mathrm{sat}}\right){\displaystyle \frac{W_{\mathrm{p},\mathrm{s}}}{H_{\mathrm{p},\mathrm{s}}}},\end{array}$$

where a and b are scaling parameters set to 0.45 and 1 as in Durán et al. (2011). Note that this change will significantly influence the rate at which dunes of different sizes grow. In particular, W_tot≲28 m outflux will be higher under scaled flux, while for larger dunes the outflux is higher under the unscaled flux.

2.3 Calving

The sizes of the flanks of a dune are coupled only through the third term in Eq. (2) above. However, when dunes become very asymmetric, this term may become larger than the other terms in Eq. 2, meaning that spanwise fluxes exceed streamwise fluxes. Since this is no longer the typical behaviour of a barchan, we take this to define a maximum asymmetry threshold for dunes in our model. This threshold occurs when the asymmetry ratio, defined as $\mathit{\gamma }=\max (W_{\mathrm{p}},W_{\mathrm{s}})/\min (W_{\mathrm{p}},W_{\mathrm{s}})$, exceeds

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{\mathit{\gamma }}_{\mathrm{c},\mathrm{shift}}\left(W_{\mathrm{p}},W_{\mathrm{s}}\right)& =\mathrm{1}+{\displaystyle \frac{\mathit{\alpha }{q}_{\mathrm{sat}}}{{\mathit{\lambda }}_{\mathrm{2}}{q}_{\mathrm{shift}}-\mathit{\alpha }{q}_{\mathrm{sat}}}},\\ & +{\displaystyle \frac{\mathrm{\Delta }{q}_{\mathrm{sat}}}{\mathrm{2}\left({\mathit{\lambda }}_{\mathrm{2}}{q}_{\mathrm{shift}}-\mathit{\alpha }{q}_{\mathrm{sat}}\right)\min \left(W_{\mathrm{p}},W_{\mathrm{s}}\right)}}.\end{array}\end{array}$$

For highly asymmetric dunes, the typical barchan geometry shown in Robson and Baas (2023) can break down if the asymmetry ratio exceeds the threshold

$$\begin{array}{}\text{(5)}& {\mathit{\gamma }}_{\mathrm{c},\mathit{\lambda }}={\displaystyle \frac{\mathrm{2}{\mathit{\lambda }}_{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{1}}}}-\mathrm{1}.\end{array}$$

Thus, combining these two terms, the maximum asymmetry of a dune with flank widths W_p and W_s is

$$\begin{array}{}\text{(6)}& {\mathit{\gamma }}_{\mathrm{c}}\left(W_{\mathrm{p}},W_{\mathrm{s}}\right)=\min \left({\mathit{\gamma }}_{\mathrm{c},\mathit{\lambda }},{\mathit{\gamma }}_{\mathrm{c},\mathrm{shift}}\left(W_{\mathrm{p}},W_{\mathrm{s}}\right)\right).\end{array}$$

When a dune exceeds this maximum asymmetry threshold, its two flanks spontaneously break apart and are reformed into symmetric barchans, a process we refer to as calving.

2.4 Collisions

One of the reasons that barchans have been studied so extensively is that they migrate rapidly at a rate governed by the sand flux conditions and the size of the bedform. We use the commonly accepted expression for the migration rate as

$$\begin{array}{}\text{(7)}& v_{\mathrm{mig}}={\displaystyle \frac{c{q}_{\mathrm{sat}}}{W_{\mathrm{p}}+W_{\mathrm{s}}+W_{\mathrm{0}}}},\end{array}$$

where c=45 (Elbelrhiti et al., 2008) relates to the speeding up of wind passing over the crest of a dune, and W₀=16.6 m (Elbelrhiti et al., 2008) ensures that the migration rate of small dunes remains finite.

From this expression, it follows that larger dunes will migrate more slowly than smaller ones, which ensures that collisions are frequent occurrences in barchan swarms. The TFABM uses a simple rule for collisions, which is, nonetheless, able to reproduce a wide range of known collision behaviours (Robson and Baas, 2023). Each dune has three associated centres of mass (CoMs), namely one for each flank and one for the combined dune. The collision algorithm is triggered whenever one of these centres of mass intersects with the footprint of another dune. When that condition is met, the following algorithm is used to determine the output of a collision:

1. The intersecting flanks merge together.

2. The volume of merged flanks is used to calculate an effective flank width W_merged.

3. A non-intersecting flank with the width W_i will also join the merged flanks if $W_{\mathrm{merged}}/W_i\le {\mathit{\gamma }}_{\mathrm{c}}(W_{\mathrm{merged}},W_i)$.

4. If $W_{\mathrm{merged}}/W_i>{\mathit{\gamma }}_{\mathrm{c}}(W_{\mathrm{merged}},W_i)$, then the non-intersecting flank forms into a separate symmetric barchan in the same manner as during calving.

5. Finally, the flanks which merged together form a barchan with a preserved centre of mass and an asymmetry ratio determined by the relative mass that was initially to the left and right of the centre of mass.

The algorithm itself can produce between one and three output dunes, although calving may occur on the following time step, such that a collision may effectively produce up to four outputs. We refer to collisions which decrease the number of dunes (i.e. have a single output) as merging or aggregation, collisions which conserve the dune number as exchange, and collisions which increase the dune number as fragmentation. For more information on the phase space of the TFABM collision rule, see Robson and Baas (2023). In particular, the supplementary material of Robson and Baas (2023) compares the ability of the collision rules applied in the different barchan ABMs to reproduce collision phase spaces from cellular automata (Katsuki et al., 2011) and subaqueous experiments (Assis and Franklin, 2020), revealing that the TFABM is able to reproduce the greatest number of collision outcomes.

2.5 Boundary conditions and dune injection

The simulation space has a downwind length of 10 km and a central strip with a width of 3 km or, in some cases, 5 km set to ensure a sufficient number of dunes were present. On either side of the central strip there were empty strips of the same width to allow for an oblique movement of the dunes. At the upwind boundary over the entire width (9 or 15 km) of the simulation space, we supply a free flux q₀=0.25 q_sat (Elbelrhiti et al., 2008). This flux is projected downwind until it is intercepted by a dune, downwind of which the free flux is zero.

In this work, we do not use periodic boundaries; thus, once dunes exit the simulation space, they are lost. To account for this, new dunes enter the simulation space at the upwind boundary. The rate of injection of new dunes is governed by assuming that just upwind of the domain, there is a number density ρ₀ of symmetric barchans each with a width $W_{\mathrm{eq}}=\mathrm{\Delta }{q}_{\mathrm{sat}}/(q_{\mathrm{0}}-\mathit{\alpha }{q}_{\mathrm{sat}})$, which corresponds to the known unstable equilibrium size under the unscaled outflux (Hersen et al., 2004; Worman et al., 2013). Dunes are injected only into the central 3 or 5 km wide strip of the simulation space to ensure that the dunes do not reach the lateral borders. The length of a time step was set to ${\mathit{\delta }}_t=\mathrm{1}/\mathrm{8}$ years ≈45 d. As we will see, the smallest dunes in the system are typically the newly created bedforms with width W_eq. With our choice of q_sat, c, and W₀ (see Eq. 7), the turnover time (time taken for a dune to migrate its own length) of barchans with width W_eq is t_turnover(W_eq)=84 d, meaning that our choice of time step is still much shorter than the shortest timescale relevant to the dunes in our simulations. The number of dunes entering the system each time step is then determined by

$$\begin{array}{}\text{(8)}& N_{\mathrm{enter}}={\displaystyle \frac{W_{\mathrm{simulation}}}{\mathrm{3}}}{\mathit{\rho }}_{\mathrm{0}}{v}_{\mathrm{mig}}\left(W_{\mathrm{eq}}\right){\mathit{\delta }}_t.\end{array}$$

3 Results

3.1 Unimodal wind

For the first series of simulations, we use a Gaussian distribution for the direction of the wind centred around 0°. We experimented with different values of the standard deviation but found that, provided the standard deviation remains relatively small (we tested up to 5°), the exact value is not important. Therefore, all of the simulations shown in this section were performed using a standard deviation of 3°. Since the wind angle is normally distributed, this standard deviation ensures that 99.7 % of the time, the angle of the wind is within ±9° of the primary mode. Under this wind regime we investigated the effect of changes in q_shift and ρ₀, as well as the differences between scaled and unscaled outflux. We also performed simulations in which the properties were allowed to stabilise before the injection density was changed in order to explore how the TFABM might be applied to investigate boundary condition change (Marvin et al., 2023).

3.1.1 Varying q_shift

As shown in Robson and Baas (2023), the rate of equilibration between the flanks, represented by the flank-balancing flux, q_shift, has a major impact on the output of collisions. This is because a higher value of q_shift means a lower maximum asymmetry ratio (see Eq. 4), thus increasing q_shift leads to a larger portion of the collision phase space, resulting in fragmentation collisions (Robson and Baas, 2023). However, the properties of the entire phase space of collisions may not be relevant in a real-world swarm, since it has been shown that barchans in swarms exhibit spatial structuring, such as preferential alignment (Elbelrhiti et al., 2008; Lima et al., 2002), which may result in collisions occurring more frequently in a particular part of the phase space, i.e. only a small range of lateral offsets. Therefore, to investigate the effects of varying q_shift, we fixed ρ₀=37 km⁻² and performed simulations with q_shift=0, 0.05 q_sat, 0.1 q_sat, and 0.15 q_sat. Each of these simulations was performed over a timescale long enough for the entire 10 km length of the simulation space to be filled and for the properties of the dunes to have stabilised; this meant simulations of >350 years in all cases. In Figs. 1 and 2, we show the final states of these simulations for unscaled and scale outfluxes respectively alongside plots showing that the population size and mean width had stabilised during the simulations. Animations showing the entire evolution of the swarms can be found in Robson (2024a).

Figure 1 The final states of the simulated swarms and plots showing that the mean dune width and population size had stabilised for the unscaled outflux, ρ₀=37 km⁻². $q_{\mathrm{shift}}/q_{\mathrm{sat}}$ = (a) 0.00, (b) 0.05, (c) 0.10, and (d) 0.15.

Figure 2 The final states of the simulated swarms and plots showing that the mean dune width and population size had stabilised for the scaled outflux, ρ₀=37 km⁻². $q_{\mathrm{shift}}/q_{\mathrm{sat}}$ = (a) 0.00, (b) 0.05, (c) 0.10, and (d) 0.15.

Regardless of the choice of outflux, lower values of q_shift lead to a smaller total number of dunes in the system, and therefore a lower overall density, while the average size of dunes in simulations with lower q_shift is greater than when q_shift is high. To explain this, we show in Fig. 3 how changing q_shift affects the frequency of the different collision types. We observe that, as q_shift increases, the frequency of merging collisions decreases, while the frequency of fragmentation increases sharply from the least frequent at q_shift=0 to the most frequent when q_shift=0.15 in the case of unscaled flux or, similar to exchange collisions, in the case of scaled flux.

Figure 3 Relative frequencies of different types of collisions: merging (circle), exchange (triangles), and fragmentation (squares) observed in simulations with ρ₀=37 km⁻². $q_{\mathrm{shift}}/q_{\mathrm{sat}}=\mathrm{0.0}$, 0.05, 0.1, and 0.15, with scaled (solid markers) and unscaled (open markers) outfluxes.

It has been found that the sizes of dunes in real-world swarms remain stable in the longitudinal dimension (Elbelrhiti et al., 2008; Durán et al., 2011); however, this is not something that previous agent-based models have been able to reproduce (Lima et al., 2002; Durán et al., 2011; Worman et al., 2013). In Fig. 4, we show that, in contrast, when using unscaled outflux, we are able to produce stabilised sizes of dunes when q_shift=0.15 q_sat, with the swarms exhibiting a trend whereby the sizes of dunes increase until around 2–3 km, after which the size distribution remains approximately steady. For the unscaled outflux and $q_{\mathrm{shift}}/q_{\mathrm{sat}}\le \mathrm{0.1}$, the stabilisation of sizes did not occur within the simulation space, though it is possible that it may have done further downwind. Such an increase in size has been observed in narrow real-world swarms (or monocorridors) but not large swarms, as simulated here (Elbelrhiti et al., 2008; Durán et al., 2011).

Figure 4 The mean width and dune density in 500 m cross sections averaged from measurements at the end of each year once the swarm properties stabilised. $q_{\mathrm{shift}}/q_{\mathrm{sat}}=\mathrm{0}$, 0.05, 0.1, and 0.15 in panels (a)–(d) respectively and unscaled dune outflux. The grey areas represent 1 standard deviation.

Figure 5 shows that when dune outflux scaled with influx, low values of q_shift (0 or 0.05 q_sat) lead to swarms in which dune size increases with downwind distance; for q_shift=0.1 q_sat, the size stabilises at some distance, beyond which it remains constant, and for q_shift=0.15, the dune size initially increases before reaching an apex and subsequently decreasing with downwind distance. In other words, it is possible to produce swarms with longitudinally homogeneous size distributions for both types of outflux, but the value of q_shift required to do so is higher for the case where the flux streaming off of the dunes is not scaled by the flux they receive.

Figure 5 The mean width and dune density in 500 m cross sections averaged from measurements at the end of each year once the swarm properties stabilised. $q_{\mathrm{shift}}/q_{\mathrm{sat}}=\mathrm{0}$, 0.05, 0.1, and 0.15 in panels (a)–(d) respectively and scaled dune outflux. The grey areas represent 1 standard deviation.

Whereas, in some instances, we observe longitudinally homogeneous dune sizes, the simulated swarms are not homogeneous in the dune number density. In all of the simulated swarms with unscaled outflux, the overall dune density shows a decreasing trend with downwind distance. On the other hand, for the scaled outflux, the dune density initially increases to a maximum and then begins to decrease in all cases. However, for q_shift=0.15, the density of the profile is approximately flat between around 2.5–6 km downwind; i.e. there is an area of around 4 km over which the density profile appears homogeneous.

We have already shown how varying q_shift affects the relative frequency of the different types of collisions; however, q_shift also controls asymmetry growth and calving in the TFABM. In Fig. 6, we show the different asymmetry distributions (the ratio of the widths of the flanks) from the scaled flux simulations with different q_shift alongside the asymmetry distribution observed across three areas of a swarm, containing a total of around 3000 dunes, in Tarfaya (Robson et al., 2022), which has a close-to-unimodal wind regime (Robson and Baas, 2024b). The data on dune morphology used to calculate the asymmetry distributions are freely available from Robson and Baas (2024a). The comparison reveals that the Tarfaya asymmetry distribution lies somewhere between the simulated distributions for q_shift=0 and 0.05 q_sat, with much lower values than those which produce longitudinally homogeneous size distributions. However, the asymmetry distributions of other real-world swarms have very different standard deviations, as shown in Fig. 7. In particular, the width of the asymmetry distribution for a swarm in Mauritania is best fit by the simulated swarm with q_shift=0.1 q_sat for which the size distribution is longitudinally homogeneous.

Figure 6 Histograms of dune asymmetries defined at the ratio of port and starboard flank widths are shown in panels (a)–(d) for respectively $q_{\mathrm{shift}}/q_{\mathrm{sat}}=\mathrm{0}$, 0.05, 0.1, and 0.15, with ρ₀=37 km⁻² and scaled flux. The dotted black lines show the distribution from a real-world swarm in Tarfaya.

Figure 7 The standard deviation of the ratio of flank widths for ($W_{\mathrm{p}}/W_{\mathrm{s}}$) for the simulated swarms and real-world swarms in Tarfaya, Mauritania, and two swarms in the northern circumpolar region of Mars. For the natural swarms, the ordering of the flank widths in the ratio was chosen to minimise the mean of the distributions.

3.1.2 Varying ρ₀

The previous section shows that varying q_shift leads to significant differences in the simulated swarms through changing the collision, asymmetry, and calving behaviour of the dunes. However, the rate of injection of new dunes at the upwind domain boundary was kept at ρ₀=37 km⁻² in those simulations. Real-world swarms, however, are observed to have a wide range of densities (Elbelrhiti et al., 2008). Figure 8 shows simulated swarms and property convergence for q_shift=0.15 q_sat and used ρ₀=12, 24, and 37 km⁻² (these correspond to 5%, 10%, and 15% of the density of Tarfaya 1 in Robson and Baas, 2024a; Robson et al., 2022). We focus on the results from the scaled outflux simulations; however, the swarms displayed similar properties in the case of unscaled outflux. As seen in Fig. 8, as the injection density increases, so too does the overall number of dunes in the swarm, while, as in the case of increasing q_shift, the mean sizes of dunes decreases. While these densities are lower than some swarms (Robson et al., 2022), there are many examples of swarms with similar densities to these (Elbelrhiti et al., 2008). It is possible to simulate swarms with a higher density using the TFABM, but the computational time increases with dune number, meaning that high-density swarms need more time to simulated.

Figure 8 The final states of simulations and stabilisation of the dune number and mean size with time for swarms with scaled outflux. q_shift=0.1 q_sat and ρ₀=12, 24, 37 km⁻² in panels (a)–(c) respectively. Note that in panels (a) and (b), the width of the swarms is 5 km, and in panel (c) the width is 3 km.

The longitudinal variation in the dune size and the number density are shown in Fig. 9. Since q_shift remained constant across the simulations, the relative frequency of the different collision types did not vary, leading to all three simulations reaching a longitudinally homogeneous size distribution. On the other hand, the downwind distance at which the dune size stabilised decreased with increasing injection density.

Figure 9 The mean width and dune density in 500 m cross sections averaged from measurements at the end of each year once the swarm properties stabilised. The grey areas represent 1 standard deviation. Scaled outflux was used for simulations with $q_{\mathrm{shift}}/q_{\mathrm{sat}}=\mathrm{0.1}$ and ρ₀=12, 24, and 37 km⁻² in panels (a)–(c) respectively.

As we have already discussed, for ρ₀=37 km⁻² the dune density exhibits a non-monotonic behaviour with downwind distance, but after a peak at around 2.5 km, the number density decreases, moving downwind. The same trend is observed for an injection density ρ₀=37 km⁻², with the peak moving downwind to around 3.5 km. However, for ρ₀=12 km⁻², we are able to produce a swarm which for which the dune density remains constant from around 4.5 km. In the same swarm, the size distribution stabilises at around 6 km, meaning that over the final 4 km of the simulation space both the dune density and dune size are homogeneous, as expected for real-world swarms. The size distribution of that swarm from 6 km (when it stabilises) is shown in Fig. 10 and can be fit with a log-normal distribution, as has been reported for barchan swarms in many studies (Elbelrhiti et al., 2008; Durán et al., 2009, 2011; Génois et al., 2013a; Robson et al., 2022).

Figure 10 The stable size distribution from 6–10 km in the scaled flux simulation with $q_{\mathrm{shift}}/q_{\mathrm{sat}}=\mathrm{0.1}$ and ρ₀=12 km⁻². The solid black line shows a fitted log-normal distribution.

3.1.3 Boundary condition change

Recent studies have begun to consider the impacts of future anthropogenic climate change on dune systems (Gunn et al., 2022; Baas and Delobel, 2022). Other studies have previously examined how cultivation (Long and Sharp, 1964), quarrying (Al-Dousari and Pye, 2005), and mining (Engel et al., 2018) can lead to increases in both aeolian activity and dune formation.

To explore how the TFABM might be used to study such human-induced changes, we performed simulations in which a swarm was allowed to form and stabilise under a particular combination of boundary conditions before these conditions were then suddenly changed.

In Fig. 11, we show the snapshots and longitudinal width profiles of a stable swarm exposed to a sudden decrease in the injection density. From the width profiles, one can see that after the change in boundary conditions, the swarm adjusts to reach a new steady state. In fact, the new steady state is the same as that of a swarm that had formed entirely under the changed boundary conditions. In other words, after a sufficiently long time, the properties of the swarm no longer preserve any information about the initial boundary conditions and are described solely by the new state. This challenges the idea that barchan swarms record palaeoenvironmental change (Hesse, 2009).

Figure 11 A simulated swarm and its longitudinal width profile at 300, 400, 500, and 600 years in panels (a)–(d) respectively. The swarm formed with q_shift=0.1 q_sat and a scaled outflux with an initial injection density of ρ₀=24 km⁻² until 356 years, after which ρ₀=12 km⁻². The dashed and dotted lines represent the results from simulations where ρ₀ was kept fixed at 24 and 24 km⁻² respectively, with the other model parameters unchanged.

Although there is no long-term memory, a swarm does not instantaneously adapt to the new conditions; instead, we observe a “wavefront” progressing through the snapshots shown in the left panels of Fig. 11. Comparing the position of this front (around 4 km in Fig. 11b, 6 km in Fig. 11c, and 10 km in Fig. 11d) with the longitudinal width profiles, we see that upwind of the front, the swarm has reached its new state, while downwind of the front, the swarm is still in its original state. Thus, there is a finite time during which the properties of a swarm might record changes in boundary conditions. In the case of a change in injection density, the time of the transition is the time it takes for the wavefront to travel the entire length of the swarm. Of course, other boundary condition changes which simultaneously impact all areas of a swarm (rather than just the upwind boundary, as for the case of injection) might have very different rates of transition to the novel state.

It has recently been suggested that changing boundary conditions in aeolian systems are associated with a temporary increase in the interaction density (Marvin et al., 2023). We observe an overall decrease in the rate of collision during the change in boundary condition simulations shown here; however, this change in the global collision rate simply reflects the fact that the total density of dunes is decreasing. On the other hand, the wavefront does appear to show a higher density at a very small range of downwind distances. If this is the case, then there would be a very localised increase in the interaction density, consistent with Marvin et al. (2023). However, due to the limited longitudinal range of the front, and the high variability in collision rate in the swarm as a whole, we have been unable to identified any localised increase in collisions.

3.2 Bimodal simulations

While the simulations in the previous section all simulated a unimodal wind, many real-world swarms are exposed to a secondary mode for some parts of the year (Robson and Baas, 2024b), which is something that we can simulate using the TFABM but which was inaccessible to earlier models. We performed simulations with ρ₀=37 km⁻² and q_shift=0.1 q_sat, where for three-quarters of the year the wind follows the dominant mode, set to a Gaussian distribution centred on 0° with standard deviation of 3°, while for the remaining season, the wind direction was set to a Gaussian distribution with the same standard deviation but centred on an angle θ_b for which we used 22.5, 45, and 67.5°. We chose these angles since barchan swarms have been observed in locations with bimodal wind regimes with acute angular separations but not for obtuse separations.

In Fig. 12, we show the final states of the simulations for the different modal angular separations with ρ₀=37 km⁻² and q_shift=0.1 q_sat alongside plots indicating that both the mean size and dune number stabilise during the course of the run. The number of dunes is typically less stable than the unimodal swarms, which suggests that the individual dunes are not so stable when the secondary mode is included. For all cases, we also observe that the starboard edge of the swarm (which is further upwind during the secondary mode) appears to have a higher density of larger dunes compared to the downwind (port edge) area, with this affect becoming more pronounced with increasing angular separation.

Figure 12 Final states and stabilisation of mean size and dune number for ρ₀=37 km⁻² and q_shift=0.1 q_sat, scaled outflux, and secondary wind mode for the last quarter of each year with angular separation. (a) 22.5°, (b) 45°, and (c) 67.5°.

Although q_shift is the major control on collision outcomes in the TFABM, it has also been shown that the relative asymmetry of colliding dunes (Robson and Baas, 2023) and direction of migration (Bo and Zheng, 2013) can also control barchan interactions. However, we did not observe significant changes in the relative frequencies of the different collision types as the angular separation increased. There was a slight decrease in the number of exchange collisions (62 % to 56 %) as the angular separation increased, accompanied with an increase from 5 % to around 13 %, while merging changed only by around 1 %. The rate of calving also increased as the angular separation increased, but overall it remained significantly lower than the rate of collisions, as predicted in mean field modelling (Robson et al., 2022).

It was shown in Robson and Baas (2023) that, using the Two-Flank Agent-Based Model, under a bimodal wind where the secondary mode has an acute angular separation from the dominant mode, asymmetry grows according to the Bagnold model (Bagnold, 1941; Bourke, 2010), whereas for obtuse angular, separation it grows according to the Tsoar model (Tsoar, 1984; Tsoar and Parteli, 2016; Bourke, 2010). Similar findings have also been reported using other types of modelling (Parteli et al., 2014; Lv et al., 2016).

The distributions of asymmetry are shown in Fig. 13. One would expect that under the Bagnold model for asymmetry growth (Bagnold, 1941), the starboard flank would grow predominantly under our bimodal winds, leading to $W_{\mathrm{p}}/W_{\mathrm{s}}<\mathrm{1}$, while under the Tsoar model we would see $W_{\mathrm{p}}/W_{\mathrm{s}}>\mathrm{1}$. In Fig. 13e), we see that in all of the bimodal runs, the average of $W_{\mathrm{p}}/W_{\mathrm{s}}$ is less than one, indicating that overall the asymmetry of dunes is growing predominantly, according to the Bagnold model (Bagnold, 1941; Bourke, 2010). We do observe, however, that the mean asymmetry varies with the angle of the secondary mode. Specifically, as the angle increases, the mean decreases from unity. However, none of the swarms become as asymmetric as those observed in Mauritania and Mars zone 1 (Robson et al., 2022; Robson and Baas, 2024a, b). The standard deviation of the distribution is not significantly affected by the change in angle and is primarily controlled by q_shift, as shown in Figs. 6 and 7.

Figure 13 Histograms of dune asymmetries defined at the ratio of port and starboard flank widths are shown in panels (a)–(d) for θ_b=0, 22.5, 45, and 67.5°, with $q_{\mathrm{sat}}/q_{\mathrm{shift}}=\mathrm{0.1}$ and scaled outflux. The dashed and dotted lines show the distributions for real-world swarms in Tarfaya and Mauritania. Panels (e) and (f) show the mean and standard deviations of the distributions shown in panels (a)–(d) compared to four real-world barchan swarms.

A feature of barchan swarms that we have not yet discusses is that a downwind dune tends to align with the horn of its upwind neighbour (Bagnold, 1941; Elbelrhiti et al., 2008). In the earliest barchan ABM, this was evidenced by peaks in an angular correlation function (Lima et al., 2002); however, such peaks have not been reported widely for real-world barchan swarms. On the other hand, Elbelrhiti et al. (2008) showed that this affect could be measured by looking at the normalised lateral offset of dunes and their nearest downwind neighbours, which showed double-peaked distributions (with peaks $\sim \pm \mathrm{1}$ indicating alignment downwind of the horns) when the downwind neighbour was smaller than the upwind dune and a flat distribution when the downwind dune was larger than its neighbour. In Fig. 14, we show the distributions of the lateral offset of dunes and their nearest downwind neighbour in each case normalised by the width of the upwind dune for each of the six barchan swarms from Robson et al. (2022) and Robson and Baas (2024a, b). The three Tarfaya swarms (Fig. 14a–c) show the same double-peaked behaviour as reported in Elbelrhiti et al. (2008), although the distribution for which the downwind neighbour is larger is not as uniform as previously reported. The Mauritanian swarm (Fig. 14d) shows an approximately uniform distribution for both cases, with only small peaks at ±1. Finally, the Martian swarms (Fig. 14e and f) show double-peaked distributions but with the peaks occurring at ±2, indicating that the barchans do not occur with any lateral overlap with their upwind neighbour.

Figure 14 The normalised lateral offset of the nearest downwind neighbours (defined as the lateral offset divided by the width of the upwind dune) in the populations of (a) Tarfaya 1, (b) Tarfaya 2, (c) Tarfaya 3, (d) Mauritania, (e) Mars 1, and (f) Mars 2 (available at https://doi.org/10.5281/zenodo.10848514) (Robson and Baas, 2024a).

In Fig. 15, we show lateral offset of the nearest downwind neighbours in the simulated swarms with θ_b=0, 22.5, 45, and 67.5°. Unlike the real-world results, we see the double-peaked shape in both cases where the downwind dune is smaller and larger than the upwind dune. Increasing the angle between the wind modes has the affect of widening the peaks in both distributions. We also see that for greater angular separation the distributions for the cases of smaller and larger dunes begin to become offset from one another, whereas for unimodal wind they are perfectly aligned.

Figure 15 The lateral offset of the nearest downwind neighbour of dunes normalised by the width of the upwind dune for scaled outflux $q_{\mathrm{shift}}/q_{\mathrm{sat}}=\mathrm{0.1}$. θ_b=0, 22.5, 45, and 67.5° in panels (a)–(d) respectively.

4 Discussion

The efficiency of the TFABM has meant that we have been able to perform a significant number of swarm-scale simulations and investigate the impact of a number of parameters on the types of swarms that are produced. Previous ABMs have not been able to produce the longitudinally homogeneous size distributions (Lima et al., 2002; Durán et al., 2011; Worman et al., 2013) which are expected for real-world swarms (Elbelrhiti et al., 2008; Durán et al., 2011). In contrast, we have shown that for several combinations of model parameters the TFABM is able to produce longitudinally homogeneous swarms.

The phenomenological parameter q_shift, which relates to the lateral motion of sediment within an asymmetric barchan and, therefore, controls the stability of asymmetric dunes, is found to be the most important factor in generating stable dune sizes through its influence on the collision dynamics. However, the types of collisions required to produce homogeneity depend upon whether the outflux of sand from the dunes is scaled with the influx that they receive. If the outflux is saturated across the horns, there is an unstable equilibrium size which should lead to runaway growth or shrinking (Hersen et al., 2004; Hersen and Douady, 2005; Worman et al., 2013). We find, in such cases, that homogeneity results from a dominance of fragmentation collision which occurs for q_shift≥0.1 q_sat. Fragmentation collisions serve to counteract the runaway growth of dunes which become larger than equilibrium, and so the homogeneity is seen as a balance between the rate of growth of the large dunes due to flux and the rate of collisions resulting in fragmentation.

Fragmentation collisions are able to reduce the size of a large dune through the creation of at least one additional new dune. However, these dunes are often too small to persists on their own, as has previously been reported using continuum simulations (Durán et al., 2005). The formation of ephemeral small dunes, whose sediment is converted to free flux when the dunes become too small to persist, is what leads to the decrease in dune density in the unscaled flux simulations. Without additional effects, such as the spontaneous creation of sand patches which evolve into dunes (Delorme et al., 2023), this decrease in density is inevitable if fragmentation collisions are required to stabilise dune size.

If, on the other hand, the outflux from dunes scales with the influx they receive then, rather than an equilibrium size for dunes, there is an equilibrium flux for which a barchan of any size would be stable. In this case, longitudinal homogeneity can be achieved if collisions typically drive dunes to a particular size. Size selection has been shown to be possible using exchange collisions only (Durán et al., 2009; Robson et al., 2021), with stronger peaks indicating the higher the relative prevalence of exchange over merging and fragmentation (Robson et al., 2022). The explains why we observe stable sizes in TFABM swarms with scaled outflux when q_shift=0.1 q_sat, since at this value exchange collisions dominate, as shown in Fig. 3. At q_shift=0.1 q_sat, over 60 % of collisions result in exchange, with ∼35 % merging and the remainder being fragmentation. This also explains why, although exchange dominance leads to a constant size distribution which can be described by a log-normal distribution (Fig. 10), the overall density typically decreases with downwind distance, since merging collisions (which decrease dune number) occur significantly more frequently than fragmentation collisions (which increase dune number). For the values of q_shift we tested, only if the density of dunes is low, which reduces the number of collisions, do we observe a swarm that is homogeneous in both size and density (Fig. 9).

Looking at Fig. 3, it is possible that between $q_{\mathrm{shift}}/q_{\mathrm{sat}}=\mathrm{0.1}$ and 0.15 there exists a point at which the rates of fragmentation and merging are equal, allowing the homogeneity of size and density to coexist at higher densities; however, we were unable to find such a point with a manual search of the parameter combinations.

In order to constrain the range of parameters in the TFABM to reproduce realistic swarms, the longitudinal size and density profiles of real-world swarms require further investigation. In Fig. 1 of Durán et al. (2011), the size profiles of four swarms in Tarfaya are shown to be longitudinally homogeneous for between 3 and 13 km, but there is no discussion of dune density. Qualitatively the swarms shown in that work display non-monotonic density over their downwind extent. Elbelrhiti et al. (2008) is the only other study which has explored longitudinal profiles of size and density over 4–5 km. Some of the profiles show sizes which remain constant over the entire distance, while others shown an initial increase in size followed by either a plateau or decrease in size (similar to many of our simulations). Similarly, the density profiles for some of the locations are flat, while others display non-monotonic behaviour, including an increase to a peak followed by a decrease, as observed in many of the swarms we presented here. Interestingly, in paragraph 8 of chapter 14 in Bagnold (1941) (p. 218 of the 1954 reprint), it was claimed that it is “common for barchan belts to terminate because the frequency of their occurrence dwindles to nothing”. In other words, the density of dunes should decrease significantly towards the end of a swarm. Bagnold (1941) also suggests that this could be mitigated if additional sediment becomes available following that supplied at the upwind boundary. Indirect evidence for the stability of dune density comes from the Tarfaya zones of Robson et al. (2022) and Robson and Baas (2024a), which could be interpreted as several areas within a large swarm and which have very similar sizes and densities to one another. However, much more work is needed to establish how barchan size and density behave in natural swarms.

Assuming that swarms are longitudinally homogeneous, at least in size, the values of q_shift which produce homogeneity are 0.15 q_sat and 0.1 q_sat in the unscaled and scaled flux cases respectively. These were the same values of q_shift which were shown to reproduce over 80 % of collisions from cellular automata (Katsuki et al., 2005), and water tank experiments (Assis and Franklin, 2020), as shown in the supplementary material of Robson and Baas (2023). Thus, it seems that collisions produced by other techniques coincide with those that produce longitudinally stable size distributions in swarms.

On the other hand, we have shown in Fig. 7 that there are significant differences in the width of asymmetry distributions in real-world swarm, and consequently, a wider range of q_shift is required to reproduce all of the different distributions in observed swarms. This suggests that it might be necessary to decouple the dynamics of collisions and asymmetry growth in the TFABM to capture a wider range of swarm-scale properties. Figure 13 also shows that, although bimodal wind regimes can drive dunes away from symmetry, we are unable to produce the extreme asymmetries observed in some locations. In order to reproduce such asymmetries, one could adjust the shape parameters in the TFABM which currently prevent dunes from becoming too asymmetric according to Eq. (5).

Finally, we have been able to reproduce the two-peak distribution of the alignment of dunes and their downwind neighbours, which is something that has been observed previously for real-world swarms (Elbelrhiti et al., 2008). However, there are some differences between the observed distributions and those produced using the TFABM; the most significant is that in the simulated swarms, there is little difference between the shape of the distribution for smaller and larger downwind dunes. It is possible that one must include the lateral diffusion of sand flux (Schwämmle and Herrmann, 2005) or the deflection of sand flux around a barchan so that the effective shadow cast by a dune only lasts a finite distance. For instance, it has been observed that erosion downwind of the bulk of a fixed model of a barchan in a wind tunnel persists only to around 20 dune heights, depending upon the wind speed and size of the dune, with deposition resuming beyond this point (Cai et al., 2021). We do observe, in Fig. 14, a slightly different behaviour than Elbelrhiti et al. (2008) in some cases (especially the Martian swarms), and so future work should look further into the alignment of barchans in swarms.

5 Conclusions

We have demonstrated that the Two-Flank Agent-Based Model is capable of generating swarms of barchans which have longitudinally homogeneous size distributions which had not been produced using previous agent-based models but which are known to be the case for some real-world barchan swarms. The model and results presented here, therefore, represent a step forward in our understanding of the behaviour of these complex swarms. We observed that, depending on whether or not outflux of sand from dunes scales with the influx, different relative frequencies of collision types are required to produce longitudinal homogeneity. However, in both cases, the necessary collision dynamics coincide with those reported in numerical and experimental studies and are produced by the TFABM with q_shift in the range 0.1–0.15 q_sat. It may be possible to test this prediction using more involved computational models or water tank experiments.

Longitudinal stability of dune density is a rarer occurrence in the simulated swarms; however, it is possible to create sparse swarms exhibiting this kind of homogeneity. Current quantitative evidence is lacking for the assessment of whether such stability is indeed a universal property of barchan swarms in nature. If it is a universal property, then it may be necessary to incorporate sand patch formation or additional sediment supply in the TFABM to produce dense homogeneous swarms.

Unlike with previous agent-based models, the Two-Flank Agent-Based Model is capable of simulating swarms under bimodal winds, which revealed that bimodality may be responsible for asymmetry distributions centred away from unity, and it also changes in the alignment of dunes with their neighbours. This adds to our understanding of how dune interaction shape the properties of swarms.

Being able to simulate more realistic swarms using an agent-based model presents many exciting opportunities for future work, including the investigation of changing boundary conditions. Preliminary results indicate that there is a finite transition time for barchan swarms, after which historical boundary conditions are no longer reflected in the properties of swarms. Future work could explore changing boundary conditions further to understand how swarms will be affected by anthropogenic change. It may also be possible to use the TFABM to train machine-learning algorithms to trace individual barchans through satellite imagery of dense swarms.