Internal sedimentary structure of linear dunes modelled with a cellular automaton

Linear dunes are the most common type of dune found on Earth and exist on several extra‐terrestrial bodies, but despite this abundance their internal stratigraphy has not been commonly agreed. A cellular automaton is deployed to simulate the development of linear dunes, starting from a flat bed, under bi‐modal oblique wind regimes of varying degrees of asymmetry. The internal stratigraphy of the linear dunes is monitored by keeping track of (buried) erosion surfaces, avalanche deposits and vertical accumulation, as well as the age of last subaerial exposure of the sediments. The simulations show the initial pattern‐coarsening of a network of small dunes into fewer larger longitudinal ridges via bedform interactions and Y‐junction dynamics. Three newly recognized types of bedform interaction are identified that relate to initial Y‐junction dynamics: longitudinal crest‐splitting, which creates free dune tips that can interact with adjacent dunes, and laterally oscillating interactions that lead to ephemeral Y‐junctions (normal or reverse). The results show that these three bedform interactions leave no persistent signatures in the stratigraphic record. However, a further three bedform interactions involving the superposition of one dune onto another – merging, cannibalizing and repulsion (known from transverse dune field dynamics) – do leave specific evidence in the internal stratigraphy of the remaining dune, a buried interaction surface at a specific inclination. The preservation potential of this interaction surface varies between the three types. After the initial pattern‐coarsening phase, the linear dunes become larger and more independent and their crest orientation follows the net resultant transport direction. The stratigraphies of mature dunes under wind regimes of differing asymmetry show that under (nearly) symmetrical winds the dune accumulates mainly vertically, with strata dipping parallel to the flanks, while under more asymmetrical wind regimes the internal stratigraphy resembles that of transverse dunes.


