A Simple Agent‐Based Model That Reproduces All Types of Barchan Interactions

We introduce a novel agent‐based model for simulating interactions between migrating barchan dunes. A new two‐flank representation of barchans allows modeling of bedform asymmetries that are intrinsic to collision dynamics but have not been explored before. Although simple compared with real‐world barchans or those in continuum and cellular automata simulations, all known barchan behaviors emerge from the rules of our model. In particular, the two mechanisms for asymmetry growth in bimodal winds are observed and qualitatively agree with existing theories. We also reproduce the emergence of calving and all types of collisions that have been reported in reductionist models, water‐tank experiments, and field observations. The computational efficiency of the new model, compared with continuum simulations, enables the simulation of large swarms of dunes while maintaining the complex phenomenology of these bedforms, some of which has been lacking in previous agent‐based models.

Plain Language Summary Barchans are naturally occurring sand dunes found in regions where the wind direction is near-constant, and the overall supply of sand is low.Because of these conditions, barchans migrate very quickly resulting in collisions between the dunes.Interactions between dunes also occur as sand streams off upwind dunes and is absorbed into downwind barchans.In this work, we present an agent-based model which treats the dunes themselves as the elementary objects, rather than the sand and airflow.Such models are capable of simulating large populations of dunes.We model barchans as comprising two flanks which can grow semi-independently.With this structure, we can replicate complex phenomena, including dune asymmetry due to varying winds, restoration of symmetry under a constant wind, and the spontaneous breakup of dunes due to strong winds in directions close to 90° from the usual.These phenomena were inaccessible to previous agent-based models of barchans.We are also able to reproduce all of the different types of collision which have been observed in lab experiments and more computationally intensive models.The new model, therefore, represents an improvement on previous agent-based models while remaining computationally cheap enough to simulate large populations.

10.1029/2023GL105182
2 of 12 & Herrmann, 2003; Worman et al., 2013).The most realistic of these models are agent-based models, shedding light on how processes such as exchange collisions, fragmentation collisions, and spontaneous calving affect size distributions and cluster formation in barchans swarms (Durán et al., 2011;Diniega et al., 2010;Génois et al., 2013;Lima et al., 2002;Parteli & Herrmann, 2003;Worman et al., 2013).Past agent-based models have been confined, however, to incomplete sets of phenomena for the dunes, resulting in unrealistic properties such  as spatially heterogeneous size distributions in the simulated swarms (Durán et al., 2011;Worman et al., 2013) or other unrealistic features such as dunes being injected at a maximal rather than minimal size (Génois et al., 2013).Existing models have furthermore considered only idealized barchans and wind regimes, limiting the range of properties that can be investigated.
In this work, we introduce a novel agent-based model in which a simple set of rules reproduce a wide range of the behaviors exhibited by real-world barchans, and those that have been simulated using microscopic models or reproduced in the lab.

Model Structure
Our model shares many similarities with previous agent-based models, including the size dependence of migration rate and the volume change due to incoming and outgoing sand fluxes.The major simplification (and advantage) of all agent-based barchan models is to treat the dunes as elementary objects (or agents) by reducing the detailed local interactions between (secondary) airflow dynamics, sand transport and dune morphology into simple rules governing the size-evolution and migration of the objects as a whole (Durán et al., 2011;Lima et al., 2002;Worman et al., 2013).
Our model differs from the previous works by using a more versatile shape of dune allowing us to account for the fact that barchans in real-world swarms are rarely symmetric (see Figure 1c) (Bourke, 2010) and variations in wind occur due to synoptic variations on daily and seasonal scales (Elbelrhiti et al., 2005).

Morphology
Locations of dunes are stored as the position of the toe.Each dune is divided into a left and right flank with widths W l and W r .The flanks adjoin along the central axis for the length of the main body of the dune,   = 1 2 1( + ) .The parameter λ 1 is the proportionality constant relating the length and width of the main body.Based upon field data (Sherman et al., 2021) we set λ 1 = 1.5.The dunes' central axes are fixed parallel with the dominant wind direction, the positive x-direction.Downwind of the body, each flank extends into a horn, with the streamwise distance from toe to the horntip proportional to the flank-width L l,r = λ 2 W l,r .From field data (Sherman et al., 2021) we set λ 2 = 2.5 The y-coordinate of the horn-tip extends half a horn-width H l,r from the respective edge of the dune.We use an existing expression for the scaling of horn and flank-widths with α = 0.05, and Δ = 4.6 m (Hersen et al., 2004;Worman et al., 2013).All of these properties are shown in Figure 1b.Finally, we assume a pyramidal shape for the bodies and horns of each flank giving flank-volumes proportional to the cube of the width (Elbelrhiti et al., 2008) where λ 3 is the height-width ratio.

Sand Flux
Barchans absorb incoming sand flux across their entire upwind margin and emit flux at a saturated rate q sat from their horns (Bagnold, 1941;Hersen et al., 2004;Lima et al., 2002).In the simulations, we set q sat to 100 m 2 year −1 .The upwind boundary of the domain receives a free flux input, q 0 < q sat , which propagates in the wind direction until it is absorbed by a dune.We fix q 0 = 0.2q sat unless otherwise stated.Since the rate of flux propagation is much higher than the dunes' migration speeds, unobstructed flux is assumed to travel from the upwind to downwind boundary within a timestep (Durán et al., 2011;Lima et al., 2002;Worman et al., 2013) set to 0.001 years, around 9 hr.
In addition to flux in the wind direction, we introduce a lateral flux q shift to account for the exchange of sediment between the flanks, perpendicular to the driving wind.This parameter may represent restructuring of the dune bulk through avalanches and reputation, sediment splash, and short-period gusts normal to the primary wind.
Since the magnitude of any turbulent eddies will depend upon the wind strength, we assume that q shift ∝ q sat with a proportionality constant of 0.2 used unless otherwise stated.This value was chosen as it enables all the interactions described below and is also consistent with typical isotropic turbulence intensity in atmospheric boundary layers (Kaimal & Finnigan, 1994).In the absence of detailed lab or field measurements of asymmetry reduction and collisional phase-spaces, we are unable to constrain this parameter at this time (see Supporting Information S1 for a more detailed discussion).
For a symmetric dune, lateral flux has no net effect but when one flank is larger than the other, q shift leads to a transfer of material from the larger to the smaller flank at a rate proportional to the difference in wind-parallel projections of the flank-lengths.Combining streamwise and spanwise fluxes, flank-volumes evolve according to where q in is the received influx, θ is the wind direction, and  W is the flank-width perpendicular to the wind.In the absence of upwind dunes, q in = q 0 .In the presence of upwind dunes q in may be modified due to absorption by and emission from other barchans.Additionally, when wind is not in the dominant direction, one flank may prevent the other from absorbing sand flux.

Maximum Asymmetry and Calving
Drivers of asymmetry include: wind variations, asymmetric sand-supply, collisions, and topography (Bourke, 2010).In our model, as a dune becomes more asymmetric, the rate at which material is lost from the larger flank to the smaller due to q shift increases (final term in Equation 3).If allowed to increase freely, eventually the material loss from the large flank due to q shift will exceed the material this flank is losing from its horns (the second term in Equation 3).In such a case, we argue the agent is no longer behaving like a barchan.Asymmetry dominance occurs when the ratio of the larger flank width to the smaller flank width is greater than a maximum asymmetry threshold γ c given by where W min is the smaller flank width.To ensure that all dunes behave as barchans, any dune that exceeds the threshold breaks apart into two dunes.We call this process calving.When calving occurs, the two flanks become separate, symmetric dunes with volumes conserved.The y-positions of the new dunes are set to the y-coordinates of the flanks' centres of mass prior to calving.The x-positions are determined by placing the smaller product of calving just downwind of the larger dune and then moving the two so that the overall centre of mass of the dunes remains unchanged.

Migration Rates
Barchan migration rates are proportional to the saturated sand flux and increase with decreasing dune size (Bagnold, 1941;Elbelrhiti et al., 2008;Hersen et al., 2004).It has also been shown that the migratory behavior is unaffected by asymmetry (Bourke, 2010).Dunes in our model migrate in the wind direction at a rate given by where c is a model parameter connected to the wind speed-up over the slope of the dunes.We set c = 50 in all simulations presented here unless otherwise stated, consistent with previous agent-based models (Durán et al., 2011) and continuum simulations (Durán et al., 2010).

Collisions
Size-dependent migration rates ensure that collisions are ubiquitous in swarms.As a smaller upwind dune catches a downwind dune, the two merge to form a complex intermediate bedform.In some cases, this reverts to a typical barchan shape, however the combined bedform may, instead, be unstable and break apart, resulting in several 10.1029/2023GL105182 5 of 12 barchans.The two-flank structure of dunes in our models enables a simple collision algorithm which produces a wide range of interaction outcomes.
There are three centres of mass (CoMs) associated with dunes in our model; one for each flank and one for the entire bedform as a whole.We define a collision as occurring when one or more of the CoMs of the upwind dune is contained within the basal area of the downwind dune meaning that dunes can overlap for some time before a collision is triggered accounting for the time taken for the intermediate stages of real-world barchan collisions.
Once the collision condition is met, the outcome is determined as follows: 1. Calculate volume V col of intersecting flanks.2. Convert to width W col using Equation 2.