INTRODUCTION
It is commonly recognized that linear dunes are the most abundant dune type on Earth (Lancaster, 1985;Craddock et al., 2015) and exist on extra-terrestrial bodies like Titan and Mars (Lorenz et al., 2006;Rubin & Hesp, 2009;Radebaugh et al., 2010). Linear dunes are associated with bimodal wind regimes, with alternating winds from oblique directions on either side of the crest. They are usually categorized as longitudinal dunes (Hesp et al., 1989), where the angle between dune crest and resultant sand transport direction is less than 15 degrees (Hunter et al., 1983), while in terms of morphology some oblique dunes can also be classed as linear (Rubin & Hunter, 1985). Compared with other dune types, they better preserve earlier sedimentary records in their dune bodies (Lancaster, 1995;Munyikwa, 2005;Telfer et al., 2010). Therefore, their chronology records have long been reconstructed with optically stimulated luminescence dating (OSL) (Duller & Augustinus, 2006;Fitzsimmons et al., 2007;Bristow et al., 2007;Stone & Thomas, 2008), and interpreted with a one-dimensional model (Telfer et al., 2010). Here a cellular automaton (CA) model for dune dynamics is developed to map the sedimentary structure of linear dunes and a few typical dune behaviours are simulated with this model.
Despite the abundance of linear dunes and many field studies and dating work, their dynamics and internal sedimentary structures are not yet resolved; for example, the relative importance of lateral migration in stratigraphic signatures. Evidence has been found to support migration (Bristow et al., 2005) as well as non-migration (Tsoar et al., 2004) theories. The debate is largely caused by incomplete understanding surrounding the lee side wind patterns, together with challenges of accessing the buried architecture, while linear dune evolution takes place over timescales that are too long to be monitored in real-time. As a result, current understanding of linear dune development and movement is incomplete. For the sedimentary structure, traditional in situ approaches, such as digging, observing exposed stratigraphy (e.g. McKee & Tibbitts, 1964;Hunter, 1977;Tsoar, 1982) and ground penetrating radar (GPR) or sonic techniques (Bristow et al., 2000(Bristow et al., , 2005Vriend et al., 2007), have been limited by depth and cannot always cover the full history of the dune from inception.
The formation of linear dunes was first attributed to roll vortices by Bagnold (1941), based on the symmetry of the dune cross-sections, as evidenced by some field observations (e.g. McKee & Tibbitts, 1964;Hanna, 1969). Since then, remote sensing (Hesp et al., 1989), flume experiments (Rubin & Ikeda, 1990) and stratigraphy analysis (Rubin & Hunter, 1985;Bristow et al., 2000Bristow et al., , 2005Bristow et al., , 2007 have revealed evidence that lateral migration has a significant role in the formation of linear dunes (e.g. Livingstone et al., 2007), largely discrediting the roll vortices hypothesis. Figure 1 shows an example of a GPR cross-section of a linear dune (Bristow et al., 2005), together with the stratigraphic interpretation identifying sets of cross-stratification and bounding surfaces.
The interpretation indicates a lateral migration and growth of the dune over time. Instead of formation through lateral shifting, linear dunes may also grow and elongate in the crest parallel direction by sand transport along the dune flanks driven by secondary airflow dynamics in the lee of the crest (Tsoar, 1982;Tsoar et al., 2004).
Through the coalescence of smaller segments, linear dune fields evolve over time towards a coarser state with fewer and larger dune ridges, and fewer junctions. This is how dune populations behave as a self-organized system. Patterning of dune fields has been measured on various planetary bodies, for example, on Earth (Thomas, 1986;Kocurek & Ewing, 2005;Derickson et al., 2008;Rachal & Dugas, 2009), on Mars (Parteli & Herrmann, 2007;) and on Titan (Savage, et al., 2014;Ewing, et al., 2015), and Day & Kocurek (2018) have found a constant spatial density of bedform interaction by analyzing dune fields on multiple planetary bodies. The patterning of dunes was also modelled with numerical algorithms focusing on defects and Y-junctions (Werner & Kocurek 1999;Diniega et al., 2010). Pattern coarsening and bedform interaction have been studied in terms of surface topographies and development, but less so with respect to internal architectures (Day & Kocurek, 2017) and mass exchange. More recently, Telfer et al. (2017) investigated the relationship between the age and pattern coarsening of linear dunes, and Brothers et al. (2017) combined GPR profiles with remote sensing images to interpret bedform interactions, identifying two types of bedform repulsion. By observing deflated stratigraphy, Day & Kocurek (2017) defined the key features of repulsing elements: a target and an impactor. Dating samples, GPR surveys and observations in the field, however, only reveal a static picture of bedding facies that require manual interpretation.
The morphological development of linear dunes has been simulated with three-dimensional (3D) numerical computer models since the work of Werner (1995), particularly with reduced-complexity approaches such as Cellular Automata (CA) and Lattice-Boltzmann Methods (LBM) (Gao et al., 2015a,b;Rozier et al., 2019;Gadal et al., 2020).
The internal stratigraphy of aeolian bedforms has been modelled using computational techniques, including geometric visualizations (Rubin, 1987), ripple stratigraphy (Forrest & Haff, 1992), grain size sorting on ripples (Anderson & Bunas, 1993) and 1D chronologies of linear dunes (Telfer et al., 2010;Leighton et al., 2014). Recently, Swanson et al. (2019) employed a numerical 2D transect model to explore internal stratigraphy and set-scale architecture of aeolian bedforms, and Zhang et al. (2012) simulated internal stratigraphy of star dunes. The full 3D sedimentary structure of whole linear dunes and dune fields, however, has not been simulated previously. The work here presents insights into the 3D internal stratigraphy of linear dunes and their bedform interactions. These results enable a clearer differentiation between types of bedform collisions as previously identified mainly from transverse dunes and 2D profiles.

Sand movement
The basis of the CA here is the dune algorithm of Werner (1995), which has been applied in many other simulation studies (Baas & Nield, 2007;Barchyn & Hugenholtz, 2012;Liu & Coulthard, 2015;Yan & Baas, 2017;Mayaud et al., 2017). In this algorithm, dune topography is composed of sand slabs, of a fixed thickness, stacked on a grid, and elevations are changed by removing or adding slabs. A sand movement cycle starts with the polling of a random cell on the grid. The top slab is removed with an erosion probability, p e , which can be set to 1 when simulating bare sand dunes. The eroded slab is moved over a fixed 'jump-length' distance, L, downwind to a new cell location. Here a deposition probability, p d , is applied to decide whether the slab is deposited. If not, the slab jumps another distance L downwind and the deposition check is repeated. The deposition probability is greater over sandy beds than over rocky surfaces to represent the differing saltation rebound properties. In addition to the fundamental slab movement rules, a shadow zone in the lee of sand piles simulates the sheltering effect, preventing local erosion and trapping all incoming slabs. After erosion or deposition, the local elevation change alters the slopes with its neighbouring cells. An angle of repose for loose sand is enforced to initiate avalanching when the angle is greater than the threshold, so as to prevent unrealistically steep slopes. The angle of repose is typically set to around 30 degrees, depending on the discretized angles available as a function of slab thickness (Narteau et al., 2009). In this study, the angle of repose is set as 27 degrees and the angle of initial yield as 35 degrees being slightly greater (Kocurek & Ewing, 2016). The two thresholds are recognized by previous studies on sand pile dynamics (Makse et al., 1997), and are dictated by the discrete slab-height of the algorithm. In an algorithm cycle, after being polled, possible endings for a slab therefore are: left remaining on its original cell (no erosion), deposition on a downwind cell, or coming to rest after being part of an avalanche event.
The Werner algorithm does not include the effects of secondary airflow dynamics, such as flow deflection or roller vortices, and is therefore not able to replicate the type of linear dune elongation via sand movement along the flanks by secondary helical flow, as is proposed by Tsoar et al. (2004).

Delineation of the internal structure
The internal stratigraphy of dunes is reconstructed by labelling slabs when they are moved or when they are exposed. A complete sand slab transport cycle contains three sedimentary processes where the elevation changes: erosion, deposition and avalanche. When a slab is removed from a cell during the erosion stage, the next slab in the stack on that cell becomes exposed to the air and is labelled as an erosional surface. The slab that is moved is subsequently labelled as a (pure) deposited slab if it is dropped on a destination cell without triggering an avalanche, or as an avalanched slab if it is involved in this process. All other slabs involved in an avalanche are also labelled as avalanched slabs. Further, in order to record a chronology, burial age is marked on each slab. Slabs moved or exposed are marked by the current iteration number, so that afterwards this iteration label reflects the time of burial. This allows entire views of chronology records to be mapped on a 3D basis horizontally throughout dune fields and vertically reaching the innermost parts of dunes.

Implementation of oblique winds
To produce linear dunes, bi-modal oblique winds are applied to the model domain. There are a few different methods for implementing oblique winds over a grid of cells, including: (i) rotating the entire domain as a continuous surface associated with resampling after each wind shift (Bishop et al., 2002;Nield & Baas, 2008); (ii) moving slabs in combinations of discrete lateral and forward steps (Werner, 1995;Nield & Baas, 2008); and (iii) allowing slab fractions to distribute over multiple cells in accordance with a wind path length (Mayaud et al., 2017). In the study here, wind directions are not changed continuously but remain at the same obliqueness during a simulation. This enables the cells to be defined as rectangular, rather than square, with a fixed 'side to side' ratio of 2 to 1 (lateral width to downwind length). Wind is then simulated as transporting slabs along the diagonals of the cells, so that √5 is the fixed jump length unit and the angle between two opposing oblique winds is 128 degrees. The rotated shadow zone fully covers the lee slope under this direction setting. This solution avoids surface resampling issues and jump length limitations, because wind directions and jump length can be easily changed by reshaping cells and adjusting gradients.

Definition of the timescale
In the bare-sand algorithm, time and space dimensions are not fixed to a specific scale. Time evolution is tracked in terms of iterations, I, defined as a number of slab erosion polling events equal to the number of cells in the domain (Baas & Nield, 2007;Eastwood et al., 2011). In order to observe effects caused by the bi-modal wind regime, a 'model year' is manually defined here as 500 iterations.

Parameters and simulation conditions
The slab height, h, is set to one tenth of the cell unit (as in : Nield & Baas, 2008;Eastwood et al., 2011) to provide the vertical resolution. Deposition probabilities are set to 0.4 for bare surface cells (mimicking a hard rock base) and 0.6 for cells with sand slabs present. The angle of initial yield is distinguished from the angle of final repose, which are seven-slab and five-slab differences between adjacent cells, respectively. The shadow zone angle is set to 15 degrees, as has been standard in previous studies. The jump length is set to the minimum unit of √5, along the cell diagonals, to enable pattern formation at the smallest scale.
Simulation scenarios start with a flat horizontal sand bed overlying a non-erodible base.
Bi-modal oblique wind regimes are applied, with both winds at the same sand carrying capacity, but with varying ratios of durations per model year. The standard domain orientation and descriptions, as shown in subsequent plan view figures, is defined as top ('north') being upwind, bottom ('south') being the resultant downwind, with the first oblique wind blowing from top-left ('north-west') to bottomright ('south-east'), and the second wind blowing from north-east to south-west.
The domain-averaged potential volumetric transport rate can be quantified as the product of the volume of a slab and its mean path length ((p e /p e )L) per iteration, transported per lateral span (Nield & Baas, 2008). If the length unit of the domain is defined as a metre, and given that a model year is defined as 500 iterations, the volumetric transport rate for the simulation conditions and parameters used in this study amounts to a cumulative equivalent of 210 m 3 m À1 yr À1 . This is on the same order of magnitude of potential sand transport rates reported for the linear dune zones of the Namib Sand Sea by Lancaster (1985).

Simulation settings
The simulations aim at exploring two aspects: (i) bedform interactionsincluding sequences naturally extracted from pattern coarsening processes, supplemented with simplified binary dune collision to provide further insight; and (ii) how different bi-modal wind regimes affect the linear dunes and their internal structures.
Three wind regimes are simulated with different ratios of two oblique winds. The two winds have identical sand carrying capacity but are applied for differing durations within a fixed model year of 500 iterations. The first regime is a perfectly balanced one where each of the two oblique winds blows for 250 iterations in turn, hence termed as N = 1 regime, where N is the wind duration ratio (Gao et al., 2015b). This regime yields a resultant sand transport that is longitudinal, i.e. parallel with the resultant crest lines, moving towards the 'south' of the domain. Along the axis normal to the crest, sand is moved back and forth in equal amounts with a net total displacement of zero. The second regime is one where the two oblique winds blow slightly unequally: for 240 and 260 iterations (N = 1.08), and the third regime consists of winds blowing for 150 and 350 iterations (N = 2.33). The second regime yields a resultant transport direction that is 1.4 degrees deflected from the longitudinal axis, while the third regime yields a transport direction that is 38 degrees deflected. These two asymmetrical regimes were selected for three reasons. First, the asymmetries of 20 iterations difference (for N = 1.08) and 200 iterations difference (N = 2.33) represent two orders of magnitude of net lateral sand transport. Second, even the slight asymmetry of only 20 iterations difference leads to significant effects on dune morphology and internal structure over an extended timescale. Third, because 'linear' here is a morphological term, dunes that are created by the N = 1.08 regime are fully longitudinal (resultant transport direction is within 5 degrees of the crest line axis), while the N = 2.33 regime generates dunes that are longitudinal, but turn oblique at later stages. These differences reflect the 'fingering' and 'bed instability' controls towards dune orientation (Courrech du Pont et al., 2014).
The simulation domain is a grid of 200 cells wide by 1800 cells long in the resultant transport direction. Because of the rectangular cell geometry this yields a domain of 400 9 1800 length units. The domain is intentionally long in the resultant transport direction in order to avoid having to use periodic boundaries, which would lead to recycling of sand and migrating dune ridges that disrupt and interfere with the development of the internal stratigraphy, thus not representing natural conditions. Instead, the very long domain is initially populated with a 30-slab thick flat sand bed covering only a spanwise band from a streamwise fetch distance of 100 to 500 cells (hence covering a square surface). This domain provides sufficient space for dune development, maturation and elongation starting from the square layer of sand without the dune tips extending beyond the edge of the grid. The lateral edges of the domain, however, are connected as periodic boundaries, so that dune ridges are preserved as they migrate laterally.
For each wind regime simulation, during the initial pattern coarsening stages the bedform interactions and Y-junction dynamics were captured and analysed. The analysis and interpretation of Y-junction dynamics was supplemented with simplified 2D simulations of bedform interactions and mass exchange between two interacting dunes, as a means to elucidate the complexity and stratigraphy of Yjunction dynamics. To investigate the impact of wind regimes on stratigraphy in more mature dunes, the authors select a specific developed ridge after the dune field has coarsened substantially and most of the early bedform interactions have concluded. This ridge is then isolated (all other sand topography is removed from the domain) and the wind regime simulation is continued.

RESULTS AND DISCUSSION
The scenario simulations yield small-scale linear dunes, of which the sedimentary structure is not significantly occupied by any low angle plinth. They migrate laterally back and forth under the alternating oblique winds, while extending at their downwind tips. The migration is driven by avalanching and the dune flanks are at the angle of repose. This is not representative of typical real-world linear dunes, where slip faces occupy only the upper parts of the flanks and the base or plinth of the dunes extend much further laterally into the interdune areas, with lower topographic slopes. The linear dunes modelled here are thus more comparable with those that migrate sideways. This model limitation is shared with other CA and LBM simulation studies that exclude secondary-flow dynamics, where dune flanks are typically at the angle of repose (Rozier et al., 2019;Gadal et al., 2020). The internal sedimentary structure of the linear dunes is jointly shaped by bedform interactions and the wind regime. At the earlier, pattern coarsening stage, collision stratigraphy is significant. As the dune field matures, dunes become more isolated with reduced connectivity and interactions, at which point the internal structure is determined more by the wind regime. Particularly in the asymmetrical wind regimes the earlier collision facies are then gradually exposed and reworked into newly buried deposition stratigraphy.

Pattern coarsening and the internal structure of Y-junctions
Pattern coarsening From the start of a simulation, small heaps of sand are initiated and develop into a small-scale bedform pattern. Interactions between the bedforms lead to a reduction in number of crest lines, growth in dune height and pattern coarsening.
In the earliest stage the initial dunes are small enough to migrate whole dune widths and be reworked within single oblique winds and any long-term internal stratigraphy is not preserved.
The longer duration component of each bi-modal regime produces a total sand displacement of: 105 m 3 m À1 for 250 iterations (N = 1), 109 m 3 m À1 for 260 iterations (N = 1.08), and 147 m 3 m À1 for 350 iterations (N = 2.33). These displacements represent dune heights of 7.3 m, 7.5 m and 8.6 m, respectively. Simulation results show that the bedforms reach these heights within two to three model years, after which their morphology and internal stratigraphy start to reflect the full bi-modal wind regime. The simulated dunes are very similar in height and spacing as the (now vegetated) linear dunes visible in Rice Valley, California, USA (N33°57´, W114°45´).
The signatures of pattern coarsening during the subsequent model years are mostly preserved in the stratigraphy, but eventually all sand in the domain is reworked over longer periods of time, due to lateral migration and bedform interactions. The sand from the original flat deposit at T = 0 is completely reworked after 10 years under the N = 2.33 regime, after 20 years under the N = 1.08 regime, and under the N = 1 regime some original deposit survives at the core of the developing linear dunes for up to 47 model years. Figure 2 shows for each wind regime a sequence of morphological development over the course of the simulation (the time steps shown are optimized for illustrating various Yjunction dynamics, see next section). For the two more symmetrical wind regimes (N = 1 and N = 1.08) the linear dunes are aligned along the longitudinal axis, with only the downwind tips deflecting along an individual wind direction. The sequence for the N = 2.33 regime, however, shows a gradual change from an early stage where crest lines are largely longitudinal (running 'north' to 'south') towards a later stage where dune ridges align along the resultant sand transport direction for this wind regime (a deflection of 38 degrees to the 'west' or 'left' of the domain). This can be attributed to a phase shift from a fundamental bed instability mode related to the dominant lateral sand transport components, to a resultant transport or 'fingering' mode related to ridge extension at the tip (Courrech du Pont et al., 2014), once the density of Y-junctions is sufficiently reduced in the domain. At this stage, the dunes are so large that their turnover time becomes longer than the wind cycle and they move independently without making stable contact with adjacent sand bodies. Bedform interaction stratigraphy -Y junction and dune tip dynamics Dune terminations and Y-junctions are key elements in bedform interactions that lead to pattern coarsening. This study identifies six fundamental types of interactions, conceptually represented in Fig. 3, that underlie pattern coarsening in linear dune fields. Examples of these dynamics are labelled in Fig. 2 and can also be observed in the simulation animations provided in the Supplementary Materials. For the six types of interactions identified, signature features of their stratigraphy, with extracted profiles shown in Fig. 4, supported with simplified 2D simulations to track sand mass exchanges, shown in Fig. 5. The core process is lateral interaction of dune ridges caused by differential migration speeds related to dune heights or parts of a dune ridge that are lower than other sections, with lower sections or dunes moving faster than larger ridges. The crest-splitting type (1) is a longitudinal process that breaks up a long dune ridge along significant variations in dune height (a 'saddle'). The lower section of the ridge loses sand to the higher parts, or in other words: the lower section moves faster in the longitudinal direction. The ridge is pinched out until the crest splits, creating two new terminations that can subsequently operate in other interaction types. Crest-splitting is also reported in transverse dunes, in terms of defect creation  linked to a transverse instability mechanism. The longitudinal nature of the crest-splitting in the simulations herein is however unique to linear dune field dynamics. The oscillating interactions (types 2 and 3 in Fig. 3) are a reflection of the bi-directional transport regime and are not observed under unidirectional conditions. Y-junctions may or may not form, depending on the extent to which the oscillating tip approaches the neighbouring ridge. The oscillating interactions do not leave significant or long preserved structures in the stratigraphy. The other three types (4, 5 and 6) may be grouped as fundamental 'reorganization' dynamics, and include merging, cannibalizing and defect repulsion, following the nomenclature established for transverse dune interactions by Kocurek et al. (2010). The bedform interactions include two processes: (i) sediment transfer over a distance between two separate dunes, without permanent contact; and (ii) superposition when two main dune bodies come into permanent contact. The oscillating interactions (types 2 and 3) occur where interdune spacing is too wide relative to the size/speed of the dunes to enable permanent contact and superposition, and instead the smaller, faster-moving dune termination dissipates by feeding all of   Fig. 4). H1 and H2 refer to the (absolute) crest elevation of the target and impactor, respectively. its sand to the neighbouring ridge. The three reorganization dynamics (types 4, 5 and 6) all involve superposition.

Sediment transfer over a distance
The first step in developing lateral sediment transfer is usually a crest split, which produces fast-moving dune terminations at either side of the gap (Fig. 2A, type 1). The stratigraphy example in Fig. 4A shows that the splitting dune ridge in Fig. 2A consists of two older cores (disregarding the left-most peak, which is an artifact of the upwind edge of the original sand deposit at time = 0). The northern core is younger in age and displays low-angled bedding planes, whereas the southern core traces back to 'original' sand at its northern tip and shows bedding planes that are inclined at slip-face angles. The upper part of the stratigraphy is a layer of active sand that has been deposited in the most recent year, with a topographic low where the older cores connect. Figure 2A shows how the sand in the topographic low moves southward, leaving a gap between two old cores and creating the crest split. The newly created terminations oscillate laterally, with two possible outcomes: (i) the terminations lose sand, sometimes via ephemeral Y-junction connections, until the entire sand body disappears (types 2 and 3); or (ii) one of the terminations permanently connects with a neighbouring dune ridge, forming a stable Yjunction, either in reversed (northern tip connects) or normal shape (southern tip connects) (types 4, 5 and 6).
The examples of interaction types 2 and 3 outlined in Fig. 2A and B show oscillating terminations intermittently connecting to a main dune ridge to the 'east'. In the asymmetrical wind regimes the more dominant 'south-west' moving transport component is able to disconnect the termination away from the ridge every year. Due to the lack of a stable and permanent body contact the oscillating interactions do not leave any significant evidence in the stratigraphic record. The oscillating termination simply acts as a sand source for the growth of the main dune ridge, and the ephemeral Y-junctions are assembled and disassembled seasonally. This may explain why the interaction structure is less recognized and documented in linear dune fields than in transverse dune fields. Wind seasonality and the related stages of dune field evolution produce features that only exist seasonally and are not preserved in the depositional record. The stratigraphy is visible in the cross-sections of Fig. 4B and C: the main dune ridge has an internal age structure and stratigraphy that is similar to a solitary dune ridge, showing 'old' sand at its core and active sand at the surface, a growth of annual accumulations on its flanks, and buried erosion surface layers that are inclined at the same angle as the dune flanks. The oscillating termination, on the other hand, is composed purely of active sand that has been moved within the most recent year and with no further internal stratigraphy. If a termination permanently connects with a neighbouring dune ridge ('ii' above) the sediment transfer progresses to superposition, whichunlike for ephemeral Yjunctionsdoes result in stratigraphy facies.

Superposition
The three bedform reorganization interactions occur after superposition starts. Examples of the different types (4, 5 and 6 in Fig. 3) are indicated in Fig. 2, with stratigraphic cross-sections at two time-steps shown in Fig. 4. To help to clarify the stratigraphy facies Fig. 5 shows a 'schematic'. They are simplified 2D uni-directional simulations of the interaction between a faster moving 'impactor' dune with a crest elevation, H2, and a 'target' dune with a crest elevation, H1, that can be related to previous studies (e.g. Kocurek & Ewing, 2005;Brothers et al., 2017). The mass exchange between interacting bodies are visualized.
Throughout a superposition process, a recognizable feature is the interaction surface with inclinations much lower than the angle of repose. In the merging cross-sections (Fig. 4D), the interaction surface is visible as a shallowangled base extending to the right of the main dune ridge over which the smaller approaching dune crest is migrating to the left. This base consists of sand of the main dune that has been kept in place by the shadow zone extending from the approaching dune crest under the dominant left-ward transport component of the bi-directional regime. The smaller dune, consisting only of active sand, climbs up the main ridge along the interaction surface, and its material is redistributed along the flanks, as indicated by Fig. 5A. The internal stratigraphy of the main ridge shows left-dipping slip-face strata, but the interaction surface remains visible as a shallowsloping erosion surface layer dipping towards the right in the stratigraphy of Fig. 4D, step 2. As the main ridge is migrating to the south-west under the asymmetrical regime, the lightly-buried interaction surface is however quickly exposed and eroded and is not preserved in the long-term. This merging type interaction occurs when the approaching impactor dune crest (to the east in the example labelled in Fig. 2D) is much smaller than the target dune ridge (to the west in this example), as is also evident in Fig. 5A.
An increasing impactor : target size ratio leads to a cannibalizing interaction. Both merging and cannibalizing results in the reduction of dune ridge number. The key differentiating point between merging and cannibalizing is whether the impactor crest can exceed the height of the target dune while the former climbs up the flank of the latter, as can be seen in comparing Fig. 5A step 2 with Fig. 5B step 2. When this happens the top of the target dune is eroded away and the impactor dune overtops the older base of the target dune that is left behind. The interface surface is relatively flat and separates the remnant slip-face stratigraphy of the lower part of the target dune below from the active and unconsolidated sediment of the cannibalizing impactor dune above. Because the interface surface extends through the whole width of the combined dune ridge it is preserved for longer than in the merging type interaction.
As the heights of the impactor and target dunes approach parity repulsion can take place (Figs 4F and 5C). The example of this interaction indicated in Fig. 2C shows that the two ridges are of nearly equal size, similar to the size ratio in Fig. 5C. In this process the impactor dune starts to climb over a low-angled plinth of the target dune, but as the impactor intercepts the dominant incoming sand drift the target dune loses volume, shrinks in size and escapes (or is repulsed) away from the impactor. As the target escapes, the base that is left in its wake thins out. Since the growing (and now slower) impactor dune migrates over this base, the resulting stratigraphy displays an interaction surface that dips down to the left in Fig. 5C, also visible as an inclined erosion surface layer in Fig. 4F step 2. This facies is only fully destroyed if the dune migrates laterally through an entire width's distance, and so this interaction surface is preserved for the longest and has more probability of being encountered in a deposition record. This may be one reason why repulsions are well-documented in transverse and barchan dune fields (e.g. Brothers et al., 2017). The interaction surfaces visible in Figs 4 and 5 may be similar to the bounding surfaces identified in the GPR profile in Fig. 1 (Bristow et al., 2005). Day & Kocurek (2017, fig. 1) provide a conceptual diagram and field evidence of the plan view stratigraphy of a defect repulsion. In comparison, this study shows in Fig. 6 the plan view stratigraphy of a horizontal slice through the repulsion described above. Both post-repulsion and ongoing repulsion are reflected in this frame, as the northern part of the target dune ridge has separated, whilst the southern part is still attached to the Y-junction. Figure 5 shows the situation just before the southern part breaks free towards the west (see Supplementary Materials for full animation sequence). The repulsion process has left a stratigraphic signature as a diagonal line of erosion surfaces at the interface between strata dipping in different directions. The cross-strata in the northern part of the impactor dune to the right dip down the slip-face flank (towards the north-west), while the strata at the Y-junction itself dip towards the south-west aligned to the resultant net transport direction under this N = 2.33 wind regime. The interface between these two sets is left behind in the stratigraphy of the impactor dune as the target dune ridge is tearing away while the Y-junction migrates towards the south. The interface between strata dipping in different directions resembles the plan view stratigraphy reported in Day & Kocurek (2017), except that here the asymmetrical oblique winds drive repulsions that migrate diagonally to the linear dune crest line, rather than perpendicular as would be the case in uni-directional transverse dune dynamics.
When comparing the internal stratigraphy of the three superposition interaction types the key distinction is the slope and extent of the interacting surface that separates the active sediment of the impactor dune from the buried slip-face strata of the target dune. In the merging type this surface slopes parallel to the upwind flank at a small depth and is reworked relatively quickly. In the cannibalizing type this surface is relatively horizontal and is buried more in the centre of the ridge, taking a longer time to be exposed and reworked. In the repulsion type the interaction surface dips downwind and lies very close to bedrock, surviving the longest.
The simulation results presented here of the internal stratigraphy of bedform interactions enables a clearer differentiation between the types of bedform collision as previously identified mainly in transverse dune fields (e.g. Kocurek et al., 2010;. In particular, the internal stratigraphy can reveal signatures of past bedform interactions, specifically from merging and cannibalizing, that are otherwise not visible from the external dune shape alone. It should be recognized, however, that in all of the above stratigraphic signatures the model algorithm does not simulate the potentially significant impact of lee-side scour which may remove the entire base of the target dune stratigraphy depending on the scale ratio with the impactor dune (Kocurek & Ewing, 2016;Cardenas et al., 2019).

Internal structure of independent linear dunes
As shown above, dune interactions leave signatures in the internal stratigraphy that persist for some time. Through pattern coarsening the dunes eventually become so large, however, that their turnover time exceeds the duration of a wind cycle, interactions become increasingly rare and the dunes are independent. Under further growth and reworking the wind regime then becomes the sole determinant for internal stratigraphy. Figure 7 shows the internal stratigraphy of an isolated dune for each wind regime. The dunes have been subject to the wind regime, after isolation, for a long enough duration to ensure complete re-working of their sediment, thus excluding any bedform collision signatures in their stratigraphy. Their orientation is dominated by tip-extension, with the terminations pointing in the resultant sand transport direction. Because the dunes have been losing sand from their tips while migrating, however, they have reduced in length (compared to before they were isolated, cf. Fig. 2). The large transverse size of the dunes, meanwhile, means that the annually transported volume of sediment is distributed over large flanks and this leads to thin layers of sand deposits that are less distinguishable than their counterparts in the earlier coarsening phase (Fig. 4). The deposition imagery is thus not suitable for displaying internal chronology, particularly for the N = 1 and N = 1.08 wind regimes with their lower net transport. Instead, in Fig. 7 chronology is visualized by tracking deposition times within periods of multiple years. For N = 1 (Fig. 7A), sand ages on a 10-year cycle are mapped, showing that the internal stratigraphy covers the most recent 35 years. In plan view, the deposition sequence shows a symmetrical, chevron-line pattern, and in the cross-section the chronology is also symmetrical, with the deposition from the last oblique wind covering the left flank. For N = 1.08 (Fig. 7B), chronology is also mapped on a 10year cycle, and the stratigraphy is shown for two distinct time steps, the first after 25 years of wind regime simulation after isolation, and the second after a further 20 years beyond this. These two time steps are shown to demonstrate how the stratigraphyplan view chronology in particulardevelops towards more symmetry over time. Figure 7B1 plan view chronology reflects a nearly unidirectional pattern, with much wider deposition bands on the left flank. After the additional 20 years of migration and development, the plan view chronology is much more symmetrical and starts to approach that of the N = 1 regime (Fig. 7A). The change of symmetry is attributed to the change of sand supply. Before the isolation, adjacent dunes trap all sand moving towards their neighbour, whilst afterwards, sand recycled by the periodical boundary feeds the dune. For the N = 2.33 regime, Fig. 7C visualizes deposition times within a two-year cyclea much short cycle than the others because this dune migrates laterally very quickly and the internal chronology extends back for only five years. The chronology shows how the deposition is entirely unidirectional, migrating strongly to the left under the much more dominant second wind of the regime.
The maps and cross-sections of within-year (seasonal) deposition time mainly reflect the back and forth reworking of sediment under the two alternating winds. Under the N = 1 and N = 1.08 regimes, most of the sand is reworked within the year and only the deposits from the beginning and end of each year are preserved because they are left buried by the following wind. Under the N = 2.33 regime the deposits are mostly from 150 to 350 iterations within the 500-iteration year, because the strongly asymmetrical wind regime preserves only deposits that are laid down during the middle part of the year.
For slab types, particularly for erosion surfaces in cross-section, the visualization is similarly affected by the small volume of net sand movement, and differentiation between layers from consecutive years is difficult, again especially for N = 1 and N = 1.08. Slab type images of cross-sections and plan views in Figs 6B and 7A apply a 10-year cycle, within which the last five years is left invisible. Thus, visible surfaces are stacked up with five years of record, and between each five-year band sits an (invisible) adjacent five-year interval. With this method the internal stratigraphy can be visualized reflecting the migration and elongation of the dune towards the south. In the cross sections, erosion surfaces dipping to the left are the buried signatures of the temporary (erosion) stoss slopes existing under the weaker wind in the regime. The evidence of erosion surfaces in Fig. 6 indicate that for the N = 1 regime, the surfaces stack-up dipping symmetrically left and right. Under the N = 1.08 regime the surfaces stack-up less symmetrically, and under the N = 2.33 regime all erosion surfaces dip down left-ward, indicating the strong lateral migration of the ridge.
In the longitudinal profiles of the dunes the green erosional surfaces are on the original annual basis, as the longitudinal transport and deposition is a continuous and more voluminous process creating sufficient spacing between layers. These longitudinal cross-sections are reminiscent of the empirical profile data collected by Telfer (2011), reporting an example of linear dune elongation using OSL dating. Comparison to GPR profiles (such as Bristow et al., 2000) is more challenging because these data are complicated by sinuosity signatures. Under N = 1 and N = 1.08, the elongation and migration signatures dominate the lower parts of the dunes. They show full sequences of annually formed erosional surfaces, between which are the grain flow strata deposited during the year. The youngest sediments are distributed around the south termination. Under the N = 1 regime the upper part of the stratigraphy is dominated by a mix of grain fall deposits and erosion surfaces, reflecting the cyclic movement of sand in two oblique directions. The stratigraphy for the N = 1.08 regime in the earlier phase is more complex and shows two distinct bands of erosion-surface and grain flow striping, suggesting a secondary surface undulation progressing along the top of the main linear ridge. At this phase, the region of mixed grain fall deposits and erosion surfaces is confined to the top layer of the dune, reflecting the fact that the dune is slowly migrating laterally under the slightly asymmetrical wind regime. At the later phase (Fig. 7B2) the dune has reworked itself significantly and the internal stratigraphy in slab types approaches that of the N = 1 regime.
The longitudinal stratigraphy of the N = 2.33 regime is distorted, because the linear dune has become strongly curved and nearly crescentshaped, under the highly asymmetrical wind regime and as it has shrunk in length, while the transect runs straight north-south through the body. The authors tried extracting a curving transect that follows the crest line, but the resulting stratigraphy is very noisy and disjointed. The concave banding in the longitudinal stratigraphy of Fig. 7C reflects the varying height of the deposits and erosion surfaces of the curving dune as it rapidly migrates laterally to the west. The clear separation between erosion surfaces (left behind under the weaker secondary wind) with thick layers of grain flow strata reflects the nearly unidirectional nature of the dune migration.
The internal stratigraphies of the linear dunes display features of two models of development. In the (near-) symmetrical wind regimes (N = 1 and N = 1.08) the cross-section stratigraphy indicates strata dipping parallel to the dune flanks and interlinking at the crestssignatures of the mainly vertical accumulation of sediment layers on each of the flanks under the alternating oblique winds. This conforms to the model proposed by Bagnold (1941) and supported by evidence of McKee & Tibbitts (1964). The visual resolution of the model results is not high enough, however, to reveal the precise details of the inter-meshing of the tops of the dipping strata at the crest points. In the strongly asymmetrical wind-regime (N = 2.33) the stratigraphic signatures of lateral migration dominate and resemble the cross-strata found in transverse dunes, much like the model results of Rubin & Hunter (1985) and shown in the GPR evidence of Bristow et al. (2007). Because the cellular automaton (CA) model does not include secondary airflow features it is not possible to simulate the type of sinuous crest lines of linear seif-dunes that could be achieved for the kind of stratigraphy proposed by Tsoar (1982).

CONCLUSIONS
The simulation results reported here reveal a variety of stratigraphic signatures as a consequence of linear dune field pattern coarsening via bedform interactions, as well as more mature independent linear dunes. During the pattern coarsening phase six types of bedform interactions were identified in terms of morphological development. Three of the types, longitudinal crest-splitting and oscillating interactions forming normal or reverse Y-junctions, are unique to linear dune field development but they do not leave much evidence in the stratigraphic record because the interactions are mediated by sediment transfers conforming to the normal erosion, transport and deposition processes. The other three types, however, involve the superposition of an impactor dune body onto a target dune body, occurring during merging, cannibalizing and repulsion. These interactions do leave a signature in the stratigraphy in the form of a buried interaction surface, and the inclination of this surface in cross-section reflects the type of superposition involved, governed by the impactor : target size ratio of the two original dune bodies.
The stratigraphy of more mature independent linear dunes conforms to two different models of accumulation: under symmetrical wind regimes the flanks build up under the alternating oblique winds, leaving strata parallel to the flanks dipping in both directions; whereas under strongly asymmetrical wind regimes the stratigraphy is dominated by lateral migration, resembling the cross-strata of transverse dunes. Even though linear dunes under asymmetrical wind regimes extend at their tips oriented along the net resultant sand transport direction, their internal stratigraphy reflects mainly lateral migration.
YL is now pursuing a PhD degree at the University of Oxford. We thank Charlie Bristow for helpful discussions during the earlier stages of development, and for his kind permission to reuse Figure 1. We thank R Ewing, M Day and an anonymous reviewer for extensive and thoughtful comments that have helped improve the text and figures.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.