Non-intersecting flanks join intersecting ones if
4. Flanks where W col /W flank > γ c (W flank ), become separate symmetric barchans in the same manner as calving.
The positions of the new dunes are determined in a similar way to calving.The dune that formed from the intersecting flanks may be asymmetric based on the relative volumes to the left and right of the CoM prior to collision.Total volume and CoM are preserved during collisions.

Isolated Dune Phenomena
The most important difference between our model and previous works (Durán et al., 2011;Génois et al., 2013;Worman et al., 2013) is the two-flank morphology of the dunes and the way that these flanks interact with one another.The two-flank structure is capable of reproducing complex phenomena observed for isolated dunes exposed to different wind regimes.

Asymmetry Growth and Calving
Several drivers of barchan asymmetry have been identified (Bourke, 2010;Tsoar & Parteli, 2016) with bimodal winds most often cited.Two conceptual models have been proposed: the Bagnold model (Bagnold, 1941) and the Tsoar model (Tsoar, 1984).Due to the complex nature of real-world systems, field evidence remains inconclusive (Bourke, 2010).However, cellular automata and continuum simulations have revealed that the effects depend upon the angular divergence between the primary and secondary winds, their relative strengths, and their duration (Lv et al., 2016;Parteli et al., 2014).
In Figures 2a-2c we show the behavior of an initially symmetric barchan subject to bimodal wind regimes, simulated by assigning a fixed probability of 50% for the secondary wind direction each timestep.Changing the probability and saturated fluxes affects the rate of asymmetry growth and calving but does not change the overall phenomenology.For acute angular wind divergences, the upwind limb of the dune under the secondary wind becomes larger than its counterpart, similar to Bagnold's conceptualization (Bagnold, 1941).An obtuse angular wind divergence leads to the growth of the downwind limb, similar to Tsoar's model (Tsoar, 1984).Similar phenomenology has been observed in continuum simulations (Figure 3

Asymmetry Reduction
The bimodal winds simulated above may represent seasonal or daily variations, or storms as in Bagnold's model (Bagnold, 1941).In the case of storms, after the short period of secondary wind direction, the now asymmetric barchan is then re-exposed to the dominant wind.Figure 2d shows the evolution of an asymmetric dune under a unidirectional wind, with asymmetry reducing and the central axis shifting laterally over time.Similar to our observations, Parteli et al. (2014) showed how an asymmetric dune reverts to a symmetric shape.The lateral shifting of the central axis is reminiscent of cellular automata simulations of dunes subject to off-centre fluxes in Katsuki and Kikuchi (2011).

Interaction Phenomena
Interactions between barchans in a swarm often play a more important role than external drivers, therefore it is vital that an agent-based model can replicate the phenomenology of those interactions.

Sand Flux Interactions
Absorption of incoming flux and emission from a dune's horns is a key mechanism for interactions between barchans.In Figure 3a we show an example where sand flux absorption by an upwind dune causes a loss of material in the downwind dune.At first, the smaller upwind dune is catching up with the downwind dune, but the downwind dune soon shrinks and accelerates, moving away from the upwind dune.This effect has been reported in experiments by Assis and Franklin (2020) and Bacik et al. (2020).
The downwind dune in the previous example was able to maintain its size toward the end of the simulation by realigning itself with the horn of the upwind barchan via lateral shifting which occurred due to asymmetry induction.The promotion of asymmetry by the upwind dune can also cause the downwind dune to calve, as shown in Figure 3b.Induced calving has been observed in nature (Elbelrhiti, 2012;Elbelrhiti et al., 2008), and in water-tank experiments (Assis & Franklin, 2020;Assis et al., 2022).The phenomenon is also similar to an interaction described in earlier experiments (Endo et al., 2004), cellular automata models (Katsuki et al., 2011) and conceptually (Bagnold, 1941).
2. Fragmentation-more dunes are produced as outputs than were involved in the collision.
3. Exchange-where material is exchanged between dunes but the total dune number does not change.
Finally, fragmentation collisions have been observed numerous times (Assis & Franklin, 2020;Assis et al., 2022;Durán et al., 2005Durán et al., , 2009Durán et al., , 2011) ) with various dynamics identified, including "breeding" and "budding" (see Figures 1b-1d in Durán et al. (2005)) and "fragmentation-exchange" (see Figure 1e in Assis and Franklin (2020)).The "split" interaction in Endo et al. ( 2004) is also similar, however in that interaction no collision takes place and so it is more similar to the induced calving we described above.Fragmentation processes differ mainly in the intermediate stages and appear similar when only the initial and end states are considered.Fragmentation collisions producing either three or four outputs have been observed (see Figures 1d and 1c in Durán et al. (2005) respectively).
The collision algorithm we described in Section 2.5 can directly result in one, two, or three outputs, while even four outputs may arise if the dune which forms from the intersecting flanks is sufficiently asymmetric that it subsequently calves.Examples of collisions producing the different numbers of outputs are shown in Figures 3c, 3e, 3f, and 3g.Some collision outputs can be asymmetric, noticeably in Figures 3f and 3h.Additionally, the second type of exchange collision described above, with dunes intersecting only at the edge of their horns, can also occur in our model since a collision is not said to occur unless the centre of the flank is contained within the other dune.This type of glancing interaction (or lateral linking) is shown in Figure 3d and closely resembles those shown in Figure 1 of Gay Jr (1999).The result is not that dissimilar to an exchange collision but the dune tracking shown in Figure 3 reveals the distinction.
The collisions described above involve symmetric dunes, however, when the colliding dunes are asymmetric, the result of the collision also depends on which flanks are larger, as shown in Figure 3h.
As well as asymmetry, the relative sizes of the colliding dunes (Assis & Franklin, 2020;Bo & Zheng, 2013;Diniega et al., 2010;Durán et al., 2005Durán et al., , 2009Durán et al., , 2011;;Endo et al., 2004;Katsuki et al., 2011) and their lateral offset (Assis & Franklin, 2020;Bo & Zheng, 2013;Durán et al., 2009Durán et al., , 2011;;Katsuki et al., 2011) are known to impact collision dynamics.To isolate these dependencies we explored the phase-space of simulations starting with two dunes in series, with wind direction θ = 0°.If long-range interactions are permitted then the initial longitudinal displacement of the dunes is an important factor in determining the output.To focus on the phase-space of the collision algorithm alone we, therefore, rescaled q sat , q 0 , and q shift by a factor ηll1 while simultaneously rescaling c by 1/η to maintain migration rates.We define the lateral offset as the distance between the toes of the dunes in the y-direction normalized by the combined width of the adjacent flanks.Since sand-flux effects are nullified, the initial streamwise separation of the dunes is not important, however, for clarity, we set the initial separation such that the downwind dune is initially 1 m downwind of the tip of the horns of the upwind dune.
In Figure 4 we show the outcomes of collisions as a function of lateral offset and initial size ratio for different values of q shift and sizes of the smaller dune.The only model parameter which influences collision type is q shift which controls the asymmetry threshold γ c .High q shift corresponds to low γ c , which means that fragmentation collisions make up a larger proportion of the phase-space in the lateral offset-size ratio plane.As q shift decreases, dunes can become more asymmetric without fragmenting and hence this type of collision becomes rare and exchange and aggregation collisions dominate.No fragmentation collisions were observed for q shift /q sat = 0.05.Previous works have given details of only a very limited region of the full collision phase-space (Durán et al., 2005(Durán et al., , 2011) ) or else, been carried out in water-tank experiments in which turbulent wake effects cannot be de-coupled from the collision processes (Assis & Franklin, 2020;Endo et al., 2004;Katsuki et al., 2011).These limitations mean that we cannot infer a value of q shift by comparing directly with previous results.However, in  Supporting Information S1 we compare the outputs of our collisions to those observed in two previous studies, one computational (Katsuki et al., 2005) and one experimental (Assis & Franklin, 2020).In both cases, our collision rule reproduces the vast majority of collisions and we also show that it performs better than the collision rules used in previous agent-based models (Durán et al., 2011;Génois et al., 2013;Lima et al., 2002;Worman et al., 2013).
Regardless of the value of q shift , merging collisions occur for dunes that are nearly perfectly aligned as, when aligned dunes collide, all four flanks intersect.One class of exchange collisions represents the same case, but where the maximum asymmetry threshold is sufficiently close to unity that the asymmetric merged dune immediately splits (this occurs only for high q shift ).Exchange collisions also occur when three flanks intersect at the time of the collision, those three flanks form one output dune while the remaining flank forms another.Fragmentation collisions most often occur when the lateral offset is high so that only two of the flanks overlap during the collision.However, for high q shift , the larger product of exchange can be sufficiently asymmetric to immediately calve, resulting in fragmentation collisions at lower values of the lateral offset.
While the dependence of the phase-space on q shift is easily interpreted, the dependence on absolute dune sizes is more complex.Dependence on size has been identified in water-tank experiments (Endo et al., 2004)and its effect may be masked by definitions of lateral offset which only depend upon the size of one of the dunes (Durán et al., 2011).In our model, absolute size affects the outcome of collisions due to the non-linear scaling of volume on width which affects whether or not the maximum asymmetry threshold is exceeded.Generally, we observed in Figure 4 that four-body fragmentation collisions are more common for collisions between larger dunes, whereas a greater proportion of three-body fragmentation occurs for smaller dunes.At high values of q shift , however, the situation is complicated by secondary fragmentation (described above), and dependence on absolute size in collision dynamics requires further investigation.

Conclusions
The agent-based model we propose here remains very simple compared to realistic but computationally expensive continuum and scale-coupled cellular automata.Despite this, we demonstrate a wide range of complex phenomena emerging from the small number of rules, including: 1. Asymmetry induced by bimodal winds according to both Bagnold's and Tsoar's hypotheses.2. Transience of asymmetry under a unimodal wind direction.
3. Lateral displacement of asymmetric dunes.4. Spontaneous fragmentation (calving) due to oblique winds. 5. Long-range repulsion due to sand flux absorption.6. Calving induced by the presence of an upwind dune.7. One, two, three, and four-output collisions, which can also cause asymmetry.8. Merging collisions predominantly for aligned dunes.9. Temporary lateral linking via the horns.
These phenomena have all previously been reported in field studies, continuum simulations, cellular automata, and laboratory experiments.Many of these behaviors could not be simulated in previous agent-based models since they focused on symmetric dunes under a unidirectional wind.Furthermore, whereas the collision rules in earlier models either omitted some known collision types or introduced many additional parameters which did not have an obvious physical interpretation, in our new model, the collision, calving, and asymmetry behaviors are controlled mainly by a single parameter, the lateral flux, q shift .There are no existing measurements of the rate of asymmetry reduction under a unidirectional wind, which would provide an explicit verification of q shift .Furthermore, the collision phase spaces mapped in existing studies do not have sufficient detail to infer a value of q shift .Future work should endeavor to investigate this parameter and better constrain the links between asymmetry and collision outcomes.
There remain some aspects of real-world swarms that cannot be studied with this model, such as variations in wind forcing within the swarm.However, the wide range of phenomena that we do observe means that the model represents a significant improvement from previous works and will allow for the study of previously inaccessible properties such as bedform asymmetry and spatial distributions within swarms containing thousands of dunes on Earth and on Mars.

Figure 1 .
Figure 1.(a) A section of a barchan swarm near Tarfaya, Morocco.Imagery from 2023 CNES Airbus courtesy of Google Earth.(b) A schematic showing the morphological characteristics of the simulated dunes.(c) An asymmetric barchan from the Tarfaya swarm in panel (a) overlaid with the asymmetric simulated dune from panel (b).Barchan imagery from 2023 CNES Airbus courtesy of Google Earth.
in Parteli et al. (2014)) and cellular automata (Figure4inLv et al. (2016)).If the angular separation, duration or strength of the secondary wind mode is large enough, we observe that the dune becomes sufficiently asymmetric to induce calving (Figure2b), as also observed in continuum simulations (Figure3of Parteli et al. (2009) and Figure 4 of Parteli et al. (2014)) under the same conditions.

Figure 2 .
Figure 2. (a) Growth of upwind limb under secondary mode with angular separation 40°.(b) Angular separation of 80° leading to growth upwind limb and eventual calving.(c) Growth of downwind limb when angular separation is 120°.(d) Reduction of asymmetry and lateral movement of central axis of asymmetric dune under unidirectional wind in the positive x-direction.

Figure 3 .
Figure 3. Different interaction phenomena observed in the model are shown with ID numbers displayed representing the order in which dunes were created.(a) An upwind dune causes a larger downwind dune to decrease in size and speed up such that a collision does not occur.(b) The presence of an upwind dune results in asymmetry growth and calving of a downwind dune.(c) An exchange collision.(d) Two dunes overlap but do not collide.(e) A merging collision.(f) A fragmentation collision resulting in three output dunes.(g) A four-output fragmentation collision.(h) Collisions between two pairs of dunes with the same total width and asymmetry but where, in one pair, the larger flanks are on the same side of the dunes and, in the other pair, are on opposite sides.Please see Supporting Information S1 for gif-animations of all simulations.

Figure 4 .
Figure 4. Collision type as a function of lateral offset and initial width ratio for different sizes of the smaller dune and values of the lateral sand flux q shift .